Functional-renormalization-group approach to strongly coupled Bose-Fermi mixtures in two dimensions

We study theoretically the phase diagram of strongly coupled two-dimensional Bose-Fermi mixtures interacting with attractive short-range potentials as a function of the particle densities. We focus on the limit where the size of the bound state between a boson and a fermion is small compared to the average interboson separation and develop a functional-renormalization-group approach that accounts for the bound-state physics arising from the extended Fr\"{o}hlich Hamiltonian. By including three-body correlations we are able to reproduce the polaron-to-molecule transition in two-dimensional Fermi gases in the extreme limit of vanishing boson density. We predict frequency- and momentum-resolved spectral functions and study the impact of three-body correlations on quasiparticle properties. At finite boson density, we find that when the bound-state energy exceeds the Fermi energy by a critical value, the fermions and bosons can form a fermionic composite with a well-defined Fermi surface. These composites constitute a Fermi sea of dressed Feshbach molecules in the case of ultracold atoms while in the case of atomically thin semiconductors a trion liquid emerges. As the boson density is increased further, the effective energy gap of the composites decreases, leading to a transition into a strongly correlated phase where polarons are hybridized with molecular degrees of freedom. We highlight the universal connection between two-dimensional semiconductors and ultracold atoms and we discuss perspectives for further exploring the rich structure of strongly coupled Bose-Fermi mixtures in these complementary systems.


I. INTRODUCTION
Ever since the theoretical explanation of conventional superconductivity as arising from the attractive interaction between electrons mediated by phonons [1,2], Bose-Fermi mixtures have been the subject of intense research. As they combine systems of different quantum statistics, their many-body behavior can be vastly different from that of the underlying bosonic or fermionic subsystems alone. Consequently, they can feature rich many-body physics ranging from superconductivity to the formation of composite bosonic or fermionic bound states similar to mesons and baryons in particle physics.
In solid-state physics, bosons typically appear as collective degrees of freedom. These may be, for instance, phonon excitations of an underlying crystalline lattice or collective excitations of the electronic system itself in the form of, e.g., magnons or plasmons. Beyond such systems, experimental progress in the fields of atomically thin semiconductors [3] and ultracold atoms [4] makes it now possible to enter a new regime of strongly coupled Bose-Fermi mixtures. Here -akin to the physics of nuclear matter-fermions and bosons appear on equal footing, both representing pointlike particle degrees of freedom.
Crucially, direct pairing between bosons and fermions is a new essential ingredient in these mixtures. Recently it was shown [5] that for such strongly coupled Bose-Fermi mixtures a description in terms of Fröhlich or Holstein models [6,7], in which fermions couple linearly to the bosonic degrees of freedom, fails. In addition, the coupling to bosons at quadratic order becomes relevant, which has to be accounted for in an extended Fröhlich Hamiltonian [5], giving rise to qualitatively new physics recently observed in experiments in cold gases [8][9][10] and Rydberg systems [11].
Various aspects of atomic, three-dimensional Bose-Fermi mixtures have been investigated theoretically using the Fröhlich model -thus disregarding the crucial quartic interaction term. This revealed a rich structure of the phase diagram ranging from polaron formation [12][13][14] and boson-induced p-wave superfluidity [15], to phonon softening and phase separation [16].
Similarly, the phase diagram of two-dimensional Bose-Fermi mixtures has been explored using the Fröhlich model. These studies were motivated in particular by exciton-electron mixtures in semiconductors, and, following initial work by Ginzburg [17], it was predicted that the system may turn superconducting [18][19][20] while other works proposed a transition to supersolidity [21,22], or that the formation of both phases might be intertwined [23].
Due to the shortcomings of the Fröhlich model and mean-field inspired approaches that neglect pairing [24][25][26][27][28][29], these initial studies missed the fact that the microscopic interaction between atoms in ultracold gases and between excitons and electrons in semiconductors is fundamentally attractive. While in cold gases interactions arise from long-range van der Waals forces, the polarization of charge-neutral excitons by electrons gives rise to attractive forces in semiconductors. Crucially, in both cases the interactions support bound states between the fermionic and bosonic particles. Consequently, as the strongly coupled regime is entered, one has to consider the extended Fröhlich Hamiltonian in order to account for the pairing to fermionic Feshbach molecules in cold arXiv:2104.14017v2 [cond-mat.quant-gas] 8 Mar 2022 atoms and exciton-electron bound states, called trions, in semiconductors.
The presence of this novel bound-state physics renders the description of strongly coupled Bose-Fermi mixtures an outstanding theoretical challenge. This is reflected by the fact that until now -except for initial studies in three dimensions [30][31][32][33][34][35][36][37]-the phase diagram of strongly coupled Bose-Fermi mixtures as function of the density of bosons n B and fermions n F , schematically shown in Fig. 1, remains unexplored. With the discovery of atomically thin transition-metal-dichalcogenides, the semiconducting class of layered van der Waals materials, the exploration of this phase diagram in two dimensions becomes particularly urgent. This is not only due to the potential of layered materials for technological applications, but also due to the possibility of realizing long-lived, stable exciton-electron mixtures that feature a striking similarity to cold atomic mixtures [38][39][40]. This universal connection, detailed by the comparison of typical scales in both systems shown in Table II below, opens the possibility to explore emerging phases in strongly interacting systems in two complementary and seemingly disparate systems, that while playing on vastly different energy and length scales, are governed by the same dimensionless system parameters.
In this work, we study theoretically the phase diagram of strongly coupled two-dimensional Bose-Fermi mixtures as a function of the boson and fermion densities.
A key theoretical challenge is that the pairing between bosons and fermions gives rise to fermionic composite particles. Due to their fermionic nature, these particles evade conventional mean-field approaches and are thus much harder to describe than their bosonic counterparts in Fermi mixtures, where they emerge as Cooper pairs or bosonic molecules. Moreover, the existence of such fermionic composites implies a phase diagram that is richer in possible phase transitions compared to the simpler Fröhlich model. Here we tackle this challenge by developing first steps towards a comprehensive functionalrenormalization-group approach that allows access to the full phase diagram of Bose-Fermi mixtures in two dimension. Our approach accounts for the bound-state physics arising from the extended Fröhlich Hamiltonian and can be systematically extended to describe the plethora of competing phases illustrated in Fig. 1.
In order to explore this phase diagram it is crucial to start from limits that allow for a controlled understanding of the physics involved. One such limit is found at extreme population imbalance where just a single boson is immersed in a fermionic bath. This so-called Fermi polaron problem already displays rich physics that has been studied extensively in three dimensions [42][43][44][45][46][47][48][49][50][51][52][53][54]. Here one finds that as the interaction between the impurity and the bath is tuned, the system undergoes a sharp transition from a polaronic to a molecular state. While in the polaron state the impurity is essentially weakly dressed by bath excitations, in the molecular state the impurity binds tightly to one fermion of the surrounding environ- At strong-coupling the system is described by the extended Fröhlich model that accounts for the formation of a two-body bound state between fermions and bosons of energy B . The limit nB = 0 (along the y-axis) defines the Fermi polaron problem, discussed in Section III, where a single bosonic impurity interacts with a fermionic bath. In this limit the impurity can either bind with a fermion into a molecule or remain unbound as a Fermi polaron. At finite boson density, discussed in Section IV, we find a transition from a molecular phase which hosts a Fermi sea of bound molecules (light gray) to a mixed phase in which a condensate of bosons hybridizes fermionic and molecular degrees of freedom (red/dark shading). With the exception of the extreme limit nF = 0 that corresponds to the Bose polaron problem (along the x-axis), as the boson density is increased beyond the regime nB nF (red/dark to blue/light shading), the phase diagram remains largely unexplored. Starting with the possibility of bipolaron formation [41], various competing phases can be conjectured based on studies of the simpler weak-coupling Fröhlich model, ranging from supersolid charge density wave states [21,23] to boson-mediated s/p-wave fermion pairing [15,16,18,19]. ment giving rise to a state that, close to the transition [52,53], is orthogonal to the polaron state.
The two-dimensional case has received attention over the last decade [55][56][57][58][59][60][61] as well. It turns out that this case is more challenging to describe due to the increased significance of quantum fluctuations in reduced dimensions. While early works based on simple variational wave functions found no polaron-to-molecule transition [55], later studies showed that this finding was in fact an artifact caused by the neglect of three-body correlations. Including these, one indeed recovers a polaron-to-molecule transition in two dimensions [56,59], a result supported by subsequent studies using a variety of Quantum Monte-Carlo (QMC) techniques [58,60,61].
As the preceding discussion shows, there are strong constraints on any approach that aims to reliably describe strongly coupled Bose-Fermi mixtures in two dimensions even on a qualitative level. First, in order to address the strong-coupling character of the problem correctly, it must be based on the extended Fröhlich model. Second, it must go beyond perturbation theory in order to describe the formation of fermionic bound states. Third, at vanishing boson density it must correctly reproduce the quantum impurity limit, which necessitates the incorporation of three-body correlations. Fourth, in order to describe the phase diagram at finite boson density n B , the approach must be able to deal with the fermionic nature of the composite particles that will experience Pauli blocking at finite density similar to baryons in atomic nuclei.
All these requirements are met by the functional renormalization group (fRG). Based on an implementation of Wilson's renormalization group idea, the fRG has been successfully applied to the study of strongly coupled systems in a broad range of areas [62][63][64][65], spanning from the asymptotic safety of quantum gravity [66][67][68] to highenergy [69,70], statistical [71][72][73] and condensed matter physics [74][75][76][77][78]. In addition to addressing the aforementioned constraints imposed by the two-dimensional polaron problem, the fRG technique developed in this work displays several other advantages. First, in contrast to variational approaches based on particle-hole excitation expansions, it provides a fully self-consistent approach that naturally includes high-order quantum fluctuations and treats polaron and molecular states on equal footing. Second, compared to conventional quantum field theory approaches the fRG includes quantum fluctuations in a coarse-grained fashion -momentum-by-momentum shell-that makes it ideally suited to treat competing ordering instabilities. Third, similar to variational techniques, the fRG can be improved systematically by using increasingly refined truncations of the underlying quantum effective action. Finally, it offers an easier access to spectral and dynamical response functions compared to Monte Carlo approaches where the analytic continuation of noisy data is required.
We demonstrate the applicability of our approach by focusing on the case where the size a B of the fermionic bound state is small compared to the average distance d ∼ n −1/2 B between bosons. Since for sufficiently shortranged attraction this bound state always exists in two dimensions [79], its binding energy B = 2 /(2µa 2 B ) (with µ the reduced mass) is the relevant interaction scale, i.e. we work in the limit ( 2 /2µ)n B / B 1. By including the full feedback of three-body correlations on the renormalization group flow, we demonstrate the correct description of the polaron-to-molecule transition in the single boson limit. In particular we predict the transition to occur at a critical dimensionless interaction strength ( F / B ) * = 1/18.78 in excellent agree-ment with state-of-the-art variational [56,59] and diagrammatic MC approaches [58,60,61].
Having thus established the limiting case of the phase diagram, we extend the renormalization group (RG) flow to finite boson density. At small dimensionless Fermi energies F / B , we find that fermionic composites build up a well-defined Fermi surface, leading to the formation of a trion liquid in the case of semiconductors and a Fermi sea of dressed Feshbach molecules in the case of ultracold atoms. As the boson density is increased, the effective energy gap of the composites decreases, leading to a transition into a strongly correlated phase where fermions are hybridized with molecular degrees of freedom. This extension of a single boson framework does not take into account the formation of higher-order bound states including more than one boson [80][81][82]. While this description might thus be missing some of the phases and states at play, recent theoretical and experimental results suggest that this simplified treatment may, however, still be sufficient to describe the physics relevant on experimental times scales [37,83].
Adapting the fRG approach to account for the full frequency-dependence of self-energies, we predict the spectral properties of the model. We find that the inclusion of three-body correlations has a strong impact on the effective masses of polarons and molecules (trions) which can be observed using state-of-the-art experimental techniques recently developed in ultracold atoms [84][85][86].
The paper is structured as follows: in Section II we introduce the strong-coupling model of Bose-Fermi mixtures and discuss the effective action formalism. Here we also introduce our fRG approach, derive the corresponding renormalization group equations and discuss how the various phases discussed in this work can be distinguished. As this section contains a detailed discussion of the used technique, readers mainly interested in the predictions of our work may proceed from the introduction to Section III and the sections thereafter. In Section III we discuss the universal connection between strongly coupled Bose-Fermi mixtures in atomically thin semiconductors and ultracold atoms. We benchmark our approach on the limiting case of a single boson embedded in a fermionic environment, obtain the ground-state energy of the system and study the evolution of correlation functions in dependence on the fermion density and interaction strength. In Section IV we turn to the case of finite boson density. We determine the phase diagram both as a function of the chemical potential and density of both species. In Section V we adapt the fRG scheme to describe the spectral functions of the model and we predict the properties of quasiparticles emerging in the theory. We conclude in Section VI, discuss perspectives for possible experimental realizations and provide an overview of open questions and promising extensions of the fRG approach introduced in the present work.

