Fermionic matter-wave quantum optics with cold-atom impurity models

Motivated by recent cold-atom realisations of matter-wave waveguide QED, we study simple fermionic impurity models and discuss fermionic analogues of several paradigmatic phenomena in quantum optics, including formation of non-trivial bound states, (matter-wave) emission dynamics, and collective dissipation. For a single impurity, we highlight interesting ground-state features, focusing in particular on real-space signatures of an emergent length scale associated with an impurity screening cloud. We also present novel non-Markovian many-body effects in the quench dynamics of single- and multiple-impurity systems, including fractional decay around the Fermi level and multi-excitation population trapping due to bound states in the continuum.

Crucially, the study of quantum emitters in structured baths has been largely restricted to single-particle physics, with only few attempts to move to the multiexcitation regime [30][31][32].Consequently, very little is known about the many-body physics of such models.In particular, the role of particle statistics (i.e. the difference between fermionic and bosonic particles) has, to our knowledge, not been thoroughly examined.In the field of condensed matter physics, fermionic impurity models are well-established and play a significant role in a wide variety of systems.For instance, individual localised fermionic impurity levels can affect transport properties of mesoscopic systems [33][34][35][36], and their presence in superconductors can lead to the emergence of localised quasi-particle states [37] or Majorana modes [38].
By contrast, such models have not received much attention from the quantum optics community, perhaps due to a perceived lack of applicability.However, recent proposals [3,41,42] and subsequent realisations [18,43] of cold-atom matter-wave analogues to traditional waveguide QED setups provide a quite natural pathway to exploring these models from a quantum optics perspective: by trading in the bosonic atoms used in these setups for fermionic ones, such experiments could unlock a new scenario of "fermionic matter-wave quantum optics".
In this work, we explore this bridge between the fields of quantum optics and condensed matter physics in the context of a system of fermionic impurities coupled to a structured bath.Specifically, we model the bath as a lattice of non-interacting spinless fermions and the impurities as additional fermionic modes coupled by local hopping to different sites of the lattice.We analyze separately the case of a single impurity and the case of multiple impurities.In either case, we consider both onedimensional (1D) and two-dimensional (2D) lattices, and we analyze in detail the thermodynamic limit, in which the lattice becomes infinite.
The single-impurity version of our model, known as the Resonant Level Model (RLM) [44], has been extensively studied in 1D, but remains largely unexplored in 2D beyond the effectively 1D description obtained by assuming that the impurity level couples only to l = 0 bath modes.We explore the ground-state features of this model, incorporating both the effects of a finite bandwidth in the bath and a 2D lattice without full rotational symmetry, and present signatures of a screening cloud in the scaling of impurity-bath correlations as a function of distance r from the impurity.Specifically, in 1D we observe a cross-over from logarithmic scaling to the characteristic r −1 scaling of correlations in a bath without an impurity, consistent with previous findings [45,46].However, we show that this scaling is universal over a much wider range of parameter values than previously believed.In 2D, we show that the long-range scaling of the correla-arXiv:2305.11610v2[quant-ph] 6 Jul 2023 tions is also the same as that found in a bath without an impurity, of the form r −3/2 , contrary to the behaviour predicted in the literature [45,46].
We also study the quench dynamics of an initial state in which the impurity modes are occupied and the bath is in its Fermi sea ground state.We use a master equation to describe the dynamics in the Markovian regime [47].To describe non-Markovian effects, we adapt the singleexcitation resolvent formalism from quantum optics to the multi-excitation sector of our model.We then discuss the circumstances under which the system relaxes to the ground state, and clarify the role of the singleparticle bound states of the system alluded to above in this many-body setting.In particular, for a single impurity, we show that fractional decay occurs not only at the band edges, but also around the Fermi level, and how the emission dynamics differs depending on the value of the Fermi level relative to the impurity on-site potential.In the multiple-impurity case, we show that the dynamics generally follow a multi-exponential decay, which is fractional in certain impurity configurations supporting bound states in the continuum (BICs), and we generalise the theory of BICs in 2D baths presented in Refs.[2,4].
Finally, we discuss the feasibility of observing some of these effects in state-of-the-art cold-atom experiments.In particular, we show that the signatures of the groundstate screening cloud can be detected in finite systems, with sizes that are readily achievable in experiments, and we suggest how to prepare the many-body ground state, using an adiabatic or a dynamical protocol.We also discuss the effects of finite temperature on the observation of non-Markovian dynamical phenomena.
The paper is structured as follows: In Sec.II, we introduce our model and develop the theoretical tools used throughout.We then discuss the physics of a single impurity (Sec.III) and multiple impurities (Sec.IV), before addressing experimental factors in Sec.V.

II. THEORETICAL FRAMEWORK
In this section, we introduce our model and the notation and conventions employed throughout this work.We briefly review the theory of free-fermion models and fermionic Gaussian states and introduce both a perturbative and an exact method for studying the dynamics of our model, based on a Markovian master equation and an extension the resolvent formalism, respectively.

A. Fermionic impurity models
We consider N spinless fermionic impurities, described by creation and annihilation operators {c † n , c n } N n=1 , coupled via local tunneling to a d-dimensional spinless fermionic bath of L d sites.The bath modes are likewise described by operators {b † j , b j } L d j=1 , and we denote the site to which the n-th impurity is coupled by j n .We restrict our attention to the 1D (d = 1) and 2D (d = 2) cases.The Hamiltonian for this system can be written as H = H S + H B + V , with ( = 1) Since the bath is translation-invariant, and it has a single site per unit cell, it can be diagonalised by a Fourier transform.
we can express where we use the shorthand notation r n ≡ r jn , with r j denoting the position of the j-th bath site.The dispersion relation in the bath takes the form and the on-site potential ∆ constitutes a detuning with respect to the lower band edge.
The Hamiltonian H conserves the number of particles.Hence, when studying the dynamics it produces for a given initial state with a well-defined number of particles, we can restrict the Hilbert space to the subspace of states with that same number of particles.In the thermodynamic limit, we instead consider the filling fraction or, equivalently, the Fermi level E F , which is defined in terms of the ground state of H B (which also conserves the number of particles) as the highest energy of the occupied bath eigenmodes.This ground state is where "FS" stands for "Fermi sea", and |vac denotes the vacuum state.A Fermi level E F ≤ 0 corresponds to an empty bath, whereas E F ≥ 2 d+1 J corresponds to a fully occupied bath.The standard Hamiltonian employed in quantum optics to describe quantum emitters coupled to structured (photonic) environments is recovered from Eq. ( 1) by replacing the impurity and bath modes with spin and boson operators, respectively.This spin/boson model and our fermionic model share the same single-particle physics, however they differ in their many-body physics due to the different particle statistics and the intrinsic non-linearity introduced by spin impurities.We therefore want to move beyond the single-excitation regime and in the rest of this section we develop the formalism that makes this possible.

B. Lattice Green's function
Throughout this work, we will repeatedly refer to the lattice Green's function of the bath, which we also call the self-energy function and which is defined in the thermodynamic limit (L → ∞) as [48] Σ For a 1D bath, an analytical expression can be obtained (see Appendix A).However, in 2D, there is no general analytical expression, except for r = (n, n) (n ∈ Z), in which case Σ(z, r) can be expressed in terms of hypergeometric functions [49].The lattice Green's function does satisfy certain recursion relations that can in principle be used to compute it for other values of r [50,51].In practice, however, we find it more convenient to reduce Eq. ( 4) to a 1D integral which can then be computed numerically in a simple manner (see Appendix A).

