Nonlinear dynamical Casimir effect at weak nonstationarity

We show that even small nonlinearities significantly affect particle production in the dynamical Casimir effect at large evolution times. To that end, we derive the effective Hamiltonian and resum leading loop corrections to the particle flux in a massless scalar field theory with time-dependent Dirichlet boundary conditions and quartic self-interaction. To perform the resummation, we assume small deviations from the equilibrium and employ a kind of rotating wave approximation. Besides that, we consider a quantum circuit analog of the dynamical Casimir effect, which is also essentially nonlinear. In both cases, loop contributions to the number of created particles are comparable to the tree-level values.


Introduction
The dynamical Casimir effect (DCE) describes particle production by nonuniformly accelerated mirrors [1][2][3][4][5][6] and aptly illustrates many prominent features of the quantum field theory. On the one hand, this effect is similar to the celebrated Hawking [7][8][9] and Unruh [10][11][12] effects. Essentially, all these phenomena are related to a change in the system's quantum state due to the interactions with strong external fields; the only difference is the nature of the interactions and external fields. On the other hand, the DCE is amenable to a plethora of theoretical and experimental approaches, most of which are concisely described in reviews [13,14]. In particular, recently, this effect was experimentally measured in a system of superconducting quantum circuits [15][16][17]. The mentioned properties make the DCE a perfect model for the study of nonstationary quantum field theory.
Most of the research on the DCE, including experimental implementations, focuses on a simplified two-dimensional model with a free massless scalar field and perfectly reflecting boundary conditions (we set = c = 1): where mirror trajectories L(t) and R(t) are stationary in the asymptotic past and future, i.e., L(t) ≈ L ± + β ± t and R(t) ≈ L ± + Λ ± + β ± t as t → ±∞. Here, Λ ± > 0 is the distance between the mirrors, and β ± is their velocity in the corresponding limits. Note that for physically meaningful trajectories, L(t) < R(t) and |L(t)| < 1, |Ṙ(t)| < 1 for all t. Usually one also makes a Lorentz boost, shifts the origin and conveniently sets L(t < 0) = 0, R(t < 0) = Λ − without loss of generality.
Due to the nonstationarity of the model (1.1), its Hamiltonian cannot be diagonalized once and forever, so the notion of a particle is different in the limits t → ±∞. In other words, a natural mode decomposition of the free quantum field φ(t, x) is different in the asymptotic past and future: φ(t, x) = n a in n f in n (t, x) + H.c. , as t → −∞, n [a out n f out n (t, x) + H.c.] , as t → +∞, (1.2) where functions f in n , f out n solve the free equations of motion and form a complete basis with respect to a proper inner product. In general, the modes and creation/annihilation operators that diagonalize the free Hamiltonian in the asymptotic past (in-region) and future (out-region) are related by a generalized Bogoliubov (or canonical) transformation [18][19][20]: with nonzero Bogoliubov coefficients α nk and β nk . Hence, the number of the out-particles (f out n modes) created during the nonstationary evolution, in the Heisenberg picture, is given by the following expression (see Appendix A for an alternative derivation): N n = in|(a out n ) † a out n |in = k |β kn | 2 + k,l (α kn α * ln + β kn β * ln ) n kl + k,l α kn β ln κ kl + k,l α * kn β * ln κ * kl .
Here, we introduce the short notations for the level density n kl = in|(a in k ) † a in l |in and anomalous quantum average (correlated pair density) κ kl = in|a in k a in l |in and schematically denote the initial quantum state (possibly mixed) as |in .
In the Gaussian approximation, i.e., in the absence of interactions and nonlinearities, particle creation in the DCE is related only to mixing of positive-and negative-frequency modes. Indeed, assume that the initial quantum state of the system coincides with the vacuum, a in n |in = 0 for all n. In this case, both n kl = 0 and κ kl = 0, so the identity (1.4) significantly simplifies: (1.5) This identity also approximately holds at low temperatures: in this case, n kl is exponentially suppressed and κ kl is zero, so they give negligible contributions to the right-hand side of (1.4). Thus, on the tree level, particle creation is fully determined by the Bogoliubov coefficient β kn . For this reason, identity (1.5) is widely considered as the essence of the DCE, and the vast majority of literature -including seminal papers [1][2][3][4][5][6], reviews [13][14][15], and textbooks [18][19][20] -is purely dedicated to the calculation of N free n and derivative quantities. Nevertheless, real-world systems are rarely free from interactions and nonlinearities. In such systems, approximation (1.5) is not applicable even if the initial values of n kl and κ kl are small. In fact, in nonstationary situations, these quantities can receive nonnegligible loop corrections 1 : where U (t, t 0 ) denotes the evolution operator in the interaction picture, and t 0 , t set the moments when the interaction term is adiabatically switched on and off. In this case, identity (1.4) does not reduce to (1.5). Furthermore, loop corrections can be significant even if interactions are apparently weak. As was recently discovered, in many nonstationary interacting systems, n kl and/or κ kl receive secularly growing loop corrections, i.e., n (n) kl ∼ (λt) an and κ (n) kl ∼ (λt) bn with some constant integer powers a n > 0 and b n > 0 in every order of the perturbation theory in λ. Such corrections indefinitely grow in the limit t → ∞ and are not suppressed even if the coupling constant goes to zero, λ → 0. For instance, the secular growth of loop corrections was observed in the de Sitter [28][29][30][31][32][33] and Rindler [34] spaces, strong electric [35][36][37] and scalar [38][39][40] fields, gravitational collapse [41], and nonstationary quantum mechanics [42,43]. A short review of the origin and consequences of the secularly growing loops can be found in [44].
The DCE was also recently shown to possess secularly growing loop corrections [45,46]. Unfortunately, this analysis was restricted to the first few loops. At the same time, due to the secular growth, high and low order loop corrections in this model are comparable at large evolution times. Therefore, a definitive conclusion about the destiny of n kl and κ kl can be made only after a resummation of the leading contributions from all loops. In other words, the number of created particles (1.4) and other observable quantities in the nonlinear DCE have physical meaning only nonperturbatively.
In this paper, we estimate a nonperturbative contribution to the number of created particles in the nonlinear scalar DCE. For illustrative purposes, in most of the paper, we consider a simple nonlinear generalization of the model (1.1) with a quartic potential: restrict ourselves to small deviations from stationarity and assume that the initial state of the field coincides with the vacuum. This approximation simplifies the calculation of the Bogoliubov coefficients and allows us to resum the perturbative series of loop corrections to n kl and κ kl . As a result, we obtain that the nonperturbative loop contributions to (1.4) in the approximation mentioned above are comparable to the tree-level approximation (1.5) at large evolution timeseven when the coupling constant λ → 0. In addition, we extend these calculations to a more physically relevant model: which naturally emerges in a Josephson metamaterial implementation of the DCE [17] due to the small nonlinearity of superconducting quantum interference devices (see Appendix C for the details). For small deviations from stationarity, i.e., for relatively small variations of the light speed v(t, x), this model is qualitatively equivalent to a simplified model (1.8). In this case, the resummed loop contribution to the created particle number is also significant at large evolution times. This paper is organized as follows. In Sec. 2, we discuss the quantization of a free scalar field and calculate the Bogoliubov coefficients in the model (1.1). In Sec. 3, we derive the effective Hamiltonian of the nonlinear model (1.8). In Sec. 4, we calculate level density and anomalous quantum average in the assumption that almost all created particles populate the fundamental (lowest-energy) mode. In this case, the effective Hamiltonian is essentially quantum mechanical. In Sec. 5, we consider a large N model, which is qualitatively equivalent to the model (1.8), and confirm that at small deviations from stationarity, the fundamental mode is the most populous one. Thus we justify the assumption of Sec. 4. In Sec. 6, we repeat these calculations for a physically motivated model (1.9). Finally, we discuss the results and conclude in Sec. 7. Besides, we explain the physical meaning of N n in Appendix A, present the full expression for the interacting Hamiltonian of the model (1.8) in Appendix B and derive the model (1.9) from the Hamiltonian of a Josephson metamaterial in Appendix C.

