Stationary discrete solitons in circuit QED

We demonstrate that stationary localized solutions (discrete solitons) exist in a one dimensional Bose-Hubbard lattices with gain and loss in the semiclassical regime. Stationary solutions, by defi- nition, are robust and do not demand for state preparation. Losses, unavoidable in experiments, are not a drawback, but a necessary ingredient for these modes to exist. The semiclassical calculations are complemented with their classical limit and dynamics based on a Gutzwiller Ansatz. We argue that circuit QED architectures are ideal platforms for realizing the physics developed here. Finally, within the input-output formalism, we explain how to experimentally access the different phases, including the solitons, of the chain.


I. INTRODUCTION
Realizations of quantum nonlinear media as ultracold atoms in optical lattices [1], ion-traps [2] or superconducting circuits [3,4] are interesting candidates for future quantum information processors. Apart from this challenging goal, they are also testbeds to explore new many body states of matter both in the classical and quantum regime [5]. Among others, solitons -localized and form preserving solutions -are a paradigmatic example of collective nonlinear solutions. In the so-called classical limit for the Bose-Hubbard (BH) model, the operators are replaced by their c-number average, obtaining the well known discrete nonlinear Schrödinger equation (DNLS) [6]. Stable exact and numerical localized solutions, discrete solitons, exist in different dimensions and topologies [7][8][9][10]. Quantum solitons have been hypothesized to exist in the BH model with and without dissipation. Theoretical predictions based on different approaches, Gutzwiller- [11] , truncated Wigner-approximations [12], Gaussian expansions [13] as well as density matrix renormalization group techniques [14] have found slowly decaying localized solutions. Therefore, none of those were stable solutions for the dynamics. Thus, quantum fluctuations seem to kill these topological solutions. Experimental realizations in the quantum realm are few. Bose-Einstein condensates confirmed the temporal existence of bright [15] and dark [16] localized modes [17]. For ions in optical traps, a proposal for the observation of solitons [18] was shortly after followed by the experimental observation [19]. Things may change if dissipation and gain coexist. In the classical limit yielding the dissipative driven DNLS (DD-DNLS) equation, localized solutions have been reported [20,21]. Furthermore, the DD-DNLS exhibits "spontaneous walking" solitons [22]; using nonlinear gain and * Corresponding author: naether@unizar.es dissipation, exact travelling discrete solitons exist as stable dynamical attractors [23]. Therefore, an open question remains in the literature. What about quantum solitons in nonlinear media with loss and gain? In our opinion, the combination of many body physics, dissipation and driving is interesting. It provides new phases to explore with non thermal but equilibrated states, as already demonstrated in the dissipative driven BH (DD-BH) model [24,25]. Besides, it establishes a links with man made realizations of lattice systems where dissipation can be an issue [5]. In the present context these novel phases could provide solitons. In this work, we will discuss the existence of stationary equations. Stationary solutions, if stable, are obtained via the dissipative dynamics no matter of the initial state (belonging to the basins of attraction). Therefore their preparation is easier and more robust. Along this work, we first argue that for the DD-BH model quantum solitons have no anti-continuous limit, i.e, the uncoupled lattice system has a unique stationary solution [26][27][28]. This is important, since the single site DD-DNLS for the same parameter regime can have different fixed points in their irreversible dynamics. This is a qualitative difference and should alert to take some care in the utilization of the DD-DNLS for finding solitons. To fix this problem, we include quantum fluctuations up to second order. In doing so, we recover the uniqueness in the stationary single-site solution. Moreover, solitons exists within this (we call it) semiclassical approximation. We discuss their stability and range of existence. We also describe other types of phases appearing when the soliton solution is unstable or absent. We complement our study with a Gutzwiller-Ansatz showing that localization is more persistent within the Gutzwiller in the range where solitons exist within the semiclassical limit. Finally, we discuss a physical realization for the DD-BH based on a circuit-QED architecture [29]. The physical support for our model is complemented with a proposal for a measurement scheme based on an input-output theory, to access the different phases using (already demonstrated) experimental capabilities. The rest of the paper is organized as follows. The next section II presents the model including dissipation and gain. Besides, we briefly describe the different theoretical approaches used: a second order (in the quantum fluctuations) expansion (SOE) and the Gutzwiller-Ansatz. We finish the section with a possible implementation in circuit QED architectures. In Sect. III we show our numerical results in both approximation schemes and we compare them against the classical DD-DNLS. We present in Sect. IV the input-output formalism for measuring the different phases and conclude with some discussion in. V.
(1) It marks a minimal model for interacting bosons in a lattice. The model (in the rotating frame of the drive) is characterized via δω (the detuning of the bare resonator frequency ω 0 from the pump frequency ω d , δω = |ω 0 − ω d |), A (the driving amplitude for this coherent external driving), U (the onsite repulsion) and J (the strength of the hopping among sites). Phenomenologically, single-particle-losses can be casted in a Gorini-Kossakowski-Sudarshan-Lindblad master equation [30] with γ −1 the time scale for the losses and {A, B} = AB + BA the anticonmutator. A pictorial and physical realization based on circuit QED is shown in fig. 1. Without loss and driving, the BH is a cornerstone in many body physics. The generalized BH, Eqs. (1) and (2) mark, then, a paradigmatic model for the study of collective phenomena with driving and dissipation. The dynamical equations for the averages a † n ...a m ≡ Tr(a † n ...a m ) are given by, i d a n dt = (δω − i γ 2 ) a n + A − J( a n+1 + a n−1 ) (3a) + U (2 a † n a n a n + a 2 n a † n − 2 a n 2 a † n ) i d a † l a n dt = − iγ a † l a n + A( a † l − a n ) + J( a † l−1 a n + a † l+1 a n − a † l a n−1 − a † l a n+1 ) (3b) + U ( a † l a † n a n a n − a † l a † l a l a n ) i d a l a n dt = (2δω − iγ) a l a n + A( a l + a n ) − J( a l−1 a n + a l+1 a n + a l a n−1 + a l a n+1 ) (3c) + U ( a † l a l a l a n + a † n a n a n a l + δ n,l a n a n ) . . .
The dots above indicate that, due to the interaction term U a †2 n a 2 n , an endless hierarchy of equations for the n-point correlators a † n ...a m is obtained. Therefore, the set needs to be cut at some order.

