Superconducting circuit boundary conditions beyond the Dynamical Casimir Effect

We study analytically the time-dependent boundary conditions of superconducting microwave circuit experiments in the high plasma frequency limit, in which the conditions are Robin-type and relate the value of the field to the spatial derivative of the field. We give an explicit solution to the field evolution for boundary condition modulations that are small in magnitude but may have arbitrary time dependence, in a formalism that applies both to a semiopen waveguide and to a closed waveguide with two independently adjustable boundaries. The correspondence between the microwave Robin boundary conditions and the mechanically-moving Dirichlet boundary conditions of the Dynamical Casimir Effect is shown to break down at high field frequencies, approximately one order of magnitude above the frequencies probed in the 2011 experiment of Wilson et al. Our results bound the parameter regime in which a microwave circuit can be used to model relativistic effects in a mechanically-moving cavity, and they show that beyond this parameter regime moving mirrors produce more particles and generate more entanglement than their non-moving microwave waveguide simulations.


Introduction
The quantum theory of relativistic fields with moving boundaries was first explored by Moore in a remarkably original paper on the quantum formulation of linearly polarized light in a one-dimensional moving cavity [1]. The primary result of this investigation was the discovery that moving mirrors in vacuum create photons. Later, motivated by developments in quantum field theory in curved spacetimes, the specialization to a single moving mirror in Minkowski spacetime was carried out by Fulling and Davies [2], who again found that non-uniformly accelerating mirrors generate radiation. These effects, in which particles are produced by moving boundaries, are generally referred to as the Dynamical Casimir Effect (DCE) or the Nonstationary Casimir Effect [3,4,5].
An experimental verification of the DCE in a system with a mechanically moving boundary has remained elusive because of the technological challenges in creating sufficiently large accelerations [3,4,5]. However, an experimental observation of a similar particle creation effect has been reported in a mechanically static semiconductor waveguide where the boundary condition on the field is modulated electrically, by a superconducting quantum-interference device (SQUID) [6]. A related experimental observation of particle creation in Josephson metamaterial has been reported in [7]. These observations open fascinating prospects for simulating on a mechanically static desktop device quantum phenomena due to motion, including entanglement generation and degradation, in a regime where the moving system would experience significant relativistic effects [8,9,10].
In this paper we address the evolution of a quantum field in waveguides of the type used in the experiment of [6] in situations where the modulation of the SQUID(s) at the end(s) of the waveguide is small in magnitude but may have arbitrary time-dependence, under the further assumption that the plasma frequency of the SQUID(s) is negligibly high compared with the frequencies where the experiment operates. More precisely, recall that the field Φ(t, x) in the experiment of [6] satisfies, under the approximations described in [11,12], the (1 + 1)-dimensional Klein-Gordon equation with the boundary condition where the cavity is at x ≤ 0, the SQUID is at x = 0, the meaning of the positive constants v, Φ 0 , C J and L 0 is as described in [12], and E J (t) can be given arbitrary time-dependence by modulating the magnetic field applied to the SQUID. In the regime where the SQUID's plasma frequency is large compared to the frequency of Φ, the time derivative term in (1.2) is negligible, and (1.2) reduces to where L eff 0 (t) = Φ 0 /(2π) 2 L 0 E J (t) −1 . We shall consider the regime in which the boundary condition (1.3) applies and the time-dependence of L eff 0 (t) is arbitrary in profile but small in magnitude. The regime where the SQUID's plasma frequency is large but not negligibly so is considered in [13,14,15].
Our interest in the boundary condition (1.3) is twofold. First, considering the condition in its own right, we give an explicit solution to the quantum dynamics to leading perturbative order in the time variation of the boundary condition, using a formalism that allows us to handle both the semiopen waveguide and a closed waveguide, and allowing for the closed waveguide the modulations at the two ends to remain independent of each other. The mathematical observation underlying our formalism is that (1.3) is a linear relation between the value of the field and the spatial derivative of the field, known as a Robin boundary condition, and when this condition is independent of time, it has a well-known role in making the spatial part of the wave operator in (1.1) self-adjoint and hence ensuring the unitarity of the time evolution [16,17,18]. A time-dependent boundary condition can then be treated perturbatively by combining the spectral methods of [16,17,18] to the techniques developed in [19,20] for mechanically moving cavities. For the semiopen waveguide our results agree with those found in [21,22,23,24] via a different formalism.
Second, we wish to address the sense in which (1.3) models a mechanically moving boundary of the DCE. For a Dirichlet mirror at the time-dependent location x = x DCE (t), the boundary condition on the field reads 3) reduces to (1.4) for field frequencies much smaller than v/|L eff 0 |, but for higher field frequencies the correspondence no longer holds. We shall see that our perturbative solution of the field evolution with the condition (1.3) indeed differs from the similar perturbative solution with the condition (1.4) for frequencies that are not much smaller than v/|L eff 0 |, both for a semiopen waveguide and a closed waveguide; in particular, the large frequency falloff properties of the solution are qualitatively different. Simulations of relativistic motion with the mechanically static semiconductor waveguide would hence need to take place in the low frequency regime where the successive approximations from (1.2) via (1.3) to (1.4) hold. Both the experiment of [6] and the proposals of [8,9,10] appear to operate within in this domain by a margin of approximately one order of magnitude.
The plan of the paper is as follows: Section 2 considers evolution under a small discontinuous change in the boundary condition (1.3) for a semiopen waveguide, and Section 3 presents the similar analysis for the closed waveguide. Evolution under small changes in the boundary condition with arbitrary time-dependence for both types of waveguides is written out in Section 4. Section 5 compares the evolution to that under Dirichlet boundary conditions at one or two mechanically moving boundaries. Section 6 presents a summary and concluding remarks.
Appendix A collects technical identities in a perturbative expansion of Bogoliubov coefficients. Appendix B treats the evolution under the Dirichlet boundary condition at one mechanically moving boundary in a small acceleration expansion in which velocities and travel distances are unrestricted, adapting to one moving boundary the treatment of mechanically rigid cavities given in [19,20]. In appendix C we derive the first order formula for the negativity measure of entanglement in the case when the modes have a continuous spectrum.

