Coherently amplifying photon production from vacuum with a dense cloud of accelerating photodetectors

An accelerating photodetector is predicted to see photons in the electromagnetic vacuum. However, the extreme accelerations required have prevented the direct experimental verification of this quantum vacuum effect. In this work, we consider many accelerating photodetectors that are contained within an electromagnetic cavity. We show that the resulting photon production from the cavity vacuum can be collectively enhanced such as to be measurable. The combined cavity-photodetectors system maps onto a parametrically driven Dicke-type model; when the detector number exceeds a certain critical value, the vacuum photon production undergoes a phase transition from a normal phase to an enhanced superradiant-like, inverted lasing phase. Such a model may be realized as a mechanical membrane with a dense concentration of optically active defects undergoing gigahertz flexural motion within a superconducting microwave cavity. We provide estimates suggesting that recent related experimental devices are close to demonstrating this inverted, vacuum photon lasing phase.


Introduction
One of the most striking consequences of the interplay between relativity and the uncertainty principle is the predicted detection of real photons from the quantum electromagnetic field vacuum by non-inertial, accelerating photodetectors. For the special case of a photodetector with uniform, linear acceleration a in Minkowski vacuum, the Unruh effect predicts that the detected photons are furthermore in a thermal state with temperature T = a/(2πck B ) [1][2][3][4][5]. However, to measure for example a 1 K thermal photon temperature, a uniform detector acceleration a = 2.47 × 10 20 m s −2 is required, which seems impossibly high for any current or planned tabletop experiments.
Given the difficulty in achieving uniform photodetector accelerations of sufficient magnitude and duration, certain experimental approaches have instead considered electrons that are accelerated nonuniformly by intense electromagnetic fields such that the trajectories are confined to a finite spatial volume. In particular, Bell and Leinass [6] interpreted the observed spin depolarization of electron(s) undergoing uniform circular motion in a storage ring in terms of effective heating by vacuum fluctuations (see also Refs. [7,8]). And Schutzhold et al. [9] showed that oscillating electrons emit both classical (Larmor) and quantum radiation, with the latter corresponding to the accelerating electrons converting quantum vacuum fluctuations into entangled photon pairs. Motivated by the goal to enhance the detection of photons from vacuum, a proposal by Scully et al. [10] locates an accelerating photodetector atom within an electromagnetic cavity; the idealized, perfectly conducting walls of the enclosing cavity effectively modify the spacetime geometry from Minkowki space to a finite volume for the electromagnetic vacuum [2]. This finite confining cavity volume then necessitates a non-uniform, i.e., oscillatory acceleration, similar to the accelerating electron proposal of Ref. [9]; the photodetector atom is envisaged for example as attached to a piezoelectrically driven, vibrating cantilever structure. Another recent relevant proposal [11] considers a rapidly rotating, initially excited atom in a cavity, such that the atom's spontaneous emission is enhanced when the sum of the rotation and atom splitting frequencies resonantly matches a cavity mode frequency.
In these investigations, the electrons and atom photodetectors accelerate nonuniformly (in either direction or magnitude) and in contrast to the Unruh effect as conventionally defined [3][4][5], the predicted photon spectrum is not expressible as a thermal distribution with temperature given solely in terms of the acceleration and fundamental constants [12,13]. Nevertheless, the effects of entangled photon pair production and photon detection from vacuum are still direct consequences of the non-inertial motion of electrons and atoms.
In recent work [14,15], we proposed close analogs of the vacuum, oscillatory photodetection effect that involve coplanar, i.e., two-dimensional (2D) superconducting microwave cavity circuit systems [16]. The photodetector is modelled alternatively as a harmonic oscillator [14] and as a qubit [15] that capacitively couple to the 2D cavity via a mechanically oscillating film bulk acoustic resonator (FBAR) [17]; by increasing the FBAR capacitor plate area, a strong detector-cavity coupling can be achieved. Furthermore, when the FBAR mechanical frequency matches the sum of the cavity mode and detector ground-first excited energy level splitting, the photon production rate from cavity vacuum is resonantly enhanced [10] so as to be detectable. The use of an FBAR introduces an actual mechanical acceleration, in contrast to previous circuit-based analogs where the accelerating photodetector is mimicked by an electromagnetically induced, time changing detector coupling [18,19].
However, a feasible way to demonstrate "genuine" oscillatory acceleration detection and production of photons from vacuum, where the photodetectors are accelerating within a three-dimensional (3D) electromagnetic cavity (as opposed to the above 2D analogues) is still lacking. Despite the resonant gain resulting from placing an accelerating atom inside a cavity [10,11], currently achievable cavity quality factors, atom-cavity coupling strengths, and mechanical oscillation acceleration magnitudes are insufficient for attaining measurable photon production from vacuum for a single accelerating detector.
In this paper, we describe a way to realize acceleration photodetection from vacuum with a scheme that involves "many" oscillating photodetectors contained within a microwave cavity (Fig. 1). Such a scheme may be realized by exploiting optically active defects that are embedded within a vibrating membrane structure, such that the defects effectively function as accelerating photodetectors. As one concrete example, we consider a gigahertz (GHz) vibrating diamond membrane containing nitrogen vacancy (NV) color centre defects [20]. We shall take advantage of a Dicke superradiance-like [21][22][23][24], vacuum amplification effect due to the collective motion of the many accelerating photodetectors. While the possibility has been mentioned before to coherently enhance photon production from vacuum using many accelerating electrons [9] or detectors [25,26], no such investigations have been carried out previously to the best of our knowledge. (Note, however, the proposal of Ref. [27] to amplify and hence detect cavity photons produced via the dynamical Casimir effect for an accelerating cavity mirror through their superradiant-enhanced interaction with a cloud of optically pumped atoms.)