C. Markovian master equation
Under a Born-Markov approximation [47,52], the bath modes can be traced out (see Appendix B) to obtain a master equation for the impurity density matrix ρ, where we have defined H S = H S + m,n J mn c † m c n and Physically, H S captures the coherent dynamics, including coherent interactions mediated by the bath, while D > and D < describe incoherent processes of fermion emission from the impurities into the bath and fermion absorption from the bath into the impurities, respectively.The coherent couplings J mn and collective absorption/emission rates Γ ≶ mn are related to the self-energy function (4) by Σ(∆ ± , r nm ) = J mn ∓ iΓ mn /2 and where x ± ≡ x±i0 + and r nm ≡ r m −r n .Note how Eq. ( 7) implies that only one of the two dissipators, D < or D > , is present for any given values of ∆ and E F .Qualitatively, the master equation dynamics are very intuitive: if ∆ lies inside the band but above the Fermi level, there will be unoccupied bath modes available at the impurity energy, allowing occupied impurities to emit a fermion at that energy.Conversely, below the Fermi level, occupied bath modes resonant with the impurities will populate vacant impurities.If ∆ lies outside the band, Γ mn vanishes (essentially as a consequence of Fermi's Golden Rule) and the dynamics are purely coherent.Notably, the associated interaction strengths J mn are identical to those of conventional two-level systems coupled to a bosonic bath.This reflects the fact that the interactions are mediated by certain localised single-particle eigenstates [53], which are of course insensitive to the particle statistics.The Markovian approximation underlying Eq. ( 5) amounts to the assumption of a memory-less bath which remains approximately in thermal equilibrium throughout the dynamics (see Appendix B).This assumption is valid only when ρ(ω)f (ω) ≈ ρ(∆)f (∆) for ω varied across an interval of width ∼ g 2 /J around ∆ [47].Here, ρ(ω) = −ImΣ(ω + , 0)/(πg 2 ) is the density of states in the bath and f (ω) is the Fermi-Dirac distribution characterising the equilibrium state of the bath, which we assume to be the zero-temperature Fermi sea state (3), f (ω) = Θ(E F − ω).This implies that Eq. ( 5) will fail to accurately describe the dynamics when ∆ lies close to E F (where f (ω) changes abruptly) or whenever ρ(ω) varies significantly around ∆.In these regimes, we observe non-Markovian effects.To describe them, we develop an exact description of the system in the next two subsections.

D. Free-fermion formalism
The Hamiltonian (1) is quadratic in the fermionic operators, reflecting the fact that the fermions are noninteracting.Below, we note some of the key properties of such free-fermion systems.

Fermionic Fock states
Consider a set of M independent fermionic creation and annihilation operators {ψ † l , ψ l } M l=1 .We assume the set to be complete in the sense that the states {|ψ l ≡ ψ † l |vac } M l=1 form an orthonormal basis of the singleparticle Hilbert space.A basis of the whole many-body Hilbert space is then given by the Fock states where n ∈ {0, 1} M , such that its l-th entry n l represents the occupation of the l-th mode in the state.It is straightforward to show that for an arbitrary quadratic (or one-body) operator [54] where e l is the unit vector with zeros everywhere except for the l-th component.
One natural choice of operator basis are the singleparticle eigenmodes {φ † l , φ l } M l=1 , which diagonalise the Hamiltonian as H = l ε l φ † l φ l .The many-body eigenstates of H are then given by Fock states of the form (8) in the basis formed by these eigenmodes.In particular, the ground state for a given Fermi level (filling) is According to Eq. ( 10), ground-state expectation values of arbitrary one-body operators therefore simply read

Time-evolved expectation values
The basis of single-particle eigenmodes {φ † l , φ l } M l=1 also lends itself to describing time evolution: the expectation value of a generic quadratic operator O at time t is where • 0 denotes an expectation value in the initial state.In studying the dynamics of our system, we will repeatedly consider the long-term average (LTA), It follows from Eq. ( 13) that for a quadratic operator, i.e. it is determined entirely by the correlations between eigenmodes in the initial state.When O t converges to a specific value as t → ∞, this limit coincides with O .

Gaussian state formalism
The theory of quadratic Hamiltonians is inextricably linked with the theory of Gaussian states, since Gaussianity is preserved under evolution with a quadratic Hamiltonian.A remarkable property of Gaussian states is the fact that any even-order correlations can be expressed entirely in terms of 2-point correlations according to Wick's Theorem.Thus, Gaussian states are characterised fully by their covariance matrix [55][56][57][58].This property has been exploited to develop a number of methods for studying fermionic condensed matter models, such as a time-dependent formalism for calculating spectral functions [59] and the Functional Determinant approach to electron transport and dynamics [60][61][62][63][64].
If the total number of particles in the system is welldefined and preserved during evolution, as is the case for the Hamiltonian (1), it is sufficient to consider the correlation matrix C(t) (or equivalently the one-particle reduced density matrix) whose elements are C lm (t) = ψ † l ψ m t .It evolves as [55,56] where U lm (t) = ψ l | e −iH T t |ψ m .If the system is timereversal invariant, as in our case, then H is real and H T = H, so that U (t) is simply the matrix representation of the single-particle time-evolution operator.Crucially, the correlation matrix formalism not only provides us with an economical way of computing real-time dynamics of large systems, but also allows us to extend some of the analytical tools typically associated with the single-excitation regime to the many-body regime, as we discuss in the next section.

E. Resolvent formalism
The resolvent formalism allows the exact calculation of the transition amplitudes between eigenstates of a Hamiltonian H 0 induced by an interaction V [65].In the context of quantum emitters coupled to reservoirs, it has proven powerful for understanding the non-Markovian features of the dynamics in the single-excitation sector.

Basic formalism
In our case, H 0 = H S + H B with single-particle eigenstates |e n ≡ c † n |vac and |k ≡ b † k |vac .For now, we will focus on the case of a single impurity and re-label |e 1 ≡ |e , although the formalism can also be applied to configurations with multiple impurities [2].The transition amplitudes A e (t) ≡ e| e −iHt |e , A k (t) ≡ k| e −iHt |e , and A qk (t) ≡ q| e −iHt |k can be computed as with the propagators defined (with Σ e (z) ≡ Σ(z, 0)) as Formally, these are the matrix elements of the resolvent of the Hamiltonian between single-particle eigenstates of H 0 [53].The origin of non-Markovian dynamical features then lies in the analytic properties of the G α (z) [2,65].

Multi-excitation quench dynamics
We now show how the dynamics of a multi-excitation state can be expressed in terms the amplitudes A α (t), thus effectively reducing the many-body problem to a (well-studied) single-particle problem.Since we are working in the Gaussian regime, we can characterize the state of the system completely by its 2-point correlations c † n c m , c † n b k , and b † k b q .Consider the initial state |ψ 0 = c † |FS , which will be relevant for the following discussion.For this state, e| e −iHt |e where we have first re-written the correlator as c † c t = ψ 0 | e iHt c † c e −iHt |ψ 0 and then applied Eq. ( 10) to this expectation value, using the fact that |ψ 0 is a fermionic Fock state of the form (8), and H conserves the total number of particles.We have also used the fact that A α (−t) = A α (t) for a time-reversal invariant model.The other 2-point correlators can be obtained analogously as Together with Eq. (17), this gives us a way to calculate the exact dynamics (semi-)analytically.More generally, the amplitudes computed via Eq.( 17) provide an expression for U (t) in a particular basis, which can then be used to compute the time-dependence of the correlations for any initial state via Eq.( 16).These formulae also provide a basis for developing various approximations.For example, one can generalize the single-pole or Wigner-Weisskopf approximation [2,65] to the many body regime, replacing Σ e (z) by Σ e (∆ + ) in the different propagators, but still using Eqs.(19) and (20) to compute the time-dependence of the two-point correlations.This generally leads to a different dynamics than that predicted by the master equation, in contrast to what happens in the single-particle case, where both approaches lead to the same dynamics.This is because the conditions for the Wigner-Weisskopf and master equation descriptions to be accurate, i.e. essentially ρ(ω) ≈ ρ(∆) and ρ(ω)f (ω) ≈ ρ(∆)f (∆), respectively, are not equivalent for a non-trivial bath particle distribution f (ω).

Long-term averages
We can also compute LTA expectation values easily from the integral expressions for the transition amplitudes.By Eq. ( 14), only the real poles of the propagators G α (z), which contribute terms to the A α (t) that do not decacy in time, will contribute to the LTA of the correlation matrix elements.From Eqs. (18), we see that G e (z) has two real poles, z = ω ± , which are solutions to while G k (z) shares the real poles of G e (z) and has the additional real pole ω k .Similarly, G qk (z) shares all poles of G k (z) and has the additional pole ω q .Thus, the LTA impurity occupation, for instance, can be computed as where the first and second terms are due to the real poles ω ± of G e (z) and the additional pole ] −1 is the residue associated with simple poles of G e (z).

