Unified Light-Matter Floquet Theory and its Application to Quantum Communication

Periodically-driven quantum systems can exhibit a plethora of intriguing non-equilibrium phenomena that can be analyzed using Floquet theory. Naturally, Floquet theory is employed to describe the dynamics of atoms interacting with intense laser fields. However, this semiclassical analysis can not account for quantum-optical phenomena that rely on the quantized nature of light. In this paper, we take a significant step to go beyond the semiclassical description of atom-photon coupled systems by unifying Floquet theory with quantum optics using the framework of full-counting statistics. This is achieved by introducing counting fields that keep track of the photonic dynamics. This formalism, which we dub ``photon-resolved Floquet theory"(PRFT), is based on two-point tomographic measurements, instead of the two-point projective measurements used in standard full-counting statistics. Strikingly, the PRFT predicts the generation of macroscopic light-matter entanglement when atoms interact with multimode electromagnetic fields, thereby leading to complete decoherence of the atomic subsystem in the basis of the Floquet states. This decoherence occurs rapidly in the optical frequency regime, but is negligible in the radio frequency regime. Our results thus pave the way for the design of efficient quantum memories and quantum operations. Finally, employing the PRFT, we propose a quantum communication protocol that can significantly outperform the state-of-art few-photon protocols by two orders of magnitude or better. The PRFT potentially leads to insights in various Floquet settings including spectroscopy, thermodynamics, quantum metrology, and quantum simulations.

Periodically-driven quantum systems can exhibit a plethora of intriguing non-equilibrium phenomena that can be analyzed using Floquet theory.Naturally, Floquet theory is employed to describe the dynamics of atoms interacting with intense laser fields.However, this semiclassical analysis can not account for quantum-optical phenomena that rely on the quantized nature of light.In this paper, we take a significant step to go beyond the semiclassical description of atom-photon coupled systems by unifying Floquet theory with quantum optics using the framework of full-counting statistics.This is achieved by introducing counting fields that keep track of the photonic dynamics.This formalism, which we dub "photon-resolved Floquet theory" (PRFT), is based on two-point tomographic measurements, instead of the two-point projective measurements used in standard full-counting statistics.Strikingly, the PRFT predicts the generation of macroscopic light-matter entanglement when atoms interact with multimode electromagnetic fields, thereby leading to complete decoherence of the atomic subsystem in the basis of the Floquet states.This decoherence occurs rapidly in the optical frequency regime, but is negligible in the radio frequency regime.Our results thus pave the way for the design of efficient quantum memories and quantum operations.Finally, employing the PRFT, we propose a quantum communication protocol that can significantly outperform the state-of-art few-photon protocols by two orders of magnitude or better.The PRFT potentially leads to insights in various Floquet settings including spectroscopy, thermodynamics, quantum metrology, and quantum simulations.

I. INTRODUCTION
In recent years, periodic driving has emerged as a powerful tool for the coherent control of many-body systems.This has led to the realization of novel quantum phases of matter like dynamical topological states [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] and discrete time crystals [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] as well as breakthroughs in applications like spectroscopy [33][34][35][36], metrology [37][38][39], and quantum simulation [40][41][42][43][44][45][46][47][48][49][50][51].These non-equilibrium quantum systems are generally analyzed using Floquet theory -a method first developed by Jon Shirley in 1965 [52].Interestingly, Floquet theory is also employed to investigate the dynamics of quantum systems interacting with a single-frequency quantum field.A particularly striking example of this is the case of atoms interacting with an intense laser field [53,54].Despite decades of extensive progress in quantum optics, it remains extremely challenging to employ a completely quantum mechanical treatment to this situation, due to the large number of photons involved [55].Often, a semiclassical approach is used instead, where the photon fields are assumed to In this paper, we take a significant step beyond the semiclassical Floquet theoretic description of lightmatter interactions by developing a framework dubbed "photon-resolved Floquet theory" (PRFT).The PRFT bridges Floquet theory and quantum optics by introducing full-counting statistics (FCS) of photons in the semiclassical description of the quantum system.Originally developed in the context of quantum optics [66][67][68], FCS is a powerful method that has been employed to study mesoscopic transport [69][70][71][72][73][74][75][76], quantum dots [77][78][79][80][81], spin chains [82,83], spontaneous photon emission [84][85][86][87], thermodynamics [88,89], and the entanglement entropy of noninteracting fermions [90][91][92].However, we must exercise caution in applying the standard FCS framework to coherently laser-driven systems described in the previous paragraphs.This is because FCS is inherently based on two-point projective measurements [93].Unfortunately, such measurements destroy the coherent photonic states, thereby rendering the Floquet description of the matter  The standard FCS extracts the statistical information based on two-point projective measurements at times t0 and t1.The counting statistics in the PRFT relies on two-point tomographic measurements of two independent batches.One of which is used to determine the photon statistics pn(t) at time t0, and the other determines the statistics of the time-evolved state t1 without being measured at t0.(b) Illustration of the light-matter entanglement in a two-mode Rabi model.As explained in details in Secs.II F and III B, the two-level system controls the photon transport between the two photonic modes, which leads to entanglement in the Floquet basis.Consequently, this entanglement gives rise to decoherence, which is incorrectly described by the standard Floquet theory, but correctly predicted by the PRFT.(c) Overview of some important implications and applications of the PRFT.
system invalid.
The PRFT developed here can track the photonic dynamics, without destroying the Floquet description of the matter system.This is achieved by introducing a framework based on two-point tomographic measurements of the photonic field, instead of the usual projective measurements.The distinction between two-point projective and two-point tomographic measurements is illustrated in Fig. 1(a).Formally, the PRFT introduces counting fields into the semiclassical equations of motion, leading to a dynamical cumulant generating function.This in turn enables us to investigate the redistribution of photons amongst the Fock states and gain a clear picture of the photonic dynamics.We note that while some recent interesting works have investigated the photonic dynamics of driven systems in Sambe space [94][95][96], these approaches scale exponentially with the number of frequency modes; the PRFT does not suffer from this limitation.The PRFT thus enables us to go beyond the previous approaches that connects Floquet theory to cavity dressed states as in Ref. [53], since those approaches do not provide correct results for the photon statistics.
Based on analytical derivations and extensive numerical benchmarking, we demonstrate that the PRFT is valid for coherent and number squeezed states with a moderate mean photon number n > 500 and photon standard deviation as low as σ = 4.The PRFT thus covers all types of moderate and highly occupied photonic fields in experiments, such as radio frequencies, microwaves and lasers, i.e., all driving fields for which standard Floquet theory is believed to be valid.In other words, the PRFT approach here assumes the same level of generic conditions as the standard Floquet theory, but is found remarkably capable of capturing the dynamics of matter system as well of its photonic counterpart.
We employ the PRFT to analyze multimode driving, and discover that this leads to macroscopic light-matter entanglement at long times due to the matter-systemcontrolled transport of photon between distinct modes.This light-matter entanglement causes a complete decoherence of the matter system in the basis of the Floquet states.This effect is depicted for a pardigmatic twomode Rabi model in Fig. 1(b).The PRFT thus provides a quantum-optical interpretation of Floquet states as the decohering basis for the matter system.As the standard Floquet theory is unable to describe this fundamental decoherence effect, it will incorrectly predict the dynamics of the matter system in general.The PRFT thus demonstrates that even in the semiclassical regime, fundamental physical implications of the standard Floquet theory are not understood and require further investigation.In this context, it is worth noting that light-matter entanglement can arise even in single-mode models due to the photon shot noise [97][98][99][100][101].However, this shotnoise induced entanglement has a much smaller impact on the photonic dynamics compared to the transportrelated entanglement revealed by the PRFT.A detailed analysis of these issues is presented in Sec.III A.
Furthermore, the transport-entanglement-related decoherence effect has far-reaching experimental consequences.
In particular, we demonstrate that the quantum-optical coherence time is reasonably short (a few ms) for typical optical fields used in experiments, but it can be very long for radio-frequency driving.This implies that the radio-frequency regime is optimal for realizing quantum memories and quantum operations.Furthermore, we argue that quantum time crystals provide a powerful platform for realizing quantum memories irrespective of the driving frequency.Intriguingly, we demonstrate that the light-matter entanglement described by the PRFT can be deployed in a quantum-communication protocol that is intrinsically robust against photon loss.In particular, we demonstrate that using coherent light, it is possible for the quantum state transfer rate to reach the 0.1 KHz regime over 500 km, thereby far exceeding the Hz regime that is predicted in current theoretical protocols [102].Our analysis thus demonstrates that the PRFT can play a pivotal role in the development of future quantum technologies.
This paper is organized as follows: In Sec.II, we introduce the theoretical framework of the PRFT.In Sec.III, we apply the PRFT to mulitmode quantum Rabi models for benchmarking, and investigate the light-matter entanglement.In Sec.IV, we discuss the experimental verification of the theory and implications for quantum memories, quantum time crystals, and other quantum applications.In Sec.V, we devise a quantum-communication protocol employing the PRFT framework.Finally, we summarize the main findings of this paper and discuss avenues for future research in Sec.VI.

II. PHOTON-RESOLVED FLOQUET THEORY
In this section, we introduce the basic ideas and main results of the PRFT.We emphasize that even though we primarily analyze Floquet systems in this paper, the formalism can also be used to analyze aperiodically driven systems.For a more detailed analysis, we refer the reader to Appendix A.

A. System
We consider the following generic Hamiltonian describing a matter system interacting with a multimode photonic field: where k denotes the different photonic modes, âk are annihilation operators quantizing these modes and H k acts on the matter system.The light-matter interaction strength is parameterized by g.The dynamics of this system can be determined by representing the photonic operators with Fock states, which are the eigenstates of the occupation operators Nk = â † k âk , i.e., Nk |n⟩ k = n |n⟩ k .However, for typical laser fields, the photonic modes are highly occupied, such that an analytical or numerical treatment becomes infeasible for more than two modes.
Alternatively, one can employ a semiclassical description of the system by assuming that it is initially in the state where α k e iφ k are coherent states of the photon operators âk with real-valued amplitudes α k > 0 and phases φ k ∈ [0, 2π), and the state |ϕ(t 0 )⟩ is the initial state of the matter system [52,54].In this semiclassical limit, we can substitute âk → α k e iφ k −iω k t in Eq. ( 1) such that we obtain the corresponding semiclassical Hamiltonian where we have introduced the effective light-matter interactions g k = gα k .This description is valid as long as the back action of the quantum system on the photonic field is negligible, i.e., if g ≪ ω k α k .Nevertheless the impact of the photonic field on the matter system can be large because of the product gα k .The semiclassical approach is thus valid for large mean occupation numbers ⟨ Nk ⟩ = α 2 k .For a single photonic frequency ω k or for commensurate frequencies (where all ω k ′ /ω k are rational numbers), the semiclassical Hamiltonian is time periodic H(t) = H(t + τ ) with a driving period τ .Under these conditions one can apply the celebrated Floquet theory to analyze the system [52].Unfortunately, this effective semiclassical description loses the microscopic information about the photonic field.The PRFT resolves the problem by introducing counting fields into the semiclassical description.
Before proceeding further, we would like to point out that the transition from the quantum Hamiltonian in Eq. ( 1) to the mean-field Hamiltonian in Eq. (3) corresponds to the transition from Fock space to Sambe space.In the latter, the photonic operators are replaced by their unbounded counterparts which can be formally derived by a Fourier transformation of the Schrödinger equation determined by the Hamiltonian in Eq. (3).It is worth noting that the photon number dependence of the matrix elements of â † k in Fock space can lead to a light-matter entanglement even for coherent states, which has been extensively investigated [97][98][99][100][101].This is a consequence of the photon-number dependent light-matter interaction g k (n k ) ∝ √ n k , which results in a photon-number dependent dynamics of the matter system.As this effect is induced by the finite photon number uncertainty of the coherent light fields, we will refer to it as photon shot-noise entanglement.We emphasize that this form of light-matter entanglement has a minor effect on the photonic dynamics.

