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

Georg Engelhardt<EMAIL_ADDRESS>Shenzhen Institute for Quantum
Science and Engineering, Southern University of Science and Technology,
Shenzhen 518055, China  International Quantum Academy, Shenzhen 518048, China
Guangdong Provincial Key Laboratory of Quantum Science and Engineering,
Southern University of Science and Technology, Shenzhen, 518055, China  Sayan
Choudhury<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Pittsburgh, Pittsburgh, PA 15260, USA  W. Vincent Liu
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Pittsburgh,
Pittsburgh, PA 15260, USA  Shenzhen Institute for Quantum Science and
Engineering, Southern University of Science and Technology, Shenzhen 518055,
China

###### Abstract

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 is dubbed “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 multi-frequency 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. Intriguingly, time crystals can
act as good quantum memories even in the optical 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
new 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] and discrete time crystals [16, 17, 18,
19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31] as well as breakthroughs
in applications like spectroscopy [32, 33, 34, 35], metrology [36, 37, 38],
and quantum simulation [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]. These
non-equilibrium quantum systems are generally analyzed using Floquet theory -
a method first developed by Jon Shirley in 1965 [51]. 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
[52, 53]. 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 [54]. In general,
a semiclassical approach is used instead, where the photon fields are assumed
to be high-energy coherent states, and their dynamics is neglected. This leads
to an effective Floquet description of the atomic dynamics, which can be
employed to engineer materials with novel emergent properties [55, 56, 57, 58,
59, 60, 61, 62, 63, 64]. Unfortunately, while this semiclassical treatment is
very powerful in modeling the matter subsystem, it fails to describe the
photonic driving field.

In this paper, we take a significant step beyond the semiclassical Floquet
theoretic description of light-matter interactions by developing a new
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 [65, 66, 67], FCS is a
powerful method that has been employed to study mesoscopic transport [68, 69,
70, 71, 72, 73, 74, 75], quantum dots [76, 77, 78, 79, 80], spin chains [81,
82], spontaneous photon emission [83, 84, 85, 86], thermodynamics [87, 88],
and the entanglement entropy of non-interacting fermions [89, 90, 91].
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 [92].
Unfortunately, such measurements destroy the coherent photonic states, thereby
rendering the Floquet description of the matter 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. 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 [93,
94, 95], these approaches scale exponentially with the number of frequency
modes; the PRFT does not suffer from this limitation.

We employ the PRFT to analyze multi-mode driving, and discover that this leads
to macroscopic light-matter entanglement at long times. This light-matter
entanglement causes a complete decoherence of the matter system in the basis
of the Floquet eigenstates. The PRFT thus provides a quantum-optical
interpretation of Floquet states as the decohering basis for the matter
system. Furthermore, this decoherence has far-reaching experimental
consequences. In particular, we demonstrate that the quantum-optical coherence
time is reasonably short (a few ms) for typical optical electromagnetic 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\,\text{kHz}$ regime over
$500\,\text{km}$, thereby far exceeding the Hz regime that is predicted in
current theoretical protocols [96]. 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 PRFT and
benchmark it for the Jaynes-Cummings and quantum Rabi models. In Sec. III, we
analyze the theoretical implications for periodically-driven quantum systems
and illustrate them for multi-mode versions of the Jaynes-Cummings and quantum
Rabi models. 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

###  II.1 Formalism

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.

System: We consider the following generic Hamiltonian describing a matter
system interacting with a multi-mode photonic field:

$H_{Q}=H_{0}+\sum_{k}\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}+\sum_{k=1}^{R}\tilde{g}H_{k}\left(\hat{a}_{k}^{\dagger}+\hat{a}_{k}\right),$
(1)