F. Fermions vs. bosons
Before proceeding, it is worth noting that the theory we have developed is also largely applicable to the bosonic counterpart of our model, where the impurity and bath modes describe non-interacting bosons (in which case the Hamiltonian is again quadratic).
Firstly, Eq. ( 10) is formally identical in the bosonic case.In fact, since the single-particle eigenspectrum of quadratic models is not sensitive to particle statistics, the dynamics generated by our Hamiltonian according to Eq. ( 16) are also identical.From the perspective of Gaussian states, this is due to the fact that the Hamiltonian preserves occupation number and thus only generates the subset of Gaussian transformations common to both the fermionic and bosonic Gaussian state families [57,58].
The crucial difference between fermions and bosons therefore lies in the possible initial states.It is the choice of a Fermi sea state in the bath-which is not exclusively fermionic per se, but which of course does not describe the ground state of a bosonic bath-which leads to the explicit form of Eqs. ( 5), (12), and ( 19)- (20).Indeed, general bosonic Fock states are not Gaussian, and therefore, while Eq. ( 16) still captures the evolution of the bosonic two-point correlations for the intial states we consider in this work, these do not in general characterise the full state.The transition amplitudes A α (t) themselves are once again independent of particle statistics, and so the notion of describing the many-body physics of the model in terms of its single-particle properties translates straightforwardly to the bosonic regime.Indeed, there also exists a bosonic version of Functional Determinant formalism mentioned above [66,67].

III. SINGLE IMPURITY
Having developed an extensive theoretical toolbox to study the Hamiltonian (1), we now apply it first to the case of a single impurity.Without loss of generality, we assume that the impurity is coupled to the bath at r 1 ≡ 0, so that the Hamiltonian becomes From a condensed-matter point of view, this Hamiltonian constitutes a minimal model for certain mesoscopic systems (semiconductor structures, e.g.quantum dots) [44].Moreover, it represents a limiting case of more complicated impurity models which include spin degrees of freedom and / or interactions between fermions: the Interacting Resonant Level Model [68], the Single-impurity Anderson Model [69], and the Kondo Model [70,71].As noted above, from the perspective of quantum optics, it shares the same single-particle physics as the spin/boson models that have been studied recently (see e.g.Ref. [2]).However, the fermionic character of the excitations leads to decidedly different many-body physics, as we demonstrate in this section.
The simplicity of our model allows us to exploit powerful techniques from quantum optics to describe the system exactly, in an unprecedented level of detail.Specifically, we firstly analyze the ground state properties of the Hamiltonian (23) in detail.Contrary to the standard approach followed when studying the RLM, we account for the full structure of the (finite) band of bath modes, rather than linearising around the Fermi level.Thereby, we are able to derive formally exact expressions for the ground-state correlations, allowing us to generalise existing results on the impurity-bath correlations in 1D [45] and to thoroughly examine these correlations in 2D.
Secondly, we consider the dynamics of the model, for an initial state in which the impurity mode is occupied and the bath is in the Fermi sea state (3).The formalism developed in the previous section allows us to study these dynamics in the non-Markovian regime exactly and without any approximations.
A. Ground-state features

Single-particle eigenstates
The single-particle spectrum of H comprises two distinct types of eigenstates: a continuum of scattering eigenstates |φ k , parametrised by quasi-momenta k ∈ BZ, with energies ω k [53,72], where |j ≡ b † j |vac , and two bound states (BS) |φ ± , with energies ω ± given by the real solutions to the pole equation ( 21), satisfying ω + > 2 d+1 J and ω − < 0 [30], The presence of two bound states is due to the divergent Lamb shift (with opposite sign) at the two band edges and is therefore a consequence of the finite bandwidth of the bath.In fact, for a finite but sufficiently large bath, coupling to the impurity mode results in two localized modes and L−1 scattering eigenstates in the band range.The many-body ground state is a product state of the form shown in Eq. ( 11), in which these single-particle eigenstates of energies up to E F are occupied.In the following, we focus on the case where the Fermi level lies within the band, 0 ≤ E F ≤ 2 d+1 J, such that only the lower bound state (LBS) and the scattering eigenstates of energies ω k < E F are occupied in the ground state.

Impurity occupation
Using the single-particle eigenstate formulae above, the ground-state impurity occupation in the thermodynamic limit takes the form If the coupling is small enough and the detuning is deep enough in the band, we can neglect the first term, which is the contribution from the LBS [2], since If we additionally extend the lower integration limit to −∞ under an infinite-band approximation [73], Thus, the impurity is either occupied or empty for ∆ E F or ∆ E F , respectively.The transition point occurs in a range of energies of width ∼ Γ F around E F − δω F .Clearly, this reasoning breaks down in the 2D bath at the Van Hove singularity (E F ≈ 4J), where Γ F diverges [48].In this case, a simplified integral expression for c † c GS can be obtained by substituting an asymptotic expression for Σ e (ω + ) as ω → 4J.However, this does not lead to a qualitative change in the dependence of c † c GS on ∆.

Impurity screening cloud
A more salient ground-state property of the singleimpurity model is an emergent length scale associated with an impurity screening cloud, which constitutes a hallmark feature of generic impurity models [40].It is reflected, for example, in the spatial dependence of the envelope of c † b j GS [45,74,75], which we denote here by F (r j ).Using the expressions for the single-particle eigenstates and Eq. ( 12), in the thermodynamic limit, The first term, coming from the LBS, decays exponentially or faster as the distance , and, as noted above, its prefactor R e (ω − ) is very small if ∆ lies sufficiently deep inside the band.The second term gives the contribution of the scattering eigenstates, and it decays algebraically, as ∼ |r j | −α , so it dominates at large distances.Previous works seem to suggest a constant exponent α = 1 for 1 ≤ d ≤ 3 [45].However, as we show in Appendix C, this algebraic decay is the same as the one of the bathbath correlations b † 0 b j FS in an impurity-free bath with the same Fermi level, and is ultimately determined by the dimensionality, and the smoothness of the Fermi surface.For a smooth Fermi surface, α = (d + 1)/2 [77].Note that this algebraic dependence is expected for critical (gapless) phases.
In 1D, if, in addition, the Fermi level is deep inside the band, we can make a wide-band approximation (see Appendix C) and neglect the LBS contribution to obtain , where E 1 is the exponential integral function, v F = 2g 2 /Γ F is the group velocity at the Fermi level, and We plot the envelope F (r j ) of the approximate correlations in Fig. 1a, together with a numerical evaluation of the exact formula (29).Remarkably, good agreement between the two is still observed when ∆ lies close to the band edge.The short-and long-range scaling of the correlations, also indicated in Fig. 1a, are dictated by This generalises what was already observed in Ref. [45] for the special case ∆ = E F : there is a transition from logarithmic to algebraic decay of the ground-state impurity-bath correlator at |r j | ≈ ξ 1D .In this sense, ξ 1D can be interpreted as the size of the screening cloud.
Owing to the lack of an analytical expression for Σ(ω k , r j ), we cannot obtain a "universal" formula, as Eq. ( 30), for the correlations in the 2D case.Nonetheless, we compute them using an efficient numerical method (detailed in Appendix C) that allows us to reach very large system sizes, and thus to observe the onset of their algebraic decay and determine the screening length for different directions, see Fig. 1b.In this way, we can map out the shape of the screening cloud in 2D.One might expect certain degree of anisotropy, especially if the Fermi surface is highly anisotropic, for example, when the Fermi level is tuned close to the Van Hove singularity.However, for the cases considered, the screening length seems to have a rather small angular dependence, irrespective of the value of the Fermi level.If we instead look at the actual values of the correlations along different directions, we do find larger differences the more anisotropic the Fermi surface is.
A similar angular dependence can be observed in the ground-state bath-bath correlator b † 0 b j GS , which in the thermodynamic limit reads, neglecting an exponentially decaying contribution from the LBS, Here, we have defined • FS = FS| • |FS , i.e. the first term captures the ground-state correlations in a bath without an impurity, while the second term represents a local deformation of the correlations due to the impurity.We plot this deformation Fig. 1c.

