In-medium bound states of two bosonic impurities in a one-dimensional Fermi gas

TL;DR

This paper explores how a one-dimensional Fermi gas with two bosonic impurities can form in-medium bound states due to induced attractive interactions, with exact and effective models analyzing various interaction regimes.

Contribution

It introduces an exactly solvable model for impurity interactions in a 1D Fermi gas and develops an effective model to study non-integrable cases with unequal masses or potentials.

Findings

Attractive impurity interactions lead to in-medium bound states.

Repulsive interactions do not produce bound states.

Effective models agree with exact solutions for certain parameters.

Abstract

We investigate the ground-state energy of a one-dimensional Fermi gas with two bosonic impurities. We consider spinless fermions with no fermion-fermion interactions. The fermion-impurity and impurity-impurity interactions are modelled with Dirac delta functions. First, we study the case where impurity and fermions have equal masses, and the impurity-impurity two-body interaction is identical to the fermion-impurity interaction, such that the system is solvable with the Bethe ansatz. For attractive interactions, we find that the energy of the impurity-impurity subsystem is below the energy of the bound state that exists without the Fermi gas. We interpret this as a manifestation of attractive boson-boson interactions induced by the fermionic medium, and refer to the impurity-impurity subsystem as an in-medium bound state. For repulsive interactions, we find no in-medium bound states. Second, we construct an effective model to describe these interactions, and compare its predictions to the exact solution. We use this effective model to study non-integrable systems with unequal masses and/or potentials. We discuss parameter regimes for which impurity-impurity attraction induced by the Fermi gas can lead to the formation of in-medium bound states made of bosons that repel each other in the absence of the Fermi gas.

Equations

$$H=H_f+H_b+V_{fb},$$

$$H_f=-\sum _{j=1}^{N_f}\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial X_j^2},$$

$$H_b=-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial Y_1^2}-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial Y_2^2}+g_{II}\delta (Y_1-Y_2),$$

$$V_{fb}=g\sum _{i,j}\delta (Y_i-X_j),i=1,2;j=1,\dots ,N_f,$$

$$h_{BA}=-\sum _{j=1}^N\frac{{\partial }^2}{\partial x_j^2}+2c\sum _{j<l}^N\delta \left(x_j-x_l\right),$$

$$\begin{array}{rl}& e^{i{k}_{j}l}=\frac{k_j-{\Lambda }_1+\frac{ic}{2}}{k_j-{\Lambda }_1-\frac{ic}{2}}\frac{k_j-{\Lambda }_2+\frac{ic}{2}}{k_j-{\Lambda }_2-\frac{ic}{2}},1\le j\le N;\\ & \prod _{j=1}^N\frac{k_j-{\Lambda }_1+\frac{ic}{2}}{k_j-{\Lambda }_1-\frac{ic}{2}}=1,\prod _{j=1}^N\frac{k_j-{\Lambda }_2+\frac{ic}{2}}{k_j-{\Lambda }_2-\frac{ic}{2}}=1;\end{array}$$

$$\begin{array}{rl}{k}_1& \simeq \sqrt{\frac{3c}{l}},k_2\simeq -\sqrt{\frac{3c}{l}},k_3\simeq 0,\\ k_j& \simeq k_j^{(0)}+\frac{2c}{k_j^{(0)}l},\text{for }4\le j\le N;\\ {\Lambda }_1& \simeq \sqrt{\frac{c}{l}},{\Lambda }_2\simeq -\sqrt{\frac{c}{l}},\end{array}$$

$$h_{BA}\to -\frac{J}{2}\sum ({\sigma }_x^i{\sigma }_x^{i+1}+{\sigma }_y^i{\sigma }_y^{i+1}),$$

$$h_{BO}=-\frac{{\partial }^2}{\partial x^2}+2c\left[\delta \left(x-\frac{r}{2}\right)+\delta \left(x+\frac{r}{2}\right)\right],$$

$${\mathcal{E}}_{static}(c)=\left(\pi +\frac{c^2}{\pi }\right)\arctan (\frac{c}{\pi })+c-\frac{c^2}{2}.$$

$$h_{\mathrm{eff}}=-\frac{m}{m_{\mathrm{eff}}}\frac{{\partial }^2}{\partial y_1^2}-\frac{m}{m_{\mathrm{eff}}}\frac{{\partial }^2}{\partial y_2^2}+W(y_1-y_2),$$

$$\frac{m_{\mathrm{eff}}}{m}\simeq 1+\frac{c^2}{{\pi }^4}+\left(\frac{2}{{\pi }^6}-\frac{1}{6{\pi }^4}\right)c^3+\left(\frac{4}{{\pi }^8}-\frac{1}{2{\pi }^6}\right)c^4.$$

$$W(y_1-y_2)=2c\delta (y_1-y_2)+V(y_1-y_2),$$

$$E_{\mathrm{eff}}\simeq -\frac{m_{\mathrm{eff}}}{8m}{\left[{\int }_{-\infty }^{\infty }W(y)\mathrm{d}y\right]}^2,$$

$$E_{\mathrm{eff}}=-\frac{c^2}{2}-\frac{c}{2}\int V(y)\mathrm{d}y+O(c^4),$$

$$V_{ZR}(y_1-y_2)=-\kappa \delta (y_1-y_2),$$

$$W_{ZR}(y_1-y_2)=\left[2c-2c^2/{\pi }^2\right]\delta (y_1-y_2)$$

$$E_{\mathrm{eff}}^{ZR}=-\frac{m_{\mathrm{eff}}}{2m}(c-c^2/{\pi }^2{)}^2$$

$$W_{BO}=2c\delta (y_1-y_2)+E_{BO}(y_1-y_2).$$

