# Few-photon Fock-state wave packet interacting with a cavity-atom system in a waveguide: Exact quantum state dynamics

###### Abstract

In the paper, we employ a wavefunction approach to investigate the evolution of a two-photon wave packet propagating in a one-dimensional waveguide coupled to the Jaynes-Cummings (JC) system. We derive and solve, both analytically and numerically, a set of equations of motion governing the quantum state of the system. That allows us to provide real-time analysis of the evolution of the wave packet two-photon joint spectrum (2PJS) and the excitation dynamics of the JC system in the course of its interaction with the two-photon pulse. We demonstrate that the 2PJS and the spectrum of the wave packet scattered from the JC system experience transformation for nonzero atom-cavity couplings. Moreover, using Schmidt decomposition, we show that the scattered photons feature frequency entanglement contrary to the incident ones which are not entangled.

## I Introduction

Propagating photons can act as mobile qubits carrying quantum information between stationary qubits, where it can be stored or subjected to quantum logical operations north (). Remote stationary qubits are interconnected via waveguides (either microwave or optical ones) ensuring almost lossless transport of photons within quantum information processing (QIP) devices. Waveguides confine light in a transverse direction forming an effective one-dimensional (1D) continuum for propagating photons. Due to the Purcell effect, this modified electromagnetic environment gives rise to enhanced light-matter interaction. Strong light-matter coupling enables fast and accurate quantum-state transfer between photons and stationary qubits, which is essential for the realization of high-fidelity single- and multi-qubit gates.

Remarkable progress in creation and manipulation of waveguide quantum electrodynamics (wQED) systems was achieved in the last couple of decades. In addition, several designs of single- and few-photon sources operating in the microwave hou (); hof (); pfdiaz (); pfaff () and optical claud (); snij (); hansch () domains were demonstrated experimentally. That makes wQED a highly promising and versatile hardware platform for implementation of QIP systems. Ongoing technological advances in this field galvanize extensive theoretical studies of wQED systems. Several wQED-based architectures for quantum computation were proposed recently zheng2013 (); paulisch (); pichler (). Diverse theoretical techniques were engaged for studying a photon propagation in various wQED setups. Below we present a brief overview of some of these approaches.

Interaction of Fock- and coherent-state wave packets with a two-level atom (2LA) in a 1D waveguide was studied within the Heisenberg picture in Refs. dom (); sto2013 (); sto2014 (); rou2016 ().

A theoretical framework based on master equations was employed in Ref. bar () to describe the evolution of arbitrary quantum systems driven by -photon wave packets. In Ref. tshi () a generalized master equation approach was developed for studying multiphoton transport in wQED systems with non-Markovian couplings.

Different quantum-field-theoretic approaches such as real-space Bethe ansatz shen2007 (); shen2007a (); shen2015 (), Lippmann-Schwinger formalism zheng (), Lehmann-Symanzik-Zimmermann reduction shi2009 (); shi2011 (); shi2013 (), various diagrammatic techniques plet (); schn (); koc (); see (), and Dyson series summation hurst () were adopted for the description of wQED systems as well.

An input-output formalism gard () was harnessed in Ref. sfan2010 () to establish a relationship between the exact few-photon S matrix of a 2LA and photonic input-output operators. This successful approach was employed in further studies koc2012 (); reph (); xu2015 (); can (); triv ().

It should be noted that the S matrix provides information only about asymptotic (long-time) behavior of the system, while no insight into its transient dynamics is given. To describe the evolution of the system state one should rely on alternative techniques, e.g., wavefunction approach which was employed for studying the propagation of few-photon wave packets in wQED systems containing one or few 2LAs hof2003 (); ychen (); nyst (); konyk (), a single-photon pulse propagation in a chain of non-identical 2LAs liao2016 (), interaction of one- wang () and two-photon pulses qhu () with a cavity-atom system, and few-photon scattering in systems with quantum delayed feedback fang2018 (); calajo (). Detailed information about the system can be obtained if one knows its exact wavefunction at the arbitrary moment of time. This can be achieved by solving a set of coupled equations of motion for probability amplitudes governing the quantum-state evolution of the system. These first-order ordinary differential equations (ODEs) can be solved either numerically using one of the elaborate computational routines, or analytically, for some cases.