B. Quench dynamics
We now consider the dynamics of the initial state |ψ 0 = c † |FS .This amounts to initialising the bath in its ground state at the desired filling, populating the impurity, and turning on a coupling g > 0 instantaneously at time t = 0. Below we describe the transient dynamics of both the impurity and the bath, as well as the steady state after the quench.

Impurity occupation
The dynamics of an initially empty bath and an excited impurity mode have been studied extensively in the spin/bosonic setting [2] and the results can be directly extrapolated to our model.In this regime, for ∆ sufficiently deep in the band of the bath, the impurity occupation displays approximately exponential decay at rate Γ 0 = −2 Im Σ e (∆ + ), however when ∆ lies close to the band edges, we observe fractional decay due to trapping of the excitation in the bound states [2,14,15].Addi-tionally, in the 2D bath, when ∆ = 4J (i.e.tuned to the Van Hove singularity) we observe algebraic decay [2,4].
For a non-zero filling, E F > 0, the dynamics depend primarily on the value of the detuning relative to the Fermi level.The master equation ( 5) predicts c † c t = e −Γ0t , when ∆ > E F , and c † c t = 1, otherwise.However, the exact dynamics reveal some non-Markovian features.In addition to the effect of the bound states, we can observe fractional decay when ∆ ≈ E F , see Fig. 2a.In 2D, when ∆ = 4J we find that already at a very low value of E F , the non-zero steady-state population induced by the Fermi level overwhelms the algebraic decay (see Fig. 2b).
If the coupling g is sufficiently small, the impurity occupation converges to a specific value as time increases, which is given by its LTA c † c , Eq. ( 22).This quantity is shown as a function of the detuning in Fig. 2c.Notably, the profile of the LTA impurity occupancy is "smeared out" with respect to the master equation prediction around the band edges and the Fermi level, reflecting the non-Markovian possibility of decay or absence of decay even where Fermi's Golden Rule would prohibit it.Essentially, the observation of fractional decay around the Fermi level can therefore be viewed as a non-Markovian version of Pauli-inhibited spontaneous decay [78].In fact, comparing the expression for c † c , Eq. ( 22), and the impurity occupation in the ground state c † c GS , Eq. ( 26), one realizes that both are the same, except for contributions attributed to the bound states.As we show in the next section, this kind of relationship between the long-term average (steady state) expectation value and the ground-state expectation value actually holds for arbitrary one-body operators.
Before proceeding, it is worth noting that, for detunings deep in the band (where R e (ω ± ) is negligible) and in the weak coupling regime, the fractional decay around the Fermi level is described by Eq. (28).Notably, this formula has been previously obtained in the 1D case as the steady-state population for the Resonant Level Model in the Keldysh formalism [73,79].Of course, the result we have obtained in Eq. ( 22) is more general, since it does not rely on a weak-coupling approximation and incorporates the effects of the band edges.

Relaxation to the ground state
Having established the relation between the steadystate and ground state impurity occupation above, we now discuss this relation for a general on-body operator O.According to Eq. ( 15), the LTA expectation value of any such operator after the quench is given by However, as we show in Appendix D, in the thermodynamic limit we find that for an initial state of the form |ψ 0 = c † ne |FS , with n e ∈ {0, 1}, correlations between degenerate scattering eigenmodes are the same in the Fermi sea and in the ground state of the impurity model, if ω k = ω q .This result is physically intuitive: for a free bath mode, the impurities constitute a potential with a finite range, whose only effect is to induce a phase shift as the particle is scattered [80].In terms of the LTA impurity occupation, Eq. ( 36) implies that, in the thermodynamic limit, For example, using Eq. ( 10), we have, for and noting that φ ± | c † c |φ ± = R e (ω ± ) we can check that Eq. ( 22) is consistent with Eq. (37).In essence, Eq. ( 37) formalizes the idea that the initial state |ψ 0 = c † |FS indeed relaxes to the ground state, except for a contribution from the BSs, which constitutes the effect of the finite bandwidth.

Bath excitation dynamics
Another interesting feature of the quench dynamics is the propagation of excitations through the bath following the decay of the impurity.This is captured by the timedependent expecation value b † j b j t .From Eqs. (20), We focus on the case where ∆ lies in the band, in which case we observe an emission of matter-waves into the bath.We also restrict our attention to the case of a 1D bath, which allows us to obtain some analytical results.In particular, we can consider the limit t → ∞ and x/t → κ.Under a stationary phase approximation [81], where , and k κ = arcsin(κ/(2J)).From this expression, we can distinguish three distinct regimes (see Fig. 3): for E F = 0, a wavefront forms and propagates at velocity ∼ v ∆ (i.e. the group velocity at the detuning) due to the fact that the momenta of the emitted matterwaves are strongly peaked around the detuning momentum k ∆ = arccos(1 − ∆/(2J)) (with ω k∆ = ∆) [82].For 0 < E F < ∆, the dynamics display qualitatively similar features to the case E F = 0, however now there is a significant occupation probability for sites around the impurity site, given by the ground state occupations.For E F > ∆, the emitter matter-waves are now predominantly confined to a new light cone, defined not by v ∆ but rather by the Fermi velocity v F , because modes with |k| < k F are no longer available for propagation.

IV. MULTIPLE IMPURITIES
We now turn to the problem of N ≥ 2 impurities.We will focus our attention on the collective decay of a "fully inverted" state, i.e. the dynamics of the total impurity occupation N e = n c † n c n starting from the initial state |ψ 0 = n c † n |FS .In the case of two-level impurities decaying collectively into an empty bath, this setting gives rise to a characteristic delayed superradiant burst of emission, due to a build-up of coherence between the impurities [29].While this setting of collective dynamics of multiple impurities coupled to a common bath is natural from the perspective of (matter-wave) quantum optics, it has, to our knowledge, not been studied in the context of condensed matter physics.In this section, we therefore discuss both the Markovian description of the dynamics as well as non-Markovian effects.We find that much of the conceptual discussion of the single-impurity dynamics can be intuitively extended to this more general case.In particular, we can again relate the steadystate of the system to its ground state and single-particle bound states in the spectrum of our model.

A. General description
In the Markovian regime, the collective decay is described by the master equation ( 5).The impurity correlations C mn (t) ≡ c † m c n t thus evolve in the same way as for any quadratic model, as shown in Eq. ( 16), but with a non-Hermitian effective Hamiltonian, which in the case ∆ > E F is given by H eff = m,n Σ mn (∆ + ) c † m c n , where we have defined the self-energy matrix Σ(z) as the matrix with elements Σ mn (z) ≡ Σ(z, r mn ).Thus, for a fully inverted state, C(0) = I N ×N , the total impurity occupation N e t = tr C(t) evolves as where • HS denotes the Hilbert-Schmidt norm.With reference to the discussion in Sec.II F, it is worth noting that Eq. ( 41) also holds for the bosonic counterpart of our model, in the case where the bath is assumed to remain approximately in its ground state, i.e. the vacuum state [56].
The exact dynamics feature a number of non-Markovian effects not captured by the master equation.For instance, the master equation neglects retardation effects in the bath.In reality, we have already shown that emitted fermions propagate though the bath at a finite velocity dictated by ∆ or E F (see Fig. 3).Accordingly, any collective dissipative effects will be delayed by the time taken for fermions emitted at one impurity site to propagate to the next.Moreover, even after this delay, the exact dynamics will only agree with the master equation for sufficiently small impurity separations.
The collective dynamics can also show fractional decay similar to the single impurity case.This can again be captured in terms of the LTA occupation, N e .As we show in Appendix D, we can generalize Eq. ( 37) for initial product states of the form |ψ 0 = |ψ e ⊗ |FS , where |ψ e is an arbitrary impurity state with a well-defined number of particles, as Here, the sums run through all single-particle bound eigenstates of H.In the multiple-impurity case, this includes not only BS with energies outside the band, but may also include localised dressed states emerging at the impurity energy ∆, where ∆ lies in the band.In recent work [53,83], a general framework has been developed for such states, which are known as bound states in the continuum (BIC) but have also been coined vacancy-like dressed states (VDS), since it can be shown that such states have zero amplitude on the sites to which emitters are coupled.This leads to a simple prescription for computing the wavefunction of these states, which we recall in Appendix E. Whether a BIC exists in the first place can be determined from the condition [53] det which is a generalisation of the pole equation ( 21) to the multiple-impurity case, for eigenstates at energy ∆.