Results
Cavity-N accelerating detectors model We first introduce our model for a cavity field interacting with N 1 accelerating two level system (TLS) photodetectors and show how the model maps by approximation onto a simpler, parametrically driven Dicke-type Hamiltonian [28,29]. Referring to Fig. 1, the N TLS detectors are assumed to be embedded within a two-dimensional membrane material system undergoing driven, small amplitude GHz flexural vibrations in the x-coordinate direction (i.e., displacements normal to the static, equilibrium membrane y-z surface), hence imparting acceleratory motion to the TLSs' centres of mass; the GHz membrane flexural motion may be actuated via piezoelectric transducers, for example. The membrane is enclosed by a microwave cavity 'box' with volume V = L x L y L z and box dimensions L i ∼ few cm. The cavity dimension in the flexural motion (x coordinate) direction is denoted as L x , and the membrane is assumed to be located at the cavity midway point x = L x /2. We consider the following starting model action for the combined cavity field-N detector system: where the detectors, classically modelled by internal, one dimensional anharmonic oscillator coordinates Q i , i = 1, . . . N , are linearly coupled to a 3 + 1 dimensional massless scalar field denoted as Φ(t, r) that models the electromagnetic field within the cavity. Equation (1) generalizes the starting action of Ref. [30] to include more than one detector, as well as assumed weakly anharmonic oscillator potential energy terms gQ 4 i /4!. The latter allows the detectors' internal degrees of freedom to be truncated to TLSs under parametric resonance conditions; while the harmonic oscillator detector model is beneficial for analytical calculations [30], it is somewhat artificial since no detector is purely harmonic, and in particular results in an unphysical parametric instability beyond a critical field-detector(s) coupling [14]. The detectors are assumed to have identical centre of mass rest frame internal harmonic oscillator frequency ω d0 , anharmonic coupling g, and coupling strength λ 0 between each detector's internal degree of freedom and the cavity field. Possible, direct interaction terms between the detectors are not considered in our model; as we show later below, coherent enhancements in the photon production from vacuum can occur as a consequence of the detectors coupling via the cavity field, provided the average spacing between neighbouring detectors is much smaller than the resonant cavity field mode wavelength.
With the membrane undergoing driven flexural oscillations at some frequency Ω m , the centre of mass of the ith detector follows the worldline z µ i (t) = (t, L x /2 + A cos (Ω m t + φ i ) , y i , z i ), where A is the detector's centre of mass acceleration amplitude and φ i its phase. The phase φ i accounts for the fact that the membrane's GHz flexural mode wavelength (∼ few µm) is much smaller than the extent of the detector distribution ( few mm), so that individual oscillating detectors may be out of phase with respect to each other depending on their relative separation in the y-z plane. The ith detector's proper time is denoted by τ i , and its transverse centre of mass coordinates are assumed to satisfy y i ∼ L y /2 and z i ∼ L z /2. The latter condition allows the simplifying approximation that the transversely distributed detectors couple equally to the resonant cavity mode. Note that we adopt the Minkowski metric sign convention η µν = diag(−1, 1, 1, 1).
The action (1) yields the following model Hamiltonian for the cavity-TLS detectors system in the rest frame of the cavity as derived in Supplementary Note 1: where we assume a single mode approximation for the cavity field with mode frequency ω c = ck c , k c = (2π/L x ) 2 + (π/L y ) 2 + (π/L z ) 2 , and we truncate the anharmonic oscillator detector state space to that spanned by the ground and first excited energy eigenstates with transition frequencyω d and with TLS-cavity mode coupling λ defined in Supplementary Note 1. Here, we consider the (n x , n y , n z ) = (2, 1, 1) cavity mode, which has a node at L x /2 (see Fig. 1) resulting in a first order dependence on oscillation amplitude A in the interaction term of Hamiltonian (2); if we were to instead use the (n x , n y , n z ) = (1, 1, 1) mode, the coupling would involve a cosine instead of a sine term and hence would be much smaller at second order in A. The cavity single mode approximation and detector two level truncation are justified under the condition of parametric resonance between the driven, flexural membrane frequency Ω m and the (renormalized) detector and cavity mode frequencies (see later below). Hamiltonian (2) accounts for relativistic time dilation for the accelerating TLSs through the presence of the reciprocal Lorentz factors dτ i /dt = 1 − ξ 2 sin 2 (Ω m t + φ i ), where ξ = Ω m A/c, and may be thought of as a relativistic, parametrically driven Dicke model [28,29].
With for example an achievable membrane flexural, microwave scale frequency Ω m ∼ 2π × 10 GHz and oscillation amplitude A ∼ 10 −10 m, we have ξ = Ω m A/c ∼ 10 −8 and hence the time dilation factors can be safely neglected for potential laboratory realizations. The sinusoidal interaction term can then be well-approximated as , and with ω c ∼ Ω m under conditions of resonance, we see that the individual TLS detector-cavity couplingλ is effectively reduced by the maximum TLS velocity-to-speed of light ratio v m /c 1. Nevertheless, from a theoretical standpoint it is still interesting to allow also for relativistic TLS accelerations in exploring photon production when starting initially with the TLSs in their ground states and the cavity in its vacuum state. As we show just below, even under conditions of extreme relativistic TLS motion, Hamiltonian (2) can be mapped onto a much simpler, time-independent Hamiltonian through a type of renormalized, rotating wave approximation (RWA), provided the TLS-cavity couplingλ is sufficiently small (see Supplementary Note 1 for the details of the mapping).
In particular, setting the phases φ i = 0, we can describe the collection of N TLSs as a single N + 1-level system viewed as a large pseudospin vector of length j = N/2, with the collective spin operators given by Expanding the time-dependent terms in Eq. (2) and keeping only terms up to second harmonics in Ω m , the TLS transition frequency is renormalized as ω d =ω d D 0 and there is an additive frequency modulatioñ ω d D 2 cos(2Ω m t), where D 0 and D 2 are ξ-dependent constants (defined in Supplementary Note 1). Imposing the resonance condition Ω m = ω c + ω d , we then transform the frequency renormalized, modulated Hamiltonian to the rotating frame via the unitary operator U RF (t) = exp iω c a † at + iJ z [ω d t +ω d D 2 sin(2Ω m t)/2Ω m ] . Applying the rotating wave approximation (RWA), we finally have that the relativistic, driven Dicke type Hamiltonian (2) can be replaced by the following much simpler approximate Hamiltonian: where the coupling constant is , with C 1 a ξ, Ω m -dependent constant (defined in Supplementary Note 1), and J 0 (z), J 1 (z) are Bessel functions of the first kind (not to be confused with the spin operators J ± ).
An actual microwave cavity will be lossy, while TLSs in their excited state will relax through photon and phonon emission, with the latter decay channel a consequence of the TLSs embedded within an elastic membrane. For the possible realization considered below, the dominant TLS relaxation channel is through phonon emission [23,24], and since as mentioned above the phonon wavelength is much smaller than the extent of the TLS distribution, we model the TLSs as coupled to approximately independent environments. We assume that the cavity-TLSs system dynamics can be described by the following Lindblad master equation: where γ c and γ d are the cavity and individual TLS energy damping rates, respectively, and the Lindblad superoperator is defined as Here, we suppose that the environment temperature is small compared to the frequencies of the cavity mode and TLSs (i.e., k B T / ω c , ω d ), and also that each TLS has approximately the same damping rate; for the possible realization with microwave cavity and detector frequencies of a few GHz, it suffices to work at dilution fridge temperatures of a few tens of mK or below in order to be in this low temperature regime. Analysis of the model The single cavity mode approximation in Hamiltonian (2) is justified provided the approximate resonance condition Ω m ≈ ω c + ω d is satisfied [14]. This resonance condition also justifies the RWA Hamiltonian (3), provided |λ| ω c , ω d . The advantage to working with the much simpler Hamiltonian (3) is that the Lindblad master equation (4) can be solved numerically for up to 20 or so TLSs, while for the full, relativistic time-dependent Hamiltonian (2), numerically solving the Lindblad equation is only possible for less than around 10 TLSs before the computational run times become excessively long. It is noteworthy in this respect that using the approximate Hamiltonian (3) in the master equation provides an accurate description of the average photon number dynamics even for relativistic TLS centre of mass motions with ξ = Ω m A/c 1 (see Supplementary Figure 1). Furthermore, for the assumed independent TLS damping, the average photon number dynamics is relatively insensitive to the phase φ i -dependencies of the TLSs' motions, justifying setting φ i = 0 above (see Supplementary Figure 1). From now on, we will use Hamiltonian (3) in the master equation (4).
One analyical approach to solving the above master equation (but with independent damping replaced by collective damping [29]-see below) begins with the Holstein-Primakoff (H-P) transformation which maps the collective spin operators J ± to the bosonic creation and annihilation operators b † , b [31]. In the large j = N/2 limit, the collective spin of the N TLSs is approximated as a single harmonic oscillator, with Hamiltonian (3) mapped to √ N λ(a † b † + ab). This coincides with the RWA Hamiltonian that was used in Ref. [14] to describe the oscillatory acceleration for a system comprising a photodetector coupled to a cavity electromagnetic field in the single mode approximation, with the detector's internal degrees of freedom modeled as a harmonic oscillator, and the detector's centre of mass oscillating at a frequency Ω m that matches the sum of the cavity frequency and detector internal frequency. This system describes a nondegenerate parametric amplification process, with cavity-detector photon pairs produced from the vacuum via mechanical pumping. Including cavity and detector damping and solving analytically the corresponding master equation for the above bosonic Hamiltonian, the average cavity photon number in the long time limit is . Note that the system exhibits a parametric instability when the effective coupling strength √ N λ exceeds the value √ γ c γ d /2. While such an instability indicates a breakdown of the above H-P derived approximation method, it does point to the existence of an enhanced vacuum photon production phase for N > N crit = γ c γ d /(4λ 2 ) in the original cavity-TLSs model dynamics given by equations (2)(3)(4), where there is no such instability [32]. In the following analysis, we will continue to use the definition N crit = γ c γ d /(4λ 2 ) also for the independent damping master equation (4); as we show in Supplementary Note 2, a cumulant expansion analysis of this equation establishes that N crit delineates different phases of photon production.
In Fig. 2a, we show for some example, illustrative parameters corresponding to N crit = 1, the dynamical behavior of the average scaled cavity photon number a † a /N starting from the cavity mode vacuum and with all the TLSs initially in their ground states at time t = 0. The solid line plots are obtained by numerically solving [33] the Lindblad master equation (4) with the RWA Hamiltonian (3); the utilized numerical method exploits the permutation symmetry of the density operator [34], which reduces the size of the Hilbert space required for the TLSs and increases the accessible value of N up to around 20, depending on the parameter choices. The dashed line plots are obtained by solving approximate equations for the non-vanishing first and second order moments derived from the master equation (4). We find dramatically improved accuracy as compared with usual cumulant expansion approximation methods [32,35,36], even for relatively small N , by instead setting certain fourth order instead of third order cumulants to zero (see Supplementary Figure 1). In particular, with each of the TLSs giving an identical contribution to the moment equations, we can replace σ z i , σ ± i σ z j (i = j) with σ z 1 , σ ± 2 σ z 1 respectively. Furthermore, utilizing the identity σ z 1 = 2σ + 1 σ − 1 − 1 and approximating that the fourth cumulants vanish, we obtain for example the following nonvanishing third moment approximation (see Supplementary Equation S21): Note that the latter approximation contains an additional term involving products of second moments as compared with the usually employed, third cumulant vanishing approximation [32,35,36].
The purpose for using the above cumulant expansion approximation is that the photon production from cavity vacuum can be investigated for N ≫ 1, relevant for possible realizations, while the numerical Lindblad solutions for smaller N are useful for validating the cumulant approximation, as well as giving a more complete picture of the quantum dynamics.
As N increases above N crit , a burst peak of cavity photon production from vacuum appears in Fig. 2a that progressively grows in magnitude, narrows, and shifts to earlier times. Furthermore, the long time limit steady state average photon number grows in magnitude. These features qualitatively resemble those of Dicke superradiance [21,22], although in the latter process the TLSs are initially prepared in their excited state and the cavity in the ground state. (Note that Ref. [37] considers the effect on superradiance of N = 6 uniformly accelerating, initially excited photodetectors, which is different from the collective enhancement effect for photon detection from vacuum considered here.) In order to get a better idea of the cavity mode-TLSs state, in Fig. 2b we show for N = 15 the time dependence of the logarithmic negativity measure of entanglement E N [38] between the cavity mode and TLS ensemble and also the cavity mode state Fano factor ( (a † a) 2 − a † a 2 )/ a † a in comparison with the average scaled cavity photon number. The entanglement grows and reaches a maximum roughly when the Fano factor is a maximum, while the entanglement is relatively small, although non-zero, in the long-time limit. This non-zero entanglement growth from vacuum indicates that the cavity photon production and TLS excitation are correlated. Such entanglement probes are essential in potential future experiments in order to distinguish correlated vacuum photon production and TLS excitation from possible radiative heating effects (for which the entanglement between the cavity mode and TLSs vanishes). The Fano factor gives partial information about the cavity mode state, in particular the photon number variance, and complements the information provided by the cavity mode Wigner function snapshots shown in Fig.  2c; these latter snapshots correspond to the instants when the Fano factor reaches its peak, subsequent trough, and long-time limit steady state. The photon production from a vacuum burst corresponds to the cavity mode Wigner function rapidly first spreading out and then forming a ring with width close to that of a coherent state (characterized by Fano factor = 1), and eventually settling into a thicker ring in the steady state due to environmental diffusion (with Fano factor > 1). The Fano factor and Wigner function plots help to characterize, for example, how close the cavity mode state is to an effective thermal state (i.e., having a centred, Gaussian Wigner function with Fano factor > 1); in contrast to the usual Unruh effect involving a uniformly accelerating detector, we see from Fig. 2c that the vacuum generated, cavity photon state with its ring-like Wigner function in the long-time limit is clearly non-thermal. This is not surprising given that the TLSs' motion is non-uniform, i.e., oscillatory. Possible implementation In order to get a more quantitative sense of the average cavity photon number dynamics scaling dependence on TLS number N , as well as model possible experimental set-ups, we must go to the large N crit limit where it is not feasible to numerically solve the master equation (4). As an alternative, we can apply the approximate first and second order moment equations (see Supplementary Equation S22) that should become increasingly accurate in the large N limit, provided N is not too close to N crit [35,36], and which can be straightforwardly solved numerically for arbitrarily large N since the number of coupled moment equations is fixed and small.
We will assume in part the parameters of the 3D microwave cavity-coupled nitrogen vacancy (NV) centre defect scheme of Refs. [23,24], which observed signatures of Dicke superradiance [21,22]. In particular, we consider a coupling strengthλ = 2π × 0.07 Hz, corresponding to an NV defect coupling via its magnetic moment to the magnetic field component of the considered cavity electromagnetic vacuum mode with frequency ω c = 2π × 3.2 GHz [24]. Assuming a diamond membrane flexural oscillation amplitude A = 10 −10 m, the nonrelativistic coupling strength is λ =λω c A/(2c) ≈ 1.5 × 10 −9 s −1 . For the NV defects, the dominant relaxation process is through spin-phonon interactions with rate γ d ≈ 2 × 10 −4 s −1 [24]. On the other hand, assuming a realizable, superconducting microwave cavity quality factor Q c = 10 6 [39], we have γ c = ω c /Q c ≈ 2 × 10 4 s −1 . With these numbers, we have N crit = γ c γ d / (2λ) 2 ≈ 4 × 10 17 ; the number N of microwave field-coupled defects in the experiment of Ref. [24] is in the range (0.36 − 1.5) × 10 16 , not far below this N crit value. Fig. 3a gives the first burst peak maximum and the long time limit steady state rescaled cavity photon number a † a /N dependencies on N/N crit . We observe clear evidence of a phase transition, with the slope of average photon number dependence on N changing sharply as N moves through N crit . In particular, in the long time steady state, we obtain the following approximate analytical expressions well below and above N crit , respectively (see Supplementary  Equation S24): a † a = N Ncrit γd γc+γd ≈ 2 × 10 −26 N ≪ 1 for N N crit , and a † a = N γ d /(2γ c ) ≈ 5 × 10 −9 N ≫ 1 for N N crit ; we can clearly see that the average cavity photon number generated from vacuum is negligible below and non-negligible above N crit . While these steady state average photon numbers scale linearly with N , from Fig. 3b we see that in contrast the first burst peak scales as N α with α ≈ 2 for N > N crit (i.e., quadratically with N to a good approximation); the peak average photon number rapidly grows in magnitude relative to the steady state value with increasing N above N crit . In Fig. 3c, we see that the delay t d in the appearance of this first burst peak starting from the initial cavity vacuum state becomes shorter as N increases, with the inverse delay t −1 d scaling linearly with N . While these observed scaling dependencies for the first burst peak of photon production from vacuum coincide with those for Dicke superradiant bursts [24], they are not properly a superradiant phase however [32,36,40]. In order to understand better the nature of this enhanced vacuum photon production phase, it is informative to apply the unitary transformation U = exp(iπJ x /2) to the master equation (4); the RWA Hamiltonian (3) then transforms to the Tavis-Cummings Hamiltonian H = λ(a † J − + aJ + ), which involves co-rotating terms only, while the Lindblad , with each TLS initial ground state transformed to its corresponding excited state. The cavity-oscillatory TLSs system is thus unitarily equivalent to a Tavis-Cummings model with an incoherent pump. While the latter model does not exhibit a Dicke superradiant transition, which requires the presence of both co-rotating and counter-rotating terms of the Hamiltonian, it nevertheless exhibits a so-called "inverted lasing" (or "counterlasing") transition [32,36,40]. For this reason, the transition to enhanced acceleration detection from vacuum may be thought of as a transition from a normal, incoherent phase below N crit to a coherent inverted lasing phase above N crit .