Many of the previous studies model a local emitter by a 2LA. It is the simplest saturable (i.e., nonlinear) quantum system since it can contain only one excitation per moment. This property was employed in Ref. rou () to derive an exact analytical solution for the problem of -photon scattering on a 2LA. Investigation of a multi-photon transport in wQED systems with more complex emitters constitutes a more intricate problem for analysis. A system composed of a cavity coupled to an atom, which is referred to as the Jaynes-Cummings (JC) system, is one of the most important model systems in quantum optics. It demonstrates a number of interesting phenomena such as photon blockade birn (); henn (); hamsen (), photon phase switching tiecke (), generation of photon bundles sanch (), etc.

The paper concerns a few-photon transport in the wQED setup with the emitter represented by the JC system. We develop a wavefunction approach to provide a fully quantum-mechanical real-time analysis of scattering of two-photon Fock-state wave packets on the JC system coupled to a 1D unidirectional (chiral) waveguide. The exact solutions of the equations of motion for the probability amplitudes determining a state of the system at the arbitrary moment of time are obtained both analytically and numerically. In the long-time limit of the system evolution when all interaction processes are over, and the scattered photons propagate in the waveguide as free excitations, we ultimately arrive at the results equivalent to those derived earlier using other approaches. We use Schmidt decomposition of the two-photon spectral distribution function (SDF) of the outgoing wave packet to demonstrate that the scattered photons are entangled.

The outline of the paper is as follows. In Sec. II a model system is described and its possible experimental incarnations are briefly discussed. In Sec. III we study the case of a single-photon ingoing pulse. Sections IV-VI contain the results regarding the case of a two-photon Fock-state ingoing pulse. In Sec. IV the equations of motion governing the quantum-state evolution of the system are derived. The dynamics of the JC system is studied in Section V. In Sec. VI we present the exact expression for the SDF of the wave packet scattered from the JC system. The properties of the spectrum of the outgoing photons are analyzed in Sec. VII. In Sec. VIII Schmidt decomposition is employed to investigate the frequency entanglement of the scattered photons. We summarize in Sec. IX. Details of derivation and solution of evolution equations are presented in Appendices.

## Ii The Model

We consider a model system consisting of a single-mode cavity coupled to an individual 2LA. A cavity mode overlaps with modes of a unidirectional 1D waveguide leading to that photons are able to leak in and out of the cavity. In order to populate the cavity and to excite the atom coupled to it, the former is driven by a few-photon wave packet propagating in a waveguide. The concept of the considered waveguide-cavity-atom system is illustrated in Fig. 1.

In practice, the prototypical waveguide-cavity-atom system considered here can be implemented within a variety of physical architectures. Those include, but not limited to photonic-crystal structures forming waveguides and cavities for photons, which interact with natural atoms birn (); tiecke (), individual semiconductor quantum dots (QDs) lohd2015 (), or embedded color centers engl2010 (); radu2017 (); tapered optical fibers coupled to whispering-gallery-mode resonators interacting with cold atoms trapped near their surface dayan (); superconducting circuit QED (cQED) setups, where nonlinear properties of the Josephson junctions (JJs) are utilized to embody atom-like systems clarke (); schm (); nori (). Remarkably, but general properties of those diverse physical systems can be described by the same mathematical model discussed below.

### The Hamiltonian

The considered waveguide-cavity-atom system is described by the Hamiltonian shen (); reph (); oeh (); chu2011 ():

(1) |

Henceforth, we set . Thus, throughout the paper, all energies are given in frequency units.