B. 1D bath
As shown in Appendix A, in 1D the coherent and dissipative couplings in Eq. ( 5) read where Γ 0 is the single-impurity decay rate.For simplicity, we can consider configurations with the emitters equidistantly placed at a separation d, such that r mn = (n−m)d.In Fig. 4a, we show the collective decay of the fully inverted state into an initially empty bath for such a configuration.While the master equation predicts the same multi-exponential decay dynamics for different values of d yielding the same couplings J mn and Γ mn , the exact dynamics display the retardation effect discussed above, and the time taken for the onset of cooperative decay intuitively increases with d.Note also that, as expected, the agreement between the exact and Markovian dynamics at later times is better for smaller values of d.For E F > 0, we also observe the anticipated fractional decay around the Fermi level (see Fig. 4b): N e t relaxes to the non-zero value (42), computed numerically.A particularly interesting 1D emitter configuration is the one where the emitter separation obeys cos(k ∆ d) = ±1.In this case, Eqs.(44) imply J mn = 0 and Γ mn = (±1) |m−n| Γ 0 , whereby the master equation ( 5) takes the form of the Dicke master equation [84], therefore we refer to this situation as the Dicke regime.In the case For reference, we also plot the master equation dynamics (black solid) and the independent decay of N = 5 impurities (grey solid).b) Fractional decay of N = 5 impurities at different Fermi levels, for ∆ = J, g = 0.4J, and d = 5, and independent decay of N = 5 impurities at same values of EF (grey solid).The LTA occupation, given by Eq (42), is also shown (grey dashed).
of atoms (two-level systems), the permutation-invariance of the Dicke equation together with the spin statistics leads to the characteristic scaling of the peak emission rate of the superradiant burst as ∼ N 2 [84].However, in our case, fermionic statistics have essentially the opposite effect: In this regime, the self-energy matrix Σ(∆ + ) can be diagonalised by a unitary transformation, whereby N e t = n e 2 Im λnt , with {λ n } N n=1 denoting the eigenvalues of Σ(∆ + ), all of which have non-positive imaginary part.In the Dicke regime, there is actually only a single non-zero eigenvalue, such that implying lim t→∞ N e t = N − 1.In other words, the dynamics only involve the emission of a single particle because the (in this case unique) bright mode can be ini- tially occupied by at most one fermion.Again, it is worth noting that Eq. ( 45), like the more general Eq. ( 41), still holds in the case where the c n are bosonic operators [42,85].However, in this case the initial state is no longer Gaussian and therefore Eq. ( 45) does not capture the full evolution of the bosonic state, which in principle includes also non-trivial higher-order correlations.The Markovian prediction N e = N − 1 from Eq. ( 45) can be re-framed in terms of population trapping in BICs.In 1D, Eq. ( 43) takes the form which implies that a BIC exists whenever two neighbouring emitters are separated by a distance d ∈ (π/k ∆ )Z.In the vacancy picture, such a pair of emitters "cuts out" a finite chain of d−1 sites, which then supports a eigenstate at energy ∆ localised between the two impurity sites (see Appendix E).This implies in particular that N − 1 such BICs emerge in the Dicke regime.For an empty bath (E F = 0), neglecting the effects of bound states outside the continuum, Eq. ( 42) then predicts a non-zero LTA occupation due to trapping in these BICs given by (see Appendix E) where we have defined A formally exact expression for the case E F > 0 can also be obtained, but is too lengthy to reproduce here.It is worth noting, however, that for N = 2, this more general version of Eq. ( 47) agrees with the value of N e that can be obtained using an extension of the resolvent formalism for the two-impurity case (see Appendix E).We show the agreement of Eq. ( 47) with the dynamics in Fig. 5a.The deviation of R from unity quantifies the retardation in the bath by the ratio between the timescale of the decay (∼ 1/Γ 0 ) and the time taken for excitations to propagate between neighbouring emitters (∼ d/2v ∆ ) [2].For larger d, the population trapped in the BICs therefore decreases (see Fig. 5).Conversely, if we neglect the retardation effects and assume R ≈ 1, Eq. ( 47) reproduces the Markovian result N e = N − 1.

C. 2D bath
In the 2D case, it is possible to find impurity configurations which "cut out" a full sub-lattice surrounded by vacancy-like sites in the same manner as in the 1D case.However, in 2D we also find a second, qualitatively distinct class of BIC, the first example of which was reported in Refs.[2,4].The first distinguishing feature of these states is that they emerge only when the impurity energy is tuned to the Van Hove singularity, i.e. ∆ = 4J, where some matrix elements of Σ(∆ + ) diverge.The second one is the fact that the corresponding impurity arrangement does not enclose a full section of the lattice.Instead, impurities must be placed at the vertices of a rectangle, whose sides are rotated 45 • with respect to the lattice axes.This kind of BIC has nonzero amplitude only in the interior of the rectangle (see Appendix E).We demonstrate the population trapping in two such configurations in Fig. 5.

V. EXPERIMENTAL CONSIDERATIONS A. Cold-atom impurity models
Our proposed matter-wave waveguide QED setup has been discussed in several previous works [3,18,[41][42][43]86, 87] therefore we will not review it in detail here.We simply recall that in these works, it was shown that by driving transitions between two atomic states which are strongly localised in isolated deep traps and delocalised across a more shallow lattice, respectively, a Hamiltonian of the form (1) can be realised.
Alkaline-earth-metal atoms (AEMAs) like ytterbium and strontium have emerged as arguably the most suitable for state-dependent optical lattice experiments of the kind that we envision [88].In fact, the setup required to model our Hamiltonian, with one state localised and one delocalised, has been identified as a formidable platform for simulating complex condensed matter models [89][90][91] and has already been succesfully implemented in a number of experiments, both using bosonic [92] and fermionic [93,94] atoms.
Contrary to our model, it is important to note that such AEMA setups typically feature a (non-perturbative) on-site interaction between atoms in different internal states [89,95], which breaks the assumption of noninteracting particles which underlies all of our results.However, these interactions can be tuned [96,97], for instance through Feschbach resonances [98,99] or by slightly displacing the traps for the two internal states from each other [99].Accordingly, we are confident that our proposed free-fermion matter-wave waveguide QED system is realisable in state-of-the-art experiments.We now want to establish in detail how we believe the main results that we have reported above can be observed in experiment.

B. Quench dynamics
Population trapping in bound states is arguably the simplest phenomenon to observe in cold-atom experiments, as it only requires initialising a few atoms in the deeply trapped internal state at selected positions in the optical lattice.Indeed, at the single-excitation level, experimental confirmation of the effect of a bound state on the single-impurity matter-wave dynamics has already been reported [18].It should therefore be possible to observe the population trapping of multiple excitations in the bound states (both outside and in the continuum) of multiple-impurity systems, as shown in Fig. 5.It should also be possible to observe the collective decay effects (multi-exponential decay) and non-Markovian dynamics (retardation) shown in Fig. 4 as the distance between the impurities is increased.
To explore fractional decay around the Fermi level for one or multiple impurities, we need to prepare the more complex initial state |ψ 0 = n c † n |FS .Of course, the true zero-temperature Fermi sea state cannot be realised in experiment.Instead, the state |ψ 0 would be characterised by the Fermi-Dirac distribution, i.e. b † k b q 0 = δ kq f (ω k ) with f (ω) = 1 + e (ω−E F )/k B T −1 at a temperature T .Such an initial state would no longer be a product state and therefore most of the discussion regarding dynamics in Sec.II would no longer apply.However, we can estimate the effect of this different initial state on our results by numerical evaluation of Eq. ( 16).For instance, Fig. 6a shows the single-impurity LTA occupation, complementary to Fig. 2c, for such a finite-temperature ini-tial state.Clearly, for sufficiently low temperatures, the LTA occupation still qualitatively resembles the zerotemperature profile.In recent experiments, temperatures T /T F ∼ 0.2 (with T F = E F /k B ) have already been reported [100,101] and the phenomenon of fractional decay around the Fermi level should therefore certainly be observable in experiment.
C. Ground-state correlations