II. MODEL
We consider a two-dimensional Bose-Fermi mixture consisting of a fermionic species ψ into which bosonic particles φ are embedded. The system is described by the microscopic action where x = (r, τ ) denotes the coordinate r and imaginary time τ ∈ [0, 1/T ]; moreover, x = 1/T 0 dτ d 2 r. In the following, we consider zero temperature, T = 0, and assume that bosons and fermions have the same mass m = m F = m B . We work in units = k B = 1, and set 2m = 1 unless indicated otherwise. The fields ψ and φ are of fermionic Grassmann and complex boson nature, respectively. The two species interact by means of an attractive contact potential of strength g < 0. The model is regularized in the ultraviolet (UV) by a momentum cutoff Λ.
The densities of both species are set by the chemical potentials µ ψ/φ [87]. At a finite fermion density n ψ (set by a chemical potential µ ψ > 0), tuning the boson chemical potential µ φ at fixed µ ψ and g, triggers a transition at a critical chemical potential µ c φ between a vacuum phase of bosons (µ φ < µ c φ ) with vanishing boson density to a phase of finite boson density n B > 0 (µ φ > µ c φ ). For strongly coupled Bose-Fermi mixtures it is crucial to allow for the possibility of the pairing of the bosons and fermions to a composite fermionic molecular (trion) state. In order to describe this bound state, it is essential to resolve the pole structure of the scattering vertex sufficiently well [44,88]. In order to achieve this in an efficient way, rather than considering the action in Eq. (1), we study a two-channel model where the interspecies interaction is mediated by a molecule field t, that describes a composite fermionic particle of mass 2m [89][90][91][92]. The action is given by Here, a boson and a fermion can be converted into the molecule (trion) t with a conversion Yukawa coupling h, and m t is the detuning energy of the molecule. In Eq. (2) we give the action in Fourier space where P = (p, ω) comprises the momentum p and the Matsubara frequency ω, and p,ω ≡ d 2 pdω. We operate in the limit where h → ∞ which universally describes both openchannel dominated Feshbach resonances in cold atoms [4] as well as electron-exciton scattering in atomically thin transition-metal dichalcogenides [38]. In this limit, t becomes a purely auxiliary Hubbard-Stratonovich field, i.e. it can be integrated out to yield back the original action (1) when h 2 /m t = −g is fulfilled [45,93].
In two dimensions, a bound state exists for any attractive interaction strength g < 0 [79]. Using a sharp UV cutoff in the Lippmann-Schwinger equation, the binding energy B is related to the parameters of the microscopic model through [55,79,94] Thus, rather than using the microscopic coupling g (or equivalently h and m t ) we can parametrize the interaction strength in terms of the experimentally measurable binding energy B of the molecule (in the case of cold atoms) or trion (in the case of 2D semiconductors), respectively. Note, in the following we will use the terms trion and molecule often interchangeably.

A. fRG formalism and effective action
The fRG is a momentum space implementation of Wilson's renormalization group. In the following, we briefly recall its principle; for a detailed discussion we refer to Refs. [62][63][64]77]. The idea behind the fRG is to build a family of theories indexed by a momentum scale k such that only quantum fluctuations above that scale are taken into account. Thus rather than treating fluctuations at all scales at once, one iteratively integrates out modes from high to low energies by smoothly lowering k from the microscopic UV scale Λ down to k = 0.
In practice this is done by adding to the action (2) an infrared regulator term which penalizes low-energy fluctuations, such that only high-energy modes contribute to the field integral. For bosons and fermions at vanishing density the lowenergy modes are located at small momenta. Thus the cutoff function R σ,k (P ) (σ = ψ, φ, t) is set to be large (wrt. k 2 ) for |p| k and negligible for |p| k. In this way, low-momentum fluctuations are suppressed while high-momentum ones are left unaffected. For fermions at a finite density, the low-energy modes are located around the Fermi surface. Accordingly, in this case R σ,k is chosen to suppress fluctuations of modes inside a momentum shell of width ∼ 2k around the Fermi surface.
Starting from the sum of S and ∆S k one then defines a scale-dependent partition function Z k , as well as a scaledependent effective action Γ k through a (modified) Legendre transform of the free energy ln Z k . The evolution or 'flow' of the effective action as the scale k is lowered is then given by the Wetterich equation [95], In the above expression, the supertrace STr denotes a summation over all momenta and frequencies, as well as the different fields, including a minus sign for fermions. Moreover, Γ k and R k represent the matrices of second functional derivatives of Γ k and ∆S k , respectively, with respect to the quantum fields [96].
Provided R k=Λ = ∞ [97] at k = Λ all fluctuations are suppressed and Γ k=Λ = S + const. [63,64]. On the other hand, for R k=0 = 0 one recovers at k = 0 the effective action of the original model, Γ k=0 = Γ. Crucially, the effective action Γ (Gibbs free energy) is the generating functional of all one-particle irreducible vertices. It thus contains all information about the exact solution of the theory and hence its determination corresponds to solving the non-relativistic, many-body Schrödinger equation.