The first three terms on the right-hand side (rhs) of Eq. (1) describe the interacting cavity-atom system and constitute the Hamiltonian of the JC model . The first term in describes a noninteracting single-mode cavity with a resonance frequency , where () is a bosonic annihilation (creation) operator of a photon in the cavity mode. The second term in the Hamiltonian describes an individual 2LA with being an energy gap between its ground and excited states. We have introduced the atomic raising and lowering operators connected to the Pauli operators as , , and . The Hamiltonian of the atom-cavity interaction, represented by the third term in the rhs of Eq. (1), is given in the rotating-wave approximation (RWA), which assumes the conditions

(2) |

are satisfied, where denotes a detuning between frequencies of an atomic transition and a cavity resonance. Parameter stands for an atom-cavity coupling strength. The criteria (2) are satisfied for a broad spectrum of experimental wQED realizations ranging from microwave superconducting circuits to photonic-crystal structures operating at optical frequencies. Nonetheless, one should note recent progress in the development of various system (see, e.g., reviews in Refs. pforndiaz (); kockum () and references therein) exhibiting a light-matter coupling strong enough to break down the RWA applicability and switch the system into the ultrastrong coupling regime.

The term in Eq. (1) describes a continuum of independent bosonic modes which models a waveguide. Operator () annihilates (creates) a photon with energy in the waveguide and obeys the equal-time commutation relation . The waveguide-cavity interaction is represented by the term , where is a waveguide-cavity coupling strength. Here we have neglected the actual frequency dependence of the cavity-waveguide coupling assuming that the frequencies of the cavity resonance and the incident photons are close to each other and much larger than the spectral bandwidth of the wave packet. This is expressed by the criteria

(3) |

where is a central frequency of the incident pulse with a spectral bandwidth . Due to the second criterion in (3) we extend boundaries of integration over the photon frequencies to . This approximation gives negligible deviation from the exact result due to the sharp localization of the wave packet spectrum in the vicinity of , which is expressed by the second criterion in (3).

The model Hamiltonian contains no terms describing the interaction of the system with external reservoirs leading to dissipation, which implies that in the paper we are focused exclusively on a unitary evolution of the system. This approximation is legitimate when the parameters of the system meet the constraint , which indicates that dissipation processes should be the slowest ones. It is also assumed that the temperature at which the system operates satisfies the criterion , where is the Boltzmann constant. This condition ensures that the average number of thermal excitations in the system is negligible . The above conditions can be realized in the modern cQED systems (see Ref. schm () and references therein for the parameters of state-of-the-art superconducting cQED setups).

One can show that the operator of the total number of excitation in the system is an integral of motion since it commutes with the Hamiltonian: , which indicates that the number of excitations in the system is conserved.

## Iii Single-photon wave packet

Let us first consider the case of the JC system driven by a single-photon wave packet. As long as we set that initially the JC system contains no excitations (the atom resides in its ground state , and the cavity field is in the vacuum state), the evolution of the system occurs entirely within the one-excitation domain of the Hilbert space of the system states since the number of excitations in the system is conserved. The time-dependent wavefunction of the system has the form

(4) |

where stands for a vacuum state of the entire system – a state with no excitations in both the waveguide and the JC system . In Eq. (4) is a single-photon SDF, stands for an amplitude of a state with the ground-state atom and the cavity field in the single-photon state, and is an amplitude of a state with the vacuum field in the cavity and the excited-state atom. Here and through the rest of the main part of the paper, we work within the Schrödinger picture with time-independent operators.

The initial (at ) state of the system is given by

(5) |

where is a state of the ingoing photon. It is given by

(6) |

with being an ingoing wave packet SDF satisfying a normalization condition .

Evolution of the single-excitation amplitudes entering Eq. (4) is governed by a set of coupled equations as follows

(7a) | |||

(7b) | |||

(7c) |

where and . The initial conditions are the following: and . The method of derivation of the above set of equations is presented in Appendix A. The exact analytical solutions of the system (7) can be obtained using the Laplace transform method (see the details in Appendix B).

Since the ingoing pulse has a finite duration and the photons inside the cavity mode have limited lifetime , in the long-time limit

(8) |