$$h_{\mathrm{eff}}=-\frac{m}{M}\frac{{\partial }^2}{\partial y_1^2}-\frac{m}{M}\frac{{\partial }^2}{\partial y_2^2}+2c_{II}\delta (y_1-y_2)+V(y_1-y_2),$$

$$c_{II}^{cr}\simeq \frac{c^2}{{\pi }^2}.$$

$$e^{i{k}_{j}l}=\frac{k_j-{\Lambda }_1+\frac{ic}{2}}{k_j-{\Lambda }_1-\frac{ic}{2}}\frac{k_j-{\Lambda }_2+\frac{ic}{2}}{k_j-{\Lambda }_2-\frac{ic}{2}},1\le j\le N;$$

$$\prod _{j=1}^N\frac{k_j-{\Lambda }_1+\frac{ic}{2}}{k_j-{\Lambda }_1-\frac{ic}{2}}=1,\prod _{j=1}^N\frac{k_j-{\Lambda }_2+\frac{ic}{2}}{k_j-{\Lambda }_2-\frac{ic}{2}}=1,$$

$$k_j=\frac{2\pi }{l}{m}_j+{\delta }_j,m_j\in \mathbb{Z}\setminus \{0\};$$

$$e^{i{k}_{j}l}=\underset{=1}{\underbrace{e^{i{m}_{j}2\pi }}}{e}^{i{\delta }_{j}l}\approx 1+i{\delta }_{j}l.$$

$$\begin{array}{r}\frac{k_j-{\Lambda }_1+\frac{ic}{2}}{k_j-{\Lambda }_1-\frac{ic}{2}}\frac{k_j-{\Lambda }_2+\frac{ic}{2}}{k_j-{\Lambda }_2-\frac{ic}{2}}\approx 1+\frac{2ic}{k_j^{(0)}},\end{array}$$

$$\begin{array}{rl}\frac{a+b}{a-b}& \approx 1+\frac{2b}{a}.\end{array}$$

$${\delta }_{j}l=\frac{2c}{k_j^{(0)}}\Rightarrow {\delta }_j=\frac{c}{\pi m_j}.$$

$$e^{-i{k}_{1}l}$$

$$\to e^{i{k}_{1}l}$$

$$e^{i{k}_{1}l}=\frac{k_1-{\Lambda }_1+\frac{ic}{2}}{k_1-{\Lambda }_1-\frac{ic}{2}}\frac{k_1-{\Lambda }_2+\frac{ic}{2}}{k_1-{\Lambda }_2-\frac{ic}{2}},$$

Full text

I Introduction

An environment with mobile impurity atoms is a cherished model system in quantum physics. It is a testbed for introducing and testing quasiparticle concepts, e.g., polarons and bipolarons, which naturally appear when studying the movement of electrons in crystals Landau and Pekar (1948); Pekar. (1963); Alexandrov and Mott (1995), 3He atoms in superfluid 4He Baym and Pethick (2008), or even protons in neutron matter Kutschera and Wójcik (1993). Nowadays, these concepts can be examined using quantum simulators based upon ultracold atoms Schirotzek et al. (2009); Nascimbène et al. (2009); Massignan et al. (2014); Hu et al. (2016); Jørgensen et al. (2016); Schmidt et al. (2018). An important topic that can be addressed with cold atom systems is the physics of impurity-impurity correlations induced by a medium Bruderer et al. (2007); Schecter and Kamenev (2014); Keiler et al. (2018); Naidon (2018); Dehkharghani et al. (2018); Camacho-Guardian and Bruun (2018); Camacho-Guardian et al. (2018); Pavlov et al. (2018); Mistakidis et al. (2019a, b); Reichert et al. (2019); DeSalvo et al. (2019). This topic is relevant for basic research, and in applications motivated by bound states of dressed electrons and their relation to high-Tc superconductors Alexandrov and Mott (1994),111 Experiments with ultracold atoms can give important insight into the physics of induced attractive potentials between dressed electrons. However, atoms interact via short-range potentials limiting the relationship between impurities in cold atoms and dressed electrons. For example, the predicted abrupt change of the mean distance between two polarons across the unbound-polarons to bipolaron transition Adamowski (1989) most probably cannot be simulated in these experiments because the Coulomb interaction is an important ingredient for observing this effect Schmickler et al. (2019)..

In this paper, we calculate the ground-state energy of a one-dimensional (1D) Fermi gas with two bosonic impurities, see Fig. 1. One-dimensional geometries typically enhance interaction effects Takada (1982) opening up the possibility of observing bound states supported by the induced attraction 222Note that any attractive potential supports at least one bound state in the 1D world Landau and Lifschitz (1977). Another feature that separates one spatial dimension from higher dimensions is the long-range tail of correlations. For example, Friedel oscillations Friedel (1958) decay as $\sim 1/r^{D}$ where $D$ is the dimension of space. These enhanced correlations may be useful to simulate many-body phenomena beyond short-range physics typical for cold atoms.

Our paper is organized as follows. We start by introducing the Hamiltonian of our 1D model in Section II. In Section III, we study impurity-impurity correlations in the limiting case of equal masses, $M=m$. All interactions are identical and parametrized by Dirac delta functions. The fermions do not interact among each each other due to the Pauli exclusion principle. The system is solvable by the Bethe ansatz, which is a common starting point for analyzing cold atom systems in 1D geometries Cazalilla et al. (2011); Guan et al. (2013); Batchelor and Foerster (2016). In Section IV, we go on to discuss effective models for describing two impurities in a medium and benchmark them against the Bethe ansatz results. Afterwards, we use the effective models to investigate non-integrable systems; our focus is on the appearance of in-medium bound states. We discuss the transition from unbound impurities to bound impurities, which can be tested in cold atom experiments, for example, by measuring the spectroscopic response or by studying the collapse dynamics DeSalvo et al. (2019) in imbalanced Bose-Fermi mixtures. Finally, Section V contains a summary of our results and an outlook.

