Polaron spectroscopy of interacting Fermi systems: insights from exact diagonalization

Immersing a mobile impurity into a many-body quantum system represents a theoretically intriguing and experimentally effective way of probing its properties.In this work, we study the polaron spectral function in various environments, within the framework of Fermi-Hubbard models. Inspired by possible realizations in cold atoms and semiconductor heterostructures, we consider different configurations for the background Fermi gas, including charge density waves, multiple Fermi seas and pair superfluids. While our calculations are performed using an exact-diagonalization approach, hence limiting our analysis to systems of few interacting Fermi particles, we identify robust spectral features supported by theoretical results. Our work provides a benchmark for computations based on mean-field approaches and reveal surprising features of polaron spectra, inspiring new theoretical investigations.


Introduction
A mobile impurity immersed in a many-body background represents a paradigmatic setting of many-body physics.Historically, Landau and Pekar were the first to discuss the renormalized mobility of an electron interacting with the phonons of a polar crystal, and they coined the emerging quasi-particle polaron [1].
Following the achievements in the precise preparation and control of ultracold atomic mixtures, there has been a renewed interest in studying polarons at a fundamental level [2][3][4][5].Indeed, immersing a single mobile impurity in a non-interacting or weakly-interacting background already represents a fully fledged many-body problem.So far, ultracold atom experiments have been able to probe the attractive and repulsive branches of standard, weakly interacting Fermi [6] or Bose [7,8] polarons, using radio-frequency spectroscopy.Such measurements are complicated by the metastability of atomic gases under three-body recombination and by the challenge of independently tuning the interactions within the background and those involving the impurity, so that polaron studies have yet to be extended to the case of strongly correlated backgrounds.It is nonetheless worth mentioning recent measurements of the so-called "magnetic polaron" in doped Hubbard anti-ferromagnets [9], of Nagaoka polarons arising from kinetic magnetism in triangular lattices [10,11], and of the dressed, quantum-statistics-dependent interactions between polarons [12].
In parallel, the last years have witnessed the thriving of 2D materials, such as graphene and transition metal dichalchogenide (TMD) semiconductors [13].Contrary to cold atoms, in TMD heterostructures experimentalists have access to few observables only, and optical polaron spectroscopy [14], where resonantly injected excitons are used to probe the state of the TMD material, is now daily used.More precisely, this method relies on the exciton being dressed by the electronic excitations, with the exciton-electron scattering properties being determined by the binding energy of the trion, which is the bound state of an exciton with an electron.Remarkably, TMD bilayers hosting moiré potentials have emerged as an ideal platform to simulate Fermi-Hubbard and extended Fermi-Hubbard physics [15,16].Transport or optical signatures of several exotic phases of matter have already been reported, including Wigner crystals [17,18], charge density waves [19,20], excitonic insulators [21][22][23], superconductivity [24,25] and anomalous quantum Hall states [26].However, the very large trion binding energy with respect to the moiré-electronic scales may have limited the number of features observable in polaron spectra, as we will further illustrate in this work.
Given the hardness of the many-body problem, polarons in strongly interacting many-body backgrounds have been theoretically considered in a few special settings only or under very rough approximations.We mention attempts in the context of fractional Chern insulators [27,28], Fermi superfluids along the BEC-BCS crossover [29,30], excitonic insulators [31], Mott and charge transfer insulators [32], kinetic magnetism [33,34], Holstein polarons in Luttinger liquids [35], together with a few related works on the dressing of optical excitations in Mott insulators [36,37] and fractional quantum Hall systems [38].Diffusion Quantum Monte Carlo can be used to study Bose [39] and Fermi polarons [40] in the intermediate correlation regime, but does not allow for the determination of the full spectral information.
In this work, we use exact diagonalization (ED) to tackle this challenging problem, exploring various configurations of the interacting Fermi gas.Specifically, we compute zeromomentum polaron spectra in lattice systems of a few interacting fermions, both in the spinless and in the spinful case, with different classes of repulsive or attractive interactions.The main drawback of ED is that this method is limited to very small system sizes.In spite of these finite-size effects, our results are qualitatively consistent when comparing between different lattices or, when possible, to other methods (e.g. using the Chevy Ansatz within a mean-field theory for the many-fermion background).In this sense, we believe that our ED approach provides solid results, which allow to benchmark approximate methods but also provide novel insights.Moreover, studying how the growth of many-body correlations scale with the system size is an interesting topic on its own [41,42].
This manuscript is structured as follows: in Section 2, we illustrate the ED method and benchmark it for an impurity immersed in a non-interacting Fermi sea.In Section 3, we discuss polaron spectra in the presence of strong charge density wave correlations for a system of spinless fermions with long range repulsion.In Section 4, we tackle the spinful case with a spin-dependent impurity-fermion interaction, and we study the fate of two-body and threebody absorption lines in the presence of contact repulsive interactions between the fermions.In Section 5, we consider polaron spectroscopy of fermionic superfluids, in the attractive Fermi-Hubbard model, both for balanced and nearly-balanced populations of the two species.Finally, we draw our conclusions and outline future directions in Section 6.

