Shaping Dynamical Casimir Photons

Temporal modulation of the quantum vacuum through fast motion of a neutral body or fast changes of its optical properties is known to promote virtual into real photons, the so-called dynamical Casimir effect. Empowering modulation protocols with spatial control could enable to shape the spectral, spatial, spin, and entanglement properties of the emitted photon pairs. Space-time quantum metasurfaces have been proposed as a platform to realize this physics via modulation of their optical properties. Here, we report the mechanical analog of this phenomenon by considering systems whose lattice structure undergoes modulation in space and in time. We develop a microscopic theory that applies both to moving mirrors with modulated surface profile and atomic array meta-mirrors with perturbed lattice configuration. Spatio-temporal modulation enables motion-induced generation of steered frequency-path entangled photon pairs in co- and cross-polarized states, as well as vortex photon pairs featuring frequency-angular momentum entanglement. The proposed space-time dynamical Casimir effect can be interpreted as an induced dynamical asymmetry in the quantum vacuum.


Introduction
The generation of photon pairs out of the quantum vacuum, the so-called dynamical Casimir effect (DCE), was originally described as a motion-induced phenomenon [1], but it can occur when any kind of temporal modulation is exerted on the vacuum to promote virtual into real photons [2][3][4][5]. Motion-induced photon generation has not been observed to date because it requires unfeasibly large modulation frequencies of a mechanical boundary, and several analog DCE systems have been demonstrated involving temporal modulation of material properties [6][7][8][9]. Still, the physics of motional dynamical Casimir effects offers interesting insights into the interplay between matter and field fluctuations in non-equilibrium systems. Motional DCE (also known as motion-induced or mechanical DCE) is typically described in a "field-centric" approach based on quantum fluctuations of the electromagnetic field supplemented with time-dependent boundary conditions. A "matter-centric" approach has been recently pursued based on microscopic models that emulate a moving dielectric mirror as a collection of accelerated dipoles that emit quantum radiation [10]. Interestingly, the two descriptions result in identical predictions for the angular emission profile of DCE photons. This duality between field-centric and matter-centric approaches also occurs in equilibrium fluctuation-induced interactions [11]. Microscopic models can be used to study the dissipative counterpart of DCE emission, namely the drag force on the moving mirror, as well as the related problem of quantum friction and associated near-field DCE emission of surface polaritons [12,13].
Modulation protocols having both temporal and spatial control can enable novel functionalities such as the generation of complex structured Casimir light. We have recently proposed an analog dynamical Casimir effect based on the concept of space-time quantum metasurfaces [14], in which the optical properties of a quantum metasurface are modulated in space and time and generate DCE photon pairs with tailored spatial profiles at giant production rates.

arXiv:2105.04510v1 [quant-ph] 10 May 2021
Here, we develop the motion-induced version of this effect based on a microscopic model of a spatio-temporally modulated mirror. We describe the mirror as a collection of dipoles whose center-of-mass coordinates are modulated in space and in time, and find that photon pairs are generated out of vacuum with tailored spatial modes determined by the spatial modulation protocol. The collection of dipoles behaves as a quantum phase array antenna that emits structured Casimir light. The same microscopic theory applies to a quantum meta-mirror comprising an atomic array with modulated lattice structure [15]. We call the proposed effect space-time DCE to distinguish it from the standard DCE process that involves temporal modulation only.

