Kondo QED: The Kondo effect and photon trapping in a two-impurity Anderson model ultra-strongly coupled to light

The Kondo effect is one of the most studied examples of strongly correlated quantum many-body physics. Another type of strongly correlated physics that has only recently been explored in detail (and become experimentally accessible) is that of ultrastrong coupling between light and matter. Here, we study a system which we denote as"Kondo QED") that combines both phenomena, consisting of a two-impurity Anderson model ultra-strongly coupled to a single-mode cavity. While presented as an abstract model, it is relevant for a range of future hybrid cavity-QED systems. Using the hierarchical equations of motion approach we show that the ultrastrong coupling of cavity photons to the electronic states (impurity) noticeably suppresses the electronic Kondo resonance due to the destruction of many-body correlations of the Kondo cloud. We observe this transfer of correlations from the Kondo cloud to the cavity by computing the entropy and mutual information of the impurity-cavity subsystems. In addition, in the weak lead-coupling limit and at zero-bias, the model exhibits a ground-state photon accumulation effect originating entirely from counter-rotating terms in the impurity-cavity interaction. Interestingly, in the strong lead-coupling limit, this accumulation is ``Kondo-enhanced'' by new transition paths opening when increasing the hybridization to the leads. This suggests a new mechanism for the generation of real photons from virtual states. We further show that the suppression of the Kondo effect is stable under broadening of the cavity resonance as a consequence of the interaction to an external bosonic continuum. Our findings pave the way for the simultaneous control of both the Kondo QED effect and a photon accumulation effect using the ultrastrong coupling of light and matter.