II Formulation

We consider two bosonic impurities interacting via a zero-range potential with $N_{f}$ spinless (fully spin-polarized) fermions. For convenience, we assume that $N_{f}$ is an odd number. This assumption does not limit the generality of our results as we are interested in the limit $N_{f}\to\infty$. Note that the bosons may, in principle, possess spin quantum numbers. However, the bosonic spin part of the wave function is not important for our discussion, since we focus on the ground state. It is only necessary that the wave function is symmetric with respect to the exchange of the spatial coordinates of the impurity atoms. The particles are confined to a ring of radius $R$, see Fig. 1. The Hamiltonian for the system reads

[TABLE]

where $H_{f}$ describes fermions:

[TABLE]

with $m$ the fermion mass. $H_{b}$ describes the impurity bosons:

[TABLE]

where $M$ is the boson mass, and $g_{II}$ is the strength of the boson-boson interaction. The interaction between fermions and bosons is written as

[TABLE]

where $g$ is the corresponding interaction strength. We solve the Schrödinger equation $H\psi=\epsilon\psi$ for the ground state for different $N_{f}$ and sizes $L=2\pi R$ of the system. Then, we extrapolate the energies to the thermodynamic limit, $N_{f},L\rightarrow\infty$, assuming a fixed density of the Fermi gas, $N_{f}/L=\rho$. For convenience, we introduce the dimensionless quantities $y_{j}=Y_{j}\rho$, $x_{j}=X_{j}\rho$, $l=L\rho$, $c_{II}=mg_{II}/(\hbar^{2}\rho)$, $c=mg/(\hbar^{2}\rho)$, and $\varepsilon=2m\epsilon/(\hbar^{2}\rho^{2})$.

Note that working with a finite number of particles allows us also to study the transition of the energy from a few- to many-body limits, an exciting research venue, which can be addressed with cold atoms Wenz et al. (2013); Zinner (2016); Rammelmüller et al. (2017a). We observe that the ground state energies of systems with $N_{f}$ of the order of 10-20 particles can be accurately described using the results derived in the thermodynamic limit.

III Solvable Limits

III.1 Bethe-Ansatz-Solvable Case

First, we consider the most symmetric case: $c_{II}=c$ and $m=M$, whose Hamiltonian we write as

[TABLE]

where $N=N_{f}+2$; we set $x_{N-1}=y_{1}$ and $x_{N}=y_{2}$ to explicitly demonstrate the particle exchange symmetry. The ground state of $h_{BA}$ with fermions at the coordinates ($x_{1},...,x_{N_{f}}$) and bosons at ($x_{N-1},x_{N}$) can be studied experimentally with SU(3)-symmetric fermions, e.g., with 173Yb Pagano et al. (2014). Indeed, the ground state of SU(3) fermions with the particle decomposition $N_{f}+1+1$ has a bosonic symmetry for the exchange of the particles in the $1+1$ subsystem. To understand this, note that: (i) the particles in the $1+1$ subsystem are distinguishable particles and, hence, there exist no apriori symmetry requirements for their exchange; (ii) the Hamiltonian $h_{BA}$ commutes with the particle exchange operator; (iii) the bosonic symmetry leads to the lowest energy. Furthermore, we expect that the ground state energies of isotopic systems with a small mass imbalance, e.g., 6Li-7Li and 39K-40K (cf. Ferrier-Barbut et al. (2014); Wu et al. (2011)), can be accurately described by Eq. (5).

The spectrum of the Hamiltonian $h_{BA}$ can be found using the Bethe ansatz (BA) Yang (1967). Let us briefly summarize this approach. For every ordering of particles (e.g., for $x_{1}<x_{2}<...<x_{N}$), the wave function is written as a sum of the plane waves $e^{i\sum_{j}k_{j}x_{j}}$. For this wave function to fulfill the boundary conditions at $x_{i}=0$, $x_{i}=l$ and $x_{i}=x_{j}$ for all $i$ and $j$, the quasi-momenta $k_{j}$ must satisfy the BA equations

[TABLE]

where the bosonic and fermionic symmetries have already been implemented Lai and Yang (1971); Batchelor et al. (2005); Imambekov and Demler (2006). $\Lambda_{1}$ and $\Lambda_{2}$ are to be determined together with the set $\{k_{j}\}$. Once the BA equations are solved, the energy of the system is determined as $\varepsilon=\sum_{j=1}^{N}k_{j}^{2}$. Note that the number of unknowns in Eqs. (6) for the ground state can be reduced to $(N_{f}+3)/2$ from $N_{f}+4$. Indeed, the total (angular) momentum must be zero in the ground state, $\sum_{j}k_{j}=0$. This together with the fact that the wave function is real makes the quasi-momenta appear in pairs (we exemplify this below). In addition, one can show that $\Lambda_{1}=-\Lambda_{2}$.

To solve Eq. (6) for the ground state, we apply Newton’s method, which requires an accurate initial estimate of $k_{j}$ and $\Lambda_{j}$. For $c\to 0$, we obtain this estimate directly from the BA equations (see Appendix A):

[TABLE]

where $k_{j}^{(0)}$ is the quasi-momenta at $c=0$. Note that in Eq. (7) the quasi-momenta are related pairwise, e.g., $k_{1}=-k_{2}$, as has already been mentioned. The only non-paired quasi-momentum is $k_{3}=0$. Estimate (7) allows us to calculate $\{k_{i}\}$ and obtain solutions for weak interactions, $|c|\ll 1$. We then follow the solutions as $c$ is changed in small steps. An initial guess for moderate interactions is obtained from a Taylor series constructed using solutions at smaller values of $|c|$, see Appendix B.