Discussion
The challenge of any experimental scheme is to exceed the threshold for coherent photon production from vacuum (equivalently the inverted lasing threshold [32,36,40]), which in terms of the TLS number N is given by the condition N > N crit = γ c γ d / (2λ) 2 . In order to reduce the size of N crit , we require microwave cavities with large quality factors, TLSs with small damping rates, and large TLS-photon coupling strengths λ. One way to enhance λ is to reduce the microwave cavity volume, for example by utilizing 2D coplanar microwave cavities [41][42][43]. However, this creates challenges for locating a sufficiently large number of TLSs within the microwave cavity region. On the other hand, while achievable λ couplings for 3D microwave cavities are a few orders of magnitude smaller than what is possible for 2D coplanar cavities, correspondingly much greater numbers of TLSs can be located within the 3D microwave cavity region, furnished by defects within a flexurally vibrating membrane. The above mentioned NV centre scheme for probing Dicke superradiance [23,24] is a promising direction, although larger cavity quality factors are required, as well as larger surface area diamond membranes that can be actuated into GHz frequency flexural motion.
While we have focused in the present investigation on a microwave cavity-GHz oscillating membrane defect scheme, it is worth considering other possible schemes as well that effectively incorporate many accelerating photodetector/emitter degrees of freedom. In particular, it would be interesting to revisit the oscillating electron-photon pair production scheme of Ref. [9], but instead consider many electrons accelerating in a strong, periodic electromagnetic field. Under conditions where the dominant wavelength of emitted photon pairs is large compared to the average electron separation, we might expect a superradiant-like coherent enhancement of the pair production rate. Alternatively, we might consider a cloud of atoms (i.e., quantum dipoles) [44], although the challenge to impart sufficient acceleration magnitude to the cloud's centre of mass would need to be addressed. Furthermore, electron-electron and dipole-dipole interactions would likely need to be taken into account; such interactions might result in quantum synchronization [44], where small initial differences in the individual electron/dipole acceleration magnitudes and phases do not influence the resulting coherent emission of radiation starting from the vacuum state. Finally, with the long-sought goal to demonstrate the original Unruh effect, it would be interesting to investigate possible collective enhancements of photodetection from vacuum for a dense cloud of detectors accelerating uniformly over a sufficiently long time interval in Minkowski vacuum, such that the detected photon spectrum is thermal to a good approximation.
In conclusion, we have proposed a way to demonstrate photon production from vacuum as a consequence of genuinely accelerating photodetectors. Our scheme involves a membrane with a dense cloud of embedded NV defects (equivalently TLS photodetectors) undergoing driven, transverse flexural vibrations within a microwave cavity, and modelled as a driven Dicke-type Hamiltonian. Under the condition of resonance, where the TLS's centre of mass acceleration frequency matches the sum of the cavity mode and internal TLS's transition frequency, the system Hamiltonian considerably simplifies via the RWA, thus allowing for an accurate analytical solution using a cumulant expansion approach. When the number N of TLS defects exceeds a critical value N crit , the system undergoes a transition to an inverted lasing phase, signalled by a 'burst' peak in the average cavity photon number that scales as N 2 , yielding significantly enhanced, collective photon production from vacuum.
While the primary motivation for the present investigation concerns the interplay between the fundamental physics of relativistic quantum field vacuum and many body quantum dynamics, the possibility to realize entangled many TLS (i.e., qubit)-microwave photon states, purely by mechanically driving a membrane structure initially in a cavity vacuum, may also find broader relevance and application in quantum information science and technology.