where $k$ denotes the different photonic modes, $\hat{a}_{k}$ are annihilation
operators quantizing the photonic modes and $H_{k}$ acts on the matter system.
The light-matter interaction strength is parameterized by $\tilde{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 $\hat{N}_{k}=\hat{a}_{k}^{\dagger}\hat{a}_{k}$, i.e.,
$\hat{N}_{k}\left|n\right>_{k}=n\left|n\right>_{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 the system is initially in the state

$\displaystyle\left|\psi(t_{0})\right>=\left|\phi(t_{0})\right>\otimes\prod_{k=1}^{R}\left|\alpha_{k}e^{i\varphi_{k}}\right>,$
(2)

where $\left|\alpha_{k}e^{i\varphi_{k}}\right>$ are coherent states of the
photon operators $\hat{a}_{k}$ with real-valued amplitudes $\alpha_{k}>0$ and
phases $\varphi_{k}\in\left[0,2\pi\right)$, and the state
$\left|\phi(t_{0})\right>$ is the initial state of the matter system [53, 51].
In this semiclassical limit, we can substitute
$\hat{a}_{k}\rightarrow\alpha_{k}e^{-i\varphi_{k}-i\omega_{k}t}$ in Eq. (1)
such that we obtain the corresponding semiclassical Hamiltonian

$\mathcal{H}(t)=H_{0}+\sum_{k=1}^{R}H_{k}g_{k}\cos(\omega_{k}t+\varphi_{k}),$
(3)

where we have introduced the effective light-matter interactions
$g_{k}=\tilde{g}\alpha_{k}$. This description is valid as long as the back
action of the quantum system on the photonic field is negligible, i.e., if
$\tilde{g}\ll\omega_{k}\alpha_{k}$. Nevertheless the impact of the photonic
field on the matter system can be large because of the product
$\tilde{g}\alpha_{k}$. The semiclassical approach is thus valid for large mean
occupation numbers $\langle\hat{N}_{k}\rangle=\alpha_{k}^{2}$.

For a single photonic frequency $\omega_{k}$ or for commensurate frequencies
(where all $\omega_{k^{\prime}}/\omega_{k}$ are rational numbers), the
semiclassical Hamiltonian is time periodic
$\mathcal{H}(t)=\mathcal{H}(t+\tau)$ with a driving period $\tau$. Under these
conditions one can apply the celebrated Floquet theory to analyze the system
[51]. 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.

Photon-resolved time-evolution: The time-evolution operator $\mathcal{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 real-valued counting fields
$\chi_{k}\in\left[0,2\pi\right)$, and define the generalized time-evolution
operator:

$\displaystyle
U_{\boldsymbol{\chi}}(t)=\exp(i\sum_{k}\chi_{k}\hat{N}_{k})U(t)\exp(-i\sum_{k}\chi_{k}\hat{N}_{k}).$
(4)

This transformation implies that the annihilation (creation) operators
transform as $\hat{a}_{k}\rightarrow\hat{a}_{k}e^{-i\chi_{k}}$
($\hat{a}^{\dagger}_{k}\rightarrow\hat{a}^{\dagger}_{k}e^{i\chi}_{k}$),
leading to the following generalized time-evolution operator in the
semiclassical limit:

$\displaystyle\mathcal{U}_{\boldsymbol{\chi}}(t)$ $\displaystyle\equiv$
$\displaystyle\mathcal{T}e^{-i\int_{0}^{t}\mathcal{H}_{\boldsymbol{\chi}}(t^{\prime})dt^{\prime}},$
(5) $\displaystyle\mathcal{H}_{\boldsymbol{\chi}}(t)$ $\displaystyle=$
$\displaystyle
H_{0}+\sum_{j}^{R}H_{k}\alpha_{k}\left(e^{i\omega_{k}t+i\chi_{k}}+e^{-i\omega_{j}t-i\chi_{k}}\right),$
(6)

where $\boldsymbol{\chi}=\left(\chi_{1},\dots,\chi_{R}\right)$ is a vector of
counting fields (see Appendix A for a detailed derivation). Based on the
generalized time-evolution operator in Eq. (5), we define the photon-resolved
time-evolution operators

$\mathcal{U}^{(\boldsymbol{m})}(t)=\frac{1}{\left(2\pi\right)^{R}}\iiint_{0}^{2\pi}d\boldsymbol{\chi}\mathcal{U}_{\boldsymbol{\chi}}(t)e^{-i\boldsymbol{m}\boldsymbol{\chi}},$
(7)

where $\boldsymbol{m}=(m_{1},\dots,m_{R})$ is a vector of photonic transition
numbers. In terms of the photon-resolved time-evolution operator, we can now
express arbitrary system observables. For example, the probability that the
photonic modes are in the Fock states $\boldsymbol{n}=(n_{1},\dots,n_{R})$ is
given as

$\displaystyle\left<\hat{P}_{\boldsymbol{n}}\right>_{t}$ $\displaystyle=$
$\displaystyle\sum_{\boldsymbol{m},\boldsymbol{m}^{\prime}}\left<\mathcal{U}^{(\boldsymbol{m})\dagger}(t)\mathcal{U}^{(\boldsymbol{m}^{\prime})}(t)\right>_{t_{0}}$
(8) $\displaystyle\times$
$\displaystyle\prod_{k}a_{n_{k}-m_{k}}^{(k)*}a_{n_{k}-m_{k}^{\prime}}^{(k)},$

where $\hat{P}_{\boldsymbol{n}}$ denotes the projector onto the Fock states
with quantum numbers $\boldsymbol{n}$, $a_{m}^{(k)}$ are the expansion
coefficients of the photonic initial state
$\left|\alpha_{k}e^{i\varphi_{k}}\right>=\sum_{m}a_{m}^{(k)}\left|m\right>_{k}$,
and
$\langle\hat{A}\rangle_{t}\equiv\langle\psi(t)\mid\hat{A}\mid\psi(t)\rangle$,
where $\left|\psi(t)\right>$ can be either a state in the matter system or the
composite light-matter system depending on the enclosed operator $\hat{A}$.
Thus, by semiclassically calculating the generalized time-evolution operators
in Eq. (5), we can evaluate genuine quantum properties of photonic
observables.

The derivations in the PRFT make the following assumptions for the expansion
coefficients $a_{m}^{(k)}$ in the Fock basis of the photonic initial states
$\sum_{m}a_{m}^{(k)}\left|m\right>_{k}$: (i) The amplitudes
$\left|a_{m}^{(k)}\right|$ vary slowly with the photon number $m$; (ii) The
phases are well defined in the sense that $\arg a_{m}^{(k)}=\varphi_{k}\cdot
m$ with constant $\varphi_{k}$. These assumptions are naturally fulfilled for
the coherent state in Eq. (2) when $\alpha_{k}$ is large, such that
$\left|a_{m}^{(k)}\right|^{2}$ obeys the Poisson distribution. Moreover, the
PRFT framework can also be extended to include other kind of initial states
(like number-squeezed states), as long as conditions (i) and (ii) are
fulfilled.

Figure 1: Schematic illustration of the fundamental difference between FCS and
the PRFT. (a) The standard FCS extracts the statistical information based on
two-point projective measurements at times $t_{0}$ and $t_{1}$. (b) The
counting statistics in the PRFT relies on two-point tomographic measurements
of two independent batches. One of which is used to determine the statistics
at time $t_{0}$, and the other determines the statistics of the time-evolved
state $t_{1}$ without being measured at $t_{0}$.

Full-Counting Statistics: While the analytical (or numerical) evaluation of
the photon-resolved time-evolution operators in Eq. (7) is an infeasible task
in many cases, it is relatively easy to compute the moments and cumulants of
the photonic modes $\hat{a}_{k}$. Importantly, the counting statistics within
the PRFT is fundamentally different from the standard FCS formalism; this is
illustrated in Fig. 1. 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 [92, 97, 68] (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 photon-number 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. [98].

Formally, the counting statistics of the photon modes can be calculated via
the cumulant- or the moment-generating functions, which are defined by

$\displaystyle K_{\boldsymbol{\chi}}(t)$ $\displaystyle=$
$\displaystyle\log\left(M_{\boldsymbol{\chi}}(t)\right),$ (9) $\displaystyle
M_{\boldsymbol{\chi}}(t)$ $\displaystyle=$
$\displaystyle\left<e^{i\sum_{k}\chi_{k}\hat{N}_{k}}\right>_{t},$ (10)

respectively. The associated $n$-th cumulant and moment of mode $k$ are
determined via

$\displaystyle\kappa_{n}^{(k)}(t)$ $\displaystyle=$
$\displaystyle\frac{d^{n}}{d\,(i\chi_{k})^{n}}K_{\boldsymbol{\chi}}(t),$
$\displaystyle m_{n}^{(k)}(t)$ $\displaystyle=$
$\displaystyle\frac{d^{n}}{d\,(i\chi_{k})^{n}}M_{\boldsymbol{\chi}}(t).$ (11)

We are interested in the time-evolution of the cumulant- and moment-generating
functions. To this end, we define the dynamical cumulant-generating function
$K_{dy,\boldsymbol{\chi}}(t_{1})$ via

$\displaystyle K_{\boldsymbol{\chi}}(t_{1})$ $\displaystyle=$ $\displaystyle
K_{dy,\boldsymbol{\chi}}(t_{1})+K_{\boldsymbol{\chi}}(t_{0}),$ (12)

where $K_{\boldsymbol{\chi}}(t_{0})$ and $K_{\boldsymbol{\chi}}(t_{1})$ can be
determined by independent tomographies at the beginning and the end of the
time evolution [98]. As shown in detail in Appendix A.3, the dynamical
cumulant-generating function can be expressed as

$K_{dy,\boldsymbol{\chi}}(t)=\log\frac{1}{2}\left<\mathcal{U}_{\boldsymbol{\varphi}}^{\dagger}(t)\mathcal{U}_{\boldsymbol{\varphi}+\boldsymbol{\chi}}(t)+\mathcal{U}_{\boldsymbol{\varphi}-\boldsymbol{\chi}}^{\dagger}(t)\mathcal{U}_{\boldsymbol{\varphi}}(t)\right>_{t_{0}},$
(13)

where $\boldsymbol{\varphi}=\left(\varphi_{1},\dots\varphi_{R}\right)$ is the
vector of phases of the photonic states. Employing Eq. (12), we can now obtain
the change of the cumulants:

$\kappa_{dy,n}^{(k)}(t)\equiv\left.\frac{d^{n}}{d(i\chi_{k})^{n}}K_{dy,\boldsymbol{\chi}}(t)\right|_{\boldsymbol{\chi}=0}=\kappa_{n}^{(k)}(t)-\kappa_{n}^{(k)}(t_{0}),$
(14)

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

$\displaystyle\Delta\langle N_{k}(t)\rangle$ $\displaystyle=$
$\displaystyle\kappa_{1}^{(k)}(t)-\kappa_{1}^{(k)}(t_{0}),$
$\displaystyle\Delta\sigma_{k}^{2}(t)$ $\displaystyle=$
$\displaystyle\kappa_{2}^{(k)}(t)-\kappa_{2}^{(k)}(t_{0})$ (15)

correspond to the change of the mean photon number $\langle N_{k}\rangle$ and
the variance $\sigma_{k}^{2}$, 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.

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 moment-generating function is
$\tilde{M}_{\boldsymbol{\chi}}(t)=\langle\mathcal{U}_{\boldsymbol{\varphi}-\boldsymbol{\chi}/2}^{\dagger}(t)\mathcal{U}_{\boldsymbol{\varphi}+\boldsymbol{\chi}/2}(t)\rangle_{t_{0}}$
according to the standard FCS, where the phases $\boldsymbol{\varphi}$ have
been evaluated heuristically. This results is in profound differences with the
predictions of the PRFT for several important quantities (like the second
cumulant), and highlights the crucial role played by two-point tomographic
measurements in the PRFT.

Probability redistribution: Using the moment-generating function of the final
state given in Eq. (10), we can calculate the probability distribution of the
Fock states at time $t$ via

$\displaystyle p_{\boldsymbol{n}}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{R}}\int_{0}^{2\pi}M_{\boldsymbol{\chi}}(t)e^{-i\boldsymbol{\chi}\cdot\boldsymbol{n}}d\boldsymbol{\chi}$
(16) $\displaystyle=$ $\displaystyle\sum_{m}q_{n-m}(t)p_{m}(t_{0}),$

where we have introduced the  quasiprobabilities  as

$q_{\boldsymbol{n}}(t)=\frac{1}{(2\pi)^{R}}\int_{0}^{2\pi}M_{dy,\boldsymbol{\chi}}(t)e^{-i\boldsymbol{\chi}\cdot\boldsymbol{n}}d\boldsymbol{\chi},$
(17)

where
$M_{dy,\boldsymbol{\chi}}(t)=\exp\left[K_{dy,\boldsymbol{\chi}}(t)\right]$.
Akin to regular probabilities, these quasiprobabilities are real valued and
fulfill

$\sum_{\boldsymbol{n}}q_{\boldsymbol{n}}=1.$ (18)

However, due to interference effects, $q_{\boldsymbol{n}}$ may also be
negative. Such negativity is a characteristic signature of nonclassical
temporal correlations [99]. We note that according to Eq. (16), the
quasiprobabilities may be interpreted as the kernel of the probability
redistribution.

###  II.2 Benchmark Calculations

We now proceed to apply the PRFT to two paradigmatic models of atom-photon
interactions in a cavity, and demonstrate that this framework is exact in the
semiclassical limit.

Figure 2: (a) Time evolution of the spin components and the first two photonic
cumulants in the quantum Rabi model. The time evolution calculated in the Fock
space (dashed) agrees perfectly with the prediction of the PRFT (solid).
Parameters are $\omega=10h_{z}$ and $g=10h_{z}$ and the initial mean photon
number is $\langle\hat{N}\rangle_{0}=1000$. (b) Quasiprobabilities $q_{\Delta
n}$ [defined in Eq. (17)] and probability redistribution $\Delta
p_{n}=p_{n}(t)-p_{n}(0)$ at time $t=1.2\cdot(2\pi/\omega)$ for
$\omega=g=h_{z}$ (top row) and $\omega=g=10h_{z}$ (bottom row). The negativity
of the quasiprobabilities is more pronounced in the latter case due to virtual
processes that arise in this regime.

Jaynes-Cummings model: A natural starting point for our analysis is to study
the photonic dynamics in a solvable model of atom-photon interactions — the
Jaynes-Cummings model. This model describes a single two-level atom
interacting with a near resonant cavity mode of the photonic field [100, 101]:

$H_{JC}=\frac{h_{z}}{2}\hat{\sigma}_{z}+\omega\hat{a}^{\dagger}\hat{a}+\tilde{g}\left(\hat{\sigma}_{+}\hat{a}+\hat{\sigma}_{-}\hat{a}^{\dagger}\right),$
(19)

where $\hat{\sigma}_{\pm}=\hat{\sigma}_{x}\pm i\hat{\sigma}_{y}$ and
$\hat{\sigma}_{z}$ represents the usual Pauli matrices for spin-$1/2$
particles. Assuming the photon field to be classical, the effective
Hamiltonian describing the atomic subsystem is:

$\mathcal{H}(t)=\frac{h_{z}}{2}\hat{\sigma}_{z}+g\left(\hat{\sigma}_{+}e^{-i\omega
t}+\hat{\sigma}_{-}e^{i\omega t}\right),$ (20)

where $g=\tilde{g}\alpha$ with $\alpha$ being the amplitude of coherent
photonic state. The generalized time-evolution operator is then given by Eq.
(5) where $\mathcal{H}_{\chi}(t)$ can be obtained from Eq. (20) by replacing
$\hat{\sigma}_{\pm}\rightarrow\hat{\sigma}_{\pm}e^{\mp i\chi}$. In the
interaction picture defined by
$U_{0}(t)=\exp\left[{-i\frac{\omega}{2}\hat{\sigma}_{z}t}\right]$, the
generalized time-evolution operator becomes

$\displaystyle\mathcal{U}_{\chi}(t)=\cos\left(Et\right){\mathbbm{1}}+i\sin\left(Et\right)\hat{\sigma}_{\chi},$
(21)

where

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\left(h_{z}-\omega\right)^{2}+16g^{2}},$
$\displaystyle\tan\theta$ $\displaystyle=$
$\displaystyle\frac{2g}{h_{z}-\omega},$ $\displaystyle\hat{\sigma}_{\chi}$
$\displaystyle=$
$\displaystyle\left[\cos\theta\hat{\sigma}_{z}+\sin\theta\left(\cos\chi\hat{\sigma}_{x}+\sin\chi\hat{\sigma}_{y}\right)\right].$

The photon-resolved time-evolution operators can be obtained by applying Eq.
(7) and read as

$\displaystyle\mathcal{U}^{(0)}(t)$ $\displaystyle=$ $\displaystyle
e^{-i\frac{\omega}{2}\sigma_{z}t}\cos\left(\frac{Et}{2}\right)\mathbbm{1}+i\sin\left(\frac{Et}{2}\right)\cos\theta\hat{\sigma}_{z},$
$\displaystyle\mathcal{U}^{(\pm 1)}(t)$ $\displaystyle=$ $\displaystyle
ie^{-i\frac{\omega}{2}\sigma_{z}t}\sin\left(\frac{Et}{2}\right)\sin\theta\hat{\sigma}_{\mp}.$
(22)

This form of the photon-resolved operators encodes the conservation of the
quantity: $N_{\rm tot}=\hat{a}^{\dagger}\hat{a}+\hat{\sigma}_{z},$ which is a
salient feature in the Jaynes-Cummings model. To illustrate this, let us
consider the initial state
$\left|\psi(t_{0})\right>=\left|\uparrow\right>\otimes\left|n\right>,$ that
gives rise to the following photonic occupations of the Fock states $n$ and
$n+1$

$\displaystyle\left<\hat{P}_{n}\right>$ $\displaystyle=$
$\displaystyle\left[\cos\left(Et\right)\right]^{2}+\left[\cos\theta\sin\left(Et\right)\right]^{2},$
$\displaystyle\left<\hat{P}_{n+1}\right>$ $\displaystyle=$
$\displaystyle\left[\sin\theta\sin\left(Et\right)\right]^{2},$ (23)

which maintains the probability
$\langle\hat{P}_{n}\rangle+\langle\hat{P}_{n+1}\rangle=1.$ Intriguingly, the
occupations in Eq. (23) are identical to the exact time evolution of the
Jaynes-Cummings model. Thus, we have determined the dynamics of a genuine
quantum model by employing semiclassical methods of the PRFT.

Rabi model: While the Jaynes-Cummings model is exactly solvable, it is valid
only when the frequency of the cavity mode is nearly resonant with the atomic
level splitting. In order to go beyond this regime, we now examine the quantum
Rabi model describing the interaction of an atom with a single-mode
electromagnetic field:

$H_{QR}=\frac{h_{z}}{2}\hat{\sigma}_{z}+\omega
a^{\dagger}a+\tilde{g}\hat{\sigma}_{x}(a^{\dagger}+a),$ (24)

where $h_{z}$ is the atomic level splitting, $\omega$ is the frequency of the
photonic mode and $\tilde{g}$ is the light-matter interaction. The Jaynes-
Cummings model emerges as the effective description of this system when
$\tilde{g}\ll\omega$, $h_{z}\ll|\omega+h_{z}|$, and $|\omega-
h_{z}|\ll|\omega+h_{z}|$ [102, 103, 104, 105]. In the semiclassical limit, the
atomic dynamics is described by
$\mathcal{H}_{R}=h_{z}\hat{\sigma}_{z}/2+2g\hat{\sigma}_{x}\cos(\omega t).$

We compare the photon statistics obtained using the PRFT with exact numerical
results obtained for the quantum Rabi model. For these calculations, we have
approximated the Poisson distribution by a Gaussian function with mean
$\overline{n}$, standard deviation $\sigma=\sqrt{\overline{n}}$, and
normalization $\mathcal{N}$. As shown in Fig. 2, we find that the PRFT
calculations reproduce the full quantum dynamics exactly. We have tested a
wide range of parameters for both $\omega$ and $g$ and have always found a
perfect agreement, even for very long times. These results clearly illustrate
that the PRFT provides a powerful semiclassical route to compute photon
counting statistics for archetypal quantum-optical models.

##  III Theoretical implications

Having established the validity of the PRFT, we now proceed to analyze systems
subjected to multi-mode driving where the photonic field is composed of
multiple commensurate photonic frequencies $\omega_{j}$. 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. Finally, we illustrate these results
for the two-mode Jaynes-Cummings and three-mode Rabi models.

###  III.1 Floquet-state analysis

According to Floquet theory, the generalized time-evolution operator can be
written as

$\mathcal{U}_{\boldsymbol{\chi}}(t)=\mathcal{U}_{k,\boldsymbol{\chi}}(t)\exp\left(-iH_{Fl,\boldsymbol{\chi}}t\right),$
(25)

where the generalized Floquet Hamiltonian can be expanded as

$H_{Fl,\boldsymbol{\chi}}=\sum_{\mu}E_{\mu,\boldsymbol{\chi}}\left|u_{\mu,\boldsymbol{\chi}}\right>\left<u_{\mu,\boldsymbol{\chi}}\right|.$
(26)

As the kick operator
$\mathcal{U}_{k,\boldsymbol{\chi}}(t)=\mathcal{U}_{k,\boldsymbol{\chi}}(t+\tau)$
is time periodic, it accounts only for periodic changes in the photon
redistribution. Thus, the asymptotic dynamics of the cumulant-generating
function in Eq. (13) is determined by the generalized quasienergies
$E_{\mu,\boldsymbol{\chi}}$. The stroboscopic dynamics of the matter system is
characterized by the Floquet states $\left|u_{\mu,\boldsymbol{\chi}}\right>$,
that generalize the common eigenstates in time-independent systems. When
expanding the initial state in the Floquet state basis
$\left|u_{\mu}\right>=\left|u_{\mu,\boldsymbol{\chi}=0}\right>$

$\left|\Psi(0)\right>=\sum_{\mu}c_{\mu}\left|u_{\mu}\right>\otimes\left|\boldsymbol{\alpha}\right>,$
(27)

where $\left|\boldsymbol{\alpha}\right>$ is the initial state of the photonic
field, the asymptotic dynamical cumulant- and moment-generating functions read
as

$\displaystyle K_{dy,\boldsymbol{\chi}}(t)$ $\displaystyle\rightarrow\log
M_{dy,\boldsymbol{\chi}}(t),$ (28) $\displaystyle M_{dy,\boldsymbol{\chi}}(t)$
$\displaystyle\rightarrow\qquad\qquad$ $\displaystyle\sum_{\mu}$
$\displaystyle\frac{\left|c_{\mu}\right|^{2}}{2}\left(e^{i\left(E_{\mu,\boldsymbol{\varphi}}-E_{\mu,\boldsymbol{\varphi}+\boldsymbol{\chi}}\right)t}+e^{i\left(E_{\mu,\boldsymbol{\varphi}-\boldsymbol{\chi}}-E_{\mu,\boldsymbol{\varphi}}\right)t}\right).$

Thus, $M_{dy,\boldsymbol{\chi}}$ is a weighted average of the dynamical
moment-generation functions of the Floquet states, with the weights given by
the expansion coefficients $\left|c_{\mu}\right|^{2}$. We can now obtain the
mean photon-number change to be:

$\Delta\langle\hat{N}_{k}\rangle=-\sum_{\mu}\left|c_{\mu}\right|^{2}\frac{dE_{\mu,\boldsymbol{\varphi}}}{d\varphi_{k}}\cdot
t=\sum_{\mu}\left|c_{\mu}\right|^{2}\Delta\langle\hat{N}_{k}\rangle\mid_{|\mu\rangle}\cdot
t,$ (29)

where we have evaluated the quasienergies at the phases of the photonic fields
$\boldsymbol{\varphi}$. Intriguingly, when preparing the system initially in
an arbitrary Floquet state $\mu$ the variance change $\Delta\sigma_{k}^{2}$
vanishes

$\left.\Delta\sigma_{k}^{2}\right|_{|\mu\rangle}=0.$ (30)

It is worth noting that the standard FCS incorrectly predicts
$\left.\Delta\sigma_{k}^{2}\right|_{\left|u_{\lambda}\right>}=\frac{d^{2}E_{\lambda,\boldsymbol{\varphi}}}{d\varphi_{k}^{2}}\cdot
t^{2}$, which is a consequence of the initial projection on the photon number
state basis.

For a superposition of Floquet states, Eqs. (17) and (28) imply that the
quasiprobabilities are given by a weighted average of Floquet-state dependent
quasiprobabilities $q_{n\mid\mu}$, i.e.,
$q_{n}=\sum_{\mu}\left|c_{\mu}\right|^{2}q_{n\mid\mu}$. Similarly, from Eq.
(16), we can infer that

$p_{n}=\sum_{\mu}\left|c_{\mu}\right|^{2}p_{n\mid\mu}\quad.$ (31)

Thus, $p_{n\mid\mu}$ defines a conditional probability distribution, and the
dynamics of the photonic state is controlled by the Floquet states.

Figure 3: (a) Schematic of the two-mode Jaynes-Cummings model. (b,c,d) Photon
distribution in mode $k=1$ for three different initial states
$\left|u_{1}\right>,\left(\left|u_{1}\right>+\left|u_{2}\right>\right)/\sqrt{2}$,
and $\left|u_{2}\right>$, where $\left|u_{\mu}\right>$ are the two Floquet
states of the semiclassical system. The time steps $t\cdot h_{z}=0,600,1200$
are depicted in brown, green, and yellow, respectively. The black solid lines
depict the simulation in Fock space, while the colored areas are the
predictions of the PRFT. (e,f,g) Quasiprobabilities defined via the dynamical
moment-generating function in Eq. (17). The initial conditions are the same as
in (b,c,d), respectively. (h,i,j) Dynamical cumulants of the photon
redistribution as a function of time for the same initial conditions as in
(b,c,d). The system parameters are $\varphi_{1}=0$ and $\varphi_{2}=\pi/2$,
and $\omega=h_{z}$, $g_{1}=g_{2}=0.2h_{z}$. The variances of the initial
photonic distributions are $\sigma_{k}^{2}=100$ for $k=1,2$.

Our results lead to a quantum-optical interpretation of Floquet states. To see
this, we consider the time-evolution 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:

$\left|\Psi(t)\right>=\sum_{\mu}c_{\mu}\left|u_{\mu}(t)\right>\otimes\left|A_{\mu}(t)\right>,$
(32)

where the mean photon number of each $\left|A_{\mu}\right>$ changes linearly
in time [see Eq. (29)], while its variance remains unchanged [see Eq. (30)].
Consequently, the photonic states $\left|A_{\mu}\right>$ will become mutually
orthogonal for sufficiently long times, and the reduced density matrix of the
matter system becomes:

$\displaystyle\rho_{matter}=\sum_{\mu}\left|c_{\mu}\right|^{2}\left|u_{\mu}(t)\right>\left<u_{\mu}(t)\right|.$
(33)

Thus, from a quantum-optical point of view, Floquet states acts as the
decohering basis. We note that this interpretation holds as long as there is a
non-vanishing photon flux between distinct photonic modes.

###  III.2 Two-mode Jaynes-Cummings model

The single-mode models in Sec. II.2 cannot exhibit large-scale transport. To
illustrate the theoretical findings in Sec. III.1, we investigate the photon
counting statistics in the two-mode generalization of the Jaynes-Cummings
model [106, 107]:

$H_{TMJC}=\frac{h_{z}}{2}\hat{\sigma}_{z}+\sum_{k=1}^{2}\omega\hat{a}_{k}^{\dagger}\hat{a}_{k}+\sum_{k=1}^{2}\tilde{g}_{k}\left(\hat{\sigma}_{+}\hat{a}_{k}+\hat{\sigma}_{-}\hat{a}_{k}^{\dagger}\right),$
(34)

where both modes have the same frequency $\omega_{k}=\omega$, and the photonic
modes are initially in coherent states
$\left|\alpha_{k}e^{i\varphi_{k}}\right>$ with real valued $\alpha_{k}$ and
$\varphi_{k}$. A schematic illustration of this model is shown in Fig. 3(a).

Moment-generating function: Following the procedure in Sec. II.2, we proceed
to calculate the moment-generating function. The semiclassical Hamiltonian
including the counting fields is given as

$\displaystyle\mathcal{H}_{\boldsymbol{\chi}}(t)=\frac{h_{z}}{2}\hat{\sigma}_{z}+\left[\hat{\sigma}_{+}e^{-i\omega
t}G(\boldsymbol{\chi})+\hat{\sigma}_{-}e^{i\omega
t}G^{*}(\boldsymbol{\chi})\right],$

where we have defined $G(\boldsymbol{\chi})=\sum_{k}g_{k}e^{-i\chi_{k}}.$ In
an interaction picture defined by
$U_{0}(t)=e^{-i\frac{\omega}{2}\hat{\sigma}_{z}t}$, the time-evolution
operator reads as

$\displaystyle\mathcal{U}_{\boldsymbol{\chi}}(t)$ $\displaystyle=$
$\displaystyle e^{-iH_{\boldsymbol{\chi}}t}$ (35) $\displaystyle=$
$\displaystyle\cos\left(E_{\boldsymbol{\chi}}t\right)\mathbbm{1}+i\sin\left(E_{\boldsymbol{\chi}}t\right)\hat{\sigma}_{\boldsymbol{\chi}},$

where
$H_{\boldsymbol{\chi}}=\frac{h_{z}}{2}\hat{\sigma}_{z}+\left(\hat{\sigma}_{+}G(\boldsymbol{\chi})+\hat{\sigma}_{-}G^{*}(\boldsymbol{\chi})\right)$
is the Floquet Hamiltonian and

$\displaystyle\hat{\sigma}_{\boldsymbol{\chi}}$ $\displaystyle=$
$\displaystyle\cos\theta_{\boldsymbol{\chi}}\hat{\sigma}_{z}+\sin\theta_{\boldsymbol{\chi}}\left(\cos\phi_{\boldsymbol{\chi}}\sigma_{x}+\sin\phi_{\boldsymbol{\chi}}\hat{\sigma}_{y}\right),$
$\displaystyle E_{\boldsymbol{\chi}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\left(h_{z}-\omega\right)^{2}+16\left|G(\boldsymbol{\chi})\right|^{2}},$
$\displaystyle\tan\theta_{\boldsymbol{\chi}}$ $\displaystyle=$
$\displaystyle\frac{2\left|G(\boldsymbol{\chi})\right|}{h_{z}-\omega},$
$\displaystyle\phi_{\boldsymbol{\chi}}$ $\displaystyle=$ $\displaystyle\text{
arg}\,G(\boldsymbol{\chi}).$ (36)

The counting-field dependent Floquet states
$\left|u_{\mu,\boldsymbol{\chi}}\right>$ are the eigenstates of
$\sigma_{\boldsymbol{\chi}}$ with quasienergies $E_{\mu,\boldsymbol{\chi}}=\pm
E_{\boldsymbol{\chi}}$. The dynamical moment-generating function then becomes

$\displaystyle 2M_{dy}(\boldsymbol{\chi},t)$ $\displaystyle=$
$\displaystyle\left[\cos\left(E_{\boldsymbol{\varphi}}t\right)\cos\left(E_{\boldsymbol{\varphi}+\boldsymbol{\chi}}t\right)+\sin\left(E_{\boldsymbol{\varphi}}t\right)\sin\left(E_{\boldsymbol{\varphi}+\boldsymbol{\chi}}t\right)\left<\hat{\sigma}_{\boldsymbol{\varphi}}\hat{\sigma}_{\boldsymbol{\varphi}+\boldsymbol{\chi}}\right>_{0}\right]$
(37) $\displaystyle+$ $\displaystyle
i\left[\cos\left(E_{\boldsymbol{\varphi}}t\right)\sin\left(E_{\boldsymbol{\varphi}+\boldsymbol{\chi}}t\right)\left<\sigma_{\boldsymbol{\varphi}+\boldsymbol{\chi}}\right>_{0}-\sin\left(E_{\boldsymbol{\varphi}}t\right)\cos\left(E_{\boldsymbol{\varphi}+\boldsymbol{\chi}}t\right)\left<\hat{\sigma}_{\boldsymbol{\varphi}}\right>_{0}\right]+(c.c.,\boldsymbol{\chi}\rightarrow-\boldsymbol{\chi})$
$\displaystyle\approx$
$\displaystyle\cos\left[\left(E_{\boldsymbol{\varphi}}-E_{\boldsymbol{\varphi}+\boldsymbol{\chi}}\right)t\right]+i\sin\left[\left(E_{\boldsymbol{\varphi}}-E_{\boldsymbol{\varphi}+\boldsymbol{\chi}}\right)t\right]\left<\hat{\sigma}_{\boldsymbol{\varphi}}\right>_{0}+(c.c.,\boldsymbol{\chi}\rightarrow-\boldsymbol{\chi}),$

where $\boldsymbol{\varphi}=\left(\varphi_{1},\dots,\varphi_{R}\right)$. In
the approximate expression, we have neglected the counting-fields in
$\hat{\sigma}_{\boldsymbol{\varphi}+\boldsymbol{\chi}}\rightarrow\hat{\sigma}_{\boldsymbol{\varphi}}$,
as it only accounts for sub-leading contributions in time.

Benchmarking: Before discussing the implications of the PRFT, we benchmark our
predictions with the exact quantum calculation. For the PRFT, we calculate the
probabilities via Eq. (16), where we obtain the quasiprobabilities by
numerical integration of Eq. (17) using the analytical expression of the
moment-generating function in Eq. (37). To simplify the simulation, we use the
conservation of excitations
$\hat{N}_{ex}=\hat{\sigma}_{z}+\sum_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}$ in
the two-mode Jaynes-Cummings model, and number-squeezed states with a large
mean photon number but a small variance (see Appendix C for more details). We
find that the PRFT perfectly agrees with the exact quantum simulation as can
be seen in Fig. 3.

Probability redistribution: We now investigate the redistribution of photons
between the two modes. Firstly, by examining the probability distribution of
mode $k=1$ [see Figs. 3(b), 3(c), and 3(d)], we observe the Floquet-state
dependence of the photonic dynamics. In Figs. 3(b) [3(d)], we observe that the
probability distribution is shifted to higher (lower) Fock states for the
initial Floquet state $\left|u_{1}\right>$ ( $\left|u_{2}\right>$), i.e.,
photons are flowing into (out of) mode $k=1$. In both cases, the Gaussian
shape of the probability distribution is conserved. This is in agreement with
Eq. (30) stating that the variance does not change for an initial Floquet
state. Finally, in Fig. 3(c), we investigate an initial superposition state of
Floquet states with balanced coefficients $c_{1}=c_{2}=1/\sqrt{2}$. According
to Eq. (31), the time-evolved probability distribution is a weighted product
of the probability to be in a Floquet state $\left|c_{\mu}\right|^{2}$ and the
conditional probabilities $p_{n\mid\mu}$. Consequently, the probability
distribution consists of two peaks, each corresponding to one of the Floquet
states $\left|u_{k}\right>$.

Quasiprobabilities: Analogous to the single-mode Rabi model studied earlier,
we observe that the quasiprobabilities in Figs. 3(e), 3(f), and 3(g) are not a
smooth function of $\Delta n_{1}$. The quasiprobabilities for the Floquet
initial state in Fig. 3(e) are spread in both directions, yet, for positive
$\Delta n_{1}$ the quasiprobabilities alternate between positive and negative
values. In the convolution in Eq. (16), these quasiprobabilities cancel each
other. In contrast, the quasiprobabilities for negative $\Delta n_{1}$ are
essentially positive. These quasiprobabilities add up constructively in the
convolution in Eq. (16), which leads to an overall shift to larger
probabilities in Fig. 3(b).

The quasiprobabilities of the Floquet state $\left|u_{2}\right>$ in Fig. 3(g)
are related to Fig. 3(j) by inversion at $\Delta n_{1}=0$, which can be
explained with a generalized chiral symmetry [108]. Finally, we observe that
the quasiprobabilities of the superposition initial state in Fig. 3(f) are
given by an average of the quasiprobabilities in Figs. 3(e), and 3(g) in
agreement with the explanations in Sec. III.1.

Cumulants: For the two-mode Jaynes-Cummings model, the cumulants for the three
different initial states are shown in Figs. 3(h), 3(i) and 3(j). For the
approximated moment-generating function in Eq. (37), the moments are:

$\displaystyle m_{2l+1}^{(1)}$ $\displaystyle=$
$\displaystyle(-1)^{l}\left(\frac{2g_{1}g_{2}}{E_{\varphi}}\sin(\varphi)t\right)^{2l+1}\left<\hat{\sigma}_{\boldsymbol{\varphi}}\right>_{0}+\mathcal{O}\left(t^{2l}\right),$
$\displaystyle m_{2l}^{(1)}$ $\displaystyle=$
$\displaystyle(-1)^{l}\left(\frac{2g_{1}g_{2}}{E_{\varphi}}\sin(\varphi)t\right)^{2l}+\mathcal{O}\left(t^{2l-1}\right),$
(38)

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
$\varphi=\varphi_{2}-\varphi_{1}$. Interestingly, only the odd moments depend
on the expectation value of $\sigma_{\boldsymbol{\varphi}}$.

When the initial state is a Floquet state, the mean and the variance changes
of the photon number distribution are:

$\displaystyle\left.\Delta\langle\hat{N}_{1}\rangle\right|_{\left|u_{\mu}\right>}$
$\displaystyle=$
$\displaystyle(-1)^{\mu}\frac{2g_{1}g_{2}}{E_{\varphi}}\sin(\varphi)t+\mathcal{O}\left(t^{0}\right),$
$\displaystyle\left.\Delta\sigma_{1}^{2}\right|_{\left|u_{\mu}\right>}$
$\displaystyle=$ $\displaystyle 0+\mathcal{O}\left(t^{0}\right).$ (39)

For $\varphi\neq\left\\{0,\pi\right\\}$, there is a net photon flux between
the photon modes. The variance vanishes up to temporal fluctuations. According
to Eq. (39), the direction and magnitude of the photon flux can be controlled
by the relative phase $\varphi$. 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). Finally,
we note that the higher-order cumulants oscillate very rapidly and do not have
a significant impact on the probability redistribution.

