Temporal boundaries in electromagnetic materials

Temporally modulated optical media are important in both abstract and applied applications, such as spacetime transformation optics, relativistic laser-plasma interactions, and dynamic metamaterials. Here we investigate the behaviour of temporal boundaries, and show that traditional approaches that assume constant dielectric properties, with loss incorporated as an imaginary part, necessarily lead to unphysical solutions. Further, although physically reasonable predictions can be recovered with a narrowband approximation, we show that appropriate models should use materials with a temporal response and dispersive behaviour.


I. INTRODUCTION
Mobile technologies are considered the biggest technology platform in history, with transformative advances occurring across all society. These are having a profound impact in diverse areas including health care, education, and industry [1]. Key to developing mobile technologies is the ability to predict electromagnetic (EM) wave propagation though systems with differing permittivities, and in particular predicting losses. In both physics [2] and engineering [3], EM loss is often expressed using the electric loss tangent ("tan δ"), the ratio of the imaginary to the real part of the permittivity. This constant permittivity model is widely used in condensed matter physics [2], and is critical to the design of technologies diverse as mobile phones, imaging systems, consumer electronics, radar, accelerators, sensors, and even microwave therapy [4]. In such situations, microscopic models involving e.g. atomic structure or quantum mechanical effects typically provide no significant advantage.
Here we prove that modelling loss using a constant complex permittivity and permeability [15,[26][27][28] is physically incompatible with a temporal boundary. Such models can lead to unphysical post-boundary solutions that grow exponentially, despite being applied to passive and lossy materials; or fields may become complex-valued despite being real-valued before. This important point, which provides an unambiguous warning to non-specialists, is not addressed in other recent work on temporal boundaries, which focus on either reflection and refraction inside a medium with parabolic dispersion [29], appropriately generalised Kramers Kronig relations [30], or the complicated effects resulting from a Lorentzian response model [31]. Our dynamic material model, in contrast to more complicated ones [30,31] is explicitly designed to provide a minimal example with simple behaviour which clearly reveals the basic physical principles relating to the treatment of temporal boundaries: the necessity of considering the dynamics of the bound current, the requirement for material-property boundary conditions, and the resulting secondary implications for frequency-domain properties such as the dispersion relations.
In this article the term "boundary" refers to an interface between two regions with different constitutive relations (CRs). In electromagnetism CRs are most simply given by a permit-tivity and permeability µ, although more complicated CRs are also allowed. In contrast to spatial boundary conditions that describe an interface between two static media, here we consider a temporal boundary at t b , where the medium has one set of CR properties before the transition (t < t b ), and different CRs afterwards (t > t b ). In our idealised transition, the CR for the entire region changes instantly; although gradual transitions are also possible [32,33]. The unphysical consequences arise irrespective of which of the two possible types of temporal boundary condition (TBC) we consider. The first "natural" TBC [34] is derived from Maxwell's equations on the assumption that all the currents are finite, which leads to the continuity of D and B. The second "fundamental" TBC treats only E and B as physical EM fields, with D and H as derived fields acting as a gauge for the current [35][36][37]; here E and B are also continuous, but a temporal-surface dipole current appears at the transition, as shown in Fig. 1.
After demonstrating and defining the problem, we present two methods for obtaining physically meaningful solutions. The first applies the constant complex CR model whilst using a narrow band approximation (NBA). This leads to a solution based on complex conjugate pairs of frequencies and refractive indices. Although partially successful, it merely hides the fact that to model lossy media correctly when there is a temporal boundary, one requires a time-dependent material response, and therefore dispersive CR, where and/or µ depend on frequency. Besides, just as a dynamic medium model requires its own fields to represent its state, in the frequency domain we see that a dispersive medium generates one or more additional modes; and this necessary information cannot be included in the standard constant complex CR model. Thus in either time or frequency, additional boundary conditions (ABCs) must be specified.
Note that proofs and additional discussion are presented in the Supplementary Material following the References.

II. LINEAR MEDIA
In linear media any EM field can be represented in the Fourier domain by a sum or integral of terms of the form exp(−iωt + ik · x), where ω ∈ C is a complex frequency, and k ∈ C 3 a complex wavevector. Since the source free Maxwell's equations (5) are linear, with J total = 0 and ρ = 0, it is sufficient to consider just a single mode with k oriented along the z-axis (alongẑ), and E alongx. Now consider a temporally dispersive medium where the CRs specify permittivity (−ω) and permeability µ 0 , and where the "−ω" is a consequence of choosing e −iωt in (1).
If ω and (−ω) are both real then (−ω) can be replaced with (ω); but this is not allowed in our following calculations, because extra care must be taken when using complex permittivity to model damping. From Maxwell (5) we obtain a dispersion relation, and define the refractive index [38][39][40]. These are where c 0 = ( 0 µ 0 ) −1/2 is the vacuum speed of light. We now ask whether these CRs correspond to a passive lossy medium, i.e. one dampened with no external energy added. Given that both ω and k can be either real or complex, there are two possibilities that are straightforward to consider. These fit into the temporally propagated and spatially propagated viewpoints respectively [41], and are: First, if ω is real and positive, we require that plane waves are spatially evanescent in the propagation direction, which implies Im (−ω) > 0.
Second, if k is real, then we need Im(ω) < 0 to damp the field, leading to the requirement that when ω 2 (−ω) is real and positive, then Im(ω) < 0.
How the fields represented by these modes, change as they cross a temporal boundary will depend on how the change in CRs is specified, and on the chosen TBCs.

III. TEMPORAL BOUNDARY CONDITIONS
Different types of electromagnetic TBC can be identified, depending on the material response [11,12,42]; here we summarise focussing on how bound currents represent the material response. First, as depicted in Fig. 1, we could identify the sudden change in the CRs at t = t b as leading to an instantaneous temporal-surface dipole current (J surf ). The total current, J total , in the medium is where J reg is the usual finite current in Maxwell's equations: Since the fields E, B, D, H may be discontinuous we write where D a and D c are the D before and after the transition and θ is the Heaviside function. Using (4) and (6) in the Maxwell-Ampère equation, we have This approach can also be used for B in the Maxwell-Faraday equation to derive a similar result. The TBC are where , etc. One option [5,15] is to set J surf = 0, to obtain the natural TBC, i.e.
[D] = 0 and Alternatively, if treating E and B as the only physical fields, we get the fundamental TBC, i.e.
[E] = 0 and which rely on (7) to calculate the conserved temporal-surface current J surf . This is an analogous approach to that describing the surface current around a permanent magnet [35].

IV. CONSTANT COMPLEX CR GIVES UNPHYSICAL RESULTS
Even a time boundary between a vacuum with = 0 and a lossy medium with = c , where c is a non real constant with Im( c ) < 0 and Re( c ) > 0, results in a failure. For simplicity, we set Pre-boundary (t < 0), we start with a single real mode Here ω a , k a are both real and positive and satisfy the vacuum dispersion relation Since the post-boundary lossy material, with ω c and k c , must have Im( c ) > 0, the t > 0 general solution is with the dispersion relation The choice of root for n c is unimportant, as both roots are included in (11). The first root Re(n c ) > 0, since Im( c ) > 0, requires Im(n c ) > 0; and the second root Re(n c ) < 0 with Im(n c ) < 0. Applying the fundamental TBC (9), just before the time boundary, we have and B y (0 − , z) = (k a /ω a )E x . Just after the time boundary Note that except for different constants, the result is the same as would be obtained using the natural TBC (8).
From (13) and (14) we see that k c = k a or k c = −k a must hold and ω c = n a ω a /n c . Because we can choose k c = k a without affecting the analysis, we can now rearrange the four unknowns to get Now we have k c = k a > 0 and c I < 0, so E increases exponentially with time despite this being a lossy medium. Clearly this is physically invalid, so the constant complex CR model has failed. This has a crucial significance for any technology relying on temporal boundaries. For example, consider an EM wave propagating in an engineered material based on acrylonitrile butadiene styrene (ABS) whose permittivity is changed at time t 0 from that of ABS ( = 2.55 + i0.007) to = 2 + i0.008; according to a constant complex permittivity model the amplitude of the EM wave is predicted to amplify exponentially, even though the material stays lossy. Further, substituting (15) into (16) we observe that in general E x (t, z) is complex valued, despite the TBC (9) and the initial field (10) being real-valued. This is because the wave equation for the medium when t > 0, from (5) and J total = 0, is i.e. not a real equation in real unknowns; a point whose significance might be missed if expecting to take the real part.

V. USING A NARROWBAND APPROXIMATION
Using a narrowband approximation we can recover the correct physical behaviour, although when solving the dispersion (2) one needs to be careful about the square root. We choose ω, k in (1) to satisfy the dispersion relation and consider all relevant frequencies. With an over-bar denoting complex conjugates, the associated solutions are created by substituting ω → ±ω, ω → ±ω, k → ±k, and k → ±k. Then, combining (2) and (3) gives c 2 0 k 2 − ω 2ñ (−ω) 2 = 0. The eight roots of this and its complex conjugate are given by (see Supplementary Material) which usesñ(−ω) =ñ(ω), a consequence of the reality condition onñ. Now, if we apply the NBA, we know that all the fields are concentrated about two modes, i.e. at ω ≈ ω 0 and ω ≈ ω 0 . Let N =ñ(−ω 0 ), so that for ω ≈ ω 0 we haveñ(−ω) ≈ N . Similarly, we also have N =ñ(ω 0 ).
To obtain a real E x (t, z) we need to add the complex conjugate. As Im(ω 0 ) = Im(−ω 0 ) we construct the total field from (19) and (22) above to give We are now in a position to reconsider our time boundary system in the context of dispersive media. Before the time boundary t < t b = 0 we have the vacuum. Assuming a narrowband initial pulse enveloping the single mode (10), where ω a = c 0 k a ∈ R, and k a > 0. Whichever TBC we assume we obtain k c = ±k a . Again the choice of root is unimportant, therefore we set k c = k a , so that k c > 0. After the time boundary, we have dispersion relations given by (19)- (22). In order to obtain the physical solutions consistent with k c > 0 we choose Re(N ) > 0, Im(N ) > 0 and the modes given by (19) and (22). For convenience we define a wave speed Supplementary Material), and the EM field remains real and damped. However, we no longer have a single differential equation (17), as we now need two (N , N ): So by using the NBA with dispersion, we can match the TBC, and obtain a physical, dampened, real-valued solution for the electric field; but the two post-boundary refractive indices required suggest that more appropriate models are necessary.

