Three-body bound states of two bosonic impurities immersed in a Fermi sea in 2D

We consider two identical impurities immersed in a Fermi sea for a broad range of masses and for both interacting and non-interacting impurities. The interaction between the particles is described through attractive zero-range potentials and the problem is solved in momentum space. The two impurities can attach to a fermion from the sea and form three-body bound states. The energy of these states increase as function of the Fermi momentum $k_F$, leading to three-body bound states below the Fermi energy. The fate of the states depends highly on two- and three-body thresholds and we find evidence of medium-induced Borromean-like states in 2D. The corrections due to particle-hole fluctuations in the Fermi sea are considered in the three-body calculations and we show that in spite of the fact that they strongly affect both the two- and three-body systems, the correction to the point at which the three-body states cease to exist is small.


Introduction
The improvement of the techniques for cooling atoms has been boosting advances in physics for several years [1,2,3,4]. A remarkable example is the finding of experimental evidences of the so-called Efimov states, that were predicted by V. Efimov in the seventies [5]. Using few-body techniques, he showed that the energy levels of three identical bosons accumulate with a geometrical separation between them when the binding energy of the pairs goes to zero. More than thirty years passed until a signature of such states was found in loss experiments in 2006 [6]. Although being a true breakthrough, this was still an indirect measurement and another nine years were necessary until this effect was found in a direct manner in experiments with 4 He [7].
This impressive technological advance also opened new avenues as the quasi-2D realization of single-species systems [8,9,10], two-component gases [11,12], mixedspecies systems [13,14] and even heteronuclear diatomic molecules [15] in ultracold atomic traps. On the theoretical side, it has been shown that three-identical weakly attractive bosons in 2D have only one two-body and two three-body bound states [16] and further studies have also shown that the Efimov effect does not happen when the system is restricted to lower dimensions [17,18]. Recent studies have focused on mixed-species two-dimensional systems [20,19,21], which, despite of not supporting the Efimov states, can have a rich energy spectrum in the case of highly mass-imbalanced systems [22,23].
While most of these calculations were done in vacuum or disregarding the medium, bound states can be drastically affected depending on the surrounding medium in which the system is immersed. If one thinks about a single impurity immersed in a bath of different atoms, its motion distorts the medium, which in turn acts back on the impurity. The final state depends on the interaction strength between the different species, where for weak enough attraction the impurity behaves as a free particle and for strong enough attraction it forms a bound state with a particle from the sea. On the other hand, for intermediate attraction the impurity is dressed by the particles from the sea and becomes a polaron, namely, a quasi-particle whose basic properties as effective mass and charge differs from the ones of the original impurity. These different properties result from the dressing of the impurity by excitations of the background medium. The original idea of polaron dates back to Landau [24] and a discussion about it in the context of ultracold atomic gases can be found in [25].
It is possible to simulate polaronic systems with cold atomic gases by mixing a few atoms of a specific type to a bath of atoms of a different kind [26] (different masses or spin states, for example) even when the system is restricted to two dimensions and the majority atoms are fermions, in which case the quasi-particle is known as a Fermi polaron. Fermi gases have been successfully trapped to quasi-2D geometries and the Fermi polaron was already observed through Radio-Frequency Spectroscopy techniques [27,28,29,30]. The fate of a single impurity immersed in a twodimensional Fermi sea has been studied in recent years both experimentally [31] and theoretically [32,33,34,35,36], where the importance of the particle-hole fluctuations in 2D systems was emphasized.
In real systems it can be very hard to have only a single particle propagating in the medium and it is therefore important to understand how the presence of extra impurities affect the system. In this work we take a first step in this direction and consider the case of two identical impurities immersed in a Fermi sea. The study applies to a broad range of masses and both for interacting and non-interacting impurities. In both cases, the two impurities can attach to a fermion from the sea and form three-body bound states. We show how the energy of the states are shifted depending on the Fermi momentum k F and also how these states decay into the two-and three-body continuum. The corrections due to the fluctuations in the sea are considered when the impurities are not allowed to interact with each other. We show that in spite of the strong influence both the twoand three-body systems, the corrections from fluctuations in the Fermi sea have small influence on the point where the three-body bound state cease to exist. This paper is organized as follows: the influence of the background environment in the state of an impurity bound to a fermion from the sea is discussed in Sec. 2 and the derivation of the equations which consider corrections to the two-body system due to the particle-hole fluctuations in the sea are presented in Sec. 3. Sec. 4 brings a discussion about the dependence of the two-and three-body thresholds in the Fermi momentum. Three-body systems composed of two identical interacting impurities and a fermion from the sea are studied in Sec. 5. Corrections to the three-body states induced by the fluctuations in the sea are discussed in Sec. 6 for two non-interacting particles, where it is also argued how the corrections would affect the three-body system when the impurities do interact. Discussion and outlook are in Sec. 7 and technical details are given in Appendix A and Appendix B.