B. Photon-resolved time evolution
The time-evolution operator U(t) corresponding to the Hamiltonian in Eq. ( 3) does not contain the information about the microscopic state of the photonic fields.In order to track the photonic dynamics, we introduce realvalued counting fields χ k ∈ [0, 2π), and define the generalized time-evolution operator, where U (t) is the time-evolution operator corresponding to the Hamiltonian in Eq. ( 1).This transformation implies that the annihilation (creation) operators transform as âk → âk e iχ k (â † k → â † k e −iχ k ), leading to the following generalized time-evolution operator in the semiclassical limit, where χ = (χ 1 , . . ., χ R ) is a vector of counting fields (see Appendix A for a detailed derivation).Based on the generalized time-evolution operator in Eq. ( 6), we define the photon-resolved time-evolution operators where m = (m 1 , . . ., m R ) is a vector of photonic transition numbers.In terms of the photon-resolved timeevolution operator, we can now express arbitrary system observables.For example, the probability that the photonic modes are in the Fock states n = (n 1 , . . ., n R ) is given as where Pn denotes the projector onto the Fock states with quantum numbers n, a m are the expansion coefficients of the photonic initial state m |m⟩ k , and ⟨ Ô⟩ t ≡ ⟨ψ(t) | Ô | ψ(t)⟩, where |ψ(t)⟩ can be either a state in the matter system or the composite light-matter system depending on the enclosed operator Ô.Thus, by semiclassically calculating the generalized time-evolution operators in Eq. ( 6), we can evaluate genuine quantum properties of photonic observables.

C. Full-counting statistics
While the analytical (or numerical) evaluation of the photon-resolved time-evolution operators in Eq. ( 8) is an infeasible task in many cases, it is relatively easy to compute the moments and cumulants of the photonic modes âk .Importantly, the counting statistics within the PRFT is fundamentally different from the standard FCS formalism.As shown in Fig. 1(a), the standard FCS framework is based on two-point projective measurements, where the state is formally projected to the Fock state basis at the beginning t 0 and at the end t 1 of each experimental run [69,93,103] (see Appendix B for more details).However, performing a projective photon number measurement at the beginning of the experiment would completely destroy the coherent photonic state.To circumvent this issue, we propose using two-point tomographic measurements, which are performed at the beginning and the end of the time evolution.The tomography is independently carried out for two batches of experimental runs with the same initial states.As illustrated in Fig. 1(b), the photon statistics is determined by photonnumber measurements at t = t 0 for the first batch, and at t = t 1 for the second batch.This alternative approach to FCS has been investigated for heat transport between incoherent baths in Ref. [104].
Formally, the counting statistics of the photon modes can be calculated via the cumulant-or the momentgenerating functions, which are defined by respectively.The associated nth cumulant and moment of mode k are determined via We are interested in the time evolution of the cumulantand moment-generating functions.To this end, we define the dynamical cumulant-generating function K dy,χ (t 1 ) via where K χ (t 0 ) and K χ (t 1 ) can be determined by independent tomographies at the beginning and the end of the time evolution [104].As shown in detail in Appendix A 3, the dynamical cumulant-generating function can be expressed as where φ = (φ 1 , . . .φ R ) is the vector of phases of the photonic states.Employing Eq. ( 13), we can now obtain the change of the cumulants, which describes the change of the photon statistics in the two-time tomographic measurement sketched in Fig. 1(b).We can see for instance that the first and second dynamical cumulants correspond to the change of the mean photon number ⟨ Nk ⟩ and the variance σ 2 k , respectively.In the rest of this paper, we will analyze Eq. ( 14) in a variety of contexts to gain a transparent picture of the photonic dynamics.The derivation of Eq. ( 14) assumes a separable initial state of the form in Eq. (2).In Appendix A 6, we explain how the PRFT can be generalized to more general initial states such as entangled states in a similar fashion.
It is worthwhile to point out that the well-defined phase in the PRFT is in contrast to the well-defined particle number in the standard FCS.This is the origin of the difference in the measurement protocols for the two cases.As shown in Appendix B the cumulant-generating function is Kχ (t) = log⟨U † φ−χ/2 (t)U φ+χ/2 (t)⟩ t0 according to the standard FCS, where the phases φ have been evaluated heuristically.This results is a profound difference with the predictions of the PRFT, and highlights the crucial role played by two-point tomographic measurements in the PRFT.We note that the first cumulant in Eq. ( 16) agrees with the standard semiclassical calculation as shown in Appendix A 7.

D. Probability redistribution
Using the moment-generating function of the final state given in Eq. (11), we can calculate the probabil-ity distribution of the Fock states at time t via where we have introduced the quasiprobabilities as where M dy,χ (t) = exp [K dy,χ (t)].Akin to regular probabilities, these quasiprobabilities are real valued and fulfill However, due to interference effects, q n may also be negative.Such negativity is a characteristic signature of nonclassical temporal correlations [105].We note that according to Eq. ( 17), the quasiprobabilities may be interpreted as the kernel of the probability redistribution.

E. Full-system state
The PRFT not only predicts the photonic probability distribution, but also the state of the matter system, from which we can even reconstruct the state of the total system.We can express Eq. ( 17) in the form where we have defined the probability operators It is easy to show that the set of probability operators fulfill Moreover, in the parameter regime and time scope in which the PRFT is valid, we have 0 ≤ Pn (t) ≤ 1, which is equivalent to 0 ≤ p n (t) ≤ 1 in Eq. (20).Thus, the operators Pn (t) define a positive-operator-valued measurement (POVM) [106], which consistently describes the photon measurement process sketched in Fig. 1.
According to the theory of quantum measurements [106,107], the reduced density matrix of the matter system after the measurement is given by where ρ(t 0 ) = |ϕ(t 0 )⟩ ⟨ϕ(t 0 )| is the initial density matrix, and ρ n denotes the reduced density matrix conditioned on the measured photon number n.As the POVM Pn (t) acts on the initial state, the time evolution operator U φ (t) has been added heuristically to account for the dynamics of the matter system.Equation ( 23) thus defines a valid set of Krauss operators describing photon-numberdependent quantum channels.If ρ(t 0 ) is a pure state, then ρ n (t) ≡ |u n ⟩ ⟨u n | will be also pure as Pn (t) is positive semidefinite.A purification of ρ M is given by where we have assigned the phases e −i(ωt−φ)•n in terms of the frequency vector ω = (ω 1 , . . ., ω R ), such that we can interpret the states |n⟩ as the Fock states of the photonic system.Tracing over the light system, we can verify that indeed tr L [|Ψ(t)⟩ ⟨Ψ(t)|] = ρ M (t).Thus, we can identify the |Ψ(t)⟩ as the state of the total light-matter system.

F. Floquet-state analysis
We now proceed to analyze systems subjected to multimode driving where the photonic field is composed of multiple commensurate photonic frequencies ω k .In this scenario, the matter subsystem is still described by a Floquet Hamiltonian in the semiclassical limit.However, we find that the PRFT leads to some important insights about such systems (like light-matter entanglement), that are beyond the reach of a standard Floquet analysis and imply a fundamental limitation on its validity.
According to Floquet theory, the generalized timeevolution operator can be written as where the generalized Floquet Hamiltonian can be expanded as As the kick operator U kick,χ (t) = U kick,χ (t + τ ) is time periodic, it accounts only for periodic changes in the photon redistribution.Thus, the asymptotic dynamics of the cumulant-generating function in Eq. ( 14) is determined by the generalized quasienergies E µ,χ .The stroboscopic dynamics of the matter system is characterized by the Floquet states |u µ,χ ⟩, that generalize the common eigenstates in time-independent systems.When expanding the initial state in the Floquet-state basis |u µ,φ ⟩ where |A(t 0 )⟩ is the initial state of the photonic field, the asymptotic dynamical cumulant-and moment-generating functions read as e i(Eµ,φ−Eµ,φ+χ)t + e i(Eµ,φ−χ−Eµ,φ)t .
This derivation is rigorously explained in Appendix A 4. Thus, M dy,χ is a weighted average of the dynamical moment-generation functions of the Floquet states, with the weights given by the expansion coefficients |c µ | 2 .We can now obtain the mean photon number change, where we have evaluated the quasienergies at the phases of the photonic fields φ, and we have denoted the meanphoton change in a specific Floquet state by ∆⟨ Nk ⟩ | µ .Intriguingly, when preparing the system initially in an arbitrary Floquet state µ, we find i.e., the variance change ∆σ 2 k vanishes.For a superposition of Floquet states, Eqs. ( 18) and (28) imply that the quasiprobabilities are given by a weighted average of Floquet-state dependent quasiprobabilities q n|µ , i.e., q n = µ |c µ | 2 q n|µ .Similarly, from Eq. ( 17), we can infer that Thus, p n|µ defines a conditional probability distribution, and the dynamics of the photonic state is controlled by the Floquet states.
Our results lead to a quantum-optical interpretation of Floquet states.To see this, we consider the timeevolution of the generic initial state in Eq. (27).According to Eq. ( 31), the light-matter state becomes entangled in the course of the time-evolution: where |u µ,φ (t)⟩ = U kick,φ (t) |u µ,φ ⟩, and the photonic wave functions for long times is given by which is in agreement with Eq. (24).Thereby, the mean photon number of each |A µ ⟩ changes linearly in time [see Eq. ( 29)], while its variance remains unchanged [see Eq. (30)].This means that the conditional probability can be approximated by where [∆n µ (t)] k = −dE µ,φ /dφ k t, which can be numerically efficiently evaluated.Consequently, the photonic states |A µ (t)⟩ will become mutually orthogonal for sufficiently long times, and the reduced density matrix of the matter system becomes, Thus, from a quantum-optical point of view, Floquet states act as the decohering basis.We note that this interpretation holds as long as there is a nonvanishing photon flux between distinct photonic modes.This effect is illustrated in Fig. 1(b) and will be further analyzed in Sec.III B for a two-mode Rabi model.We emphasize that a large-scale photon flux cannot develop when the quasienergy does not depend on the counting field, as we see from Eq. ( 28).This situation occurs for single-mode systems, where the single counting field can be transformed away by a counting-fielddependent shift in time t → χ/ω, which can be seen in Eq. (7).Consequently, the transport-induced lightmatter entanglement can strictly appear only in two-or higher-mode systems.We recall however that there are always two polarization modes of light, such that this decoherence effect discussed here has practical importance.
At this point, it is instructive to compare our results with Ref. [53], which has computed photonic observables by considering a static phase for photonic coherent states.However, the static-phase assumption leads to unphysical predictions such as the diverging higher-order moments and cumulants and, thus, is not suitable to quantitatively capture the dynamics of the photon field or the decoherence of the matter system.It is worth mentioning that the PRFT formally operates in the Sambe space rather than in Fock space.Thus, the PRFT cannot account for the shot-noise induced entanglement discussed in Sec.II A. However, as we demonstrate in Sec.III A, this effect has a minor influence on the photonic dynamics.