The situation changes dramatically when the initial state is an balanced
superposition of two Floquet states $\left|\phi_{sp}\right>$. In this case the
first two cumulants become

$\displaystyle\left.\Delta\langle
N_{1}\rangle\right|_{\left|\phi_{sp}\right>}$ $\displaystyle=$
$\displaystyle\mathcal{O}\left(t^{0}\right),$
$\displaystyle\left.\Delta\sigma_{1}^{2}\right|_{\left|\phi_{sp}\right>}$
$\displaystyle=$
$\displaystyle\left(\frac{2g_{1}g_{2}}{E_{\varphi}}\sin(\varphi)\right)^{2}t^{2}+\mathcal{O}\left(t^{1}\right),$
(40)

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. (40) is equal to the squared mean change in Eq. (39). This time
dependence is easily understood from Fig. (3)(c), as the rapid variance
increase is a consequence of the linearly growing distance between the two
Floquet-state-dependent peaks.

Figure 4: Purity as an entanglement measure of the light-matter 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
$\sigma^{2}=\text{Var}\;\hat{N}_{1}(0)$, 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 and the parameters
are the same as in Fig. 3.

Light-matter entanglement: To quantify the light-matter entanglement we
investigate the purity:

$\mathcal{P}=\text{Tr}\left[\rho_{M}^{2}\right],$ (41)

where $\rho_{M}$ is the reduced density matrix of the atomic system. The state
is separable (maximally entangled) when $\mathcal{P}=1$ ($\mathcal{P}=1/2$).
We examining the purities for an initial Floquet state and a superposition of
Floquet states employing both exact numerical simulations as well as the PRFT.
The latter predicts

$\displaystyle\mathcal{P}$ $\displaystyle=$
$\displaystyle\frac{1+\left(\left|c_{1}\right|^{2}-\left|c_{2}\right|^{2}\right)^{2}}{2}+2\left|c_{1}c_{2}\right|^{2}e^{-\frac{\left(E^{\prime}_{2,\varphi}-E^{\prime}_{1,\varphi}\right)^{2}t^{2}}{2\sigma^{2}}},$

where, $c_{1}$ and $c_{2}$ are the expansion coefficients of the initial state
in the Floquet basis, and $E^{\prime}_{\mu,\varphi}=dE_{\mu,\varphi}/d\varphi$
(see Appendix D). Our results are shown in Fig. 4.

We find that for an initial Floquet state (e.g., $c_{1}=1$ and $c_{2}=0$), the
purity stays close to $\mathcal{P}\approx 1$ [see Fig. 4(a)] in accoradnce
with Eq. (LABEL:eq:timeDependentPurity). On the other hand, for the balanced
superposition state ($\left|c_{1}\right|^{2}=\left|c_{2}\right|^{2}=1/2$ and
$c_{2}=0$), the purity rapidly decays to $\mathcal{P}\approx 1/2$ [see Fig.
4(b)], showing maximal light-matter entanglement. According to Eq.
(LABEL:eq:timeDependentPurity), the purity degrades faster for a smaller
initial variance $\sigma^{2}$. We recall that the derivation of the PRFT has
assumed a broad initial photon distribution, i.e., large $\sigma$.

###  III.3 Three-mode Rabi model

Figure 5: Analysis of the photonic statistics in the three-mode Rabi model in
Eq. (43). (a) Change of the mean photonic occupation of mode $k=3$. (b) Change
of the corresponding variance. Parameters are $h_{z}=2.1\omega_{1}$,
$\omega_{2}=2\omega_{1}$,$\omega_{3}=3\omega_{1}$, and
$g_{k}=\alpha_{k}\tilde{g}_{k}=\omega_{1}$,
$\varphi_{1}=\phi_{3}=2\varphi_{2}=0.5$.

To illustrate that the light-matter entanglement explained in Sec. III.1 is a
generic effect, we now investigate a three-mode Rabi model

$H_{3MR}=\frac{h_{z}}{2}\hat{\sigma}_{z}+\sum_{k=1}^{3}\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k}+\sum_{k=1}^{3}\tilde{g}_{k}\hat{\sigma}_{x}\left(\hat{a}_{k}+\hat{a}_{k}^{\dagger}\right),$
(43)