B. Truncation schemes
While the flow equation (5) is exact, it is, in most practical cases, impossible to solve without resorting to approximations. A standard strategy is to propose an Ansatz for the flowing effective action Γ k . When dealing with fermions, it is customary to expand in the powers of the fields in a so-called vertex expansion [77]. Following this strategy we choose the Ansatz for the field-dependent part of the effective action We perform an additional gradient expansion by neglecting a possibly emerging momentum dependence of the Yukawa coupling h k via vertex corrections. Each field σ = ψ, φ, t carries renormalized flowing single-particle Green's functions whose momentum dependence is approximated within the gradient expansion as parametrized by inverse quasiparticle weights A σ,k and detunings m σ,k . Note that for the boson field φ we have absorbed the dependence on the chemical potential µ φ into the definition of the detuning m φ,k for convenience. The Ansatz (6) incorporates in detail two-body correlations between the bosons and fermions. In particular, it describes well the pairing correlations between the particles which is essential to enter the strong-coupling regime. As a short-hand we refer to the effective flowing action (6) as the 'two-body truncation'. The two-body truncation has been used successfully to study the Fermi polaron problem in three space dimensions [53,98,99]. In two space dimensions, however, quantum fluctuations are stronger and previous works [55,56] have established that higher-order correlations must be taken into account to describe the ground state of the system. Indeed, as we shall see in Section III, the two-body truncation is not sufficient to describe the polaron-to-molecule transition.
Consequently, we extend the Ansatz for the effective action to a 'three-body truncation'. To this end we add a term to the two-body truncation that accounts for the build up of three-body correlations during the RG flow: The additional term proportional to the contact coupling λ k describes the scattering between composite molecules and fermions, and thus, by virtue of the tree-level diagram depicted in Fig. 2(a), it accounts effectively for the emergence of three-body correlations in the system. Let us briefly comment on the validity of the gradient expansion used for both truncations (6) and (10). In the single-boson limit, we expect the low-energy excitations of the boson φ and the composite particle t to be at small momenta and we may thus expand the momentumdependence of their propagators in a power series about p = 0, ω = 0. For the fermions, on the other hand, we expect the most relevant excitations to lie around the Fermi surface. We thus expand their propagator about p 2 = F , ω = 0.
As we extend our calculation to a finite boson density we retain the expansion around p = 0 for the molecules as we will find that their phase appears in a regime of the phase diagram where n B n F . Thus the Fermi energy of molecules always remains small. Moreover, we employ a gradient expansion that neglects effective mass corrections as these are not expected to be crucial to correctly capture the qualitative physics of the phase diagram (except for large mass ratios m F m B [100], a regime not considered in this work).
In the quantum impurity limit, the vanishing of the boson density implies that the properties of the fermionic Green's function are not affected by interactions; i.e. the propagator in Eq. (7) with A ψ,k = 1 and m ψ,k = 0 is exact. This can also be verified explicitly from the flow equations derived further below [cf. Eqs. (15) to (18) and (E2)]. At finite boson density we neglect the renormalization of the fermionic propagators since throughout this work we will remain in the regime of density ratios n B n F . While the truncation in Eq. (10) can be improved systematically, e.g., by considering higher-order correlations or a more involved momentum dependence of the propagators or the vertices, the model in Eq. (10) is sufficient to accurately describe the intricate quantum impurity limit, as shown in Section III. In particular, even though h k and λ k have no momentum dependence, the Bose-Fermi scattering T -matrix, as described by the exchange treelevel diagram shown in Fig. 2(b), acquires a momentum dependence due to the dynamic field t that is sufficient to describe accurately the Bose-Fermi scattering at the relevant energy scales.
We note that at finite boson density our truncation does not account for the possible formation of bound states between two or more bosons and a single fermion [80][81][82]. In a realistic experimental setting, where the system will be prepared adiabatically, the formation of these higher-order bound states requires several bosons to be located in the close vicinity of the fermions. Since we focus here, however, on the regime where the boson density is significantly smaller than the fermion density, n B n F , the probability to find such configurations will be small. As a result, compared to the time scale of Fermi polaron or molecule formation, the formation of higher-order bound states will be sup-pressed, enabling the observation of the phase diagram studied in this work on transient time scales. Nevertheless, while recent results suggest that this treatment is appropriate [37,83], the framework used in this work can be extended to feature bound states between two bosons and a fermion; both in the vacuum limit as well as at finite density (for details see Appendix E). This highlights that this study provides only an initial step in the exploration of this phase diagram which, given sufficiently stable bound states, may feature an even richer structure.

C. Regulators
For the regulators R σ,k we use sharp cutoff functions [77], defined so that the regulated inverse flowing propagators appearing on the rhs. of the flow equation (5) acquire the simple form Here, F,k = µ ψ − m ψ,k /A ψ,k is the Fermi energy of the fermionic species ψ. For the fermions ψ the regulator suppresses fluctuations at momenta in a shell of width 2k around the, in principle, flowing Fermi-surface of the bath [101]. Even though the molecule is a fermion as well and thus may develop a Fermi surface at finite boson density, we regulate it about zero momentum as all phases considered in this work appear in the regime n B n F . The choice of sharp cutoff functions has several advantages [102]. Foremost, it allows for an analytic derivation of the flow equations. In addition, it facilitates the comparison to previous FRG studies [53,98,99] as well as to self-consistent diagrammatic approximations that display a similar mathematical structure [5].

D. Flow Equations
We now turn to the explicit derivation of the RG equations [103,104] of all running coupling constants. For the three-body truncation Γ 3,k (Γ 2,k is a subset obtained by setting λ k ≡ 0 in all flow equations) all vertices can be expressed in terms of the six running couplings A σ,k , m σ,k , h k and λ k . Following the prescription detailed in Appendix A 1, the flow equations are obtained from appropriate functional derivatives of the Wetterich equation. Their diagrammatic representation is shown in Fig. 3 and in terms of the flowing Green's functions they read and In these expressions∂ k stands for the derivative with respect to the k dependence of the regulator only, i.e.
As discussed in Appendix A 2, from Eqs. (15) to (17) the flow equations of the couplings A σ,k , m σ,k are obtained by projection onto the momentum dependencies given in Eqs. (7) to (9).

E. RG initial conditions
The initial conditions for the flows are obtained by setting Γ k=Λ = S +const. First, we discuss the UV initial conditions for m t and A t which are obtained as follows. The two-body problem of a single boson scattering with a single fermion can be solved exactly. The resulting initial condition for m t,k=Λ is given by Eq. (3). To arrive at this expression one may recognize that in the twobody problem the molecule is the ground state. As such it has to be a gapless degree of freedom in the infrared, i.e. m t,k=0 = 0. Moreover, µ ψ and µ φ must be set to negative values µ vac ψ and µ vac φ to yield a vanishing density of either species. In addition, µ vac ψ + µ vac φ = − B has to be fulfilled to ensure that the energy cost to create two particles from the vacuum to form a bound state is given by the molecular binding energy B . Using these conditions, together with the fact that in the two-body problem neither the boson and fermion propagators G ψ,k and G φ,k nor the vertices h k and λ k renormalize, the flow of m t,k can be solved analytically to yield Eq. (3) (for more details we refer to Appendix B).
We work in the limit of large h k=Λ which ensures t to be purely an auxiliary field and we use A t,k=Λ = 1. Furthermore, we set λ k=Λ = 0 as it does not appear in the classical action in Eq. (2). The initial condition for the field renormalization of the boson field is naturally given by A φ,Λ = 1, and the UV value of its detuning is set by the boson chemical potential, m φ,Λ = −µ φ . Finally, since we will study only phases at small ratios n B n F we can assume that the fermion field is not renormalized, i.e. A ψ,k = 1 and m ψ,k = 0 throughout the RG flow.

F. Chemical potentials and distinction of phases
The numerical integration of the flow equations yields the physical value of the propagators and interaction vertices at the infrared scale k = 0. Depending on their properties we can distinguish various states and phases of the strongly coupled Bose-Fermi mixture, summarized in Table I. In the single-boson limit, yet at finite fermion density, we distinguish two states: a molecular state in which the boson is paired into a composite particle, and a polaron state where the boson is dressed by fluctuations of majority fermion particles. At finite boson density, we distinguish two phases: a molecular phase, where all bosons are paired into fermionic molecules, n t > 0 and n φ = 0, and a mixed phase where molecules and unpaired polarons coexist [30]. In the mixed phase, n φ > 0, so that the condensate of bosons creates a bilinear coupling in the effective action ∼ h √ n φ (t * ψ + c.c.) leading to a hybridization of the fermions with the molecular degree of freedom. This means that no purely polaronic phase with n t = 0 and n φ > 0 is possible. In the limit of n B → 0 the mixed phase connects to the polaron state, whereas the molecular phase connects to the molecular state.
In order to differentiate between these states and phases we consider the different densities defined by integrals proportional to p,ω G σ,k=0 (p, ω). These densities are nonzero only when poles of G σ,k=0 (p, ω) lie in the upper half of the complex ω-frequency plane. Hence, from the location of poles, manifest in the energy gaps of the particle in the infrared, we can determine whether the corresponding densities vanish.
Specifically, the boson vacuum corresponds to a finite excitation gap for both the boson and the molecule, m φ,k=0 , m t,k=0 > 0. Likewise, in the single-boson limit the ground state has to be gapless while the excited state is gapped since this limit marks the boundary between the boson vacuum and the many-boson regime.
At finite boson density, the molecular phase corresponds to m φ,k=0 > 0 and m t,k=0 < 0, i.e. molecules feature a Fermi surface determined by their Fermi energy −m t,k=0 /A t,k=0 . For the mixed phase, the situation is more subtle. Our Ansatz does not allow for the description of a condensate at finite boson density that could be accounted for, e.g., by shifting the φ-field expectation value by a coherent state transformation. However, it is still possible to predict whether a boson condensate forms. Indeed, a necessary condition for the existence of a φ-condensate is that for some 0 ≤ k ≤ Λ the boson gap m φ,k vanishes [105]. In that case, even though we are unable to further pursue the RG flow, we identify the phase to be the mixed phase.
In this mixed phase the bilinear term mentioned above leads to a mixing of the fermionic and the molecular propagators. Consequently, these propagators share the same pole structure and the corresponding species are thus populated simultaneously. As a result all three particle species are present in this phase. This implies that in our model a regime populated exclusively by majority fermions and condensed minority bosons is possible only in the single-boson limit at n B = 0.