the system reaches the steady state: the cavity field is in the vacuum state, the atom resides in its ground state, while the scattered photon propagates in the waveguide as a free excitation. Applying the condition (8) to solutions of Eq. (7) one obtains the final state of the system , where

(9) |

is the state of the scattered photon. Here stands for a single-photon frequency-dependent phase shift induced by the JC system. It is determined as shen (); reph (); oeh ():

(10) |

The complex single-photon resonances of the open (i.e. coupled to the waveguide) JC system are given by:

(11) |

where . In the resonant case () one has and . For the strong-coupling regime of the JC system () one can employ an approximation .

## Iv Two-photon problem: Quantum-state evolution

### iv.1 The wavefunction

Here we consider the evolution of the waveguide-JC system when the incident wave packet is in the two-photon Fock state. Due to conservation of the excitation number the dynamics of the system is restricted exclusively by the two-excitation domain of the Hilbert space of the system states. The wavefunction of the system at arbitrary moment is represented by a superposition of two-excitation states:

(12) |

The first term in Eq. (12) describes a state with both photons in the waveguide. The time-dependent two-photon SDF exhibits symmetry to permutation of photon frequencies because of the bosonic nature of photons. A quantity is referred to as a two-photon joint spectrum (2PJS) and determines a distribution of a joint probability of finding a pair of photons with frequencies and in the waveguide. The probability amplitudes correspond to the states of the system with the waveguide and the JC system both containing a single excitation. The amplitudes correspond to the states with the waveguide void of photons and the JC system containing two excitations. Subscripts and specify the state of the atom (ground or excited).

The initial state of the system in the case of the two-photon drive is with the wavefunction of the ingoing photons given by

(13) |

In the above expression stands for a two-photon SDF of the incident wave packet obeying the condition .

It is assumed that the incident two-photon wave packet is composed of a pair of indistinguishable photons with SDFs given by . In this case, factorizes into the product of single-photon SDFs as

(14) |

Therefore, the initial photonic state is separable and can be expressed in the form

(15) |

Thus, following the arguments by Rohde et. al roh (), here we deal with an ingoing wave packet in the two-photon multimode Fock state. In general, multimode Fock states are a special case of a broader class of multimode multiphoton states.

### iv.2 Evolution equations

Equation of motion governing the evolution of the two-photon SDF reads as follows

(16) |

with the initial condition . The rest of the probability amplitudes in Eq. (12) have zero initial conditions and obey the evolution equations as follows:

(17a) | |||

(17b) | |||

(17c) | |||

(17d) |

The details of derivation of Eqs. (16) and (17) are given in Appendix A. We solve the above system of ODEs governing evolution of the amplitudes both numerically and analytically. For numerical solution we use NDSolve function of the Mathematica system. The analytical solutions are obtained using the Laplace transform method (see derivations in Appendix B).

We proceed to investigation of evolution of the 2PJS of the wave packet in the course of its interaction with the JC system. For calculations we model the single-photon SDF of the ingoing wave packet by the Lorentzian function given by

(18) |

where parameter denotes a moment of time when the incident wave packet arrives at the cavity site. Hereinafter we set . For modeled by a Lorentzian shape determined by Eq. (18) one obtains:

(19) |

where is the Heaviside function.

The animation of 2PJS evolution can be found in Supplemental Material suppl (), while the series of snapshots for the specific moments of time and different parameters of the system are demonstrated in Fig. 2. These computations reveal that the 2PJS of the wave packet experiences a transformation in the course of its interaction with the JC system. Another important detail is that for , which corresponds to the setup with the cavity either empty or decoupled from the atom, the 2PJS ultimately regains its initial shape for times satisfying the condition (8). This result will be rigorously explained below in the paper. On the contrary, for the case of non-zero interaction () between the cavity and the atom, the 2PJS of the scattered wave packets [for times satifying the criterion (8)] are transformed compared to the initial 2PJS. As it will be demonstrated later in the paper (see Sec. VIII), unlike the ingoing photons, which are not entangled [see Eq. (13)], for the scattered photons are in the frequency-entangled state.