with distinct commensurate photonic frequencies $\omega_{k}$. For
concreteness, we choose $\omega_{2}=2\omega_{1}$ and $\omega_{3}=3\omega_{1}$.
In this case, the semiclassical Hamiltonian is a periodically-driven Rabi
model with period $2\pi/\omega_{1}$, and we compute the photon-mean and
photon-variance changes using the PRFT.

Following the analysis for the two-mode Jaynes-Cummings model, we consider the
three initial conditions $\left|\phi(0)\right>=\left|u_{1}\right>$,
$\left|\phi(0)\right>=\left|u_{2}\right>$, and
$\left|\phi(0)\right>=\frac{1}{\sqrt{2}}\left(\left|u_{1}\right>+\left|u_{2}\right>\right)$,
where the $\left|u_{\mu}\right>$ are the Floquet states of the matter system,
i.e., the eigenstates of Eq. (26). As shown in Fig. 5, we find that the
photonic mean and variance change for the three-mode Rabi model is
qualitatively similar to that of the two-mode Jaynes-Cummings model [cf.
corresponding curves in Fig. 3 (b), 3(c) and 3(d)]. In particular, when the
initial state is the Floquet state $\left|u_{1}\right>$
($\left|u_{2}\right>$), the mean photon number grows (decreases) linearly in
time, while the variance remains almost unchanged. In contrast, the mean
photon change is close to zero for the superposition state, while the variance
increases rapidly. Thus analogous to the two-mode Jaynes Cummings model, this
effect can be interpreted as the generation of a maximally entangled light-
matter state.

##  IV Experimental and technological implications

In this section, we discuss the decoherence related to the light-matter
entanglement introduced in Sec. III.1 for coherently driven system. The effect
can be detected in experiments and impacts the design of quantum memories and
quantum operations.

###  IV.1 Quantum-optical coherence time

To estimate the quantum-optical coherence time for periodically-driven
systems, we consider the mean photon-number changes
$\left.\Delta\langle\hat{N}_{k}\rangle\right|_{\left|u_{\mu}\right>}=-E^{\prime}_{\mu}(\varphi_{k})t$
for two distinct Floquet states $\mu_{1}$ and $\mu_{2}$, where
$E^{\prime}_{\mu}(\varphi)=dE_{\mu}(\boldsymbol{\varphi})/d\varphi$. 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.,

$t_{c}\cdot\left|E^{\prime}_{\mu_{1}}(\varphi_{k})-E^{\prime}_{\mu_{2}}(\varphi_{k})\right|=\sigma_{k},$
(44)