Field quantization and Bogoliubov coefficients
Consider a quantization of the free model (1.1) with L(t < 0) = 0 and R(t < 0) = Λ − , i.e., initially resting mirrors. As usual, we expand the operator of the quantized field in the mode functions: Here, operators a in n and (a in n ) † satisfy the standard bosonic commutation relations, and mode functions f in n are expressed in terms of auxiliary functions G(z) and F (z): which solve the generalized Moore's equations [1]: and comply the initial conditions G(z ≤ Λ − ) = F (z ≤ 0) = z/Λ − . Note that due to stationarity of the motion of mirrors in the asymptotic future, functions G(z) and F (z) periodically grow as z → +∞: G(z + ∆z G ) = G(z) + 2 and F (z + ∆z F ) = F (z) + 2 with ∆z G = 2Λ + 1−β + and ∆z F = 2Λ + 1+β + . This property can be easily inferred from the geometric method for constructing modes [46][47][48][49]. Moreover, in this limit, the first Moore's equation implies a simple relation: Let us return to mode functions (2.2). First, they can be straightforwardly shown to solve the free equation of motion and form a complete orthonormal basis with respect to the Klein-Gordon inner product: (2.5) Second, these functions have definite positive frequency with respect to the Killing vector ξ = ∂ t : and diagonalize the free Hamiltonian at the past infinity: Due to these reasons, functions (2.2) are usually referred to as in-modes. At the same time, at the future infinity in-modes contain both positive and negative frequencies: where we change to the coordinatest = γ + (t − β + x),x = γ + (x − β + t − L + ) and denoteΛ = γ + Λ + , γ + = 1/ 1 − β 2 + for shortness. The coefficients in the decomposition (2.8) are nothing but the Bogoliubov coefficients determined by the inner products (2.5): of in-modes (2.2) and out-modes: For arbitrary mirror motion, functions G(z), F (z) and the Bogoliubov coefficients are very difficult to compute. Nevertheless, this computation significantly simplifies if we restrict ourselves to weak deviations from stationarity, i.e., consider trajectories of the form L(t) = l(t), R(t) = Λ − + r(t), 1, and exclude motions that induce constructive interference. In this case, G(z) and inverse function g(y) = G −1 (y) are approximately linear: G(z) ≈ 2 ∆z G z + O( ), g(y) ≈ ∆z G 2 y + O( ). Keeping in mind that these functions periodically grow with z or y, we expand their derivatives into the Fourier series: and find high order Fourier coefficients to be small, G n ∼ and g n ∼ for all n = 0. In addition, these coefficients are symmetric, G −n = G * n and g −n = g * n , since the original functions are real. These properties allow us to expand complex exponents from (2.2) into simple plane waves: Combining this identity with the relation (2.3), changing to the (t,x) coordinates and calculating inner products (2.9), we straightforwardly obtain the Bogoliubov coefficients (we single out the Kronecker delta δ n,k to make the -expansion explicit): (2.14) We emphasize that these approximate identities are valid only for relatively small frequencies, n 1/ , where particle wavelength is much larger than the characteristic mirror displacement, λ ∼Λ/πn Λ . We also expect that at high frequencies, interactions between the field and the mirror can be approximated by a quasistationary process (at least if we assume no constructive interference). Hence, the corrections to the stationary Bogoliubov coefficients (i.e., identity transformation) in this case are approximately zero, α n =k ≈ 0 and β nk ≈ 0 for n 1/ or k 1/ . We remind that identities (2.14) are valid only for weakly nonstationary motions. Notable examples of such a motion include "broken" trajectories with a small final velocity : 15) or resonant oscillations with destructive interference 2 , e.g.: (2.16) 2 We emphasize that this is a very special case of resonant oscillations; in fact, the interference is constructive for the majority of oscillation amplitudes, frequencies, and dephasing angles [50,51].