Microscopic model for space-time motional DCE
A microscopic model for a spatio-temporally modulated dielectric mirror consists of an array of N multi-level atoms, each of which undergoes its own accelerated trajectory R j (t) and couples to the electromagnetic field via the dipolar interaction. The trajectories are driven by some external agent and in this work we consider oscillatory trajectories typically employed in studies of the dynamical Casimir effect. We note, however, that our theory can be extended to arbitrary trajectories. The spatial modulation is incorporated via a synthetic phase Φ(R j ) that is imprinted by a temporal delay of the oscillation of each atom. We thus consider where R j is the static position of the j-th atom (j = 1, ..., N), ∆ is the oscillation amplitude, and Ω the oscillation frequency. A conceptual representation of the spatio-temporally modulated atoms is shown in Fig. 1. For example, a traveling-wave modulation on the x − y plane corresponds to Φ(R j ) = β · r j with r j = (x j , y j ) and β = (β x , β y ) a momentum "kick" on the same plane. A spinning-wave modulation around an axis orthogonal to the mirror's plane corresponds to Φ(R j ) = ϕ j , where is an integer denoting an imprinted angular momentum and ϕ j is the azimuthal coordinate of each atom. Note that under this modulation the atoms are not rotating, just the temporal dephasing of their oscillatory motion along the z direction is spinning around the z− axis. The traveling-wave modulation could be implemented by an acoustic perturbation that launches a plane-wave on the surface of the mirror. In the case of a meta-mirror formed by an atomic array monolayer trapped in a three-dimensional deep optical lattice, the modulation could be accomplished by shaking the vertical lattice and using spatial light modulators to set a phase shift among the vertical trapping fields on different atoms to imprint the linear synthetic phase. The dipole array interacts with the electromagnetic field via the Hamiltonian (SI units used throughout) where d j (t) is the electric dipole operator of atom j, E(R, t) and B(R, t) are the electric and magnetic fields both evaluated at the atomic positions, and v j (t) = −ẑΩ∆ sin(Ωt − Φ(R j )) is the (non-relativistic) velocity of each atom. In this paper the atomic internal degrees of freedom and the electromagnetic field are treated as quantum dynamical variables, while the atomic center-of-mass is a classical prescribed motion. The first term in the Hamiltonian is the usual electric dipole interaction and the second is the so-called Röntgen interaction. We note that the Röntgen term must be taken into account in the DCE far-field emission of photons, as we study here. For the DCE problem of emission in the near-field, e.g., pairs of surface polaritons induced on a surface by a close-by moving atom, it can be neglected [12,13]. The non-zero matrix elements of the electric dipole operator d j (t) in the interaction picture are where |g j is the ground state of atom j, |e j is an excited state,η is a real unit vector denoting the orientation of the j-th dipole, and d eg is the matrix element. Note that all atoms have the same d eg as we assume identical atoms. When the atoms are isotropic, the sum over orientations of any given atom gives ∑η j (η j ) a (η j ) b = δ ab . We consider that the atoms are sufficiently spaced within the array to neglect photon multi-scattering among different atoms, rendering the evolution of different atoms identical irrespective of their location within the array. The electromagnetic field in the interaction picture is then simply given by the free field, that we expand in a set of modes in the usual way Here, V is a quantization volume, E γ (R) is the spatial mode, ω γ > 0 is the mode frequency, and a † γ , a γ are creation and annihilation operators of photons in mode γ. The specific choice of modes is determined by the symmetries of the synthetic phase, as we discuss later in the paper.