General framework
This work explores the problem of an impurity immersed in a system of spin 1/2 fermions, described by a Hamiltonian of the form Here H f f denotes the interacting fermionic background, described by a generalized Fermi-Hubbard Hamiltonian where c † jσ denotes the creation operator of a fermion of spin σ ∈ {↑, ↓} at lattice site j = (x, y), t f is the hopping constant and V i j denotes the matrix elements of a general (possibly long range) interaction between the fermions.
The kinetic energy for free impurities is H I = −t im 〈i, j〉 a † i a j , where a † j denotes the creation operator of the impurity and t im its hopping rate.In this work, we will limit ourselves to contact impurity-fermion interactions Notice that the coupling constant U σ depends on the spin of the fermion.This is a natural choice both in atomic systems, where Feshbach resonances are spin-dependent [43], and in solid-state spectroscopy, where the trion binding energy depends on the polarization of the probe excitons [44][45][46] (e.g. in molybdenum based TMD monolayers, spin-triplet trions are typically unbound).
Describing a Fermi polaron in a non-interacting Fermi sea is already a fully fledged manybody problem, which cannot be solved exactly [2][3][4].Here we commit ourselves to the study of models of interacting fermions, such that accessing the ground state | f 0 〉 of H f f already represents a formidable task.We address this problem using ED, which consists in expressing the Hamiltonian as a sparse matrix in the real space basis [47][48][49][50], and in using the Lanczos method [51] to determine the ground state.
ED gives also access to the ground state |Ψ〉 of the full impurity-fermion Hamiltonian H.However, the most natural observable in both ultracold atom [52] and solid-state experiments [14] is provided by the spectral function of the impurity.Our main goal is thus to compute polaron spectra, which experimentally corresponds to the following protocol: first, the many-fermion system is prepared in its ground state in the absence of the impurity, or the fermion-impurity interaction is switched off by preparing the impurity in some hyperfine level; then, the impurity is resonantly injected, or its internal state is resonantly flipped to a hyperfine level characterized by a sizable fermion-impurity interaction; the final observable is the frequency-resolved absorption curve of the resonant excitation.Mathematically, the polaron spectral function is defined as where |Ψ 0 〉 is the ground state of the fermion-impurity decoupled Hamiltonian H f f + H I , with energy E 0 , while H is the full Hamiltonian with finite fermion-impurity coupling.Since generally |Ψ 0 〉 = a † k=0 | f 0 〉, one can make a link with optical spectroscopy in TMDs, where an exciton is injected with very small momentum.The spectral function is also related via Fourier transform to the overlap S(t) = 〈Ψ 0 |e −i(H−E 0 )t |Ψ 0 〉 following a quench of the impurity-fermion interaction.Notice that the spectral function (4) satisfies the normalization In plotting the spectra, the lines are artificially broadened by replacing ω → ω + iγ.
A crucial technical remark concerns the fact that, even though the full spectrum of H cannot be obtained for large sparse matrices, the spectral function can still be reliably and efficiently obtained [48,51,53].The trick consists in constructing the Krylov space of dimension M K by applying M K times the full Hamiltonian H to the decoupled ground state |Ψ 0 〉.The Hamiltonian in this space is represented by a tridiagonal matrix, for which the resolvent of the first entry can be conveniently computed by recursion and expressed as a continued fraction.It can be proven that this approach captures exactly the first 2M K + 1 moments of the spectral function.
In this work, we restrict ourselves to nearest-neighbour hopping without complex phases and we consider two-dimensional square or triangular lattices.We will be limited to very small system sizes such that finite-size effects represent a major concern; however, in the following, we will show and argue that one can still extract solid and useful qualitative insights using this approach.In our philosophy, this method naturally needs to be complemented with other approaches, such as Chevy Ansatz computations [54,55].For instance, Section 5 focuses on a polaron in a Fermi superfluid, which is a good example of how ED can confirm a highly nontrivial feature previously reported using mean-field calculations complemented by the Chevy Ansatz [30,31], an approach reviewed in Appendix D.
An insidious consequence of finite-size effects is that the ground states |Ψ 0 〉 and/or |Ψ〉 may not be invariant under the action of the spatial symmetries of the model.For instance, the ground state of the non-interacting system with 2 identical fermions (i.e.N ↑ = 2, N ↓ = 0, V i j = 0) consists of a superposition of states in the form c † k=0↑ c † k̸ =0↑ |0〉, which is not a zero momentum eigenstate of the total momentum.In the following, we restrict ourselves to the sector of zero total momentum and require that |Ψ 0 〉 = a † k=0 | f 0 〉, i.e. that the fermion and impurity momentum before quenching the impurity-fermion interaction are separately zero: this imposes a strong constraint on the number of particles that one can fit into the system.In the case of crystalline phases, one should also respect the commensurability of the crystal with respect to the size of the system.