Effective Hamiltonian
The full Hamiltonian of the nonlinear DCE, model (1.8), is very complex. Therefore, we need a reasonable approximation to simplify the calculations and estimate the number of created particles. First, we assume small deviations from stationarity, i.e., β nk ∼ α n =k ∼ O( ) with 1. In the leading order in , this approximation forbids almost all processes that violate the energy conservation law. Second, we consider small coupling constants and large evolution times 4 , λ → 0, t → ∞, λΛ + t = const, i.e., keep only the leading secularly growing contributions to the evolution operator. This limit is similar to the rotating-wave approximation in quantum optics [62][63][64][65] and suppresses rapidly oscillating parts of the effective Hamiltonian. We also note that in this limit, the Bogoliubov coefficients are approximately constant. Both these approximations are inspired by a quantum mechanical toy model of nonstationary particle production [43,66,67].
First of all, let us calculate the free Hamiltonian in the suggested limit. Substituting mode decomposition (2.8) into identity (2.7), rotating to the (t,x) coordinates, integrating overx and neglecting rapidly oscillating and constant contributions, we obtain the following approximate identity: as t → +∞. Here, we denote the frequency of the f out n mode as ω n = πn/Λ. In the last line, we also substitute the approximate Bogoliubov coefficients (2.14) and keep only two leading terms in the -expansion.
Note that the Hamiltonian of the massive variant of the model (1.1) also has the form (3.1) with a frequency ω n = πñ Λ 2 + m 2 ≈ πñ Λ +Λ 2πn m 2 , as m → 0. At the same time, the Moore's 3 In comparison, the direct calculation of β nk for "broken" trajectories (2.15) yields the following expression: which confirms the validity of approximations (2.14). The apparently troubling phase factor e iπ k is significant only at very large mode numbers, k ∼ 1/ , where its prefactor is itself proportional to O 2 .
approach [1], which we used to calculate the Bogoliubov coefficients in Sec. 2, does not work in the massive model due to the lack of conformal invariance. Hence, in the massive theory, approximation (2.14) and the second identity in (3.1) are valid only at relatively small evolution times,t 1/m 2Λ . Now we employ the same method to calculate the interacting Hamiltonian at the future infinity: where H int and H (1) int are the normal-ordered quartic interaction terms proportional to 0 and 1 , respectively: a l a m a n + g k+l+m+n a k a l a m a n + H.c. , (3.4) and δH free is the normal-ordered quadratic term, which can be absorbed into the renormalization of frequency in the free Hamiltonian (3.1): In the last identity, we single out the correction to the physical mass (which is, essentially, the one-loop correction): and introduce the ultraviolet cutoff n UV . In what follows, we assume that this correction is canceled by the corresponding tree-level counterterm 5 . Furthermore, we believe that this assumption does not affect our calculations. First, they are devoted to infrared rather than ultraviolet contributions. Second, we are interested in the evolution of the quantum state, which is not directly related to the mass of the field. For these reasons, in the following sections, we assume zero physical mass. Note that the Kronecker deltas in the H (0) int establish energy conservation in the scattering (a † a † aa) and decay (a † aaa and a † a † a † a) processes. This is a consequence of the emergent time translation symmetry in the limit → 0. Due to the same reason, there are no substantial loop corrections in this limit. Conversely, higher-order parts of the interacting Hamiltonian violate the energy conservation law 6 and allow usually forbidden processes. For instance, the aaaa term in the H (1) int describes a simultaneous production of four correlated particles. However, in what follows, we restrict ourselves to the leading nontrivial order in , i.e., assume that the "forbidden" processes occur only once or twice during the evolution of the system. We also remind that in this section, we truncate the -expansion of the interacting Hamiltonian (3.2) after the first order term. The full expression, which is expressed through the exact Bogoliubov coefficients and do not require the limit → 0, can be found in Appendix B.