Two-photon emission rate
We assume that the initial state of the N-atom system plus electromagnetic field is all the (identical) atoms in their ground state and the field in vacuum, i.e., |ψ(0) = |{g}; vac , where we denote |{g} = |g 1 , g 2 , ..., g N the multi-atom ground state. The time-evolved state in the interaction picture to second order in the dipolar couplings d eg is The first term is a correction to the amplitude of the initial state due to the coupling. The second term corresponds to one photon emitted and the atomic array being in an equal superposition of one atom excited in level e j and all the rest in the ground state. The third, fourth, and fifth terms are two-photon contributions in which two, one and no atom is excited, respectively. We denote |e j , {g} j = |g 1 , ..., g j−1 , e j , g j+1 , ..., g N the state in which one atom is excited, and |e j , e j , {g} j,j = |g 1 , ..., g j−1 , e j , g j+1 , ..., g j −1 , e j , g j +1 , ..., g N the state in which two atoms are excited. We make the crucial assumption that the modulation frequency is much smaller than any atomic transition frequency from the ground state, Ω ω eg . In this approximation, the second, third, and fourth terms give transition amplitudes that are non-resonant and their contributions to the process of one-and two-photon generation is negligible and are discarded in the following. In this regime, the DCE process is solely given by the fifth term, for which photon pairs are generated and all atoms remain in their ground state.
One can compute the transition amplitude c γ 1 ,γ 2 (t) in two ways: (a) use second-order time-dependent perturbation theory based on the bilinear Hamiltonian H(t) that depends on field and atomic degrees of freedom, or (b) use first-order time-dependent perturbation theory based on an effective Hamiltonian H e f f (t) that depends quadratically on field degrees of freedom and contains the information about the atoms through their ground-state polarizability. Approach (b) was used in [10] for the case of a single atom and no synthetic phase: the atom-field interaction is written in the atom's comoving frame as the standard static dipolar interaction, it is then re-written as an effective Hamiltonian that traces over the atom's internal degrees of freedom rendering the interaction quadratic in the field and depending on the atomic ground-state dynamic polarizability, and finally the Hamiltonian is Lorentz boosted to the lab frame to get H e f f (t). Although approach (b) is technically simple, it has the drawback that it obscures the joint atom-field dynamics. Furthermore, for multiple atoms and non-zero synthetic phases it requires to introduce multiple comoving frames because atoms in the array have different instantaneous velocities. For these reasons, in this work we prefer to follow the physically more transparent albeit more cumbersome approach (a).
The transition amplitude is given by the usual expression in second-order time-dependent perturbation theory that sums over intermediate virtual states |e j , {g} j ; γ . In the limit when the relevant field wavelengths are much larger than the atomic displacements, ω∆/c 1, we can approximate the Röntgen term of the Hamiltonian evaluating the fields at the static atomic positions, so that it has a time-dependency ∼ sin(Ωt − Φ(R j )) arising from the velocity. Also, the dipolar term in the Hamiltonian can be expanded around the equilibrium positions, giving a time-dependency ∼ cos(Ωt − Φ(R j )) from the spatial gradient of the field. After these approximations, the concatenated time integrals in Eq. (6) can be readily evaluated as they are reduced to integrals of products of simple harmonic functions. The final result is a sum of various terms that oscillate in time. We drop the co-rotating term e i(ω 1 +ω 2 +Ω)t because it averages out to zero. Counter-rotating terms of the form (±ω 1,2 ± Ω ∓ ω eg ) −1 × e i(±ω 1,2 ±Ω∓ω eg )t and (±ω 1,2 ± ω eg ) −1 × e i(±ω 1,2 ±ω eg )t could potentially produce a resonant effect if their denominator vanished. However, under the assumption that the modulation frequency is much smaller than the atomic transition frequencies from the ground state, Ω ω eg , this cannot occur (note that DCE photons necessarily have frequencies ω 1,2 ≤ Ω per energy conservation). The only counter-rotating term contained in the transition amplitude that can produce a resonant effect is (ω 1 + ω 2 − Ω) −1 × e i(ω 1 +ω 2 −Ω)t . Upon using all these considerations, the transition amplitude results where we already performed the sum over orientations of the isotropic dipoles. Here, Note that the factor involving the summation over excited states is the ground state static polarizability of an isotropic atom, α 0 = ∑ e (2d 2 eg /3 0h ω eg ). As expected, this expression for the transition amplitude is identical to the one that results from the effective Hamiltonian approach and first-order time-dependent perturbation theory. As already mentioned, the merit of the presented approach is that it highlights atom-field dynamical processes that remain obscured in the approach that starts out from polarizabilities. The rate of production of photon pairs in modes γ 1 , γ 2 is obtained taking the long time limit r which is to be understood within time-dependent perturbation theory: t is typically not longer than a fraction of the relevant atomic life times. Finally, The Dirac delta ensures energy conservation in the DCE process and the last factor contains other conservation laws that depend on the synthetic phase, as described below. This factor has the form of an array form factor akin to classical antenna theory: Each atom emits with a phase e −iΦ(R j ) weighted by W γ 1 ,γ 2 (R j ). Interestingly, the atomic array is like as driven quantum antenna emitting DCE photon pairs.