## V JC system dynamics

For the analysis of the excitation dynamics of the JC system we use Eq. (12). The average photon number in the cavity is given by

(20) |

The atomic excited state population is determined as

(21) |

Moreover, using the exact wavefunction (12) one can calculate the probability to find the JC system in the state with a given number of excitations (). The probabilities of finding one and two excitations in the JC system at an instant are determined as

(22) |

and

(23) |

Calculations of the cavity and atomic excitation dynamics presented in Fig. 3 show that shorter ingoing pulses (larger ratio ) results in more efficient excitation of both the cavity and the 2LA. The cavity and atomic populations exhibit the oscillations arising due to exchange of excitations between the 2LA and the cavity. For times exceeding the duration of the ingoing pulse the cavity and atomic populations decay exponentially as .

Computations of and shown in Fig. 4 reveal that the probability of finding the JC system in the two-excitation state is considerably reduced compared to the empty cavity () case (see Figs. 4c and 4d). This is explained as follows. Coupling the cavity mode to the 2LA, which is a saturable (i.e. nonlinear) system due to its ability to emit or absorb only one photon at a moment, introduces nonlinearity into the composite JC system. This nonlinearity gives rise to the photon blockade effect birn (); henn () when the presence of one excitation in the system prevents the appearance of the second one.

Since the probability of the JC system to reside in the two-excitation state is inhibited due to the atom-induced photon blockade, the frequency of oscillations of the cavity and 2LA populations is approximately the vacuum Rabi frequency of the JC system . A slight mismatch between these frequencies is due to the contribution of two-excitation processes and cavity-waveguide coupling.

## Vi Scattered two-photon state

Besides the transient quantum-state dynamics of the system, the asymptotic (long-time) solutions are of considerable interest as well. For times satisfying the condition (8), one arrives at the result that only the term which corresponds to the free evolution of the scattered two-photon field in the waveguide survives in the wavefunction of the system expressed by Eq. (12). The remaining terms vanish indicating that the JC system is completely depopulated and rests in its ground state . The final state of the system is then given by , where stands for the scattered photons wavefunction. It is expressed as

(24) |

where is the two-photon SDF of the scattered wave packet given by (see derivation in Appendix B):

(25) |

where has the form

(26) |

The result expressed by Eq. (25) can be equivalently derived using the two-photon S matrix of the JC system obtained, e.g., in Refs. shi2011 (); reph (); oeh ().

The term in Eq. (25) corresponds to the uncorrelated elastic scattering of two photons when the energy of each one is conserved. The term describes inelastic scattering when the energies of individual photons are not conserved. The complex two-photon resonances of the open JC system are given by

(27) |

### Cavity decoupled from atom

Let us consider a limiting case of the JC system when the cavity is empty or decoupled from the atom (). In this case, the inelastic term vanishes in Eq. (25), and the two-photon SDF of the scattered wave packet reduces to

(28) |

where is the SDF of the scattered photon with the single-photon phase shift given by . Thereby, the photons acquire only phase shifts in the course of scattering from the cavity which does not contain an atom. Using Eq. (24) one has

(29) |

where . The above expression indicates that for the state of the scattered photons is separable.

## Vii Spectrum of scattered photons

Spectrum of the scattered photons is determined as sto2013 (); sto2014 (), and characterizes the density of photons with frequency in the scattered wave packet. The average number of photons in the frequency band is given by . Integration of over the entire frequency range gives the average scattered photon number. Since the number of excitation in the system is conserved, one has , where is the spectrum of the ingoing wave packet. Using Eq. (24) one arrives at the expression for as follows:

(31) |

It follows from Eqs. (31) and (30) that for the setup with the cavity decoupled from the atom one has implying that for the spectrum of the scattered photons is identical to that of the ingoing ones.