G. Error analysis
To describe the deviation from the exact time evolution quantitatively, we specify the initial state to be where |ϕ(t 0 )⟩ is the initial state of the matter system, and the photonic state of mode k is given by Thereby, n k is the mean photon number, σ k is the standard deviation, φ k is the mean phase, and N is a normalization factor.For σ k < √ n k (σ k > √ n k ) the system is in a number (phase) squeezed state, while for σ k = √ n k it is in a coherent state.As discussed in details in Appendix A 5, the deviation of the probabilities p n in Eq. ( 17) from the exact ones p (Ex) n scales as where F [x] denotes the scaling function, which depends on the set of ratios n k of all photonic modes k.This shows that the PRFT performs well for photonic states with a large standard deviation σ k and large mean photon number n k .Large σ k makes sure that the phase is well defined, which is reflected by the first argument of the function F in Eq. (38).A large n k guarantees that the matrix elements of the photonic operators âk do not depend on the photon number, which is described by the second and third argument of F.
Clearly, the latter source of error is absent when considering the Floquet theory in Sambe space [introduced in Eq. ( 4)].All requirements are naturally fulfilled for coherent states with Intriguingly, the PRFT makes accurate predictions, even when σ k is small; we discuss this in Sec.III B.
The linear change of the mean photon number in Eq. ( 29) is a consequence of the vanishing photon-number dependence of the matrix elements ⟨n This establishes a 'translational invariance' in Fock space.In each unit of time, the matter system can pump a certain number of photons from one driving mode into another independent of the photon number.When the mean photon number change reaches the order n k , the translational invariance in Fock space is lost, and the PRFT theory breaks down.

III. APPLICATIONS
We apply the PRFT to three versions of the quantum Rabi model, and demonstrate that the framework is accurate in the semiclassical limit.The Hamiltonian describing the system is where the two-level system is described by the common Pauli matrices σα with α = {x, y, z}, and h z denotes level splitting of the two-state model.We denote the eigenstate of σz with eigenvalue 1 (−1) by |↑⟩ (|↓⟩).The photonic system operators and parameters have been described in Sec.II A.  To compare the photonic probability distributions predicted by the numerical exact quantum calculation p and by the PRFT p n k of mode k, we employ the trace distance We quantify the entanglement of the light and matter systems in terms of the purity of the two-level system P ≡ tr ρ 2 M , where ρ M is the reduced density matrix of the two-level system.For a pure state P = 1, while for a completely mixed, i.e., maximally entangled, state P = 0.5.

A. Rabi model
First, we investigate the paradigmatic single-mode quantum Rabi model with R = 1.In the semiclassical limit, the atomic dynamics is described by H R = h z σz /2 + 2g 1 σx cos(ω 1 t).In the following discussions, we neglect the index k = 1.Since there is no possibility for large scale photon transport in the single-mode Rabi model, the error will be mainly determined by the photon shot-noise of the initial states.As we will see, this error is small as compared to the transport dynamics introduced in Sec.II F and discussed in Sec.III B. In Appendix C 1, we also benchmark the PRFT against the dynamics in the Sambe space [introduced in Eq. ( 4)], where the accuracy of the PRFT is even further improved.
Time evolution.In Fig. 2(a), we depict the dynamics for n = 10 5 photons and two different σ = √ n (red, dash-dotted) and σ = 0.1 √ n (yellow, solid), i.e., for a coherent and a number squeezed photonic state.The PRFT results are depicted by a black dashed line.We observe that the PRFT results of the spin observables ⟨σ α ⟩ with α ∈ {x, y, z} agree perfectly to the numerical results for the coherent and number-squeezed initial conditions.Likewise, we do not observe differences in the mean photon number change ∆ N .For the photon variance change ∆σ, we observe that the PRFT prediction agrees reasonable well with the exact calculation for the number squeezed state for short times.Yet, the PRFT strongly deviates for the coherent photonic initial state.These findings are in agreement with the error analysis in Eq. (38), which shows that the error scales with σ k / √ n k due to the photon shot noise.Yet, we note that deviations in the probability distributions are heavily enhanced by the definition of the variance, which disproportional weights photon numbers away from the mean value with ∝ (n − n) 2 .Intriguingly, the PRFT agrees perfectly with the numerical results when the dynamics is simulated in the Sambe space instead of the Fock space as demonstrated in Appendix C 1.For these reasons, we continue to analyze the photon statistics in terms of the trace distance defined in Eq. (40).Moreover, we mention that the observed variance values in Fig. 2 are small compared to the rapidly diverging variance, which can appear in multimode systems due to photon transport.
Photon probabilities.In Fig. 2(b), we benchmark the photon probability distribution in the weak-(upper row), intermediate-(middle row), and strong-coupling (bottom row) regimes for number-squeezed (left) and the coherent (right) initial states, respectively.To test the validity of the PRFT, we choose a scaled time t σ ≡ (2π/g)6σ, which agrees with the analytically predicted validity according to the first argument in Eq. (38).The scaled time t σ is a sufficient time scale to observe the transport entanglement effect explained in Sec.II F and demonstrated in Sec.III B. In the weak-and the intermediate-coupling regime, we observe that the PRFT agrees well with the exact numerical calculations.For the coherent photonic state, we observe some minor deviations exhibiting an oscillating dependence with photon number, which will be explained below.
Error scaling.In Fig. 2(c), we investigate the error quantified by Eq. ( 40) in the weak-, intermediate-and strong-coupling regimes at the scaled time t σ .Thereby, we investigate the error for three different level splittings h z .In the left column, we investigate the error as a function of σ for n = 10 5 .In the weak-and intermediatecoupling regimes we find that the error increases with σ in agreement with the second and third argument in Eq. (38).For very small σ ≈ 0.05 √ n, we also observe a rapid error increase, which is due to the first argument in Eq. (38).In contrast, the error exhibits an overall decaying behavior as a function of σ in the strong-coupling regime, suggesting that the first argument in Eq. ( 38) plays a more prominent role.
In the right column of Fig. 2(c) we investigate the error as a function of n(t 0 ) while simultaneously scaling σ = 0.3 n(t 0 ).In the weak-and intermediate-coupling regime, we find that the error scales approximately as d 1 ∝ n(t 0 ) −0.5 ∝ σ/n, i.e., according to the second argument in Eq. (38).While it might appear that the error scales as d 1 ∝ n(t 0 ) −1 ∝ 1/σ 2 in the strong-coupling regime, we interpret this as a consequence of numerical fluctuations, which are caused by the σ-dependent evolution time t σ .Similar observations can be also found for other ratios of σ/ n(t 0 ) and matter initial states (not shown).
Overall, we find that the error is smaller than d 1 ⪅ 0.01 for the Rabi model at resonance h z ≈ ω.Away from the resonance condition h z > ω, we find that the error even decreases.We explain this by a reduced shot-noise entanglement, which occurs most efficient at resonance.Analytical analysis.When |h z − ω| ≪ g, we can neglect the counter-rotating terms σ + â † , σ − â and the Hamiltonian reduces to the Jaynes-Cummings model In this case, the photonresolved time evolution operators in Eq. ( 8) are, where is the energy of the excited eigenstate (see Appendix C 2 for a detailed derivation).This form of the photon-resolved operators encodes the conservation of the quantity, N tot = â † â + σz , which is a salient feature in the Jaynes-Cummings model.To illustrate this, let us consider the initial state |ψ(t 0 )⟩ = |↑⟩ ⊗ |n⟩ , that gives rise to the following photonic occupations of the Fock states n and n + 1 which maintains the probability ⟨ Pn ⟩ + ⟨ Pn+1 ⟩ = 1.Intriguingly, the occupations in Eq. ( 42) are identical to the exact time evolution of the quantum Jaynes-Cummings model.Thus, we have determined the dynamics of a genuine quantum model by employing semiclassical methods of the PRFT.Analysis of Eq. ( 42) explains the minor oscillating deviations in the weak-and intermediate-coupling regimes in Fig. 2(b) for the coherent photonic state.The PRFT assumes a fix Rabi frequency g = g√ n.Taking a more microscopic perspective, each initial Fock state |n⟩ determines its specific oscillation frequency g(n) = g√ n, such that the probabilities in Fig. 2(b) oscillate slower (faster) for photon numbers n below (above) the mean photon number, leading to the observed minor derivations.This interpretation is further underpinned by the improved accuracy of the PRFT in Sambe space, which is analyzed in Appendix C 1.