III. QUANTUM IMPURITY LIMIT: SINGLE BOSON IN A FERMI SEA
We first apply our approach to the limiting case of the Bose-Fermi phase diagram where an individual boson is immersed in a bath of fermions. This limit defines the so-called Fermi polaron problem, and its solution determines the phase diagram along the y-axis of Fig. 1. In order to reach this single-boson limit, the boson chemical potential is tuned to the critical value µ φ = µ c φ that separates the boson vacuum (µ φ < µ c φ ) from the phase of a finite boson density (µ φ > µ c φ ); see Table I.

A. Fermi polaron problem in ultracold atoms and atomically thin semiconductors
The nature of the ground state of the Fermi polaron problem universally depends on the ratio of the two relevant energy scales of the problem: the kinetic energy, represented by F , and the interaction energy, set by B . While F / B can, in theory, be tuned by adjusting either F or B , in experiments it depends on the physical system which parameter is accessible for easy tunability.
The two main systems in which strongly coupled Bose-Fermi mixtures can be realized today are ultracold atoms and atomically thin semiconducting transitionmetal dichalcogenides (TMD). To support the following discussion, in Table II we summarize key parameters and quantities describing the universal connection between these systems.
In monolayer TMD, B represents the trion binding energy which is typically fixed [38,39,[106][107][108]. However, by electrostatically doping the system with charge carriers, the Fermi energy F is easily adjusted and F / B can thus be tuned. In cold atoms the situation is reversed. Here, the binding energy B can be tuned using Feshbach resonances, while adjusting the Fermi energy over a wide range of values is challenging. As a result, in cold atoms the Fermi energy F is the natural unit and, correspondingly, the spectrum of the system is expressed as a function of the dimensionless energy E/ F and interaction strength B / F . In contrast, in TMD the binding energy B provides the appropriate unit, and the spectrum is expressed as a function of E/ B and F / B .
Of course, physics does not depend on the chosen units. It is, however, still instructive to compare spectra for both sets of units, as the choice of units reflects the experimental protocols employed to observe the physics of Fermi polarons: in TMD using gate-doping of F and in

B. Quasiparticle energies
In order to obtain the spectrum of the Fermi polaron problem we first determine the ground-state energy of the system, set by µ c φ , the critical energy needed to bring a boson from the vacuum. The procedure is summarized in Fig. 4: when the polaron is the ground state, m φ,k=0 = 0, and the polaron energy is given by E pol = µ c φ . In this 'polaron regime' the molecular state is an excited state whose energy is determined from the pole of its Green's function which yields E mol = E pol + m t,k=0 /A t,k=0 . In turn, in the 'molecular regime' the molecule is the ground state. Here, m t,k=0 = 0, and the molecule energy is given by E mol = µ c φ , while the polaron is an excited state with an energy gap E pol = E mol + m φ,k=0 /A φ,k=0 .
In Fig. 5 we show the polaron and the molecular energy as obtained from the two-and three-body truncations. The spectrum of the Fermi polaron problem is shown both in units convenient for cold atoms [ Fig. 5(a)] as well as 2D materials [ Fig. 5(b)]. The comparison of (a) and (b) demonstrates that despite the fact that both panels contain fully redundant information, they yet represent seemingly different behaviour which is, however, solely due to the different choice of units.
In Fig. 5 the results obtained from the two-body truncation (6) are shown as dashed lines. This truncation takes into account a similar set of diagrams as a nonself-consistent T -matrix approach [57] which, in turn, is equivalent to a variational Chevy approach [44,55]. By contrast to the aforementioned approaches our fRG is self-consistent. As expected from these approaches, we find that the two-body truncation is not sufficient to generate a polaron-to-molecule transition.
Instead we find that the inclusion of irreducible threebody correlations is crucial, which is consistent with diagrammatic MC [60] and higher-order variational approaches [56,59]. We find that the inclusion of the three-    Table III.
Similar to previous field-theoretical or variational approaches [55][56][57]59], we do not include all possible twobody correlations and focus on the effect of pairing correlations. Further two-body correlations can, for instance, be generated by the re-emergence of the four-point vertex ∼ γψ * ψφ * φ. One may justify the exclusion of this vertex by an analogy to BEC superconductivity. There the vertex γ accounts for induced interactions in the particle-hole channel, leading to a contribution similar to the Gorkov corrections to BCS superconductivity [109][110][111]. In the BCS case, it leads to an effective shift of the inverse dimensionless interaction strength that appears in the gap equation determining T c /T F . Based on this analogy, we expect that such terms will not establish a new polaron-to-molecule transition, but rather only shift the location of an already present transition. Thus we concur with previous studies that it is three-body corre-lations that are essential to establish the formation of a phase of trions in strongly coupled Bose-Fermi mixtures [112].
We note that at low Fermi energies, we find a weak non-monotonous behaviour of the polaron energy in the dependence on F / B . Such a behavior is not present in works using variational [55][56][57]59] or MC approaches [60,61]. As discussed in Appendix C, we attribute this effect to the limited resolution of the frequency-and momentum-dependence of the vertex functions in both our truncations. This effect is, however, not relevant for our study of the Bose-Fermi phase diagram which depends only on the relative energy gaps between the polaron and molecular state and not on their respective absolute values.
have scaled both vertices by powers of h that reflect the scaling of the vertices with the molecular wave function renormalization A −1 t,k=0 yielding results independent of h in the contact-interaction limit at h → ∞.
Atom-molecule scattering.-The vertex λ k describes the scattering between the composite fermionic molecules and the excess fermions in the system. During the RG flow, λ k evolves from λ k=Λ = 0 in the UV to a negative value in the infrared at k = 0. Thusλ yields an attractive contribution, shown in Fig. 2(c), to the overall atom-molecule scattering amplitude that has an additional, significant contribution from the tree-level φ-exchange diagram depicted in Fig. 2(d). Fig. 6 shows the absolute value of the scattering vertex in the three-body limit (dashed orange line) where it takes the valueλ =λ (3B) = − F / B , for details see Appendix D. Thus the vertex scales proportional to the square of the size of the molecular bound state a B ∝ 1/ B . The solid orange line shows the result for λ in the polaron problem. At small fermion density the molecule is the ground state. In this 'molecular regime' the density of fermions is so low that the average interfermion spacing greatly exceeds the molecular size a B . Thus the atom-molecular scattering vertex is essentially unaffected by the presence of the fermionic medium, and λ follows the three-body resultλ (3B) .
As F / B is increased we observe a suppression of the atom-molecule scattering vertex. We attribute this effect to two contributing factors. First, the molecule becomes an excited state beyond the critical interaction ( F / B ) * . In this case the molecule is gapped and within our FRG approach which projects vertex functions on vanishing external vertex frequencies and momenta (see Appendix A),λ is thus suppressed by the molecular energy gap. More importantly, however, as the Fermi energy becomes larger than B , F / B > 1, the size of the bound state starts to exceed the typical inter-fermion distance. As a consequence, in-medium effects come into play leading to significant modifications ofλ. Indeed, these corrections become so strong thatλ starts to decrease at even larger values of F / B . Molecular gap.-The dimensionless molecular gap m t = m t,k=0 /h 2 k=0 is shown as a blue line in Fig. 6. For interaction strengths F / B < ( F / B ) * where the molecule is the ground state, the molecule is gapless, m t = 0. Beyond the transition the molecule becomes an excited state and we find that its gap vanishes linearly as The corresponding crossing of the molecular and the polaron state can also be interpreted as leading to an effective Feshbach resonance in the polaronfermion scattering where the tree-level diagram shown in Fig. 2(b), evaluated on-mass-shell, diverges. The associated polaron-fermion scattering length changes sign at the transition, with a positive value signaling the existence of a fermionic bound state.
In turn, within a single-channel theory that is formulated purely in terms of the 'atomic fields' ψ and φ, the divergence of the effective polaron-fermion scattering ver-tex ∼ h 2 /P t signals the instability towards a phase of fermionic bound states. In this language, entering this phase at finite boson density would necessarily require the introduction of the emergent fermionic composite states. Finally we note that in Fig. 6 we show only results from the three-body truncation Γ 3,k since in the two-body truncation Γ 2,k , the vertex λ k = 0 vanishes by definition throughout the RG flow. Moreover, since no polaron-to-molecule transition is present in this simpler truncation,m t always remains finite.

IV. BOSE-FERMI MIXTURE AT FINITE BOSON DENSITY
We now turn to the mixture regime, where a finite density of bosons interacts with a bath of majority fermions. As discussed in Section II B, within our truncations we can identify two phases: a molecular phase, where all bosons are bound into molecules, and a mixed phase where molecules are hybridized with majority fermions and coexist with a condensate of polarons.
While we can describe the molecular phase directly, we can not fully access the regime in which a condensate of polarons exists since this would require to explicitly include the condensate and thus an effective potential for the bosonic field. However, we can still determine the critical system parameters at which the system becomes unstable towards condensation. Indeed, the associated phase boundary is determined by the vanishing of the scale-dependent boson gap m φ,k /A φ,k at the end of the RG flow.
For large values of the boson chemical potential µ φ , the underlying assumption n B n F is no longer valid. When this condition breaks down, we thus terminate the fRG flow. While this does not define a phase, we dub this part of the phase diagram the 'stopped flow region', further discussed below.