Figure 5 demonstrates the dependence of the scattered photons spectrum on the central frequency of the ingoing wave packet . Calculations reveal that the most pronounced modification of the scattered wave packet spectrum occurs when lies in the vicinity of the JC single-photon resonances . The effect of atom-cavity coupling strength on the scattered spectrum for is shown in Fig. 6. One can see that increased atom-cavity coupling leads to stronger modification of the scattered photons spectrum.

If one substitutes Eq. (25) into Eq. (31) the scattered photons spectrum can be represented as

(32) |

In the above expression stands for the contribution to the scattered photons spectrum arising purely from the inelastic scattering determined by

(33) |

The term

(34) |

arises due to interference of components of the outgoing wave packet produced by elastic and inelastic scattering from the JC system. Contributions of all these terms to the scattered photons spectrum for the different parameters of the system are shown in Fig. 7. Calculations reveal that takes negative values producing the “dip” in the shape of the scattered photons spectrum seen in Fig. 6.

## Viii Frequency entanglement

As it was rigorously shown in the Sec. VI when the two-photon wave packet interacts with the cavity decoupled from the atom the scattered photons acquire only phase shifts while their SDF remains factorable. The state of the scattered photons is separable as that of the ingoing ones. The situation changes starkly when one switches on the coupling between the cavity and the atom. The JC nonlinearity gives rise to the effective photon-photon interaction which results in strongly-correlated scattering (photon blockade) discussed in Sec. V. According to Ref. xu2013 (), the photon-photon interaction leads to frequency entanglement. In our setup the latter arises due to the presence of the “inelastic” term in the SDF of the scattered wave packet.

The frequency entanglement of the scattered photons can be quantified using the Schmidt decomposition of their SDF fed ():

(35) |

which represents the SDF of the outgoing two-photon wave packet as a weighted sum of products of single-photon SDFs. The decomposition coefficients satisfy . The Schmidt-mode single-photon SDFs and form a complete set of orthonormal functions: and . The Schmidt coefficients and SDFs are determined via solution of the eigenvalue problem par (); wal ():

(36) |

where the integral kernels and are given by

(37) |

Since the SDF obeys the permutation symmetry , this leads to which, in turn, gives . We tackle the eigenvalue problem (36) numerically. For this purpose we discretize into a uniform grid on a square domain of frequencies ranging around .

Using the decomposition (35) along with Eq. (24) one represents the outgoing photons state as a superposition of two-photon separable states:

(38) |

where .

If Schmidt coefficients are all zero except for one, this indicates that the two-photon SDF is factorable and the photons are not entangled. Otherwise, we deal with the frequency-entangled photon pair. As a measure of entanglement, we also employ the von Neumann (entanglement) entropy . For entangled states , while zero entanglement entropy indicates that the quantum state is separable. The von Neumann entropy is expressed in terms of Schmidt coefficients as benn ()

(39) |

It follows from Eq. (28) that for the setup with the empty cavity () the SDF of the outgoing photons is factorable and there is only one nonzero Schmidt coefficient indicating that there is no frequency entanglement. On the contrary, for the cavity interacting with the atom () the outgoing photons exhibit frequency entanglement since and as presented in Fig. 8.

## Ix Summary

In this section, we summarize our results. We studied the scattering of two-photon Fock-state wave packets on the JC system employing a wavefunction approach. Using the interrelation between the Heisenberg and Schrödinger representations (see Appendix A) we derived sets of coupled equations of motions governing the evolution of the probability amplitudes which determine the state of the system at the arbitrary moment of time. Using the exact solutions of these equations of motion we tracked the evolution of a 2PJS of a wave packet in the process of its interaction with the JC system. The excitation dynamics of the latter was studies as well. We determined the probabilities to find the JC system containing one and two excitations at the arbitrary instant. We found that the probability of finding two excitations in the JC system is inhibited compared to the case of the empty cavity which indicates the photon blockade induced by the 2LA nonlinearity.

In the long-time limit, when the system reaches its steady state, we derived the exact expressions for the SDF of the scattered two-photon wave packet. This result agrees with that obtained using the exact S matrix of the JC system derived in the literature. We also analyzed the scattered photons spectrum. Our calculations reveal that it differs significantly from the incident photons spectrum when the central frequency of the ingoing two-photon wave packet is close to one of the JC single-photon resonances.