Bound and virtual states of the two-body system
We consider an atom labeled a of mass m a interacting with a fermion b of mass m b that is on top of a background Fermi sea of momentum k F . The interaction is described for the one term separable potential V = λ |χ χ|, it is attractive (λ < 0) and has the form factor g(p) = p| χ = 1 because we use a Dirac delta potential. The matrix elements of the two-body T-matrix [37,19] under the influence of the Pauli blocking effects are calculated as where p is the relative momentum, q is the center-of-mass (CM) momentum of the pair with respect to the Fermi sea, the reduced mass is m ab = m a m b /(m a + m b ) and E 2 is the internal (or binding) energy of the pair. The Θ−function ensures that the fermion is not inside the Fermi sea. The divergence in (1) is treated with the subtraction method [38], where the strength of the potential λ is connected to the binding energy of the pair |E ab | Notice that if the background Fermi sea is absent, k F = 0 and (2) recovers the two-body T-matrix elements of an ab system in vacuum, namely It is possible to identify three different thresholds in (2). As discussed in [34] for m a = m b , two of them come from the branch cut of the square root and are given by The third one arises from the branch cut of the logarithm and reads E 0 th = 0 for q ≥ ma m ab k F . The two-body system only supports bound states when (i) there is a pole in the two-body T-matrix and (ii) the energy where the pole appears is below the lowest threshold, namely E 2 < E − th or E 2 < E 0 th for q ≥ ma m ab k F . Setting (2) equal to zero, the two-body energies corresponding to bound states are described by the expression which tends to the energy of the pair in vacuum when the Fermi sea is not relevant, i.e., E 2 → −|E ab | for q/k F → ∞. This expression agrees with the ones calculated using many-body techniques, variational methods etc. [32,33,34,39,36]. The internal two-body energy given in (4) always satisfies the condition (ii), but surprisingly it does not satisfy the condition (i) for all possible combinations of momenta q and k F . There is a strong competition between binding, the Fermi sea and the CM momentum such that the pair is unbound for k F > 8m ab |E ab | when the value of these parameters fulfill the relation Although the analytical expression (4) is a solution of τ −1 (E 2 ) = 0, numerical calculations from (1) also do not indicate two-body bound states in the region given in (5), which has its boundaries framed by the vertical arrows in figure 1(b). The absence of bound states in the region defined by (5) was already pointed out in [34] and can also be seen in [40] for the spin-balanced system and the remaining question is to know what happens to the state in this region. It is shown in figure 1 that for k F > 8m ab |E ab | the state touches the lowest threshold (E − th ) in two points, strongly suggesting that it enters through the two-body cut in one point and comes back in another one. This cut is determined exclusively by the square-root inside the logarithm in (2) and the logarithm itself does not contribute to the cut since its argument is positive and non-zero throughout the regions where E 2 < E − th and E 2 < E 0 th for q ≥ ma m ab k F . The analytic continuation of the matrix elements (2) is labeled τ −1 * and found by flipping the sign in front of the square-root inside the logarithm in (2), since the logarithm itself does not contribute to the cut. Next, this expression is used to calculate E 2 as a solution of τ −1 * (E 2 ) = 0, which gives precisely (4) as the result. Notice that now the solution still satisfies condition (ii), but only satisfies condition (i) in the region given in (5) and not outside this region as before. These two-body states have entered through the cut and are in the second energy-sheet, therefore they are virtual states. These states are shown in figure 1(b).   The matrix elements of the two-body T-matrix as given in (2) may lead to inconsistent results in the scattering region (E 2 > E − th ), where the imaginary part of the energy becomes important. For instance, setting E 2 → E 2 + iǫ in (2), expanding the terms inside the square root and collecting energies and ǫ's together makes it possible to see that the imaginary part becomes null in some cases. If one disregards for a moment the small imaginary term that should appear inside the square roots and considers that the imaginary contribution due to these terms come only from the arguments inside it, then three cases have to be considered, namely, th . When the energy crosses from E 2 < E − th to E 2 ≥ E − th , the analytic extension of the square root has to be properly taken in account. In other words, its imaginary part has to follow the sign of the iǫ term that appears outside the square roots in (2), meaning that 1/2m ab (k F − m ab /m a q) 2 The same argument is taken in the transition from E − th ≤ E 2 ≤ E + th to E 2 ≥ E + th and therefore, the matrix elements are written as with g(E 2 , k F , q) given by for for E 2 > E + th . A comparison between the absolute value of the matrix elements with the analytic extension in (6), with E 2 → E 2 + iǫ in (2) and numerically calculated from (1) is shown in figure 2. The analytic extension becomes more relevant as E 2 approaches E + th . For E 2 < E + th the analytic form in (2) nicely reproduces the numerical calculation. The correction introduced by the analytic extension can be seen in the spectral molecular function, defined as A mol (E, q) = −2ℑ τ (E, q). Figures 3(a) and 3(c) show the spectral molecular function respectively for E ab E F = 0.1 and 2.5 calculated with (2), where a small imaginary part is added to the two-body energy. On the other hand, figures 3(b) and 3(d) are calculated with (6), where the analytic extension of the matrix elements has been properly taken into account. Notice that the use of the proper expression slightly changes the contrast in the figures. Figures 3(a) and 3(c) have been obtained in [34] and in order to straightforwardly compare figure 3 to the previous work, calculations here were made replacing the binding energy in (2) and (6) by the total two-body energy measured in terms of the Fermi energy E F plus the chemical potential of the impurity (µ a ), specifically an. extension pure numeric Figure 2. Comparison between the absolute value of the matrix elements with the analytic extension in (6), with E 2 → E 2 + iǫ in (2) and numerically calculated from  . Spectral molecular function A mol (E, q) for E ab EF = 0.1 (top) and E ab EF = 2.5 (bottom). The left-hand-side is the result of calculations with (2) and is the same as presented in [34]. The right-hand-side is calculated with (6), where the analytic extension of the matrix elements has been properly taken into account. Notice how the contrast change from left to right and that the corrected expression for the matrix elements is more relevant in the large binding limit.