A. Zeroth order: The DNLS equation
The simplest approximation is the so-called classical limit, consisting in replacing operators by their averages: a †2 n a n → |ϕ n | 2 ϕ n , with ϕ n = a n . The approximation can be understood as the zeroth order cumulant expansion in the quantum fluctuations. In doing so, the equation for the first moments (3a) forms already a closed set. The resulting equations are the celebrated DNLS equations, in this case, with driving and dissipation: (4) For this set of equations [20,22], apart from dark soliton and kinks, there also exist bright solitons for a defocusing nonlinearity U = −1, on which we will focus along this work. The main question that we tackle is, if the solutions found in the DNLS survive the inclusion of quantum fluctuations.

B. Second order expansion
Let us consider now equations (3) up to second order of correlations. In doing so, we rewritê a n =â n − a n + a n =: δâ n + a n .
Neglecting terms with O(δâ 3 ) 0, any n-correlator can be written in terms of 2-point correlators: Above j = k = l = m and A j can be any annihilation a n (creation a † n ) operator. Consequently, equations (3a), (3b) and (3c) form a closed set which stands for a second order expansion (SOE). Its numerical solution provides us the results in the subsection III B. It is worth emphasizing that, instead of the SOE, Gaussian expansions as the Hartree-Fock-Bogoliubov (HFB) or higher order terms might also be considered [13,31,32]. However, in the parameter regime explored, our SOE approaches better the exact result for the single site case [26].