Semiopen waveguide
In this section we discuss a semiopen waveguide under a sudden change in the Robin boundary condition (1.3). Subsection 2.1 establishes the notation and reviews the known properties of a time-independent Robin boundary condition. The sudden change is implemented in subsection 2.2, and small sudden changes about the Robin boundary condition relevant for the experiment of [6] and about the Dirichlet boundary condition are discussed respectively in subsections 2.3 and 2.4.

Static boundary condition
We adopt units in which the phase velocity v in (1.1) is set to unity. We may hence think of the field as a real scalar field φ on a (1 + 1)-dimensional Minkowski spacetime, with the global Minkowski coordinates (t, x) and the metric ds 2 = −dt 2 + dx 2 . We take the boundary to be at x = 0 and the field to live in the half 0 ≤ x < ∞.
We consider the massive Klein-Gordon field equation where µ ≥ 0 is the mass. The massless special case µ = 0 reduces to (1.1). We keep here µ general, in part because setting µ = 0 does not significantly simplify the analysis, but also in part in view of prospective future comparisons with mechanically moving cavities in situations where transverse dimensions may generate a positive µ by Fourier decomposition [19,20,25]. We introduce at x = 0 the Robin family of boundary conditions where D ∈ R ∪ {∞} is a constant independent of t. The special case D = 0 gives the Dirichlet boundary condition, φ(t, 0) = 0, the special case D = ∞ gives the Neumann boundary condition, φ (t, 0) = 0, and all other values of D mix φ(t, 0) and φ (t, 0). Any choice for D makes −∂ 2 x self-adjoint [16,17,18] and would hence yield a consistent quantum theory of a non-relativistic particle on the half-line, although the choice D = 0 may in that context be considered less fine-tuned than the others [26]. In the present context of a relativistic quantum field theory, we consider only the values of D that make the spectrum of −∂ 2 x + µ 2 strictly positive. For µ = 0 this means D ∈ (−∞, 0] ∪ {∞}, and for µ > 0 it means D ∈ (−∞, 0] ∪ {∞} ∪ (µ −1 , ∞). A negative eigenvalue of −∂ 2 x + µ 2 would give a tachyonic instability, and the zero eigenvalue that occurs when µ > 0 and D = µ −1 would give a zero mode and hence a theory without a Fock vacuum.
The field is quantised in the usual fashion. A continuum of mode solutions that are positive frequency with respect to ∂ t are where k > 0, ω = k 2 + µ 2 , δ is determined as a function of k from tan δ = −kD , (2.4) and we may fix the phase by choosing |δ| < π/2 for D ∈ R and δ = π/2 for D = ∞. When µ > 0 and µ −1 < D < ∞, there is in addition a discrete ground state mode, given by Writing the Klein-Gordon inner product in the conventions of [27] as where the overline denotes complex conjugation, the continuum modes φ k are Diracorthonormal, and when the discrete ground state φ g (2.5) exists, it is normalised and orthogonal to the continuum modes. The Fock space is built on the vacuum that is annihilated by the annihilation operators associated with the continuum modes φ k and with the discrete mode φ g when the latter exists.