We solve Eq. (6) for a sequence of $N_{f}$ and extrapolate to the thermodynamic limit. To this end, we subtract from the energy the zero-interaction offset and fit the difference with $\varepsilon(c)-\varepsilon(0)=\varepsilon_{\infty}+A_{1}/N+A_{2}/N^{2}$, where $\varepsilon_{\infty},A_{1}$ and $A_{2}$ are fitting parameters. It is straightforward to argue for the form of the fitting function, $\varepsilon_{\infty}+A_{1}/N+A_{2}/N^{2}$, in the case of strong interactions ($c\to\pm\infty$) for which the energies can be calculated using a non-interacting Fermi gas. We do not attempt to validate the fitting function for finite values of $c$, since we observe that the form of the function is not important for our analysis (see Appendix C).

To investigate induced correlations in the thermodynamic limit, we introduce the “in-medium binding energy” $E=\varepsilon_{\infty}-2\mathcal{E}$, where $\mathcal{E}$ is the energy gain for immersing one impurity in a Fermi gas McGuire (1965, 1966) (see Appendix D). The quantity $2\mathcal{E}$ describes the energy of two non-correlated impurities. $E$ is presented in Fig. 2. Note that $E$ cannot be positive, since any induced correlations between impurities must vanish when they are far apart. If $E=0$ the two impurities are not correlated, in general, they are infinitely far from each other. In this case, we say that there is no in-medium bound state, whereas if $E<0$ then there is at least one. Next, we analyze cases with $c<0$ and $c>0$ separately.

Repulsive interactions, $c>0$. We calculate the energy for $c\lesssim 2$ and find that $E=0$ (within numerical accuracy), which means that there are no in-medium bound states. In principle, an in-medium bound state is possible for two impurities that repel each other without the Fermi gas. This can happen if the energy cost of embedding an impurity in the medium is larger than the cost of bringing the impurities close together. Our calculations show that an in-medium bound state does not appear when $c_{II}=c>0$. For $c\to 0$ it can easily be understood that $E=0$, since the impurity-impurity interaction without the Fermi gas scales as $c$ (see Eqs. (3) and (5)), whereas the induced impurity-impurity interaction is expected to scale as $c^{2}$ (see Sec. IV). Therefore, the interaction volume, i.e., the space integral of the effective impurity-impurity interaction, is necessarily a small positive quantity for $c\to 0$, which does not allow for the existence of a bound state Landau and Lifschitz (1977); Simon (1976).

In the limit of strong interactions some extra insight can also be gained. For $c\to\infty$ the important degrees of freedom are spins Volosniev et al. (2014); Deuretzbacher et al. (2014); Volosniev et al. (2015); Massignan et al. (2015), which allows one to map the Hamiltonian $h_{BA}$ onto an XX spin chain Deuretzbacher et al. (2017) with constant coefficients,

[TABLE]

where $\sigma_{x}^{i}$ and $\sigma_{y}^{i}$ are the Pauli matrices acting on the spin at site $i$; $J$ is an exchange coefficient proportional to $1/c$, see Volosniev (2017) for the derivation in a homogeneous environment. The system in Fig. 1 with $c\to\infty$ is then identical to a linear spin chain with two spin deviates (magnons) for which a bound state is not expected Hodgson and Parkinson (1984).

Attractive interactions, $c<0$. Figure 2 shows that for $c<0$ there is an in-medium bound state whose energy is below the ground state energy of the Hamiltonian that describes two bosons without the Fermi gas, i.e., $H_{b}$ from Eq. (3). This lowering of the energy is a manifestation of the induced impurity-impurity correlations, which we interpret in Sec. IV using an effective impurity-impurity potential. Since we are interested in interactions mediated by the Fermi gas, we do not discuss in this paper the limit $c\to-\infty$, which implies the formation of a tightly bound trimer (in analogy to trions in SU(3) symmetric Fermi gases Guan et al. (2013)). This trimer is supported by the fundamental (not induced) interaction, and, hence, is beyond the scope of this paper. For the sake of discussion, we present the energy of a trimer in Fig. 2. We expect that this energy, $-2c^{2}$, determines the behavior of the system as $c\to-\infty$. In the rest of the paper, we only explore $c\gtrsim-2$.

III.2 Two Static Impurities

Before we discuss effective models that describe two mobile impurities in a sea of fermions, we consider two static impurities $M\to\infty$ – another analytically solvable limit of the Hamiltonian (1). The solution is obtained by solving the one-body problem: one particle in a ring with a potential due to the two impurities fixed at $-r/2$ and $r/2$

[TABLE]

where the subscript emphasizes the connection to the Born-Oppenheimer (BO) approximation, which will be employed below. The spectrum of $h_{BO}$ depends on the distance $r$. We calculate this dependence only for attractive interactions, i.e., $c<0$; the repulsive case can be calculated in a similar manner. To obtain the ground state energy of the Fermi gas, $\varepsilon_{BO}(c,r)$, we add the energies of the lowest $N_{f}$ eigenstates of the Hamiltonian $h_{BO}$ (see Appendix E). The thermodynamic limit is calculated by extrapolating the results for systems with different values of $N_{f}$ and $l$ and a fixed ratio $N_{f}/l$. We observe that already for $N_{f}\simeq 19$ the energy for small values of $r$ can be used (for our discussion) as the energy in the thermodynamic limit. The solution for $r\to\infty$ can be obtained by fitting to the known form of the tail Recati et al. (2005); Fuchs et al. (2007).

Figure 3 illustrates the energy $E_{BO}(r)=\varepsilon_{BO}(c,r)-\varepsilon_{BO}(0,r)-2\mathcal{E}_{static}$ for $c=-0.5$ (we assume that $g_{II}=0$ for the sake of discussion). $\mathcal{E}_{static}$ is the energy gain for immersing a single static impurity in a Fermi gas Parisi and Giorgini (2017)