Self-energy correction to the two-body system and Polaron
A way of taking into account fluctuations in the background medium is by introducing the self-energy of the impurity, which considers particle-hole excitations on top of the Fermi sea. The full propagator for the self-energy of the impurity relates to the bare one through the Dyson equation In the self-consistent way, the transition operator and the self-energy obey a set of non-linear coupled equations. This problem is handled by solving this set of non-linear coupled equations iteratively, which is equivalent to a Ward-Luttinger approach [41]. The zero-order transition operator is the bare one, given in (6) and the zero-order selfenergy is a function of the bare transition operator. The first-order transition operator is a function of the zero-order self-energy and so on (see Appendix A for the details). This procedure leads to the set of coupled equations where τ −1 0 (E 2 , q) is given in (6), m a is the mas of the impurity and m b the mass of the fermion.
The problem of one single impurity on top of a Fermi sea has been successfully solved with non-self-consistent methods (see, e.g. [34,36,39,42,43]). However, we found the self-consistent one more suitable to handle in the investigation of the threebody problem of two impurities immersed in a Fermi sea, whose results are presented in the following sections. This choice is supported since the self-consistent method describes well the two-body system when particle-hole fluctuations are taken into account. For instance, calculations done with the zero-order expression of the selfenergy (11 with n = 0) for m a = m b reproduce previous results in the literature. Using (11) and (12) we found that the energy of the attractive and repulsive branches of the polaron, which are solutions from E = q 2 / (2m a ) + ℜ [Σ 0 a (E, q)], as well as the weights of the branches ( a (E, q))]} −1 ) and the polaronic spectral function are the same as found in [39,34]. In [36,42] the molecular regime, namely the region where the energy of the molecule under the effect of particle-hole fluctuations from the medium (dressed molecule) is lower than the polaron one, calculated with non-self-consistent methods happens from large and negative η up respectively to -0.97 or -0.95, above which the polaron energy will always be below the molecular energy. We define This transition between regimes is often called polaron-molecule transition and we found a first-order (n = 0 in (11) and (12) ) polaronmolecule transition to occur around η = −0.54, which is below the results obtained with variational [36] and diagrammatic Monte Carlo [42] calculations, but is close to the η ≈ −0.60 found in the 2D limit of quasi-two-dimensional highly polarized Fermi gases [43]. The different behavior found between our method and previous studies in the small window −0.96(±0.01) < η < −0.54 [36,42] or −0.60 < η < −0.54 [43] is not important for our purpose of using the self-consistent approach in the three-body calculations. The actual point of crossing is not the relevant information since experiments [31] have shown that the crossing of the energies of the polaron and the molecule as function of η happens at a very shallow angle and that the two-body energy spectrum is not affected by the polaron-molecule transition (see [31] for details).
In view of that, it is more important that the model reproduces this transition than the actual point of occurrence. Solving (11) and (12) with n = 0 we find that the energy spectrum of the polaron and the dressed molecule are very similar to each other and that the crossing between the energies does happen at a shallow angle, as shown in figure 4. Here a difference between the self-consistent method and the non-self-consistent one can be seem. Comparing figure 4 to [43], we note that the polaronic spectrum is exactly the same, while the spectrum of the dressed molecule has a small quantitative difference. However, as explained above, since our model captures the correct behavior of the two-body system, this small quantitative difference is not a problem. As an extra argument for using the self-consistent method in our analysis, notice that our three-body results when particle-hole fluctuations are taken into account lie mostly in η > −0.54 (see figures 8 and 10), from where the energy of the polaron is always lower than the molecular one and it is the same as calculated with non-self-consistent methods [34,36,39,42,43].
It was found in [42], with diagrammatic Monte Carlo calculation, that corrections beyond one particle-hole in the medium does not contribute significantly to the energy spectrum of the polaron and the dressed molecule. In view of that, the convergence of the iterative solution of (11) and (12) is studied. Σ n a (E a , 0) and τ −1 n+1 (E 2 , 0) are shown as function of E 2 /E ab for different η's, where the strong binding-limit corresponds large negative η and the weak binding-limit to large positive η.
Results for large binding (η = −2) are shown in figure 5(a). Notice that the selfenergy converges very fast as Σ 1 a (E 2 ) and Σ 2 a (E 2 ) are almost indistinguishable from each other and both of them differ just slightly from Σ 0 a (E 2 ). As expected, the self-energy correction is almost negligible in this limit, which is confirmed by noticing that τ −1 0 (E 2 ), τ −1 1 (E 2 ) and τ −1 2 (E 2 ) are practically identical. Increasing k F and moving towards the small binding limit, a few iterations are still needed for the self-energy to converge. For η = 0, for instance, the difference from Σ 2 a (E 2 ) to Σ 1 a (E 2 ) is clearly much smaller than the difference from Σ 1 a (E 2 ) to Σ 0 a (E 2 ), as shown in figure 5(b). The difference is that now the self-energy influences the two-body system, rendering it more bound. There is a big difference between the pole position of τ −1 0 (E 2 ) and τ −1 1 (E 2 ), but this difference is already much smaller from τ −1 The scenario keeps this pattern for η > 0 and figures showing details are not presented here. The results in figure 5 provide our argument for working with τ −1 1 (E 2 ) and Σ 0 a (E 2 ) in the three-body calculations below. The computational cost of going beyond this are considerable and we do not expect neither qualitative nor large quantitative changes.
One great simplification comes by noticing that the function Σ is almost independent of the angles in its argument, namely where The agreement between the both sides of (13) is in general better than 0.1%. The approximation in (13)