Sudden change in the boundary condition
Suppose that for t < 0 we use the field modes as introduced above and for t > 0 we use a similar set of field modes with D replaced by D . Denoting the new continuum modes byφ k and the new discrete mode byφ g , we may match the two sets of modes at t = 0 in the notation of [27] as where we have included the lower left subscript o in the Bogoliubov coefficients o α and o β to indicate that the change in the boundary condition is sudden. Expressions for the Bogoliubov coefficients can be found by taking inner products of (2.8) with the untilded modes and their complex conjugates [27]. For the continuumto-continuum Bogoliubov coefficients, we find where δ and δ are determined by and P in (2.9a) denotes the principal value in the integration over k at k = k . In the special cases D = 0 and D = 0 the formulas (2.9) are understood in their welldefined limiting sense. The continuum-to-discrete and discrete-to-continuum Bogoliubov coefficients will not be needed below and we omit the formulas. We note that while o β k k (2.9b) is a function for all k and k , o α k k (2.9a) has distributional support at k = k because of the Dirac delta in the first term and the principal value integral in the second term.

Small sudden change: far from Dirichlet
We now consider the case where D and D are negative and close to each other. This is the situation relevant for the the waveguide experiment of [6]. (Recall that in our conventions the field φ lives at x ≥ 0, while the field Φ in (1.1)-(1.3) lives at x ≤ 0.) We write where Λ is a positive constant of dimension length and η is small. The mode expansion contains no discrete mode for any value of µ. Expanding the Bogoliubov coefficients (2.9) in η, we find As a consistency check, we have verified that the expansion (2.12) is consistent with the identities satisfied by the Bogoliubov coefficients, collected in Appendix A, in the sense that the linear order identities (A.4a) hold and the off-diagonal part of the quadratic order identities (A.4b) holds. We are not aware of reasons to suspect inconsistencies in the diagonal part of (A.4b) but we have not undertaken the distributional analysis to examine this.

Small sudden change: near Dirichlet
As a second regime of interest, we perturb around the Dirichlet boundary condition. This case is not directly relevant for the experiment of [6], but we shall see in Section 5 that it exhibits close theoretical similarity with a mechanically moving boundary.
We set D = 0 and assume D to be close to 0. If D is negative, the result may be obtained from (2.12) by writing ηΛ = −b and letting Λ → 0 while b remains finite and negative but small: then D = b < 0, and the Bogoliubov coefficients are given by If D is positive, the t > 0 theory has a tachyonic instability (for all positive D if µ = 0 and for 0 < D < µ −1 if µ > 0) because of the negative eigenvalue of −∂ 2 x + µ 2 ; however, the tachyon is nonperturbative in D , and we have verified that just ignoring the tachyon and proceeding directly from (2.9) leads again to (2.13), where now D = b > 0.
As a consistency check, we have verified that (2.13) satisfies the linear order Bogoliubov identities (A.4a) and the off-diagonal part of the quadratic order Bogoliubov identities (A.4b), regardless the sign of b.

Closed waveguide
In this section we adapt the analysis of Section 2 to a cavity waveguide that is closed at both ends. For simplicity, we treat only the massless field, µ = 0.

Static boundary condition
We follow the notation of Section 2, setting µ = 0. We place the cavity at 0 ≤ x ≤ L, where the positive constant L is the length of the cavity.
We write the static Robin boundary conditions at the ends of the cavity as where D 1 and D 2 are constants independent of t, taking values in R ∪ {∞}. Any choice for D 1 and D 2 makes −∂ 2 x self-adjoint [16,17,18]. To avoid instabilities and zero modes, we assume initially D 1 and D 2 to be such that the spectrum of −∂ 2 x is strictly positive. We shall however see in subsection 3.3 that a perturbative treatment in D 1 and D 2 remains consistent even in the presence of a nonperturbative tachyon, on a par with what happened for the semiopen waveguide in subsection 2.4.
The field is again quantised in the usual fashion. The Klein-Gordon inner product is as in (2.6) but integrated over x ∈ [0, L]. The equation for the eigenvalues can be written down using (3.1).