C. Gutzwiller Ansatz
To consider higher on-site correlations, we will compare the results of SOE with the time evolution of a density matrix using a Gutzwiller Ansatz [11,32] (which assumes a factorized form for the density matrix): with site-dependent density matrices ρ n . Using that tr( n ) = 1 (tr(∂ t n ) = 0), we obtain the quantum nonlinear master equation set: Though several systems may be modeled by means of a BH model with losses and external driving as in Eqs. (1) and (2), we fix our attention on circuit QED architectures [24]. The latter seem to be an ideal platform to study many body physics [5]. In this subsection we argue that the fundamental blocks for simulating (1) have been already experimentally demonstrated. The first ingredient is having nonlinear resonators. For that, we think about recent experiments where coplanar waveguide resonators are interrupted by a Josephson Junction (JJ) [Cf. Fig. 1c]. The JJ provides the nonlinearity through the term E J cos(2π/Φ 0 δφ) in the effective action. Here E J is the Josephson energy, Φ 0 the flux quanta and δφ the jump in the flux at both sides of the junction [33,34]. In Ref. 35, the authors measured nonlinear resonators that can be modeled within the Hamiltonian (after expansion of the cosine): with ω 0 ∼ = 6GHz and U ∼ = −700KHz. Therefore, by choosing pumps with driving frequencies detuned from the ω 0 on the KHz − MHz regime different U/δω in (1) can be simulated. Finally, higher order terms can be safely discarded, U /U ∼ = 10 −3 . An intercavity coupling as in (1): has been already measured in a wide range of values for J, even reaching values of J/ω 0 ∼ = 0.2 [36]. Moreover, a tunable coupling J has been measured [37]. Highly reproducibility in the resonator bare frequencies, a necessary ingredient for building many body arrays, has been also achieved [38]. Finally, a measurement scheme is mandatory. Here, we rely on the field tomography techniques developed in the circuit QED community [39][40][41]. As explained in section IV, measuring field-field correlators is sufficient for accessing the different phases of (1), including the solitonic solutions. A possible architecture is depicted in figure 1 c). Inspired on Ref. 29 we envision a one-dimensional array of nonlinear cavities: superconducting resonators interrupted by a JJ. The design is such, that the coupling can be tuned by locally approaching the resonators [36]. The measurement can be accomplished by an auxiliar transmission line that couples the array and where an input field is impinged and the ouput is measured as we will explain in sec. IV. Therefore, the simulation and measurement may be possible within the technological state of the art.

III. RESULTS
We summarize here our numerical findings. We first review the classical DNLS limit. Then, we report on our quantum results, both for the SOE and Gutzwiller-Ansatz.

A. DNLS: solitons with anti-continuous limit
When considering the DNLS one common procedure for finding localized modes is as follows. Imagine that two stable solutions exist for the single site case, say ϕ and ϕ . Then, at zero hopping (J = 0), the solution ϕ m = ϕ , ϕ n = ϕ ∀n = m is a stable localized solution. This example corresponds to a soliton with a localized amplitude at one site (m). This (trivial) soliton can be used as starting point to find solutions by turning on the hopping, J = 0.
In the nonlinear jargon, the zero hopping case is named as anti-continuous limit and [6,7,42] a variety of numerical continuation techniques from zero to non-zero hopping have been used as, for example, the Newton-Raphson method. This procedure is not restricted to conservative settings, but also works very well in the driven dissipative classical models mentioned before. So was used in [20,22] to characterize the localized modes for the DD-DNLS. The results of this numerical continuation for the DD-DNLS (4) will be compared to SOE in the following.