Benchmarking Fermi polarons
In this subsection, we consider the simplest scenario of an impurity immersed in a noninteracting Fermi sea of spin-polarized fermions (N ↓ = 0 for definiteness).
The fermion-impurity binding energy in the vacuum E B > 0 is related to the coupling U ↑ appearing in H f I by the Lippmann-Schwinger equation [56] 1 with L x , L y the number of lattice sites in the two independent spatial directions and ε f k , ε I k denoting the free fermion and impurity dispersions, respectively.We define the Fermi energy through the formula E F = 2πħ h2 m f n, which holds in continuous infinite systems, where n = N ↑ /A is the density.Here we defined the area of the system A, where, assuming unit distance between adjacent lattice sites, A = L x L y for a square lattice and A = 32 L x L y for triangular lattices; furthermore, for the effective fermion mass m f , one can use the correspondence ħ h 2 2m f = t f on a square lattice and ħ h 2 2m f = 3 the case of a 6 × 6 triangular lattice with N ↑ = 4 (right panel).The red lines mark the vacuum impurity-fermion binding energy −E B , while the magenta AP and RP labels indicate the attractive and repulsive polaron branches, respectively.In these plots we used a small linewidth γ = 0.1E F and a logarithmic color scale to highlight the fine structure of the molecule-hole continuum, which depends on the details of the computation.We devote Appendix A to a more detailed comparison between different lattice sizes and number of particles, as well as to showing convergence in the Krylov space dimension M K .

Charge density waves
In this section, we will be interested in the physics of fermions with repulsive finite-range interactions.For definiteness, we will consider spin-polarized systems with Coulomb repulsion , where V 0 is the coupling constant and r i denotes the position of the i-th lattice site.Given the underlying lattice structure and the finite filling factors used in the following, we expect our results to qualitatively hold for generic finite-range repulsive potentials.
Polaron spectroscopy has already been used in experiments to detect Wigner crystals in TMD monolayers [17] and bilayers [18] without a sizable moiré potential.The main signature of the presence of the Wigner crystal is a weak umklapp line, which departs from the repulsive polaron peak with a splitting that scales with doping1 like ∝ , with m X the mass of the exciton and a M the lattice constant of the Wigner crystal.These experiments can be modeled assuming that the Wigner crystal is a static external potential that scatters the probe exciton.
While Wigner crystals are formed in continuum systems and spontaneously break a continuous translational symmetry, our ED approach is limited to finite-size lattice models, and we will refer to states with strong crystalline fluctuations as "charge density waves".Because of the finite size, the discrete translational symmetry cannot be spontaneously broken and one needs to inspect the density-density correlator 〈 f 0 |n i n j | f 0 〉, with n i = c † i↑ c i↑ , in order to monitor the crossover from a Fermi liquid to a charge density wave.
Impurities in lattice systems in the presence of strong repulsion have been considered as well in the TMD literature.For instance, the umklapp peak was used to reveal the onset of charge incompressibility of correlated electrons in moiré bilayers [59], even without the breaking of the discrete translational symmetry.The optical signatures of generalized trions in a moiré system displaying charge density waves at fractional filling were very recently reported in [60].
We computed the polaron spectrum for a 6 × 6 square lattice with N ↑ = 4 fermions, a number commensurate with the formation of a crystal of lattice constant a C DW = 3a 0 .We fix the strength of fermion-impurity interactions according to E B /E F ≃ 2.2 and scan a broad range of values of V 0 ; as it is usually done for Wigner crystals, we express the strength of the interactions as κ = V 0 a C DW E F , i.e. the ratio of the typical interaction to kinetic energy.The results are shown in Fig. 2. In panel (a) the polaron spectrum is depicted as a function of κ.The color mesh is in linear scale, while in panel (b) we plot exactly the same data in log scale.In panels (c) and (d), computed along the vertical dotted slices of panel (a), we display the density-density correlator 〈 f 0 |n (0,0) n (x, y) | f 0 〉 at small and large κ, respectively.These two plots illustrate the building up of crystalline correlations along the crossover from Fermi liquid to charge density wave.
Let us now focus on the features visible in panels (a) and (b).To start with, the attractive polaron (AP) experiences a redshift for increasing κ.We attribute this feature to the renor- malization of the fermion mass, resulting in a larger binding energy; in particular, in the limit V 0 → ∞ all the fermions are perfectly correlated, so the impurity binds to an object of total mass N ↑ m f .Then, we have highlighted by the pink dashed line the umklapp peak [17], well defined at large V 0 .This can be thought of as a scattering state of the impurity with the fermionic crystal and occurs at an energy essentially determined by the folding of the free impurity dispersion in the first Brillouin zone of the crystal, i.e. in a square lattice at wavevector 2π/a C DW .In Appendix B we have indeed verified that this peak behaves like expected.
Finally, the most intriguing and original feature is the avoided crossing occurring for comparable values of the binding energy E B and the effective interaction strength V 0 /a C DW .Here we speculate on the origin of this effect, which could be further studied and confirmed in the future using the Chevy Ansatz.In a Hartree-Fock picture, the formation of a charge density wave is described in terms of Bloch waves on top of the Hartree-Fock potential, with the same periodicity as the charge density wave and determined self-consistently by filling the lowest Bloch band.When V 0 is increased, the gap between the lowest and second band scales ap- proximatively linearly; moreover, the low-lying bands also become flatter and flatter.It seems then likely that the repulsive polaron (RP) hybridizes with the state made of the bound impurity plus an excitation of the neutral mode corresponding to the interband excitation of one Hartree-Fock quasi-particle; the dashed-dotted line in panel (b) is a tentative guide to the eye for the energy of this candidate state, which blueshifts linearly with V 0 .Notice that a similar feature was observed for a homobilayer model with local repulsion, in a ladder computation using the electron Green's function computed in DMFT, and was attributed to doublons [32].An even more characteristic feature is observed on a triangular lattice, as we show in Fig. 3(a).Understanding this peculiar behavior is also the subject of ongoing efforts.
A piece of evidence that these results are not just an artifact of finite-size effects comes from the study of one-dimensional chains, for which it is possible to reach a reasonably large number of sites L = 22 with N ↑ = 11 fermions.In 1D, we define the Fermi energy as E F = 2πt f n 2 .As one can see in Fig. 3(b), the avoided crossing is also found in such a long chain.Interestingly, no umklapp peak is found in this 1D case.
To conclude this section, we comment on the relation of our results with existing experiments in TMDs.The avoided crossing, in particular, has not been observed in the charge density waves and Wigner crystals reported so far in the literature.This is due to the fact that the trion Bohr radius is much smaller than a moiré cell, or in other words, the trion binding energy is large compared to the electronic many-body gap.The avoided crossing occurs when the two energies are comparable, and this could be made possible by engineering trions with a smaller binding energy in heterostructures with strong moiré potential.Such trions may be accessible in multilayer systems, with the excitons and the fermions mainly localized in two different layers.