[TABLE]

The quantity $E_{BO}$ has a deep minimum at $r=0$ given by $\mathcal{E}_{static}(2c)-2\mathcal{E}_{static}(c)$ and an oscillatory tail. To derive this limiting value, note that when both impurities are at $r=0$, they act as a single impurity with the strength $2c$. For $c\to 0$ the tail can be written simply as $c^{2}\cos(Ar+\delta)/r$, where $A$ and $\delta$ are constants. This tail can be obtained from Friedel oscillations Friedel (1958).333An alternative derivation is given in Refs. Recati et al. (2005); Fuchs et al. (2007). These oscillations determine the density of the Fermi gas at the position of the second impurity, provided that the first impurity is separated by the distance $r$. This density in turn determines the energy of the system, according to first order perturbation theory. It is worthwhile noting that the dependence of $E_{BO}$ on $r$ is observable. It can, in principle, be probed in cold atom experiments by spectroscopy Recati et al. (2005).

IV Effective model

IV.1 Bethe-ansatz-solvable system

The ground state of a one-dimensional Fermi gas ($N_{f}\to\infty$) with a single interacting impurity is orthogonal (for all non-zero values of interaction strengths) to the ground state of the corresponding non-interacting system. This phenomenon is related to the Anderson orthogonality catastrophe Anderson (1967). For the SU(2) symmetric case it can be conveniently studied using the BA equations Castella and Zotos (1993). This orthogonality reduces the applicability of the polaron picture. For example, the dynamics after a sudden change of parameters cannot be captured by a single quasiparticle, instead it requires a continuum of states. Still, the notion of the effective mass and self-energy can be used to describe the low-energy spectrum of a Fermi gas with one impurity McGuire (1965, 1966); Castella and Zotos (1993); Giraud and Combescot (2009); Mao et al. (2016), suggesting the use of the following Hamiltonian to model the binding energy for the system of two impurities

[TABLE]

where $m_{\mathrm{eff}}(c)$ is the effective mass of the impurity whose analytical form is presented in Refs. McGuire (1966); Castella and Zotos (1993); Giraud and Combescot (2009), so that the Hamiltonian $h_{\mathrm{eff}}$ with $W=0$ correctly describes the low-energy spectrum of two non-interacting impurities. For simplicity, we work with the expansion of the effective mass truncated at the order $O(c^{4})$,

[TABLE]

The terms of the order $O(c^{5})$ can be neglected for interactions considered here, $c\in(-2,0)$. This equation shows that the mass does not increase by more than a few percent for the considered parameters.

The function $W$ in Eq. (11) describes the impurity-impurity correlations. We choose it as

[TABLE]

where the first term is the interaction between impurities without the Fermi gas; the second term is an effective interaction mediated by the environment. Note that the exact shape of $V$ is not required for our purpose. We are interested in the weak and moderate interaction regimes for which the knowledge of a few integrated properties of $V$ is sufficient, e.g., only $\int V(y)\mathrm{d}y$ determines the binding energy for $c\to 0$. Indeed, the ground state energy $E_{\mathrm{eff}}$ of $h_{\mathrm{eff}}$ for weak interactions reads Landau and Lifschitz (1977); Simon (1976)

[TABLE]

where the integration is over the whole system in the thermodynamic limit. This expression can be expanded as

[TABLE]

where we assume that $\int V\mathrm{d}y$ scales as $c^{2}$ (see below) to estimate the neglected pieces. Note that the renormalization of mass enters in $O(c^{4})$ meaning that this effect may be disregarded for $c\to 0$ when compared to the effect of the two-body effective interaction. The first term in Eq. (15) is the ground state energy of two particles without the Fermi gas. The integral $\int V\mathrm{d}y$ must be negative to ensure that the energy of the in-medium bound state is below $-c^{2}/2$ as in Fig. 2. Therefore, the effective interaction must be overall attractive. Let us discuss the two (arguably) simplest approximations that can be used to appropriately choose $V$.

Zero-range Potential. The most basic form of $V$ in Eq. (13) consistent with the first two terms of the expansion (15) is the zero-range (ZR) potential

[TABLE]

where $\kappa\equiv|\int V\mathrm{d}y|$. This potential can be used to reproduce low-energy properties of two impurities when higher order terms in Eq. (15) are not important. The parameter $\kappa$ can be obtained, for example, from a single-phonon exchange Schecter and Kamenev (2014), in which case $\kappa=2c^{2}/\pi^{2}\simeq 0.202c^{2}$ for $c\to 0$. If the potential

[TABLE]

is used in Eq. (11) then a single bound state with the energy

[TABLE]

is produced. This effective model captures qualitatively the exact results, see Fig. 4 (top). We show the ground state energies of $h_{\mathrm{eff}}$ with $m_{\mathrm{eff}}$ from Eq. (12) and with $m_{\mathrm{eff}}=m$, to illustrate that the mass renormalization leads only to a marginal correction for the considered values of $c$.

Induced interaction from the Born-Openheimer approximation. The potential $V$ in Eq. (13) can also be derived in the Born-Oppenheimer approximation, where it is assumed that the Fermi gas follows the impurity adiabatically. The potential in this case is simply the energy in Fig. 3. This approximation must be accurate if the impurities are very heavy. For mobile impurities this approximation is useful if the impurities move slowly in comparison to the Fermi velocity, which defines the dispersion of a sound mode in a Fermi gas.