Linear synthetic phase
In the case Φ(R) = β · r, we quantize the electromagnetic field using plane-wave modes labelled by γ = {K, λ}, where K = (k, k z ) is the wave-vector with k its projection on the x − y plane, k z the normal projection, and λ the polarization state. The field modes are E is the dispersion relation. In the continuum limit for the momenta ∑ K → V(2π) −3 dK, the two-photon emission rate into modes K 1 , whereW withK 1,2 = K 2 /K 1,2 and K 1,2 = ω 1,2 /c. Note thatW K 1 ,λ 1 ;K 2 ,λ 2 =W K 2 ,λ 2 ;K 1 ,λ 1 . The rate is equal to the one for a single oscillating atom (first line in Eq. (10)) multiplied by a multi-atom correction (second line) that has the form of a product of two array form factors. In principle, each atom emits independently of its neighbors, resulting in incoherent emission of DCE photon pairs that should scale linearly in the number of atoms in the array. However, since EM quantum fluctuations have all possible wavelengths, there are large-wavelength fluctuations that coherently couple to all atoms and coherent emission à la super-radiance should be possible. In this case, the rate should scale as the square of the number of atoms, which we now show that is indeed the case. To get a close and simple expression for the form factors, we assume atoms are arranged into a finite size cubic array with inter-atom distance d. We choose the coordinate system so that the static positions of the atoms are R j = dm xx + dm yŷ + dm zẑ (1 ≤ m i ≤ N i with dN i = L i the size of the array in each direction). We compute the modulus square of each of the summations in Eq. (10) and express the result as AF 1 (k 1 , k 2 ; β) AF 2 (k 1z , k 2z ), where the array form factors are Both array form factors are periodic functions of their arguments. In the N x , N y N z limit, which mimics a finite-width large-area slab, the first array factor is approximately equal to a two-dimensional comb of sharp peaks, all of equal height and proportional to the square of the total number of in-plane atoms, (N x N y ) 2 . This indicates that in the large-N limit the emission is coherent. Figure 2 shows how the momentum kick controls the directionality of the emitted DCE pairs. In the figure, the emission direction of one of the photons is fixed to be vertical (k 2 = 0) and we trace over its polarization state. The emission lobes of the other photon are plotted for different values of β. For each β, maximal emission is when k 1 = β. The emission rate for TM photons is always larger than for TE photons.
For an infinite atomic array in the x − y plane and finite along the z direction, i.e., mimicking a finite-width infinite-area slab, one can calculate the summation over in-plane atomic positions taking advantage of periodicity. In this case we use the lattice summation identity where the sum over q corresponds to reciprocal momenta. Here, A is the (infinite) area of the slab and n S = N x N y /A = 1/d 2 is the number surface density of atoms. When all relevant wavelengths are much larger than the inter-atomic distance d, one can take the continuous limit approximation in which only q = 0 survives the sum over q. Non-zero values of the reciprocal momentum correspond to high-order diffraction modes that are evanescent and do not contribute to the generation of DCE photons. Then we obtain the in-plane linear momentum conservation condition The time-evolved quantum state can be written as and corresponds to a color-path entangled superposition.
In the absence of kick (β = 0), the angular spectrum of the emitted photons is always symmetric with respect to the normal of the slab, k 1 = −k 2 , which is what happens in the standard DCE problem of a rigidly oscillating mirror [10]. In stark contrast, when β = 0, the traveling-wave modulation generates directed ripples on the mirror and produces photons that are emitted asymmetrically. This steered DCE emission can be interpreted as a modulationinduced asymmetry of the quantum vacuum.
The DCE emission rate is obtained by integrating out one of the photons in the emitted pair, r (β) . For the finite-width infinite-area atomic array, the inplane momentum k 2 is fixed by the momentum conservation condition k 2 = β − k 1 , and the out-of-plane momentum is also fixed by the dispersion relation and the energy conservation condition, k 2z = ζ 2 [(Ω − ω 1 ) 2 /c 2 − |β − k 1 | 2 ] 1/2 , where ζ 2 = ±1 gives the two possible emission directions normal to the array (z > 0 or z < 0). Hence, the momentum integration above can be done straightforwardly and only the summation over polarization states remains. Note that propagative DCE photons can be emitted only when k 1z = ζ 1 (ω 2 1 /c 2 − |k 1 | 2 ) 1/2 and k 2z are real. The spectral photon emission rate per unit area for photons with polarization λ 1 , in-plane momentum in the interval (k 1 , k 1 + dk 1 ), out-of-plane direction ζ 1 , and frequency in the interval (ω 1 , ω 1 + dω 1 ), takes the form dΓ (β) where f (β) We perform the calculations in the linear polarization basis. The emitted photons can have both TE or both TM polarization, and emission of photons with non-collinear in-plane momenta can be cross-polarized becausẽ In space-time DCE this can be non-zero since k 2 = β − k 1 , while cross-polarized emission is not possible in standard DCE because photons have collinear in-plane momenta, k 1 = −k 2 . Figure 3 depicts f (β) ζ,λ (k, ω) for the case of an infinite monolayer, λ = TE polarization and ζ = 1 (by symmetry, the plots for ζ = −1 are identical). Polar density plots are shown both for the high-frequency (ω 1 > Ω/2) and low-frequency (ω 2 < Ω/2) photons in the emitted pair. In the absence of kick, the high-frequency photon can be emitted in any azimuthal direction but it has a maximum polar angle of emission (panel a, area to the right of the central vertical solid line), while no such a constraint exists for the low-frequency photon (panel b, area to the right of the central vertical solid line). As the magnitude of the momentum kick β increases, the distributions undergo intricate changes. The region of allowed emission for the first photon gets deformed when the kick is non-zero and at a critical value cβ = Ω − ω 1 (between panels d and e) an "island" of emission appears surrounded by a sea of forbidden emission directions (shaded areas). The island drifts to higher polar angles until it touches the grazing emission line when cβ = 2ω 1 − Ω (between panels e and f). The island starts to shrink in size (panels f-j), and finally at β max = Ω/c it collapses to a point (after panel j, not shown) and the photon is only emitted parallel to the kick. Far-field emission above that value of the kick is not possible. Regarding the second photon, its emission distribution remains mostly unperturbed until at cβ = Ω − 2ω 2 two areas of forbidden emission appear at large polar angles and opposite to the kick direction (between panels e and f). The forbidden region grows until it engulfs its allowed emission region and a second island forms at cβ = Ω − ω 2 (between panels h and i). Finally, it ends up being emitted at a grazing angle but in a direction anti-parallel to the kick (after panel j, not shown). The modulation also excites hybrid entangled pairs composed of one photon and one evanescent surface wave (shaded areas), and when β > β max only evanescent modes are created on the atomic monolayer and subsequently decay via non-radiative loss mechanisms.
The DCE spectral emission rates for photons with polarization λ are obtained by integrating over all emission directions and depends only on the modulus of the kick. The TE and TM spectral emission rates are shown in Fig. 4 as a function of frequency and momentum kick. At zero momentum kick, they are both symmetric with respect to ω = Ω/2 and the rate of TE emission is smaller than TM emission (note the different vertical scales in the plots). The two rates are identical to those derived in [10] for the standard DCE problem in the absence of synthetic phase. At non-zero kick, they both become slightly asymmetric, with the peak of the TE (TM) moving to frequencies larger (smaller) than ω = Ω/2. The origin of the asymmetry is the non-zero cross-polarized emission. Indeed, the TE spectral emission rate has contributions from two terms,W 2 K 1 ,TE;K 2 ,TE andW 2 K 1 ,TE;K 2 ,TM . The first term is symmetric around ω = Ω/2 because upon interchanging the frequency, momentum, and spin indices one getsW 2 K 1 ,TE;K 2 ,TE = W 2 K 2 ,TE;K 1 ,TE . However, the second term is not symmetric because the same interchange gives W 2 K 1 ,TE;K 2 ,TM =W 2 K 2 ,TM;K 2 ,TE =W 2 K 2 ,TE;K 2 ,TM . The same argument applies to the TM spectral emission rate. The spectral emission rate summed over polarizations is symmetric, though.
The total rate per polarization is obtained by integrating over frequency As shown in Fig. 5, the TM rate shows a monotonous decay to zero at β max = Ω/c, while the TE rate decays non-monotonically. Furthermore, Γ TE (β) < Γ TM (β) and their difference diminishes as the kick grows. The figure also shows the emission rate for circularly polarized photons, which is the same for right-and left-circular polarization, Γ R (β) = Γ L (β), and is initially flat and then decays monotonically to zero. As follows from the figure, Γ TE (β) < Γ R/L (β) < Γ TM (β) and the total rate verifies Γ(β) = Γ TE (β) + Γ TM (β) = 2Γ R/L (β).