Reduction to quantum mechanics
As explained in the previous section, the limit → 0 establishes an approximate energy conservation, i.e., restricts the energy exchange with the external world. Hence, we expect this limit to imply a small number of created particles. Furthermore, we can expect that created particles mainly populate the most easily excited, lowest-energy mode with n = 1. In this case, we can ignore all modes with n > 1 and reduce the interacting Hamiltonian (3.2) to a quantum-mechanical one: where we denote a ≡ a 1 for shortness. Recall that Fourier coefficients g 2 and g 4 , which contain the information about the motion of mirrors (see Eq. (2.11)), are small, g 2 ∼ g 4 ∼ , in the limit → 0. Keeping in mind the normal-ordered structure of this Hamiltonian and assuming the vacuum initial state, |in = |0 , we straightforwardly obtain the evolved quantum state in the limit λ → 0, t → +∞ and → 0: (4. 2) The leading contribution to Eq. (4.2) is ensured by multiple scattering of virtual particles (i.e., by powers of the a † a † aa term). Note that the T -ordered exponential is resolved into an ordinary exponential because in the limit in question, time-varying parts of the effective Hamiltonian are suppressed by powers of λ, so the identity [H(t 1 ), H(t 2 )] ≈ 0 approximately holds for large evolution times. In other words, T -exponential is resolved because we want to keep only leading secularly growing terms in its expansion. Identity (4.2) readily implies the approximate level density and anomalous quantum average: 3) The first and second terms on the right-hand side of (4.5) describe the tree-level and loop-level contributions, respectively. Note that the loop contribution is always positive and has the same order in as the tree-level approximation. Moreover, at very large times,t 1/λΛ, the oscillating part can be replaced with the average value: (4.6) Therefore, nonlinearities enhance the number of particles created in the fundamental mode and measured at large evolution times. Furthermore, the loop contribution significantly surpasses the tree-level approximation in the case |g 4 | |g 2 |, e.g., in resonant oscillations (2.16) with q = 4. We emphasize that the loop corrections to the number of created particles are related to a change in the initial quantum state, which becomes possible due to the violation of the energy conservation law in a background field. The leading correction, equation (4.2), describes a state with four correlated particles. Note that the full series in (4.2) is obtained by a unitary evolution from a pure state, although unitarity is violated in any fixed order in .