B. SOE: Localization without anti-continous limit
In the quantum regime, the continuation from the zero hopping (anti-continuous limit) can not be used. The reason is, that the single site of (1) and (2) has a unique solution [26]. Therefore, if quantum solitons exists they do not have an anti-continuous limit. Without the possibility of finding solitons by continuation, an educated guess is to try the search within the parameter regime where they exist in the classical DNLS limit [20,22]. As anticipated, through this work we use the approximate treatment SOE. Therefore, for consistency, we must check that SOE has a unique solution in the single site case. In figure 2 we delimit this consistency region. Darkgray areas denote multivalued SOE solutions of (3) and therefore forbidden regions. For comparison we also draw (in light-gray) the forbidden space for the DD-DNLS. The red star marks a point with a unique solution for the SOE but with proven solutions in the DNLS. It seems a good point to start with. The concrete parameters are U = −1 A = γ = 2 and -for a better comparison with the DNLS [20,22]-detuning of δω = 3 + 2J. The only free parameter remaining is the coupling J. Those parameters are used along the text.
In general, steady-state solutions can be obtained by simply integrating the dynamics for (3) up to sufficiently long times such that a stationary dynamics is reached [43]. By construction, only stable solutions are found. Besides, there is no necessity of fine-tuned state preparation. Finally, this method provides not only steady-state solutions, but also time-periodic modes. Another possibility, which also finds unstable modes, is using the corresponding algebraic set of equations for d t . . . = 0 with e.g. a Newton-Raphson scheme. Unfortunately, this method has no guarantee to converge and could be used successfully only for specific parameters (see the unstable soliton mentioned below). We use the long-time dynamics for all stable modes presented here. Examples of different solutions assuming periodic boundary conditions are plotted in fig. 3. In figure 3a the dynamics for a stationary ripple mode is depicted. In Fig. 3b) and c) examples for the time evolution of stationary and periodic localized modes are plotted respectively. To visualize the physical mechanism yielding these solutions we choose to plot the mean amplitudes 0.5 (| a n (t min ) | + | a n (t max ) |) vs. sites n and J in Fig.  3d. This averaging is recommendable to better illustrate the localized character of the oscillatory mode, in all other regions of steady-state modes it does not change the picture. The figure shows that, for vanishing and small J, the homogeneous mode is the only stable solution. It becomes unstable at J 1 0.1, a symmetrybreaking bifurcation not present in the DNLS limit. For small, but finite J > 0.1 the ripple modes with one site having a higher amplitude than its two neighbors dominate the dynamics. For the spatial periodicity of these modes there is a certain dependency on the number of sites, as can be seen in fig. 3(d) with n = 25 leading to a defect in the right bottom corner, whereas for n = 30 in fig. 3 (b) no such defect can be found. As the coupling increases, repeated bifurcations into modes with different periodicity can be observed as the extension of the maxima grows and less peaks can be accommodated within the lattice. Finally, this leads to only one central localized mode in fig. 3(d) at J 2 appearing dynamically as the steady state. Increasing J, the stationary soliton starts to oscillate at J 3 and finally relaxes to the homogeneous mode for J 4 5.38. Thus, the SOE does exhibits a train of symmetry breaking bifurcations towards more and more localization, a behavior not found in the DD-DNLS.
We have seen that the classical DNLS is a particular limit of the quantum model in which the quantum fluctuations are neglected. If the quantum corrections were negligible, both the SOE and the DNLS would produce similar results. As shown in Fig. 3(e) this is not the case. There we plot the dependence of the center site amplitude | a c (t evol ) | on J after the dynamics settled into a steady-state or periodic state of SOE ((3) (red) within the evolution time t evol = 100. Please note, that all possible phases are shown in 3(e) and the value of | a c | does not necessarily indicate that there is a difference to neighbouring sites. The arrows point to the bifurcation points J i . When the dynamics is determined by the periodic mode we also show the maximum and minimum of | a c (t) | in green and blue. The algebraic stationary and unstable localized mode is shown with a blue dashed line. For comparison we show the results of the classical limit (4) in black, also exhibiting a periodic mode, but for higher J. The classical amplitude is nearly twice as high as the SOE value, but the main difference is in the classes of solutions found. Whereas the soliton mode is stable in the classical DNLS limit from J = 0 up to the appearance of the periodic solution at J cl , for the SOE limit there is no anti-continuous limit. At J 1 the ripples appear and persist for J 1 < J < J 2 . The bifurcation into the periodic solution is located at J 3 < J cl ; as well as the high-coupling homogenous mode at J 4 . The symbols in the SOE families denote the examples shown in 3(a)-(c).

C. Gutzwiller Ansatz
We complement the SOE with the dynamics within the Gutzwiller Ansatz (7). As initial conditions, we use the homogeneous steady state for the corresponding value of J at all sites but the center, where we assume a coherent state with higher n c (t = 0) 8.8. The Fock-space per lattice site ρ i is truncated to a maximum of 15 excitations in a lattice of 15 sites. Within this parameter space, we were able to find a region at J 2, where the initially localized distribution survives much longer. Even though we can compare the SOE and the Gutzwiller Ansatz only qualitatively, this corresponds to the regime of stable soliton modes in the SOE limit. In fig. 4(a) Fig. 4(c) shows, that the initial value decays abruptly to n c 3, from that point on the decay becomes very slow. This indicate the existence of a weakly unstable localized mode. For the example shown in fig. 4(b), the initial decay is equally fast, furthermore some oscillations can be observed, which is a reminiscence to the existence of the periodic mode in the SOE limit for these parameter values. Additionally, hints of these oscillations can be seen in ∆n c,1 for J 3.5 . We could not find any indications of an ripples modes, since it presents higher-order inter-site correlations neglected within the Gutzwiller-Ansatz, as we will show in the following section.

IV. IN-OUT MECHANISM
It still remains the question of how to extract the information stored in our discrete array of cavities. In order to do this, we will follow the input-output formalism [44]. The basic idea here is to make our dissipative system interact with the electromagnetic field (EM), in the form of a transmission line (TL) [Cf. figure 1 b) and c)]. The EM field can be decomposed into two contributions: the input, or the radiation impinging onto the system, and the output, the sum of a reflected plus a radiated component [See Fig.1(b)]. The latter is determined by the system and its interaction with the TL. Therefore, measuring the output (and comparing it with the input) we can infer information on the system dynamics. We will briefly sketch the main steps to derive the inputoutput relations. Details on the calculations can be found in Appendix A (see also the Supplementary Material of Ref. [45]). The system we want to acquire information from is an open system described by the Hamiltonian (1) and the master equation (2 figure 1 c). The total Hamiltonian, including the TL and the interaction of our system with it, reads H = H open +H TL +H int . The TL interacts directly with the cavities and it can be viewed as a one dimensional EM field. In second quantization, the EM field can be described as a collection of harmonic oscillators. As it is depicted in fig. 1 b), we should consider carefully the direction of propagation of excitations in the TL. For this task we will introduce the EM field operators: l(p) (l † (p)) which destroys (creates) a photon with momentum p propagating to the left and r(p) (r † (p)) which destroys (creates) a photon with momentum p propagating to the right. In terms of these and for a linear dispersion relation (ω p = v|p|), the Hamiltonian of the EM field in the TL reads (11) with v the speed of light in the TL. Finally, the interaction H int considers the most general type of coupling in a solid-state device. It consists of an inductive part (flux interaction) and a capacitive contribution (charge interaction). The interactions are weak and point-like, happening at the position of every cavity, yielding where we have introduced a plane wave expansion for the cavity operators a n (t) = 1 √ N k e ikxn/d a k (t) (13) d being the lattice spacing of our cavity array with N sites. We are now able to solve the Heisenberg equations for the left and right operators. The idea is to relate the contributions of all momenta before the interaction (from an initial time t 0 → −∞) and after the interaction (up to a time t 1 → ∞). As it is shown in appendix A, for every cavity momentum k we only have significant contributions from those momenta p in a narrow region around the former. For the right operators, we call this momentum interval +k and for the left ones −k [cf. eq. (A10)]. Taking this into account, we introduce the (right) input operator r +k in (t) = 1 √ 2π +k dp e −ivp(t−t0) r 0 (p) (14) for times t > t 0 (r 0 denotes the r operator at t = t 0 ) and the (right) output operator r +k out (t) = 1 √ 2π +k dp e −ivp(t−t1) r 1 (p) for times t < t 1 (r 1 denotes the r operator at t = t 1 ). And similarly for the l operators. Thus, the Heisenberg equations lead to the following input-output relations with Γ a constant characterizing the strength of the TLnonlinear cavity coupling, and a k defined as In the presence of the TL the equations of motion for the a k (t) operators differ from those obtained from (2). The TL plays now the role of a second environment for the system described by (1). However, under the same approximations which led to (2) (Markov approximation), we immediately see that the role of the TL is to renormalize the decay rates γ. Namely, to add a contribution proportional to g 2 to them. In addition, the input field will renormalize the driving field amplitude (with strength of order g). Within relations (16) and (17) we can map output fieldfield correlations to cavity modes correlations. For example, using relation (16) r +k † out (t 1 )r +k out (t 2 ) gives us information on the system two-point correlator a † k a k provided that we only impinge a signal from the left Similar relations hold for other two-point correlations. Therefore, two time correlations in the output field map to the system two point correlations in momentum space. The first are experimentally accessible, the latter as we will show now provide sufficient information for distinguishing the different phases described throughout this work. We introduce the connected correlator, It is worth emphasizing that the above is identically zero in the classical limit. In fig. 5 we plot f a k , a † k for different stationary SOE solutions. The homogeneous solution (J = 0) is shown 5 a). The ripples (J = 0.3) are drawn in 5 b) with the same parameters as in fig. 3 a). Localized solutions, static and oscillatory are given in 5 c) and d) respectively. Their real space counterparts are given in 3 b) and c). The momentum space for each phase are clearly distinguishable from each other and show, that the input-output mechanism presents a perfect measurement scheme to prove the existence and stability of localized modes.