VI. DYNAMIC MATERIAL MODELS
It is not possible to both implement a physically consistent time boundary with a constant complex c , and escape the requirement for the NBA in (25). This means we must instead use an explicitly causal [43] dynamic response model for the medium, thus providing fully dispersive CR. Since a response model adds extra field(s) to describe the material response, the coupled EM-material system will both have more than two modes and need ABCs. Essentially, defining a material response model immediately creates a demand for ABCs -they are simply the boundary conditions on the auxilliary fields (such as polarization P or bound current J b ) that one has decided to use.
A minimal but sufficient material response can be given by  b /(χ0λ); but we show 1 2 EJ b to increase visibility. Total energy is conserved for t < 0; but just after the time boundary, EE (and hence EEM) rapidly reduce as excitation is transferred to Jb, where it is strongly damped. The energy lost though Jb depends on how strongly it is driven by E until the new dynamic near-equilibrium is reached. Direct damping from the lossy ω3 solution is shown by the rapid and significant decay of EJ b just after the transition. so that in the steady state P = χ 0 λ −1 E. In this model, as long as the loss λ is large compared to the field frequency, with the desired change in permittivity ∆ = χ 0 /λ held fixed, it indeed responds as if it were a medium of complex constant CR. This model also requires a boundary condition for the dielectric polarization P field, namely [P ] = 0.
However, it is significant that the polarization field (26) is in fact derived from a bound current J b =Ṗ = −λP + χ 0 E. This J b is driven by, and hence indirectly applies loss to the field E; and our simulation results shown in figure 2 use this J b approach. Used after the time boundary (t > t b ), this medium has a dispersive behaviour given by where χ 0 > 0 and λ > 0; and the steady state is reached when ω/λ → 0. The resulting cubic dispersion relation is This equation in ω and fixed k 2 produces three solutions as opposed to the two given in (11) or (24). Since (28) has real coefficients when written as a polynomial in (iω), and as a cubic has at most two non-real roots for iω, the three solutions consist of a complex conjugate pair and a real valued one. The pair correspond to counter propagating waves (modes) with the same damping, with the other solution (mode) being nonpropagating and purely damped (see figure 2(a)). At a time boundary where k ∈ R and k > 0, we have three outgoing modes, requiring three boundary conditions. Two TBC are given by either the natural (8) or fundamental (9) TBC (7); the ABC of course is just the [P ] = 0 evident from the dynamic model (26). This ABC is analogous to the Pekar ABC [44], which are required when an EM wave passes into a spatially dispersive medium [45,46].
Simulation results are shown in figure 2(b,c,d), demonstrating the system behaviour as it passes the boundary. Despite the excellent match to a medium with constant complex CR before and sufficiently far after the boundary transition, it does not exhibit the unphysical behaviour of prescribed constant complex CRs. Instead, just after t = t b = 0 in figures 2(c,d) we can see a rapid rebalancing as E and J b (i.e.Ṗ ) adjust to the recently changed χ 0 .
Note also the overshoot in E J just after the boundary, which can be attributed to the dampened ω 3 . These occur on a timescale set by λ, and incur an E dependent energy loss. From a dynamic perspective, system appears to have two loss processes, although both are in fact different manifestations of the same λ loss term. In a steady-state medium, i.e. away from the boundary, the losses are gradual and proportional to ω/λ. This counter-intuitive dependence on the 1/λ is because the loss depends on the mismatch between the ideal J b (or polarization P ) and the actual value; but for larger λ values, the mismatch becomes smaller. In contrast, as the system transitions across the time boundary, the sudden change in χ 0 means that the E-dependent mismatch can suddenly become very large, and this causes equally large and rapid losses as J b (or P ) re-synchronise to the electric field E. However, note that in the special case where E is zero at the boundary, the mismatch remains small and no significant rapid loss takes place.
A key further point of interest is whether there are any sideeffects if the model parameters λ and χ 0 change with time. On the basis of (26) this does not appear to be the case, since the only time derivative term is that on the LHS, applied to P . However, to check this properly we need to reformulate the model so that it is explicitly based on bound currents, which are the true microscopic physical property. Taking the time derivative of (26) and rearranging (see Appendix E) leads to the model equatioṅ where K = J b − χ 0 E is an offset for the bound current J b . This construction has ensured that there is only one time derivative applied to one field quantity (i.e. K), so that it retains an unambiguously causal form [43]; it is in fact this equation we integrate, along with Maxwell's equations, to get the results shown in Fig. 1(b-d) We can see from the equation that the value of K is going to be continually trying to catch up to the present value (albeit scaled) of E, on a timescale set by λ. Since the introduction of K removes any dependence on the time derivative of χ 0 , a temporal boundary in the χ 0 value is implemented simply by changing the parameter χ 0 as we integrate step by step. However, extra care need to be taken if λ changes at a temporal boundary; since its time derivativė λ also appears. This would either be specified as part of the simulation environment defining χ 0 , λ,λ; or we might add an auxilliary equation for λ if it had its own temporal dynamics. This model (i.e. (26), or (29)), being defined by a temporal differential equation, is necessarily explicitly causal [43]. It works as intended, i.e. to create a near-constant permittivity, when (i) the desired positive permittivity shift χ 0 > 0, and when (ii) χ 0 /λ 0 3. In this regime the polarization (or microscopic polarization current) adiabatically follows the phase of electric field, and so accurately models the desired effective permittivity; indeed the effective loss at low frequencies is proportional to ω/λ, i.e. increasing the polarization current loss λ actually reduces the effective damping in the CW limit. However, if the first condition does not hold, a positive polarization can have a negative energy; thus the amplitudes of E and P can increase without limit whilst still conserving energy. Alternatively, if the parameters change so that the second condition starts to fail, the three modes -two electromagnetic and one polarization -become ever more strongly coupled and eventually exhibit a complicated dynamics (and dispersion) not relevant to our presentation here.