Large N generalization
Now let us show that the fundamental mode is indeed the most populous one. To that end, consider the O(N )-symmetric, large N generalization of the model (1.8): where i = 1, · · ·, N with N 1, and we assume the summation over the repeated upper indices. The effective Hamiltonian of this model in the limit λ → 0, t → ∞ and → 0 coincides with the Hamiltonian (3.2) after the appropriate change in the operator products: 3a † k a † l a m a n → (a i k ) † (a i l ) † a j m a j n + (a i k ) † (a j l ) † a i m a j n + (a i k ) † (a j l ) † a j m a i n , a † k a l a m a n → (a i k ) † a i l a j m a j n , a k a l a m a n → a i k a i l a j m a j n , 3a † k a l → (N + 2)(a i k ) † a i l , 3a k a l → (N + 2)a i k a i l , and coupling constant, λ → λ/N . Similarly to the original Hamiltonian (3.2), the main contribution to the generalized O(N )-symmetric Hamiltonian is ensured by elastic scattering processes 7 , and the leading relevant correction to this Hamiltonian is given by the production of four correlated quadruplets. Due to this reason, we believe that the qualitative behavior of both models are approximately the same, at least at small deviations from stationarity. At the same time, the evolved quantum state in the model (5.1) is straightforwardly calculated: where A p,q (t) = k,l,m,n δ p+q,k+l+m+n klm δ q,n e −iτ C k,l,m,n + iτ C k,l,m,n − 1 C 2 k,l,m,n  .7), we obtain the future asymptotics of the resummed quantum averages: (5.8) The first term in Eq. (5.8) describes the tree-level contribution, N free n ; the second and third terms appear only in the nonlinear, interacting theory.
Let us also explicitly evaluate N n (t) for a particular physically meaningful mirror motion. As an example of such a motion, we consider resonant oscillations with q = 2 frequency and destructive interference (2.16). We approximately calculate the particle number for relatively small, λΛt 1, and large, λΛt n UV / log(n UV ), evolution times, where n UV defines the effective ultraviolet cutoff. At intermediate times, these asymptotics are connected by a smooth curve, which can be calculated numerically (Fig. 1).
At relatively small evolution times, λΛt 1, the oscillating functions in (5.8) can be expanded in a series. The next-to-the-leading term in this series, which originates from the two-loop correction to the quantum averages, quadratically grows with time (compare with [45,46]): k,l,m,p |g n+k | 2 g 2 0 δ n,l+m+p + δ k,l+m+p (n + k)lmp ∼ 2 n N − 2τ 2 log 2 n n 2 + · · · , for 1 n 1 . In the last line, we discard the overall factor and fit the loop contribution, N loop n = N n − N free n , from the numerically calculated dependence (Fig. 2a). We emphasize that the total number of created particles is always positive, although the loop contribution is negative for n ≥ 4 and small evolution times (see Fig. 1).
At large evolution times, λΛt n UV / log(n UV ), we can replace the oscillating functions with their average values: n,k C 2 l,m,p,q − 2δ k,q + 2δ n,q C 2 l,m,p,q ∼ 2 n N + C log 2 n n + · · · , for 1 n 1 , where we again discard the overall factor and fit the loop contribution from the numerically calculated dependence (Fig. 2b). The factor C in front of the loop contribution is positive and slightly changes with the ultraviolet cutoff, C ∼ log log(n UV ). For reasonable values of n UV , this factor has the order of C ∼ 100. It is worth mentioning that Eqs. (5.9) and (5.10) are valid only for not very high mode numbers, n 1/ , where the Bogoliubov coefficients can be approximated by Eqs. (2.14). However, we emphasize that for every n 1/ loop contribution to the created particle number have the same order in as the tree-level expression.
We also note that the number of created particles is proportional to a small factor, 2 , and goes to zero at high frequencies, N n → 0 as n → ∞. This qualitatively validates the low-particle approximation of Sec. 4.

