Vacuum radiation and frequency-mixing in linear light-matter systems

Recent progress in photonics has led to a renewed interest in time-varying media that change on timescales comparable to the optical wave oscillation time. However, these studies typically overlook the role of material dispersion that will necessarily imply a delayed temporal response or, stated alternatively, a memory effect. We investigate the influence of the medium memory on a specific effect, i.e. the excitation of quantum vacuum radiation due to the temporal modulation. We construct a framework which reduces the problem to single-particle quantum mechanics, which we then use to study the quantum vacuum radiation. We find that the memory changes the vacuum emission properties drastically: Frequencies mix, something typically associated with nonlinear processes, despite the system being completely linear. Indeed, this effect is related to the parametric resonances of the light-matter system, and to the parametric driving of the system by frequencies present locally in the drive but not in its spectrum.


Introduction
Light experiences dispersion as it passes through an optical medium, such as the glass in your window or the water in your glass, and different frequencies appear to be travelling at different rates. On a quantum level, the vacuum inside the glass is different from the vacuum outside it. In light of recent studies that return to the problem of time-dependent media [1][2][3][4][5][6][7][8], it is worth asking if dispersion plays an additional, non-trivial, role also in a medium whose properties changes with time.
Optical dispersion is of course well known, and is accurately described by the theory of macroscopic electrodynamics [9,10], in which one ignores the microscopic make-up of the medium, replacing the chain of absorption and re-emission processes of the constituents (from which the dispersion originates) with a phenomenological frequency-dependent permittivity ε. This greatly simplifies the problem on a classical level, but introduces some difficulties when attempting to quantise the theory, as the Lagrangian and Hamiltonian of the theory becomes ill-defined. Consequently, many different approaches have been pursued (a good review of which can be found in [11,12] and references therein). On a conceptual level these issues have now largely been resolved by introducing phenomenological microscopic degrees of freedom, often in the manner proposed by Hopfield [13], see for instance the work of Huttner and Barnett [14] or Philbin [15]. Such phenomenological microscopic degrees of freedom usually consists of a simplified version of medium constituent dynamics (i.e. microscopic detail), which nonetheless gives the same result at a macroscopic scale. The presence of the medium directly affects the electromagnetic vacuum, leading to Casimir-Polder forces, as is discussed thoroughly in [16]. Calculations can however become complex within these frameworks, especially when introducing timedependencies to the medium.
In this work we will study the temporal modulation of bulk media, or fibre-like scenarios, at multiple frequencies. Specifically, our aim is to develop models of quantum vacuum radiation relevant for experiments such as [62] and [63], due to recent progress in photonics as well as experimental interest. In particular, we examine the role of dispersion, especially with regards to the temporally delayed response, in the production of photons from the vacuum state due to the medium time-dependence. Interestingly, different physics transpire depending on whether it is the light-matter coupling (such as Rabi frequency) or the resonance frequency that is temporally modulated. As we will discuss, the former acts similarly to direct driving whereas the latter, on which we will focus, is a type of parametric driving. Whilst the spectrum of vacuum radiation is qualitatively similar in both cases to first order in the size of the refractive index modulation, this is no longer the case at higher orders (also discussed in [64]). We will therefore focus our attention on non-trivial second order effects, which is a topic of increasing importance with the advent of so-called epsilon-near-zero materials [62,[65][66][67][68] where changes to the refractive index in time can be in the order of unity. We will employ a model for macroscopic electromagnetism where the microscopic degrees of freedom is treated phenomenologically, in the spirit of Hopfield [13], and similar to [14,59,69,70]. This model allows to fully account for dispersion and memory effects. As a result we uncover a frequency mixing mechanism that modifies the spectrum of the emitted photon pairs.
Usually quantum vacuum radiation is emitted when the sum of two polariton frequencies match the frequencies contained within the spectrum of the modulation [51,59]. In our case, we modulate the resonance frequency at ν 1 and ν 2 , and the spectrum is thus strongly peaked around these frequencies. However, multiplefrequency modulations will form an interference pattern in the time domain, which oscillates at frequencies outside the spectrum. Question then becomes whether or not energy can be absorbed by this. Interestingly, we find that frequency-mixed photons appear when the sum of two light-matter quasiparticle frequencies match ν 1 , ν 2 or n n  | | 1 2 . The latter are indeed the beating frequencies. This not only provides a physical manifestation of time-dependent media but also provides an additional route for the detection of photons in a background free environment (i.e. at frequencies that are displaced from those of any input fields). Whilst frequency mixing is usually connected to nonlinear processes, here the underlying assumption is that the medium response is at all times linear. Instead, the mixing phenomenon is related to a parametric response of a coupled system. In particular, we find that energy can be absorbed from the modulation interference pattern precisely because of the time-delayed response of the medium. In this process, energy is absorbed from the wave oscillating at ν 1 , stored until a (anti)quanta of energy is absorbed by the second wave oscillating at (−)ν 2 (or vice-versa). The total energy n n  | | 1 2 is then emitted in the form of a polariton pair. The latter is related to the 'superoscillations' studied in [64], and the 'bichromatic' driving briefly mention in [54] can be seen as a special case of this.
The manuscript is structured as follows: In section 2, we define a microscopic phenomenological action for the light-matter system, whose classical equation of motion results in a common type of dispersion relation. We then define polariton branches and quantise using a path integral formalism in section 3. Transition amplitudes for temporally modulated media are then discussed in section 4, an in-depth example of which we treat in section 5. Discussion of the methods and concluding remarks are then presented in section 6.

The model and effective action
It is well-known that the dispersive response of the medium complicates calculations. The origin of this complexity is the two distinct types of time dynamics at interplay: optical parameters that change with time, as well as the time-delayed response of the medium. The time-delayed response is directly connected to dispersion, as the rate at which the medium constituents absorb and re-emit light depends on the frequency. Such frequency-dependence of the response implies, by necessity, that the Hamiltonian/Lagrangian is nonlocal in time. The medium is therefore characterised, in the time-domain, by a memory kernel connecting past events with the present [9, 10, 71]. In the context of macroscopic electromagnetism, a time-dependent medium is introduced by allowing a model parameter, such as the resonance frequency, to change with time. The resulting time-dependent permittivity is then described by a memory kernel which changes non-trivially with time.
In this work, we will model the optical medium as a set of harmonic oscillators R i with natural oscillation frequencies Ω i respectively, at a spatial density of ρ. As we will see shortly, these oscillation frequencies will act as the resonance frequencies of the medium. Note, we will use units such that c=ÿ=ò 0 =1 for notational simplicity. Coupling this to electromagnetism by dipole terms, quantified by dipolar coupling strengths q i , yields the action where A and j are the vector and scalar potentials respectively. In Coulomb gauge, i.e. when  = · A 0, the equation of motion for the scalar potential is and substituting in j from equation (2) leads to Since the last term is just quadratic in the oscillator fields R i , we can re-diagonalise. Furthermore, this can be done without impacting the form of the action, since the oscillator parameters ρ, Ω i and q i are all phenomenological, i.e. chosen to fit experimental data. We therefore arrive at an action describing the electromagnetic vector potential A coupled to a set of oscillators R i by a dipole term, where the latter phenomenologically take into account the microscopic details of the matter degree of freedom, given by where R i is the position of each oscillator in its potential well and ρ is the density of oscillators. This action is inspired by the Hopfield models employed in [13,14,59,69,70]. In the case of constant W ºW ( ) t x, i i , we find that equation (3) leads to a dispersion relation for the electric field of the familiar Sellmeier form, as is widely adopted in the optics literature [72]. In other words, the above action is a suitable starting point for modelling any dielectric where absorption is negligible. From this we can see that a time-dependent Ω i induces temporal changes in the refractive index. We note that it is however also possible to create a time-dependent medium through a coupling strength q i that depends on time. This has been studied in various scenarios, referred to as a time-dependent Rabi frequency, and we will not delve deeply into this scenario here. Already at this stage we can see, from the above action (equation (3)), that such a time-dependence will act more akin to a direct driving force than a parametric drive. Similarly to [69], we want to compute an effective action for the photons, [ ] S A eff , by integrating out the oscillator degree of freedom. Schematically, we do this by computing , as we are not interested in the dynamics of R i . In this path integral we integrate over each possible configuration of the oscillator position R i as a function of time that fulfils the stated boundary conditions, as defined in [73]. However, as the coupling in is linear, it is easy to show that the quantum fluctuation of R i does not affect A and is contained in the normalisation constant  (here set to unity) [73]. Therefore, performing this path integral for R i with the above boundary conditions is equivalent to solving the classical equation of motions for R i driven byq A i [74]. This can be done by method of Green's functions, that is, by solving , , 0 . We can link this to usual optics parameters by noting that the medium response function, commonly denoted as c ¢ ( ) t t x, , [71], is given by , , In other words, substituting yields the effective action for photons. In order to make this more tangible, let us also expand the vector potential in the polarisation vectors = å l , defined with respect to some reference vector p such that e λ ·p=0. We should note here that since A is completely transverse, so are the oscillators R i , and thus satisfy the Coulomb 'gauge' condition  = · R 0 i . Through this, we find the effective action where Δ i is the oscillator propagator given in equation (5), with r = g q i i 2 2 being the effective plasma frequencies for each resonance respectively.
Since the two polarisations de-couple, we will from here on drop the λ subscript for notational simplicity, and work only with the scalar quantity ( ) A t x, . This is so far general, and we have specified neither the space nor the time-dependence of W ( ) t x, i 2 . In the next section, we will consider the case of a static but inhomogeneous set of oscillators, such that W º W ( ) . The spatial dependence will be taken into account by expanding in an appropriate set of normal modes ( ) u x k , the exact form of which depends on the physical situation. Let us start with the equation of motion for the vector potential from equation (6), under the assumption that the oscillator frequency is time-independent. This is given by The goal is now to expand the vector potential in a set of normal modes . The form of the mode functions ( ) u x k depends on the physical scenario. We thus look for functions that satisfy with appropriate boundary conditions for the situation. In this work, we will focus on bulk media, and as such there is no spatial dependence on the oscillator frequency (W º W ( ) x i i ). A brief aside into a fibre-like scenario can be found in appendix A. There are nonetheless multiple ways of expanding in terms of normal modes for bulk media.

Plane waves
For bulk media, a natural choice of normal modes are the momentum modes Here we find the dispersion relation given by equation (4).

Paraxial waves
In most experimental scenarios however, the simple plane waves are not accessible, and are instead replaced by structured paraxial beams. Let us once again consider a homogeneous bulk medium where W º W ( ) x i i 2 2 , but where we restrict equation (7) to the paraxial limit, with the z-direction chosen to be the propagation direction.
In other words, let with k being the momentum in the z-direction and  q k 2 1 2 2 , where r is the transverse plane coordinates and q its associated momentum. Similarly to [75], we then find that k must follow the dispersion relation of equation (4), and r ( ) u k satisfies the paraxial wave equation where 2 is the transverse Laplacian. Solutions include the familiar Laguerre-Gaussian modes and Hermite-Gaussian modes [76].

Quantisation
It is often noted that constructing a quantum field theory reduces to quantising an infinite set of harmonic oscillators [73,77,78], one for each (commonly continuous) position/momentum. By expanding a field in terms of suitable normal modes however, one can take this seriously and proceed to quantise each normal mode. This is usually done within the canonical quantisation scheme, but we will here use a path integral language. Whilst this quantisation technique is known for plane waves in vacuum [73,79], it is not commonly employed for computations, nor has it been generalised for dispersion. As we shall show here however, we find this technique particularly suitable for tackling the type of problems addressed by macroscopic quantum electrodynamics.
Let us start by re-writing the effective action of equation (6) in the frequency domain, yielding , is a real quantity. Here we integrated by parts on the -term and used that D ¢ ) D k, depends on the particular normal modes used (see appendix A). As an example, for plane waves this reduces to equation (4).
The solution of the classical equations of motion for each normal mode takes the form of , . This defines the quasiparticles of the system. In other words, by 0for ω as a function of normal mode label k, we find N quasiparticle branches. These are usually referred to as polaritons. The exact number of polariton branches depends on the explicit form of w ( ) D k, . We will label these branches by the subscript α, and an example can be seen in figure 1(a). Inspired by [69], we can do the following field transformation in order to define a polariton action. Note, this transformation is always well-defined as the w , diverges at the same points and at the same rate as w w a ( ) k 2 2 goes to zero. The action of equation (8) is nonlocal in time (i.e. dispersive) in field-coordinates w ( ) A k : by this field-transformation, we trade nonlocality in time for nonlocality in space. This simplifies the quantisation procedure whilst still taking dispersion into account, as dispersion is now implicit in the definition of the polariton fields and their dependence on the momentum mode k. Temporal nonlocalities in quantum theory can be treated, but usually at a computational cost since one must now define a (commonly) infinite set of conjugate momenta (see for instance the discussions in [80,81]). A spatial nonlocality on the other hand, which here means that the polariton frequencies w a ( ) k contain terms of higher order than k 2 , is straightforward to tackle since we will treat each momentum mode k independently.
Written in field-coordinates A α , and after transforming back into temporal space, the action is that of a set of complex harmonic oscillator This is the action with which we will work. From now on, we will be working with the dynamics of single normal modes k, so we will drop the sum over k and corresponding identifier in order to simplify notation. Also, we will drop the index α on all but the mode frequency w a for the same reason. Although this is a field theory, in terms of normal modes, all the usual techniques from single-particle quantum mechanics apply. This can be quantised in the manner most familiar to the reader. In this work, we choose a path integral method as it allows for a straightforward definition of time-nonlocal perturbation theory.
As is usual in path integral quantisation, we want to add the driving terms * JA and * J A to the action for future use. These driving terms physically originate from free currents in the system, i.e. the movements of free charges. We will however not consider physical driving here, but use the driving terms for computational purposes. After simplifying the notation and adding the driving, the action takes the form Let us now proceed by calculating the polariton transition amplitude where we have the boundary condition A(t i )=A i and A(t f )=A f . Here we calculate the probability amplitude for a polariton in branch α, normal mode k and polarisation λ, starting with field amplitude A i at time t i and transitioning to field amplitude A f at time t f 3 . First we note that the quantum fluctuations decouple from the classical dynamics, as the action in equation (11) is quadratic in the fields. As a consequence, the transition amplitude factorises as ( cl is the classical action and the pre-factor ( ) T is determined by the quantum fluctuations η. We here define the quantum fluctuation by splitting the field into classical and quantum components h Explicitly, this pre-factor is given by We calculate the classical action using the equation of motion along with the boundary conditions at t i and t f . Finally, we find the transition amplitude is that of a complex driven simple harmonic oscillator. See appendix B for a detailed calculation. As usual, this expression contains all information required for computations.

Connecting polaritons and photons
The field transformation in equation (9) also has a physical interpretation. In doing this, we project the photon field in terms of polariton fields. The expansion coefficients, a generalisation of the co-called Hopfield coefficients [13], are given by  w a ( ) k . As we are always quadratic in the fields, it is convenient to define the squared coefficients   w º a a ( ) k k 2 . These are given by Here we have used the fact that the polaritons live on-shell (i.e. ω = ω α ), and that 2 . It can easily be shown that    a 0 1 k , and corresponds physically to a factor describing the degree to which the polariton is 'photon-like'. In other words, in spectral regions where w a  k, this factor is close to unity, and vice versa. In figure 1(b), an example of this can be seen.
We should note that in order to go from polariton observables to photon observables, the field transformation in equation (9) needs to be undone. In general, integral expressions will come with factors of  a k when transforming from polariton to photon degrees of freedom, although in the actual path integral it can be absorbed into the normalisation.

Transition amplitudes
A time-dependent medium can generally change the number of polaritons in the system: Quanta can be excited from the vacuum [44], whose accompanied spectrum is of interest, and like-wise polaritons can be absorbed into the vacuum. The former is the vacuum radiation. Each process has a transition amplitude a ¬ ( mn pq f i k , denoting a transition from a (pq)-state with p+q polaritons at time t i into a (mn)-state with m+n polaritons at time t f , whose absolute square gives the associated probability. Here we will first consider this general situation. We will once again drop the k identifier to simplify notation, unless otherwise stated. Throughout this, we will use a quantisation box of volume  , as is standard (see [77]), and the normal modes used take the form Also, we should note that these transition amplitudes are the polariton Fock space propagators. However, we will first take a detour into a system where driving is present, as this links directly to a time-dependent medium in a perturbative setting.

Generating functionals
Let us first consider a driven medium, whose amplitudes will act later as generating functionals when considering time-dependent media perturbatively. We will first calculate vacuum persistence amplitude for notational simplicity. This is given by the Gaussian integral , , , where * = d A dAdA 2 and Ψ 00 (A) is the time-independent version of the groundstate wavefunction seen in appendix C. Note that in this A is a complex variable, and not a function. Computing this yields cos .
As can be expected, this is simply the generalisation of the vacuum persistence amplitude in [73] to the case of a complex harmonic oscillator.
However, this calculation becomes increasingly complex for higher energy states, and we will therefore use a trick similar to what is done in appendix C in order to derive the Fock wavefunctionals Ψ mn . That is, we use the wavefunctionals and calculate the transition amplitude This amplitude can be seen as a generating functional of sort. If we expand f a and f b in terms of the Fock wavefunctionals Ψ mn , we find that In this way, we find the transition amplitudes The explicit form of F(b, a) J can be found in appendix D. This captures all processes possible. We can be a bit more explicit and ask ourselves what is the amplitude of exciting two polaritons back-to-back from the vacuum into mode k in branch α: where we have ignored global phases.

Time-dependent medium
Let us now turn our attention to time-dependent media. In particular, let us consider a homogeneous medium with a weak space-and time-dependent resonance frequencies, i.e. W ºW where  | | f 1. We can then perturbatively construct the oscillator propagators (equation (5)) in orders of | | f :  The 0th-order is simply the usual propagator, leading to the dispersion relation in equation (4). We can also relate the shift in oscillator frequency to the change in refractive index through where ò i is the characteristic amplitude of f i . As we will see below, we can perform the same field transformation as before (equation (9)) and we arrive at the polariton action in equation (10). The higher order propagators translate into (perturbative) potentials for the polaritons. We cannot however trade the temporal nonlocality of the higher order potentials for additional spatial nonlocality, because these terms are not diagonal in frequency space (i.e. w w ¢ ¹ -). Therefore these will be temporally nonlocal two-time potentials, also in the polariton field coordinates. For clarity, let us proceed step-by-step.
Before transforming to the polariton fields we have the effective action where we have let w w ¢  -¢ in the second integral, and we have defined the auxiliary propagator Applying the polariton transformation in equation (9) to the above action, and simplifying the notation, yields where we sum over momenta ¢ k and polariton branches a¢. Also, we have now defined the polariton projected i t i t kk k k k k with  w a ( ) k being the polariton projection operator defined in equation (9). This takes the form of a complex harmonic oscillator, along with an additional two-time harmonic potential , , 2 0 k k k k k , which connects the normal mode at k with the one at ¢ k . This latter term we will treat perturbatively, which is done by computing There are two separate sectors here: either polaritons are excited from the vacuum into the same polariton branch, or two separate ones. We will treat these two sectors separately for clarity, and will be referred to as intrabranch and interbranch vacuum radiation respectively. In both cases, we are interested in the probability amplitude of exciting a polariton pair back-to-back, as illustrated in figure 2.

Intrabranch vacuum radiation
Let us first consider the case when a a ¢ = in the perturbative potential of equation (20). We can then compute the necessary functional derivatives to equation (14). where we have expanded to second order for consistency, and considered scattering states where = -¥ t i and = ¥ t f . We have also ignored the overall phases. Substituting the auxiliary propagator s aa kk in terms of oscillator propagators, equations (19) and (17), yields the final result:  where  is the volume of the medium. It is worth noting that in this process, where we consider two polaritons are emitted back-to-back, the medium modulation f i doesn't contribute with any additional momentum. Thus it is the homogeneous part of the modulation that is sampled. This is expected, as a pair of back-to-back polaritons automatically conserve momentum. Secondly, we are mostly interested in a periodically modulated medium, i.e. the dynamical Casimir effect, and therefore the zero frequency response is very small. Hence we can safely ignore the second line, which is proportional to˜( ) f 0, 0 i . As for the vacuum radiation spectrum, we have two separate mechanisms here. One is a direct emission that depends only on the spectrum of the modulationf , this is the first term, whereas the second term explicitly depends on past events due to the integral over axillary frequency w¢. This latter term allows for vacuum radiation resonances outside the spectrum of the modulation.

Interbranch vacuum radiation
For interbranch vacuum radiation we will first consider a slightly different driven amplitude, since in this case the two polaritons are distinguishable (at separate frequencies). Instead of equation (14), we must take the product of exciting one polariton into each of the branches. Thus we have where we have added α or a¢ identifiers for clarity, and made sure that the process conserves momentum by involving a k and a -k polariton respectively. We can now substitute equation (24) into the perturbative procedure in equation (21) where we have already neglected terms that would involve a factor of˜( ) f 0, 0 i , for the same reason as for the intrabranch polaritons. Finally, we find the probability amplitude  We will return to this amplitude shortly, but it is once again worth noting that these interbranch processes open up the possibility for a variety of frequency mixing processes, as generally ω α and w a¢ are at different frequencies.
The spectrum of vacuum radiation depends directly on the spectrum of the modulation f i , but due to the integral over w¢ in the second line, also frequencies outside is possible.

Correlators
As a quick aside, it is worth mentioning that correlators can be calculated with relative ease. This is done by applying the appropriate number of additional functional derivatives with respect to J to the transition amplitude, before setting J=0. For instance, we can calculate the field-field correlator related to transitioning from vacuum to two back-to-back polaritons by where Yñ | is the ground state ñ |0 propagated with the time-modulated kernel, and S 1 is the action seen in the exponential when calculating the perturbative transition amplitudes in equation (21). We should note that it would here make sense not to consider a transition from the vacuum at = -¥ t i to an excited state = ¥ t f , but rather from = - , and track the evolution of correlations as t = -T t increases. However for the sake of brevity, we will not further discuss correlators in this work.

Frequency mixing of vacuum radiation
In this section, we will explore further the dispersion-induced mixing processes mentioned briefly in the end of the last section. Specifically, let us consider a two-frequency time-dependence with the Fourier transform 2 .
i i In both the intrabranch (equation (23)) and interbranch (equation (27)) sectors, there is an integral over an auxiliary frequency w¢. We can evaluate this mixing integral in the large τ-limit (i.e. modulating for many periods), yielding For simplicity, we will here modulate the mth-resonance of the medium ( d = i i m ) for some large time t w a  1 . Now, recall equation (16), where we can relate the size of the modulation  i to the change in the refractive index δn. It is however more convenient to work with changes to the permittivity ε=n 2 , which we will denote δε: the two are related through de w w d w = -( ) ( ) ( ) n n 2 , and thus For this type of modulation, we see that the intrabranch amplitude can be re-written as It is worth noting that the interbranch resonances are suppressed in general, as they require both branches to be photon-like simultaneously (so that   a a¢  0 k k ). Consequently, the last two term of equation (32) can safely be neglected, as they furthermore contribute at the next order in perturbation theory. These additional vacuum radiation resonances are nonetheless possible. We can now calculate the total excitation probability by + ¬ ¬ | | G G 11 00 intra 11 00 inter 2 4 . Let us at this point specify the medium as fused silica (as in figure 1), and as we are usually interested in optical frequencies, we will modulate the first ultraviolet (Ω 2 ) medium resonance only. Specifically, we let ν 1 =Ω 2 /5 and ν 2 =Ω 2 /6, and choose ò such that δn ; 10 −3 (small but standard for fused silica). The associated probability spectrum can be seen in figure 3(a), where solid and chequered shading denotes an intraand interbranch processes respectively. The polariton branches of interest are shown in figure 3(b), along with the relevant modulation terms.
As can be seen, the temporal modulation provides the energy to resonantly connect a polariton branch with some antipolariton branch, which causes polaritons to be emitted from the vacuum state. Only the ω 1 and ω 2 branches are at a comparable scale to the modulation frequency (∝Ω 2 ), and are thus the only ones into which vacuum radiation is emitted. There are nonetheless several different possibilities, where the modulation energy will match either 2ω 2 (intrabranch) or ω 2 +ω 1 (interbranch). This opens up for the possibility of frequencymixed vacuum radiation, where the frequency of emitted vacuum radiation is given by a combination of the frequencies present in the system. Note however, that both polaritons in any given pair will however oscillate at the same frequency when measured outside the optical medium, as they are at the same wavelength, leading to a measured spectrum such as the one seen in figure 3.
Starting with the intrabranch resonances, we see the two expected dynamical Casimir-like resonances, that is, ω 2 =ν 1,2 /2. However, we also see resonances at ω 2 =ν 1,2 as well as at the mixed frequency ω 2 =(ν 1 +ν 2 )/2. The difference frequency would become relevant when w n n = - , which is in the far The polariton branch of interest as a function of vacuum wavelength. An excitation process always involves a polariton-antipolariton pair, the latter having negative frequency. The time-modulation then provides the energy connecting the two branches (denoted by coloured arrows). (c) Illustration of the possible mixing processes at second order in perturbation theory. 4 Note that the total probability density for emission is given by Also, there is no need for renormalising this integral, as we are considering differences between the occupation in each state, not the total occupation number in each.
infrared and is ignored here. Furthermore, the interbranch resonances contribute also, when ω 2 +ω 1 =ν 1 and ω 2 +ω 1 =ν 2 , denoted as chequered shading with yellow and red solid line, respectively, in figure 3. This is reminiscent of nonlinear processes, where sum and difference frequency generation is commonplace [71]. However, the system studied in this work is by assumption linear. In fact, these resonances in the spectrum of emitted vacuum radiation has much in common with the resonances of classical parametric oscillators. It is known that a stand-alone parametric oscillator with oscillator frequency Ω has a primary resonance at n W = 2 if driven at frequency ν, and several sub-harmonic resonances at Ω=ν and 3ν/2, and so on, where the strength of each resonance down the line is significantly weaker than the last [17]. Also, coupled parametric oscillators has been shown to exhibit a variety of combination (i.e. frequency-mixing) resonances [18,19,85], closely connected to the interbranch processes discussed here.
However, this does not explain resonances of the form w w n n The type of frequency mixing is of a different nature than the 'combination' parametric resonances. Instead, the mixing relates to the parametric driving of the system by the beating pattern formed by the different components of the modulation (in the time domain), which oscillates at frequencies outside its spectrum. An example of this is the 'superoscillations' studied in [64], but is in this case of much familiar origin: the two waves, n t cos 1 and n t cos 2 , beat at (ν 1 +ν 2 )/2 and (ν 1 −ν 2 )/2. The system is however unable to absorb the energy represented by this beating pattern directly. Rather it is a two-stage (virtual) process where one quanta of energy is absorbed by the first modulation wave ( n t cos 1 ), which is stored, while the second modulation wave ( n t cos 2 ) either adds or removes another quanta of energy from the system. Note, the removing of energy comes from the absorption of an anti-quanta of the modulation wave. The total energy is then emitted in the form of two polaritons.
It is worth pointing out that had we instead chosen to temporally modulate the light-matter coupling strengths g i instead of the oscillator frequencies Ω i , it is easy to see that we would not get modifications to the oscillator propagator seen in equation (15), the last line of which is responsible for the time-nonlocal integral in equation (28). Instead this would act similarly to a driving force. We would indeed also find quantum vacuum radiation in this scenario, and to  ( ) look very similar, fulfilling the condition w w n + = a a ¢ r 1,2 for some integer r [50,51]. However, we expect the contribution from both frequencies of the drive ( n n  | | 1 2 ) to disappear in this case (at least to the same order in ò).

Discussion and conclusion
In conclusion, we have studied quantum vacuum radiation excited by temporal changes to the resonance frequency of an optical medium. In particular, we have examined how the dispersive response affects the spectrum of emitted photons. We studied this with bulk media in mind, and specified fused silica as an example. We found that the delayed temporal response of the medium, responsible for dispersion, introduces frequency mixing to the system. The spectrum of emitted photons then takes on a character reminiscent of nonlinear optics, where both sum and difference frequency emission is possible.
Specifically, we showed that photons are emitted when the sum of two polariton branch frequencies match a combination of modulation frequencies. This we found led to several quantum vacuum radiation resonances, including n n  | | 2 1 2 as well as the usual dynamical Casimir-like emission at ν 1 /2 and ν 2 /2, when modulating the medium at frequencies ν 1 and ν 2 . We note that the system is by assumption linear, as to not confuse this with a nonlinear phenomenon. We found instead that there are two separate, linear, mechanisms by which frequencies can mix, related either to the energy emission process or the energy absorption process, or a combination thereof.
The mixing of polariton branch frequencies is a consequence of the nature of coupled systems having multiple modes of oscillation, which in this case are the polariton branches. In the most simple case when the optical medium only has a single resonance frequency Ω, the two modes oscillate at frequencies 4 . It follows that any excitation in the system, and hence emitted vacuum radiation, must consist of some combination of ω + -polaritons and ω − -polaritons. In the case of fused silica, there are further branches, whose algebraic form is considerably more complicated, but the physics is the same. On the other hand, the mixing of drive frequencies (i.e. the sum/difference frequency peaks) has a more subtle origin, and is connected to the time-delayed response of the medium to changes in its resonance frequencies. We find that when modulated at multiple frequencies, say ν 1 and ν 2 simultaneously, the medium can absorb energy from the beating pattern formed between the two waves. This process relies on a time-delayed response to changes in the resonance frequency, as the medium must first absorb one quanta of energy from one drive (say ν 1 ), and at a later time absorb a (anti)quanta of energy (−)ν 2 from the second drive. The total energy of n n  | | 1 2 is then emitted in form of polaritons, and as such energy conservation requires that w w n n + =  a a ¢ 1 2 . In order to study this, we used a microscopic phenomenological model for electromagnetism in an optical medium with a generic Sellmeier dispersion relation. We quantised this using a path integral formalism. No approximations were made with regards to the delayed response, and dispersion was therefore fully taken into account. Within this framework, we induced a time-dependent change to the refractive index n by weakly perturbing the resonance frequencies of medium. The model is however extendible to include also temporal changes to other parameters of the optical medium, such as the density and dipolar coupling strengths. It is worth noting that this model relates most readily to experiments in bulk media, such as in [62,63,68], rather than the typical cavity set-up where polariton physics is more commonly discussed [50][51][52][53][54].
The origin of the time-dependent resonance frequencies has not been mentioned explicitly in this work, but has been kept general. Nonetheless, the results are directly applicable to experiments in which the temporal changes to the resonance frequency originates from the quadratic Stark shift (as discussed in [70]), i.e. g W  W + E i i pump 2 , for some strong electric field E pump . Whilst this mechanism does introduce an actual nonlinearity to the system, we want to highlight that this nonlinearity affects the pump beam only, and the physics of the quantum vacuum discussed here is at all times linear, especially since typical vacuum electric field fluctuations are exceedingly weak. Indeed, this is the same line of reasoning as some recent discussion of the overlap between nonlinear optics and Casimir-Polder physics [86]. Therefore, the framework is applicable to experiments with strong electric fields propagating in bulk, or structured, media, such as the fibre experiment in [63]. In the context of these bulk media with a strong pump pulse, we expect the mixed-frequency quantum vacuum radiation discussed in this work to be readily observable. Whilst the mixing is indeed a second order effect, the fact that it allows us to shift the frequency of the vacuum radiation to ranges with better detector efficiencies, such as the optical to infra-red regime [87], greatly improves the observability of quantum vacuum radiation. Suppose that the fused silica slab in figure 3 is a thin film of roughly m 100 m thickness, and considerably larger than the pump laser spot size A spot in the transverse direction. We can estimate the number of photon pairs emitted per unit angle dθ as . This radiation would be emitted in the orthogonal direction to the pump beam (i.e. the transverse plane). This is the emission per pulse, so a repetition rate of 1 MHz yields roughly three photon pairs per second, which is measurable with current technology [87], given that they can be out-coupled from the medium (an experimental challenge but not impossible). Importantly, this frequency mixing is off-set from any other frequency of the system, and is therefore unlikely to be filtered away (a common problem for quantum vacuum radiation).
In addition, the dispersion also allows you to choose to work at frequencies where the physics is sensitive to small changes to the optical parameters, such as close to the point where the group velocity dispersion is close to zero (a common point of interest in fibre optics [88]). In fact, we would argue that this is indeed the mechanics of photon pair production in [63], albeit this requires further analysis that is outside the scope of this work.
Another experiment that relies on this mechanism is described in [7], where the refractive index of a thinfilm epsilon-near-zero metamaterial is changed rapidly in time, building on experiments performed in [62,68]. In light of the present results, additional physics can be expected associated with the linear frequency mixing mechanisms. This work suggests that the probability of emission for mixed frequency vacuum radiation to be d µ( ) n 2 , where δn is the absolute change of the refractive index. A back of the envelope calculation for the conditions in epsilon-near-zero materials (where δn ; 0.9) suggests near unity probability of emitting quantum vacuum radiation where ∼20 % of the photons emitted would be frequency mixed. Further study is required however, since this is clearly not a perturbative change to the refractive index, andabsorption cannot always be neglected. This present work does nonetheless indicate that rich physics can be explored in the spectrum of emitted quantum vacuum radiation, especially in experiments with large changes to the refractive index.
Finally, we note that we expect this vacuum radiation mixing phenomenon to be rather general, occurring in any temporally modulated system that has delayed temporal responses, and we note also that it is related to the parametric resonances of the system.