Small sudden change: far from Dirichlet
We consider the case where D 1 < 0 and D 2 > 0, and there is a sudden but small change in their values. This models a waveguide whose each end terminates at a SQUID as in the experiment of [6].
For t < 0, we set D 1 = −κ 1 L and D 2 = κ 2 L, where κ 1 and κ 2 are positive dimensionless constants. With this boundary condition −∂ 2 x is positive definite. The mode functions are q goes over the positive solutions to tan δ q = κ 1 q, and we choose the phase so that 0 < δ q < π/2. There is exactly one q in each interval mπ < q < (m + 1)π, m = 0, 1, 2, . . .. The mode functions are normalised to (φ q , φ q ) = δ qq . For t > 0, we set D 1 = (−κ 1 + η 1 )L and D 2 = (κ 2 + η 2 )L, where η 1 and η 2 are dimensionless constants, assumed to be small. We work perturbatively in η 1 and η 2 , setting both of them to be proportional to a formal expansion parameter η which at the end is set to unity. The mode functions are proportional to e −ikt sin(kx + δ) where k and δ are determined from (3.1). We label the mode functions by the positive solutions to such that k = p/L+O(η), we denote them byφ p , and we choose their phase to agree with that of (3.2) in the zeroth perturbative order. We then find that the eigenfrequencies are given by The Bogoliubov coefficients are now defined by matching the modes at t = 0 as Proceeding as in [19,25], we find where the map q → ϕ q labels the consecutive solutions to (3.4) by consecutive integers. The expressions for the order η 2 terms in (3.6) and (3.8) are lengthy and we suppress them here. As a consistency check, the linear terms in (3.8) satisfy the linear order Bogoliubov identities (A.4a). We are not aware of reasons to suspect inconsistencies in the quadratic order identities (A.4b) but examining these would require a nontrivial evaluation of the left-hand side in (A.4b) and we have not carried out this evaluation.

Small sudden change: near Dirichlet
We consider also a perturbation around the Dirichlet boundary condition. If η 1 < 0 and η 2 > 0, the result may be obtained from (3.8) simply by taking the limit κ 1 → 0 and κ 2 → 0. The t < 0 mode functions are given by where n = 1, 2, . . .. Using positive integers to label both the t < 0 mode functions and the t > 0 mode functions, we find that the t > 0 eigenfrequencies are given by where m = 1, 2, . . ., and the Bogoliubov coefficients are given by where we have now displayed also the order η 2 terms. If η 1 ≥ 0 and/or η 2 ≤ 0, the t > 0 theory may have tachyonic instabilities or zero modes; however, both of these are nonperturbative, and we have verified that setting D 1 = D 2 = 0 for t < 0, D 1 = η 1 L and D 2 = η 2 L for t > 0, ignoring any tachyons or zero modes, and working directly from (3.1) and (3.7), yields (3.10) and (3.11) regardless the signs of η 1 and η 2 . As a consistency check, the expressions in (3.11) can be verified to satisfy the perturbative Bogoliubov identities (A.4). The elements of the matrix square on the left-hand side of (A.4a) are given by absolutely convergent sums that can be evaluated by residue techniques.

Boundary condition with arbitrary time-dependence
When the boundary condition has arbitrary time dependence but the variations remain so small in magnitude that first-order perturbation theory suffices, the evolution of the field can be obtained by composing the sudden changes of Sections 2 and 3 and passing to the limit [20]. We discuss first the semiopen waveguide and then the closed waveguide.

Semiopen waveguide
For the semiopen waveguide of Section 2, we consider a boundary condition of the form (2.2) where D may change in time but only within the interval t 0 ≤ t ≤ t f .

Far from Dirichlet
Consider first the far-from Dirichlet case. We write D = −Λ 1 + η(t) , where Λ is a positive constant and the function η(t) is vanishing outside the interval t 0 ≤ t ≤ t f and satisfies |η(t)| 1.
At t < t 0 , we introduce early time basis modes that are as in (2.3) but with the replacement e −iωt → e −iω(t−t 0 ) . At t > t f , we similarly introduce late time basis modes that are as in (2.3) but with the replacement e −iωt → e −iω(t−t f ) . Let α k k and β k k be the coefficients in the Bogoliubov transformation from the early time modes to the late time modes. Working perturbatively in η, we may proceed as in [20], and the outcome can be read off from formulas (6) and (7) therein. We find When µ = 0, (4.1) and (4.2) reduce to what was found in [21,22,23,24] via a different formalism.

Near Dirichlet
In the near-Dirichlet case, we take D = b(t), where the function b(t) is vanishing outside the interval t 0 ≤ t ≤ t f . Proceeding as above, we find

Closed waveguide
For the closed waveguide of Section 3, we consider a boundary condition of the form (3.1) where D 1 and D 2 may change in time but only within the interval t 0 ≤ t ≤ t f . We may proceed as above. The only new aspect is that the method of [20] needs to be generalised to accommodate the linear term that appears in the frequencies (3.6) and (3.10).