Nonlinear DCE in a Josephson metamaterial
Now let us apply the approach of Secs. 2-4 to another model of the nonlinear DCE, which has a clear experimental implementation [17]. Namely, consider the Hamiltonian of a Josepshon metamaterial, i.e., an array of superconducting quantum interference devices (dc-SQUIDs) in a varying external magnetic field: The effective cavity length Λ, speed of light v(t, x), and coupling constant λ are specified by the parameters of the metamaterial (see Appendix C for the details). Moreover, variations of the external magnetic field directly translate into the variations of the light speed v(t, x). We assume that these variations are small, nonresonant, and vanish at the past and future infinity, i.e., v(t, x) = v ∞ + ṽ(t, x), where v ∞ = const, 1 andṽ(t, x) → 0 as t → ±∞. In other words, we consider the Hamiltonian (6.1) for small deviations from stationarity.
First of all, we split the Hamiltonian (6.1) into the free (quadratic) and interacting (quartic) parts and quantize the free model similarly to Sec. 2: Here, a † n and a n are the standard bosonic creation and annihilation operators; in-modes f in n (t, x) solve the corresponding free equations of motion, form the complete orthonormal basis with respect to the Klein-Gordon inner product and have definite positive frequency at the past infinity 8 : At the future infinity, in-modes are rewritten in terms of the out-modes: , as t → +∞, (6.4) where Bogoliubov coefficients α nk and β nk are determined by the variations of the speed of light. If these variations are small and nonresonant, Bogoliubov coefficients are close to the identity, α nk ≈ δ n,k + α nk n/k and β nk ≈ β nk n/k (compare with (2.14)).
Then we substitute the in-modes into the interacting Hamiltonian, expand it to the first order in , take the integral over dx and neglect rapidly oscillating terms: klmn 4δ k,l+m+n a † k a l a m a n + (3δ k+l,m+n + 6δ k,m δ l,n ) a † k a † l a m a n + H.c. , k,l,m,n √ klmn 4β k,l+m+n a k a l a m a n + f (1,3) k,l,m,n a † k a l a m a n + f (2,2) k,l,m,n a † k a † l a m a n + H.c. , as t → +∞. Here, f (1,3) k,l,m,n and f (2,2) k,l,m,n are some constant O(1) tensors that depend on the Bogoliubov coefficients; we do not need their explicit form to determine the leading correction to the number of created quasiparticles in the limit 1. In fact, in this limit, the leading loop contributions to the level density (1.6) and anomalous quantum average (1.7) are approximated by the following expressions (we set t 0 = 0 for shortness): Keeping in mind that at small deviations from the stationarity, created quasiparticles mainly populate the lowest-energy mode (compare with Secs. 4 and 5), we truncate the summations at n UV = 1 for a rough estimate of the effective Hamiltonian: We emphasize that at large evolution times, t t * ∼ Λ 3 /λ , the loop contribution is comparable to the tree-level value, N free 1 ≈ 2 β 1,1 2 . For realistic parameters of the Josephson metamaterial 9 , this time is approximately t * ∼ 10 −5 s. In other words, small intrinsic nonlinearities of the dc-SCQUIDs are negligible only at relatively small times, t t * ; at larger times, the measured number of created particles will deviate from the tree-level value 10 .