Thresholds of the three-body system
We now treat the problem of two identical impurities of mass m a immersed in a Fermi sea background and interacting with a fermion of mass m b that belongs to the sea.
We assume that the total angular momentum is zero for both interacting and noninteracting impurities. As for the two-body system, the interaction between particles is assumed to be attractive s−wave zero-range potentials and now the center-of-mass (CM) of the three-body systems is considered at rest in the frame of the Fermi sea. The coupled homogeneous integral equations for the bound state of an aab system are derived in detail in [19]. Here we have a Fermi sea and need to take Pauli blocking into account. We do so by inserting appropriate Θ-functions in the coupled integral equations. The Faddeev components describing the system of two impurities immersed in the Fermi sea are given by where the three-body reduced masses are m aa,b = 2m a m b /(2m a + m b ) and m ab,a = m a (m a + m b )/(2m a + m b ). The ab and aa transition amplitudes are respectively given in (2) and (3) and are calculated for energies of the corresponding subsystems within the aab system. Before seeing how the presence of the Fermi sea affects the three-body states, it is important to understand how the two-and three-body thresholds move with the Fermi momentum. As the a particles are not directly affected by the Fermi sea, the transition operator of this subsystem is given in (3). The two-body breakup due to the interacting identical a particles (E th aa ), where aab → aa + b, is defined as the minimum energy which satisfies τ −1 aa E th aa − q 2 /2m aa,b = 0. The assumption that the three-body system is at rest in the Fermi sea frame implies that the total momentum of the pair is restricted by the momentum of the fermion (b particle), i.e., q ≥ k F . The momentum which minimizes the energy is k F and E th aa moves with k F as In 2D we set E aa = 0 when the a particles are not interacting with each other [44,45].
In this case the two-body transition operator in (14) diverges and decouples the homogeneous integral equations (14) and (15). Besides, the aa + b decay channel is not available and E th aa need not be considered in the three-body calculations. In the same way, the two-body breakup due to the interaction between the ab particles (E th ab ), where aab → ab + a, is defined as the minimum energy which satisfies τ −1 ab E th ab − q 2 /2m ab,a = 0. As the b particle does interact with the Fermi sea, the transition operator is given in (2) and the expression which satisfies τ −1 ab (E 2 ) = 0 is presented in (4). Replacing E 2 → E th ab − q 2 2m ab,a in this equation allows us to calculate the minimum three-body energy that gives a pole in the two-body T-matrix. Notice that there is a subtlety in this part since in cases where k F ≥ 8m ab |E ab |, the minimum energy which defines the two-body breakup does not correspond to an ab bound state (see figure 1(b)). When the two-body subsystem is virtual (unbound), the three-body system will only be stable while its energy is below the lowest two-body threshold. This condition is found by calculating the expression which minimizes the energy in the expression for E − th when E − th → E th ab − q 2 2m ab,a . Collecting everything, the dependence of E th ab with k F reads where k * = (2m 2 a |E ab |/m ab,a ) 1/2 . (17) is illustrated in figure 6. The three-body system has to be below either the ab subsystem threshold or the two-body continuum when the first one is not available.