Impurity interacting with two Fermi seas
While polaron spectra of molybdenum-based TMDs are essentially consistent with the Fermi polaron picture [14,61], experiments in WSe 2 monolayers display a much richer structure [44,62].At small electron doping, two attractive polaron lines are visible together with the repulsive polaron; at high densities, instead, all the oscillator strength is taken by a low energy line, which redshifts with increasing density.This difference originates from the fact that two identical electrons and a hole will generally not bind into a trion, as a consequence of the anti-symmetrization of the wavefunction [44][45][46].In MoS 2 and MoSe 2 monolayers exciton formation involves the lowest conduction bands of quantum numbers (K, ↑) and (−K, ↓), so that right circularly polarized radiation can excite only the (K, ↑) transition and the only visible trion will be the one formed from this exciton and an electron in (−K, ↓).Here, ±K denote the corners of the first Brillouin zone and the arrows the physical out-of-plane spins.In WSe 2 monolayers, instead, the bright conduction bands (K, ↑) and (−K, ↓) lie higher in energy than the dark (K, ↓) and (−K, ↑) lower conduction bands.This means that the dark conduction bands get doped, while the electron forming an exciton comes from the bright conduction bands and can form a trion with either the (K, ↓) and (−K, ↑) electron, which is distinguishable by valley or spin.
This argument explains qualitatively the presence of two attractive polaron resonances.Some theoretical understanding of the crossover to a single line red-shifting with doping has been provided in [62,63] based on variational wavefunction calculations.The first step is to realize that the attractive polaron is not really associated with a trion (exciton bound to an electron), but rather with a tetron, i.e. an exciton bound to an electron-pair formed out of the conduction band Fermi sea.At small doping, the phase space of the hole is limited to very small momenta and the tetron is basically a trion very loosely bound to a hole.In practice, to locate the two attractive polaron lines at small doping, one can compute the binding energy of the singlet and triplet trion.However, when the exciton can bind to two Fermi seas and at large doping, a so-called "hexciton" can be formed from the exciton and a particle-hole pair out of each of the two Fermi seas.The high doping condition ensures that the particle-hole pair is small and akin to a neutral object; otherwise, if the holes were highly delocalized, the exciton could not bind to two charged electrons.
This rich physics can, in principle, also be accessed in ultracold atom setups with three species [64,65]; however, to our knowledge, no experiment has ever investigated the spectral properties of a polaron in two Fermi seas to this date.
Here, we apply our ED method to an impurity immersed in a background with two species of fermions.For the sake of dealing with only one parameter to tune the strength of repulsion, we consider the case of contact repulsive interactions with coupling V i j = V ↑↓ δ i j .We have checked that one obtains similar results in the presence of Coulomb terms, both inter-and intra-spin (contact interactions, instead, are only effective for fermions of different spin, due to the Pauli principle).To observe two attractive polaron peaks, we allow for spin imbalance of the contact fermion-impurity attraction U ↑ ̸ = U ↓ .In order to plot our results, we define the and the binding energy scale E B from the energy of the molecule2 formed by the impurity and a fermion with coupling U, according to 1 . Since we cannot continuously change the density, we change the ratio E F /E B by varying U.
Exact diagonalization results are plotted in Fig. 4 for a system of N ↑ = N ↓ = 4 fermions in a 4×4 triangular lattice (plots for the square lattice are very similar, see Appendix C).We choose a small spin asymmetry of 3 for the fermion-impurity interaction and we increase the fermion-fermion contact repulsion V ↑↓ through panels (a-c).Apart from the repulsive polaron peak, one can spot two attractive polaron resonances, denoted respectively AP ↑ and AP ↓ in panel (a), which are associated with spin up and spin down two-body bound states and which are predominant at small E F /E B ; and a lower energy resonance AP ↑↓ , corresponding to the three-body bound state of the impurity with both the spin up and spin down particles.The identification of these lines is supported by the calculation of the two-and three-body binding energies in vacuum, which recover the many-body lines in the limit of small E F /E B (cyan lines in panel (a)).
In the absence of the repulsive interactions between the fermions, V ↑↓ = 0, all these three lines are simultaneously visible, with a slow oscillator strength transfer to the three body resonance for increasing E F /E B .We explain this effect by noting that, for the three-body line, the oscillator strength should scale with the square of the density, while for the usual two-body attractive polarons it scales linearly in the density.In other words, at large E B the three-body state has a very small radius and its wavefunction is very different from the ground state without impurity.
When the repulsive interaction V ↑↓ is turned on, there seems to be an avoided crossing between the lowest two-body line AP ↑ and the three-body line AP ↑↓ .This comes from the fact that, at first order in perturbation theory, the three-body state blueshifts linearly with V ↑↓ , since the spin up and down fermions are bound together by the interaction mediated via the impurity.Interestingly, panel (b) reminds of the experimental data from WSe 2 monolayers [62].The splitting of the avoided crossing decreases with V ↑↓ , see panel (c).Notice that, since the vertical axis is in units of E B , the redshift of the lines with E F /E B is due to the well-known shift of the attractive polaron with density, which is linear for small E F /E B [66].
Our method is somehow complementary to the variational computation of [62,63], which is effectively a few-body computation, where the many-body background is traced out and enters via an effective static screening potential.Also, the variational method is well suited to compute the ground state energy and oscillator strength, but it is difficult to get information about the full spectrum using that method.This makes our ED approach particularly interesting to study the crossover between the two-body and three-body lines as a function of repulsive interactions.