Measurement schemes
Another relevant question is the feasibility of measuring single-impurity ground-state correlations like the ones shown in Fig. 1.In particular, we have discussed how the real-space structure of the impurity screening cloud can be inferred from the correlator c † b j GS .The experimental challenge of measuring this correlator directly can be circumvented by exploiting the Gaussianity of our model: Wick's Theorem implies and both terms on the right hand side of this expression can be obtained through site-resolved density measurements of the kind utilised in well-established coldatom microscope experiments [102][103][104][105].While we find that the phase of the correlations c † b j cannot be reconstructed even from higher-order correlations, the magnitude | c † b j | is sufficient to characterise the screening.
On the other hand, the full ground-state correlations, including also their complex phase, could be measured using an alternative approach similar to the photoemission spectroscopy scheme proposed in Ref. [106]: by introducing local modulations of the bath lattice at the impurity site and another site j, a small number of fermions (less than one) can be resonantly excited into an auxiliary, detection lattice, which is initially empty, and whose chemical potential is offset from the original lattice by an energy much larger than the interlattice tunelling.The collective momentum-space modes of this auxiliary lattice right after the modulation (at t = 0) can then be expressed in terms of the original bath operators as bk = b 0 + e ik•rj b j .After a sufficiently long period of ballistic expansion in the detection lattice, the occupation of these modes can be measured using a standard time-of-flight protocol [107] allowing us to infer both | b † 0 b j GS | and ϕ j from the obtained data for (many) different values of k.While the above discussion focuses on bath-bath correlations, it stands to reason that a similar method could also be used to measure the phase of c † b j GS , by coupling only the impurity site (state) and not the bath site to which it is coupled to the auxiliary detection lattice.