B. Two-mode Rabi model
We found for the R = 1 Rabi model that the shotnoise entanglement causes a small error for the photonic probabilities (d 1 ⪅ 0.01 in Fig. 2), which is even vanishing in the thermodynamic limit n → ∞, σ → ∞ or when working in the Sambe space defined in Eq. ( 4).As a more advanced example, that still allows for numerical benchmark calculations, we consider now the two-mode quantum Rabi model in Eq. ( 39) with R = 2.This model allows for a large-scale photon transport between the two photonic modes, which leads to a more prominent photonic dynamics as in the R = 1 model, and a light-matter entanglement effect, that persists in the thermodynamic limit.
The initial state is given in Eq. (36) for n k = 5000 and three different photonic distribution widths σ k ∈ {2, 3, 4} for both modes k = 1, 2. Larger values for σ k cannot be numerically simulated.According to the error scaling in Eq. ( 38), the PRFT requires large values of σ k .However, we find that the PRFT is already very accurate for σ k = 4 and improves quickly.In Appendix C 1, we also benchmark the PRFT for various mean photon numbers n k , and in the dephasing, adiabatic, and high-frequency driving regimes.
As explained in Sec.II F, the time evolution sensitively depends on the initial state.For this reason, we carry out benchmarking for an initial Floquet state |ϕ(t 0 )⟩ = |u µ ⟩, and an initial spin-up state |ϕ(t 0 )⟩ = |↑⟩ in Figs. 3 and  4, respectively.In the first column of Figs. 3 and 4 we depict the variational distance d 1 defined in Eq. (40).The second column depicts the probability distribution at selected times.The third column shows the expectation value ⟨σ y ⟩ of the two-level system, while the forth column depicts the purity of the matter system.
For short simulation times t < 2π/ω, we use Eq. ( 17) to determine the photon probability distribution, while for longer times t > 2π/ω we employ Eq. ( 34).To calculate ⟨σ y ⟩ and P within the PRFT, we use Eq. ( 23) for short simulation times.For long times, we employ the reduced density matrix of the state in Eq. (32).For the spin observables, we also depict the predictions of the standard Floquet theory.In the following, we discuss the performance of the PRFT in the weak-, intermediate-, and strong-light-matter-coupling regimes.
Weak coupling.In the weak-coupling regime for an initial Floquet state depicted in Fig. 3(a), we observe that the photonic probability distribution is shifted to smaller photon numbers n 1 with increasing time.For small σ 1 = 2, the compact initial photon distribution significantly diffuses with time.However, when increasing to σ 1 = 4, we already find a very good agreement with the exact probabilities.The improvement with increasing σ 1 can be also clearly seen in the trace distance d 1 .This improvement with increasing σ is according to the first argument in Eq. ( 38), while the two other arguments have a negligible influence on the dynamics.An excellent agreement of both the standard Floquet theory and PRFT to the exact calculation can be also observed for ⟨σ y ⟩ and the purity.
As we investigate here the resonant system with h z = ω 1 = ω 2 in the weak-coupling regime, we can apply the rotating-wave approximation, and investigate the corresponding two-mode Jaynes-Cummings model.As explained in Sec.II F, the asymptotic dynamics is defined by the counting-field depend quasienergies, which in this case can be exactly calculated and read as with µ = 1, 2, where we have defined G(χ) = k=1,2 g k e −iχ k .The corresponding Floquet states are given by √ 2 |u µ,χ ⟩ = e iϕχ |↓⟩ ± e −iϕχ |↑⟩ with ϕ χ = arg G(χ).Interestingly, the quasienergies depend only on the difference χ = χ 1 − χ 2 .Evaluating the first dynamical cumulant, we obtain where φ = φ 1 − φ 2 is the phase difference of the two photonic states.For φ ̸ = {0, π}, there is a net photon flux between the photon modes.The variance change ∆σ 1 vanishes according to Eq. (30).The current vanishes for φ = {0, π} due to the symmetry of the initial condition.In Fig. 3(a), we consider Floquet state µ = 1, such that photons are transported from mode k = 1 to k = 2, as can be inferred from Eq. ( 44).Likewise, photons would be transported from k = 1 to k = 2 when the system is initialized in Floquet state µ = 2. Figure 4(a) depicts the same as Fig. 3(a) but for the initial state which is here a balanced superposition of the two Floquet states.The overall accuracy is similar for both initial conditions.According to Eq. ( 29), the photon flow is controlled by the Floquet state, such that the photon redistribution in the long-time limit becomes entangled with the two-level system.This effect can be clearly observed for σ 1 = 4, where the left peak is entangled with the Floquet state |u 1,φ ⟩, while the right peak is entangled with |u 2,φ ⟩ [see also Fig. 1(c)].
The light-matter entanglement leads to decoherence as is clearly manifested in the time evolution of ⟨σ y ⟩.It is noteworthy that the PRFT agrees almost perfectly with the exact quantum calculation, while the standard Floquet theory strongly deviates.In particular, the PRFT accurately predicts the evolution of an initial pure state (with purity P = 1) to a maximally entangled state (purity P = 0.5).In contrast, the purity remains 1 in standard Floquet theory.
Intermediate coupling.In Fig. 3(b) we depict the dynamics for an initial Floquet state in the intermediate coupling regime for h z = ω = g k = g.As in the weakcoupling regime, we observe that the PRFT calculation rapidly approaches the exact probabilities for increasing σ 1 .The mean photon number n 1 (t) increases linearly, while the width σ 1 (t) stays almost constant, which can be seen for σ 1 = 4.We note that even though the difference between exact numerics and the PRFT is larger than in the weak coupling case, the convergence can be clearly anticipated even for small σ 1 .In general, large couplings g result in more Fourier components in the periodic dynamics of the Floquet states, that require an initial probability distribution with larger σ 1 to be smoothed out (see also Appendix A 5).
Furthermore, similar to the dynamics in Fig. 4(a), we observe that the probability distribution eventually splits into two peaks in Fig. 4(b), each being entangled with a Floquet state.The height of each peak is thereby deter- 20 10 0 10 20 20 10 0 10 20 n 1 ) as a function of time.The second column shows the photonic probability distributions for selected times t1g = 4, t2g = 8, and t3g = 20 as labeled in the panels.The third column depicts the expectation value of ⟨σy⟩ as a function of time.The fourth column shows the purity.(a), (b), and (c) depict the dynamics in the weak-, intermediate-, and strong-light-matter-coupling regime, respectively.Colored lines depict the numeric quantum calculations.We choose a symmetric coupling of both modes to the two-level system g k = g and n k (t0) = 5000 for k = 1, 2. Other parameters are specified in the figure .mined by the amplitude of the expansion coefficients in the Floquet basis according to Eq. (31).
Strong coupling.In Figs.3(c) and 4(c) we analyze the dynamics in the strong-coupling regime.As in the weak-and intermediate-coupling regime, we observe that the PRFT calculations for both the photonic and spin observables approach the exact calculation for increasing σ 1 .In contrast, the standard Floquet theory clearly fails 20 10 0 10 20 20 10 0 10 20 to reproduce the correct dynamics of the two-level system.
Due to numerical limitations, the probability distribution is depicted only for times t smaller than the driving period t i ≪ τ = 2π/ω in Fig. 3(c).As the findings in Sec.II F are only valid at stroboscopic times t = nτ , we do not observe a steady photon flux towards a higher (lower) photon number as in Fig. 3 (a) [Fig.3(b)].However, we expect that both peaks will merge at stroboscopic times t = nτ at a mean photon number as pre-dicted by Eq. (29).Similarly, we cannot associate the two peaks with the two Floquet states in Fig. 4(c) for the depicted times t j ≪ τ .

C. Three-mode Rabi model
To illustrate that the light-matter entanglement explained in Sec.II F is a generic effect, we now investigate a three-mode Rabi model with distinct commensurate photonic frequencies ω k .For concreteness, we choose ω 2 = 2ω 1 and ω 3 = 3ω 1 .In this case, the semiclassical Hamiltonian is a periodically-driven Rabi model with period 2π/ω 1 , and we compute the photon-mean and photon-variance changes using the PRFT.Noteworthy, a simulation of the full quantum Rabi model in Eq. ( 39) is numerically hardly tractable due to the three photonic modes.When representing each mode with m states, the Hilbert space has dimension D = 2×m 3 , i.e., D = 2×10 6 states for a moderate m = 100.Similar to the analysis for the two-mode Rabi model, we consider the three initial conditions |ϕ(t 0 )⟩ = |u 1,φ ⟩, |ϕ(t 0 )⟩ = |u 2,φ ⟩, and , where the |u µ ⟩ are the Floquet states of the two-level system, i.e., the eigenstates of Eq. ( 26) for χ = 0.
In Fig. 5, we depict the photonic mean and variance change of mode k = 3 for each initial condition.When the initial state is the Floquet state |u 1 ⟩ (|u 2 ⟩), the mean photon number grows (decreases) linearly in time, while the variance remains almost unchanged.This is analog to the probability distributions in Fig. 3 (a) in the twomode model, where the variance remains constant for a Floquet initial state.
In contrast, the mean photon number change is close to zero for the superposition state, while the variance increases rapidly.This corresponds to the probability distributions for an initial superposition in Fig. 4(a).Since we have chosen a balanced superposition of the two Floquet states with c 1 = c 2 = 1/ √ 2, the mean photon num-ber remains unchanged according to Eq. ( 29).However, as the Floquet state µ = 1 linearly increases the photon number, while Floquet state µ = 2 simultaneously linearly decreases the photon number, the variance diverges quadratically.

IV. EXPERIMENTAL AND TECHNOLOGICAL IMPLICATIONS
In this section, we discuss the light-matter induced decoherence for coherently driven systems.This effect can be detected in experiments and has a significant impact on the design of quantum memories and quantum operations.The following analysis therefore focuses on the transport-induced entanglement and neglects the shotnoise entanglement.

A. Quantum-optical coherence time
To estimate the quantum-optical coherence time for periodically driven systems, we consider the mean photon-number changes ∆⟨ Nk ⟩ µ = −E ′ µ (φ k )t for two distinct Floquet states µ 1 and µ 2 , where E ′ µ (φ k ) = dE µ (φ)/dφ k .The system is completely decohered at time t c when the difference in the mean photon number for these two Floquet states exceeds the width of the photon distribution in Fock space, i.e., where σ 2 k is the initial variance of photon mode k.We recall that the variance does not change with time for Floquet states as shown in Sec.II F. To estimate σ k in terms of physical quantities, we distinguish two cases: (i) a closed light-matter system, where the photonic field is confined in a cavity; (ii) an externally driven quantum system, where the photonic field is a traveling wave.Each case gives a different scaling behavior for t c .
Closed light-matter systems: For simplicity, we consider coherent states, whose mean and variance are equal The mean photon number in a cavity mode is given by n k = ϵ 0 E 2 V /(2ℏω k ) where E is the electromagnetic field, V is the cavity volume, and ϵ 0 is the dielectric constant.This implies that the quantum-optical coherence time is given as We note that the coherence time in Eq. ( 46) is an approximation, since the above arguments have assumed that the photon modes decohere independently.
Externally-driven quantum systems: For the experimentally more relevant situation, in which the matter system is driven by a traveling wave, the estimate for the coherence time has to be modified.Here we establish a connection to a cavity setup to get an estimate for the coherence time.Consider that the light field is a pulse of duration t p with central frequency ω k and spectral width ∆ω ∝ 1/t p .For long times t p , the spectral width vanishes and we assign all pulse photons to the central frequency.The mean photon number in a pulse is n k = P (ω k )t p , where P (ω k ) is the power of the electromagnetic field at frequency ω k .Now, we model the pulse as a cavity mode with initial occupation n k .Consequently, we find that the quantum-optical coherence time is, We note that since each photonic mode contributes to the decoherence effect, the coherence time in Eq. ( 47) can be considered as an upper bound.