The function $E_{BO}(y_{1}-y_{2})$ decays as $1/|y_{1}-y_{2}|$, however, it leads to an effectively short-range potential due to the oscillatory tail. For example, $\int E_{BO}(y)\mathrm{d}y$ is well-defined. We calculate that $|\int E_{BO}(y)\mathrm{d}y|\simeq 0.22c^{2}$ for $c\to 0$. This is slightly larger than $\kappa\simeq 0.202c^{2}$ for the zero-range potential discussed above. Even though, the long-range tail is not expected for integrable systems Schecter and Kamenev (2014), the potential $E_{BO}$ performs as well as the zero-range potential, confirming that only integral properties of $V$ are important. Figure 4 (center) gives the binding energy calculated using the potential

[TABLE]

Figure 4 (bottom) shows the quantity $E+c^{2}/2$ to single out the effect of the induced interaction. The center and bottom panels of Fig. 4 demonstrate that the $E_{BO}$ can be used to qualitatively analyze in-medium bound states. We note that the difference between the filled and empty circles in Fig. 4 (bottom) for $|c|<0.5$ can be fit with $-0.036c^{2}-0.1873c^{3}$, which demonstrates the importance of terms with $c^{n}$, $n>2$. To reduce this disagreement between the exact results and the effective model one can include coulpings between eigenstates of $h_{BO}$ due to the motion of the impurities. We leave this discussion for future studies.

To calculate the data in Fig. 4 (center, bottom), we exactly diagonalize the Hamiltonian $h_{\mathrm{eff}}$ from Eq. (11) with $W_{BO}$: We use the eigenstates of the effective Hamiltonian with $W_{BO}=0$ as a basis to write $h_{\mathrm{eff}}$ as a matrix and diagonalize this matrix after truncation. The energy is found by fitting to the form $E_{K}=E_{\mathrm{eff}}+D/K$, where $K$ is the number of used basis states and $D,E_{\mathrm{eff}}$ are fitting parameters. This slow $1/K$-convergence is expected for zero-range interactions Volosniev and Hammer (2017).

Other potentials. The effective potential can also be calculated using other approximation schemes. For example, trial wave functions Adamowski (1989); Vansant et al. (1994); Dehkharghani et al. (2018) can be used. We do not discuss those approaches here. However, we note that different methods may lead to different shapes of the effective potential. This is not surprising, since the effective potential is not an observable quantity for mobile impurities.

IV.2 Non-integrable cases

Motivated by the accuracy of the effective model (11) for the most symmetric case, we extend this model to study appearance of in-medium bound states in non-integrable systems, i.e., $c_{II}\neq c$ and/or $m\neq M$. We write the corresponding effective Hamiltonian as

[TABLE]

where we use $m_{\mathrm{eff}}=M$ for simplicity. This approximation relies on the observation that the mass renormalization is not important for the Bethe-ansatz-integrable case for weak interactions. $h_{\mathrm{eff}}$ supports a bound state for all $c_{II}<0$ because the induced interaction is attractive. For $c_{II}>0$ the bound states appear only if $c_{II}<c_{II}^{cr}$ for which the repulsive impurity-impurity interaction in the absence of fermions is overtaken by the attractive interaction mediated by the Fermi gas.

We first consider the case when $c_{II}\neq c$ and $m=M$, in which the particle exchange symmetry is broken by the interaction term 444For strong interactions, this case can still be studied with integrable Heisenberg Hamiltonians using the mapping onto a spin chain Volosniev et al. (2014); Deuretzbacher et al. (2014); Volosniev et al. (2015); Massignan et al. (2015). We do not discuss this limit here, instead we focus on the regimes with weak and moderate interactions. The kinetic energy is still symmetric with respect to the exchange of two particles. We calculate binding energies using the zero-range effective potential and the potential from the adiabatic approximation, see Fig. 5 (top). Both potentials lead to similar results for $c_{II}=0$ and small values of $c$. Note that in this regime the one-phonon exchange potential gives a leading contribution in $c$. For larger values of $c$ the ZR potential must be corrected.

For $c_{II}=0.2$ the ZR potential predicts that the in-medium bound state is formed at smaller values of $c$ in comparison to the BO potential, see Fig. 5 (top). The difference between $c_{II}^{cr}$ in the two methods is, however, marginal and one can use the ZR potential to derive the critical value for the appearance of the bound state:

[TABLE]

The transition from an overall repulsive to an overall attractive induced interaction can be studied in Bose-Fermi mixtures by looking, for example, at the collapse dynamics DeSalvo et al. (2019): a Bose gas can be stable only if its particles repel each other. One could also study spectroscopically the energy needed to break a bound state by transfering the system into a non-interacting state. Let us estimate the in-medium binding energy for $m=M$. For $c=-2$ and $c_{II}=0$, we have $E\simeq-0.06$ (see Fig. 5 (center)), which for 6Li atoms with $\rho=3/(\mu m)$ translates into $\simeq 22$ nK $\times k_{B}$, where $k_{B}$ is the Boltzmann constant. This means that in-medium bound states can be observed only at ultracold temperatures.

The in-medium binding energy can be increased if the impurities are heavy, cf. Eq. (14) with $m_{\mathrm{eff}}=M$. To explicitly show this, we consider $c_{II}\neq c$ and $m\neq M$. In this case both the interaction and kinetic energies are not symmetric with respect to the exchange of two particles. For the sake of discussion, we first use the BO potential to calculate the energies of the system with $M=2m$, see Fig. 5 (center). The behavior of the energy resembles that for $m=M$, but, as expected, the overall energy scale is now larger. It is worthwhile noting that the critical value $c^{cr}_{II}$ obtained using the ZR approximation does not depend on the mass $M$. The dependence of $c^{cr}_{II}$ on $M$ is also negligible for the BO potential. Therefore, Eq. (21) can be used to predict the appearance of bound states for different masses of the impurity.