Far from Dirichlet
In the far-from Dirichlet case, we write D 1 = −κ 1 + η 1 (t) L and D 2 = κ 2 + η 2 (t) L, where κ 1 and κ 2 are positive constants and the functions η 1 (t) and η 2 (t) are vanishing outside the interval t 0 ≤ t ≤ t f and satisfy |η 1 (t)| 1 and |η 2 (t)| 1. Indexing the mode functions in the notation of subsection (3.2), writing and proceeding as above, we find that the coefficients in the Bogoliubov transformation from the early time modes to the late time modes are given by

Near Dirichlet
In the near-Dirichlet case, we write D 1 = η 1 (t)L and D 2 = η 2 (t)L. Indexing the past and future mode functions by positive integers in the notation of subsection (3.3), and writing ω m = πm/L , ω n = πn/L , (4.8) we find

Comparison
We are now ready to compare the evolution under the time-dependent Robin boundary condition to the evolution under the Dirichlet condition at a mechanically moving boundary.

Semiopen waveguide
For the semiopen waveguide, the evolution of a field with mass µ ≥ 0 under the timedependent Robin boundary condition was given in subsection 4.1. The evolution of a massless field under a Dirichlet condition at a mechanically moving boundary is given in Appendix B in terms of the acceleration of the boundary, in a small acceleration approximation that allows the velocity and the travel distance to remain arbitrary and overlaps with the DCE literature results [3,4,5] in the common domain of validity. Comparing (4.3)-(4.4) and (B.11)-(B.12), we see that the massless field with the mechanically moving boundary can be simulated to the leading order in perturbation theory by the µ = 0 near-Dirichlet Robin boundary condition provided we choose b(t) so that a(τ ) = ∂ 2 τ b(τ ), where a(τ ) is the proper acceleration of the boundary as a function of its proper time τ , and the modulation starts and ends so gently that both b andḃ vanish. This is precisely the relation one would have expected from the low frequency equivalence between the Robin boundary condition (1.3) and the mechanically moving Dirichlet boundary condition (1.4). The simulation is reliable in the frequency range where the first-order perturbation theory results on both the Robin side and on the mechanical side remain reliable; we shall not attempt to quantify this range more precisely, but for given a(τ ) and b(τ ) the range is unlikely to include arbitrarily high frequencies.
Comparing further (4.3)-(4.4) and (B.11)-(B.12) with (4.1)-(4.2), we see that the massless field with the mechanically moving boundary can be simulated by the µ = 0 far-from-Dirichlet Robin boundary condition provided we choose η(t) so that a(τ ) = −Λ∂ 2 τ η(τ ), the frequencies are much smaller than Λ −1 , and the modulation starts and ends so gently that both η andη vanish. Again, this is precisely the outcome one would have expected from (1.3) and (1.4). When the frequencies are not much smaller than Λ −1 , the evolution with the Robin boundary condition differs from the evolution with the mechanically moving boundary because the square root factors in (4.2) differ from unity.
Investigating these differences further, recall that the beta coefficients are directly related to the total photon production number by: Taking η(t) to be a sinusoidal function: with positive constant ε 1 and driving frequency ω d , the integrals in (4.2) can be performed exactly. The dominant part of the photon flux occurs for frequencies below the driving frequency. We therefore define the photon flux density, n(k), as a function of the reduced frequency k := k /ω d , by 1/ω d . In both cases the spectrum has in the range 0 < k < 1 a distinctive parabolic shape that is qualitatively characteristic of the DCE, and for Λ 1/ω d there is good agreement with the zero temperature curve of Figure 8 in [12].
For a quantitative comparison with the DCE, Figure 1 presents also the photon flux density for the mechanically moving Dirichlet boundary, obtained from (B.11)-(B.12) with the matching a(τ ) = −Λ∂ 2 τ η(τ ). When the effective length, Λ, is much smaller than the effective wavelength of the driving frequency (Λ 1/ω d ), the moving boundary and Robin-boundary systems produce almost identical spectral functions, even though the modulation used for the plot has nonvanishingη at the start and end. On the other hand, when the effective length is larger than the effective wavelength of the driving frequency (Λ 1/ω d ), the spectral functions from the moving boundary system and the Robin-boundary systems differ: the radiation from the moving boundary has a higher intensity, and the spectrum from the Robin boundary condition decays more rapidly at high frequencies. Both of these behaviours result from the square root factors which suppress the beta coefficients in (4.2).
One of the signatures of the DCE radiation is that emitted photons are entangled. We can quantify the amount of entanglement produced using the negativity [28,29] which is defined as minus the sum of the negative eigenvalues of the partially transposed state. We show in appendix C that for perturbative Bogoliubov coefficients of the type (4.1)-(4.2) the leading term in the negativity between sharply peaked field modes k and k is given by ∆k|B k k | where ∆k is the spectral line-width of the frequencies.
For the driving function given in equation (5.2), |B k k | is highest for values of k + k near the driving frequency. Therefore, we define ∆ω = (k − k )/2, and we expand the frequencies near half driving-frequency as: The results for |B k k | as a function of ∆ω/ω d are shown in Fig. 2. These plots are qualitatively similar to those of Fig. 14 in [12] for the "without resonator" line. We again observe good agreement when Λ 1/ω d between moving mirrors and non-moving microwave waveguides. Also, when Λ 1/ω d we find that the moving mirrors generate more entanglement than their associated non-moving microwave waveguide simulations. (Left) In the small effective length limit Λ 1/ω d the two systems display near identical behavior. (Right) In the large effective length limit Λ 1/ω d the moving system produces more entanglement than the non-moving system. We use the same parameters as in Fig. 1.