Discussion and Conclusion
In this work, we have shown that nonlinearities (i.e., interactions between the quantum fields) significantly affect particle creation in a weakly nonstationary scalar DCE. We emphasize that at large evolution times, the correction to the particle number is comparable to the tree-level value even if the coupling constant is small, i.e., nonlinearities are apparently weak. First of all, we calculated the Bogoliubov coefficients and derived the effective Hamiltonian of the quantized field. Then we employed this Hamiltonian to evaluate the nonperturbative (in the coupling constant λ) correction to the tree-level number of created particles. We also discussed the effective Hamiltonian of a Josephson metamaterial, where DCE can be measured experimentally. In both cases, created particles mainly populate the fundamental (lowest-energy) mode. More importantly, in both cases tree-level and loop-level contributions to the number of created particles have the same order in the small parameter , which determines the magnitude of deviations from stationarity.
We emphasize that in this paper, we work in the interaction picture. In other words, we separate the free, Eq. (3.1), and interacting, Eq. (3.2), Hamiltonians and assume that we exactly know how the free Hamiltonian evolves the field operator φ(t). This knowledge is contained in the Bogoliubov coefficients. Then we approximately calculate the final quantum state, i.e., evolve the initial quantum state with the interacting Hamiltonian. This information is contained in the resummed level density and anomalous quantum average. Note that the quantum state continues evolving even when the Bogoliubov coefficients stop changing, i.e., free evolution is essentially finished. This approach is inspired by similar problems from the high energy physics [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44].
Note that we restricted ourselves to the calculation of the number of created particles, Eq. (1.4) with the resummed level density (1.6) and anomalous quantum average (1.7). Nevertheless, our analysis can be straightforwardly applied to the calculation of the resummed Keldysh propagator, which is also expressed through the quantum averages (1.6) and (1.7). Using this propagator, one can easily obtain other observable quantities (e.g., see Appendix A).
Besides, we emphasize that the summation of loop corrections in Secs. 4-6 corresponds to the solution of a reduced system of Schwinger-Dyson equations on exact correlation functions written in the Schwinger-Keldysh technique. Namely, for this summation, we single out only "tadpole" and "bubble" contributions to propagators and vertices, respectively. In Secs. 4 and 6, such a reduction is achieved by employing the limit → 0 and λ → 0, t → ∞; in the model of Sec. 5, it naturally emerges in the large N limit. Then, as soon as the exact Keldysh propagator is calculated, quantum averages and created particle number are unambiguously restored. More 9 In the experiment [17], authors created a Josephson metamaterial with the following parameters: length Λ = 4 mm, number of nodes N = 250, Josephson critical current I c ∼ 10 −5 A, and the effective speed of light in the unperturbed circuit v ∞ ∼ 0.5c, with c being the speed of light in the vacuum. 10 We remind that we resum only the leading secularly growing corrections in the limit λ → 0, t → ∞, λt = const, i.e., keep only the terms of the form (λt) k with an integer k > 0. The subleading corrections has the form λ m t n with some integers m > n > 0; at finite λ and t, these contributions are suppressed by the powers of λ. For realistic parameters of the model, see footnote 9 and paper [17], these subleading terms become significant only at the time scale t * * ∼ Λ 5 v ∞ /λ 2 2 ∼ 1 s, which is much larger that t * . details on the relation between the different approaches to the summation of loop corrections to the created particle number can be found in [43], where a simplified model of the DCE with a similar diagrammatics is discussed.
Our work can be extended in several possible directions. First, we expect that the effects studied in this paper can be measured in a quantum circuit analog of the DCE [15][16][17] (see also [68][69][70][71][72][73][74]). Namely, we expect that at large evolution times, t t * ∼ 10 −5 s, the number of the quasiparticles created in a quantum circuit will deviate from the tree-level value due to the intrinsic nonlinearity of Josephson junctions (see Sec. 6). In other words, this means that the Gaussian approximation used in [15][16][17] for a theoretical estimate of N is valid only at relatively small times.
However, we remind that in this paper, we considered only small deviations from the stationarity. At the same time, feasible experimental implementations of models (1.8) or (1.9) require a resonant pumping of energy into the system, i.e., strongly deviate from the stationarity. In this case, Bogoliubov coefficients and number of created particles rapidly grow with time already at the tree level [5,[13][14][15][50][51][52]. Therefore, it is promising to extend the results of the present paper to strong deviations from stationarity, including the resonant pumping of energy into the system. In the large N limit, such an extension can be possibly implemented using the Schwinger-Keldysh technique similarly to papers [28,36,37,43]. We expect that in strongly nonstationary systems, nonlinearities also significantly affect particle production through the generation of quantum averages n kl and κ kl because these averages are eventually multiplied by large Bogoliubov coefficients and amplified on a par with vacuum fluctuations, see Eq. (1.4) (also compare with the simplified model of [43] where calculations for strong deviations from stationarity were performed).
Second, the expression for the effective Hamiltonian is easily generalized to other realistic models of the DCE, including massive quantum fields, semitransparent mirrors, and resonant oscillations with constructive interference. In general, the strategy is similar to Secs. 2 and 3: one calculates the Bogoliubov coefficients for the free field and substitutes them into the Hamiltonian (B.1).
Finally, two-dimensional O(N ) and CP N models on a finite interval are known to possess nontrivial classical solutions with a negative energy [79][80][81][82][83][84][85]. In these models, the naive vacuum, a n |in = 0, is unstable at some values of the model parameters. In principle, nonequilibrium techniques (similar to those used in Sec. 5) might confirm this instability and explicitly trace the evolution from the initial quantum state to the true vacuum. An example of such a calculation in a different setup can be found in [40]. Furthermore, it is promising to generalize the predictions of [79][80][81][82][83][84][85][86], which were obtained in a stationary approximation, to time-dependent boundary conditions. The calculations of Sec. 5 may be considered as a starting point for such generalizations.
A Alternative derivation of the created particle number Correlation functions are the most fundamental objects in the nonstationary quantum field theory [44]. First, they have a clear physical meaning even in time-dependent backgrounds where the notion of a particle is ill-defined. Second, all directly observable quantities are derived from correlation functions. For example, the exact Hamiltonian of the nonlinear DCE (1.8) can be obtained from the exact Keldysh propagator: where we denote x = (t, x) for shortness. At the tree level, the Keldysh propagator is defined as an anticommutator of quantum fields: where we employ the mode decomposition (2.1) and use short notations for the level density and anomalous quantum average. In the interacting theory, the Keldysh propagator preserves the form (A.2) with exact quantum averages (1.6) and (1.7) if the difference between the times is much smaller than their average, t 1 − t 2 (t 1 + t 2 )/2. This allows us to estimate the future asymptotics of the exact Hamiltonian in the interacting theory: δ k,l + n kl (α km α * lm + β km β * lm ) + κ kl (α km β lm + β km α lm ) + H.c. + O λΛ , (A.3) where we take the limit t → +∞, substitute the asymptotics of the in-modes (2.8), rotate to the (t,x) coordinates and integrate overx. Note that the contribution of the nonlinear part is negligible if we consider small coupling constants, λ 1/Λ 2 . Finally, using the properties of the Bogoliubov coefficients: and symmetry of the level density, n * kl = n lk , and anomalous quantum average, κ kl = κ lk , we obtain the following expression for the exact Hamiltonian: where N m is given by equation (1.4) with level density (1.6) and anomalous quantum average (1.7). The divergent sum ∞ m=1 πm 2Λ describes the standard static Casimir energy and does not depend on the mirrors motion at intermediate times. The second sum describes excitations above stationarity. This confirms that N m has the meaning of the particle number created in a linear or nonlinear DCE at large evolution times.
Here, we introduce the short notation for the symmetrized quantities, e.g., A (m,n) = 1 2! (A mn + A nm ). We also neglect rapidly oscillating and constant terms that do not contribute to (1.6) or (1.7) in the limit in question. In other words, we average the Hamiltonian, H int (t) −→ 1 T t+T t H int (t )dt , and neglect the subleading contributions in the limit T → +∞.
It is worth mentioning that the Kronecker deltas of the form δ i+j,k+l appear only in a massless theory. In a massive two-dimensional theory, these deltas are multiplied by the time-dependent function, sin[(ω i +ω j −ω k −ω i+j−k )T ] (ω i +ω j −ω k −ω i+j−k )T , where ω k = (πk/Λ) 2 + m 2 . This function goes to one only if i = k, i = l, or m → 0; otherwise, it vanishes at large evolution times, T 1/m 2Λ . This behavior illustrates the well-known fact that elastic scattering of massive particles in a flat two-dimensional spacetime reduces to a simple momentum exchange.
Finally, substituting the approximate Bogoliubov coefficients (2.14) into this identity and discarding the subleading in terms, we straightforwardly obtain Eq. (3.2).

C The Hamiltonian of a Josephson metamaterial
Following [74] and [17], we straightforwardly obtain the Hamiltonian of a Josephson metamaterial that consists of an array of dc-SQUIDs (Fig 3): Here, Φ n (t) = t −∞ V n (t )dt is the flux variable associated with the n-th node, Q n is the charge at that node, Φ 0 = π /e is the magnetic flux quantum, C 0 is the capacitance, and E J is the Figure 3: A schematic representation of the Josephson metamaterial [17], which models the DCE in a cavity with a varying speed of light. The square with a cross represents a Josephson junction.
Josephson energy of each SQUID. In an external magnetic field, the Josephson energy varies, E J (Φ ext ) = 2E J (0) cos πΦext Φ 0 (we assume that both Josephson junctions of the SQUID have approximately the same properties). The total length of the array is Λ, and the distance between the adjacent nodes is fixed, δx = Λ/N . Note that the array is grounded at both ends, which models the Dirichlet boundary conditions for propagating excitations. Now let us derive the continuum limit of this Hamiltonian. To that end, we introduce the capacitance and charge per unit length, c 0 = C 0 /δx and Q(t, x n ) = Q n (t)/δx, use the identity Q n (t) = C 0 ∂ t Φ n (t) and expand the cosine to the fourth order:

(C.2)
For shortness, we also introduce the inductance per unit length, l 0 (t, x) = In the last identity of (C.2), we also redefine the flux variable, φ(t, x) = √ c 0 Φ(t, x). Note that in practice, variations of the external magnetic field are usually spatially uniform, Φ ext (t, x) = Φ ext (t). However, we consider arbitrary functions Φ ext (t, x) for generality. Now it is easy to see that in the naive limit λ → 0, which is justified at small evolution times, t Λ 3 / λ, Josephson metamaterial effectively reproduces the standard model of the scalar DCE, Eq. (1.1). Otherwise, Hamiltonian (C.2) implies the equations of motion (1.9).