The result in
Finally, the three-body breakup (E th 3B ), where aab → a + a + b, is defined as the minimum energy which renders one of the denominators on the right-hand-side of (14) and (15) zero. For k F = 0 this energy is E th 3B = 0 and a finite k F increases the breakup Bound state for k F > 8m ab |E ab | E − th for k F > 8m ab |E ab | Bound state for k F < 8m ab |E ab | Figure 6. Two-body total energy, E 2 /|E ab |, as function of q/ m a |E ab |. For k F < 8m ab |E ab | the relevant threshold for the three-body aab system is the bound ab subsystem, while for k F > 8m ab |E ab | it is the unbound subsystem.
value. Searching for the minimum energy which makes the denominator vanish on the first term of both (14) and (15), we have that one of the momenta q or k is allowed to vary in the range [0, ∞], while the Θ−function restrict the other one to [k F , ∞]. In both cases, the minimum energy that makes the denominator zero is E th 3B = k 2 F /2m ab . In the missing term, the second one on the right-hand-side of (15), both momenta are allowed in the range [0, ∞], but the Θ−function restricts this whole term to contribute only when q + k ≥ k F (see (B.7a)). The dependence of E th 3B on k F is therefore which is always below k 2 F /2m ab . Equation (18) thus defined the three-body threshold for the breakup aab → a + a + b.

Interacting impurities
Each impurity particle, a, interacts with a fermion from the sea by a zero-range potential characterized by a two-body binding energy E ab = (2m ab a 2 2D ) −1 , forming three-body bound states aab. These states can exist even when the impurities are not interacting with each other. In both the interacting (E aa = 0) and non-interacting (E aa = 0) cases, the three-body energies increase with increasing Fermi momentum k F . Since the two-and three-body thresholds are also moving up as k F increases, three-body bound states are found below the Fermi energy even with positive energies. As an example, we show in figure 7 the case where the two impurities have the same mass as the fermion (m a = m b ). The two impurities can, in principle, interact with each other with any energy. We set E aa = E ab , which for identical masses and k F = 0 leads to the well-known case of just one two-body state and two three-body bound states with energies E (0) [16,19]. As seen in figure 7, the energies and thresholds increase with increasing k F . For E aa = E ab , the decay channel aab → a + ab is always closed. The three-body threshold is highest for k F = 0 and moves up at the smallest rate. It goes below E th aa when k F / m a |E ab | > √ 2. It is interesting to note that the first excited state decays into an atom and a dimer when k F / m a |E ab | ≈ 1.1, while the ground state is breaking up into three free atoms for k F / m a |E ab | ≈ 5.1.
In view of what was discussed in the previous section, for k F ≥ k * = (2m 2 a |E ab |/m ab,a ) 1/2 three-body states of the two interacting impurities, which in 2D forms a L=0 bound state, and a fermion from the sea are supported even when two of the subsystems are in the virtual state, as shown in figure 7.  Figure 7. Three-body energy, E 3 /|E ab |, as function of the Fermi momentum, k F / m a |E ab |, for m a = m b and E aa = E ab The ↓ marks the point k F = k * . For k F > k * the two ab-subsystems are unbound.
The results in figure 7 indicate that cold atomic experiments may be able to deepen the knowledge of 2D three-body systems by looking at two impurities in the presence of a Fermi sea with which the impurities interact. The Fermi momentum is given by density of the sea, which is a controllable parameter. Then, it is possible to change it and measure the number of atoms lost in the trap, which must have a peak when the three-body system is close to either the two-or three-body continuum [6]. Different bound states may decay into different systems, as seen in figure 7 and a larger number of bound states can be reached in highly asymmetric mass systems (m b /m a ≪ 1) [22].
In order to give some insight for experiments, it is possible to systematically extend the procedure used in the investigation of the symmetric mass case and study how the three-body states vanish for a large range of mass-ratios. The final result is summarized in a mass versus k F diagram, shown in figure 8 for E aa = E ab . The Arabic numerals indicate the number of bound states in each region, where 0 means that no bound state was found, 1 indicates that only the ground state is available and so on. The central black-line divides the plot in two main regions where the three-body states go into either a dimer plus atom or three-atom continuum. Taking as example two 133 Cs atoms immersed in a sea of 6 Li fermions, the mass-ratio is m b /m a = 0.045 and four bound states are present for both E aa = 0 and E aa = E ab when k F = 0 [46,21]. It is shown in figure 8 that the four states for E aa = E ab disappears at η ≈ 1.96 (ground), η ≈ 0.59 (first), η ≈ 0.058 (second) and η ≈ −1.42 (third), respectively. Furthermore, notice that the three deeper states decay into three atoms while the highest one decays into atom plus dimer. The result shown in figure 8 is, to the best of our knowledge, the first one to consider the problem of two impurities immersed in a Fermi sea in 2D when the impurities are allowed to interact with each other. These results were achieved without taking into account particle-hole fluctuations in the Fermi sea. We now consider fluctuation corrections starting with two non-interacting impurities.