Closed waveguide
For the closed waveguide, the evolution of a massless field under the time-dependent Robin boundary condition was given in subsection 4.2. The evolution of a massless field in a moving, mechanically rigid cavity of proper length L with Dirichlet boundary conditions is given by [19,20] m and n are positive integers, ω m and ω n are given by (4.8), τ is the proper time at the centre of the cavity and a(τ ) is the proper acceleration at the centre of the cavity. Comparing (4.9)-(4.10) and (5.6)-(5.7), we see that the massless field in the rigidly moving cavity can be simulated to the leading order in perturbation theory by the near-Dirichlet Robin boundary condition provided we choose η 1 (τ ) = η 2 (τ ) and a(τ ) = L∂ 2 τ η 1 (τ ), and the modulation starts and ends so gently that both η andη vanish. Again, this is precisely the relation one would have expected from (1.3) and (1.4), and the simulation is reliable in the frequency range where the first-order perturbation theory results on both sides are reliable.
We anticipate that the simulation can be extended to a moving cavity that is not rigid, in the small amplitude regime commonly considered in the DCE literature [3,4,5], by equating Lη 1 (t) (respectively Lη 2 (t)) to the variation in the position of the left (right) boundary. We have however not examined this question systematically.

Conclusions
In this paper we have analysed the evolution of a quantum field in 1+1 dimensions under time-dependent Robin boundary conditions that occur in superconducting microwave circuit experiments in the high plasma frequency limit. We solved the evolution explicitly to linear order in the time variation of the Robin boundary condition, in a formalism that allowed us to handle both a semiopen waveguide and a closed waveguide, and for the latter allowing the boundary conditions at the two ends to be varied independently. For frequencies that are much smaller than the effective inverse length associated with the SQUID(s) at the end(s) of the waveguide, we verified that a suitable modulation of the SQUID(s) allows the waveguide to simulate a Dirichlet boundary condition at a mechanically moving boundary of the DCE, even in the regime where the mechanical motion is relativistic; both the experiment reported in [6] and the experimental proposals of [8,9,10] appear to operate within in this domain by a margin of approximately one order of magnitude. For higher frequencies the waveguide still exhibits particle creation and mode mixing, but these can no longer be quantitatively matched to those of the DCE, and in particular the large frequency falloff properties in the evolution are qualitatively different. These features in the large frequency Bogoliubov coefficients result in differing particle emission spectrum for the moving and non-moving systems when the effective length is larger than the inverse driving frequency L 1/ω d . In this limit, mechanically moving boundary radiation can be characterised as having a larger total flux and a less steep falloff at high frequencies compared to radiation from the static waveguide with time-varying Robin boundary conditions.
While the analogy between a moving Dirichlet boundary and a time-varying Robin boundary condition is useful for simulation purposes it is important to keep in mind that the physical systems corresponding to these two situations are different and therefore can lead to different outcomes. On the one hand, our results support proposals to simulate in a mechanically static waveguide quantum phenomena due to motion, including entanglement generation and degradation, even in a regime where the mechanically moving system experiences significant relativistic effects [8,9,10]. On the other hand, our results demonstrate that the interpretation of a waveguide experiment in terms of the simulation of the DCE is possible only in certain parameter ranges. Within the Robin boundary condition approximation that we have analysed, the range of DCE interpretation is determined just by the effective inverse length scale set by the SQUID(s) at the end(s) of the waveguide. It remains a subject to future work to establish the range of DCE interpretation when all relevant effects beyond the Robin boundary condition approximation are considered [11,12], including the effects due to a large but finite SQUID plasma frequency analysed in [13,14,15].