We employed Schmidt decomposition and the von Neumann entropy as a measure of entanglement of the outgoing photons. We found that the photons scattered from the JC constitute a frequency-entangled photon pair or a biphoton. The waveguide-JC system can be exploited as a tunable deterministic source of itinerant frequency-entangled biphotons. This source can be built on the state-of-the-art superconducting cQED architecture.

In the paper, we focus on the case of the ingoing photons without spectral entanglement, thereby, investigation of the local quantum system response to the frequency-entangled few-photon states is of interest. The approach employed in the paper can be used for studying more complex quantum systems, for example, multilevel emitters, coupled-resonator arrays, optomechanical systems. Moreover, one can go beyond the RWA employed in the paper for the description of atom-cavity coupling and consider multi-photon scattering on a quantum system operating in the regime of ultrastrong light-matter interaction. Consideration of those systems and regimes constitute possible directions for the follow-up studies.

###### Acknowledgements.

The author thanks O. O. Chumak for reading of the manuscript and useful comments, and A. Sokolov and S. Lukyanets for fruitful discussions.## Appendix A Derivation of evolution equations for the probability amplitudes

Here we present derivations of the equations of motion which govern the evolution of the single- and two-excitation probability amplitudes. For that purpose, in this Appendix we work within the Heisenberg picture and start with the derivation of the evolution equations for the operators introduced earlier in the paper. The model Hamiltonian given by Eq. (1) generates the equation of motion for the waveguide variable as follows

(40) |

which can be formally integrated as

(41) |

where stands for the annihilation operator of a photon propagating in the waveguide as a free excitation.

Equation of motion for the cavity photon annihilation operator reads as

(42) |

Eliminating the reservoir variables by substituting Eq. (41) into Eq. (42) and integrating over one arrives at the equation of the form

(43) |

where we have introduced an operator . Using Eq. (41) one can show that the commutation relations hold

(44) |

These commutators are routinely utilized below for derivation of equations of motion for the probability amplitudes.

The Heisenberg equation for the atomic lowering operator takes the form

(45) |

Equations of motion for the creation operators are derived by the Hermitian conjugation of the corresponding equations for the annihilation operators.

Using Eq. (4) one can formally express the single-excitation probability amplitudes in the form as follows:

(46) |

Employing the property of the vacuum state along with the relations and , where denotes some quantum-mechanical operator, one can transfer the time dependence from the wavefunction onto the operators which gives . By doing so, one obtains the following representation for the probability amplitudes:

(47) |

Using this representation along with the Heisenberg equations for the corresponding operators and taking into account that and one immediately arrives at the set of equations of motion (7).

Following the same receipt, we represent the two-excitation amplitudes entering the wavefunction (12) as

(48) |

Employing the Heisenberg equations (40), (43) and (45), the commutators (44) along with the properties and one obtains equations of motion for the two-excitation amplitudes (17).

This scheme can be directly applied to the -excitation problem. When the JC system is driven by the -photon Fock-state wave packet, one arrives at the set of coupled equations of motions, where -excitation amplitudes depend on -excitation amplitudes and so on down to -excitation amplitudes governed by Eq. (7). This hierarchy of evolution equations can be then solved numerically with the arbitrary precision using one of the available ODE solvers.

## Appendix B Solution of evolution equations using the Laplace transform

### b.1 Single-excitation amplitudes

The Laplace transform turns the system of ODEs (7) governing the evolution of the single-excitation amplitudes into a set of algebraic equations as follows

(49a) | |||

(49b) | |||

(49c) |

The inverse Laplace transform of the solutions of these equations gives

(50a) | |||

(50b) | |||

(50c) |

Taking the long-time limit (8) in Eqs. (50) and integrating over one obtains

(51) |

Employing the relation in the square brackets in the rhs of the expression for one obtains