Ground-state preparation
The more significant experimental challenge is the preparation of the many-body ground state itself.A wellestablished approach in this context is to start from a different, simpler many-body state and adiabatically tune the parameters of the system to transform this state into the desired ground state [97,[108][109][110].In our case, this would involve preparing the Fermi sea state and adiabatically turning on the coupling g.Crucially, the preparation time depends on the size of the energy gap between the ground state and the first excited state during the protocol [108], which for our model is given by ∼ πv F /L d .For large systems, it therefore becomes unfeasible to adiabatically prepare the ground state.This is exacerbated by the fact that the overlap between any two ground states of impurity models with different coupling strengths g decays with the system size, a phenomenon known as the orthogonality catastrophe [111,112].However, the cross-over from logarithmic to algebraic scaling of the correlations can already be observed in comparatively small finite systems, which can be adiabatically prepared in reasonable time compared to experimental timescales (see Fig. 6b).
An alternative way to observe the impurity screening cloud is given by Eqs. ( 37) and ( 42): In a finite system with ∆ sufficiently deeply in the band, the initial state |ψ 0 = c † |FS relaxes to a state sufficiently resembling the ground state with the same number of excitations which displays the desired impurity screening effect.

VI. CONCLUSION
In this work, we have studied the physics of a spinless fermionic system consisting of one or more impurities coupled to a single-band bath in 1D and 2D.By extending some of the tools originally developed in quantum optics to address the properties of quantum emitters coupled to structured baths, we have been able to incorporate many effects that have typically been neglected in condensed-matter studies of such systems -for example, the effect of band edges and impurity-bath bound states.Conversely, from the perspective of quantum optics, our work provides an insight into the effect of manybody physics and fermionic particle statistics on typical quantum-optical phenomena, such as particle emission, fractional decay, effective bath-mediated interactions, and collective effects such as sub-and superradiance.
Specifically, for a single impurity, we presented exact expressions for the two-point correlation functions which fully capture the ground state.This allowed us to characterise in detail the real-space signatures of an impurity screening cloud in our model.In 1D, we showed that the impurity-bath correlations obey a universal scaling law over a much larger parameter range than previously expected.More generally, in arbitrary dimension d, we demonstrated that the impurity-bath correlations have the same long-range algebraic scaling as the bath-bath correlations in a bath without an impurity, as ∼ |r| −(d+1)/2 for |r| → ∞.
Furthermore, we studied the quench dynamics of an initial state with the impurity populated and a Fermi sea state in the bath.We unveiled qualitatively different emission of matter waves from the impurity into the bath depending on the value of the impurity energy relative to the Fermi level, and we derived a formula capturing the steady state after emission, which clarifies the effect of single-particle bound states on the dynamics and characterises when the system relaxes to its ground state.For multiple impurities, we showed that the quench dynamics of an analogous initial state, in which all impurities are occupied and the bath is in a Fermi sea state, typically display multi-exponential decay of excitations into the bath.We also characterized the steady state of the system and once again highlighted the role of single-particle bound states, which, for certain impurity configurations, can also emerge at energies in the continuum.In 1D, we related the fractional decay in the Dicke regime (i.e. the regime with only a single collective "bright" mode) to the existence of such BICs, which trap the excitations in the impurity modes.In 2D, we also discussed population trapping in BICs, generalizing examples of BICs previously reported in the literature.
Importantly, some of these effects, such as the signatures of single-impurity screening cloud, fractional decay due to Fermi statistics, and population trapping in BICs, could be observed with state-of-the-art cold atom experiments, which allow for the exploration of regimes that are typically inaccessible in more conventional solid-state setups.
Moreover, our proposed scenario of fermionic matterwave quantum optics opens the door to a number of interesting avenues of further research.While we have focused so far on spinless non-interacting models, we expect that including spin and interactions will reveal even more intriguing phenomena intrinsically linked to the fermionic character of the excitations.In particular, the scenario of collective dissipative dynamics of interacting fermionic impurities constitutes a potentially interesting intersection between many-body physics and quantum optics.In the superlative, it could prove fruitful to study the quantum optics of Kondo-type impurities coupled to structured reservoirs on the one hand, and emitters in strongly-correlated baths on the other hand.However, moving beyond our simple models could certainly constitute a considerable computational challenge requiring a more sophisticated theoretical toolbox than we have presented here.

VII. ACKNOWLEDGEMENTS
We are grateful to S. Blatt for helpful discussions and advice.B.W., M.B., and J.I.C. acknowledge funding from the projects FermiQP and EQUAHUMO of the Bildungsministerium für Bildung und Forschung (BMBF), as well as from the Munich Center for Quantum Science and Technology (MCQST), funded by the Deutsche Forschungsgemeinschaft (DFG) under Germany's Excellence Strategy (EXC2111 -390814868).ED acknowledges support from the ARO Grant No. W911NF-16-1-0361 and the SNSF project 200021 212899.

VIII. CODE AVAILABILITY
The computer codes used to obtain the results presented in this work are publicly available at https://doi.org/10.5281/zenodo.8096145and β out = e −iφ .Now, Eq. (A3b) implies φ = arccos(1 − ∆/(2J)) ≡ k ∆ .Substituting back into Eq.(A4), As can be seen from the first equality, the single-emitter self energy Σ e (∆ + ) ≡ Σ(∆ + , 0) is purely imaginary for values of ∆ within the band range, so in the second equality we express the general self-energy function in terms of the single-impurity decay rate Γ 0 .This formula allows us to evaluate the coherent couplings and collective decay rates appearing in the master equation ( 5) in a straightforward manner, leading to Eqs. (44).

Relation to 1D integral in 2D
In the 2D case, where r = (r x , r y ), the self-energy integral can be simplified by integrating over one momentumspace component in the manner outlined above to obtain where β in (k, z) and β out (k, z) denote the roots of the polynomial p(β; k, z) = β 2 + (z/J − 4 + 2 cos k)β + 1, with |β in | < 1 and |β out | > 1.While the remaining integration cannot in general be performed analytically, it does make the evaluation of Σ(z, r) less numerically demanding than performing the full 2D integral.
Appendix B: Fermionic master equation Here, we derive the fermionic master equation, Eq. ( 5), by formally tracing out the bath modes.The starting point is the time-local, Born-Markov master equation [52] expressed here in the interaction picture, O I (t) = e i(H S +H B )t Oe −i(H S +H B )t , and where we take ρ B = |FS FS|.This equation can be derived formally using projection operators [47].Qualitatively, the underlying assumption is that the bath correlation time is short compared to the characteristic timescale of the impurity-bath interaction g, implying that the bath remains approximately in thermal equilibrium despite its coupling to the impurity [47,52].Expanding the nested commutators we can express where we have defined With the explicit form of H S , H B , and V given in Eqs. ( 1) and ( 2), it is easy to see that Substituting this into T 1 (t, t ) and T 2 (t, t ), using the fact that for the particular bath state considered tr The integrals in time can be performed substituting ∆ → ∆ ± i , and then taking the limit → 0, and similarly, (B9) In the last step, we have used the Sokhotski-Plemelj theorem, expressing the end result as the sum of two contributions: one coming from the principal value of the integral, J Last, it is a matter of algebra to show that ρI = −i m, n Going back to the Schrödinger picture, noting also that J < mn + J > mn = J mn , this then yields the result presented in the main text.In the case d = 1, Eq. ( 29) can be simplified using the expression for the self-energy function introduced in Appendix A. It then reads, neglecting the bound state contribution in anticipation of the wide-band approximation, Under a wide-band approximation, we linearise the dispersion relation around the Fermi level as , where v F denotes the group velocity at the Fermi level.We also assume Σ e (ω where For large k F , the integrand is highly oscillatory and hence we can extend the lower limit of integration to −∞.This allows us to evaluate the integral analytically as where Ci(z) and Si(z) are the cosine and sine integrals, respectively [113].In addition, we find that f (z 0 |r j |) is approximately independent of the argument of z 0 , so that we can replace z 0 |r j | by |r j |/ξ 1D with ξ 1D = |z 0 | −1 , as defined in Eq. (31).Noting that Re[e ik Neglecting the LBS contribution, Eq. ( 29) can be cast in the form of a one-dimensional Fourier transform.To do so, one has to perform a change of variables, expressing the integral in terms of the momentum parallel and perpendicular to r j .More concretely, in 2D, for where T is a conformal linear transformation, Thus, the first term in Eq. ( 29) becomes where we have defined the one-dimensional integral Here, r ≡ (m 2 + n 2 )r and Ω k ≡ ω T −1 (k) .Note that the original integration region, BZ, is the square with vertices 2π × {(0, 0), (1, 0), (0, 1), ( 1 where p(y; z, x) is the 2m-order polynomial p(y; z, x) = x n y 2m + x −m y m+n + (z/J − 4)y m + x m y m−n + x −n with roots y = y j .Note that if y j is a simple root, the residue is simply given by Res (1/p(y; z, x), y j ) = 1/p (y j ; z, x).Thus, the second term in Eq. ( 29) can be written as with a new auxiliary one-dimensional integral defined as Finally, we arrive at the compact expression which follows from the fact that the correlations are real, since the total Hamiltonian is real, and For the special case of (m, n) = (1, 0), we can find an analytical expression for σ(z, k), which reads In practice, for general r = (m, n) we compute Re f n (k) numerically, and then evaluate Eq. (C14) using a Fast Fourier Transform algorithm.Note that the same approach can be used to compute the bath-bath correlations (34) in an efficient way, and this is how we obtain Fig. 1c.

b. Long-range scaling
The long-range behavior of the Fourier transform is determined by the singularities of the function to be transformed.The singularities of Re f 1 (k) arise because the integrand in its definition, Eq. (C7), is discontinuous at the Fermi surface, Ω k = E F .Similarly, the integrand in the definition of Re f 2 (k) (after performing the integral in the q-variables) has a logarithmic divergence at the Fermi surface.Thus, if the Fermi surface is smooth, after integrating along the perpendicular momentum k , which results in a long-range scaling of the correlations as ∼ r −3/2 [81].In fact, the singularities are of the same type as the ones one would obtain following a similar procedure for the bath-bath correlations in an impurity-free bath, which can be computed as [77] b A similar analysis can be performed for higher dimensions, arriving at the same conclusion.
Appendix D: Thermodynamic LTA A general eigenstate of the Hamiltonian (1) with energy E can be written as where α = 1, . . ., D(E) labels the degenerate eigenstates at energy E and D(E) denotes the associated degeneracy.The associated eigenvalue equation where we have defined a α = (a α1 , . . ., a αN ) T and g k = gL d/2 (e −ik•r1 , . . ., e −ik•r N ) T .By the second equation, and substituting this solution into the first equation, where we have introduced Note that this definition of Σ(z) is consistent with the definition of the self-energy matrix in Sec.IV of the main text.Thus, the eigenenergies E are solutions to where G(z) ≡ [z − ∆ − Σ(z)] −1 , and the amplitudes on the impurities a α , which fully determine the associated eigenstate through Eq. (D3), can be chosen as unit vectors in the kernel of G −1 (E).It is straightforward to show that the norm of the state is given by where Σ (z) ≡ ∂ z Σ(z).Similarly, the degenerate eigenstates obey the orthogonality condition for α = β.Note that in the thermodynamic limit, the sum in Eq. (D5) diverges for an eigenenergy E within the band of the bath, since the bath energies ω k form a continuum and so for some k in the sum, ω k → E, causing the denominator of the corresponding term to vanish.The norm of the eigenstates with energies in the continuum may diverge or not in the thermodynamic limit, depending on the specific values of a α .The first kind correspond to scattering eigenstates, while the second correspond to bound states in the continuum (BIC), as discussed in the main text.In the following discussion, we disregard these states and focus only on the unbounded eigenstates.
We now consider the correlations of degenerate eigenstates in the initial state |ψ 0 = c † |FS .We find that where we have dropped the argument E in the denominator, and where Σ < (z) is defined analogously to Eq. (D5), with the sum restricted to k ∈ BZ such that ω k < E F .We want to understand the behaviour of this correlator in the thermodynamic limit.As noted above, the denominator in this expression will diverge.If E > E F , then the divergent term which dominates Σ(E) is not included in Σ < (E), and therefore the numerator in Eq. (D9) remains finite and, accordingly, the correlator vanishes.On the other hand, if E < E F , the numerator is dominated by the divergent term and hence a ), this implies that the correlator vanishes for α = β.Conversely, for α = β, the numerator and denominator diverge in the same way, so that finally for unbounded eigenstates in the thermodynamic limit, which is what we set out to show.
the finite chain are standing waves.The particular standing wave eigenmode at energy ∆ thus forms the localised bath component of a BIC, and by Eq. (E1), this BIC only has amplitudes on the two impurities which "cut out" the finite chain from the lattice.In particular, this implies that in the Dicke regime there emerge N − 1 such BICs, whose explicit form can be calculated, after normalisation, as [114] |φ  (48).
We now focus on the case of the LTA total impurity occupation, for which O = X in Eq. (E2).Using the wavefunction (E4), we can compute where ∓ refers to cos(k ∆ d) = ±1.Since S is tridiagonal and Topelitz, its eigensystem can be obtained analytically and reads [115] λ n = 1 + R cos nπ N , (E6a) with n, m = 1, . . ., N − 1. Notably, V = V T = V † and this transformation will diagonalise any tridiagonal Topelitz matrix with equal values on the off-diagonals.In particular, this includes X, which has eigenvalues To draw a connection to previous work, we can confirm that our explicit analysis of population trapping in terms of BICs is consistent with the resolvent approach used to describe this effect in the single-excitation sector [2,4].For two impurities, we consider (anti)symmetrised impurity operators c ± = (c 1 ± c 2 )/ √ 2, since for a timereversal symmetric bath (ω k = ω −k ) these couple to two sets of orthogonal bath modes with the same dispersion relation ω k as the original bath.Specifically, for emitters at a distance d from each other, H = H + + H − with which satsify {b † k,α , b q,β } = δ αβ δ kq for k ∈ [0, π).Since [H + , H − ] = 0, the two-impurity problem can thus be mapped onto two single-impurity problems.
To extend the resolvent formalism to each of these problems, we define the states |e ± ≡ c † ± |vac and |k ± ≡ b † k,± |vac and the transition amplitudes A ± (t) ≡ e ± | e −iHt |e ± and A k,± (t) ≡ k ± | e −iHt |e ± , which can be calculated from an integral of the form (17) from propagators G ± (z) and G k,± (z) identical to G e (z) and G k (z) in Eqs.(18) except for the substitution Σ e (z) → Σ ± (z) ≡ Σ e (z) ± Σ(z, d) and g → g k,± .By a derivation analogous to the one for the single-impurity model, the collective decay of the fully inverted state is described by (E13) The long-term limit of this expression is determined by the real poles of G ± (z) and G k,± (z).This again includes BS outside the energy band, given by the real solutions to the pole equations z − ∆ − Σ ± (z) = 0, however we neglect their contribution here.This leaves only the pole ω k of G k,± (z), as well as a new pole at z = ∆, which we can show emerges for G ∓ (z) when cos(k ∆ d) = ±1.The associated residue is given by R as defined in Eq. ( 48) [2], such that the LTA occupation is finally given by (E14) Given Eq. (E10), it is straightforward to see that this agrees with Eq. (E2) in the case where N = 2.

Van Hove singularity (2D)
When searching for BICs in the 2D bath with impurities detuned to the Van Hove singularity, we are faced with the challenge that the entries of the self-energy matrix Σ(∆) will, in general, diverge at ∆ = 4J.We should therefore understand Eq. ( 43) as the statement that Σ(∆) has at least one eigenvalue which vanishes in the limit ∆ → 4J.However, rather than trying to establish which constraints this condition imposes on the impurity positions in general, we note that there is a more intuitive approach for constructing BICs at the Van Hove singularity.
For such a BIC, it is easy to see that the eigenvalue equation H • B |Ψ = ∆ |Ψ amounts to the condition j,j ψ j = 0 ∀ j = j n .

(E15)
We can then use Eq.(E15) to iteratively construct a simple, real wavefunction for the BIC supported by this configuration: starting at some site j adjacent to a vacancy site, we assign ψ j = +1 or ψ j = −1 arbitrarily on this site.To ensure that Eq. (E15) holds for the nearest neighbours of this site, we must then assign amplitudes ψ j = ±1 to its next-nearest neighbours j .We then proceed to enforce Eq. (E15) to those sites in turn, until we arrive ultimately at a checkerboard pattern, where every other site in the rectangle has amplitude ψ j = +1 (which we label j ∈ P ) or ψ j = −1 (j ∈ M ).After normalisation, we are left with a BIC wavefunction where we have defined With reference to Refs.[2,4], we note that the above analysis applies for instance to the diamond configuration of four impurities at sites r 1 = (0, −m), r 2 = (−m, 0), r 3 = (0, m), and r 4 = (m, 0) (shown for m = 3 in Fig. 7a).In this case, |P | + |M | = m 2 and therefore R V = 1 + (gm) 2 /(2J) 2 −1 .Noting also that σ = (1, ±1, 1, ±1), where + (−) refers to the case of odd (even) m.It is then easy to show that the subradiant decay in this configuration reported in Refs.[2,4] is consistent with population trapping in a BIC of the form (E16).
More generally, our analysis can also be applied to extended configurations giving rise to multiple BICs of the form (E16).We focus in particular on the two examples shown in Fig. 5, with N = 6 and N = 8 impurities, respectively, whose wavefunctions can be constructed shown in Fig. 7.It is straightforward to show that these configurations support multiple BICs |φ n with n = 1, 2 and n = 1, 2, 3, respectively.With access to the full BIC wavefunction, we are able to construct the matrices S, O, X and Y (k) explicitly, allowing us to calculate LTA expectation values from Eq. (E2).

FIG. 1 .
FIG.1.Single-impurity ground-state correlations.a) Scaling of the 1D ground-state impurity-bath correlation envelope at fixed EF = 2J and g = 0.2J for varied ∆, computed using the wide-band formula (30) (black solid line) and by numerical integration of Eq. (29) (coloured markers).The short-and long-range scalings(32) are also shown (red dashed lines).b) 2D ground-state impurity-bath correlation for two different ∆ = EF and g = 0.2J, computed along two different directions as described in Appendix C. For reference, two directions are indicated in the left pane of c).We also highlight the algebraic scaling of correlations.c) Deformation of the ground-state bath-bath correlations as given by Eq. (34), for g = 0.2J and various EF = ∆.