Spinning synthetic phase
In this section we briefly discuss the case of a spinning synthetic phase Φ(R) = ϕ. Because of the symmetry properties of the phase, one needs to quantize the electromagnetic field with angular momentum. Usually, this is done in the paraxial approximation in terms of Laguerre-Gauss modes that carry orbital angular momentum [16]. However, for DCE this is not appropriate because the emitted photon pairs are non-paraxial, and a more general approach is required. We follow the quantization scheme based on vector-Bessel modes, that form a complete basis for the electromagnetic field with angular momentum and do not require any paraxial approximation [17]. In this scheme, neither orbital angular momentum (OAM) nor spin angular momentum (SAM) are good quantum numbers, only their sum is. Vector-Bessel modes are labelled by γ = {k, k z , η, m}, where k ≥ 0 is the transverse linear momentum, k z is the axial linear momentum, η = ±1 is associated with the sign of the transverse spin s t = ±hck z /ω, and m is the total angular momentum (OAM plus SAM). The dispersion relation ω 2 /c 2 = k 2 + k 2 z = K 2 does not depend on m or η. In cylindrical coordinates, the vector-Bessel mode is E where J m (x) are Bessel functions. Note these modes are non-diffracting (constant amplitude along the z direction), have a topological vortex singularity along the z-axis with a phase wrapping equal to 2πm, and decay along the radial direction.
The two-photon generation rate is where we took the continuum limit ∑ k,k z → V ∞ −∞ dk z ∞ 0 kdk. The function T γ 1 ,γ 2 (ρ j ) results from plugging the vector Bessel mode into Eq. (9) and has a cumbersome expression that we do not write here (it will be written below after performing the summation over j). Emitted photon pairs can have the same (η 1 η 2 = 1) or opposite (η 1 η 2 = −1) signs of transverse spin. In contrast to the rate for the linear synthetic phase Eq. (10), it is not possible to express the rate Eq. (22) as that of a single atom multiplied by a multi-atom correction. The underlying reason is the nontrivial topology of the imprinted modulation, which is ill-defined for a single atom.
To perform the sum over atoms, we assume they are arranged into a cylindrical geometry of radius R and height L z , forming parallel disks along the z direction (inter-disk distance d) and concentric rings on the orthogonal plane (inter-ring separation r). The axis of the cylinder coincides with the axis of spinning modulation axis (z-direction). The second line in Eq. (22) is the product of two array form factors: one is the same AF 2 discussed in the previous section, and the other is a new AF 3 that depends on all quantum numbers k 1 , k 2 ; k 1z , k 2z ; η 1 , η 2 ; m 1 , m 2 of the two photons. We compute AF 3 by writing the radial coordinates of the atoms as ρ j = rn, with 0 ≤ n ≤ N R and rN R = R, and the azimuthal coordinates on the n−th ring as ϕ j (n) = (2π/N(n))p with p = 0, 1, ..., N(n) − 1. We assume the number of atoms on a given ring is proportional to the ring perimeter, hence N(n) = ξn with ξ a constant. In the limit ξ 1, the atoms form a quasi-continuum on the azimuthal direction and the summation over the azimuthal index p gives | ∑ n,p T(ρ j )e i(m 1 +m 2 − )ϕ j | 2 = δ m 1 +m 2 − | ∑ n N(n)T(rn)| 2 . Hence, the spinning modulation stirs the quantum vacuum and generates vortex photon pairs whose angular momentum must add up to the imprinted topological charge : in agreement with angular momentum conservation. For = 0 the generated photons in a pair have opposite twist. The two-photon state is frequency-angular momentum entangled, The remaining summation over the radial coordinates is performed in the limit of a continuum distribution of atoms, r → 0 and N R → ∞ such that rN R = R is constant. We define x = nr and approximate the sum as an integral over the continuous variable x. The array form factor takes the final form where n S is the number of atoms per unit of disk area. The total emission rate from the modulated atomic array is The spectral weight function f (ω) sums the product For = 0, the integrals H (±) and H (0) are one of Lommel's integrals and have a closed form, but the other three integrals do not. In the limit R → ∞, all the integrals are of the Weber-Schafheitlin form, ∞ 0 dyy q J µ (κy)J ν (κ y), with exponent q = 1 or q = 0. The above expressions allow to calculate f (ω, m) in space-time motion-induced DCE under spinning modulation. Direct inspection of f (ω, m) shows that it satisfies the condition f (ω, m) = f (Ω − ω, − m). This identity states the simple fact that, as photons are emitted in pairs satisfying energy and angular momentum conservation, the emission rate of a photon with given frequency and angular momentum must be the same as the emission rate at the complementary frequency and angular momentum. Also, one intuitively expects that one of the emitted photons should most probably be emitted with an angular momentum equal to the driving, i.e., m = . Because of angular momentum conservation, its companion in the pair is most probably emitted with m = 0. This intuition was confirmed in [14] for the case of DCE photons emitted from a quantum metasurface with spinning modulation of its optical properties.
Finally, we briefly discuss the structure of the energy density and Poynting vector in space-time DCE with spinning synthetic phase. Their expectation values on the evolved quantum state have a single vortex singularity along the z-axis for = 0. The reason why there is a single vortex line and not two (one per emitted photon) is due to the nature of an expectation value: it corresponds to summing over infinite many realizations of DCE two-photon emission events, and as such it is not surprising that individual dual vortex events average out and lead to a single vortex aligned with the modulation axis. Since there is only a single preferred axis in the problem, symmetry considerations explain why dual vortices cannot occur. A more formal justification is the non-diffracting nature of the vector-Bessel modes employed in the quantization scheme. However, a single experimental realization will have dual vortices along random directions of emission.