A Perturbative Bogoliubov identities
In this appendix we record identities satisfied by a perturbative expansion of the Bogoliubov coefficients. These identities are used in the main text for consistency checks of the perturbative treatment.
Let α and β denote the matrices of a Bogoliubov transformation in the notation of [27], with the matrix elements α jk and β jk . The indices may be discrete or continuous; in the latter case the matrix is understood as the kernel of an integral operator and may include a distributional part. By construction, the matrices satisfy the Bogoliubov identities where I stands for the identity matrix. Suppose now that α and β have the formal power series expansions where is a real-valued expansion parameter. Substituting (A.2) in (A.1) and collecting terms order by order, order 0 is identically satisfied while orders 1 and 2 give (A.3b) When α and β are real, (A.3) simplifies to 1 :

B Accelerated boundary in Minkowski spacetime
In this appendix we consider a massless scalar field in (1 + 1)-dimensional Minkowski spacetime subject to the Dirichlet boundary condition at one accelerated boundary, in the limit where the acceleration is treated perturbatively but may have arbitrary timedependence, and the velocity and travel distance of the boundary remain arbitrary. We follow the methods that were developed for a mechanically rigid accelerated cavity in [19,20]. The results overlap with those in the DCE literature [3,4,5] for small amplitude oscillations in the common domain of validity.

B.1 Inertial boundary to uniformly accelerated boundary
We work with a massless scalar field φ in (1 + 1)-dimensional Minkowski spacetime, in the notation of the main text: the metric is ds 2 = −dt 2 + dx 2 , and the wave equation is (2.1) with µ = 0. For t < 0, we take the field to live in the half-space x ≥ a −1 , where a is a positive constant of dimension inverse length, and we adopt at x = a −1 the Dirichlet boundary condition. We adopt the basis of mode functions where k > 0. φ k has the positive frequency k with respect to ∂ t , and the normalisation in the Klein-Gordon inner product is (φ k , φ k ) = δ(k − k ). For t ≥ 0, we make the boundary follow the uniformly-accelerated worldline x = √ t 2 + a −2 . The proper acceleration of the boundary is equal to a. The field lives in the region x ≥ √ t 2 + a −2 , and we take the field to satisfy the Dirichlet boundary condition at the accelerated boundary. We adopt the basis of mode functions where k > 0. Φ k has the positive frequency k/a with respect to the boost Killing vector x∂ t + t∂ x , and it has the positive frequency k with respect to the proper time of the boundary. The normalisation in the Klein-Gordon inner product is At the junction t = 0, we match the two sets of modes by the Bogoliubov transformation From the inner products that give the Bogoliubov coefficients [27], we find

B.2 Small acceleration expansion
We wish to find the asymptotic form of o α k k and o β k k (B.4) when the acceleration a of the boundary is small compared with both k and k . The small a expansion of o β k k may be obtained by applying to (B.4b) the method of repeated integration by parts [30]. The result is The small a expansion of o α k k is more involved since we expect the coefficients in this expansion to be no longer functions but distributions, the leading term being δ(k − k ). We shall not give a rigorous treatment but proceed heuristically as follows.
Starting from the integral in (B.4a), we replace k in the integrand by k + i where > 0. If is considered fixed, the asymptotic small a expansion can be obtained the method of repeated integration by parts [30], with the result Each of the three terms shown in (B.6) has a well-defined distributional limit as → 0, as follows from the identity and its derivatives. Assuming that the → 0 limit commutes with the small a expansion, we obtain where the distributions G 1 and G 2 are given by Both G 1 and G 2 are hence representable by a function except at the coincidence limit, and we may write As a consistency check, the expansions in (B.5) and (B.8) satisfy the linear order Bogoliubov identities (A.4a). The quadratic order Bogoliubov identities (A.4b) would require a distributional treatment and we shall not analyse them here.

B.3 Arbitrarily-accelerated boundary
Let τ denote the proper time of the boundary and a(τ ) the proper acceleration of the boundary, such that positive (negative) values of a(τ ) mean acceleration towards increasing (decreasing) x. We assume that a(τ ) is vanishing outside the interval τ 0 ≤ τ ≤ τ f and non-negative within this interval, and we assume that a(τ ) remains much smaller than the frequencies to be considered.
We define the early (respectively late) time modes by (B.1) in the inertial frame in which the boundary is at rest in the early (late) times. Proceeding as in [19,20], or in Section 4 of the present paper, we find that the Bogoliubov coefficients from the early time modes to the late time modes are and we omit the examination of the distributional part ofÂ k k at k = k . The above treatment assumes a(τ ) to be non-negative. We shall not examine the validity of (B.12) when a(τ ) may take negative values.