B. Experimental verification
The decoherence effect discussed above provides a route to verify the PRFT without measuring photon statistics, which may be a challenging task.We explain our approach for the two-mode Rabi model for illustration.The key idea is to determine the decoherence effect by purity measurements of the atomic system.For this, it is crucial to isolate the quantum-optical decoherence from other decoherence sources.We achieve this by varying the amplitudes of the coherent states α k and the lightmatter interactions gk such that the effective parameters g k = gk α k remain constant.In doing so, the spin system experience the same semiclassical driving fields [see Eq. (C5)] and, thus, is subject to the same decoherence sources other than the quantum-optical decoherence.As for a coherent driving field σ k ∝ α k ∝ P (ω k ) 1/2 , the quantum-optical decoherence dynamics can be accessed by measuring the purity decay for various driving-field powers P (ω k ).
To get an estimate for the coherence time in the optical regime, we evaluate Eq. ( 47) for the two-mode Rabi model, where the photonic frequency is assumed to be ω = 400 THz.The quasienergy difference in the denominator in Eq. ( 47) can be approximated by a typical Rabi frequency of atomic systems Ω R = 40 MHz, generated by a laser with power P (ω) = 10 µW.For these parameters, the quantum-optical coherence time is t c = 5 ms, which is comparable to the duration of typical cold-atom experiments [108][109][110], but significantly shorter than quantum information storage times achieved with trapped ions [111].
We obtain the order of magnitude of the coherence time for radio frequencies using the following parameters: ω = 10 MHz, Ω R = 10 kHz, and P (ω) = 10 W, which are typical in current experiments [112].In this regime, the quantum-optical coherence time is very long (t c ≈ 3×10 18 s).This is a consequence of the high photon occupation of radio-frequency modes for realistic experimental parameters.Thus, in the radio-frequency regime, the electromagnetic field can be considered as completely classical, and the decoherence can be neglected.

C. Quantum memories and quantum operations
The quantum-optical decoherence in driven systems can have a significant impact on the design of quantum memories and quantum operations.To achieve long quantum information storage times, sophisticated control protocols have been developed, that typically involve time-periodic electromagnetic fields.Using dynamical decoupling, quantum information could be conserved for more than six hours in rare-earth atoms embedded in a crystal structure [112], and more than 50 s in trapped ions [111].In both cases, the quantum information is stored in the hyper-fine levels of the ground-state manifold, which are energetically separated in the radiofrequency regime.
Depending on the control protocol, the quantumoptical coherence time in Eq. ( 46) influences the performance of quantum memories.The quantitative considerations in Sec.IV B suggest that driving protocols involving optical-frequencies should be avoided, while radio-frequency control protocols are optimal.We note that quantum information storage and retrieval protocols often employ optical frequencies [113].Even though the pulse duration in these cases may be rather short, t p < 1 µs, an inappropriately adjusted pulse sequence might lead to a degradation of the quantum information via the quantum-optical decoherence effect.Similar considerations also apply to other quantum operations based on periodic driving, such as two-qubit gates.As quantum error correction usually requires high fidelities > 99% [114], even a minor quantum-optical decoherence can have a significant effect on quantum operations.Furthermore, inspection of the coherence time in Eq. (47) reveals that the coherence time can be enhanced by choosing control protocols for which the difference of quasienergies are not sensitive to the driving phases.
As the coherence time depends on the specific system, no general statements can be made about how to eliminate the quantum-optical decoherence effect.An alternate intriguing approach to mitigate this decoherence would be to employ quantum time crystals [19,20,24].Discrete time crystals, that are driven by an external driving field with period τ = 2π/ω, exhibit subharmonic response with a frequency ω/n, where n > 1 is an integer.This intriguing subharmonic response arises from the structure of the eigenspectrum which is composed of Floquet eigenstates that are separated by a quasienergy of ω/n.Since time crystals are robust to generic perturbations, the difference of quasienergies would have little dependence on the driving phases, leading to stable quantum memories.

V. QUANTUM COMMUNICATION
The light-matter entanglement discussed in the previous section can be employed in a quantum communication protocol that is robust against photon loss.To this end, we consider that Alice and Bob-the two participants in the communication-successively carry out the light-matter entanglement process described for the two-mode Rabi model and postselect the measurement results.We now proceed to explain how Greenberger-Horne-Zeilinger (GHZ) states can be employed to speedup the light-matter entanglement generation, before delving into the details of the communication protocol in Sec.V B.

A. Rapid generation of light-matter entanglement
Typically, the light-matter interaction between a single atom and the light-field is relatively weak, which slows down the generation of maximally-entangled light-matter states as depicted in Fig. 1(b).To enhance this effect, we again employ the setup in Fig. 1(b), but with N A atoms.The corresponding Hamiltonian reads as where the Pauli operators σ(j) z , σ(j) + , σ(j) − act on the atoms j = 1, . . ., N A .This is a rotating-wave approximation version of the Hamiltonian in Eq. ( 39) with many twolevel systems, i.e., the Tavis-Cummings model, in which the many-body interaction is mediated via the quantized electromagnetic field.As in Sec.III B, we consider two photonic modes k = 1, 2 which are initially in coherent states α k e iφ k .
The PRFT is a powerful tool to analyze this system when the atom numbers are large and exact numerical calculations are very expensive.Upon introducing the counting fields χ 1 and χ 2 , the corresponding semiclassical Hamiltonian reads as where H (j) (t) denotes the semiclassical Hamiltonian of atom j.The atoms are formally decoupled in the semiclassical description, yet, the interaction is still implicitly encoded in the counting fields.The quasienergy of the total system can be written as the sum of the quasienergies of the individual atoms where E µj ,χ is the quasienergy of atom j in Floquet state µ j .We have introduced the vector notation µ = (µ 1 , . . ., µ N A ), that contain the quantum numbers µ j of the N A atoms.The corresponding Floquet states read as where u (j) µj ,χ denotes the Floquet state of atom j with quantum number µ j .
To drastically enhance the number of photons transported from mode k = 1 to k = 2 (or vice versa), we prepare the system in either of the Floquet states characterized by the N A -component quantum numbers µ (1) = (1, . . ., 1) or µ (2) = (2, . . ., 2), i.e., the state in which all atoms j are in the same Floquet state µ j = α ∈ {1, 2}.The quasienergy of the atom ensemble in either of these states is then According to Sec.II F, this leads to an enhancement of the photon transport proportional to the atom number.
To be more specific, we choose h z = ω and the driving phases φ 1 = π/4, φ 2 = −π/4.In this case, the Floquet states of the matter system |u µ ⟩ for each atom are , and the Floquet states of the atom ensemble in Eq. ( 51) By constraining the system to the collective states |− − −⟩ and |+ + +⟩, we restrict the model to an effective two-level system with a renormalized quasienergy, for which the findings in the Jaynes-Cummings model in Appendix C 3 are valid.To create light-matter entanglement in the basis |+ + +⟩ and |− − −⟩, we use the initial condition where |GHZ⟩ is the celebrated GHZ state.This state is defined by a superposition |GHZ⟩ = (|0 0 0⟩ + |1 1 1⟩) / √ 2, where |0 0 0⟩ ≡ |0 . . .0⟩ and |1 1 1⟩ ≡ |1 . . .1⟩.The Hadamard gate ÛH = exp [−iπσ y /4] that locally rotates the state of each atom, is independently applied to all atoms.For later purpose, we illustrate the light-matter entanglement generation in Fig. 6(a) showing a superposition of states, where either driving field is enhanced and the other is reduced.We emphasize that even though we take the Tavis-Cummings model as an example, this enhancement effect is valid for general Floquet systems according to the PRFT.

B. Remote entanglement generation
The most crucial task of the quantum communication protocol is the generation of remote entanglement between two atomic ensembles possessed by Alice and Bob.The dashed lines sketch the light paths.After Bob's lightmatter entanglement, the coherent fields are measured (gray half circles).A vanishing signal difference heralds success of the protocol.The three insets show the photon distributions pn before and after Alice's entanglement process, and after Bob's process.(c) Analysis of the which-path information.After interaction with Bob's GHZ state, the system is in a superposition of four states.In two of which, the two output fields have changed their amplitude which reveals the path of quantum information and results in failing of the protocol.In the other two states, the which-path information remains hidden and results in success of the protocol.
Quantum state transfer can then be carried out via quantum teleportation [107].As schematically sketched in Fig. 6(b), this is achieved by repeating the light-matter entanglement process using the same driving fields.The three steps of the protocol are as follows: 1. State preparation: For an efficient light-matter entanglement generation, we assume that both Alice and Bob have prepared a GHZ state such that the initial state Alice and Bob carry out local Hadamard gates such that their states become where X = A,B.The choice of the basis |− − −⟩ and |+ + +⟩ is thereby determined by the driving phases, which we assume to be φ 1 = −π/4 and φ 2 = π/4.

2.
Light-matter interaction: Alice impinges two coherent light beams onto her atom ensemble leading to the generation of a light-matter entangled state according to the explanations in Secs.II F and III B: where |A + ⟩ (|A − ⟩) denotes a photonic state with enhanced (diminished) amplitude.The output light fields are transmitted to Bob, where they interact with Bob's atom ensemble.The resulting state can be written as in the basis of collective excitations |0 0 0⟩ X and |1 1 1⟩ X .
The projective measurement can be implemented by measuring the intensity difference of both output fields.If the difference is close to zero, then the two atoms have been successfully entangled.Otherwise, the process has failed and must be repeated.
The working principle of the protocol is based on the which-path information, that is illustrated in Fig. 6(c).After interacting with Alice ensemble, either of the output fields will be enhanced.This effect can be either repeated or reversed after interaction with Bob's ensemble.If the field k = 1 is two times enhanced, the intensity measurement reveals that it is in the photonic state |A 2+ ⟩.Consequently, Alice's and Bob's ensemble both are in the state |+ + +⟩ and thus not entangled.A similar reasoning applies to the state |A 2− ⟩.If the effect is reversed and the intensities of the output modes are equal, the which-path information of the photonic fields remains hidden, such that both ensembles preserve their uncertainty and become entangled.From Fig. 6(c) we find that the success probability is 50%.We note that the photon probability distribution of the coherent states decay exponentially from the mean value such that the states |A 2+ ⟩, |A 0 ⟩, and |A 2+ ⟩ can be distinguish with only a small error probability in the intensity difference measurement.