A. Phase diagram as a function of chemical potential
In Fig. 7 we present the phase diagram of the Bose-Fermi mixture for both Γ 2,k and Γ 3,k [Eqs. (6) and (10)] at a fixed Fermi energy F , as function of B and µ φ . In the three-body truncation Γ 3,k [ Fig. 7(a)] a molecular phase forms at finite boson density in the interaction regime where the molecule is the ground state of the quantum impurity limit discussed in Section III.
In fact, the ground-state energy of the quantum impurity limit determines the chemical potential µ φ = µ c φ ( B , F ) that separates the vacuum of bosons from the mixed phase or the phase of a finite density of molecules. Along this phase boundary the system undergoes a transition from a polaronic to a molecular ground state.
In the interaction regime B / F > ( B / F ) * , increasing the boson chemical potential starting from values µ φ < µ c φ leads to a boson-vacuum-to-molecule transition as µ φ crosses the critical chemical potential. Directly on the critical line one enters the quantum impurity regime and a single molecule forms [115]. Increasing µ φ beyond µ c φ one enters the molecular phase where a finite density of bosons, all bound into molecules, exists. In this phase m t,k=0 /A t,k=0 < 0, and the molecules acquire a Fermi surface. Tuning µ φ further to larger values one reaches the phase boundary to the mixed phase. Here, a finite density of molecules coexists with gapless boson particles.
For B / F < ( B / F ) * there is no molecular phase and one transitions directly from the boson vacuum to the mixed phase. As the flow is terminated at a finite RG scale k end once the boson becomes gapless m φ,k=k end /A φ,k=k end = 0, at the boundaries to the molecular phase and to the boson vacuum phase the boson turns gapless at the end of the flow at k end = 0. Moving further into the phase from these boundaries the value of k end at which the flow is terminated increases.
When the flow is stopped in the mixed phase, the molecules might have already formed a molecular Fermi level during the course of the RG flow. This is indicated by the dashed gray line in Fig. 7. Above this line the molecule has developed a Fermi surface when the flow ends or is terminated at k = k end . Below the line the molecule has remained gapped. As expected, for B / F > ( B / F ) * this line parametrizes the bosonvacuum-to-molecule transition. For B / F < ( B / F ) * on the other hand, it bisects the mixed phase. These regions then correspond to phases of a single Fermi sea (boson vacuum), two Fermi seas (molecular phase), two Fermi seas with a bosonic condensate (mixed phase above the gray dashed line) and a bosonic condensate with only a single Fermi sea (mixed phase below the gray dashed line) as discussed in Refs. [30,116].
Increasing the bosonic chemical potential µ φ further within the mixed phase, the bosonic density increases until eventually the molecular Fermi wave vector becomes larger than the fermionic Fermi wave vector (−2m t,k /A t,k > F ). Within this regime, the bosonic density has become comparable to the fermionic density. This means that it is no longer justified to neglect the renormalization of the fermionic Green's function and to disregard higher-order correlations along with subdominant interaction channels. As we expect that in this case our truncation no longer renders an appropriate description of the system, we terminate the flow at finite scale k end once −2m t,k /A t,k > F . When this happens during the RG flow, a molecular Fermi sea has already formed while the bosons are still gapped m φ,k /A φ,k > 0. This 'stopped flow region' (mint in Fig. 7) occurs after the bosonic chemical potential has been tuned well into the mixed phase. We therefore expect that in the stopped flow region, close to the boundary to the mixed phase, the system would still be in a mixed phase, if one were to continue the flow.
Within the two-body truncation [see Fig. 7(b)] it is unsurprising to see that no molecular phase forms at finite

B. Phase diagram as a function of density
In the previous subsection results were given as a function of chemical potential. Experimentally it is, however, often simpler to determine the density of particles instead of their chemical potential. Thus, to make direct connection to experiments, it is useful to also consider the phase diagram as a function of particle densities. Since in the effective action formalism employed in this work, the chemical potentials are the parameters of the theory, the canonically conjugate densities have to be computed explicitly.
In principle, the fermion and boson densities can be determined directly from the two-point Green's functions. Within the derivative expansion and two-channel model an alternative approach is, however, more convenient.
Here one makes use of the fact that the densities are connected to the derivative of the effective potential U evaluated at the equilibrium field configuration σ eq by the standard relation Here n B and n F , respectively, denote the total density of bosons and fermions in the system, including those bound into molecules. The effective potential U (σ eq ), in turn, is obtained from the derivative-free part of the infrared effective action evaluated at the field expectation values U (σ eq ) = Γ k=0 [σ eq ]/(V /T ) where for the considered phases σ eq = (ψ eq , φ eq , t eq ) = 0. In the absence of approximations, determining the densities from the effective potential or from the Green's functions are equivalent methods, as follows from the Luttinger theorem [30,117]. Within our fRG scheme we, however, expect it to be computed more accurately using the flow of U than using the flow of G σ as this approach relies on lower-order vertices.
In the fRG, the effective action is promoted to a flowing effective action that depends on the RG scale k. Accordingly, it is convenient to define corresponding scaledependent densities n F/B,k and to determine the densities of the systems from their value at the end of the RG flow. The resulting density values are then associated with the corresponding phases. Since there is no polaron-to-molecule transition for Γ 2,k , in the following we discuss only results obtained in the three-body truncation.
The flow equation of the effective potential is obtained by evaluating the Wetterich equation (5) at vanishing fields, where ξ φ = 1 for bosons and ξ σ = −1 for fermions (ψ and t). Due to the pole structure of the integrand, Eq. (21) can be simplified further (for details see Appendix F) to Here, the step functions Θ σ,k (p) originate from the sharp regulators in the flow equations [see Eqs. (12) to (14)].
For the bosonic and molecular field they are defined as Θ t,k (p) = Θ φ,k (p) = Θ(p 2 − k 2 ), while for the fermionic field Θ ψ,k (p) is defined in the following in Eq. (23). As the scheme described in Section II does not feature a renormalization of the majority propagator it is evident from Eq. (22) that, within that approximation, the fermions do not contribute to the flow of the effective potential U k . Consequently, from the integration of Eq. (22) the density of fermions would not be calculated accurately since the depletion of majority carriers, resulting from fermions being bound into molecules, is not taken into account.
In order to take this effect into account, we deriveseparate from the flow of the Green's functions of the bosons, molecules and the interaction vertices-a flow equation for the propagator of the majority species, that does not feed back into any flow other than that of the effective potential. Since the majority fermions have a finite density already at the start of the RG flow, we regulate the fermions around their flowing Fermi level [101]. Accordingly, the step function in the first line in Eq. (22) is given by To derive the flow equations of A ψ,k and m ψ,k , we evaluate the RG flow of the associated vertex function at external frequency and momentum (p 2 , ω) = ( F , 0), i.e., we perform the gradient expansion around the bare Fermi surface of the majority species. This flow is then used to determine the effective potential U , and, in turn, the boson and fermion densities through Eq. (20). In order to reproduce the majority carrier density F /4π in the UV with regard to Eq. (20), the initial condition for the density flow is given by the mean-field result U k=Λ (σ eq ) = − 2 F /8π. In Fig. 8 we show the resulting phase diagram of the system as a function of the boson and fermion density. It can be regarded as the counterpart of Fig. 7(a), expressed in different variables. Specifically, to obtain Fig. 8, for the combinations of boson chemical potential µ φ and interaction strength B / F that lie in the molecular phase we computed the corresponding values of n B / B and n F / B . For combinations that lie inside the mixed phase or the stopped flow region we can not compute the boson and fermion density as the flow is terminated at finite k end . In Fig. 8 we thus identify density combinations outside the molecular phase as being in the mixed phase [118].
In Fig. 8, the single-boson limit discussed in Section III corresponds to the y-axis at n B / B = 0, and the polaron-to-molecule phase transition occurs at n F / B = ( F / B ) * /4π = 0.00424. As the boson density is increased, the mixed phase becomes favorable, i.e., the maximal density of fermions for which all bosons are bound into molecules decreases. We find that there is also a minimal fermion density required to enter the molecular phase. Below that critical value one again enters the mixed regime.