C Entanglement formula for continuous spectra
In this appendix we derive the formula for the perturbative approximation to the bipartite mode entanglement of the field, as measured by the entanglement negativity, in the case when the field solutions have continuous spectra. These generalise the negativity formulas found for field solutions with discrete eigenvalues in [31]. The field is prepared initially in the vacuum state and is subjected to an evolution which can be described by Bogoliubov transformations that are assumed to take the general form: where η 1, k and k are continuous real-valued parameters,Â k k andB k k are of linear order in η and c(k) is a phase taking the general form, c(k) = e if (k) , where f is a real-valued function of k. Furthermore, we assume thatÂ k k andB k k satisfy the conditions:Â where the star denotes complex conjugation. Equations (C.2) are satisfied by the evolution in the semiopen waveguide (4.1) and by the evolution with the accelerated mirror (B.12).
The added difficulty of systems with continuous spectra is that the states are Dirac δ-normalised which can lead to apparent infinities in formulas for the entanglement if not correctly handled. The infinities are an artefact of working with infinitely precise frequencies which in practice is not possible -there is always a spectral line-width, ∆k, associated with the measurement of a frequency. Experiments which measure bipartite entanglement of the field do so by measuring two distinct frequencies each with a small line-width. For simplicity, we will assume that the modes measured are uniform wavepackets of frequencies having spectral line-width ∆k and central frequencies k p , k p respectively. We define: otherwise.
We will also assume that the measured frequencies are sufficiently separated such that there is no overlap of their supports in frequency space: Letâ(k) be the annihilation operator associated with the frequency k. After the evolution the annihilation operators are transformed by the Bogoliubov transfomations according to:â The quadrature operators associated with these frequencies are: where we use the conventions of [32]. As already alluded these frequencies are not measured precisely rather the actual measurement occurs over a small band of frequencies. We can define the quadrature operators associated with a wavepacket centred at the frequency k p by: x kp = f kp (k)x(k)dk, (C.7a) p kp = f kp (k)p(k)dk. (C.7b) A short calculation shows that the quadrature operators satisfy the commutation relations: x i ,p j = 2iδ ij , (C.8a) x i ,x j = p i ,p j = 0, (C.8b) where i, j ∈ {k p , k p }.
Since the Bogoliubov transformations are linear and the vacuum state is a Gaussian state, it follows that the final state of the field will also be a Gaussian state. Gaussian states are completely characterised by their first and second statistical moments. It is the second moments which are important for determining the amount of entanglement in a Gaussian state. The second moments can be arranged into a covariance matrix: we define R = (x kp ,p kp ,x k p ,p k p ), then the covariance matrix can be defined by: where curly braces denote anti-commutator and the covariance matrix is normalised such that the covariance matrix of the vacuum state is the identity matrix. Expectation values are to be taken with respect to the initial state which in our case is the vacuum state. It is easy to see that the second term in (C.9) vanishes. We define the kernels: where superscript T indicates matrix transposition and σ a ≡ f kp (k)f kp (k )S(k, k )dkdk , (C.13a) Using relations (C.1) and (C.2), the matrix S(k, k ) simplifies to S(k, k ) = δ(k − k ) + S 1 (k, k ) + O(η 2 ), (C.14) where the matrix elements of S 1 (k, k ) are given by: We can also write the covariance matrix as a Taylor expansion in η: σ = σ 0 + σ 1 + O(η 2 ), (C. 16) and with equation (C.14) it is easily seen that σ 0 is the 4 × 4 identity matrix.
The symplectic values ofσ can be found from the absolute values of the eigenvalues of the matrix Ωσ, where Ω = diag(ω, ω) and The zeroth order eigenvalues are ±i and both eigenvalues have double degeneracy. It is therefore necessary to use degenerate perturbation theory to determine the eigenvalues to the linear order. The linear order corrections to the eigenvalues are found to be: I, (C. 19) where I ≡ dl dl dk dk f kp (l)f k p (l )f kp (k)f k p (k )c(l )c(k) β ll β k k . (C.20) The symplectic eigenvalues are therefore: and the negativity is: whereν s is the smallest of the two symplectic values. Let us now assume that the wavepackets are very sharply peaked ∆k/k p 1 and ∆k/k p 1, then the integrals in equation (C.20) can be approximated by: (C. 23) and the negativity (to first order in η) simplifies to: This is the continuous spectrum formula for the entanglement negativity and holds for non-integer values of the frequencies. It is equivalent to the formula for discrete frequency modes in the case when the modes have opposite parity, i.e., when (k + k ) is odd [31].