I. INTRODUCTION
Understanding the properties of strongly correlated open quantum systems remains one of the significant challenges in quantum many-body physics, with applications in quantum computation [1], machine learning [2], quantum optics [3,4], and condensed matter physics [5].The Kondo resonance, arising from the strong quantum correlations formed between magnetic impurities and the surrounding electrons, has not only provided a testing ground for fundamental theories, but also for quantitative comparisons between theory and experiments [6].
Quantum dots (or single molecules [7,8]) are often used as controllable impurities and can be engineered to manifest the Kondo effect [9][10][11][12][13][14].In addition, they are promising for a range of technological applications, like single-electron transistors [7,[15][16][17].Importantly, for our purpose, it has been demonstrated experimentally that both charge and spin degrees of freedom can be coupled to microwave photons [18][19][20][21].So far, electronic systems (ESs), like quantum dots or impurity spin, operating in the Kondo regime in concert with electron-photon in-teractions [22,23] have offered a way to non-invasively probe quantum correlations in fermionic many-body systems through the phases and amplitudes of the photonic signals [24,25].In these studies the electronic properties are largely unaffected by the cavity photons due to the weak electron-photon interaction.
At the same time, it has been shown that light-matter coupling can be tuned to be on the same order of magnitude as the bare frequencies of the isolated systems [26][27][28].In this ultra-strong coupling (USC) regime, virtual processes which simultaneously create or annihilate both light and matter excitations (usually neglected by the rotating wave approximation in the weak and strong coupling regimes [26,27]) become important.Interestingly, these processes, enabled by the so-called counterrotating terms in the Hamiltonian, induce an hybridization between light and matter even in the ground state, which becomes dressed by virtual excitations [29].This allows for the emergence of counterintuitive phenomena in various fields, such as quantum optics [29][30][31], transport [32,33], chemistry [34][35][36], and condensed matter [37,38].USC can also be realized in the context of open quantum systems [39].Additional novel potential applications have been explored in relation to quantum information processing [40,41], quantum memories [42], quantum plasmonics [43], and quantum thermodynamics [44,45].
A. Kondo QED: An overview In this section, we start with an intuitive overview and explanation of our results.In this work, we provide a qualitative description of how Kondo resonance behaves in the presence of light-matter interactions, which we term Kondo Quantum Electrodynamics (QED).In particular, we are interested in what occurs when those interactions are allowed to be of order of other system energies (the USC regime).
We investigate two scenarios: (i) a single-mode cavity, and (ii) a bosonic continuum, both ultra-strongly coupled to a two-impurity electronic system (ES).In both cases the ES is sandwiched between two fermionic environments that are designed to interact solely with the lower energy impurity |g [see Fig. 1 (a)].
Each impurity can have the following four electronic configurations: a vacant impurity with no electrons, two states where the impurity is occupied by an electron with either spin up or spin down, and a state featuring double occupancy (with both a spin-up and a spin-down electron).
When the coupling between the ES and the leads is large, strong correlations are established between the impurity and the electrons in the leads at the Fermi energy, allowing them to form spin-antisymmetric states.The 'bare' basis of single-and double-electron occupancy described above are then insufficient to describe this regime.In addition, these correlations cause the emergence of a zero frequency peak in the density of states (DOS) of the electronic system, as shown in Fig. 1 (b).This peak is a well-known spectroscopic signature of the Kondo effect [46][47][48][49], and it has been experimentally detected [50].
Because of the many-body nature of the fermionic leads, one must resort to sophisticated numerical methods to describe this parameter regime.One such method, the hierarchical equation of motion (HEOM) approach, allows us to calculate the density of states (DOS) in a numerically exact manner, accounting for strongly non-Markovian bath correlations without resorting to perturbative approximations [51][52][53][54][55].This means that the zero-frequency Kondo peak, which describes the many-body correlations between electrons in the lead and the electron system (ES), can be accurately captured by the HEOM method.
In order to characterize the impact of the light-matter interactions in the USC limit on the Kondo effect, which will also hybridize with the impurity states, we employ the HEOM method in two ways: first, by explicitly including the light degrees of freedom in the system part of the model (for the cavity limit), and second, by including them in the so-called auxiliary HEOM degrees of freedom (for the continuum limit).From this method we calculate the DOS and find that the hybridization between light and matter (emerging as a consequence of the ultrastrong coupling of the cavity field with the ES) can reduce the electronic Kondo resonance, by suppressing the correlations with the Kondo cloud.This result can be understood as a competition between ES-cavity and ES-leads hybridizaton.Stronger correlations between the system and the cavity intuitively decrease the electron availability to form delocalized states with the leads, thereby modulating the Kondo effect.
Furthermore, in the Ultrastrong Coupling (USC) regime, the presence of quantum-fluctuations in the lightmatter ground state [26] enables new transition paths which result in a steady-state photon accumulation effect in the cavity.This intriguing effect manifests under both weak and strong matter-lead coupling conditions.
When the system and the leads interact weakly, this photon accumulation effect is enhanced as the lightmatter coupling increases.This can be intuitively understood using an effective master equation valid in the limit of small system-lead-coupling.This clearly demonstrates how, in this weak-coupling limit, increasing lightmatter coupling leads to increased virtual photon accumulation due to larger ground-state light-matter hybridization [32].The trapping (conversion from virtual to real) of those photons occurs due to an inherent lightmatter decoupling mechanism which activates when both impurities are occupied.
However when the lead-system interaction enters the non-pertubative Kondo regime, a different trend is observed.We first see a larger overall photon accumulation magnitude, but with a counter-intuitive suppression of the photon accumulation effect taking place as the lightmatter coupling is increased.This suppression can be traced back to the cavity-induced decoupling effect which reduces the importance of higher-order transitions enabled by the combination of system-lead and light-matter hybridization.
In addition, another surprising feature is that despite being triggered by the counter-rotating terms in the light-matter Hamiltonian, the photon-accumulation effect appears for relatively weak light-matter couplings and, as mentioned above, its magnitude can be "Kondoenhanced" in the strong matter-lead coupling regime.This feature may provide a new way to indirectly detect the presence of ground-state virtual photons in the Rabi Hamiltonian [56][57][58][59].
In parallel, we also note that the time required to accumulate a photon in the cavity remains sensitive to the intensity of the light-matter and matter-lead interaction strengths.
The rest of this article is organized as follows.In section II.A, we start by introducing our primary cavitybased model.In Sec.III, we lay out our principal findings, specifically highlighting the suppression of the Kondo peak provoked by the ultrastrong coupling to the single-mode cavity, and the photon accumulation effect which is reciprocally influenced by the Kondo effect.In Sec.IV we consider a continuum bosonic bath rather than a single-mode cavity, and show how it impacts the Kondo resonance.In Section V, we wrap up with a conclusion of our findings, their potential applications, and directions for future research.For those interested in more detail, we have included an Appendix.It covers the step-by-step derivation of the DOS using HEOM, an examination of the convergence of the DOS, a derivation of a modified master equation for degenerate systems, a description of the dressed states of the Kondo QED system, and an explanation of the reabsorption of the diamagnetic term in the light-matter interaction.

II. THE MODEL
A. Kondo QED: A two-level ES ultrastrongly coupled to cavity photons Here, we consider a two-level ES ultrastrongly coupled to a single-mode cavity and sandwiched between left and right leads as shown in Fig. 1(a).Such an ultrastrongly coupled ES-cavity can be practically implemented by various near-future cavity QED and circuit-QED setups, such as hybrid superconducting circuits [60], semiconductor quantum wells coupled to a microcavity circuits [61], molecular excitons coupled to a metal-clad microcavity [62], and hybrid solid-state architectures [24,[63][64][65], especially quantum-dot based systems [18,66,67].We model this setup by the Hamiltonian (with = 1 throughout) where the system Hamiltonian describes the ES (H e ), the single-mode cavity (H c ), and their interactions (H ec ).Here, the Hamiltonian of the electronic system is given by where d † nσ creates an electron at level n = g, e with energy n .Here, nnσ = d † nσ d nσ is the electronic number operator with spin σ.The Coulomb repulsion energy U n represents a non-linear effect and requires both spin up ↑ and spin down ↓ electrons to occupy the same state.The two-level ES is near resonance with the fundamental frequency (ω c ) of the single mode a inside the cavity with Hamiltonian In addition, the light-matter coupling between the electronic system and the cavity (c) photons, known as the light-matter coupling, is described by where the coupling constant g ct can originate from a purely transverse (t) engineered interaction [68][69][70][71].It is important to note that we assume g ct and ω c have been implicitly renormalized by the A 2 term [26], as demonstrated in Appendix E. To study the influence of the light-matter USC on the Kondo effect, we further assume a high-quality microwave cavity so that the dissipation to its bosonic environment can be neglected.The leads (labeled by α) are electronic reservoirs described by the Hamiltonian α,k , where c † α,k creates a fermion (f) in the state k of the lead α.Importantly, the electrons in the leads are assumed to couple to only the lowest level |g of the ES.Hence, the interactions between the ES and two separate leads can be characterized by the interaction Hamiltonian The interaction between the electronic system (ES) and the fermionic (f) leads can be fully characterized by the Lorentzian spectral density where Γ α represents the coupling strength between the system and the α-lead with bandwidth W f and chemical potential µ α .

A. Suppression of the Kondo peak
We now compute the DOS of the electronic system.This is a key quantity in describing the Kondo effect [48] and can be engineered to improve electronic device performance.The DOS of this ultra-strongly coupled ES-cavity system can be calculated as in Eq. (A1) using a parity-dependent HEOM based on a recent canonical derivation of the influence superoperator [72,73], see Appendix A for more details.
Here, to have a better resolution of the Kondo effect, we restricted our analysis to the DOS of the lowest state A |g (ω).We further set a large repulsion energy as U n = 15ω c , to avoid any overlap between the Kondo peak and other resonances when increasing the transverse coupling.By varying the cavity coupling from g ct = 0.6ω c to g ct = ω c (deep in the USC regime), the Kondo peak gradually disappears as shown in Fig. 1(d).
We note that, to optimize the memory requirement of the simulation, we truncated the Fock space to three photons.While this is typically insufficient to achieve convergent USC effects, an increase in the truncation only slightly affects the zero-frequency component of the DOS in the Kondo regime.The convergence properties of the Remarkably, in order to observe a noticeable impact on the zero-frequency component of the DOS in the Kondo regime, the strength of the transverse coupling g ct has to be in the deep-strong coupling regime, i.e., it has to be comparable to the cavity resonant frequency ω c .As illustrated in Fig. 2, when we decrease the lead coupling from Γ α = ω c to 0.6ω c , the suppression of the Kondo peak, A |g (0), only occurs as g ct approaches ω c , not at the reduced value of Γ α .Similarly, as the cavity resonant frequency ω c increases, even up to twice its original value, we see a reduction in the Kondo suppression effect, correlated with leaving the USC regime into the strong coupling range.Consequently, maintaining the cavity coupling within the deep ultrastrong coupling regime is also a crucial prerequisite for observing a noticeable Kondo suppression.
Additionally, the suppression of the Kondo effect shows that a potential reduction of the correlations between the system and the leads occurs in the USC light-matter coupling regime.To clarify this, we analyzed the von Neumann entropy of the steady-state reduced density operator ρ ec (∞) = Tr f [ρ T (∞)] by partially tracing over the Hilbert space of the fermionic (f) leads As shown in Fig. 3(a), increasing the light-matter coupling g ct results in reducing the entropy of the system, indicating a decoupling of the ES from the leads.In Fig. 3(b) we also show the von Neumann entropy of the cavity and ES alone, where Combined with the mutual information see Fig. 3(c), these quantities show that correlations between ES and cavity increase when the g ct increases.This is consistent with the mentioned suppression of the Kondo peak in the DOS, providing more evidence that the USC light-matter coupling decouples the ES from the leads (akin to the decoupling seen in the tunnel-coupled two-impurity Anderson model [74]).

B. The photon accumulation
In analyzing the behavior of the cavity photons in the Kondo regime we see that, surprisingly, the average photon number (N ph ) in the cavity is large for small coupling and decreases with increasing light-matter coupling g ct , as shown in Fig. 2(d), which is the opposite of what one might naively expect in USC physics.
To better analyze this effect and its relationship to the Kondo suppression in more detail, it is useful to consider the weak ES-lead coupling regime.In this non-Kondo regime, shown in Fig. 4, we can employ a Born-Markov quantum master equation (BMME) to describe the influence of the leads on the ES-cavity system.This approach uses a Lindblad dissipator written in terms of a decomposition on the different eigenstates |ϕ i of H s with eigenenergies I , i.e., where [•, •] ({•, •}) denotes the commutator (anticommutator).Note that the density operator in this Born-Markov quantum master equation is allowed to contain both even and odd parity, i.e., Here, we assume that there is no coherence in the initial system states (otherwise, the Born-Markov quantum master equation should be modified to take into account degenerate energy levels in the eigenoperator decomposition).The details of the derivations are shown in Appendix C. The Born-Markov quantum master equation provides information on the transitions between different system eigenstates which happen at the rates written in terms of the spectral density of leads J fα and Fermi-Dirac distribution n fα = 1/(e βω + 1), with β = (k B T ) −1 (k B = 1).Assuming an initially empty electronic system, an electron will rapidly enter the system due to the higher potential of the leads relative to the impurity energies.Due to the USC between cavity and ES, this electron can enter the ground (|G ± with (a) N ph ≈ 0.5) and higher (|ϕ ± l=9,10 with N ph ≈ 1.5) photondressed states as shown in Fig. 4(a).
Importantly, these states contain components in which the ES is excited to the impurity 2 (state |e ) and virtual photons are present in the cavity, due to the counterrotating terms in the light-matter interaction.These components allow for a non-zero rate to two-electron states, where both impurity 1 (state |g ) and 2 (state |e ) are occupied, and physical photons populate the cavity, through the paths where |s represents the state empty of photons and electrons.Here, the stationary dressed states can be approximately written as and Note that |↑ (↓), ↑ (↓), n c represents the uncoupled eigenstate containing the electrons with spin up ↑ (spin down ↓) configuration singly occupying the lower |g and higher |e energy levels.Here, f i(j) k is the corresponding amplitude of each uncoupled eigenstate.In this sense, |ϕ i=1,...,4 and |ϕ j=5,...,8 mainly contain around 1 and 3 photons, but, because of the double electron occupation, they are uncoupled from the cavity, see Appendix D for more details.
As shown in Fig. 3(d), in the long time limit of the non-Kondo regime, photons accumulate in the cavity, even for arbitrarily small light-matter coupling in a time which depends on g ct .A larger g ct enhances this ground-state photon accumulation rate as shown in Fig. 4(c).Importantly, this photon is not virtual, and will eventually decay into the electromagnetic environment, allowing for a potential observation of this effect.
As we increase the ES-lead coupling to reach the Kondo regime, this strong ES-lead coupling can allow for higherorder transitions to transient states, such as |G ↑↓ with double occupation in the lower state and |2 with no electrons but two photons, leading to the dressed states |ϕ 11 and |ϕ k=13,14,15 as illustrated in Fig. 4(a).Meanwhile, compared to the non-Kondo regime, such strong ES-lead coupling in the Kondo regime can drastically enhance the ground-state photon accumulation rate, as displayed in Fig. 4(c).As expected, in the Kondo regime, the dynamics of the corresponding low energy states cannot be de-scribed by the Born-Markov master equation, as shown in Fig. 5(a) and Fig. 5(b).This master equation also fails to model the non-perturbative effects causing transitions to higher excited states such as |G ↑↓ , |2 , and |ϕ k=11,14,15 , which did not play a role in the non-Kondo regime, as seen in Fig. 4(a) and Fig. 5(c).
In the long-time limit, shown in Fig. 5(d), single photon dressed states (|ϕ i=1,...,4 ) dominate the steady-state occupation in the non-Kondo regime.The photon occupation of these states increases as we increase g ct (due to an increase in the expected photon number of intermediate transient states) giving rise to the increase of N ph , see Fig. 3(d).However, in the Kondo-regime, at weaker light-matter couplings, we see larger average populations because of access to new transition paths involving threephoton dressed states (|ϕ j=5,...,8 ) and even higher energy dressed states (|ϕ ± k=9,...,15 ), resulting in Kondo-enhanced dressed photon accumulation.At the same time, these pathways are suppressed as we increase the light-matter coupling, thereby isolating the ES-cavity system from the leads (also causing the suppression of the Kondo peak) and reducing the photon accumulation effect, as shown in Fig. 3 (d).

C. Two-impurity Anderson model ultrastrongly coupled to a bosonic continuum
In this section, we generalize the previous analysis based on the resonant interaction to single-mode cavities to allow for a coupling to a continuum of environmental modes.This setting is usually introduced to investigate a richer domain of phenomena such as hybridization with the bath [39], unraveling of multiple excitation bound states [75], harvesting single photons from the vacuum [76], developing robust long-distance entanglement protocols [77], and also to promote advancements in the design of quantum computing and sensing devices [78].Thus, in addition to single-mode cavities, we also explore systems comprised of a two-level ES ultrastrongly coupled to a bosonic continuum and positioned between two leads, as illustrated in Fig 6(a).
The Hamiltonian describing the bosonic continuum (b) fields b j with energy ω b,j is denoted as The interaction between the ES and the bosonic continuum can be characterized by Assuming the coupling is of the transverse type (t), as referred to in Eq.A11, we recast the interaction Hamiltonian into the correlation function given by Eq.A22.Subsequently, we characterize the spectral density of the bosonic continuum using the Drude-Lorentz model, where λ represents the coupling strength between the electronic subsystem (ES) and the bosonic continuum, which functions as a bosonic reservoir with a bandwidth of W b .Assume that W b = ω c as depicted in Fig. 6(b), the bosonic continuum can still have a high amount of energy in the particular mode, which can simulate the singlemode cavity in the previous example but with broadening energy distribution contributed by other modes.This assumption offers an advantage in comparing the effects of the coupled cavity between its distinctive single and multimode characteristics.
To examine the influence of the ultrastrongly coupled bosonic continuum on the Kondo effect, we determine the DOS employing the hybrid HEOM approach, which encompasses both fermionic and bosonic hierarchy (refer to Eq. (A33)).
For the numerical implementation, we truncate the bosonic hierarchy to N b = 4 and set the number of exponents in the bosonic correlation to m max = 10 (see Eq.A24).The convergence properties of the entire DOS concerning the truncation of N b and m max are further shown in Appendix B. Upon implementing a coupling strength of λ = 0.25ω c for the bosonic continuum on ES, there is virtually no impact on the Kondo peak in the DOS, as demonstrated in Fig. 7(a).Moreover, when we increase λ to 0.5ω c , a noticeable suppression of the Kondo effect becomes apparent.
To compare these results with the single-mode cavity case depicted in Fig. 1(b), it is important to introduce a renormalization of the parameter λ to account for the broad nature of the system-bath coupling described by Eq. ( 21).This could be done by mapping, in specific parameter regimes, the overdamped spectral density in Eq. ( 21) to its underdamped version [79,80].
Alternatively, it is also possible to follow a more intuitive route and consider that the continuum can be approximately replaced by a single effective ancillary mode whose coupling strength to the system scales as the residue √ λω c of the spectral density J b (ω) at iω c , see Eq. (E18) in [81].Using this scaling, we observe a remarkable similarity in the Kondo peak suppression (as a function of the effective, normalized coupling strengths to both the bath g ct /ω c and λ/ω c ) between the cases of a single-mode and a continuum bath.These results suggest that, even in the presence of broadening, the resonant interaction to an environmental bosonic mode continues to play a dominant role in suppressing the Kondo correlation between electrons in the ES and the leads.

IV. CONCLUSION
In summary, we have shown that the ultrastrong lightmatter interaction can suppress the Kondo screening and  simultaneously allow for a Kondo-enhanced steady-state photon-trapping effect (via counter-rotating terms and virtual transitions), which is reduced due to lead isolation as the light-matter coupling is increased.
For the latter, our results indicate that, in the Kondo regime, an increase of the ultrastrong coupling to the cavity does not necessarily imply stronger light-matter hybridization.While increasing the lead coupling allows more electrons to participate in the light-matter hybridization, simultaneously increasing the ultrastrong coupling to the cavity generates competition between cavity-induced delocalization across the electronic system (ES) and ES-lead coupling-induced delocalization of electrons between the ES and leads.This competition gives rise to the counterintuitive decoupling effect between the system and the leads which, in turn, further reduces the photon accumulation effect and the suppression of the Kondo correlation.Therefore, our work offers a comprehensive perspective to gain both qualitative and quantitative understanding of the complex interplay between light-matter coupling and Kondo physics at the nanoscale.
We note that any physical cavity inevitably interacts with its own electromagnetic environment resulting in photon loss [82].At the same time, the choice of different gauges needs to be done very carefully [83,84].More general conditions, including non-zero photonic decay, will be considered in future work.Furthermore, the USC Kondo-photon interaction considered in this work can also be combined with different impurity configurations, e.g., coupling of both impurities to leads, additional spinorbit coupling [85], arbitrary mixing longitudinal and transverse couplings of the cavity photons [70], or RKKY interaction [74] to generalize the possible physics observable in this interesting regime, and allow for other potential ways to tailor the competition between the formation of electron-photon dressed states and many-body entangled states in the Kondo effect.The density of states (DOS) can be written in a compact form as [46,49] in terms of the retarded G R nσ (t) and the advanced G A nσ (t) which depend on the system correlation functions as Here, the Heaviside function Θ(t) ensures that causality is properly accounted for.The correlation functions of the system are given by where G e/o (t){ô e/o } represents the propagator with arbitrary parity symmetry depending on whether it is applied on an even ôe or odd ôo parity operator.Here, the density operator (ρ) referred to the total system (T) is considered to be in the long-time limit where the initial density operator is evolved into a correlated steady-state from an initially uncorrelated condition which depends on the thermal (th) equilibrium state of the fermionic (f) bath (ρ th f ), the thermal equilibrium state of the bosonic (b) bath (ρ th b ), and the initial state of the system (s) ρ s (0).The propagator G e/o (t){ρ T (0)} describing the behavior of the quantum system can be obtained by solving the hierarchical equations of motion (HEOM) as follows.
Without loss of generality, we consider a N -level electronic system (ES) coupled to n α metallic leads.In addition, the ES is coupled to a bosonic field, which, generically, can represent a cavity mode [28,32], phonon mode [86][87][88] or surface plasmon [89,90].We consider the Hamiltonian, where H s = H e +H c +H ec is the Hamiltonian of the system (s) including the electronic (e) system (H e ), cavity (c) field (H c ), and their interaction (H ec ).Here, H b represents the Hamiltonian of the bosonic (b) bath.We further assume that the degrees of freedom of the system can also be coupled to their bosonic environment H sb .The Hamiltonian of the N -level ES can be written as in terms of the bare energies n and Coulomb interactions U n .The interaction Hamiltonian between the ES and leads can therefore be written as Additionally, the Hamiltonian describing the cavity (c) fields a with energy ω c is denoted as ) the interaction between the ES and the cavity field with coupling strength g c can be characterized by where Q σ represents the fermionic interaction operator of the ES.The form of Q σ depends on the specific coupling type of the cavity field.For example, one may have longitudinal or transverse coupling [69,70].The specific form of the interaction operator Q σ we consider here is a transverse (t) coupling g ct .Moreover, the couplings of the bosonic (b) environment to the interior degrees of freedom of the system (s) with coupling strength g sb,j (j labels the mode number) can be modeled by where b † j (b j ) represents the creation (annihilation) operator of the bosonic environment, Note that V s refers to the Hermitian ES-interaction operators acting on fermionic degrees of freedom.In the fermionic case, V s must have even parity to be compatible with charge conservation.For the system-bath interactions, we can then write down the interaction Hamiltonian in the interaction picture as In this frame, the system (s) density operator (ρ) rotates (denoted by the tilde) as ρs = U † s (t)ρ s (t)U s (t), where U s (t) = e iHst .In order to derive the HEOM of the system, we begin with using the Liouville-von Neumann equation in the interaction frame which can be integrated to obtain the formal solution By iteratively replacing ρs (t 1 ) with ρs (t) in Eq. (A15), one obtains the Dyson series of the von Neumann equation in terms of the time ordering superoperator T where Here, the two hats ( •) refers to superoperator.For the reduced Dyson series to be solvable, we apply the canonical approach in Ref. [72] which does not require any path integral.The formal expression of Eq. (A16) can be written as where Ĝ(t)[•] is a superoperator which propagates the even-(p = +) or odd-parity (p = −) density operator and can be used to calculate the DOS.Its explicit form [72] is given by in terms of the following fermionic superoperator and the following bosonic superoperator [39,54] where p = ∓ represents the projection on the even or odd sector.Here, the [•, •] − and [•, •] + denote the commutator and anticommutator, respectively.For simplicity, we define ν to denote the presence (ν = 1) or absence (ν = −1) of a Hermitian conjugation.Moreover, we define ν = −ν.As it can be seen from the previous expressions, the effects of the fermionic and bosonic environments on the system are completely encoded in the correlation functions which, in the fermionic case, depend on the spectral density Analogously, in the bosonic case, they depend on the spectral density J b (ω) = 2π j g 2 sb,j δ(ω − ω j ) and the Bose-Einstein distribution n eq b (ω where k B is the Boltzmann constant and T fα (T b ) represents the absolute temperature of the α-fermionic (bosonic) bath.Non-zero chemical potential (µ α = 0) in the α-fermionic bath can account for non-equilibrium physics.
To proceed in the derivation of the HEOM, the bath correlation functions will be expressed as a sum of exponential terms, which allows to define an iterative procedure.Specifically, based on some spectral decomposition schemes, such as the Matsubara spectral decomposition [91] or the Padé spectral decomposition [92], the correlation functions of both fermionic and bosonic environments can be written as a sum of exponentials as and where τ = t 1 − t 2 .However, in order to obtain a closed form for the HEOM, the bosonic correlation function has to be further decomposed into its real and imaginary part unless χ m = χ * m .Here, N R refers to the number of exponentials used to obtain C R b (τ ).Similarily, for C I b (τ ).By plugging the fermionic, [Eq.(A23)], and bosonic correlation functions [Eq.(A25) and Eq.(A26)] in Eq. (A18) and taking the time derivative, one can obtain the explicit form of the superoperators Âν nσ (t), Bν αlnσ (t), Pm (t) and K(t) which leads to where Ps [•] = P s [•]P s represents the parity superoperator of the system with Here the symbols j and q involve the multi-index {α, l, n, σ, ν} and {u, m}, respectively.Therefore, we can sequentially define the fermionic and bosonic part of the first-tier auxiliary superoperator propagator as Ĝ(1,0) and Ĝ(0,1) To obtain the higher-tier auxiliary superoperator propagators, we must repeatedly take the time derivative of the different-tier auxiliary superoperator propagators.First, the derivative of the first-fermionic-tier and first-bosonic-tier auxiliary superoperator propagators can be expressed, respectively, as and The influence of different cavity photon numbers on the Kondo effect is subtle.In our simulations we found that a truncation in the Fock space at N ph = 3 photons achieves a good convergence accuracy at the price of a reasonable computational cost.
In the bosonic continuum scenario, we also evaluate the convergence of the DOS using the same approach.We examine the difference ∆N b(i,i+1) between the DOSs of the ith tier and its neighboring tier (i + 1)th in the bosonic HEOM, and the difference ∆m j,j+1 between the jth and its neighboring (j + 1)th Padé approximants of the bosonic continuum correlation function in the bosonic HEOM.The extremely small deviation values of ∆N b(4,3) and ∆m (10,9) indicate a high level of convergence when N b = 4 and m max = 10 as shown in Fig. 11.

Appendix C: Modified master equations for degenerate systems
To model a degenerate system in the USC regime, we here use the fermionic influence superoperator to derive a modified master equation, which takes the degenerate energy levels into account in its form for the Lindblad dissipators.We first decompose all operators into the eigenbasis |ϕ k and |ϕ l of the system Hilbert space in the interaction picture and define the operator A ν σ,kl (ω) as where k and l are the eigenenergies corresponding to the states |ϕ k and |ϕ l , respectively.With the definition of τ = t 1 − t 2 , the expression in Eq. (A17) for the related density operator ρs (t) can be rewritten as where the fermionic influence superoperator is given by We now define ∆ω = ω − ω as the difference between the eigenenergies ω and ω.The nonsecular terms in terms of exp(±i∆ωt 1 ) can be neglected due to their fast oscillations when ω = ω, Hence, by applying this secular approximation on the influence superoperator, one obtains   Here, the information about degenerate transitions is encoded in the decomposed operator, A ν σ (ω) = kl A ν σ,kl (ω).By neglecting the Lamb shift due to its small value in Eq. (C5), one can obtain the modified master equation in Lindblad form as The superindex eq refers to equilibrium.In this section we have derived the modified master equations for degenerate system via the canonical approach [72].

Figure 1 :
Figure 1: (a) An electronic system, described by two states having energy g = −6ωc and e = −5ωc, coupled to a high quality factor single-mode cavity with strength gct and sandwiched between left and right leads with coupling strengths ΓL and ΓR, respectively.On the right, we represent the four possible electronic states of each impurity |g and |e : (bottom) two states where the impurity is occupied by an electron with either spin up or spin down, (middle) an empty state with no electrons, and (top) a double-occupied state with both spinup and spin-down electrons.(b) Illustration of the Kondo regime, characterized by strong correlations between system electrons in |g and the leads electrons at the Fermi level.As a consequence, a Kondo peak at the Fermi level (zero frequency) in the DOS appears which requires non-perturbative methods, such as the HEOM approach.(c) Suppression of the Kondo peak around zero frequency for the equilibrium DOSs (µL = µR = 0) of the ES due to increasing the transverse light-matter coupling gct from 0 (black dotted) to ωc (purple dash-dot) at temperature T = 0.5ωc.The Coulomb repulsion between spin-up and spin-down electrons is set to Un = 15ωc.The coupling strengths to both leads are ΓL = ΓR = Γ = ωc with bandwidth W f = 10ωc.The truncation of the Fock space dimension, of the HEOM tiers, and of the Padé series, are set to N b = 3, Nc = 3, and u = 5, respectively.

Figure 2 :
Figure 2: (a) The zero-frequency component A |g (0) of the DOS as a function of the transverse coupling strength gct, with Γ = 0.6ωc (red star markers) and Γ = ωc (blue square markers).(b) The zero-frequency component of the DOS, A |g (0), plotted against the cavity frequency, ωc.The suppression of the Kondo peak diminishes as ωc increases.

Figure 3 :
Figure 3: Correlations between subsystems.As a function of the light-matter coupling strength gct, we show (a) the Von Neumann Entropy (Sec) (b) the cavity (Sc) and ES (Se) entropy, (c) the mutual information (Iec), and (d) the average photon number (N ph ).The scaling of the correlations depends on the coupling strength to the leads (Kondo regime in red and non-Kondo regime in blue).
Figure 4: (a) On the left, we show an energy diagram for the uncoupled system (gct = 0, Γ = 0).The cavity states |i describe i photons while the ES states (|e and |g ) can be occupied by electrons with arbitrary spin configuration.As the USC of cavity to ES occurs, the tunneling of an electron from the leads causes a transition from the empty state |s to the intermediate transient states |G ± and |ϕ ± l=9,10 (red arrow).Shortly afterwards, another entering electron will participate in this process, further inducing the transition to the stationary dressed states |ϕi=1,...,4 (blue arrow) and |ϕj=5,...,8 (blue dashed line arrow).By increasing the ES-leads couplings, there is a transition from the non-Kondo to the Kondo regime where the dressed states with higher energy, i.e., |G ↑↓ and |ϕ k=11,...,15 , can be excited.(b) The short-time dynamics of different dressed states obtained by solving the HEOM are shown with markers while the solid-lines represent the solutions of the Born-Markov quantum master equation (BMME).The two methods are in good agreement in the non-Kondo regime (Γ = 0.01ωc).(c) Evolution time needed to accumulate an expected number of photons equal to one as a function of the light-matter coupling gct in the non-Kondo (Γ = 0.01ωc) and the Kondo (Γ = ωc) regimes.

Figure 5 :
Figure 5: Short-time quantum dynamics showing the failure of the BMME (solid curves) in comparison with the HEOM (markers) in the Kondo regime.Due to the Born-Markov approximation, the BMME can be seen to underestimate the populations of |G ± and |ϕ ± l=9,10 in (a) but to overestimate the populations of |ϕi=1,...,4 and |ϕj=5,...,8 in (b).(c) shows that the populations of dressed states with higher energy can only be obtained by solving a non-perturbative method such as the HEOM.(d)shows the populations of stationary dressed states as a function of gct in both the non-Kondo and Kondo regime.

Figure 6 :
Figure 6: (a) An electronic system, characterized by two states with energies g = −6ωc and e = −5ωc, is coupled to the jth mode in a bosonic continuum with strength gct,j.Simultaneously, this electronic system is situated between left and right leads, with coupling strengths ΓL and ΓR, respectively.(b) The spectral density of the bosonic continuum is characterized by the Drude-Lorentz model with its peak at ω = ωc.
We acknowledge Stephen Hughes and Xiao Zheng for helpful suggestions and discussions.N.L. acknowledges partial support from JST PRESTO through Grant No. JPMJPR18GC, and the Information Systems Division, RIKEN, for the use of their facilities.M.C. acknowledges support from NSFC (Grants No. 12050410264 and No. 11935012) and NSAF (Grant No. U1930403).F.N. is supported in part by: Nippon Telegraph and Telephone Corporation (NTT) Research, the Japan Science and Technology Agency (JST) [via the Quantum Leap Flagship Program (Q-LEAP), and the Moonshot R&D

Figure 7 :
Figure 7: (a) The DOS of the ES ultrastrongly coupled to a bosonic continuum with λ = 0, λ = 0.25ωc, and λ = 0.5ωc is represented by black dotted, light blue dashed, and red solid curves, respectively.At λ = 0.5ωc, the decrease in the zero-frequency DOS peak, A |g (0), indicates the suppression of the Kondo effect.(b) The height of the zero-frequency DOS peak, A |g (0), is plotted as a function of the coupling strength.Blue-square and red-star markers represent the coupling of the ES to a single-mode with coupling strength gct and to a bosonic continuum with rescaled coupling strength √ λωc, respectively.The Kondo effect is suppressed as the effective coupling to the bosonic environment is increased.
Appendix A: Canonical derivation of the density of states

Figure 9 :
Figure 9: The deviations ∆l of DOSs computed in lmax = 2 from lmax = 3, lmax = 3 from lmax = 4 and lmax = 4 from lmax = 5.The difference between five terms and four terms ∆l (5,4) , four terms and three terms ∆l (4,3) , three terms and two terms ∆l (3,2) of the Padé expansion for the DOS results are presented in red solid, blue dashed, and cyan dotdashed curves, respectively.It also shows that the results are convergent with lmax = 5. (a)(c) and (b)(d) correspond to gct = 0.6ωc and ωc, respectively.Meanwhile, (c) and (d) display magnified plots near the zero-frequency DOS.

Figure 10 :Figure 11 :
Figure 10: The differences ∆N ph between N ph = 3 and N ph = 2 (∆N ph(3,2) ) as well as N ph = 2 and N ph = 1 (∆N ph(2,1) ) for the DOSs are presented in red solid, blue dashed curves, respectively.It shows that the result is convergent around ω/ωc = 0 with N ph = 3. (a)(c) and (b)(d) correspond to gct = 0.6ωc and gct = ωc, respectively.Furthermore, (c) and (d) showcase the magnified plots in the vicinity of the zero-frequency DOS.