Data availability
The codes and generated data that were used to produce the figures in this paper are available by request from the corresponding author.   In the following, we derive the RWA Hamiltonian (3) starting from the relativistic cavity field-coupled-N detector action (1).
The starting model action for the cavity field-N detector system is given by where the ith detector's worldline is z µ i (t) = (t, L x /2+A cos (Ω m t + φ i ) , L y /2, L z /2), and τ i denotes the ith detector's proper time. With the relation S = dtL, the system Lagrangian in the laboratory frame is where the Lorentz factor dτ i /dt = 1 − ξ 2 sin 2 (Ω m t + φ i ) with ξ = Ω m A/c. Performing the Legendre transformation on the Lagrangian (S2), we obtain the following cavity field-detector system Hamiltonian: where Π(t, r) =Φ(t, r)/c 2 and P i (t) = m 0Qi (t)dt/dτ i are the momenta conjugate to Φ(t, r) and Q i (t), respectively. The first term in Eq. (S3) is the free cavity scalar field Hamiltonian, the second term is the free detector Hamiltonian, and the third term is the interaction Hamiltonian. Quantizing by replacing the position and momentum coordinates with their corresponding operators, we assume weak anharmonic potential energy terms for the detector degrees of freedom, so that the free detector energy eigenvalues E n , n = 0, 1, 2, . . . , can be approximated using first order perturbation theory. The approximate energy eigenstates for the free detectors are then harmonic oscillator Fock states: |n 0 , n 1 , . . . , n N , n i = 0, 1, 2, .... With the parametric resonance condition applied between a given detector's centre of mass oscillation frequency, cavity mode frequency, and transition frequency between the ground and excited detector's energy levels (see below), we can truncate the Hilbert space of each detector to its two lowest energy eigenstates |n i , n i = 0, 1, associated with the ground and first excited energy eigenvalues E 0 , E 1 , respectively. Therefore, definingω d = (E 1 − E 0 )/ , the free detectors' position operators and free Hamiltonian are expressed in terms of Pauli operators: We consider a 3D cavity box with side lengths L x , L y , and L z , and impose the Dirichlet boundary conditions Φ(t, 0, y, z) = Φ(t, L x , y, z) = 0, Φ(t, x, 0, z) = Φ(t, x, L y , z) = 0, and Φ(t, x, y, 0) = Φ(t, x, y, L z ) = 0. The cavity quantum field operator can be decomposed in terms of classical normal mode solutions and associated creation and annihilation operators a n , a † n as follows: sin (k nx x) sin k ny y sin (k nz z) a n (t) + a † n (t) , with k ni = n i π/L i , i = x, y, z. Under the resonance condition (see below), the cavity field is truncated to the single mode with vector label n = (2, 1, 1). The corresponding mode frequency is ω c = ck c , k c = (2π/L x ) 2 + (π/L y ) 2 + (π/L z ) 2 . Substituting (S4-S6) into Hamiltonian (S3) and dropping the subscript 2, 1, 1 on a ( †) 2,1,1 for notational convenience, we then obtain the single-mode, relativistic parametrically driven Dicke Hamiltonian (2): whereλ = √ 2 λ 0 / m 0 ω d0 cL x L y L z | k 2,1,1 |. Hamiltonian (S7) expressed in terms of J z , J ± , ξ = Ω m A/c, and with the detectors' (TLSs') phases all set to zero (φ i = 0), is as follows: where the reciprocal Lorentz factor is dτ /dt = 1 − ξ 2 sin 2 (Ω m t). Applying the Fourier series expansion to dτ /dt and the Jacobi-Anger expansion to the sin [ω c ξ cos(Ω m t)/Ω m ] term, we obtain sin where J 2n+1 (z) is a Bessel function of the first kind. Keeping only terms up to second harmonics in Ω m , Eqs. (S9) and (S10) become approximately where the ξ dependent D 0 and D 2 coefficients can be read off from Eq. (S9) and the ξ, Ω m dependent coefficient C 1 = 2J 1 (ω c ξ/Ω m ) can be read off from Eq. (S10). Substituting Eqs. (S11), (S12) into Hamiltonian Eq. (S8), we obtain where ω d =ω d D 0 is the renormalized detector oscillator frequency. Transforming to the rotating frame via the unitary operator U RF (t) = exp iω c a † at + iJ z [ω d t +ω d D 2 sin(2Ω m t)/2Ω m ] , the cavity mode and detector annihilation operators pick up time-dependent phase terms as follows: The system Hamiltonian (S13) then becomes in the interaction picture: where B =ω d D 2 /2Ω m < 1. Making use of the Jacobi-Anger expansion again such that e ±iB sin(2Ωmt) ≈ J 0 (B) ± 2iJ 1 (B) sin(2Ω m t), (S16) and substituting Eq. (S16) into Eq. (S15), we arrive at the following expression for the system Hamiltonian: Imposing the parametric resonance condition Ω m = ω c + ω d and combining the cos(Ω m t) term with the first two terms within the braces, we retain the time-independent terms and drop the oscillating terms at integer multiples of Ω m (RWA). As a result, we recover the approximate time independent Hamiltonian (3): where we have dropped the subscript I and the renormalized coupling constant λ = (S18)]. The plots are obtained by numerically solving for some example parameter values the Lindblad master equation (4) using the Quantum Toolbox in Python (QuTiP) [1]. From the figure, we can see that the Lindblad master equation (4) with RWA Hamiltonian (3) [Eq. (S18)] gives a good approximation to the average photon number dynamics for the parametrically driven Dicke Hamiltonian (2) [Eq. (S7)], even for relativistic motion with ξ ∼ 1. Furthermore, the dynamics is relatively insensitive to the individual TLS's phases φ i . Under nonrelativistic conditions ξ 1 relevant for possible experimental realizations, we have D 0 ≈ 1, D 2 ≈ 0 (with dτ /dt ≈ 1), and C 1 ≈ ω c ξ/Ω m , so that the renormalized coupling is λ ≈λω c ξ/(2Ω m ); the Lindblad master equation (4) with RWA Hamiltonian (3) [Eq. (S18)] gives an even more accurate approximation to the average photon number dynamics for the parametrically driven Dicke Hamiltonian Hamiltonian (2) [Eq. (S7)] in the nonrelativistic regime than for the relativistic regime.

SUPPLEMENTARY NOTE 2: CUMULANT EXPANSION
In the following, we derive from the Lindblad master equation (4), approximate equations for the second moments of the cavity and TLSs degrees of freedom via a cumulant expansion.
We start from the dynamics of the cavity-TLS system as described by the Lindblad master equation: where the Lindblad superoperators are defined by L A [ρ] = AρA † − 1 2 A † Aρ − 1 2 ρA † A. Consider the unitary transformation G = exp[iθ(a † a − J z + N/2)], which leaves Hamiltonian (3) [Eq. (S18)] invariant. The master equation for the transformed system density matrixρ = GρG † is given bẏ where we have applied the transformation relations G † aG = e iθ a and G † σ − i G = e −iθ σ − i . Assume that the system state is initialized as a product state of the individual ground states of the isolated cavity and N isolated TLSs: |ψ(0) = |0 c ⊗ |−, −, . . . , − s . The fact thatρ(0) = ρ(0) and the transformed equation (S20) coincides with Eq. (S19), indicates thatρ(t) = ρ(t) for the whole time range: G is a symmetry of the cavity-TLS system and initial state. As a consequence, only operators that are invariant under the G transformation give non-vanishing expectation values, and thus for the first and second moments, we need consider only the following non-zero moments: σ z 1 , a † a , aσ − 1 , a † σ + 1 , and σ + 1 σ − 2 ; all TLSs give the same moment values (a consequence of having the same assumed coupling and damping rates), which allows us to replace σ z(+,−) i and σ + i σ − j (i = j) with σ z(+,−) 1 and σ + 1 σ − 2 , respectively. We have verified numerically for a few N TLSs using QuTiP [1] that the following noninvariant moments that would otherwise appear in the moment equations vanish, i.e., a 2 = aJ + = J + J + = 0.
With Eq. (S19), the dynamical differential equations for the second moments can be obtained through a cumulant approximation [2]. These equations also include the nonvanishing third moment terms a ( †) σ ± 2 σ z 1 and a † aσ z 1 . The usual way to approximate these third moments is to set to zero the corresponding third order cumulants and obtain an approximate expression involving a sum of products of first and second order moments. However, we instead utilize an alternative cumulant approximation that is found to be more accurate. To proceed, we first rewrite σ z 1 through the identity σ z 1 = 2σ + 1 σ − 1 − 1, and substitute into the third moments to obtain the sum of a fourth moment term and one second moment term (note that alternative substitutions such as σ z 1 = 1 − 2σ − 1 σ + 1 give poorer approximations than the former, normal ordered choice). By approximating the fourth moments through setting the fourth cumulants to zero, we obtain the following improved third-moment approximations: The resulting approximate differential equations for the nonvanishing second moments and σ z 1 are d dt a † a = −γ c a † a + 2iN λ aσ − 1 ,