Discussion
As mentioned in the introduction, motional DCE involving temporal modulation of a boundary has not been observed to date because experimentally feasible mechanical modulation frequencies are insufficient for generating a large amount of photons per unit of time. For example, the rate of photon creation from an oscillating perfectly-reflecting mirror of area A is Γ DCE = AΩ 5 ∆ 2 /15(2π) 2 c 4 [18], which gives ∼ 10 −21 photons s −1 for a mirror of 1 cm 2 area, modulation frequency Ω/2π ∼ 1 MHz, and modulation amplitude ∆ ∼ 100 nm. This is less than one photon every 10 years. An alternative approach to moving the whole mirror is to make its surface oscillate by launching acoustic waves, a scheme that is related to the proposed space-time motional DCE. However, usual materials can bear maximal relative deformations δ max = ∆ max /(v s /ω w ) ∼ 10 −2 (v s is the speed of sound and ω w is the frequency of the acoustic wave) and have maximal velocities of the boundary v max ∼ δ max v s ∼ 50 m/s, which again leads to negligibly small photo-production rates [2]. The same limitation occurs in space-time motional DCE. According to Fig. 5, the maximal rate happens for zero synthetic phase and is We estimate it for experimentally feasible parameters of a spatio-temporally modulated metamirror consisting of an array of 87 Rb atoms (ground state polarizability α 0 = 5.9 × 10 −28 m 3 [19]) loaded into a 2D square optical lattice and vertically trapped by another shaking optical lattice. For a lattice constant d = 532nm, unit filling, N ∼ 200 atoms (A ∼ 50µm 2 , n S ∼ 4µm −2 ) [20], modulation amplitude ∆ ∼ 100 nm, and modulation frequency Ω/2π ∼ 10 kHz [21,22], the maximal rate is Γ max ∼ 10 −77 photons s −1 . Again, unrealistic large modulation frequencies would be required to get measurable DCE rates. In order to face these severe limitations of both standard and space-time motional DCE, analog DCE set-ups based on modulation of optical properties are required. For standard DCE, the experiment [6] based on a superconducting microwave coplanar waveguide is an analog of a one-dimensional DCE mirror, with the rate given by Γ = (Ω/12π)(v e f f /c) 2 , where v e f f is the effective velocity of the mirror. For the experimental parameters Ω/2π = 11 GHz and v e f f /c = 0.05, the estimated rate is Γ ∼ 10 6 photons s −1 . For space-time DCE, a possible analog of a DCE mirror moving in three dimensions is a set-up consisting of a graphene-disk metasurface whose electro-optical properties are spatio-temporally modulated. All-optical modulation of the Fermi energy at THz frequencies should enable giant rates Γ ∼ 10 12 photons s −1 from centimeter-sized metasurfaces [14].