C. Quantum state transfer rate
The major obstacle in quantum communication is photon loss.A typical damping rate of optical fibers is γ ≈ 0.051 1/km, leading to a loss of more than 99% of photons after 100 km.This heavily limits the reach of quantum state transfer protocols based on few photons.Current theoretical transfer protocols predict transmission rates of up to 1 Hz over 500 km [102].These protocols typically employ quantum repeaters strategically placed between the transmission endpoints.
The coherent-light protocol introduced in Sec.V B is naturally robust against photon loss.The information of Alice's qubit is encoded as an enhanced or diminished light amplitude during transmission [see inset in Fig. 6(b)].When a photon is lost, it is hardly possible to determine from which transmission peak it originated.Still, while the which-path information is preserved, photon loss has a detrimental effect, as it leads to a broadening of the photonic probability distribution.When the broadening exceeds the distance of the two peaks, i.e., ⟨∆ N1 ⟩ ≈ ∆σ 1 the quantum information is lost.
Along the same arguments for photon loss, the quantum information is also robust against classical amplification.To compensate the photon loss, we assume that it will be amplified with rate γ Amp = γ, such the mean photon number is conserved during transmission.Yet, the amplification will lead to a broadening of the probability distributions.Modeling photon loss and amplification as independent Poissonian processes with rate γ, the width of each peak increases by ∆σ 1 = 2γdP t p /(ℏω), where P is the power of the transmitted pulse, t p is the pulse duration, and d the transfer distance.We recall from Sec. IV A that the separation between the two peaks is given as In atomic systems, the quasienergy splitting can be associated with the Rabi frequency Ω R .Using the GHZ amplification in Eq. ( 52), the probability peaks are separated by The quantum state transfer rate can be thus estimated as where we introduced S ≡ σ/ √ n as the ratio of the photonic distribution widths of a number-squeezed state σ and a coherent state √ n.For instance, we assume a photonic frequency ω = 400 THz and a typical Rabi frequency of atomic systems Ω R = 40 MHz, corresponding to a laser power P = 10 µW.Recent experiments of the Lukin group have successfully created GHZ states with N A = 12 Rydberg atoms [115].When assuming a coherent state with S = 1, the transfer rate is 122 Hz over 500 km, thus exceeding typical few-photon protocol by two orders of magnitude.

D. Discussion
Even though strongly idealized, the proposed protocol merits a thorough discussion.As quantum information is often stored in the ground-state manifold of atoms, that are coupled via Raman transitions, a more realistic modeling in terms of a three-or-more level system is required.We emphasize that the proposed protocol is not restricted to the Tavis-Cummings model used here for illustration, as the light-matter entanglement is a generic effect appearing in all Floquet systems as predicted by the PRFT.
There is a series of points that can be discussed independent of the concrete physical implementation.Along with future progress in quantum control, the suggested protocol has enormous development potential : • Atom number.The quantum state transfer rate in Eq. ( 57) scales quadratically in the number of atoms N A .It has been estimated that for faulttolerant quantum computation N A = 1000 physical qubits in a highly entangled state are required.This would increase the transfer rate by a factor of 10 4 compared to the estimate in Sec.V C.
• Transmission power.The transfer rate is inversely proportional to the laser power P .When Alice deamplifies the beam after interaction with her GHZ state prior to transmission by a factor of α, the transfer rate can be increased by that factor α. Bob must amplify the received signal before interaction with his GHZ state.Noiseless deamplification and amplification can be implemented with additional atom ensembles which are prepared in Floquet states as discussed in Sec.II F and shown in Figs.1(b) and 3(a).This might enhance the transfer rate by an additional factor of α = 10.
• Number squeezing.According to Eq. ( 57), the transfer frequency scales linearly with the number squeezing parameter S.Moreover, a moderate S also minimizes the shot-noise entanglement discussed in Sec.III A. We consider a number squeezing of S = 0.1 for the following estimation.
Taking these points into account, the transfer rate will be on the order of f ≈ 100MHz, and thus commercially relevant.Compared to few-photon quantum protocols, the coherent-state protocol comes with a serious of advantages: • Distance dependence.The transfer rate scales inverse proportional with the distance d.This is a more favorable scaling than for few-photon protocols whose quantum state transfer rate typically decrease exponentially due to intensive postselection in quantum repeaters [102].
• Simple implementation.The implementation is based on the light-matter entangling process.It does not require sophisticated encoding and decoding schemes to protect the quantum information that are experimentally and computationally challenging.
• On-demand light fields.Presently, the absence of on-demand single photon sources pose a major challenge to several quantum communication protocols.This problem is circumvented in the proposed protocol, as coherent pulses of light can be easily produced and controlled.
• Photon detection.The efficiency of single-photon detectors strongly influences the transfer rate of few-photon protocols.As the success flag in the coherent-light protocol is determined by measuring the intensity difference [c.f., Fig. 6(b)], a standard photon multiplier will be sufficient.
Compared with few-photon quantum communication protocols, the coherent-state entanglement protocol has two major disadvantages, though these can be overcome: • GHZ generation.To establish an efficient protocol, Alice and Bob must generate high-fidelity GHZ states, that might be technically challenging for large atom numbers.However, quantum information will be always stored in an encoded form.
While the GHZ state is the basis of Shor's quantum repetition code [116], other quantum error correction codes are based on graph states, that are generalizations of the GHZ states [114].Thus, the generation of GHZ states will only lead to negligible overhead.
• Noise and decoherence.As the quantum information is encoded in the photon number, the transmitted photonic state will be sensitive to phase noise and decoherence related to the operators â † k âk for k = 1, 2. Fortunately, as the ratio of ∆⟨ Nk ⟩ and ⟨ Nk ⟩ is on the order 10 −5 , the environment can learn only little about the quantum state, such that the which-path information is not leaked.To further enhance the protection of the quantum information, dynamical decoupling can be employed by periodically switching |− − −⟩ X ↔ |+ + +⟩ X for X = A, B, leading to a periodic variation of ∆⟨ Nk ⟩.
The entanglement generation based on coherent light can be interpreted as a physical encoding of the quantum information, which is in contrast to the logical encoding of quantum information typically deployed in quantum communication protocols [114].Assuming a pulse length of t p = 20 ms, the two peaks are separated by ∆⟨ N1 ⟩ ≈ 10 7 photons which corresponds to a Hamming distance of log 2 (∆⟨ N1 ⟩) ≈ 24 bits.As the total number of photons is on the order of ⟨ N1 ⟩ = P t/ℏω ≈ 10 12 corresponding to 40 qubits, the physical encoding presented here is thus comparable to a [ [40,1,24]] quantum error correction code.

A. Summary
In this subsection, we compare the PRFT with other well-established methods for analyzing light-matter systems.The PRFT combines important features of established frameworks, while avoiding their shortcomings: Floquet theory.The PRFT introduces counting fields into the semiclassical equation of the light-matter system to track the quantum dynamics of the photonic driving field, thereby making an important advancement to the framework of Floquet theory.Crucially, the PRFT defines a quantum channel for the dynamics of the driven matter system, that describes the decoherence induced by the light-matter interaction in the Floquet basis.This inherently quantum effect is completely neglected in the standard Floquet theory, which treats the matter subsystem as an effectively closed quantum system.We note that while this eventual decoherence may be anticipated from other semiclassical techniques [53], a quantitative calculation is generally beyond the reach of these methods.Our investigations clearly demonstrate that Floquet theory suffers from fundamental limitations in describing light-matter-coupled systems.The PRFT provides a semiclassical approach to address these issues, and emphasizes the need to carefully investigate the standard Floquet theory even in parameter regimes, in which it has been thought to be valid.The PRFT renders the Floquet theory as an open quantum system framework, by providing a microscopic derivation of the Kraus operators.
We emphasize that the PRFT has the same computational complexity as the standard Floquet theory, since it requires the integration of semiclassical equations.As a consequence, the PRFT has significant computational advantages over Sambe space methods that investigate photonic dynamics by effectively requantizing the semiclassical driving field [94][95][96].Finally, we note that, unlike other approaches, the PRFT can distinguish between modes with commensurate frequencies, thereby extending its reach over a wide class of driven systems.The phase representation approach in Ref. [53] derives photonic observables by formally considering the photon phase as a dynamical variable.However, this approach becomes problematic when specifying to coherent photonic states, for which the phase has been considered as static leading to unphysical prediction such as diverging higher-order moments and cumulants.For this reason, the method in Ref. [53] comes short to describe the light-matter entanglement, which takes a prominent role within the PRFT.
Quantum optics.An extremely appealing feature of the PRFT is that it makes accurate predictions about photonic observables, even though it only relies on semiclassical equations.Our framework thus provides a drastic computational advantage over established methods in quantum optics.In particular, we note that while it is extremely difficult to numerically investigate systems with more than two photon modes using traditional methods such as phase-space frameworks [55,117,118], the numerical effort is independent of the mode number in the PRFT.Intriguingly, the PRFT predicts the generation of light-matter entangled states in general multimode driven systems, which is essentially infeasible to describe with standard quantum optical methods.This entanglement effect is a fundamental consequence of the photon transport between different modes, which is controlled by the matter quantum system.It is thus distinct from the entanglement resulting from the photon-number uncertainty (i.e., the shot-noise) of a single mode in a coherent state [97][98][99][100][101].Our analysis has shown that the transport induced entanglement has a profound impact on the photon dynamics in contrast to the more wellstudied shot-noise induced entanglement.Furthermore, the PRFT is nonperturbative in nature, and thus capable of making predictions beyond methods from nonlinear spectroscopy [119].
Before proceeding further, we note that there are several straightforward applications of the PRFT that could not be addressed in this work.For instance, in the benchmarking of the quantum Rabi model, it has been assumed that the photon field is switched on in a quantum quench.A more realistic situation is a smooth switching of the photon density, that could be modeled as a coherent superposition in continuum mode photon field, or by a time-dependent light-matter coupling g(t).In all of these cases, the overall framework remains valid.
Full-counting statistics.The PRFT provides a stochastic description of photonic fields in terms of probabilities, moments and cumulants in a manner analogous to the standard FCS [69,120].The standard FCS is based on two-point projective measurements, and can describe spontaneous photon emission [66][67][68].However, the coherences in the photon number basis are formally destroyed in this method, thereby leading to wrong predictions for photon modes in coherent states.In contrast, the photon counting statistics in the PRFT is obtained by two-point tomographic measurements, that are compatible with coherent states.The formalism provides an exact expression for the change of the dynamical moment-and cumulant-generating functions, that describe the change in the photonic statistics of the photon modes.Finally, the quasiprobabilities capture the redistribution of the initial probability distribution and therefore describe the photonic dynamics.