Finally, we consider two bosonic 133Cs atoms in a fermionic gas of 6Li as in the experiment of Ref. DeSalvo et al. (2019), see Fig. 5 (bottom). We use both the Born-Oppenheimer potential as well as the zero-range potential. Note however that for $M\gg m$ the former is expected to perform better than the latter. For a Li-Cs mixture the energy scale is larger than that for lighter impurities and the bound states should be observable at much higher temperatures.

V Summary/Outlook

We investigate the problem of two bosonic impurities in a spin-polarized Fermi gas. First, we consider the ground state energy of the system in the Bethe-ansatz-solvable case, i.e., equal masses of fermions and impurities, $m=M$, as well as equal impurity-impurity and impurity-fermion interactions, $c_{II}=c$. We calculate the ground state energy and show that there are attractive impurity-impurity interactions induced by the fermionic medium. In the next step, we discuss an effective model for the induced interactions and compare its predictions to the exact results. We use two effective potentials to define the effective Hamiltonian for the impurity system: a zero-range potential matched to single-phonon exchange and an adiabatic potential for heavy impurities derived in the Born-Oppenheimer approximation. Both potentials are able to approximate the exact results from the Bethe ansatz. The difference between the two model potentials for $m=M$ allows us to estimate the errors and the breakdown of our effective Hamiltonian. For the Bethe-ansatz-integrable case the difference between the results derived using the two effective potentials is marginal, which argues in favor of using them for qualitative analysis of Fermi gases with impurities.

The success of the effective model in the integrable case motivates our use of the effective model to study non-integrable systems characterized by relaxing at least one of the two conditions $m=M$ and $c_{II}=c$. For repulsive impurity-impurity interactions without the Fermi gas, $c_{II}>0$, we predict that the induced interaction overcomes the repulsion if the impurity-fermion interaction satisfies $c<-\pi\sqrt{c_{II}}$, leading to an in-medium bound state. The binding energies are larger for heavier impurities such that the observation of in-medium bound states in heavy-light mixtures appears more promising.

Our findings show that the Bethe-ansatz-solvable models provide a playground for investigating induced interactions. In the future it will be interesting to use the Bethe ansatz equations (6) to investigate spatial correlations of two impurities, which will allow us to test an effective model beyond the simple energy comparison presented here. Further studies of non-integrable systems are also needed. The non-integrability due to $c_{II}\neq c$ and $m\neq M$ has been briefly discussed here. For cold atoms it is important also to consider trap effects, which break the integrability and change the properties of the system Zinner (2016); Sowiński and García-March .

It will be interesting to extend the present study to fermionic impurities. It is known that two fermionic impurities in the SU(2) case do not have an in-medium bound state Flicker and Lieb (1967). However, if the impurities are very heavy then a bound state must exist: The BO potential in Fig. 3 has attractive regions and unlike the zero-range potential of Eq. (16) has a finite range. Since it becomes exact as $M\to\infty$, the system must support a bound state in this limit. This prediction can be explored using numerical many-body methods that are able to deal with mass-imbalanced systems such as the complex Langevin approach Rammelmüller et al. (2017b, 2018).

Finally, we note that there is still a limited number of exact numerical results on two impurities in Bose gases. The present work gives insight into the properties of impurities in strongly-interacting Bose gases Girardeau (1960). However, further work should be done to understand impurities in weakly-interacting Bose gases. Since these systems cannot be studied using the Bethe ansatz, one has to employ other methods Grusdt et al. (2017); Pastukhov (2017); Dehkharghani et al. (2018); Volosniev and Hammer (2017); Parisi and Giorgini (2017); Kain and Ling (2018); Mistakidis et al. (2019c); Paeckel et al. (2019); Schollwöck (2005); Gull et al. (2011); Pasek and Orso (2019).

Acknowledgements.

We thank Lukas Rammelmüller and Simos Mistakidis for useful discussions and comments on the manuscript. This work has been supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project numbers 413495248 – VO 2437/1-1 and 279384907 – SFB 1245, by the Bundesministerium für Bildung und Forschung (BMBF) through contract 05P18RDFN1, and by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 754411.

Appendix A Weak coupling expansion

In this appendix, we derive a weak coupling expansion of the BA equations,

[TABLE]

where $N$ is an odd number. First, we consider the quasi-momenta $k_{j}$ ($j=3,...,N$) that satisfy $k_{j}(c=0)\neq 0$. For $c=0$, these quasi-momenta are multiples of $2\pi/l$. When $c\neq 0$, we write them as

[TABLE]

where $\delta_{j}$ is small. Inserting Eq. (24) into the left-hand-side of Eq. (22) leads to

[TABLE]

We write the right-hand-side of Eq. (22) as

[TABLE]

where terms proportional to $c^{n}$ with $n>1$ are neglected. Also it is assumed that $k_{j}^{(0)}\gg\delta_{j}-\Lambda_{1}$ and $k_{j}^{(0)}\gg\delta_{j}-\Lambda_{2}$. This assumption is valid, since the $\Lambda$’s lie in between the first three quasi-momenta, which are all close to zero. To derive Eq. (26), we use that for $a\gg b$

[TABLE]

With Eqs. (25) and (26), we obtain

[TABLE]

Now we investigate the quasi-momenta $k_{j}$ that vanish at $c=0$. As discussed in the main text, for the ground state, $k_{1}=-k_{2}$ and $k_{3}=0$. To show that $\Lambda_{1}=-\Lambda_{2}$, we consider Eq. (22) for $k_{1}$ and $k_{2}$. We use $k_{1}=-k_{2}$ in the equation for $k_{2}$:

[TABLE]

From this equation and the equation for $k_{1}$:

[TABLE]

we obtain that $\Lambda_{1}=-\Lambda_{2}$. Using the equation for $k_{1}$, we derive

[TABLE]