where $\sigma_{k}^{2}$ 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. III.1. To estimate $\sigma_{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. Either 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 $\overline{n}_{k}=\sigma_{k}^{2}$. The mean
photon number in a cavity mode is given by
$\overline{n}_{k}=\epsilon_{0}E^{2}\cdot V/(2\hbar\omega_{k})$ where $E$ is
the electromagnetic field, $V$ is the cavity volume, and $\epsilon_{0}$ is the
dielectric constant. Putting everything together, the quantum-optical
coherence time is given as

$t_{c}=\min_{k,\mu_{1},\mu_{2}}\frac{\sqrt{\frac{\epsilon_{0}E^{2}\cdot
V}{2\hbar\omega_{k}}}}{\left|E^{\prime}_{\mu_{1}}(\varphi_{j})-E^{\prime}_{\mu_{2}}(\varphi_{j})\right|}.$
(45)

We note that the coherence time in Eq. (45) 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 $\omega_{k}$ and spectral width $\Delta\omega\propto 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
$\overline{n}_{k}=P(\omega_{k})\cdot t_{p}$, where $P(\omega_{k})$ is the
power of the electromagnetic field at frequency $\omega_{k}$. Now, we model
the pulse as a cavity mode with initial occupation $\overline{n}_{k}$.
Consequently, we find that the quantum-optical coherence time is:

$t_{c}=\min_{k,\mu,\mu^{\prime}}\frac{P(\omega_{k})}{\hbar\omega_{k}\left|E^{\prime}_{\mu_{1}}(\varphi_{k})-E^{\prime}_{\mu_{2}}(\varphi_{k})\right|^{2}}.$
(46)

We note that since each photonic mode contributes to the decoherence effect,
the coherence time in Eq. (46) can be considered as an upper bound.

###  IV.2 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 Jaynes-Cumming 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. To this end, we use that the
purity of the system decreases more rapidly for driving fields with smaller
intensity fluctuations, i.e, smaller initial photon number variance [see Eq.
(45)]. For the two-mode Jaynes-Cummings model, this effect can be seen in Fig.
(4)(b). To isolate the quantum-optical decoherence, we vary the amplitudes of
the coherent states $\alpha_{k}$ and the light-matter interactions
$\tilde{g}_{k}$ such that the effective parameters
$g_{k}=\tilde{g}_{k}\alpha_{k}$ remain constant. In doing so, the spin system
experience the same semiclassical driving fields [see Eq. (III.2)] and, thus,
is subject to the same decoherence sources other than the quantum-optical
decoherence. As for a coherent driving field
$\sigma_{k}\propto\alpha_{k}\propto P(\omega_{k})^{1/2}$, the quantum-optical
decoherence dynamics can be accessed by measuring the purity decay for various
driving-field power $P(\omega_{k})$.

To get an estimate for the coherence time in the optical regime, we evaluate
Eq. (46) for the two-mode Rabi frequency, where the photonic frequency is
assumed to be $\omega=400\,\text{THz}$. The quasienergy difference in the
denominator in Eq. (46) can be approximated by a typical Rabi frequency of
atomic systems $\Omega_{R}=40\,\text{MHz}$, corresponding to a laser power
$P(\omega)=10\,\mu\text{W}$. For these parameters, the quantum-optical
coherence time is $t_{c}=5\,\text{ms}$, which is comparable to the duration of
typical cold-atom experiments [109, 110, 111], but significantly shorter than
quantum information storage times achieved with trapped ions [112].

We obtain the order of magnitude of the coherence time for radio frequencies
using the following parameters: $\omega=10\,\text{MHz}$,
$\Omega_{R}=10\,\text{kHz}$, and $P(\omega)=10\,\text{W}$, which are typical
in current experiments [113]. In this regime, the quantum-optical coherence
time is very long ($t_{c}\approx 3\cdot 10^{18}\,\text{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.

###  IV.3 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 [113], and
more than $50s$ in trapped ions [112]. In both cases, the quantum information
is stored in the hyper-fine levels of the ground-state manifold, which are
energetically separated in the radio-frequency regime.

Depending on the control protocol, the quantum-optical coherence time in Eq.
(45) influences the performance of quantum memories. The quantitative
considerations in Sec. IV.2 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 [114]. Even though the pulse duration in
these cases may be rather short, $t_{p}<1\;\mu\text{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 requires high fidelities $>99\%$ [115],
even a minor quantum-optical decoherence can have a significant effect on
quantum operations. Furthermore, inspection of the coherence time in Eq. (46)
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.

###  IV.4 Quantum time crystals

Recently, quantum time crystals have gained significant attention as a
potential resource for quantum memories [116, 117, 118]. Discrete time
crystals, that are driven by an external driving field with period
$\tau=2\pi/\omega$, exhibits subharmonic response with a frequency $\omega/n$,
where $n$ 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 $\omega/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 Jaynes-Cummings model and post-select the
measurement results. We now proceed to explain how Greenberger-Horne-Zeilinger
(GHZ) states can be employed to speed-up the light-matter entanglement
generation, before delving into the details of the communication protocol in
Sec. V.2.

###  V.1 Speeding up the light-matter entanglement

Typically, the light-matter interaction between a single atom and the light-
field is relative weak, which slows down the generation of maximally-entangled
light-matter states as depicted in Fig. 3(c). To enhance the effect, we again
harness the setup in Fig. 3(a), but with $N_{A}$ atoms instead of a single
one. The corresponding Hamiltonian reads as

$H=\sum_{j}\frac{h_{z}}{2}\hat{\sigma}_{z}^{(j)}+\tilde{g}\sum_{j,k}\hat{\sigma}_{+}^{(j)}\hat{a}_{k}+\hat{\sigma}_{-}^{(j)}\hat{a}_{k}^{\dagger}+\sum_{k}\omega\hat{a}_{k}^{\dagger}\hat{a}_{k},$
(47)

where the Pauli operators
$\hat{\sigma}_{z}^{(j)},\hat{\sigma}_{+}^{(j)},\hat{\sigma}_{-}^{(j)}$ act on
the atoms $j=1,\dots,N_{A}$. As in Sec. III.2, we consider two photonic modes
$k=1,2$ which are initially in coherent states
$\left|\alpha_{k}e^{i\varphi_{k}}\right>$. This system is a two-mode Tavis-
Cummings model in which the many-body interaction is mediated via the
quantized electromagnetic field.

Since this model is difficult to analyze for large atom numbers, we proceed
with a semiclassical treatment. Upon introducing the counting fields
$\chi_{1}$ and $\chi_{2}$, the corresponding semiclassical Hamiltonian reads
as

$\displaystyle\mathcal{H}(t)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N_{A}}\mathcal{H}^{(j)}(t),$
$\displaystyle\mathcal{H}^{(j)}(t)$ $\displaystyle=$
$\displaystyle\frac{h_{z}}{2}\sigma_{z}^{(j)}+\sum_{k=1}^{2}\sigma_{+}^{(j)}g_{k}e^{-i\omega
t-i\chi_{k}}+\text{h.c.},$ (48)

where $\mathcal{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. Each
$\mathcal{H}^{(j)}(t)$ is identical to the Jaynes-Cummings model in Eq.
(III.2), and the quasienergies of the total system can be written as the sum
of the quasienergies of the individual atoms

$E_{\boldsymbol{\mu},\boldsymbol{\chi}}^{(N_{A})}=\sum_{j=1}^{N_{A}}E_{\mu_{j},\boldsymbol{\chi}},$
(49)

where $E_{\mu_{j},\boldsymbol{\chi}}$ is the quasienergy of atom $j$ in
Floquet state $\mu_{j}$. We have introduced the vector notation
$\boldsymbol{\mu}=\left(\mu_{1},\dots,\mu_{N_{A}}\right)$, that contain the
quantum numbers $\mu_{j}$ of the $N_{A}$ atoms. The corresponding Floquet
states read as

$\left|\Psi_{\boldsymbol{\mu},\boldsymbol{\chi}}\right>=\bigotimes_{j}\left|u^{(j)}_{\mu_{j},\boldsymbol{\chi}}\right>,$
(50)

where the $\left|u^{(j)}_{\mu_{j},\boldsymbol{\chi}}\right>$ denote the
Floquet state of atom $j$ with quantum number $\mu_{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
$\boldsymbol{\mu}^{(1)}=(1,\dots,1)$ or $\boldsymbol{\mu}^{(2)}=(2,\dots,2)$,
i.e., the state in which all atoms $j$ are in the same Floquet state
$\mu_{j}=\alpha\in\left\\{1,2\right\\}$. The quasienergy of the atom ensemble
in either of these states is then

$E_{\boldsymbol{\mu}^{(\alpha)},\boldsymbol{\chi}}=N_{A}E_{\alpha,\boldsymbol{\chi}}.$
(51)

According to Sec. III.1, this leads to an enhancement of the photon transport
proportional to the atom number. To be more specific, we choose $h_{z}=\omega$
and the driving phases $\varphi_{1}=\pi/4,\varphi_{2}=-\pi/4$. In this case,
the Floquet states of the matter system $\left|u_{\mu}\right>$ for the Jaynes-
Cummings model are
$\left|u_{1}\right>=\left|-\right>=\left(\left|0\right>-\left|1\right>\right)/\sqrt{2}$
and
$\left|u_{2}\right>=\left|+\right>=\left(\left|0\right>+\left|1\right>\right)/\sqrt{2}$,
and the Floquet states of the atom ensemble in Eq. (50) are
$\left|\Psi_{\boldsymbol{\mu}^{(1)}}\right>=\left|-\cdots-\right>\equiv\left|\boldsymbol{-}\right>$
and
$\left|\Psi_{\boldsymbol{\mu}^{(2)}}\right>=\left|+\cdots+\right>\equiv\left|\boldsymbol{+}\right>$.

By constraining the system to the collective states
$\left|\boldsymbol{-}\right>$ and $\left|\boldsymbol{+}\right>$, we restrict
the model to an effective two-level system with a renormalized quasienergy,
for which the findings in the Jaynes-Cummings model in Sec. III.2 are valid.
To create light-matter entanglement in the basis $\left|\boldsymbol{+}\right>$
and $\left|\boldsymbol{-}\right>$, we use the initial condition

$\left|\Psi(0)\right>=\hat{U}_{H}^{N_{A}}\left|\text{GHZ}\right>,$ (52)

where $\left|\text{GHZ}\right>$ is the celebrated GHZ state. This state is
defined by a superposition
$\left|\text{GHZ}\right>=\left(\left|\boldsymbol{0}\right>+\left|\boldsymbol{1}\right>\right)/\sqrt{2}$,
where $\left|\boldsymbol{0}\right>\equiv\left|0\dots 0\right>$ and
$\left|\boldsymbol{1}\right>\equiv\left|1\dots 1\right>$. The Hadamard gate
$\hat{U}_{H}=\exp\left[-i\pi\hat{\sigma}_{y}/4\right]$, 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 diminished. We emphasize that even though we take the Jaynes-
Cummings model as an example, this enhancement effect is valid for general
Floquet systems according to the PRFT.

Figure 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.2. 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 light-mater 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 $p_{n}$ 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.

###  V.2 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. Quantum state transfer can then be carried out via quantum teleportation
[119]. 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. 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
$\left|\Psi_{AB}(0)\right>=\left|\text{GHZ}_{A}\right>\otimes\left|\text{GHZ}_{B}\right>$
is separable. Alice and Bob carry out local Hadamard gates such that their
states become

$U_{H,X}^{N_{A}}\left|\text{GHZ}\right>_{X}=\left|\boldsymbol{-}\right>_{X}+\left|\boldsymbol{+}\right>_{X},$
(53)

where $X=A,B$. The choice of the basis $\left|\boldsymbol{-}\right>$ and
$\left|\boldsymbol{+}\right>$ is thereby determined by the driving phases,
which we assume to be $\varphi_{1}=-\pi/4$ and $\varphi_{2}=\pi/4$.

  2. 2.

Light-matter interaction: Alice impinges two coherent light beams onto her
atom ensemble leading to the generation of a light-matter entangled state:

$\left(\left|\boldsymbol{-}\right>_{A}\left|A_{-}\right>+\left|\boldsymbol{+}\right>_{A}\left|A_{+}\right>\right)\left(\left|\boldsymbol{-}\right>_{B}+\left|\boldsymbol{+}\right>_{B}\right),$
(54)

where $\left|A_{+}\right>$ ($\left|A_{-}\right>$) 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

$\displaystyle\left|\boldsymbol{-}\right>_{A}\left|\boldsymbol{-}\right>_{B}\left|A_{2-}\right>+\left|\boldsymbol{+}\right>_{A}\left|\boldsymbol{-}\right>_{B}\left|A_{0}\right>$
$\displaystyle+\left|\boldsymbol{-}\right>_{A}\left|\boldsymbol{+}\right>_{B}\left|A_{0}\right>+\left|\boldsymbol{+}\right>_{A}\left|\boldsymbol{+}\right>_{B}\left|A_{2+}\right>,$

where the photonic states $\left|A_{2-}\right>$, $\left|A_{0}\right>$,
$\left|A_{2+}\right>$ 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 states of the protocol is sketched in Fig. 6(b).

  3. 3.

Measurement and postselection: Bob makes a projective measurement defined by
the operator $\left|A_{0}\right>\left<A_{0}\right|$. 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

$\left|\boldsymbol{0}\right>_{A}\left|\boldsymbol{1}\right>_{B}+\left|\boldsymbol{1}\right>_{B}\left|\boldsymbol{0}\right>_{A}$
(55)

in the basis of collective excitations $\left|\boldsymbol{0}\right>_{X}$ and
$\left|\boldsymbol{1}\right>_{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 $\left|A_{2+}\right>$. Consequently, Alice’s and Bob’s ensemble
both are in the state $\left|\boldsymbol{+}\right>$ and thus not entangled. A
similar reasoning applies to the state $\left|A_{2-}\right>$. 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\%$.

###  V.3 Quantum state transfer rate

The major obstacle in quantum communication is photon loss. A typical damping
rate of optical fibers is $\gamma\approx 0.051\,1/\text{km}$, leading to a
loss of more than $99\%$ of photons after $100\,\text{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\,\text{Hz}$ over $500\,\text{km}$ [96]. These protocols typically employ
quantum repeaters strategically placed between the transmission endpoints.

The coherent-light protocol introduced in Sec. V.2 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., $\langle\Delta\hat{N}_{1}\rangle\approx\Delta\sigma_{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 $\gamma_{Amp}=\gamma$, such the
mean-photon number is conservered 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 $\gamma$,
the width of each peak increases by $\Delta\sigma_{1}=\sqrt{2\gamma d\cdot
Pt_{p}/(\hbar\omega)}$, where $P$ is the power of the transmitted pulse,
$t_{p}$ is the pulse duration, and $d$ the transfer distance. We recall from
Sec. III.2 that the distance between the two peaks is given as
$\langle\Delta\hat{N}_{1}\rangle=|E_{1,\varphi}^{\prime}-E_{0,\varphi}^{\prime}|\cdot
t$. In atomic systems, the quasienergy splitting can be associated with the
Rabi frequency $\Omega_{R}$. Using the GHZ amplification in Eq. (51), the
probability peaks are separated by
$\langle\Delta\hat{N}_{1}\rangle=N_{A}\Omega_{R}\cdot t_{p}$. The quantum
state transfer rate can be thus estimated as

$f=1/t_{p}=\frac{N_{A}^{2}\Omega_{R}^{2}\hbar\omega}{2\gamma P}\frac{1}{d}.$
(56)

For instance, we assume a photonic frequency $\omega=400\,\text{THz}$ and a
typical Rabi frequency of atomic systems $\Omega_{R}=40\,\text{MHz}$,
corresponding to a laser power $P=10\,\mu\text{W}$. Recent experiments of the
Lukin group have successfully created GHZ states with $N_{A}=12$ Rydberg atoms
[120]. For these numbers, the transfer rate is $122\text{Hz}$ over
$500\text{km}$, thus exceeding typical few photon protocol by two orders of
magnitude.

###  V.4 Discussion

Even though strongly idealized, the proposed protocol merits a thorough
discussion. As quantum information is typically 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 Jaynes-Cummings model used
here for illustration, as the light-matter entanglement is a generic effect
appearing in all Flouqet 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. (56) scales quadratically
in the number of atoms $N_{A}$. It has been estimated that for fault-tolerant
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.3.

  * •

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 $\alpha$, the transfer rate can be
increased by that factor $\alpha$. Bob must amplify the received signal before
interaction with his GHZ state. Noiseless de-amplification and amplification
can be implemented with additional atom ensembles which are prepared in
Floquet states as discussed in Sec. III.2 and shown in Figs. 3(b) and 3(d).
This might enhance the transfer rate by an additional factor of $\alpha=10$.

Taking these points into account, the transfer rate will be on the order of
$f\approx 10\text{MHz}$, 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 [96].

  * •

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: At current, there do not exist on-demand single-photon
sources. 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 most
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 [121], other quantum error correction codes are based on graph states,
that are generalizations of the GHZ states [115]. 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 $\hat{a}_{k}^{\dagger}\hat{a}_{k}$ for
$k=1,2$. Fortunately, as the ratio of $\langle\Delta\hat{N}_{k}\rangle$ and
$\langle\hat{N}_{k}\rangle$ 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
$\left|\boldsymbol{-}\right>_{X}\leftrightarrow\left|\boldsymbol{+}\right>_{X}$
for $X=A,B$, leading to a periodic variation of
$\langle\Delta\hat{N}_{k}\rangle$.

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 [115]. Assuming a pulse length of
$t_{p}=20\,\text{ms}$, the two peaks are separated by
$\langle\Delta\hat{N}_{1}\rangle\approx 10^{7}$ photons which corresponds to a
Hamming distance of $\log_{2}\langle\Delta\hat{N}_{1}\rangle\approx 24$ bits.
As the total number of photons is on the order of
$\langle\hat{N}_{1}\rangle=Pt/\hbar\omega\approx 10^{12}$ corresponding to
$40$ qubits, the physical encoding presented here is thus comparable to a
$\left[\left[40,1,24\right]\right]$ quantum error correction code.

##  VI Summary and outlook

In this paper, we have taken a significant foundational step forward in the
development of Floquet theory by conceiving a framework — the “Photon-resolved
Floquet theory” (PRFT) — which unifies quantum optics and Floquet theory via
FCS. In this section, we first summarize the salient characteristics of the
PRFT and then discuss some important predictions that emerge from this
framework as well as an outlook for future work.

###  VI.1 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 is based on semiclassical equations and thus shares
the same computational complexity as Floquet theory. The treatment of the
matter subsystem is equivalent in both frameworks. However, 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. We
emphasize that the PRFT has significant computational advantages over Sambe
space methods that investigate photonic dynamics by effectively re-quantizing
the semiclassical driving field [93, 94, 95]. Moreover, unlike other
approaches, the PRFT can distinguish between modes with commensurate
frequencies, thereby extending its reach over a wide class of driven systems.

Quantum optics: An extremely appealing feature of the PRFT is that it makes
exact 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 [54, 122, 123], 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 multi-mode driven systems, which is
essentially infeasible with standard quantum optical methods. Furthermore, the
PRFT is non-perturbative in nature, and thus capable of making predictions
beyond methods from non-linear spectroscopy [124].

Full-counting statistics: Analogous to the standard FCS, the PRFT provides a
stochastic description of photonic fields in terms of probabilities, moments
and cumulants [68, 125]. However, since the standard FCS is based on two-point
projective measurements, which allows to correctly describe spontaneous photon
emission [65, 66, 67] , 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.

###  VI.2 Discussion and Outlook

The PRFT has far-reaching implications for quantum science and technologies.
In particular, the PRFT predicts that multi-mode driven systems undergo
complete decoherence at long times, with the Floquet states acting as the
decohering basis. This quantum-optical decoherence has serious consequences
for quantum memories and quantum operations. As quantum memories usually use
sophisticated control protocols based on periodic driving, 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 ($\approx$ 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 [126], heat engines [87] and interacting spin chains in a
cavity [127, 128, 129]. Finally, we note that the PRFT can also find
applications in areas unrelated to light-matter interactions: the moment-
generating function in Eq. (10) is formally equivalent to the position
operator studied in the context of macroscopic polarization of crystal
structures [130] and topological insulators [131, 132, 11]. A modified version
of the PRFT can thus make predictions about the electronic dipole fluctuations
in semiconductors.

## Acknowledgments

This work is supported by the Guangdong Provincial Key Laboratory (Grant
No.2019B121203002), the AFOSR Grant No. FA9550-16-1-0006, the MURI-ARO Grant
No. W911NF17-1-0323 through UC Santa Barbara, and the Shanghai Municipal
Science and Technology Major Project through the Shanghai Research Center for
Quantum Sciences (Grant No. 2019SHZDZX01). The authors acknowledge inspiring
discussions with Victor Manuel Bastidas, Georgious Sivilogou, and Zhigang Wu.

##  Appendix A Photon-resolved Floquet theory: Derivations

In this Appendix, we study the PRFT in more details. In the following
derivations, we focus on the special case of a single counting field to
enhance the readability. The generalization to multi-mode photon fields works
along the same lines.

###  A.1 Unraveling of the quantum dynamics

To reveal the quantum properties of the Floquet theory, we start from the
underlying quantum Hamiltonian in Eq. (1) with a single photonic mode $k=1$.
We unravel the corresponding time-evolution operator in a perturbation picture
defined by $H_{S}=H_{0}+\omega\hat{a}^{\dagger}\hat{a}$. In doing so, we find

$\displaystyle U(t)$ $\displaystyle=$ $\displaystyle e^{-iH_{Q}t}$ (57)
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\cdots\int_{0}^{t}dt_{n}\int_{0}^{t_{2}}dt_{1}U^{(n)}(t),$

where we have defined the perturbation-resolved propagation operator

$\displaystyle U^{(n)}(t)=e^{-iH_{S}(t-t_{n})}\cdots
e^{-iH_{S}(t_{2}-t_{1})}H_{1}\left(\hat{a}^{\dagger}+\hat{a}\right)e^{-iH_{S}t_{1}},$

in which the term $H_{1}\left(\hat{a}^{\dagger}+\hat{a}\right)$ appears $n$
times. Next, we introduce the counting field $\chi$ into the expansion in Eq.
(57) by replacing $\hat{a}\rightarrow\hat{a}e^{-i\chi}$ and
$\hat{a}^{\dagger}\rightarrow\hat{a}^{\dagger}e^{i\chi}$. This generalizes the
time-evolution operator $U(t)\rightarrow U_{\chi}(t)$, which can be formally
written as

$U_{\chi}(t)=e^{i\chi\hat{N}}U(t)e^{-i\chi\hat{N}}.$ (59)

We can expand the generalized time-evolution operator as

$\displaystyle U_{\chi}(t)$ $\displaystyle=$
$\displaystyle\sum_{m=\infty}^{\infty}U^{(m)}(t)e^{i\chi m},$ (60)

where the photon-resolved time-evolution operators $U^{(m)}$ contain all terms
in which the difference of the number of creators $m_{C}$ and the number of
annihilators $m_{A}$ is $m=m_{C}-m_{A}$. Each $U^{(m)}(t)$ can be thus
represented as a polynomial of $\hat{a}$ and $\hat{a}^{\dagger}$. The non-
Hermitian $U^{(m)}$ can be obtained from $U_{\chi}(t)$ by performing a Fourier
transformation with respect to $\chi$, i.e.,

$U^{(m)}(t)=\frac{1}{2\pi}\int_{0}^{2\pi}d\chi U_{\chi}(t)e^{-im\chi}.$ (61)

Using these photon-resolved time-evolution operators, we can express
expectation values of arbitrary observables $\hat{O}$ (can be either matter-
like, photonic, or mixed observables) as

$\left<\psi(t)\right|\hat{O}\left|\psi(t)\right>=\sum_{m,m^{\prime}}\left<\psi(t_{0})\right|U^{(m)\dagger}(t)\hat{O}U^{(m^{\prime})}(t)\left|\psi(t_{0})\right>.$
(62)

We assume a product state as initial state

$\displaystyle\left|\psi(t_{0})\right>=\left|\phi(t_{0})\right>\otimes\sum_{n=1}^{\infty}a_{n}\left|n\right>,$
(63)

where the expansion coefficients $a_{n}$ shall fulfill the conditions (i) and
(ii) introduced in Sec. II.1. To illustrate the application of the PRFT
presented in Sec. A.2 below, we consider the projector onto a particular Fock
state

$\hat{P}_{n}=\left|n\right>\left<n\right|,$ (64)

but more complicated operators can be treated accordingly. For this specific
operator, the expectation value can be written in terms of the photon-resolved
time-evolution operators as

$\displaystyle\langle\hat{P}_{n}\rangle$
$\displaystyle\equiv\left<\psi(t)\right|\hat{P}_{n}\left|\psi(t)\right>$ (65)
$\displaystyle=\sum_{m,n}\left<\left[U^{(n-m)\dagger}(t)\right]_{m,n}\left[U^{(n-m)}(t)\right]_{m,n}\right>_{t_{0}}\left|a_{m}\right|^{2},$

where we have defined the reduced photon-resolved propagation matrices by

$\left[U^{(m)\dagger}(t)\right]_{n_{1},n_{2}}\equiv\left<n_{1}\right|U^{(m)\dagger}\left|n_{2}\right>,$
(66)

which acts on the states in the matter system. The expectation value in Eq.
(65) is taken with respect to the matter initial state
$\left|\phi(t_{0})\right>$.

###  A.2 Transition to the Photon-resolved Floquet theory

The derivations in Appendix A.1 have been carried out in the Fock space. The
connection of the photon-resolved quantum time-evolution and the PRFT is
established via Eq. (60). We recall that $U_{\chi}(t)$ can be represented as a
polynomial of $\hat{a}^{\dagger}$ and $\hat{a}$ operators. As in the standard
semiclassical approximation, we now replace

$e^{i\omega\hat{a}^{\dagger}\hat{a}t}\hat{a}e^{-i\omega\hat{a}^{\dagger}\hat{a}t}\rightarrow\alpha
e^{-i\omega t-i\varphi},$ (67)

where $\alpha\gg 1$ and $\varphi$ are the amplitude and the phase of a
coherent state. More generally, this replacement is also valid for a state
which fulfills the conditions (i) and (ii) in Sec. II.1. We denote the
resulting operator as $\mathcal{U}_{\chi}(t)$, which acts on the matter
system. The transition to the PRFT is readily done by realizing that

$\mathcal{U}_{\chi}(t)\stackrel{{\scriptstyle!}}{{=}}\mathcal{T}e^{-i\int_{0}^{t}\mathcal{H}_{\chi}(t^{\prime})dt^{\prime}},$
(68)

where $\mathcal{T}$ is the time-ordering operator and the generalized time-
periodic Hamiltonian on the r.h.s. is defined as

$\mathcal{H}_{\chi}(t)=H_{0}+H_{1}\alpha\left(e^{i\omega t+i\chi}+e^{-i\omega
t-i\chi}\right).$ (69)

The calculation of $\mathcal{U}_{\chi}(t)$ can be thus performed by
analytically or numerically solving the time-dependent Schrödinger equation
with $\mathcal{H}_{\chi}(t)$ for all $\chi\in\left[0,2\pi\right)$. The final
step in the transition is done by evaluating

$\mathcal{U}^{(m)}(t)=\frac{1}{2\pi}\int_{0}^{2\pi}d\chi\mathcal{U}_{\chi}(t)e^{-im\chi}$
(70)

and replacing

$\left[U^{(m)\dagger}(t)\right]_{n_{1},n_{2}}\rightarrow\mathcal{U}^{(m)}(t)$
(71)

in Eq. (65). This replacement is well justified as the matrix elements of the
terms $\hat{a}^{{}^{\dagger}n}\hat{a}^{m}$ depend only weakly on the order of
the operators for large photon numbers. The expectation value of $\hat{P}_{n}$
in the PRFT thus reads as

$\left<\hat{P}_{n}\right>=\sum_{m_{1},m_{2}}\left<\mathcal{U}^{(m_{2})\dagger}(t)\mathcal{U}^{(m_{1})}(t)\right>_{t_{0}}a_{n-m_{1}}^{*}a_{n-m_{2}},$
(72)

where the expectation value is taken in the matter initial state.

###  A.3 Full-Counting Statistics

In the following, we derive the dynamical cumulant-generating function of the
photon field given in Eq. (13). 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

$\displaystyle K_{\chi}(t)$ $\displaystyle=$
$\displaystyle\log\left[M_{\chi}(t)\right],$ $\displaystyle M_{\chi}(t)$
$\displaystyle=$ $\displaystyle\left<e^{i\chi\hat{N}}\right>_{t},$ (73)

where the time evolution is calculated with the full quantum Hamiltonian in
Eq. (3). The moment-generating function in Eq. (73) can be rewritten as

$\displaystyle M_{\chi}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left<U^{\dagger}(t)U_{\chi}(t)e^{i\chi\hat{N}}+e^{i\chi\hat{N}}U_{-\chi}^{\dagger}(t)U(t)\right>_{t_{0}},$

where $U_{\chi}(t)$ is the generalized time-evolution operator in Eq. (59).
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_{\chi}^{*}(t)=M_{-\chi}(t)$ is maintained in the semiclassical limit, which
guarantees that all moments and cumulants are real valued.

Using now the photon-resolved time-evolution operator, we can expand Eq.
(LABEL:eq:tra:symmetrization) in terms of photon processes such that

$\displaystyle M_{\chi}(t)$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\sum_{n,m_{1},m_{2}}\left<\phi(t_{0})\right|\mathcal{U}^{(m_{1})\dagger}(t)e^{i\chi
m_{2}}\mathcal{U}^{(m_{2})}(t)e^{i\chi\left(n-m_{2}\right)}\left|\phi(t_{0})\right>\times
a_{n-m_{1}}^{*}a_{n-m_{2}}+\left(c.c.,\chi\rightarrow-\chi\right)$ (75)
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{n,m_{1},m_{2}}\left<\phi(t_{0})\right|\mathcal{U}^{(m_{1})\dagger}(t)\mathcal{U}^{(m_{2})}(t)e^{i\chi
n}\left|\phi(t_{0})\right>\times
a_{n-m_{1}+m_{2}}^{*}a_{n}+\left(c.c.,\chi\rightarrow-\chi\right).$

In this expression, we have already carried out the semiclassical replacement
of the photon-resolved time-evolution operators introduced in Eq. (71). Now we
use the assumptions in Sec. II.1 that the modulus of the expansion
coefficients $a_{n}$ varies only very slow as a function of photon number, and
their phases behave as $\arg a_{n}=n\cdot\varphi$. Therefore, for large $n$,
we can safely replace

$\displaystyle a_{n-m_{2}+m_{1}}\rightarrow
a_{n}e^{i\left(n-m_{2}+m_{1}\right)\varphi}.$ (76)

This replacement has been benchmarked in Secs. II.2 and III.2, where we have
compared numerical simulations of the quantum Rabi model and the multi-mode
Jaynes-Cummings model with the predictions of the PRFT. Some explanations
about the relation $M_{\chi}(t)=M_{-\chi}^{*}(t)$ are in order. For the
replacement in Eq. (76), this relation is not anymore microscopically
fulfilled. However, the relation becomes valid again after taking the
semiclassical limit in the following derivations.

Next, we insert the definition of the photon-resolved time-evolution operator
in Eq. (61), such that

$\displaystyle M_{\chi}(t)$ $\displaystyle=$
$\displaystyle\iint_{0}^{2\pi}d\frac{\chi_{1}d\chi_{2}}{\left(2\pi\right)^{2}}\sum_{m_{1},m_{2}}\left<e^{-im_{1}\left(\chi_{1}+\varphi\right)}\mathcal{U}_{-\chi_{1}}^{\dagger}(t)\mathcal{U}_{\chi_{2}}(t)e^{im_{2}\left(\chi-\chi_{2}+\varphi\right)}\right>_{t_{0}}\times\sum_{n}e^{i\chi
n}a_{n}a_{n}+\left(c.c.,\chi\rightarrow-\chi\right).$

We now use the identity $\frac{1}{2\pi}\sum_{m}e^{i\chi m}=\delta(x)$ to
simplify this expression. Identifying the second factor in Eq.
(LABEL:eq:der:momGenFkt) as $M_{\chi}(0)$, we finally obtain

$\displaystyle M_{\chi}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left<\mathcal{U}_{\varphi}^{\dagger}(t)\mathcal{U}_{\varphi+\chi}(t)+\mathcal{U}_{\varphi-\chi}^{\dagger}(t)\mathcal{U}_{\varphi}(t)\right>_{t_{0}}\cdot
M_{\chi}(t_{0}).$

In this form, the semiclassical moment-generating function fulfills the
correct transformation properties under inversion of the counting field
$\chi\rightarrow-\chi$. A non-symmetric representation, e.g., when only taking
the first summand in Eq. (LABEL:eq:res:momenGenFunction) does violate this
basic property. Defining the dynamical moment-generation function as

$M_{dy,\chi}(t)\equiv\frac{1}{2}\left<\mathcal{U}_{\varphi}^{\dagger}(t)\mathcal{U}_{\varphi+\chi}(t)+\mathcal{U}_{\varphi-\chi}^{\dagger}(t)\mathcal{U}_{\varphi}(t)\right>_{t_{0}},$
(79)

using the relation $K_{\chi}(t)=\log M_{\chi}(t)$, and recalling the
definition of the dynamical cumulant-generation function in Eq. (12),

$\displaystyle K_{\chi}(t)$ $\displaystyle\equiv$ $\displaystyle
K_{dy}(\chi,t)+K_{\chi}(t_{0}),$ (80)

we find that $K_{dy}(\chi,t)=\log M_{dy,\chi}(t)$, which is the expression
given in Eq. (13).

###  A.4 Relation to classical driving

Here, we show that the cumulant-generating function in Eq. (12) reproduces the
standard semiclassical definition of the energy current as considered in,
e.g., Refs. [93, 94]. We first consider the cases of a single driving mode.
Using the definition of the first dynamical cumulant in Eq. (11), we find

$\displaystyle\kappa_{dy,1}(t)$ $\displaystyle=$
$\displaystyle\Delta\left<\hat{N}(t)\right>_{t_{0}}$ $\displaystyle=$
$\displaystyle\frac{-i}{2}\left.\left<\mathcal{U}_{\varphi}(t)\frac{d}{d\chi}\mathcal{U}_{\varphi+\chi}(t)\right>_{t_{0}}\right|_{\chi\rightarrow
0}+\textit{ c.c.}$

To evaluate the derivative, we use the identity

$\displaystyle\mathcal{U}_{\chi}(t)=\mathcal{U}_{\chi}(t,t_{0})=\mathcal{U}_{0}(t+\chi/\omega,t_{0}-\chi/\omega),$
(82)

from which we can easily show that

$\displaystyle
i\frac{d}{d\chi}\mathcal{U}_{\chi}(t)=\frac{1}{\omega}\left[\mathcal{H}(t)\mathcal{U}_{\chi}(t)-\mathcal{U}_{\chi}(t)\mathcal{H}(t_{0})\right].$
(83)

Consequently,

$\displaystyle\omega\kappa_{dy,1}(t)$ $\displaystyle=$
$\displaystyle\left<\mathcal{H}(t)\right>_{t}-\left<\mathcal{H}(0)\right>_{t_{0}}.$
(84)

Deriving with respect to time, we readily find the energy current from the
matter system to the photonic mode

$\displaystyle\omega\frac{d}{dt}\kappa_{dy,1}(t)$ $\displaystyle=$
$\displaystyle\left<\frac{d}{dt}\mathcal{H}(t)\right>_{t0}\equiv I(t),$ (85)

which is the expression of the semiclassical energy current commonly used in
the literature.

The generalization of Eq. (85) 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^{\prime}\neq k$ in the semiclassical Hamiltonian in Eq. (3), such that

$\mathcal{H}(t)=\mathcal{H}_{k}(t)+H_{0}+\sum_{k^{\prime}}H_{k^{\prime}}\left(\hat{a}_{k^{\prime}}+\hat{a}_{k^{\prime}}^{\dagger}\right),$
(86)

where $\mathcal{H}_{k}(t)=H_{k}\cos(\omega_{k}t+\varphi_{k})$. This is just
the Hamiltonian in Eq. (3) with a complicated $H_{0}$ and a single
semiclassical photonic mode. Using the same reasoning as before, we can
directly conclude that the photon flux into mode $k$ is given as

$\displaystyle
I_{k}(t)\equiv\frac{d}{dt}\omega_{k}\left<\Delta\hat{N}_{k}\right>$
$\displaystyle=$ $\displaystyle\left<\frac{d}{dt}\mathcal{H}_{k}(t)\right>.$
(87)

This implies that the photon-number change in mode $k$ is captured by the
operator

$\Delta\hat{N}_{k}=\int_{t_{0}}^{t}\frac{d}{dt}\mathcal{H}_{k}(t)dt,$ (88)

which is thus an operator in space and time. Clearly, an analysis of the
transport dynamics in terms of $\Delta\hat{N}_{k}$ does not offer the
convenience of the PRFT.

##  Appendix B Standard Full-Counting Statistics

Here we review the derivation of the standard FCS [68, 125]. We assume that
the system is described by Hamiltonian $H_{Q}$ in Eq. (1). We want to count
the change of photons in photonic mode $\hat{a}$ between times $t_{0}$ and
$t_{1}$. The system is assumed to be initially in state $\rho(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
$\hat{P}_{n_{0}}=\left|n_{0}\right>\left<n_{0}\right|$, where
$\left|n_{0}\right>$ denote the Fock states of $\hat{a}$, to determine the
initial photon number $n_{0}$. The resulting density operator of the light-
matter system is given as

$\displaystyle\rho_{pr}(t_{0})$ $\displaystyle=$
$\displaystyle\sum_{n_{0}}\hat{P}_{n_{0}}\rho(t_{0})\hat{P}_{n_{0}}$ (89)
$\displaystyle=$ $\displaystyle\sum_{n_{0}}p_{n_{0}}(t_{0})\rho_{pr,n_{0}},$

where $p_{n_{0}}(t_{0})$ denotes the probabilities to find $n_{0}$ photons.
The density matrix $\rho_{pr,n_{0}}$ is the system state conditioned on the
measurement outcome $n_{0}$.

We define $p_{n\mid n_{0}}(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

$\displaystyle p_{n\mid n_{0}}(t)$ $\displaystyle=$
$\displaystyle\text{Tr}\left[\rho_{pr,n_{0}}(t)\sum_{m}\delta_{n,m}\hat{P}_{m}\right]$
(90) $\displaystyle=$
$\displaystyle\text{Tr}\left[\rho_{pr,n_{0}}(t)\frac{1}{2\pi}\int_{-\pi}^{\pi}d\chi
e^{-i\chi n}\sum_{m}e^{i\chi m}\hat{P}_{m}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}d\chi e^{-i\chi
n}\text{Tr}\left[\rho_{pr,n_{0}}(t)e^{i\chi\hat{N}}\right],$

where $\hat{N}=\hat{a}^{\dagger}\hat{a}$. The time-evolved density matrix is
given as $\rho_{pr,n_{0}}(t)=\hat{U}(t)\rho_{pr,n_{0}}\hat{U}^{\dagger}(t)$.
For the density matrix conditioned on the first projective measurement, we can
replace

$\displaystyle\rho_{pr,n_{0}}=\rho_{pr,n_{0}}e^{i\chi\hat{n}_{0}}e^{-i\chi\hat{N}}.$
(91)

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

$\displaystyle p_{n\mid n_{0}}$ $\displaystyle=$
$\displaystyle\int_{-\pi}^{\pi}\frac{d\chi}{2\pi}e^{i\chi(n-n_{0})}\text{Tr}\left[\hat{U}\rho_{pr,n_{0}}e^{-i\chi\hat{N}}\hat{U}^{\dagger}e^{i\chi\hat{N}}\right].$

We are interested in the probability distribution of the photon number change
$\Delta n=n-n_{0}$, that we denote as $\overline{p}_{\Delta n}$. This
distribution is the average of the conditional probability distribution
$p_{n_{0}+\Delta n\mid n_{0}}$ weighted by the initial probability
distribution, i.e.,

$\displaystyle\overline{p}_{\Delta n}$ $\displaystyle\equiv$
$\displaystyle\sum_{n_{0}}p_{n_{0}+\Delta n\mid n_{0}}(t)p_{n_{0}}(t_{0})$
(92) $\displaystyle=$ $\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}d\chi
e^{i\chi\Delta
n}\text{Tr}\left[\hat{U}(t)\rho_{pr}(t_{0})e^{-i\chi\hat{N}}\hat{U}^{\dagger}(t)e^{i\chi\hat{N}}\right]$
$\displaystyle=$ $\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}d\chi
e^{i\chi\Delta n}M_{dy,\chi}^{(pr)}(t).$

In the last step, we have introduced the dynamical moment-generating function
of the two-point projective measurement, that can be alternatively expressed
as

$M_{dy,\chi}^{(pr)}(t)=\text{Tr}\left[\hat{U}_{\chi/2}(t)\rho_{pr}(t_{0})\hat{U}_{-\chi/2}^{\dagger}(t)\right]$
(93)

in terms of the generalized time-evolution operator

$\displaystyle\hat{U}_{\chi}(t)$ $\displaystyle=$ $\displaystyle
e^{i\chi\hat{N}}\hat{U}(t)e^{-i\chi\hat{N}}$ (94) $\displaystyle=$
$\displaystyle e^{-iH_{Q,\chi}t}.$

In the second line, we have expressed the time-evolution operator in terms of
the Hamiltonian

$H_{Q,\chi}=e^{i\chi\hat{N}}H_{Q}e^{-i\chi\hat{N}}.$ (95)

We emphasize that the projective measurement in Eq. (89) destroys coherences
in the photon basis and thus modifies the initial state. If the initial state
$\rho(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. (89). When one blindly replaces
$H_{Q,\chi}\rightarrow\mathcal{H}_{\chi}(t)$ in Eq. (94) with the
semiclassical Hamiltonian $\mathcal{H}_{\chi}(t)$ to mimic a coherent field,
the resulting expression is not equivalent to the dynamical cumulant
generating function of the PRFT in Eq. (13) and thus makes wrong predictions
about the photon statistics.

##  Appendix C Dynamics in Fock space

In this appendix, we provide details about the numerical implementation of the
Rabi and the two-mode Jaynes-Cummings models in Fock space. The numerical
results are compared to the predictions of the PRFT and depicted in Figs. 2,
3, and 7.

###  C.1 Fock space for large photon numbers

To reduce the numerical effort, we evaluate the matrix elements of the
creation and annihilation operators in the semiclassical limit, e.g., when
acting on photonic states with high occupation numbers as specified in Sec.
II.1. Thus, we use that

$\displaystyle\hat{a}^{\dagger}\left(\sum_{m}a_{m}\left|m\right>\right)$
$\displaystyle=$
$\displaystyle\sum_{n}\sqrt{n+1}\left|n+1\right>\left<n\right|\left(\sum_{m}a_{m}\left|m\right>\right)$
(96) $\displaystyle\rightarrow$
$\displaystyle\sum_{n}\alpha\left|n+1\right>\left<n\right|\left(\sum_{m}a_{m}\left|m\right>\right),$

which is justified, if the coefficient $a_{n}$ determining the dynamics are
only significant for $n$ within the range
$\left[\alpha^{2}-x\cdot\alpha,\alpha^{2}+x\cdot\alpha\right]$, where
$\alpha\gg 1$ and $x>0$. To see this, let us consider the ratio of the matrix
elements of the creation operator $C_{n_{i}}=\sqrt{n_{i}+1}$ for
$\left|n_{i}-\alpha^{2}\right|<x\sqrt{N}$. When representing
$n_{i}=\alpha^{2}+m_{i}-1$, we find that

$\displaystyle\lim_{\alpha\rightarrow\infty}\frac{C_{n_{1}}}{C_{n_{2}}}$
$\displaystyle=$ $\displaystyle\frac{\sqrt{N+m_{1}}}{\sqrt{N+m_{2}}}=1,$ (97)

for all allowed $m_{1},m_{2}$, clearly showing that the replacement in Eq.
(96) becomes exact. We note that this replacement is equivalent to solving the
time-dependent Schrödinger equation in Sambe space.

###  C.2 Numerical simulation of the two-mode Jaynes-Cummings model

It is easy to see that the excitation number operator
$\hat{N}_{ex}=\hat{\sigma}_{z}+\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2}$
is an integral of motion of the two-mode Jaynes-Cummings model in Eq. (34). To
take advantage of this symmetry, we introduce the basis states

$\displaystyle\left|s,n\right>_{N_{ex}}\equiv\left|s\right>_{z}\otimes\left|n\right>_{1}\otimes\left|N_{ex}-n-s\right>_{2},$
(98)

where $\left|s\right>_{z}$ with $s=\pm 1$ are the spin states of the atomic
system and $\left|n\right>_{k}$ denote the Fock states of the photonic modes
$k=1,2$. The matrix elements of the two-mode James-Cummings model in this
basis read as

${}_{N_{ex}}\left<s,n\right|\hat{H}\left|s^{\prime},n^{\prime}\right>_{N_{ex}}$
$\displaystyle=$
$\displaystyle\begin{cases}\omega_{0}\sigma+\omega(N_{ex}-s)&s=s^{\prime},\;n=n^{\prime}\\\
\tilde{g}_{1}\sqrt{n^{\prime}+1}&n=n^{\prime}+1,\;s=s^{\prime}-1\\\
\tilde{g}_{1}\sqrt{n^{\prime}}&n=n^{\prime}-1,\;s=s^{\prime}+1\\\
\tilde{g}_{2}\sqrt{N_{ex}-n^{\prime}-s+1}&n=n^{\prime},\;s=s^{\prime}-1\\\
\tilde{g}_{2}\sqrt{N_{ex}-n^{\prime}-s}&n=n^{\prime},\;s=s^{\prime}+1\\\
\end{cases},$ (99)

where the $n$ dependence of the creation and annihilation operator is still
maintained. Evaluating the matrix elements in the semiclassical limit by
replacing all $\sqrt{n}\rightarrow\alpha_{k}$ and introducing the rescaled
light-matter couplings $\tilde{g}_{k}\alpha_{k}\rightarrow g_{k}$, the
Hamiltonian becomes translationally invariant with respect to $n$.

For coherent states, the initial state can be expanded in the basis Eq. (98)
as

$\displaystyle\left|\psi(0)\right>$ $\displaystyle=$
$\displaystyle\left|\phi\right>\otimes\left|\alpha_{1}e^{i\varphi_{1}}\right>\otimes\left|\alpha_{2}e^{i\varphi_{2}}\right>$
(100) $\displaystyle=$ $\displaystyle\left|\phi\right>\otimes
e^{-\alpha_{1}}\sum_{n_{1}}\frac{\alpha_{1}^{n_{1}}}{n_{1}!}e^{i\varphi_{1}n_{1}}\left|n_{1}\right>$
$\displaystyle\otimes
e^{-\alpha_{2}}\sum_{n_{2}}\frac{\alpha_{2}^{n_{2}}}{n_{2}!}e^{i\varphi_{2}n_{2}}\left|n_{2}\right>$
$\displaystyle=$
$\displaystyle\sum_{s,n_{1},n_{2}}b_{s}e^{-\alpha_{1}-\alpha_{2}}\frac{\alpha_{1}^{n_{1}}}{n_{1}!}\frac{\alpha_{2}^{n_{2}}}{n_{2}!}e^{i\varphi_{1}n_{1}+i\varphi_{2}n_{2}}\left|s,n_{1},n_{2}\right>$
$\displaystyle=$
$\displaystyle\sum_{s,N,n}b_{s}e^{-\alpha_{1}-\alpha_{2}}\frac{\alpha_{1}^{n}}{n!}\frac{\alpha_{2}^{N_{ex}-n}}{(N_{ex}-n)!}$
$\displaystyle\times
e^{i\varphi_{1}n_{1}+i\varphi_{2}n_{2}}\left|s,n\right>_{N}$
$\displaystyle\approx$
$\displaystyle\frac{1}{\mathcal{N}}\sum_{N_{ex}}\sum_{s,n}b_{s}e^{\frac{-(n-\alpha_{1}^{2})^{2}}{4\sigma_{1}^{2}}-\frac{(N_{ex}-n-\alpha_{2})^{2}}{4\sigma_{2}^{2}}}$
$\displaystyle\times$ $\displaystyle
e^{i\varphi_{1}n+i\varphi_{2}(N_{ex}-n)}\left|s,n\right>_{N_{ex}}$
$\displaystyle=$ $\displaystyle\sum_{N_{ex}}\left|N_{ex}\right>,$

where the $\left|N_{ex}\right>$ are unnormalized, but mutually orthogonal
states. In the third equality, we have expanded the spin system in the
$\hat{\sigma}_{z}$ basis as
$\left|\phi\right>=b_{+1}\left|+1\right>_{z}+b_{-1}\left|-1\right>_{z}$. In
the last step, we have used the approximation that the Poission function can
be replaced by a Gaussian function with mean $\alpha_{k}^{2}$ and the standard
deviation $\sigma_{k}=\alpha_{k}$ for large $\alpha_{k}\gg 1$. Thereby,
$\mathcal{N}$ guarantees a proper normalization. The dynamics can be now
solved independently for each $N_{ex}$ using the Hamiltonian Eq. (99), which
is independent of $N_{ex}$ in the semiclassical limit.

###  C.3 Benchmark calculations in the Jaynes-Cummings model

In the benchmark simulations, we assume number-squeezed states
$\sigma_{k}\ll\alpha_{k}$ in Eq. (100) with $\sigma_{k}=\sqrt{1000}\approx
30$. The phases are $\varphi_{1}=0.25\pi$ and $\varphi_{2}=0$. As for large
photon numbers the Hamiltonian in the semiclassical limit is independent of
$n$ and depends only on the products $g_{k}=\tilde{g}_{k}\alpha_{k}$, we do
not specify the $\alpha_{k}$ but only the effective light-matter interaction
parameters $g_{j}$. This is a demanding initial condition for verifying the
semiclassical limit approximation in Eq. (76), as for such number-squeezed
states the modulus of the expansion coefficients $a_{n}$ changes faster than
for coherent states. In Figs. 3 and 7, we compare the time evolution predicted
by the PRFT with the time evolution in Fock space for different light-matter
couplings $g_{2}$ and photon energies $\omega$. We find an excellent agreement
in all calculations despite the strong number squeezing.

Figure 7: Benchmark calculations in the Jaynes-Cummings model for various
parameters. The black colored curves depict the prediction of the PRFT and the
solid lines depict the simulation in Fock space. The change of the
probabilities $\Delta p_{n_{1}}=p_{n_{1}}(t)-p_{n_{1}}(0)$ in photon mode
$j=1$ are depicted for three different time steps $h_{z}\cdot t=0,300,600$.
The light-matter interaction is $g_{1}=0.05h_{z}$, and the initial state is
characterized by $\alpha_{k}^{2}=1000$, $\varphi_{1}=0.25\pi$ and
$\varphi_{2}=0$. The parameters for $g_{2}$ and $\omega$ are indicated in the
figure.

##  Appendix D Calculation of the purity in the two-mode Jaynes-Cummings model

In this Appendix, we apply the PRFT to calculate the purity in the Jaynes-
Cummings model as a function of time for the light-matter state

$\left|\Psi\right>=\alpha\left|u_{1}\right>\left|A_{1}\right>+\beta\left|u_{2}\right>\left|A_{2}\right>,$
(101)

where $\left|\alpha\right|^{2}+\left|\beta\right|^{2}$. The matter states
$\left|u_{\mu}\right>$ represent Floquet states. The photonic states
$\left|A_{\mu}\right>$ can have an arbitrary overlap. Because of the
excitation-number conservation in the Jaynes-Cummings model explained in
Appendix C.2, the photonic state can be expanded in orthogonal subspaces with
fixed $N_{ex}$. The Floquet-state conditioned probability to be in a
particular state Eq. (98) specified by $N_{ex}$ and $n$ is
$p_{N_{ex},n\mid\mu}=\left|\left<N_{ex},n\mid A_{\mu}\right>\right|^{2}$. As
$n=n_{1}$ by definition of the states in Eq. (98), the mean of the
distribution $p_{n}=\sum_{N_{ex}}p_{N_{ex},n\mid\mu}$ increases linearly with
time $\kappa_{dy,1}^{(\mu)}(t)=-E_{\mu}^{\prime}(\varphi)t$, while the
variance stays constant according to the explanations in Sec. III.1. From the
investigation of the two-mode Jaynes-Cummings model in Sec. III.2, we know
that higher-order cumulants grow only in sub-leading order and have thus an
minor impact on the probability distribution. Thus, according to PRFT, the
probability distribution approximately reads as
$p_{N_{ex},n\mid\mu}(t)=p_{N_{ex}}\cdot
e^{(n-\overline{n}_{\mu}(t))/2\sigma^{2}}\sqrt{\pi\sigma}$, where $\sigma^{2}$
is the initial variance,
$\overline{n}_{\mu}(t)=\overline{n}_{\mu}(0)+\kappa_{dy,1}^{(\mu)}(t)$ is the
time-dependent mean, and $p_{N_{ex}}$ is determined by the projection of the
initial state onto the subspace with excitation number $N_{ex}$. The
multiplicative form of $p_{N_{ex},n\mid\mu}(t)$ is a consequence of the
$N_{ex}$ independence of the Hamiltonian [see explanations in Appendix C.2]
and the initial state in Eq. (100).

To calculate the purity, we evaluate the reduced density matrix of
$\rho=\left|\Psi\right>\left<\Psi\right|$, i.e.,

$\displaystyle\rho_{M}$ $\displaystyle=$
$\displaystyle\text{Tr}_{L}\left(\rho\right)$ (104) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\left|\alpha\right|^{2}&\alpha^{*}\beta
c\\\ \beta^{*}\alpha c^{*}&\left|\beta\right|^{2}\end{array}\right),$

where $c=\left<A_{1}\mid A_{2}\right>$ is the overlap of the two photonic
states. For the following calculation, we use that the phases in the Fock
state expansion of both states are equal, meaning that $\arg\left<N_{ex},n\mid
A_{1}\right>=\arg\left<N_{ex},n\mid A_{2}\right>$, which is consistent with
the PRFT. In this case, the overlap $c$ is completely determined by the
conditional probabilities $p_{N_{ex},n\mid\mu}(t)$. Considering the photon
number as a continuous variable, we can then evaluate

$\displaystyle c$ $\displaystyle=$
$\displaystyle\sum_{N_{ex}}\frac{1}{\sigma\sqrt{\pi}}\int
dn\;p_{N_{ex}}e^{-\frac{(n-\overline{n}_{1})^{2}}{4\sigma^{2}}}e^{-\frac{(n-\overline{n}_{2})^{2}}{4\sigma^{2}}}$
(105) $\displaystyle=$ $\displaystyle\frac{1}{\sigma\sqrt{\pi}}\int
dne^{-\frac{2n^{2}-2n\left(n_{1}+n_{2}\right)+\overline{n}_{1}^{2}+\overline{n}_{2}^{2}}{4\sigma^{2}}}$
$\displaystyle=$ $\displaystyle\frac{1}{\sigma\sqrt{\pi}}\int
dne^{-\frac{n^{2}-n\left(n_{1}+n_{2}\right)+\frac{1}{2}\left(\overline{n}_{1}^{2}+\overline{n}_{2}^{2}\right)}{2\sigma^{2}}}$
$\displaystyle=$ $\displaystyle\frac{1}{\sigma\sqrt{\pi}}\int
dne^{-\frac{\left[n-\frac{1}{2}\left(\overline{n}_{1}+\overline{n}_{2}\right)\right]^{2}}{4\sigma^{2}}}e^{-\frac{\left(\overline{n}_{1}-\overline{n}_{2}\right)^{2}}{8\sigma^{2}}}$
$\displaystyle=$ $\displaystyle
e^{-\frac{\left(\overline{n}_{1}-\overline{n}_{2}\right)^{2}}{8\sigma^{2}}},$

where the time-dependent mean difference is given as
$\left[\overline{n}_{1}(t)-\overline{n}_{2}(t)\right]=\left[E_{2}^{\prime}(\varphi)-E_{1}^{\prime}(\varphi)\right]t$.

The purity is related to the eigenvalues $p_{r}$ of the reduced density matrix
via $\mathcal{P}=\text{Tr}\left[\rho_{M}^{2}\right]=p_{1}^{2}+p_{2}^{2}$,
where the eigenvalues are explicitly given as

$p_{1/2}=\frac{1}{2}\pm\frac{1}{2}\sqrt{\left(\left|\alpha\right|^{2}-\left|\beta\right|^{2}\right)^{2}+4\left|\alpha\right|^{2}\left|\beta\right|^{2}\left|c\right|^{2}}.$
(106)

In doing so, the purity finally becomes

$\displaystyle\mathcal{P}$ $\displaystyle=$
$\displaystyle\frac{1+\left(\left|\alpha\right|^{2}-\left|\beta\right|^{2}\right)^{2}}{2}+2\left|\alpha\right|^{2}\left|\beta\right|^{2}\left|c\right|^{2},$
(107)

which is the expression in Eq. (LABEL:eq:timeDependentPurity) after an
appropriate parameter replacement.

## References

  * Rudner  _et al._ [2013] M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin, Anomalous edge states and the bulk-edge correspondence for periodically driven two-dimensional systems, Phys. Rev. X 3, 031005 (2013).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Photonic Floquet topological insulators, Nature 496, 196 (2013).
  * Lindner  _et al._ [2011] N. H. Lindner, G. Refael, and V. Galitski, Floquet topological insulator in semiconductor quantum wells, Nature Physics 7, 490 (2011).
  * Roy and Harper [2017] R. Roy and F. Harper, Floquet topological phases with symmetry in all dimensions, Phys. Rev. B 95, 195128 (2017).
  * Lababidi  _et al._ [2014] M. Lababidi, I. I. Satija, and E. Zhao, Counter-propagating edge modes and topological phases of a kicked quantum hall system, Phys. Rev. Lett. 112, 026805 (2014).
  * Rudner and Lindner [2020] M. S. Rudner and N. H. Lindner, Band structure engineering and non-equilibrium dynamics in Floquet topological insulators, Nature Reviews Physics 2, 229 (2020).
  * Wintersperger  _et al._ [2020] K. Wintersperger, C. Braun, F. N. Ünal, A. Eckardt, M. D. Liberto, N. Goldman, I. Bloch, and M. Aidelsburger, Realization of an anomalous Floquet topological system with ultracold atoms, Nature Physics 16, 1058 (2020).
  * Maczewsky  _et al._ [2017] L. J. Maczewsky, J. M. Zeuner, S. Nolte, and A. Szameit, Observation of photonic anomalous Floquet topological insulators, Nature communications 8, 1 (2017).
  * Engelhardt  _et al._ [2016] G. Engelhardt, M. Benito, G. Platero, and T. Brandes, Topological instabilities in ac-driven bosonic systems, Phys. Rev. Lett. 117, 045302 (2016).
  * Engelhardt  _et al._ [2017] G. Engelhardt, M. Benito, G. Platero, G. Schaller, and T. Brandes, Random-walk topological transition revealed via electron counting, Phys. Rev. B 96, 241404 (2017).
  * Hu  _et al._ [2020] H. Hu, B. Huang, E. Zhao, and W. V. Liu, Dynamical singularities of floquet higher-order topological insulators, Phys. Rev. Lett. 124, 057001 (2020).
  * Huang and Liu [2020] B. Huang and W. V. Liu, Floquet higher-order topological insulators with anomalous dynamical polarization, Phys. Rev. Lett. 124, 216601 (2020).
  * Pyrialakos  _et al._ [2022] G. G. Pyrialakos, J. Beck, M. Heinrich, L. J. Maczewsky, N. V. Kantartzis, M. Khajavikhan, A. Szameit, and D. N. Christodoulides, Bimorphic Floquet topological insulators, Nature Materials 21, 634 (2022).
  * Potter  _et al._ [2016] A. C. Potter, T. Morimoto, and A. Vishwanath, Classification of interacting topological Floquet phases in one dimension, Phys. Rev. X 6, 041001 (2016).
  * Engelhardt  _et al._ [2013] G. Engelhardt, V. M. Bastidas, C. Emary, and T. Brandes, ac-driven quantum phase transition in the Lipkin-Meshkov-Glick model, Phys. Rev. E 87, 052110 (2013).
  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: A review, Rep. Prog. Phys. 81, 016401 (2017).
  * Sacha [2020] K. Sacha,  _Time Crystals_ , Springer Series on Atomic, Optical, and Plasma Physics, Vol. 114 (Springer, Cham, Switzerland, 2020).
  * Else  _et al._ [2020] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Discrete time crystals, Annu. Rev. Condens. Matter Phys. 11, 467 (2020).
  * Khemani  _et al._ [2019] V. Khemani, R. Moessner, and S. L. Sondhi, A brief history of time crystals, arXiv:1910.10745 (2019).
  * Sacha [2015] K. Sacha, Modeling spontaneous breaking of time-translation symmetry, Phys. Rev. A 91, 033617 (2015).
  * von Keyserlingk  _et al._ [2016] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in Floquet systems, Phys. Rev. B 94, 085112 (2016).
  * Khemani  _et al._ [2016a] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Phys. Rev. Lett. 116, 250401 (2016a).
  * Yao  _et al._ [2017a] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Discrete time crystals: Rigidity, criticality, and realizations, Phys. Rev. Lett. 118, 030401 (2017a).
  * Else  _et al._ [2016a] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117, 090402 (2016a).
  * Zhang  _et al._ [2017] J. Zhang, P. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Observation of a discrete time crystal, Nature 543, 217 (2017).
  * Huang  _et al._ [2018] B. Huang, Y.-H. Wu, and W. V. Liu, Clean Floquet time crystals: models and realizations in cold atoms, Phys. Rev. Lett. 120, 110603 (2018).
  * Choudhury [2021] S. Choudhury, Route to extend the lifetime of a discrete time crystal in a finite spin chain without disorder, Atoms 9, 25 (2021).
  * Choi  _et al._ [2017a] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, C. von Keyserlingk, N. Y. Yao, E. Demler, and M. D. Lukin, Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017a).
  * Russomanno  _et al._ [2017] A. Russomanno, F. Iemini, M. Dalmonte, and R. Fazio, Floquet time crystal in the Lipkin-Meshkov-Glick model, Phys. Rev. B 95, 214307 (2017).
  * Yang and Cai [2021] X. Yang and Z. Cai, Dynamical transitions and critical behavior between discrete time crystal phases, Phys. Rev. Lett. 126, 020602 (2021).
  * Peña  _et al._ [2021] R. Peña, V. M. Bastidas, F. Torres, W. J. Munro, and G. Romero, Fractional resonances and prethermal states in Floquet systems, arXiv e-prints , arXiv:2111.06949 (2021), arXiv:2111.06949 [quant-ph] .
  * Engelhardt  _et al._ [2019] G. Engelhardt, G. Platero, and J. Cao, Discontinuities in Driven Spin-Boson Systems due to Coherent Destruction of Tunneling: Breakdown of the Floquet-Gibbs Distribution, Phys. Rev. Lett. 123, 120602 (2019).
  * Gu and Franco [2018] B. Gu and I. Franco, Optical absorption properties of laser-driven matter, Phys. Rev. A 98, 063412 (2018).
  * Chen  _et al._ [2021] M.-B. Chen, B.-C. Wang, S. Kohler, Y. Kang, T. Lin, S.-S. Gu, H.-O. Li, G.-C. Guo, X. Hu, H.-W. Jiang, G. Cao, and G.-P. Guo, Floquet state depletion in ac-driven circuit qed, Phys. Rev. B 103, 205428 (2021).
  * Yan  _et al._ [2019] Y. Yan, Z. Lü, J. Y. Luo, and H. Zheng, Role of generalized parity in the symmetry of the fluorescence spectrum from two-level systems under periodic frequency modulation, Phys. Rev. A 100, 013823 (2019).
  * Lyu  _et al._ [2020] C. Lyu, S. Choudhury, C. Lv, Y. Yan, and Q. Zhou, Eternal discrete time crystal beating the Heisenberg limit, Phys. Rev. Research 2, 033070 (2020).
  * Zhang  _et al._ [2021] K. Zhang, W. Liang, P. Meystre, and W. Zhang, Quantum metrology with Bloch oscillations in Floquet phase space, arXiv preprint arXiv:2111.10506 (2021).
  * Choi  _et al._ [2017b] S. Choi, N. Y. Yao, and M. D. Lukin, Quantum metrology based on strongly correlated matter, arXiv preprint arXiv:1801.00042 (2017b).
  * Schweizer  _et al._ [2019] C. Schweizer, F. Grusdt, M. Berngruber, L. Barbiero, E. Demler, N. Goldman, I. Bloch, and M. Aidelsburger, Floquet approach to $\mathbb{Z}_{2}$ lattice gauge theories with ultracold atoms in optical lattices, Nature Physics 15, 1168 (2019).
  * Görg  _et al._ [2019] F. Görg, K. Sandholzer, J. Minguzzi, R. Desbuquois, M. Messer, and T. Esslinger, Realization of density-dependent Peierls phases to engineer quantized gauge fields coupled to ultracold matter, Nature Physics 15, 1161 (2019).
  * Arnal  _et al._ [2020] M. Arnal, G. Chatelain, M. Martinez, N. Dupont, O. Giraud, D. Ullmo, B. Georgeot, G. Lemarié, J. Billy, and D. Guéry-Odelin, Chaos-assisted tunneling resonances in a synthetic Floquet superlattice, Science advances 6, eabc4886 (2020).
  * Weitenberg and Simonet [2021] C. Weitenberg and J. Simonet, Tailoring quantum gases by Floquet engineering, Nature Physics 17, 1342 (2021).
  * Choudhury and Mueller [2014] S. Choudhury and E. J. Mueller, Stability of a Floquet Bose-Einstein condensate in a one-dimensional optical lattice, Phys. Rev. A 90, 013621 (2014).
  * Choudhury and Mueller [2015a] S. Choudhury and E. J. Mueller, Stability of a bose-einstein condensate in a driven optical lattice: Crossover between weak and tight transverse confinement, Phys. Rev. A 92, 063639 (2015a).
  * Choudhury and Mueller [2015b] S. Choudhury and E. J. Mueller, Transverse collisional instabilities of a Bose-Einstein condensate in a driven one-dimensional lattice, Phys. Rev. A 91, 023624 (2015b).
  * Geier  _et al._ [2021] S. Geier, N. Thaicharoen, C. Hainaut, T. Franz, A. Salzinger, A. Tebben, D. Grimshandl, G. Zürn, and M. Weidemüller, Floquet Hamiltonian engineering of an isolated many-body spin system, Science 374, 1149 (2021).
  * Sträter  _et al._ [2016] C. Sträter, S. C. Srivastava, and A. Eckardt, Floquet realization and signatures of one-dimensional anyons in an optical lattice, Phys. Rev. Lett. 117, 205303 (2016).
  * Meinert  _et al._ [2016] F. Meinert, M. J. Mark, K. Lauber, A. J. Daley, and H.-C. Nägerl, Floquet engineering of correlated tunneling in the Bose-Hubbard model with ultracold atoms, Phys. Rev. Lett. 116, 205301 (2016).
  * Clark  _et al._ [2018] L. W. Clark, B. M. Anderson, L. Feng, A. Gaj, K. Levin, and C. Chin, Observation of density-dependent gauge fields in a Bose-Einstein condensate based on micromotion control in a shaken two-dimensional lattice, Phys. Rev. Lett. 121, 030402 (2018).
  * Görg  _et al._ [2018] F. Görg, M. Messer, K. Sandholzer, G. Jotzu, R. Desbuquois, and T. Esslinger, Enhancement and sign change of magnetic correlations in a driven quantum many-body system, Nature 553, 481 (2018).
  * Shirley [1965] J. H. Shirley, Solution of the Schrödinger equation with a Hamiltonian periodic in time, Phys. Rev. 138, B979 (1965).
  * Guérin  _et al._ [1997] S. Guérin, F. Monti, J.-M. Dupont, and H. R. Jauslin, On the relation between cavity-dressed states, Floquet states, RWA and semiclassical models, J. Phys. A: Math. Gen. 30, 7193 (1997).
  * Sentef  _et al._ [2020] M. A. Sentef, J. Li, F. Künzel, and M. Eckstein, Quantum to classical crossover of Floquet engineering in correlated quantum systems, Phys. Rev. Research 2, 033033 (2020).
  * L. and M. [2008] M. L. and W. M.,  _Optical coherence and quantum optics_ (Cambridge University Press, New York, 2008).
  * Hübener  _et al._ [2021] H. Hübener, U. De Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, Engineering quantum materials with chiral optical cavities, Nature materials 20, 438 (2021).
  * Wang  _et al._ [2019] X. Wang, E. Ronca, and M. A. Sentef, Cavity quantum electrodynamical Chern insulator: Towards light-induced quantized anomalous Hall effect in graphene, Phys. Rev. B 99, 235156 (2019).
  * Schlawin  _et al._ [2022] F. Schlawin, D. M. Kennes, and M. A. Sentef, Cavity quantum materials, Applied Physics Reviews 9, 011312 (2022).
  * Li and Eckstein [2020] J. Li and M. Eckstein, Manipulating intertwined orders in solids with quantum light, Phys. Rev. Lett. 125, 217402 (2020).
  * Kiffner  _et al._ [2019] M. Kiffner, J. R. Coulthard, F. Schlawin, A. Ardavan, and D. Jaksch, Manipulating quantum materials with quantum light, Phys. Rev. B 99, 085116 (2019).
  * Gulácsi and Dóra [2015] B. Gulácsi and B. Dóra, From Floquet to Dicke: quantum spin Hall insulator interacting with quantum light, Phys. Rev. Lett. 115, 160402 (2015).
  * Schäfer  _et al._ [2018] C. Schäfer, M. Ruggenthaler, and A. Rubio, Ab initio nonrelativistic quantum electrodynamics: Bridging quantum chemistry and quantum optics from weak to strong coupling, Phys. Rev. A 98, 043801 (2018).
  * Eckhardt  _et al._ [2022] C. J. Eckhardt, G. Passetti, M. Othman, C. Karrasch, F. Cavaliere, M. A. Sentef, and D. M. Kennes, Quantum Floquet engineering with an exactly solvable tight-binding chain in a cavity, Communications Physics 5, 1 (2022).
  * Rokaj  _et al._ [2021] V. Rokaj, M. Penz, M. A. Sentef, M. Ruggenthaler, and A. Rubio, Polaritonic Hofstadter butterfly and cavity-control of the quantized Hall conductance, arXiv preprint arXiv:2109.15075 (2021).
  * Lloyd-Hughes  _et al._ [2021] J. Lloyd-Hughes, P. M. Oppeneer, T. P. Dos Santos, A. Schleife, S. Meng, M. A. Sentef, M. Ruggenthaler, A. Rubio, I. Radu, M. Murnane, _et al._ , The 2021 ultrafast spectroscopic probes of condensed matter roadmap, Journal of Physics: Condensed Matter 33, 353001 (2021).
  * Mandel and Wolf [1995] L. Mandel and E. Wolf,  _Optical Coherence and Quantum Optics_ (Cambridge University Press, 1995).
  * Mandel [1979] L. Mandel, Sub-Poissonian photon statistics in resonance fluorescence, Optics letters 4, 205 (1979).
  * Cook [1981] R. J. Cook, Photon number statistics in resonance fluorescence, Phys. Rev. A 23, 1243 (1981).
  * Levitov  _et al._ [1996] L. S. Levitov, H. Lee, and G. B. Lesovik, Electron counting statistics and coherent states of electric current, Journal of Mathematical Physics 37, 4845 (1996).
  * Brandes [2008] T. Brandes, Waiting times and noise in single particle transport, Ann. Phys. (Berlin) 17, 477 (2008).
  * Böhling  _et al._ [2018] S. Böhling, G. Engelhardt, G. Platero, and G. Schaller, Thermoelectric performance of topological boundary modes, Phys. Rev. B 98, 035132 (2018).
  * Xue  _et al._ [2015] H.-B. Xue, H.-J. Jiao, J.-Q. Liang, and W.-M. Liu, Non-Markovian full counting statistics in quantum dot molecules, Scientific Reports 5, 8978 (2015).
  * Engelhardt and Cao [2019] G. Engelhardt and J. Cao, Tuning the Aharonov-Bohm effect with dephasing in nonequilibrium transport, Phys. Rev. B 99, 075436 (2019).
  * You  _et al._ [2019] J.-S. You, R. Schmidt, D. A. Ivanov, M. Knap, and E. Demler, Atomtronics with a spin: Statistics of spin transport and nonequilibrium orthogonality catastrophe in cold quantum gases, Phys. Rev. B 99, 214505 (2019).
  * Dasenbrook  _et al._ [2014] D. Dasenbrook, C. Flindt, and M. Büttiker, Floquet theory of electron waiting times in quantum-coherent conductors, Phys. Rev. Lett. 112, 146801 (2014).
  * Yadalam and Harbola [2016] H. K. Yadalam and U. Harbola, Statistics of an adiabatic charge pump, Phys. Rev. B 93, 035312 (2016).
  * Groth  _et al._ [2006] C. W. Groth, B. Michaelis, and C. Beenakker, Counting statistics of coherent population trapping in quantum dots, Phys. Rev. B 74, 125315 (2006).
  * Schaller  _et al._ [2018] G. Schaller, J. Cerrillo, G. Engelhardt, and P. Strasberg, Electronic Maxwell demon in the coherent strong-coupling regime, Phys. Rev. B 97, 195104 (2018).
  * Kleinherbers  _et al._ [2018] E. Kleinherbers, P. Stegmann, and J. König, Revealing attractive electron–electron interaction in a quantum dot by full counting statistics, New Journal of Physics 20, 073023 (2018).
  * Benito  _et al._ [2016] M. Benito, M. Niklas, and S. Kohler, Full-counting statistics of time-dependent conductors, Phys. Rev. B 94, 195433 (2016).