C. Mean-field model
Remarkably, a simple mean-field-inspired argument can provide an approximate phase diagram of the model: in the single-boson limit, the polaron is a gapped excitation in the molecular regime. It has a gap ∆E = E pol − E pol which is a function of F / B , or equivalently n F / B (equal to F /4π B along the y-axis in Fig. 8). This gap was determined numerically in Section III where we found, (24) reflecting that the energy gap vanishes at the polaron-tomolecule transition and attains a value proportional to F in the strong-binding, low-density limit.
In our mean-field model of the molecular phase, the interactions are taken into account by considering the effective Hamiltonian where ε k = k 2 /2m. Even though H M F is quadratic in the fields, this effective model goes beyond naive meanfield as ∆E incorporates the non-trivial solution of the polaron problem obtained through our fRG scheme in Section III. The polaronic, mixed phase appears when it is energetically unfavorable to bind into molecules, i.e. when the Fermi energy F,t of the molecules is larger than the gap ∆E. When this condition is reached the polarons start to form a condensate, as described previously. For a molecular Fermi energy below the gap ∆E, the ground state of the mean-field model (25) is given by separate Fermi seas of densities n ψ = F /4π and n t = F,t /2π for the fermionic and molecular sectors, respectively. Hence, in the molecular phase, the total bosonic and fermionic densities are given by n B = n φ + n t = n t and n F = n ψ + n t = n t + F /4π. The mean-field transition line below which the molecular state is favored is thus parametrized by This mean-field phase boundary is shown as a dashed line in Fig. 8. While the mean-field picture is oversimplified and does not correctly capture the quantitative renormalization effects beyond the vacuum-to-molecule transition, it correctly captures the qualitative nature of the structure of the phase diagram. The phase boundary, by construction, reaches the y-axis at the polaronto-molecule transition and approaches the origin at an angle of about n F /n B ≈ 2.22 which directly follows from the behaviour of the polaron gap ∆E(( F / B ) → 0) ≈ 0.41 F .

V. QUASIPARTICLE PROPERTIES OF POLARONS AND MOLECULES IN THE QUANTUM IMPURITY LIMIT
The calculations presented in Sections III and IV only yield information about ground state properties of the system. In order to extract spectral information such as dispersion relations, particle lifetimes, effective masses or higher-lying excited states, however, the spectral functions need to be computed.
The spectral functions are obtained from the Green's functions G ψ,φ,t by analytic continuation of the Matsubara frequencies iω → Ω + i0 + which yields the retarded Green's functions G R ψ,φ,t (p, Ω). From this, the momentum-and frequency-resolved spectral functions are obtained as Two difficulties arise when determining the spectral function within the fRG. First, an analytic continuation has to be performed, either at the level of the flow equations [119][120][121] or the final output of the RG flow in the infrared [53]. Second, in order to capture non-trivial spectral functions one needs the full momentum-and frequency-dependence of the propagator, which the gradient expansion employed in Sections III and IV does not provide. A solution to the latter difficulty can be found, e.g., by the direct implementation of fully frequencyand momentum-resolved Green's functions [53] or in the BMW scheme [72,122], which also yields a full momentum-and imaginary frequency-dependence of the propagators. Both these approaches, however, do not resolve the analytic continuation issue. For this reason, we implement here a method developed in nuclear physics [121,[123][124][125] which was recently applied to the polaron problem in three dimensions [98]. In the following we shall refer to this method as the frequency-and momentum-resolved scheme (FMR).
In FMR, the flow equations [Eq. (5) and Eqs. (15) to (18) and (E2)] are analytically continued to real frequencies. In order to achieve that, rather than projecting the flow equation onto the gradient expansion parameters, we retain the full momentum-and frequencydependence of the single-particle Green's functions on the lhs. of the flow equations, while we keep the gradient expansion for the two-body [Eq. (6)] and three-body truncation [Eq. (10)] on the rhs. of the equations. This enables us to perform the loop integration over imaginary frequencies analytically. In turn, this allows us to perform the analytic continuation to real frequency to obtain direct access to the retarded Green's functions. From that we evaluate the single-particle spectral function using Eq. (28); for further details we refer to Appendix G. We remark that, when applying a non-self-consistent implementation of FMR -in which only bare quantities appear on the rhs. of the flow equations-to the spectral function of the molecule, the differential equation sys- Mol.  tem yields the same results as a corresponding T -matrix resummation [57] (see Appendix G 2). Polaron spectral function.-The polaron spectral function obtained using FMR is shown for different interaction strengths in Fig. 9. Subfigures (a), (c) and (d) are obtained in the three-body truncation Γ 3,k . Subfigure (b) shows the result from the two-body truncation Γ 2,k in order to highlight the effect of the inclusion of irreducible three-body correlations.
The polaron spectral functions show the same qualitative behavior as the corresponding spectra in 3D [53,98]. Two quasiparticle peaks -the attractive and the repulsive polaron-can be discerned, and a molecule-hole continuum in between these dominant excitations is visible. The attractive polaron is the ground state in Fig. 9 (a), (b), (c), and thus is a gapless excitation. In contrast, in Fig. 9(d), the ground-state is a molecule, and thus a small gap at p = 0 can be seen. Generally, at finite but small momenta the attractive polaron is a welldefined quasiparticle with an interaction-dependent effective mass which, along with the effective masses of the repulsive polaron and the molecule, is shown in Table IV. For larger momenta, the attractive polaron peak eventually merges with the molecule-hole continuum, such that it is no longer a well-defined quasiparticle.
The repulsive polaron appears at energies above the scattering threshold (indicated by the dashed horizontal lines in Fig. 9) as a narrow peak, indicating a long quasi-particle life-time for the interaction strengths shown. Consistent with Ref. [57] we find that as B / F decreases, the repulsive polaron gradually disappears. Moreover, while at small interaction strength ( Fig. 9(a)) the repulsive polaron eventually merges with the moleculehole continuum at finite momentum, at larger interaction strength the repulsive polaron peak remains distinct from the molecule-hole continuum at any momentum and thus keeps a long life-time at high momenta.
As evident from the comparison of Fig. 9(b) and (c), the inclusion of the irreducible three-body correlations moves the molecule-hole continuum to lower energies. This has the effect that the dispersion relation of the attractive polaron becomes flatter, increasing the polaron effective mass compared to the two-body truncation (see Table IV). Furthermore, its quasiparticle peak joins the continuum at lower momenta. In Table V the energy of the repulsive polaron is shown relative to the groundstate energy. For the repulsive polaron the inclusion of three-body correlations has the effect of slightly altering its effective mass and of lowering its energy relative to the scattering threshold. For a fermionic impurity this indicates a reduced tendency towards itinerant Stoner ferromagnetism [126].
Molecular spectral function.-In Fig. 10, the molecular spectral function is shown for different interaction strengths B / F . Here the subfigures (a), (c) and (e) in the left column are obtained in the three-  body truncation while (b), (d) and (f) in the right column result from the two-body truncation. It can be seen that a general feature of this spectral function is spectral weight that appears above a parabola centered around p = k F and that is defined by the frequency Ω p = (p − k F ) 2 + m φ,k=0 /A φ,k=0 as derived in Appendix G 3. The quasiparticle peak of the molecule follows a distorted dispersion relation which, in the strongbinding limit, tends to a free-molecule dispersion relation. Dependent on the interaction strength, at low momenta the molecular quasiparticle peak lies outside of the particle-particle continuum and joins the continuum at finite momenta just to leave it again at higher momenta. More specifically, at low B / F , the quasiparticle peak joins the continuum at a low momentum, which increases with interaction strength B / F . Likewise, the momentum at which the peak leaves the continuum again increases with B / F as well. Similar to a non-self-consistent T -matrix resummation, in our approach the molecular quasiparticle peak has a vanishing width when it is not embedded in the continuum. This can be seen analytically by inspecting the flow equation of the two-point function G R t (see Appendix G 3 for details). Apart from the structure originating from the parabola-shaped particle-particle continuum and the quasiparticle peak, further structure exists within the parabola that originates from contributions in the RG flow where the Feynman diagrams are evaluated close to their poles (see Appendix G 3).
The minimal energy of the parabola Ω p is equal to the renormalized energy gap of the polaron, indicating a close relationship between the polaron at p = 0 and the molecule at p = k F , supporting the argument that both of these states overlap with the actual groundstate of the system and possibly with each other [52]. This finding can also be understood conceptually in a meanfield picture where a bosonic minority particle at p = 0 along with a majority fermions at the Fermi surface can be interpreted as either a polaron at p = 0 or a molecule at p = k F (previously noted by Cui [127]). Note that, because the particle-particle continuum in the molecular spectrum is shifted due to the renormalization of the boson gap, this effect can not be captured in a non-selfconsistent approximation such as employed in Ref. [57]. In such an approximation spectral peaks distinct from the continuum are present, that in our implementation are a part of the continuum.
Within the spectral functions obtained using Γ 2,k , the quasiparticle peak at p = 0 -located at approximately m t,k=0 /A t,k=0 -is always at a finite energy whereas using Γ 3,k it is moved closer to Ω = 0 and eventually attains Ω = 0 past the polaron-to-molecule transition. At the same time, the minimum of the parabola, given by m φ,k=0 /A φ,k=0 , detaches from Ω = 0 as the polaron is no longer the ground state. Hence using Γ 3,k the effective mass (see Table IV) of the molecule, which is negative at small B / F , diverges with increasing B / F and eventually becomes positive at an interaction strength before the polaron-to-molecule transition. Beyond the transition the molecule is gapless at p = 0 and its effective mass is positive. Using Γ 2,k , increasing B / F makes the molecule dispersion flatter leading to an increasingly negative effective mass.