Fermi superfluids
In this section, we will be dealing with polaron spectroscopy of Fermi superfluids, with the pairing of spin up and down fermions.The fermionic Hamiltonian we consider is essentially an attractive Hubbard model, with contact inter-spin interactions V i j = δ i j V ↑↓ , where V ↑↓ < 0 is related to the binding energy E pair of a fermionic pair in vacuum via the Lippmann-Schwinger equation 1 As to the fermion-impurity interaction, we allow for spin-dependent interaction U ↑ ̸ = U ↓ .As in the previous section, we define U = and the binding energy E B from the energy of the molecule formed by the impurity and a fermion with coupling U.In other words, E B and U are related by 1 U = 1 The motivation to study this setting comes from future experiments in both ultracold atoms and solid-state devices.On the ultracold atom side, the tunability of interactions via Feshbach resonances has allowed to observe the BEC-BCS crossover [67].More recently, there has been some effort in implementing three-species Fermi mixtures [64,65], also inspired by analogies with the SU(3) group relevant in high-energy physics for the theory of the strong interaction.Therefore, polaron experiments in Fermi superfluids will hopefully be realized in the near future.On the theory side, a Chevy Ansatz study of polaron formation in 3D Fermi superfluids with spin-symmetric U ↑ = U ↓ interactions has been performed in [29], while in the fully spin asymmetric case U ↓ = 0 full spectra were reported in [30], in both 2D and 3D.
On the solid-state side, instead, we were particularly inspired by the possibility of optically probing the presence of pairing in putative excitonic insulator states found in recent transport experiments in TMD heterostructures [21][22][23]68].In this case, a 2D Chevy Ansatz study was performed in [31], where the two fermionic species correspond to two different layers and fully pseudo-spin asymmetric interactions were considered.Remarkably, unbalanced electron-hole mixtures were very recently demonstrated [68][69][70].

Spin-balanced fermion populations
We used our ED approach to compute polaron spectra in a Fermi superfluid background.The results are reported in Fig. 5 for a 4 × 4 triangular lattice with N ↑ = N ↓ = 4.As detailed in Appendix C, similar results were obtained for for N ↑ = N ↓ = 3 fermions on a 5×5 triangular or square lattice.In Fig. 5, U σ is different but fixed in each panel, while E pair is scanned.Panels (a-c) correspond respectively to In other words, the interaction is fully spin asymmetric in (a), while it is symmetric in (c), and panel (b) lies in between.The solid cyan line, labeled E 3 , represents the energy of the three-body bound state in vacuum, while the dashed-dotted line E * 3 stands for the first three-body excited state.The lowest line can be identified as the attractive polaron corresponding to the 3-body bound state.On the BCS side it is mainly the interaction with the impurity to provide the attraction, while one rather has a Fermi pair bound to the impurity on the BEC side.This follows from the fact that in the small E pair limit one basically has a Fermi polaron, while for large E pair a Bose polaron description is adequate.Unexpectedly, in the spectra this occurs not just as a redshift of a single line, but for large E pair /E F the attractive polaron line broadens and develops a secondary peak on the high energy side.Also, there is a certain sensitivity to the finite size of the system, since the secondary peak appears on the lower energy side in the instances of Appendix C. Interestingly, this feature is completely missing in the Chevy Ansatz calculations detailed in Appendix D and it will be interesting to further investigate this in the future.We speculate that it may be due to the involvement of Goldstone-like excitations, which are not captured in the Chevy approach.Another explanation may be provided by some sort of molecule-polaron transition [71].A more precise characterization is left as an open problem, but this phenomenon will constitute a good motivation and benchmark for improving on the current Chevy framework.
The other non-trivial feature visible in the spectrum of panel (a) is the avoided crossing on top of the repulsive branch for E pair ∼ E F .This is perfectly consistent with the Chevy Ansatz prediction of [30,31], which we adapt in Appendix D to our lattice setting.The wavefunction of the state shifting with E pair was shown to have 2s symmetry, suggesting its relation to the Higgs mode of the superfluid on the BCS side and to the pair 2s excited state on the BEC side [30,31].We remark that it makes a very strong argument that two completely different methods (i.e.Chevy Ansatz and ED) show the same feature, since the former is obtained by a mean-field BCS theory approximation followed by a variational restriction to the space of one quasi-particle pair excitations, while the latter is an exact method suffering from finite-size effects.