Conclusions
The space-time dynamical Casimir effect offers a novel degree of control over photoproduction from the quantum vacuum. A space-dependent synthetic phase distribution imprinted on the temporal modulation protocol of optical properties or geometrical boundaries enables the generation of photon pairs with tailored spatial modes and entanglement properties. Travelling-wave modulations, like parallel ripples propagating on the surface of a mirror, generate steered photon pairs that are frequency-path entangled. Spinning-wave modulations, like twisting ripples on the mirror, produce vortex photon pairs featuring frequency-angular momentum entanglement. The synthetic phase can be reconfigured ondemand by changing the modulation protocol, allowing to modify the nature of entanglement between the generated photon pairs.
We end the paper with a conjecture. As discussed above, the synthetic phase distribution Φ(r) is imprinted onto space-time quantum metasurfaces via temporal delay of the modulation signal on different meta-atoms [14]. The meta-atoms can all have the same geometry, as the case of the atomic array meta-mirror, or the surface can even have no meta-atoms (unstructured), as is the case of a flat mirror. In both scenarios the spatial profile of the emitted DCE photons is controlled by the synthetic phase imprinted by the modulation. On the other hand, a phase distribution Ψ(r) can be imprinted onto a static surface by decorating it with meta-atoms having judiciously designed geometrical parameters, so-called gradient metasurfaces. For example, a blazed grating mirror or a meta-mirror with resonators of varying size. The phase distribution controls the spatial profile of light reflected from the static metasurface. An interesting question is: what is the nature of DCE emission when such gradient metasurfaces are set in motion along a direction normal to the metasurface plane? We conjecture that when Ψ(r) = Φ(r), the emitted Casimir light from the moving gradient metasurface has the same properties as the Casimir light generated by the space-time modulation on an unstructured surface. For example, the blazed grating mirror moving orthogonal to the grating plane should produce steered frequency-path entangled photon pairs, just as a flat mirror under a spatio-temporal modulation with a linear synthetic phase would. A test of this conjecture will require to calculate DCE emission from moving surfaces decorated with complex nanostructures, an overwhelmingly challenging analytical calculation and a very hard numerical task on its own. Were the conjecture proven to be correct, it would allow much simpler evaluation of DCE emission from moving structured bodies by simply considering unstructured bodies modulated with the appropriate synthetic phase.
Author Contributions: All authors equally contributed and have read and agreed to the published version of the manuscript.
Funding: This work was supported by the DARPA QUEST and LANL LDRD programs.

Data Availability Statement:
The data supporting this research is available from the corresponding author upon reasonable request.