Next, we consider Eq. (23). The sum of the quasi-momenta $\{k_{3},...,k_{N}\}$ is zero, thus, Eq. (23) can be rewritten as

[TABLE]

With $k_{1}=-k_{2}$, $k_{3}=0$, this equation reads

[TABLE]

We rewrite Eq. (32)

[TABLE]

which leads to $k_{1}=\sqrt{3c/l}$ and $\Lambda_{1}=\sqrt{c/l}$.

Appendix B Numerical method to solve the BA equations

In this appendix we give a detailed explanation of the numerical method that we employ to solve the BA equations. For simplicity, we consider here the case of 5 particles, three fermions and two bosons. To solve the BA equations, we employ Newton’s method. Initial conditions for which are established differently for different regions of $c$. We exemplify different regions in Figs. 6 and 7, which show the calculated values of $k_{i}$. The quasi-momenta $k_{1}$ and $k_{2}$ are purely imaginary, they are shown in Fig. 6; the quasi-momenta $k_{3},\leavevmode\nobreak\ k_{4}$ and $k_{5}$ are purely real, they are presented in Fig. 7. For better visibility, we do not show every value of $k_{i}$ that we calculate, only every fourth point. The three numbered areas with different types of shading refer to different methods we use to establish the initial conditions for Newton’s method.

For weak interactions (area 1) we use the weak coupling expansion from Eq. (7), as an initial guess for Newton’s method. As can be seen in Fig. 6, for $c\lesssim-0.5$ the deviation between Eq. (7) and the exact solution becomes significant. We need another approach to calculate $k_{i}$ for stronger interactions.

In the second area, we use the solution at the previous step as an initial guess for Newton’s method (see the green arrows in Fig. 7). At about $c\approx-1.5$, this method requires us to compute points lying very close to each other, which can be time consuming. So once more, we change our strategy.

To come up with an initial guess for Newton’s method in the third region, we extrapolate the previously calculated solutions by using a polynomial fit with order 3. The red crosses in Fig. 7) exemplify a fitting range used to calculate the fitting function (dotted,blue curve). At $c\approx-2$ the fit function does not represent the exact solution well. Therefore, the fitting process must be repeated using the calculated solutions in the range $c\in(-2,-1.5)$ (not shown here for brevity).

Appendix C Thermodynamic extrapolation

To extrapolate the calculated energies $\varepsilon(c)-\varepsilon(0)$ to the thermodynamic limit, we shall employ the two fit functions:

[TABLE]

To illustrate the fits, we show in Figs. 8 and 9 the exact energies as functions of $N$ for two different interaction strengths together with the corresponding fits. Both functions (36) and (37) appear to represent the data well. They also produce similar results for $N\to\infty$. The values of $\varepsilon_{\infty}$ from the two fits differ only in the third digit, implying that the precise knowledge of the convergence pattern as $N\to\infty$ is not needed for the considered parameters.

Appendix D One impurity

Here we briefly review how to derive the ground state energy of a Fermi gas with a single impurity atom. This system has already been investigated McGuire (1965). We use this well-known set-up to test our numerical approach. The system is desribed the Hamiltonian

[TABLE]

where the coordinates $x_{1},...,x_{N-1}$ are the positions of the fermions, and $x_{N}$ is reserved for the impurity. The corresponding Bethe ansatz equations are given by Yang (1967)

[TABLE]

where $k^{(1)}_{j}$ is the $j$th quasi-momentum (we use the superscript to emphasize that we work with a single impurity here), and $\Lambda$ is one additional variable. We consider $N$ to be even. Once the BA equation are solved, the energy can be calculated as $\varepsilon^{(1)}=\sum_{j}\left(k^{(1)}_{j}\right)^{2}$.

We solve the BA equations with Newton’s method as already explained in the main part. For small $c$ the weak coupling expansion of the BA equations Oelkers et al. (2006) is used as an initial guess

[TABLE]

where $m_{j}\in\mathbb{Z}$ determine the quasi-momenta of the particles for zero interaction. The shift due to the small interaction is given by the terms proportional to $c$ and $\sqrt{c}$.

To extrapolate the result to the thermodynamic limit we use $\varepsilon^{(1)}(c)-\varepsilon^{(1)}(0)=\mathcal{E}+\alpha/N^{\beta}$, where $\mathcal{E},\alpha$ and $\beta$ are fit parameters. We show the result in Fig. 10. Our result fits the analytic expression quite well. The relative difference, shown in the inset, is negligible for our purposes. We present also the result for 14 particles to demonstrate that only a handful of particles are needed to simulate the ground state properties of an infinite Fermi gas with an impurity in a laboratory.

Appendix E Two static impurities

Here we calculate the eigenenergies of the Hamiltonian

[TABLE]

To this end, we divide the space into three parts: $-l/2<x<-r/2$, $-r/2<x<r/2$, and $r/2<x<l/2$. For each part we write the wave function as

[TABLE]

where $k$ is the wave number. The wave function must obey the “delta-potential boundary conditions” at $x=\pm r/2$, and periodic boundary conditions at $x=\pm l$. We divide the solutions into parity-symmetric and parity-antisymmetric ones. The energy $k^{2}$ must be represented by a real number, which means that $k$ can be either real or imaginary. These two possibilities are referred to as the “scattering states” and the “bound states”, respectively. For “bound states”, we write $k=i\kappa$, where $\kappa$ is real. The equations that determine $k$ for each class of solutions are written below.

1. Symmetric “bound states”:

[TABLE]

2. Antisymmetric “bound states”:

[TABLE]

3. Symmetric “scattering states”:

[TABLE]

4. Antisymmetric “scattering states”:

[TABLE]

To solve the equations, a genetic algorithm first finds approximate solutions for $k,\kappa$. These are then used in the Newton’s iteration method as an initial guess. We calculate as many energy levels as particles we consider.