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

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.

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 [27] opening up the possibility of observing bound states supported by the induced attraction 2 . Another feature that separates one spatial * dhuber@theorie.ikp.physik.tu-darmstadt.de † hammer@theorie.ikp.physik.tu-darmstadt.de ‡ artem.volosniev@ist.ac.at 1 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 [25] most probably cannot be simulated in these experiments because the Coulomb interaction is an important ingredient for observing this effect [26]. 2 Note that any attractive potential supports at least one bound state in the 1D world [28] dimension from higher dimensions is the long-range tail of correlations. For example, Friedel oscillations [29] decay as ∼ 1/r D where D is the dimension of space. These enhanced correlations may be useful to simulate manybody 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 [30][31][32]. 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 [23] 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 spinpolarized) 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 arXiv:1908.02483v2 [cond-mat.quant-gas] 17 Dec 2019 bosons, M fermions, m R Figure 1. An illustration of the system: Two bosonic impurities in a one-dimensional Fermi gas. Periodic boundary conditions are employed, i.e., the system lives on a ring of radius R. The mass of an impurity (fermion) is denoted by M (m).
limit N f → ∞. 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 where H f describes fermions: with m the fermion mass. H b describes the impurity bosons: 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 where g is the corresponding interaction strength. We solve the Schrödinger equation Hψ = ψ for the ground state for different N f and sizes L = 2πR of the system. Then, we extrapolate the energies to the thermodynamic limit, N f , L → ∞, assuming a fixed density of the Fermi gas, N f /L = ρ. For convenience, we introduce the dimensionless quantities y j = Y j ρ, x j = X j ρ, l = Lρ, c II = mg II /( 2 ρ), c = mg/( 2 ρ), and ε = 2m /( 2 ρ 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 [33][34][35]. 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.

A. Bethe-Ansatz-Solvable Case
First, we consider the most symmetric case: c II = c and m = M , whose Hamiltonian we write as 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 173 Yb [36]. 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., 6 Li-7 Li and 39 K-40 K (cf. [37,38]), can be accurately described by Eq. (5).
The spectrum of the Hamiltonian h BA can be found using the Bethe ansatz (BA) [39]. 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 j kj xj . 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 where the bosonic and fermionic symmetries have already been implemented [40][41][42]. Λ 1 and Λ 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 ε = N j=1 k 2 j . 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, 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 Λ 1 = −Λ 2 .
To solve Eq. (6) for the ground state, we apply Newton's method, which requires an accurate initial estimate of k j and Λ j . For c → 0, we obtain this estimate directly from the BA equations (see Appendix A): where k (0) j 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| 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 ε(c) − ε(0) = ε ∞ + A 1 /N + A 2 /N 2 , where ε ∞ , A 1 and A 2 are fitting parameters. It is straightforward to argue for the form of the fitting function, ε ∞ + A 1 /N + A 2 /N 2 , in the case of strong interactions (c → ±∞) 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 = ε ∞ − 2E, where E is the energy gain for immersing one impurity in a Fermi gas [43,44] (see Appendix D). The quantity 2E 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 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 → 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  Figure 2. The (blue) dots show the in-medium binding energy of two bosonic impurities in a Fermi gas for the Bethe-ansatzintegrable case, the size of the dots can be used as an estimate for an error from the fitting procedure used to obtain the thermodynamic limit (see Appendix C). The solid (blue) curve is added to guide the eye. The dashed (orange) curve shows the binding energy of two bosons without the Fermi gas. We interpret the gap between the dashed curve and the dots as a manifestation of attractive impurity-impurity interactions mediated by the Fermi gas. For comparison, we also present the ground state energy of two bosons and one fermion [45], see the solid (black) curve.
interaction volume, i.e., the space integral of the effective impurity-impurity interaction, is necessarily a small positive quantity for c → 0, which does not allow for the existence of a bound state [28,46]. In the limit of strong interactions some extra insight can also be gained. For c → ∞ the important degrees of freedom are spins [47][48][49][50], which allows one to map the Hamiltonian h BA onto an XX spin chain [51] with constant coefficients, where σ i x and σ i y are the Pauli matrices acting on the spin at site i; J is an exchange coefficient proportional to 1/c, see [52] for the derivation in a homogeneous environment. The system in Fig. 1 with c → ∞ is then identical to a linear spin chain with two spin deviates (magnons) for which a bound state is not expected [53].
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 → −∞, which implies the formation of a tightly bound trimer (in analogy to trions in SU(3) symmetric Fermi gases [31]). 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 → −∞. In the rest of the paper, we only explore c −2.

B. Two Static Impurities
Before we discuss effective models that describe two mobile impurities in a sea of fermions, we consider two static impurities M → ∞ -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 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, ε 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 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 → ∞ can be obtained by fitting to the known form of the tail [54,55]. Figure 3 illustrates the energy E BO (r) = ε BO (c, r) − ε BO (0, r) − 2E static for c = −0.5 (we assume that g II = 0 for the sake of discussion). E static is the energy gain for immersing a single static impurity in a Fermi gas [56] The quantity E BO has a deep minimum at r = 0 given by E static (2c) − 2E 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 → 0 the tail can be written simply as c 2 cos(Ar + δ)/r, where A and δ are constants. This tail can be obtained from Friedel oscillations [29]. 3 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 [54].

IV. EFFECTIVE MODEL
A. Bethe-ansatz-solvable system The ground state of a one-dimensional Fermi gas (N f → ∞) 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 [57]. For the SU(2) symmetric case it can be conveniently studied using the BA equations [58]. 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 selfenergy can be used to describe the low-energy spectrum of a Fermi gas with one impurity [43,44,[58][59][60], suggesting the use of the following Hamiltonian to model the binding energy for the system of two impurities where m eff (c) is the effective mass of the impurity whose analytical form is presented in Refs. [44,58,59], so that the Hamiltonian h 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 ), The terms of the order O(c 5 ) can be neglected for interactions considered here, c ∈ (−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 impurityimpurity correlations. We choose it as 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 V (y)dy determines the binding energy for c → 0. Indeed, the ground state energy E eff of h eff for weak interactions reads [28,46] where the integration is over the whole system in the thermodynamic limit. This expression can be expanded as where we assume that V dy 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 → 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 V dy 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 where κ ≡ | V dy|. 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 κ can be obtained, for example, from a single-phonon exchange [13], in which case κ = 2c 2 /π 2 0.202c 2 for c → 0. If the potential is used in Eq. (11) then a single bound state with the energy is produced. This effective model captures qualitatively the exact results, see Fig. 4 (top). We show the ground state energies of h eff with m eff from Eq. (12) and with m 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, E BO (y)dy is welldefined. We calculate that | E BO (y)dy| 0.22c 2 for c → 0. This is slightly larger than κ 0.202c 2 for the zero-range potential discussed above. Even though, the long-range tail is not expected for integrable systems [13], 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 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 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 eff as a matrix and diagonalize this matrix after truncation. The energy is found by fitting to the form E K = E eff + D/K, where K is the number of used basis states and D, E eff are fitting parameters. This slow 1/K-convergence is expected for zero-range interactions [61].
Other potentials. The effective potential can also be calculated using other approximation schemes. For example, trial wave functions [16,25,62] 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.

B. 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 nonintegrable systems, i.e., c II = c and/or m = M . We   write the corresponding effective Hamiltonian as where we use m 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 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 cr II 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 = c and m = M , in which the particle exchange symmetry is broken by the interaction term 4 . 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 inmedium bound state is formed at smaller values of c in comparison to the BO potential, see Fig. 5 (top). The difference between c cr II 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: 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 [23]: 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 −0.06 (see Fig. 5 (center)), which for 6 Li atoms with ρ = 3/(µm) translates into 22 nK ×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 eff = M . To explicitly show this, we consider c II = c and m = 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 133 Cs atoms in a fermionic gas of 6 Li as in the experiment of Ref. [23], see Fig. 5 (bottom). We use both the Born-Oppenheimer potential as well as the zero-range potential. Note however that for M 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-ansatzsolvable 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-ansatzintegrable 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 < −π √ 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 = c and m = 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 [34,63].
It will be interesting to extend the present study to fermionic impurities. It is known that two fermionic im-purities in the SU(2) case do not have an in-medium bound state [64]. 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 → ∞, 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 [65,66].
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 [67]. 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 [16,56,61,[68][69][70][71][72][73][74][75]. In this appendix, we derive a weak coupling expansion of the BA equations, where N is an odd number. First, we consider the quasimomenta k j (j = 3, ..., N ) that satisfy k j (c = 0) = 0. For c = 0, these quasi-momenta are multiples of 2π/l. When c = 0, we write them as where δ j is small. Inserting Eq. (A3) into the left-handside of Eq. (A1) leads to e ikj l = e imj 2π =1 e iδj l ≈ 1 + iδ j l .
We write the right-hand-side of Eq. (A1) as where terms proportional to c n with n > 1 are neglected. Also it is assumed that k (0) j δ j −Λ 1 and k (0) j δ j −Λ 2 . This assumption is valid, since the Λ's lie in between the first three quasi-momenta, which are all close to zero. To derive Eq. (A5), we use that for a b With Eqs. (A4) and (A5), we obtain 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 Λ 1 = −Λ 2 , we consider Eq. (A1) for k 1 and k 2 . We use k 1 = −k 2 in the equation for k 2 : From this equation and the equation for k 1 : we obtain that Λ 1 = −Λ 2 . Using the equation for k 1 , we derive Next, we consider Eq. (A2). The sum of the quasimomenta {k 3 , ..., k N } is zero, thus, Eq. (A2) can be rewritten as With k 1 = −k 2 , k 3 = 0, this equation reads We rewrite Eq. (A11) which leads to k 1 = 3c/l and Λ 1 = c/l. 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 quasimomenta k 1 and k 2 are purely imaginary, they are shown in Fig. 6; the quasi-momenta k 3 , 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 −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 ≈ −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  solid curve is the weak-coupling expansion. The area 2 is the region where we use a solution at the previous point (see the green arrows) as an initial guess for Newton's method. To establish an initial guess for the quasi-momenta in the area 3, we construct a Taylor series using a number of already calculated points. The Taylor series constructed upon the points shown with crosses is presented as a (blue) dotted curve.
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 ≈ −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 ∈ (−2, −1.5) (not shown here for brevity).

Appendix C: Thermodynamic extrapolation
To extrapolate the calculated energies ε(c) − ε(0) to the thermodynamic limit, we shall employ the two fit functions: 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 (C1) and (C2) appear to represent the data well. They also produce similar results for N → ∞. The values of ε ∞ from the two fits differ only in the third digit, implying that the precise knowledge of the convergence pattern as N → ∞ 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 [43]. We use this well-known set-up to test our numerical approach. The system is desribed the Hamiltonian 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 [39] e ik (1) where k (1) j is the jth quasi-momentum (we use the superscript to emphasize that we work with a single impurity here), and Λ is one additional variable. We consider N to be even. Once the BA equation are solved, the energy can be calculated as ε (1) 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 [76] is used as where m j ∈ 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 √ c.
To extrapolate the result to the thermodynamic limit we use ε (1) (c) − ε (1) (0) = E + α/N β , where E, α and β 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. The analytic result for the thermodynamic limit [43] is shown by the orange line. In addition the total energy for a system consisting of N = 14 particles is shown by the blue dots. Inset: The red circles display the relative difference (x − y)/x, where x is our numeric result and y the analytic expression from Ref. [43]. The curve is added to guide the eye.
where k is the wave number. The wave function must obey the "delta-potential boundary conditions" at x = ±r/2, and periodic boundary conditions at x = ±l. We divide the solutions into parity-symmetric and parityantisymmetric 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κ, where κ is real. The equations that determine k for each class of solutions are written below.
2) Antisymmetric "bound states": To solve the equations, a genetic algorithm first finds approximate solutions for k, κ. These are then used in the Newton's iteration method as an initial guess. We calculate as many energy levels as particles we consider.