Spin-unbalanced fermion populations
Within ED one can also consider the case N ↑ ̸ = N ↓ , i.e. a homogeneous mixture with spinunbalanced fermion numbers.This may be achieved in both ultracold mixtures and excitonic insulator setups, even though instabilities towards Fulde-Ferrell-Larkin-Ovchinnikov states or phase separation may strongly limit the available parameter region [72,73].
In Fig. 6, we present spectra for a 4×4 square lattice with N ↑ = 4, N ↓ = 3.We fix E B ≃ 0.5E F (where the Fermi energy is computed for spin up) and vary E pair .Panels (a-c) correspond to different fermion-impurity interactions, namely We hereby highlight some interesting features.First of all, the three-body attractive polaron line redshifts and slightly broadens for increasing E pair /E F , in analogy with the results shown in Fig. 5. Second, analogously to Fig. 5, an avoided crossing is well visible on the repulsive line in the panel (a) or (d) corresponding to U ↓ = 0 or U ↑ = 0, respectively.Finally, in the BEC limit one has three tightly bound fermionic pairs plus one spin up unpaired fermion (for this choice of particle numbers N ↑ = 4, N ↓ = 3).Then, for U ↑ < 0 the impurity can bind to either a pair or the unpaired spin up fermion, so that the attractive lines labeled AP ↑ are visible in panels (a-c), where we have also added the cyan dotted boxes as a guide for the eye.For U ↑ = 0, instead, the impurity can bind only to a pair, giving rise to only one attractive line, see panel (d).