B. Discussion and outlook
The PRFT has far-reaching implications for quantum science and technologies.In particular, the quantum optical decoherence predicted by the PRFT has serious consequences for quantum memories and quantum operations.
Quantum memories often use sophisticated Floquet control protocols, where the control fields can unintentionally induce quantum-optical decoherence.In the optical frequency regime, the quantum-optical coherence time predicted by the PRFT is reasonably short (≈ ms), which is orders of magnitude smaller than the targeted storage time of quantum information.This analysis suggests that optical frequencies should be avoided in quantum memories and quantum operations, as even a small decoherence is detrimental to maintain the high fidelity required for quantum error correction.Moreover, we have argued that quantum time crystals are ideal candidates for quantum memories.Finally, we have employed the PRFT to propose a quantum communication protocol that is robust against photon loss.Our proposal employs coherent photonic states instead of single photons to establish remote entanglement.The robustness of this protocol originates from the fact that the width of the photonic probability distribution increases only with the square root of the number of lost photons.Consequently, the quantum state transfer rate scales inversely with distance, thereby outperforming few-photon protocols based on quantum repeaters, that typically decay exponentially with distance.We plan to investigate detailed implementations of this protocol in future work.
The PRFT developed in this paper can potentially have a significant impact on various research directions.Some promising applications include thermodynamics [121], heat engines [88,122,123] and quantum phase transitions in interacting spin chains in a cavity [124][125][126], and control of many-body localization [127,128] in the presence of external driving.Thereby, the PRFT has a significant computational advantage over quantum optical methods, where accurate numerical calculations are extremely challenging in this regime.An extension to open quantum systems can also clarify the compatibility with the standard FCS.We further speculate that a suitable development of the PRFT will have important implications in spectroscopy and metrology.The PRFT can be applied in the analysis of highly occupied Fock-state lattices, which have been shown to exhibit an intriguing collapse and revival dynamics [129].Methods developed for noise suppression in electron transport can be also combined with the PRFT to control the photonic counting statistics [130].
where α ≫ 1 and φ are the amplitude and the phase of the photon field.We denote the resulting operator as U χ (t), which acts on the matter system.The transition to the PRFT is readily done by realizing that where T is the time-ordering operator and the generalized time-periodic Hamiltonian on the right-hand side is defined as The calculation of U χ (t) can be thus performed by analytically or numerically solving the time-dependent Schrödinger equation with H χ (t) for all χ ∈ [0, 2π).The final step in the transition is done by evaluating and replacing in Eq. (A9).This replacement is well justified as the matrix elements of the terms such as â † n âm depend only weakly on the photon number for large n 1 , n 2 .The expectation value of Pn in the PRFT thus reads as Pn = m1,m2 where the expectation value is taken in the matter initial state.

Full-counting statistics
In the following, we derive the dynamical cumulantgenerating function of the photon field given in Eq. ( 14).As the physical background and interpretation has been already explained in Sec.II, we focus here on the merely technical details.For simplicity, we count only the photons in mode k = 1.Nevertheless, we still implicitly allow for other photon modes, which are not explicitly counted.In this case, the cumulant-and moment-generating functions are given as where the time evolution is calculated with the full quantum Hamiltonian in Eq. ( 3).The moment-generating function in Eq. (A17) can be rewritten as where U χ (t) is the generalized time-evolution operator in Eq. (A3).Here, the moment-generating function is presented in a symmetric way.At this stage, an unsymmetrical representation would be also correct, however, this will get problematic when taking the semiclassical limit later.In doing so, we make sure that the essential transformation property M * χ (t) = M −χ (t) is maintained in the semiclassical limit, which guarantees that all moments and cumulants are real valued.
We assume that the initial state is separable where |ϕ(t 0 )⟩ is the initial state of the matter system.We can expand the initial state of the light field in the Fock basis Using now the photon-resolved time-evolution operator in Eq. (A5), we can expand Eq. (A18) in terms of photon processes such that In this expression, we have already carried out the semiclassical replacement of the photon-resolved timeevolution operators introduced in Eq. (A15).The validity of this replacement will be analyzed in Appendix A 5.
To make progress, we apply a Fourier transform to the photonic expansion coefficients Moreover, using the representation of the photon-resolved time-evolution operators in Eq. (A14), the moment-generating function reads as Albeit general, Eq. ( A23) is inconvenient to evaluate both analytically and numerically.Motivated by the abundant use of lasers and other coherent electromagnetic fields in experiments, we focus in the following on Gaussian photonic states, for which the expansion coefficients are given by where n is considered as a continuous variable, n ≫ 1 and σ denote the mean photon number and width, respectively, and φ is the phase.For σ 2 = n, the state is a coherent state, while for σ 2 < n (σ 2 > n) it is denoted as a number-squeezed (phase-squeezed) state.More general photonic states, such as multimode squeezed states and light-matter entangled states will be discussed in Appendix A 6. Expressed as a function of phase, the coefficients in Eq. (A22) read as e −(φ1−φ)σ 2 (φ1−φ) e i(φ1−φ)n .
For the first term, evaluation of the Gaussian integral gives where in the last equality we have identified the moment-generating function at time t 0 , i.e., M χ (t 0 ) = n a * n a n e −inχ .Evaluating the Gaussian integrals of the second and third terms in Eq. (A27), we obtain and

Error analysis
Here we analyze the error of the semiclassical momentgenerating function in Eq. (A31).To allow for a quantitative statement about the error, we consider the generic photonic state in Eq. (36).To simplify the notation, we consider here a single photon mode, while the generalization to a multimode system works along the same lines.
Informally speaking, the derivations in the PRFT make the following assumptions for the expansion coefficients a To quantify the error, we have to investigate two approximations: (i) the semiclassical replacement of the photon-resolved time-evolution operators in Eq. (A21); (ii) the higher-order contributions M (l≥1) χ in Eq. (A27).
(i) Semiclassical replacement.In the semiclassical replacement in Eq. (A21), we have assumed that the matrix elements of the photonic operator a are independent of the photon number during the time evolution.To estimate the validity of this approximation, we consider the ratio of matrix elements for the photon numbers n = N and n = N + κ dy,1 t + σ, where κ dy,1 t describes the change of the mean photon number and σ is the initial photon number standard deviation.As the ratio scales as the semiclassical replacement is correct as long as the ratio κ dy,1 t/n is small.As κ dy,1 (t) is of the order of light-matter interaction times time, we conclude that the semiclassical approximation gives an error of the order gt/n as indicated in Eq. (A31).Moreover, the standard deviation σ is required to be small compared to the mean photon number n to ensure the photon number independence of the matrix elements.This analysis thus justifies the first and second error scalings in Eq. (A31).
(ii) Expansion contributions.Here, we examine the magnitude of the term M (1) χ , which contributes the lowest order correction to the semiclassical momentgenerating function in Eq. (A31).The analysis of the terms M (l>1) χ works along the same lines and gives the same estimate.Transforming the sum over the photon number n into an integral, we find To make progress, we represent the auxiliary function F t,χ defined in Eq. (A26) as an exponential F t,χ (φ) = e ift,χ(φ) . (A45) In doing so, the first-order correction of the momentgenerating function can be written as where f ′ t (φ) and f ′′ t (φ) denote the first and second derivatives of f t,χ (φ) with respect to the counting field at χ = 0. Using Eqs.(A44) and (A45), we can evaluate the contribution of M (1) χ to the probability distribution to be The terms O(χ 3 ) in the exponent in Eq. (A46) generate terms of order O 1 σ 4 .Inspection of the probabilities p (1) n reveals that they are small if σ 2 ≪ f ′ t (φ) and σ 2 ≪ f ′′ t (φ)t.Both f ′ t (φ) and f ′′ t (φ) are defined via the logarithm of F t,χ (φ), which is defined in Eq. (A26).As the time evolution operators are an exponential of the Hamiltonian, we conclude that both f ′ t (φ) and f ′′ t (φ) scale with the product of the light-matter coupling g and time t.Consequently, we can estimate that the error magnitude scales as p , we find the same error scaling, such that we can conclude i.e., all terms l ≥ 1 can be neglected in the large σ limit.
As the moment-generating function can be expressed as the Fourier transformation of the probabilities, we arrive at the third error scaling given in Eq. (A31).

Generalizations
The expression of the moment-generating function in Eq. (A31) can be generalized to more general initial states along the same lines as in Appendix A 3. As the notation is tedious, we state only the final result here.
We consider a generic light-matter initial state of the form where we make no further assumption about the matter initial state |ϕ λ (t 0 )⟩.The c λ are expansion coefficients.
In contrast to the initial state in Eq. (A19) considered before, we here allow for an entangled initial state.The photonic states are Gaussian states and parameterized as where the expansion coefficients in the Fock basis |n⟩ can be written as (A51) Thereby, n λ = (n λ,1 , . . ., n λ,R ) is the vector of the initial mean photon numbers, Σ λ is a Hermitian matrix describing the corresponding standard deviation, and φ λ = (φ λ,1 , . . ., φ λ,R ) denotes the phases of the R photon modes.This parametrization covers coherent photonic states, number-squeezed photonic states, phasesqueezed photonic states, and multimode entangled photonic states.A suitable linear combination of the Gaussian states also allows for the description of Fock states.However, this would be numerically very expensive and eradicate the simplicity of the PRFT.
Generalizing the derivations in Appendix A 3 with regard to the initial state in Eq. (A49), we finally arrive at the generic moment-generating function where the initial moment-generating function is given by The consequences of the initial light-matter entanglement in Eq. (A49) may lead to new dynamical effects, whose analysis would exceed the scope of this paper.