VI. CONCLUSION
We investigated the phase diagram of strongly coupled Bose-Fermi mixtures in two dimensions. In order to make progress in the exploration of this complex phase diagram it is important to establish limits that can be understood controllably. To this end we focused on the regime of fermion-dominated population-imbalance which, in the extreme imbalance limit, connects to the Fermi polaron problem where a single bosonic impurity interacts with a Fermi sea. The opposite limit of a fermionic impurity coupled to a Bose-Einstein condensate corresponds to the Bose polaron problem which features qualitatively different physics. Already this asymmetry reflects the impact the interplay of different particle statistics has on the phase diagram away from the extreme population imbalanced limits.
In order to approach the problem we employed a functional-renormalization-group approach that allows to systematically incorporate high-order correlation functions. This enables us to reproduce the polaron-tomolecule transition in the single-boson limit which is a necessary condition for any theoretical approach that aims to describe this strong-coupling phase diagram. In contrast to the simpler three-dimensional case [42,49,53], we showed that three-body correlations have to be included to describe the polaron-to-molecule transition in two dimensions and we obtain excellent agreement with ab-initio approaches [60] that can be applied in the quantum impurity limit.
Using the fRG we extended the analysis to finite boson densities. There, depending on the boson and fermion densities (or equivalently their chemical potentials), we observed two phases: a fermionic liquid with two Fermi seas in which all bosons are bound into molecules, themselves immersed in a majority Fermi sea, and a hybridized liquid in which the condensation of bosons leads to a mixing of the fermionic and molecular sectors [30].
This hybridization and the associated mixing are not a result of the Hubbard-Stratonovich field used in our two-channel model, but they occur equally in atomic single-channel models whenever scattering vertices between fermions and bosons develop a pole in presence of a boson condensate. In this regard, the phase diagram away from the molecular phase at n B n F shares a remarkable similarity to the Bose polaron problem that describes the opposite limit of few fermions immersed in a Bose condensate, where the same hybridization mechanism leads to a crossover between the polaron and molecule instead of a transition [5,10].
Naively, one may suspect that in a mixture of bosons and fermions as many particles as possible are bound into fermionic bound states in order to maximize attractive potential energy. This, however, does not take into account the properties of the system in two ways. First, this argument neglects the fermionic nature of the bound states which leads to the formation of a molecular Fermi energy, representing a kinetic energy cost. As a result, when the bosonic density is increased, the molecular Fermi energy eventually exceeds the energy of the lowest-lying polaron state and the system enters the mixed phase.
Second, the argument misses the fact that already in the limit of a vanishingly small boson density the formation of a bound state competes with the formation of a polaron state in which a single boson interacts col-lectively with a large number of surrounding fermionic bath particles [39]. For a fixed interaction strength, the polaron state can thus profit more efficiently from an increased density of bath particles. Vice versa, as the bath density is lowered polaron dressing looses efficiency so that eventually the composite bound state becomes the new ground state (in absence of Coulomb interactions).
While the fRG approach employed in this work provides nontrivial insights into the phase diagram of the Fermi-Bose mixture, the approximations used are insufficient to explore the phase diagram in its whole richness. First hints to a plethora of exciting phenomena can already be inferred from numerous quasiparticle features of the single-particle spectral functions uncovered using the FMR scheme in Section V, ranging from non-trivial effective mass renormalization and the non-monotonous dispersion of molecules, to incoherent parts in the spectral function reflecting quasiparticle instability.
Indeed, for a more accurate description of such features it would be necessary to go beyond the gradient expansion we impose on our Ansatz and instead allowing for an arbitrary momentum and frequency dependence of vertex functions. While such a treatment has been used in the three-dimensional case [53], it remains challenging to implement numerically. Preliminary results [128], however, suggest that including the full momentum-and frequency-dependence indeed cures the spurious non-monotonous behaviour of the polaron energy discussed in Section III. Far from being only of quantitative importance, such a fully momentum-and frequency-resolved approach could give new qualitative insight into the phase diagram, e.g. by allowing for the description of transitions to non-trivial molecular Fermi surface topology [129] akin to Fulde-Ferrell-Larkin-Ovchinnikov phases in BCS superconductors [130,131].
We did not include Bose-Einstein condensation in our formalism. Its explicit inclusion would allow for the study of subregions of the mixed phase in which a bosonic condensate is accompanied by molecular or fermionic Fermi seas. Additionally, the presence of a condensate will require the incorporation of a repulsive Bose-Bose interaction to ensure the mechanical stability of the condensate. Since the bosons are strongly coupled to the fermions a strong renormalization of the boson-boson interaction has to be expected which may enhance or suppress the stability of Bose-Einstein condensation. While fermionic self-energy corrections are expected to play a subdominant role in the limit of strong population imbalance n F n B , for a study of the phase diagram away from this limit these also become an essential ingredient and may lead to striking effects such as boson-mediated pwave pairing at sufficient interaction strength [16].
The question as to which vertices (i.e. correlation functions) to include in more refined approximations of our fRG scheme is dependent on the type of phases one may expect to govern the Bose-Fermi mixture phase diagram away from the strongly-imbalanced limit -see the introductory Fig. 1. Quite generally, and similar to vari-ational techniques, in field theoretical approaches the range of phases one can discern is limited by the variety of -potentially competing-channels taken into account in the renormalization procedure. In this regard, the strongly coupled Bose-Fermi mixtures present a vast testbed to develop comprehensive theoretical approaches to competing order where a manifold of scenarios and phases may unfold, including: phase separation between the fermionic species in case of repulsive effective interactions, competing bipolaron and trion formation, boson-mediated s-or p-wave pairing of fermions, fermioninduced phonon softening that may result in supersolidity, higher-order pairing mechanisms such as bosonmediated Cooper binding of trions and phases of Efimovtype states that may condense depending on their statistics.
Moreover, as discussed in Appendix E, the formation of bound states containing several bosons may be considered. However, in ultracold quantum gases these higherbody bound states are usually subject to rapid decay to deeply bound states. The competition between such dissipative multi-particle losses and the formation of manybody phases is an intriguing perspective for future studies, posing a significant theoretical challenge that requires extension beyond equilibrium theory.
Another compelling question is what the impact of Coulomb interactions between the fermionic degrees of freedom may be. These long-range interactions will ultimately impose limits on the universal connection between strongly coupled Bose-Fermi mixtures in atomically thin semiconductors and ultracold atoms (see Table II). Coulomb interaction can be expected to play a key role in particular at low doping where screening becomes increasingly ineffective. Taking Coulomb interactions into account may indeed suppress the formation of well-defined electronic and molecular Fermi surfaces and instead lead to qualitatively different physics even in the limit of extreme population imbalance B / F , where understanding the interplay of Coulomb interaction, favoring Wigner crystallization, and boson-mediated Fermi-Fermi interactions, remains an open challenge.
Considering the myriad of open questions, the full exploration of the phase diagram of two-dimensional Bose-Fermi mixtures remains a formidable task. Due to the strong-coupling nature of the problem, uncovering the possible in-and out-of-equilibrium phases and phenomena will ultimately require a concerted effort between theory and experiment. Starting from limiting cases, such as considered in this work, that can be controllably understood and combining ab initio approaches with experimental observations will be key to tackle this outstanding challenge and can lead to new insight into effective descriptions of strongly coupled many-body quantum systems.

Gradient expansion parameters
The parameters of the gradient expansion of the twopoint functions G −1 φ,ψ,t,k are given by This prescription can then be used in Eqs. (15) to (18) and (E2) to determine the flow of the parameters. In Appendix H we provide the explicit form of the flow equations.

Appendix B: Initial conditions of the flow
The initial conditions for the flow are obtained by setting Γ k=Λ = S + const. This implies that m ψ,k=Λ = 0, m φ,k=Λ = −µ φ , m t,k=Λ = m t and A ψ,k=Λ = A φ,k=Λ = A t,k=Λ = 1. The initial conditions for the interaction vertices are given by λ k=Λ = 0 and h k=Λ = h.
The initial condition for the detuning of the molecule m t is obtained from the physical renormalization condition that in the two-body problem a bound state of energy B forms between the boson and the fermion species. Within an fRG approach, this condition is ensured by, first, setting µ ψ and µ φ to sufficiently negative values so that either species has a zero density. Furthermore, we set µ ψ + µ φ = − B . This ensures that the energetic cost to bring up a particle of both species from vacuum is given by the binding energy. Finally, the condition for the bound-state formation is given by m t,k=0 = 0 and m t,k ≥ 0. This condition guarantees that the chemical potentials are tuned correctly to the boundary between the vacuum state and the state comprised of a molecule submersed in vacuum.
The flow of the three-body vertex does not have to be taken into account, since λ k does not feed back into the solution of the two-body problem. Similarly, in the twobody problem the flow equations of G −1 ψ,k (P ) and G −1 φ,k (P ) evaluate to zero because the poles of their propagators in Eqs. (15) and (16) lie in the same half of the complex plane and thus their frequency contour integrals evaluate to zero. Physically, this is because neither of the particle species has a finite density which would be required to generate a renormalization of the particle selfenergies by particle-hole fluctuations. The molecule on the other hand is renormalized by a particle-particle diagram and thus does not require a finite density of bosons or fermions.
As the three body-vertex is not relevant in the twobody problem, the Yukawa term h does not renormal- ize. After evaluation and projection of Eq. (17), the flow equations of m t,k and A t,k are therefore given by and Using m t,k=Λ = m t,k=0 + Λ 0 dk(∂ k m t,k ) and m t,k=0 = 0 this reproduces Eq. (3). Note that since in this few-body calculation no Fermi surfaces are present, the regulators are proportional to Θ(|p| − k) for all particles involved.