VII. BEYOND THE MINIMAL DYNAMIC MODEL
The minimal model used above performs well, and has the considerable advantage of having a simple behaviour, thus clarifying the general principles for handling time boundaries. However, it is not typical of the material models used in practical situations which are often based on Drude or Lorentz oscillators that consider the polarization P . Consequently, we now consider a more general situation by using the result that any causal response can be expressed as a sum of Lorentz responses [47]. Working in the time domain, we have D = 0 E + s P s . With s max being the number of oscillators, we have then s max second order equations: Here χ s is the coupling, λ s is the damping, and α s the natural frequency of the oscillator. Note that the time derivative in the second LH term acts on the product λ s P s , not just on P s . Since the Lorentian oscillators follow second order dynamical equations, they will therefore each require two extra boundary conditions, for a total of 2s max ABCs. The natural ABCs for the s-th oscillator are To see this we substitute P (t, x) = θ(t b − t) P a s (t, x) + θ(t − t b ) P c s (t, x) and likewise for λ s , α s and χ s , into (30) and as-sume (30) holds on both sides of the boundary. This gives The coefficients of δ (t b − t) and δ(t b − t) must both be zero. This gives (31). Observe that if [λ s ] = 0 then (31) reduces to [P s ] = 0 and [Ṗ s ] = 0. This result explains why in (30) that λ s is inside the time derivative in the second LH term. If it had been outside the derivative, then if [P s ] = 0, the alternative λ s (P s )˙would contain the product δ(t − t b )θ(t − t b ), which is not defined mathematically. An alternative method for setting out these ABCs would be to factorise (30) into two first order pieces; and it is also possible to follow either type of analysis for a bound current J b [48].
We can also calculate the number of ABCs required in the frequency domain, using the case of no boundary, or when the system is very far from the boundary. Here the parameters χ s , λ s can be treated as constants, where so that (30) gives us the dispersion Expanding out the dispersion relation resulting from (2) then leads to a polynomial in ω of degree (2s max + 2), and hence (2s max + 2) modes. Matching these modes across a boundary between different materials will then require (2s max + 2) boundary conditions, 2s max of which are ABCs.