Conclusion
To summarize, we have adapted the exact diagonalization method to the computation of polaron spectra on top of different many-body backgrounds, described by extended Fermi-Hubbard models.Varying the range and the sign of the interactions between the fermions and the polarization of the system, we have predicted that charge density waves, multiple Fermi seas and Fermi superfluids display polaron spectra with distinctive features.This supports polaron spectroscopy as an effective way of probing many-body systems.
As already mentioned, cold atoms offer a promising platform for the exploration of the results presented in this work.Mobile impurities can be injected in atomic quantum gases, by exploiting controlled transfer between different atomic internal states [4].In this context, both inter-species and intra-species interactions can be tuned by means of Feshbach resonances [43,[74][75][76]; atomic mixtures with controllable imbalance can be prepared [77]; lattices of various geometries can be designed [78]; and polaron spectra can be obtained through RF spectroscopy [4,79] and (momentum-resolved) Raman spectroscopy [80].By manipulating small ensembles of atoms [41,42], these systems could explore the extrapolation between the ED results presented in this work and their thermodynamic limit.Optical lattices can also be engineered, even though, as we have shown in Appendix C, the presence of a lattice is not an essential physical ingredient to obtain our results (with the exception of Section 3, the lattice just makes it technically possible to perform ED).Possible limitations of these systems concern the contact nature of the inter-particle interactions and the existence of three-body recombination, which limits the measurement time and can be kept into account only phenomenologically in our theoretical calculations, via the width of the lines.An interesting perspective concerns polarons in dipolar gases exhibiting finite-range (dipole-dipole) interactions [39,81], as well as long-range interactions mediated by an additional Bose gas, which was recently proposed as a strategy to realize crystalline states [82].
Regarding solid-state realizations, intriguing quantum many-body phases, including the ones treated in this paper, have already been realized in moiré TMD heterostructures [15][16][17][18][19][20][21][22][23][24][25][26], and optical spectroscopy is a rather straightforward and well-established technique in this context [14].The main limitation of polaron spectroscopy in current platforms concerns the fact that the binding energy of the trion is not easily tunable, but is generally of the order of 20-30 meV, hence larger than the typical electronic many-body gaps of few meV's.In particular, this explains why the avoided crossing in the repulsive branch predicted in Section 3 was not reported in the Wigner crystal experiments of Refs.[17,19,59].However, we expect that trions with reduced binding energy may become accessible in the coming years, e.g. by a careful engineering of screening or building on the bilayer Feschbach resonance [83].
On the technical side, ED is limited to very small system sizes.Despite this limitation, comparisons with other approximate methods or experiments (whenever possible) suggest that our ED approach yields qualitatively solid results.An important theoretical puzzle, which is worth investigating in the future with complementary methods, concerns the description of the attractive polaron branch displayed in Fig. 5, in the case of the fermionic superfluid.Moreover, the ED results shown in Figs. 2 and 3 have pushed us towards developing a theory of polarons in charge density waves, an investigation which will be published soon.Other future efforts will be directed to computing polaron spectra in unconventional settings using the ED method, including topological baths and Rabi driven impurities [84][85][86].In particular, the attractive polaron looks like a very stable feature, while the particle-hole continuum is strongly dependent on the details of the simulation.The repulsive polaron is in between, in the sense that it is roughly the same through different numerical setups, but with some variability.As stated in the main text, our philosophy is that this moderate variability can actually be useful: if some feature (like the avoided crossing in the RP branch showed in Fig. 2 or the broadening in the AP of Fig. 5) is shared by different lattice shapes and sizes, it is most likely a physical result, because artifacts of lattice discretization and finite-size effects depend strongly on the numerical details.
In Fig. 7 we compare the spectrum of the polarized Fermi polaron at E B = E F for different numerical setups.We also show convergence with the dimension M K of the Krylov space by comparing with the spectral function obtained from the full diagonalization of the Hamiltonian.(In this particular case the Hilbert space has dimension N H = 18424 and so it is feasible, in the other scenarios investigated in the paper a full diagonalization is usually impossible.The largest Hilbert space considered in this work is for Fig. 5 and has dimension N H = 5290000.) In black dashes we also compare with a Chevy computation for a continuum, fully polarized non-interacting Fermi sea |F S〉.More specifically, the Chevy Ansatz leads to the equations of motion is a very short-range fermion-impurity interaction with a Gaussian shape of size ℓ W = 0.2k −1 F and depth W 0 chosen to yield E B = E F .The equations of motion are diagonalized on a polar grid with a momentum cutoff of k cut = 15k F (see Supplementary material of [31]).While all the spectral features have a very similar location in energy, it turns out that in ED the RP and molecule-hole continuum have a much large spectral weight.It is not clear whether this discrepancy is an ED finite-size artifact, an effect due to the microscopic lattice (so that in principle the two models are different), or, more intriguingly, a shortcoming of the Chevy approximation.For the RP branch in particular, multiple particle-hole scattering channels may be relevant.

B Zooming in on the umklapp peak
In this paragraph we report numerical evidence that the spectral feature labeled as umklapp peak in Section 3 behaves as expected.In particular, the umklapp peak energy is uniquely determined by the dispersion of the impurity and the lattice constant of the CDW, but will be (at first order) independent of the fermion-impurity interaction.In Fig. 8(a) we display cuts of the spectral function for N ↑ = 4 strongly repulsive fermions in a 6×6 triangular lattice and three different E B 's.While the attractive polaron branch is strongly affected, the location of the umklapp peak is mostly unchanged.In Fig. 8(b) we instead compare three different lattices at the same E B ≃ 1.8E F .Respectively, the single particle Bloch bands in an infinitesimal potential with the periodicity of the CDW suggests a RP-umklapp splitting of 6t I for the 6 × 6 triangular lattice, 4t I for the 4 × 4 square, and 3t I for the 6 × 6 square, compatible with the spectra.Due to these large splittings (arising from the small CDW lattice constant), the oscillator strength of the umklapp peak is quite small, e.g. its height is 5% for the 4 × 4 square case.
We have also verified that the umklapp peaks do not shift with t f or V 0 (not shown).

C Results for square lattice
As already discussed in Sec. 2 and in Appendix A the polaron spectra on top of a noninteracting Fermi sea are, to a very good approximation, independent of the microscopic lattice.In this Appendix, we show that this is also the case for the repulsive spinful Fermi fluid described in Sec. 4 and the superfluid of Sec. 5. Since in the main text we reported results on the triangular lattice, here in Figs. 9, 10 and 11 we show the analogous computations for the square lattice.As the reader can see, the agreement is excellent, the only qualitative discrepancy being for the superfluid, in the behavior of the AP line at large E pair /E F , as already mentioned in the main text Sec. 5, where in Fig. 10 there is a secondary peak highlighted by the pink arrow, at lower energy than the main AP peak.This is also the case for the triangular 5 × 5 lattice with N ↑ = N ↓ = 3 calculation reported in Fig. 12.While this finite-size effect prevents us from making final statements on the precise structure of the AP peak, we can reasonably expect that the presence of some broadening in the AP branch or of a close-by secondary peak is a physical result; this is a novel feature which deserves further investigation and the development of beyond Chevy theoretical tools.

D Chevy Ansatz on top of BCS mean-field state
For comparison with the Fermi superfluid results in Sec. 5 of the main text, in this Appendix we compute spectra using the generalized Chevy Ansatz on top of the BCS mean-field state [29][30][31].For notational simplicity, in the following we introduce the rescaled couplings g = V ↑↓ /L x L y and g σ = U σ /L x L y .The BCS state |BC S〉 is defined as the vacuum (i.e.γ σ |BC S〉 = 0) of the quasi-particles where where the last three terms describe respectively up-down, impurity-up and impurity-down scattering processes.The spectral weight is given by the |ψ 0 | 2 component of each eigenvector.
Differently from previous works, where (nonlinear) polar grids where used to achieve UV convergence, here we are directly interested in a lattice model.This automatically takes care of any UV issue, so that we can use Bravais grids for the momenta.As a consequence, the three interaction processes (up-down, impurity-up and impurity-down scattering) can be dealt with on an equal footing, and we can treat not only the U ↑ = U ↓ [29] and U ↓ = 0 [30,31] cases, but also all the intermediate ones.= 2, 1, 0, respectively.A triangular 11×11 lattice was used, with an average filling of one fourth, and impurityfermion interactions of strength given by E B = 0.5E F .This plot shows remarkable analogies with its ED analogue Fig. 5.
The results are reported in Fig. 13.There is an excellent qualitative agreement in the repulsive branch phenomenology with the ED calculation.What is instead missing are the secondary contributions to the AP line at large E pair /E F , which may suggest an origin from Goldstone fluctuations, which are neglected in the Chevy subspace.A way of taking care of these processes in the polaron problem is an open theoretical challenge, and will be the subject of future research efforts.

Figure 1 :
Figure 1: Polaron spectra on top of a spin-polarized non-interacting Fermi sea.Two ED calculations are compared, on a 7 × 7 square lattice with N ↑ = 3 fermions (left panel) and on a 6 × 6 triangular lattice with N ↑ = 4 (right panel).The red line represents −E B , while AP and RP indicate the attractive and repulsive polaron branches, respectively.The spectral function intensity is color-coded in a logarithmic scale.

Figure 2 :
Figure 2: Polaron spectroscopy of the crossover from non-interacting Fermi sea to charge density wave, for fixed E B /E F ≃ 2.2.(a) Polaron spectrum for increasing values of the interactions.Horizontal pink dashes indicate the position of the very weak umklapp peak (visible only from panel (b)).The vertical orange and red dots correspond to panels (c) and (d), respectively.(b) Same as (a) but in logarithmic color scale, in order to highlight the weaker features.The black dash-dots are a guide to the eye for the line that blueshifts with V 0 and generates the avoided crossing with the repulsive polaron branch RP.AP indicates the attractive polaron.Panels (c) and (d) report the density-density correlator 〈 f 0 |n (0,0) n (x, y) | f 0 〉 at small and large V 0 , respectively.The lattice is a 6 × 6 square with periodic boundary conditions and N ↑ = 4 fermions.

Figure 3 :
Figure 3: Polaron spectrum for increasing values of Coulomb repulsion on a triangular lattice (left) and on a chain (right), respectively.Pink dashes indicate the position of the umklapp peak (not visible on this scale), which is absent in the 1D system.

Figure 4 :
Figure 4: Polaron spectra for a spinful Fermi system with a slightly spin-asymmetric binding energy to the impurity.On the x-axis E B is decreased and panels (a-c) correspond to increasing contact repulsive interactions, V ↑↓ = 0, 2E F , 4E F respectively.The central panel is particularly reminescent of spectral measurements in WSe 2 monolayers.The main lines are labeled in panel (a): RP denotes the repulsive branch, AP ↑ (AP ↓ ) is the attractive polaron associated with the 2-body bound state between the impurity and the spin up (down) fermion, while AP ↑↓ comes from the 3-body bound state.The cyan dotted, dashed and solid lines at small E F /E B indicate, respectively, the AP ↓ , AP ↑ and three-body bound states in vacuum.Computation for a 4 × 4 triangular lattice with N ↑ = N ↓ = 4 fermions.

Figure 5 :
Figure 5: Polaron spectra for the Fermi superfluid, with E pair /E F tuning along the BEC-BCS crossover.Panels (a-c) correspond to increasing spin-symmetry of the fermion-impurity interactions, namely U ↑ −U ↓ U = 2, 1, 0. In the three panels we set E B /E F ≃ 0.3, 0.5, 0.5 respectively and N ↑ = N ↓ = 4 on a 4 × 4 triangular lattice.The cyan lines stand for the energies of the two lowest three-body bound states in vacuum, the ground state E 3 (solid line) and the first excited E * 3 (dash-dotted).The pink arrows highlight the presence of a weaker peak below the main AP line at large E pair /E F .

Figure 6 :
Figure 6: Polaron spectra for spin-density unbalanced Fermi superfluid, with E pair /E F tuning along the BEC-BCS crossover.Panels (a-d) correspond to U ↑ −U ↓ U = 2, 1, 0, −2.In the three panels we fix E B ≃ 0.5E F and unbalanced polarization N ↑ = 4, N ↓ = 3 on a 4 × 4 square lattice.The lines dubbed AP ↑ visible in the deep BEC regime of panels (a-c) come from the binding of the impurity with the extra unpaired spin up fermion, and are highlighted by the cyan dotted boxes.

Figure 7 :
Figure 7: (a,b) Cut of Fermi polaron spectra at E B = E F , for different number of particles, system sizes and shapes.The black dashes correspond to a continuum Chevy calculation.Panel (c) shows convergence in the Krylov space dimension.

Figure 8 :
Figure 8: Umklapp peak phenomenology.(a) Spectral function for N ↑ = 4 in a 6 × 6 triangular lattice and three different fermion-impurity interaction constants.(b) Same for different microscopic lattices, with E B ≃ E F .In all cases, the strongly repulsive regime V 0 /t I ∼ 100 is considered, for which one has a very rigid CDW.

Figure 13 :
Figure 13: Polaron spectra for a Fermi superfluid from a Chevy calculation on top of a BCS ground state.The three panels correspond to U ↑ −U ↓ U