Appendix C: The polaron energy within the gradient expansion scheme
In this appendix we discuss the weak non-monotonous behavior of the polaron energy as a function of F / B . Generally, the polaron energy lies approximately within a range of ± F around the value of − B . For small values of F / B it is thus not surprising to see that E → − B . Previous calculations [55][56][57][59][60][61] indicate that for all values of F / B the value of the polaron energy should lie below − B , in disagreement with the results shown in Fig. 5. This discrepancy highlights one of the major shortcomings of fRG, namely the dependence on regulators and on the truncation scheme.
To analyze this finding in detail in Fig. 11 we show the polaron and the molecule energies using different truncation and regulator schemes. As one can see, the Γ 2,k truncation (solid line) presented also in Fig. 5 results in polaron energies above − B . If, however, the same truncation is used and the regulators are changed such that the renormalization group flow consists of two steps, where in the first step only the molecule and in the second step only the minority particle is allowed to flow, this results (dash-dotted) in polaron energies strictly below − B . Within this scheme, however, the resulting molecule energy lies higher than before. Effectively, by treating the molecule and polaron on different footing (i.e. by treating them in different steps of the fRG) we have improved the polaron energy at the cost of a higher-lying molecular energy. Interestingly, this twostep calculation is closely related to the results obtained within the variational approach in Ref. [55] and the ladder resummations performed in Ref. [57] (crosses). If the full frequency-and momentum-resolved T -matrix in these two approaches is replaced by a gradient expansion of the T -matrix, the resulting method is equivalent to the two-step fRG. The results for this modified variational/diagrammatic calculation are shown as dots and coincide with the two-step calculation as expected. A similar equivalence of the FMR scheme is discussed in Appendix G 2.
Within the Γ 3,k calculation presented in Fig. 11 (dashed) and also in Fig. 5 the polaron energy lies again above − B for small F / B . If, however, the flow of h k is turned off (dotted) the polaron energy lies once again strictly below − B . These observations illustrate the dependence of the absolute values of the energy on the regulators and truncation employed. For example, the flow of the Yukawa vertex h k has a significant impact on the polaron energy, which is likely due to its point-like projection.
Although the relative deviations of these energies are only of the order of a few percent, we do not expect that the used fRG schemes are a reliable method of determining the absolute energy of the polaron and the molecule. Most of the variational approaches, however, do not consider the polaron and the molecule on an equal footing and therefore can produce ambiguous results when one considers transitions which depend on relative energy differences between the emergent quasiparticles. We thus believe that, by treating the polaron and the molecule on equal footing within a unified renormalization approach, the fRG scheme captures the qualitative physics correctly and can thus make qualitative predictions about transition in the quantum many-body system.

Appendix D: Three-body vertex in the vacuum three-body limit
To determine the value of the three-body vertex in the limit where two ψ-particles and a single φ-particle are present, we solve the flow equations under the initial conditions of the two-body problem discussed in Appendix B, and additionally take into account the flow of λ k . Since λ k corresponds to the on-mass-shell scattering of a molecule and a quasi-free excess fermion, we supplement the two-body initial conditions by setting the fermionic chemical potential to a small, negative value As in the two-body case, G −1 φ,k and G −1 ψ,k do not flow and as a result neither does h k . Consequentially, the flow of m t,k and A t,k is not influenced by the flow of λ k such that according to Eqs. (B1) and (B2) and similarly The flow equation of λ k is then given by leading to λ k=0 = −h 2 / B for large values of h.

Appendix E: Bose-Bose-Fermi coupling in the three-body limit and at finite density
The truncations considered in the main text neglect the emergence of a Bose-Bose-Fermi (BBF) coupling (and other higher-order couplings). In this Appendix we seek to explore the relevance of this coupling. From a physical standpoint, unlike the Fermi-Fermi-Bose (FFB) coupling λ k , the BBF coupling does not suffer from Pauli blocking and may thus be considerably stronger, potentially resulting in the formation of bound states containing more than one boson.
First, we study the BBF coupling in the limit where two φ bosons and a single ψ fermion are present. We define the corresponding coupling vertex τ k as Expanding Eq. (21) we obtain that ∂ k U k (σ eq ) has the following structure The integrand in the first term in Eq. (F1) does not have a pole in the frequency domain as the Θ σ,k -functions are frequency-independent. Stemming from the construction of the quantum field theory and the convergence factor of e iω0 + , this integral thus evaluates to zero. The second term, in contrast, possesses a pole within G σ,k and therefore does not vanish, yielding Eq. (22). Note that because the integrand only falls off fast enough due to the convergence factor, these ω contour integrals need to be closed within the upper half of the complex plane. Consequently, the second term in Eq. (22) always vanishes. In order to yield a finite value it would require the polaron to develop a finite density which we do not allow for within our phase identification scheme.

Appendix G: Frequency-and momentum-resolved spectral function
Here we provide the explicit flow equations used to obtain the frequency-and momentum-resolved spectral functions in Section V. Furthermore we comment on the analytical structure of these equations and how it relates to the structure of the particle-hole/particle-particle continua visible in Figs. 9 and 10 and the lifetimes of the molecule and the polaron. Finally, we show the close correspondence of this method to T -matrix approximation schemes.

Frequency-and momentum-resolved flow equations
In order to compute the frequency-and momentumresolved spectral functions within the FMR scheme, for a given value of F , µ φ and B in a first step the flow of the expansion parameters is computed as detailed in Section III and Appendix A 2. In a second step the solutions of the flow equations for the different gradient expansion parameters are plugged into the rhs. of the flow equations given in Eqs. (15) to (18) and (E2). This time, however, the flow equations are considered for arbitrary external momentum and frequency. Next, the Matsubara integration is performed as usual and the complex frequency ω of P = (p, ω) is continued to the real frequency axis iω → Ω + i0 + . For every external frequency and momentum the flow equations of the retarded Green's function can then be computed and yield

Equivalence to a non-self-consistent T -matrix resummation
In this subsection we show the close correspondence between the FMR scheme (Section V) and diagrammatic ladder approximations. More specifically, we show that a non-self-consistent implementation of the FMR method exactly corresponds to the result obtained for the molecule in non-self-consistent T -matrix resummation as presented in Ref. [57].
Using only bare quantities on the rhs. of the flow equation and performing the frequency integration in the quantum impurity limit, the flow of the retarded inverse molecule propagator reads Here we used that within this approximation λ k = 0, h k = h, A φ,k = 1, m φ,k = −µ φ , A ψ,k = 1 and m ψ,k = − F . Note that, since we only use bare quantities on the rhs., we have∂ k = ∂ k . Thus we can perform the k-integration analytically and obtain reproducing the molecular results presented in Ref. [57]. Furthermore, similar analysis shows that performing a modified non-self-consistent two-step fRG of the FMR scheme also reproduces the polaron results presented in Ref. [57]. In such a two-step approach the molecular propagator is renormalized in the first step as described in Eq. (G4), and in the second step the minority propagator is renormalized as prescribed by Eq. (15). In this second step, on the rhs. the coupling constants along with the majority propagator appear in their bare form and the molecular propagator with its full frequency-and momentum-dependence obtained in the first RG step is used instead of a gradient expansion. The polaron energy resulting from this calculation is shown as crosses in Fig. 11. It is worth noting, however, that as a starting point for the second step one may also perform a gradient expansion of the molecular propagator of the form and still obtain similar results (dash-dotted lines and dot markers in Fig. 11). This then directly corresponds to a version of the FMR scheme used to obtain spectral functions in which the renormalization of the molecule and the minority is divided into two consecutive steps while retaining the gradient expansion on the rhs. of the flow equations.

Analytical structure of the FMR flow equations
In the following we analyze the analytical structure of the flow of the retarded inverse Green's function of the molecule and how it is reflected in the spectral functions shown in Fig. 10. A similar analysis can be performed on the retarded inverse Green's function of the polaron as well.
Within the FMR scheme of analytic continuation, the retarded self-energy can only pick up a non-vanishing imaginary part in the limit of i0 + → 0 if during the flow one integrates over a pole caused by iω → Ω. In that case we have encountered a pole in the flow that is only avoided by the use of a retarded frequency and the selfenergy picks up an imaginary part that is non-vanishing for all i0 + .
Contributions to the spectral function defined in Eq. (28) can have two different origins. Either the Green's function picks up an imaginary part in the course of the flow as described above, or the inverse Green's function tends to i0 + resulting in a sharp excitation feature in the spectral function. In the former case the corresponding states are part of a particle-particle continuum of states with a finite lifetime, whereas in the latter case the corresponding excitations have an infinite lifetime. Inspecting the second term of Eq. (G2), we see that it causes the self-energy to develop an imaginary part if during the flow while p 2 +2|p|k cos(θ)− F > 0. For p 2 < F the minimal frequency for which this can occur is given by where we made use of the fact that m φ,k /A φ,k decreases monotonically during the flow. In turn, for p 2 > F this frequency is given by Ω > min,1 = min k,0≤k≤ (G8) Analogously, the minimal frequency for which the third term of Eq. (G2) leads to an imaginary part is given by Numerically we find that this is solved by for the interaction strengths studied here.
In Fig. 13, the spectral function from Fig. 10(f) is shown along with the minimal frequencies Ω min,1 and Ω min,2 . As it can be seen, these frequencies determine the onset of the particle-particle continua. Furthermore, as the molecule peak at low and high momenta lies outside the boundaries of the continua, the corresponding excitations possess an infinite lifetime within this renormalization scheme.