VIII. CONCLUSION
We have shown that even simple time boundaries in optics cannot be described by the standard "constant complex per-mittivity" model. Only with a dynamic or dispersive model of the propagation medium can physical results be predicted. This conclusion is supported by NBA calculations and time domain simulations, and is easily generalized to other wave systems. It has significant implications for the design and modelling of both experiments and applications of future technologies based on temporal boundary phenomena. Our conclusion may be particularly relevant to the propagation of relativistically intense electromagnetic waves in plasma, the propagation of relativistic plasma waves in laser wakefield accelerators in the bubble regime, and for relativistic ionisation fronts and relativistically induced transparency in laser-solid interactions.
It is arguably unsurprising that a good model of a time boundary requires a model that can admit non-trivial time dependence, i.e. either a time domain response model or a frequency domain dispersion. Time boundaries are a temporal, dynamic phenomena, and need to be treated as such.

ACKNOWLEDGMENTS
Thus we rearrange (E2), and substitute for P using a rearranged (E1), to get where the changes (effects) on the LHS from the first RHS term are determined solely by the known present values of K and E. If the system parameter λ is time dependent (i.e. λ(t)), then see also that we need to know ∂ t λ as well as λ; although such a specification is not required for χ 0 . This is particularly relevant in the case of an (otherwise constant) medium with an abrupt change at a time boundary. From (E6) we can see that the value of K is going to be continually trying to catch up to the present value (albeit scaled) of E, on a timescale set by λ. Clearly, for this model to function as intended we will want the λ timescale to be faster than the largest significant frequency component of E.
This means that this minimal dynamical model for the standard "constant permittivity" assumption which will match the target real-valued permittivity in the large λ limit, with the concommitant introduction of a loss that gets ever smaller in the ideal large λ limit.