FIG. 3 .
FIG. 3. Propagation of excitations in a 1D bath.Evolution of the excitation probability b † j bj t (arb.units) during the decay of a single impurity at fixed k∆ = π/4 and g/J = 0.4 for an initially empty bath (left) and for kF = π/10 (middle) and kF = 9π/10 (right), corresponding with ∆ > EF and ∆ < EF , respectively.

FIG. 4 .
FIG. 4. Collective decay.a) Collective decay of N = 5 equidistant impurities into an initially empty 1D bath at ∆ = J and g = 0.2J for varied d associated with the same value of cos(k∆d).For reference, we also plot the master equation dynamics (black solid) and the independent decay of N = 5 impurities (grey solid).b) Fractional decay of N = 5 impurities at different Fermi levels, for ∆ = J, g = 0.4J, and d = 5, and independent decay of N = 5 impurities at same values of EF (grey solid).The LTA occupation, given by Eq(42), is also shown (grey dashed).

FIG. 5 .
FIG. 5. BIC population trapping.a) Collective decay of N = 5 impurities in the Dicke regime at ∆ = J, EF = 0, and g = 0.2J for various d.For reference, we also plot the master equation dynamics (black solid) and the independent decay Ne = N e −Γ 0 t of N impurities (grey solid), as well as the LTA value (47) (grey dashed).b) Collective decay at ∆ = 4J and g = 0.2J for the impurity configurations shown in the insets, where black dots indicate a site to which an impurity is coupled and grey shaded areas indicate the localised wavefunctions of the BICs, labelled by |φn .We also plot the LTA occupations calculated in Appendix E (grey dashed).

FIG. 6 .
FIG.6.Experimental considerations.a) LTA impurity occupation for different initial bath states at temperature T , at EF = 1.5J and g = 0.4J, as a function of ∆.For reference, we also reproduce the T = 0 case from Fig.2c.As in Fig.2c, grey shaded areas indicate the regions 0 < ∆ < EF (light) and EF < ∆ < 4J (dark).b) Adiabatic ground-state preparation for a single-impurity 1D system with ∆ = EF = 2J and different system sizes, with the coupling linearly turned on from g = 0 to g = 0.4J over a time τ .For these parameters, ξ1D = 20 in the target ground state.We plot the infidelity 1 − F , where F denotes the overlap between the adiabatically varied state and the target ground state, as a function of τ .We also plot the impurity-bath correlations for a finite system (L = 200) after 200 tunneling times in the protocol, showing the cross-over from logarithmic to algebraic decay (inset).

≶
mn , and another one stemming from the Dirac delta distribution, Γ ≶ mn /2.The fact that the dispersion relation is symmetric under the change k → −k guarantees that both J , comparing Eqs.(B8) and (B9) with the definition of the self-energy function, Eq. (4), it becomes clear that Eq. (7) holds.Putting everything together, ρI = m, n

F
|rj | f (|r j |/ξ 1D )] ≤ |f (|r j |/ξ 1D )|, we can simply define the envelope of the correlations as F (|r j |) = |f (|r j |)|.2.2D bath a. Efficient numerical computation , 1)}, whereas T (BZ) is a scaled and rotated version of it.Since the integrand is periodic (invariant under translations by reciprocal lattice vectors) we could have considered instead a different integration region, BZ , the square with vertices 2π × {(0, 0), (m, n), (−n, m), (m − n, m + n)}, which has an area (m 2 + n 2 ) times that of the original BZ.The transformed integration region T (BZ ) coincides with BZ.Similarly, we can do the same change of variables to express the self-energy function asΣ(z, r j ) = g k .(C8)Assuming m ≥ n > 0, calling x = e ik and y = e ik ⊥ , we can integrate out k ⊥ using the residue integration technique, obtaining Σ(z, r j ) = where n = 1, . . ., N −1, |e ± n = (|e n ±|e n+1 )/ √ 2 with ± referring to the case where cos(k ∆ d) = ∓1, and R defined in Eq.

FIG. 7 .
FIG.7.BIC wavefunctions for the extended emitter configurations from Fig.5b, uniquely defined up to a global sign.Black dots indicate sites to which the impurities are coupled (i.e.vacancy sites) and we highlight a positive (negative) nonzero amplitude ψj on a bath site j in red (blue).