Standard classical derivation
Here, we show that the cumulant-generating function in Eq. ( 13) reproduces the standard semiclassical definition of the energy current as considered in, e.g., Refs.[94,95].We first consider the cases of a single driving mode.Using the definition of the first dynamical cumulant in Eq. ( 12), we find + c.c.
To evaluate the derivative, we use the identity from which we can easily show that Consequently, Deriving with respect to time, we readily find the energy current from the matter system to the photon mode which is the expression of the semiclassical energy current commonly used in the literature.The generalization of Eq. (A58) to multiple modes can be easily performed by requantizing the photon modes, which are not counted: Assuming we are interested in the counting statistics of mode k, we only quantize the modes k ′ ̸ = k in the semiclassical Hamiltonian in Eq. ( 3), such that where H k (t) = gH k cos(ω k t + φ k ).This is just the Hamiltonian in Eq. ( 3) with a complicated H 0 and a single semiclassical photon mode.Using the same reasoning as before, we can directly conclude that the photon flux into mode k is given as This implies that the photon-number change in mode k is captured by the operator which is thus an operator in space and time.Clearly, an analysis of the transport dynamics in terms of ∆ Nk does not offer the convenience of the PRFT.
photons in photonic mode â between times t 0 and t 1 .
The system is assumed to be initially in state ρ(t 0 ).The photon number change shall be determined by two-point projective measurements as sketched in Fig. 1(a).At time t 0 , we perform a projective measurement defined by the projector Pn0 = |n 0 ⟩ ⟨n 0 |, where |n 0 ⟩ denote the Fock states of â, to determine the initial photon number n 0 .The resulting density operator of the light-matter system is given as where p n0 (t 0 ) denotes the probabilities to find n 0 photons.The density matrix ρ pr,n0 is the system state conditioned on the measurement outcome n 0 .We define p n|n0 (t 1 ) as the probability distribution to measure n photons by a projective measurement at time t 1 , given that there have been n 0 photons at time t 0 .This conditional probability distribution can be formally written as where N = â † â.The time-evolved density matrix is given as ρ pr,n0 (t) = Û (t)ρ pr,n0 Û † (t).For the density matrix conditioned on the first projective measurement, we can replace ρ pr,n0 = ρ pr,n0 e iχn0 e −iχ N . (B3) We note that this step is the reason why we have applied a projective measurement at time t 0 .In doing so, the conditional probability distribution can be written as We are interested in the probability distribution of the photon number change ∆n = n − n 0 , that we denote as p ∆n .This distribution is the average of the conditional probability distribution p n0+∆n|n0 weighted by the initial probability distribution, i.e., In the second equality, we have expressed the timeevolution operator in terms of the Hamiltonian We emphasize that the projective measurement in Eq. (B1) destroys coherences in the photon basis and thus modifies the initial state.If the initial state ρ(t 0 ) is already diagonal in the photon basis, such as the vacuum state or thermal states, it remains unchanged and the standard FCS makes correct predictions.For this reasons, the standard FCS correctly describes spontaneous photon emission.
In contrast, coherences in the number basis are destroyed due to the projection at time t 0 in Eq. (B1).When one blindly replaces H Q,χ → H χ (t) in Eq. (B6) with the semiclassical Hamiltonian H χ (t) to mimic a coherent field, the resulting expression is not equivalent to the dynamical cumulant generating function of the PRFT in Eq. ( 14) and thus makes wrong predictions about the photon statistics.In Figs.7-9, we depict more benchmark calculations of the PRFT for the single-mode and two-mode Rabi models.Here we shortly discuss the agreement of the PRFT to the exact quantum calculation, while the detailed interpretation of the physical effects exceeds the scope of this paper.
Sambe space simulation.In Fig. 7, we investigate the accuracy of the PRFT for the same parameters as in Fig. 2, but with the photonic subsystem represented in the Sambe space [introduced in Eq. ( 4)] instead of the Fock space.In doing so, we can assess the impact of the photon number dependence of the operators âk and â † k on the accuracy.Here we shortly discuss the major differences between Figs. 2 and 7: The variance change in Fig. 7(a) predicted by PRFT agrees now perfectly to the numerical simulation.Note that the dynamics for σ 2 = n(t 0 ) lies exactly under the other two curves.The probability distributions in Fig. 7(b) of the PRFT and the numerical simulations agree perfectly to each other.The minor oscillations of the numerical simulation in Fig. 2(b) for σ 2 = n(t 0 ) have completely disappeared, as they are a consequence of the photon number dependence of âk and â † k .The error in Fig. 7(c) is significantly reduced compared to Fig. 2(c), especially for larger σ values, demonstrating that the error is caused by the photon number dependence of âk and â † k .Moreover, the error Fig. 7(c) decreases significantly faster as function of n(t 0 ) than in Fig. 2(c).We note that the error is bounded by d 1 > 10 −5 -10 −6 due to the computer precision in the numerical simulation.
Mean photon number.In Fig. 8 we investigate the PRFT for three initial mean photon numbers n k (t 0 ) = 50, n k (t 0 ) = 500, and n k (t 0 ) = 5000 for the two photon modes k = 1, 2. The standard derivations are equally σ k = σ = 4.We observe that the PRFT time evolutions n k (t 0 ) = 500 and n k (t 0 ) = 5000 agree well to the exact numerical ones.We thus conclude that the error is mainly determined by the small σ.
Mixed benchmarking.In Fig. 9 we carry out more benchmark calculations for different photon standard deviations σ k = σ for k = 1, 2 in various parameters regimes, namely, in the dephasing, adiabatic, and highfrequency driving regime.We observe an overall precise agreement of the photonic probability distribution predicted by the PRFT to the exact numerical calculation in the full quantum model.Minor deviations can be explained by numerical fluctuations.Interestingly, in the high-frequency driving regime ω k ≫ h z , g k in Fig. 9(c), we find a perfect agreement of the PRFT and exact quantum calculations.However, in this regime there is no photon flux between the photon modes as the two-level system cannot be resonantly excited.

Jaynes-Cummings model
Here we study the photonic dynamics in a solvable model of atom-photon interactions -the Jaynes-Cummings model.This model describes a single twolevel atom interacting with a near-resonant cavity mode of the photonic field [131,132]: where σ± = σx ± iσ y and σz represents the usual Pauli matrices for spin-1/2 particles.Assuming the photon field to be classical, the corresponding semiclassical Hamiltonian describing the atomic subsystem is: where g = gα with α being the amplitude of coherent photonic state.The generalized time-evolution operator is then given by Eq. ( 6) where H χ (t) can be obtained from Eq. (C2) by replacing σ± → σ± e ±iχ .In the interaction picture defined by U 0 (t) = exp −i  The photon-resolved time-evolution operators can be obtained by applying Eq. ( 8) and are given in Eq. (41).

Two-mode Jaynes-Cummings model
The two-mode generalization of the Jaynes-Cummings model [133,134] allows for an analytical calculation of  where φ = (φ 1 , φ 2 ).In the approximate expression, we have neglected the counting fields in σφ+χ → σφ , as it only accounts for sub-leading contributions in time.
Cumulants.For the approximated moment-generating function in Eq. (C7), the moments are: 2l = (−1) l 2g 1 g 2 E φ sin(φ)t 2l which grow as t 2l+1 and t 2l in time, respectively.In this form, we can study the dependence of the transport processes on the relative phase φ = φ 2 − φ 1 .Interestingly, only the odd moments depend on the expectation value of σ φ .When the initial state is a Floquet state, the mean and the variance changes of the photon number distribution are: For φ ̸ = {0, π}, there is a net photon flux between the photon modes.The variance vanishes up to minor temporal fluctuations.According to Eq. (C9), the direction and magnitude of the photon flux can be controlled by the relative phase φ.In agreement with the analysis of the probability redistribution, the photon flux can be controlled by the initial state.The variance change vanishes exactly as predicted by Eq. (30).
The situation changes dramatically when the initial state is a balanced superposition of two Floquet states |ϕ sp ⟩ = (|u which shows that the mean photon flow between the two modes approximately vanishes, while the variance increases quadratically in time.The variance change in Eq. (C10) is equal to the squared mean change in Eq. (C9).This time dependence is easily understood from Fig. 4(a), as the rapid variance increase is a consequence of the linearly growing distance between the two Floquet-state-dependent peaks.
Light-matter entanglement.Here, we apply the PRFT to calculate the purity in the Jaynes-Cummings model in order to describe the light-matter entanglement as a function of time.According to the PRFT and Eq. ( 32), the time-evolved state can be written as |Ψ(t)⟩ = c 1 e −iE1,φt |u 1,φ ⟩ |A 1 ⟩ + c 2 e −iE2,φt |u 2,φ ⟩ |A 2 ⟩ , (C11) with the in general non-orthogonal photonic states |A µ ⟩.The photon probability distribution conditioned on the Floquet state approximately reads as p n k |µ (t) = e (n k −n k,µ (t))/2σ 2 /( √ πσ), where σ is the width, and n k,µ (t) = n k (0)+κ PRFT.In this case, the overlap υ is completely determined by the conditional probabilities p n k |µ (t).Considering the photon number as a continuous variable, we can

FIG. 1 .
Photon distribution Extending validity of Floquet theory Macroscopic light-matter entanglement Decoherence in Floquet basis Efficient quantum communication protocol

FIG. 2 .
FIG. 2. (a) Time evolution of the spin components and the first two photonic cumulants in the quantum Rabi model in the intermediate-coupling regime ω = hz = g.The initial condition is given in Eq. (36) with |ϕ(t0)⟩ = (|↓⟩ + |↑⟩)/ √ 2 for n(t0) = 10 5 photons, as well as for σ = n(t0) (red, dash-dotted) and σ = 0.1 n(t0) (yellow,solid).(b) Probability distribution at the scaled time tσ ≡ (2π/g)6σ for the same coherent (right) and number-squeezed (left) initial states as in (a) and for the resonance condition hz = ω.(c) Analysis of the trace distance in Eq. (40) as a measure of the error at the scaled time tσ as a function of noise σ (left) and mean photon number n(t0) with scaled σ = 0.3 n(t0) and tσ (right).

FIG. 3 .
FIG. 3. Benchmark calculations of the PRFT for an initial Floquet state and various photon number widths σ1 = σ2 in the two-mode Rabi model.The first column depicts the trace distance d1(pn 1 , p

FIG. 6 .
FIG.6.Remote entanglement generation protocol.(a) Illustration of the local light-matter entanglement process.After the interaction of two coherent light fields at a GHZ state, the final state is a superposition of states where either coherent field is enhanced (thick arrow) while the other is reduced (thin arrow).(b) Entanglement protocol outlined in Sec.V B. After light-matter entanglement with Alice's GHZ state, the light is transmitted to Bob, where it interacts with Bob's GHZ state.The dashed lines sketch the light paths.After Bob's lightmatter entanglement, the coherent fields are measured (gray half circles).A vanishing signal difference heralds success of the protocol.The three insets show the photon distributions pn before and after Alice's entanglement process, and after Bob's process.(c) Analysis of the which-path information.After interaction with Bob's GHZ state, the system is in a superposition of four states.In two of which, the two output fields have changed their amplitude which reveals the path of quantum information and results in failing of the protocol.In the other two states, the which-path information remains hidden and results in success of the protocol.
where the photonic states |A 2− ⟩, |A 0 ⟩, |A 2+ ⟩ are close to coherent states with amplitudes smaller, comparable and larger compared to the initial coherent states.The change of amplitude refers to mode k = 1, while mode k = 2 will conversely have a larger, comparable or smaller amplitude.The photon distribution of the modes at different stages of the protocol is sketched in Fig.6(b).

3 .
Measurement and postselection: Bob makes a projective measurement defined by the operator |A 0 ⟩ ⟨A 0 |.If the measurement is successful, Alice and Bob carry out Hadmard gates (and other local operations to correct for phases accumulated during the light-matter interactions).Finally, Alice and Bob hold a share of the entangled state

m
in the Fock basis of the photonic initial states m a (k) m |m⟩ k : The amplitudes a (k) m vary slowly with the photon number m; The phases are well defined in the sense that arg a (k) m = φ k m with constant φ k .These assumptions are naturally fulfilled for the coherent state in Eq. (2) when α k is large, such that a

σ 2 .
Carrying out a similar analysis for M (l>1) χ Appendix C: Details to the application of the model calculations 1. Benchmark calculations

FIG. 8 . 1 ,
FIG.8.Same as in Fig.3, but for different initial mean photon numbers n1 = n2 as indicated in the panels.

FIG. 9 .
FIG. 9. Same as in Fig.3, but in different parameter regimes as indicated in the panels.

| 2 c * 1 c 2 υ c * 2 c 1 υ * |c 2 | 2 ,
FIG. 10.Purity as an entanglement measure of the lightmatter state in the two-mode Jaynes-Cummings model as a function of time: (a) shows the purity for an initial Floquet state for different initial photon-number variances σ 2 = Var N1(t0), and (b) shows the purity for a balanced superposition of Floquet states as initial state.The black and red lines depict the quantum simulation and PRFT respectively.Overall parameters are hz = ω = 10g and n k (t0) = 10 6 for k = 1, 2.