Self-energy corrections to the three-body system
Here we study the three-body system for two non-interacting impurities and then we see how the fluctuations in the Fermi sea affect the results. Taking E aa = 0 leads to f b ( q) = 0 in (14) and (15), simplifying the calculation. Furthermore, for m a = m b and k F = 0 only one bound state is supported [46] as shown in figure 9. Without considering the particle-hole fluctuations in the Fermi sea, the behavior of the three-body energy and of the thresholds is similar to the case where the impurities are allowed to interact. For k F / m a |E ab | > 0, the three-body threshold, which is the highest for k F = 0, increases at the smallest rate and for k F / m a |E ab | > 1.155 it becomes lower than the ab threshold, meaning that the decay channel aab → a + ab is not available any more for the three-body system. When k F / m a |E ab | ≈ 1.7 the three-body state touches the continuum and the system decays into three free atoms.
As in the case of interacting impurities, for k F ≥ k * = (2m 2 a |E ab |/m ab,a ) 1/2 the three-body ground state is supported even when the two ab-subsystems are in a virtual state. Since the last pair is also unbound (non-interacting impurities), the three-body system is bound when the three subsystems are unbound. The behavior of the ground state is similar to what has been seen in so-called Borromean systems in 2D [47,48].
In the present case we have the background medium in the form of a Fermi sea which is what makes this behavior possible. One may therefore consider this an example of a medium-induced Borromean state in 2D.
Next we introduce the self-energy correction in the three-body equations. Since the two impurities are not interacting, only the second term on the right-hand-side of (15) survives. It is important to emphasize at this point that the bosons in our problem are not from a Bose gas, but isolated impurities, or we can say that the average distance with a third boson from the gas is very large compared to any size scales in the problem (vanishing limit of the density). In this idealized situation, the connected three-body T-matrix can in principle contribute only to the two-boson T-matrix with a term proportional to the fermion density. If in addition, we take into account that the Bose gas density is negligible according to our statement of the physical situation, the contribution of the connected three-body T-matrix to the self-energy will be negligible and only the non-connected term of the T-matrix with the fermion and boson T-matrix will be relevant. Indeed, this term was already taken into account in our evaluation of the self-energy of the impurity. In a more general situation, e.g., in the presence of a Bose gas, other contributions to the impurity self-energy should be included, and the feedback from the full three-body T-matrix should be relevant as the boson density increases. In this case, one has to include corrections to the self-energy beyond the contribution of the disconnected part of the three-body transition matrix.
That said, the first correction comes in the two-body sector where the transition operator in (2) is replaced by the one in (12) with n = 0. As we discussed before, the convergence of the two-body T-matrix with n is very fast and the most significant correction is seen in figure 5 to be from τ −1 0 to τ −1 1 , justifying the choice of n = 0 in (12). In order words, the real gain of going to higher n is not worth the extra time spent in the calculation. Then, considering the dressed propagator of the impurity, the integral equation is written as with τ −1 1 given in (12). The numerical solution of (19) gives the corrected three-body energy and we need to understand again how the thresholds and energies depends on the Fermi momentum. However, since the self-energy correction is being considered, both the two-and three-body corrected thresholds are calculated numerically.
As shown in figure 5, the inclusion of the self-energy in the propagator of the impurity renders the molecule more bound, which lowers the ab threshold in the threebody calculation (see figure 9). The same happens to the three-body threshold, which starts from zero when k F = 0 and goes negative for k F > 0. The three-body state is also more bound and its energy as function of k F increases slower than in the previous cases, where corrections due to fluctuations of the medium were not considered. The three-body energy and the thresholds with and without the self-energy correction are compared in figure 9. The final result is that, although rendering both the two-and three-body bound states more bound, the inclusion of the self-energy makes the states disappear at a slightly smaller k F , since the changing rates of the thresholds are strongly affected.
The extension of the result in figure 9 for another mass ratios gives the diagram shown in figure 10, where the Arabic numeral indicate the number of bound states in each region. The continuous lines are from calculations without fluctuations in the Fermi sea and the discrete points show how the lines move when the self-energy is considered. Notice that the correction is always small within the range of the Fermi momentum and masses considered, which covers typically experimentally accessible cold atomic gas systems. Such correction is expected to be smaller in the case of interacting impurities, where the extra attraction would render the effects of the Fermi sea less relevant.
aab → a + a + b aab → ab + a 0 1 2 3+ Figure 10. Bound state diagram for E aa = 0 as function of the interaction strength η. The Arabic numerals indicate the number of bound states in each region. The "+" sign in 3+ indicates that more than three bounds states can be found in that region. The central black-line divides the plot in two main regions where the decay channel is either aab → a + ab or aab → a + a + b.