V. DISCUSSION
The phases for the Bose-Hubbard model with gain and loss have been investigated within a semi-classical approach. It has been argued that, in the zero hopping case, the unique solution is the homogeneous one. A translational-invariant broken symmetry solution (ripples) appears when the hopping term reaches some critical value. Increasing the hopping, the extension of the ripples grows and their periodicity decreases in a second bifurcation. Eventually, the discrete periodicity disappears and only one maximum remains, the stable localized mode. Passing from static to periodic (in time), this mode finally becomes unstable at higher J, transiting to a homogenous solution. Those successive symmetry-breaking transitions (from homogenous, to discrete periodic, to localized static and periodic modes, and back to homogenous) mark novel phases without counterpart in the Hamiltonian limit (zero dissipation, zero gain) of the Bose-Hubbard where the well known Mott-superfluid transition has been largely described. Apart from the interest in finding novel matter phases in the many body phase diagram, artificial systems with driving and dissipation present also a natural way of observing localization. Our calculations rely on a semiclassical approximation. We have complemented them with a Gutzwiller Ansatz where the on-site dynamics is expected to be more accurate but inter-site correlations are poorly described. Nevertheless, the regions where solitons exists in the semiclassical regime present long lived localized solutions within the Gutzwiller-Ansatz. Therefore, we expect that the semiclassical phases have some traces in the full, not yet explored, quantum dynamics.
The richness of phases presented here may be a motivation for future works considering the full quantum aspects of the model. In this line, our proposal within circuit QED presents a quantum simulator for going beyond the theory presented here.
where H open accounts for the open system: BH model + driving + environment. For simplicity, we will refer to the latter as our system. H TL describes the EM field propagating through the TL. In momentum space it is given by where we are assuming a linear dispersion relation ω p = v|p|. Our cQED proposal involves impinging a signal into the system and gathering information about it by means of the reflected and transmitted components of the former. For this task, it results helpful to decompose the EM field operators in: l † (p) (l(p)) which creates (annihilates) an excitation with momentum p propagating to the left with velocity v (the speed of light in the TL) and similarly, r † (p) (r(p)) which creates (annihilates) an excitation with momentum p propagating to the right with velocity v. Finally, H int is the interaction Hamiltonian Here we are considering a generic coupling between a system operator and the EM potentialÂ(x) In [45] following the lumped circuit element description of a TL, we decomposed the interaction into capacitive Rectangle approximation of (A8) (full line). This considers that the function is different from zero only in the region between its first two-zeros (symmetric around k/d) and that the height is constant and equal to N . and inductive contributions. The latter are encoded in the coupling function g n (x). We now introduce a plane wave expansion for the cavity operators a n a n (t) = 1 where d is the lattice spacing of our cavity array with N sites. Assuming a rotating wave approximation (RWA) regime we can rewrite (A3) as The TL only couples to the system at the position of the cavities, therefore, g n (x) = gδ(x − x n ) and we have where we have replaced p in the exponentials for ±ω p /v. We can approximate the sum by a rectangle of height N and width 2δ = 2π/N centered at k (being zero elsewhere) (see Fig. 6). Similarly, n e i(ωp/v+k/d)xn will be replaced by a rectangle centered in −k. Therefore, we can rewrite (A7) as where the integration intervals, following our previous approximation, are ± k ≡ [(±k − δ)/d, (±k + δ)/d] . (A10) We are now able to write the Heisenberg equation of motion for the EM field operators. Following (A1), (A2) and (A9) they arė r(p) = −ivp r(p) + vg N 2πω p a k (t) (A11) l(p) = +ivp l(p) + vg N 2πω p a k (t) (A12) Here we have included the explicit time dependence of the momentum cavity operators a k . Integrating (A11) from t 0 to t (t 0 < t) yields r(p, t) = e −ivp(t−t0) r 0 (p) + vg N 2πω p t t0 dt e −ivp(t−t ) a k (t ) (A13) where r 0 denotes the r operator at time t = t 0 . Notice that the former equations include a continuous momentum p and a discrete one k. Recall from our approximation to the sum (A8) that for any given k only the momenta p in a very narrow region of width δ (centered in k) contribute to our expressions. We now integrate (A13) over this momentum interval (+k for the right operators) dp r(p, t) = √ 2πr +k in (t) with the input operator r +k in defined as (following [44]) We have included the super index +k to stress that this operator sums the momentum contributions in a narrow band around k (positive for the right operators). The input operator takes into account the free evolution of all right-propagating EM field modes before the interaction with the system. Therefore, it acts as a driving field in the equations of motion of the cavity operators. In the continuum limit (δ → 0) and for t 0 → −∞, (A14) yields dp r(p, t) = √ 2π r +k as δ → 0, we make use of the following limit: lim z→0 sin(z)/z = 1 for z = δv(t − t )/d. This holds reasonably well for N 10 as it is shown in Fig. 7. In addition, we have introduced In a similar way, we can integrate (A11) from t to t 1 (t < t 1 ) and define a corresponding output operator.
r +k out (t) = 1 √ 2π +k dp e −ivp(t−t1) r 1 (p) (A18) We will find that the input and output operators are related by r +k out (t) = r +k in (t) + where we have chosen conveniently g = Γ/2π.