Discussion and outlook
We have considered the 2D problem of two identical atomic impurities (either bosonic or with distinct internal states), immersed in a background Fermi sea and interacting with a fermion on top of it in 2D. The interactions were modeled by attractive zero-range potentials and the Faddeev decomposition was used to write the homogeneous coupled integral equations for the three-body bound state.
The problem of just one impurity propagating in the Fermi sea was first considered and we have reviewed some well-known results in the literature [32,33,34,36,39], specifically the two-body T-matrix in the medium (2) and the energy of the molecule as function of its total momentum (4). The binding energy of the molecular state was studied previously as a function of its mass ratio, total momentum, and Fermi energy.
For some values of these parameters the pole of the T-matrix, which represents a bound state, enters through the cut of the two-body continuum and increasing the momentum of the molecule it returns from the cut and become a bound state again [34,40]. We have shown here that this state appears in the second energy sheet and becomes a virtual state, as shown in figure 1(b). Replacing E 2 → E 2 + iǫ in order to consider the two-body T-matrix in the medium (2) in the scattering region leads to inconsistent results in some cases (see figure 2). We have shown that the analytic extension of these matrix elements, as presented in (6), gives the right solution in this case.
Particle-hole fluctuations were then considered and the self-energy of the impurity propagating in the sea was self-consistently calculated. The non-linear couped equations for the transition operator and the self-energy (11) and (12) were solved iteratively and we showed that the convergence of both quantities with the number of iterations was very fast. Although this method has some disadvantages, we found it suitable to use mainly in the three-body calculation. This choice is supported since the self-consistent method employed in this work correctly describes the two-body system when particlehole fluctuations are taken into account. We found that the energy of the attractive and repulsive branches of the polaron, as well as their weights and the polaronic spectral function are the same as found in [39,34,43]. The polaron energy is the relevant information for η > −0.54, where most of the effects of the medium on the three-body system are calculated (see figures 8 and 10).
The formalism used here allowed us to study the complex problem of two interacting impurities propagating in the Fermi sea. Interestingly, three-body states of the two interacting impurities and a fermion of the sea are supported even when two of the subsystems are in a virtual state. Importantly, the fate of the three-body state depends strongly on the two-body systems. This feature can be used to identify the three-body states of two impurities immersed in a Fermi sea in measurements of cold atoms lost from a trap.
The complexity of the integral equations when the impurities are allowed to interact makes the inclusion of the self-energy correction very hard to be implemented. However, if the impurities are not interacting, the integral equations simplify and the effect of the correction can be studied. First of all, the three-body bound states are still allowed when the impurity-fermion subsystems are in a virtual state. Since we have two noninteracting impurities, the three-body systems are bound when the three two-body subsystems are unbound, which can be interpreted as a medium-induced Borromean state in 2D. The fate of the three-body states once more strongly depends on the twobody subsystems. The inclusion of the self-energy of the impurity drastically affects the behavior of the subsystems and of the three-body states, as it can be seen in figure 9, but the superposition of all the effects together leads to just a small correction in the final results, as shown in figure 10. We expect that the correction would be smaller in the case of interacting impurities, since the effects of the Fermi sea would decrease as there is more attraction in the system.
Three-body bound states were also found when a single impurity is immersed in a Fermi sea [36]. The different particle can be in a polaronic state or bind one or two fermions from the sea. It is interesting to understand what scenario is more favorable when another impurity is brought into the game. Would each impurity bind to two fermions or the two impurities bind to one fermion? Valuable information can also be gained by considering two polarons interacting with each other as a four-body system of two fermions and two impurities, using techniques similar to the ones employed in [49]. Another interesting direction for future investigations would be to study how the presence of the Fermi sea would affect the momentum distribution of the threebody systems by calculating the measurable two-and three-body contacts through an extension of the techniques described in [50]. Another interesting point is the understanding of how the signature of the Efimov effect in 3D three-body systems under the influence of a Fermi sea [51] would disappear and connect to the result presented here through a change in dimensionality, similarly to what was done for three-identical bosons in [52,53].
where the total energy and total momentum are respectively E = E a +E b and q = p a + p b . The integral I(E, q) reads and the expression for the self-energy is found to be The expressions for the transition amplitude and the self-energy, given respectively in (A.1) and (A.3), form a system of coupled equations. We solve this system iteratively by firstly setting the zeroth-order of the self-energy in (A.3) as where the zeroth-order transition operator τ 0 is the one given in (6) and is independent of Σ a . Then, Σ 0 a is introduced back in (A.1), leading to the first-order operator τ 1 , which requires (A.2) to be solved. Performing the integral on the energy in (A.2) is not trivial and an analysis of the poles has to be done. In the complex plane of E ′ b , the poles of I(E, q) are The pole in (A.5) is always on the lower half complex plane of E ′ b , as the only imaginary part comes from the small contribution iǫ and, for Σ a = 0, the pole in (A.6) is on the upper half plane. The self-energy is complex and contributes to the location of the pole. The study of its behavior and the results for the polaron spectral function in [34] show that ℑ(Σ a ) ≤ 0 for any (k, q, E, E ′ b ), ensuring that the pole in (A.6) is on the upper-half complex plane of E ′ b even for Σ a = 0. It is also necessary to study the analyticity of τ and Σ a as function of E and E ′ b . The coupled equations (A.2) and (A.3) indicates that both τ and Σ a must be analytic in the same region. From (1) is possible to see that τ , and consequently Σ, are analytic on the upper half complex plane of E and therefore on the lower half complex plane of E b . As the pole in (A.5) is in the semi-plane where the functions τ and Σ a are analytic, the contour is closed through the lower half complex plane of E ′ b (counterclockwise). The integral in (A.2) is calculated using the residues theorem and leads to the transition operator (see (A.1)) where λ −1 is related to the energy of the pair in the vacuum (see (1)). The set of decoupled equations for the self-energy and the transition operator are Σ n a (E a , p a ) = where τ 0 is given in (6) and is independent of Σ a .

Appendix B. Angular part of the integral Equation
Since we are interested in states with total angular momentum zero the spectator functions in (14) and (15) are independent of the angle between the momenta q and k. The angular dependence of the first term on the right-hand-side of these equations comes only from the denominator, as both the spectator (f ( q) ≡ f (q)) and the Θ−functions are angle-independent. However, the Θ−function in the second term of (15) does depend on the angle, which makes the integration of such term not as simple as the previous two. The trick to solve this integral is to transfer the information of the Θ−function to the limits of the angular integral. The angular dependence in the Θ−function on the second term of (15) is k 2 + q 2 + 2kq cos θ ≥ k 2 F , (B.1) which defines the angle θ 1 as the lowest limit of θ in (B.1) as with 0 ≤ θ 1 ≤ π. Changing integration limit from 2π to θ 1 , the angular integral in the second term of (15) reads where the main branch of tan θ has to be considered. Notice that for θ 1 = π, (B.3) gives the result for the angular integration of the first term on the right-hand-side of (14) and (15). The dependence of the Θ− function on the angle is eliminated using that | cos θ 1 | ≤

The result is
After the angular integration, the homogeneous coupled integral (14) and (15) read where the function g(q, k, k k ) reads