Temporal boundaries in electromagnetic materials

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, sensors, accelerators, and even microwave therapy [4]. In such situations, microscopic models involving e.g. atomic structure or quantum mechanical effects typically provide no significant advantage.

As new materials are developed for use in novel devices, it becomes vital that EM loss and propagation are correctly predicted. A very exciting new concept is time dependent media, where abrupt changes in permittivity can create a temporal boundary. Boundaries play a key role in many physical models; they provide initial and final states in dynamical systems, constrain analytic solutions in confined systems, and represent transitions between different modes of operation. These concepts date back to the 1950s, when Morgenthaler [5] showed that a temporal change of the permittivity produces both forward and backward propagating waves; a result echoed in directional formulations for wave propagation [6–10]. Temporal boundaries [5, 11–13] act as time-reversing mirrors in acoustics [14]; in EM an instantaneous time mirror with a sign-change in permittivity has been predicted [15] to cause field amplification. Other examples include dynamically configurable systems [16, 17], spacetime transformation devices [18–21], time crystals [22–24], 'field patterns' [25], and in laser driven plasma where relativistic changes can lead to a spatial and temporally varying permittivity.

Here we prove that modelling loss using a constant complex permittivity and permeability [15, 26–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 permittivity 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–37]; here E and B are also continuous, but a temporal-surface dipole current appears at the transition, as shown in figure 1.

Zoom In Zoom Out Reset image size

After demonstrating and defining the problem, we present two methods for obtaining physically meaningful solutions. The first applies the constant complex CR model while 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 appendices.

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 + i k ⋅ x ), where $\omega \in \mathbb{C}$ is a complex frequency, and $\boldsymbol{k}\in {\mathbb{C}}^{3}$ a complex wavevector. Since the source free Maxwell's equation (5) are linear, with J _total = 0 and ρ = 0, it is sufficient to consider just a single mode

$$\begin{equation}{E}_{x}(t,z)={E}_{0}\enspace \mathrm{exp}(-i\omega t+ikz)\quad \text{where}\enspace k,\omega \in \mathbb{C},\end{equation} \tag{ 1 }$$

with k oriented along the z-axis (along $\hat{\boldsymbol{z}}$), and E along $\hat{\boldsymbol{x}}$.

Now consider a temporally dispersive medium where the CRs specify permittivity $\tilde {{\epsilon}}(-\omega )$ and permeability μ₀, where the '−ω' is a consequence of choosing e^−iωt in (1). If ω and $\tilde {{\epsilon}}(-\omega )$ are both real then $\tilde {{\epsilon}}(-\omega )$ can be replaced with $\tilde {{\epsilon}}(\omega )$; 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–40]. These are

$$\begin{equation}{k}^{2}-{\omega }^{2}\enspace {\mu }_{0}\tilde {{\epsilon}}(-\omega )=0,\end{equation} \tag{ 2 }$$

$$\begin{equation}\tilde {n}{(-\omega )}^{2}={c}_{0}^{2}{\mu }_{0}\tilde {{\epsilon}}(-\omega ),\end{equation} \tag{ 3 }$$

where ${c}_{0}={({{\epsilon}}_{0}{\mu }_{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 $\text{Im}\left(\tilde {{\epsilon}}(-\omega )\right){ >}0$.

Second, if k is real, then we need Im(ω) < 0 to damp the field, leading to the requirement that when ${\omega }^{2}\tilde {{\epsilon}}(-\omega )$ 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.

Different types of EM TBC can be identified, depending on the material response [11, 12, 44]; here we summarise focussing on how bound currents represent the material response. First, as depicted in figure 1, we could identify the sudden change in the CRs at t = t_b as leading to an instantaneous temporal-surface dipole current (${\boldsymbol{\mathcal{J}}}_{ \text{surf}}$). The total current, J _total, in the medium is

$$\begin{equation}{\boldsymbol{J}}_{ \text{total}}={\boldsymbol{J}}_{ \text{reg}}-\delta (t-{t}_{\text{b}}){\boldsymbol{\mathcal{J}}}_{ \text{surf}},\end{equation} \tag{ 4 }$$

where J _reg is the usual finite current in Maxwell's equations:

$$\begin{equation}\begin{gathered}\nabla \cdot \boldsymbol{B}=0,\quad \nabla \cdot \boldsymbol{D}=\rho \\ \nabla {\times}\boldsymbol{E}+{\partial }_{t}\boldsymbol{B}=0\quad \text{and}\quad \nabla {\times}\boldsymbol{H}-{\partial }_{t}\boldsymbol{D}={\boldsymbol{J}}_{ \text{total}}.\end{gathered}\end{equation} \tag{ 5 }$$

Since the fields E , B , D , H may be discontinuous we write

$$\begin{equation}\boldsymbol{D}(t,\boldsymbol{x})=\theta ({t}_{\text{b}}-t){\boldsymbol{D}}_{\text{a}}(t,\boldsymbol{x})+\theta (t-{t}_{\text{b}}){\boldsymbol{D}}_{\text{c}}(t,\boldsymbol{x}),\end{equation} \tag{ 6 }$$

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

$$\begin{align}{\boldsymbol{J}}_{ \text{total}}& ={\boldsymbol{J}}_{ \text{reg}}-\delta (t-{t}_{\text{b}}){\boldsymbol{\mathcal{J}}}_{ \text{surf}}=\nabla {\times}\boldsymbol{H}-{\partial }_{t}\boldsymbol{D}\\ & =\theta ({t}_{\text{b}}-t)\nabla {\times}{\boldsymbol{H}}_{\text{a}}(t,\boldsymbol{x})+\theta (t-{t}_{\text{b}})\nabla {\times}{\boldsymbol{H}}_{\text{c}}(t,\boldsymbol{x})\\ & \quad -\theta ({t}_{\text{b}}-t){\partial }_{t}{\boldsymbol{D}}_{\text{a}}(t,\boldsymbol{x})-\theta (t-{t}_{\text{b}}){\partial }_{t}{\boldsymbol{D}}_{\text{c}}(t,\boldsymbol{x})\\ & \quad -\delta (t-{t}_{\text{b}})\left\{{\boldsymbol{D}}_{\text{c}}(t,\boldsymbol{x})-{\boldsymbol{D}}_{\text{a}}(t,\boldsymbol{x})\right\}.\end{align}$$

This approach can also be used for B in the Maxwell–Faraday equation to derive a similar result. The TBC are

$$\begin{equation}[\boldsymbol{D}]={\boldsymbol{\mathcal{J}}}_{ \text{surf}}\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad [\boldsymbol{B}]=0,\end{equation} \tag{ 7 }$$

where [ D ] = D _c(t_b, x ) − D _a(t_b, x ), etc. One option [5, 15] is to set ${\boldsymbol{\mathcal{J}}}_{ \text{surf}}=0$, to obtain the natural TBC, i.e.

$$\begin{equation}[\boldsymbol{D}]=0\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad [\boldsymbol{B}]=0.\end{equation} \tag{ 8 }$$

Alternatively, if treating E and B as the only physical fields, we get the fundamental TBC, i.e.

$$\begin{equation}[\boldsymbol{E}]=0\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad [\boldsymbol{B}]=0,\end{equation} \tag{ 9 }$$

which rely on (7) to calculate the conserved temporal-surface current ${\boldsymbol{\mathcal{J}}}_{ \text{surf}}$. This is an analogous approach to that describing the surface current around a permanent magnet [35].

Even a time boundary between a vacuum with = ₀ 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 t_b = 0, so that D (t, x ) = ₀ E (t, x ) for t < 0, and D (t, x ) = _c E (t, x ) for t > 0; and then choose a field polarization so that $\boldsymbol{E}={E}_{x}(t,z)\hat{\boldsymbol{x}}$, and $\boldsymbol{B}={B}_{y}(t,z)\hat{\boldsymbol{y}}$.

Pre-boundary (t < 0), we start with a single real mode

$$\begin{equation}{E}_{x}(t,z)={E}_{0}\enspace \mathrm{cos}({\omega }_{\text{a}}t-{k}_{\text{a}}z),\end{equation} \tag{ 10 }$$

with ${E}_{0}\in \mathbb{R}$ and B_y (t, z) = E_x /c₀. Here ω_a, k_a are both real and positive and satisfy the vacuum dispersion relation ${c}_{0}^{2}{k}_{\text{a}}^{2}={\omega }_{\text{a}}^{2}$. Since the post-boundary lossy material, with ω_c and k_c, must have Im(_c) > 0, the t > 0 general solution is

$$\begin{equation}\begin{aligned}{E}_{x}(t,z)& =\left({g}_{+}^{-}\enspace {\text{e}}^{-\text{i}{\omega }_{\text{c}}t+i{k}_{\text{c}}z}+{g}_{-}^{-}\enspace {\text{e}}^{-\text{i}{\omega }_{\text{c}}t-i{k}_{\text{c}}z}+{g}_{-}^{+}\enspace {\text{e}}^{\text{i}{\omega }_{\text{c}}t-i{k}_{\text{c}}z}+{g}_{+}^{+}\enspace {\text{e}}^{\text{i}{\omega }_{\text{c}}t+i{k}_{\text{c}}z}\right),\\ {B}_{y}(t,z)& =\frac{{k}_{\text{c}}}{{\omega }_{\text{c}}}\left({g}_{+}^{-}\enspace {\text{e}}^{-\text{i}{\omega }_{\text{c}}t+i{k}_{\text{c}}z}-{g}_{-}^{-}\enspace {\text{e}}^{-\text{i}{\omega }_{\text{c}}t-i{k}_{\text{c}}z}+{g}_{-}^{+}\enspace {\text{e}}^{\text{i}{\omega }_{\text{c}}t-i{k}_{\text{c}}z}-{g}_{+}^{+}\enspace {\text{e}}^{\text{i}{\omega }_{\text{c}}t+i{k}_{\text{c}}z}\right),\end{aligned}\end{equation} \tag{ 11 }$$

with the dispersion relation

$$\begin{equation}{k}_{\text{c}}{c}_{0}={n}_{\text{c}}{\omega }_{\text{c}}\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\quad {n}_{\text{c}}={c}_{0}{({{\epsilon}}_{\text{c}}{\mu }_{0})}^{1/2}.\end{equation} \tag{ 12 }$$

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

$$\begin{equation}{E}_{x}({0}^{-},z)=\frac{1}{2}{E}_{0}\left({\text{e}}^{\text{i}{k}_{\text{a}}z}+{\text{e}}^{-\text{i}{k}_{\text{a}}z}\right),\end{equation} \tag{ 13 }$$

and B_y (0⁻, z) = (k_a/ω_a)E_x . Just after the time boundary

$$\begin{equation}\begin{aligned}{E}_{x}({0}^{+},z)& =({g}_{-}^{+}+{g}_{-}^{-}){\text{e}}^{-\text{i}{k}_{\text{c}}z}+({g}_{+}^{+}+{g}_{+}^{-}){\text{e}}^{\text{i}{k}_{\text{c}}z},\\ {B}_{y}({0}^{+},z)& =\frac{{k}_{\text{c}}}{{\omega }_{\text{c}}} \left[({g}_{-}^{+}-{g}_{-}^{-}){\text{e}}^{-\text{i}{k}_{\text{c}}z}+({g}_{+}^{-}-{g}_{+}^{+}){\text{e}}^{\text{i}{k}_{\text{c}}z}\right] .\end{aligned}\end{equation} \tag{ 14 }$$

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 to get

$$\begin{equation}{g}_{-}^{+}={g}_{+}^{-}=\frac{1}{4}{E}_{0}(1+{n}_{\text{c}})\quad \text{and}\quad {g}_{-}^{-}={g}_{+}^{+}=\frac{1}{4}{E}_{0}\left(1-{n}_{\text{c}}\right).\end{equation} \tag{ 15 }$$

In general, for t > 0 all the ${g}_{{\pm}}^{{\pm}}$ coefficients of E_x (t, z) are non zero. Let c_c = c_R + ic_I = c₀/n_c then c_R > 0 and c_I < 0, and expand (11) to yield

$$\begin{equation}{E}_{x}(t,z)={g}_{+}^{-}\enspace {\text{e}}^{\text{i}{k}_{\text{c}}(-{c}_{\text{R}}t+z)}{\text{e}}^{{k}_{\text{c}}{c}_{\text{I}}t}+{g}_{-}^{-}\enspace {\text{e}}^{\text{i}{k}_{\text{c}}(-{c}_{\text{R}}t-z)}{\text{e}}^{{k}_{\text{c}}{c}_{\text{I}}t}+{g}_{-}^{+}\enspace {\text{e}}^{\text{i}{k}_{\text{c}}({c}_{\text{R}}t-z)}{\text{e}}^{-{k}_{\text{c}}{c}_{\text{I}}t}+{g}_{+}^{+}\enspace {\text{e}}^{\text{i}{k}_{\text{c}}({c}_{\text{R}}t+z)}{\text{e}}^{-{k}_{\text{c}}{c}_{\text{I}}t}.\end{equation} \tag{ 16 }$$

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₀ 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. Furthermore, 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

$$\begin{equation}{\nabla }^{2}\boldsymbol{E}-{{\epsilon}}_{\text{c}}\ddot {\boldsymbol{E}}=0,\end{equation} \tag{ 17 }$$

i.e. not a real equation in real unknowns; a point whose significance might be missed if expecting to take the real part.

Using a NBA 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

$$\begin{equation}k\enspace {c}_{0}-\tilde {n}(-\omega )\omega =0,\end{equation} \tag{ 18 }$$

and consider all relevant frequencies. With an over-bar denoting complex conjugates, the associated solutions are created by substituting ω → ±ω, $\omega \to {\pm}\overline{\omega }$, k → ±k, and $k\to {\pm}\overline{k}$. Then, combining (2) and (3) gives ${c}_{0}^{2}{k}^{2}-{\omega }^{2}\tilde {n}{(-\omega )}^{2}=0$. The eight roots of this and its complex conjugate are given by

$$\begin{equation}{E}_{x}=\mathrm{exp}(-i\omega t{\pm}ikz)\quad \text{satisfies}\;\;\tilde {n}(-\omega )\enspace \omega ={\pm}{c}_{0}k,\end{equation} \tag{ 19 }$$

$$\begin{equation}{E}_{x}=\mathrm{exp}(i\omega t\mp ikz)\quad \text{satisfies}\;\;\tilde {n}(\omega )\enspace \omega ={\pm}{c}_{0}k,\end{equation} \tag{ 20 }$$

$$\begin{equation}{E}_{x}=\mathrm{exp}(-i\overline{\omega }t{\pm}i\overline{k}z)\quad \text{satisfies}\;\;\overline{\tilde {n}(\omega )}\overline{\omega }={\pm}{c}_{0}\overline{k},\end{equation} \tag{ 21 }$$

$$\begin{equation}{E}_{x}=\mathrm{exp}(i\overline{\omega }t\mp i\overline{k}z)\quad \text{satisfies}\;\;\tilde {n}(\overline{\omega })\overline{\omega }={\pm}{c}_{0}\overline{k},\end{equation} \tag{ 22 }$$

(see appendix C) which uses $\tilde {n}(-\overline{\omega })=\overline{\tilde {n}(\omega )}$, a consequence of the reality condition on $\tilde {n}$.

Now, if we apply the NBA, we know that all the fields are concentrated about two modes, i.e. at ω ≈ ω₀ and $\omega \approx {\overline{\omega }}_{0}$. Let $N=\tilde {n}(-{\omega }_{0})$, so that for ω ≈ ω₀ we have $\tilde {n}(-\omega )\approx N$. Similarly, we also have $\overline{N}=~{n}({\omega }_{0})$.

To obtain a real E_x (t, z) we need to add the complex conjugate. As $\text{Im}({\omega }_{0})=\text{Im}(-{\overline{\omega }}_{0})$ we construct the total field from (19) and (22) above to give

$$\begin{equation}{E}_{x}(t,z)={g}_{+}^{-}\enspace {\text{e}}^{-\text{i}\omega t+ikz}+{g}_{-}^{-}\enspace {\text{e}}^{-\text{i}\omega t-ikz}+{g}_{-}^{+}\enspace {\text{e}}^{\text{i}\overline{\omega }t-i\overline{k}z}+{g}_{+}^{+}\enspace {\text{e}}^{\text{i}\overline{\omega }t+i\overline{k}z}.\end{equation} \tag{ 23 }$$

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 ${\omega }_{\text{a}}={c}_{0}{k}_{\text{a}}\in \mathbb{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 C = C_R + iC_I = c₀/N, where ${C}_{\text{R}},{C}_{\text{I}}\in \mathbb{R}$, C_R > 0 and C_I < 0. From (23) we have,

$$\begin{equation}{E}_{x}(t,z)={E}_{0}\enspace {\text{e}}^{k{C}_{\text{I}}t}\left(\frac{N+1}{N+\overline{N}}\enspace {\text{e}}^{-\text{i}k({C}_{\text{R}}t-z)}+\frac{\overline{N}+1}{N+\overline{N}}\enspace {\text{e}}^{-\text{i}k({C}_{\text{R}}t+z)}\right),\end{equation} \tag{ 24 }$$

where k = k_c > 0 (see appendix), and the EM field remains real and damped. However, we no longer have a single differential equation (17), as we now need two (N, $\overline{N}$):

$$\begin{equation}\begin{cases}{{c}_{0}}^{2}{\nabla }^{2}\boldsymbol{E}-{N}^{2}\ddot {\boldsymbol{E}}=0\quad & \;\text{for}\;\text{frequencies}\;\omega \approx {\omega }_{0}\enspace \text{in}\enspace \boldsymbol{E},\\ {{c}_{0}}^{2}{\nabla }^{2}\boldsymbol{E}-{\overline{N}}^{2}\ddot {\boldsymbol{E}}=0\quad & \text{for}\;\text{frequencies}\enspace {\omega}\approx \overline{{\omega }_{0}}\enspace \text{in}\enspace \boldsymbol{E}.\end{cases}\end{equation} \tag{ 25 }$$

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.

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 [42] 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

$$\begin{equation}\boldsymbol{D}={{\epsilon}}_{0}\boldsymbol{E}+\boldsymbol{P}\quad \text{where}\quad \dot {\boldsymbol{P}}=-\lambda \boldsymbol{P}+{\chi }_{0}\boldsymbol{E},\end{equation} \tag{ 26 }$$

so that in the steady state P = χ₀ λ⁻¹ E . In this model, as long as the loss λ is large compared to the field frequency, with the desired change in permittivity Δ = χ₀/λ 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 ${\boldsymbol{J}}_{\text{b}}=\dot {\boldsymbol{P}}=-\lambda \boldsymbol{P}+{\chi }_{0}\boldsymbol{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

$$\begin{equation}\tilde {{\epsilon}}(-\omega )={{\epsilon}}_{0}+\frac{{\chi }_{0}}{\lambda -i\omega },\end{equation} \tag{ 27 }$$

where χ₀ > 0 and λ > 0; and the steady state is reached when ω/λ → 0. The resulting cubic dispersion relation is

$$\begin{equation}{{\epsilon}}_{0}{\mu }_{0}\left(\lambda -i\omega \right){\omega }^{2}+{\mu }_{0}{\chi }_{0}{\omega }^{2}-{k}^{2}\left(\lambda -i\omega \right)=0.\end{equation} \tag{ 28 }$$

This equation in ω and fixed k² 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 non-propagating and purely damped (see figure 2(a)).

Zoom In Zoom Out Reset image size

At a time boundary where $k\in \mathbb{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 [45], which are required when an EM wave passes into a spatially dispersive medium [46, 47].

Simulation results are shown in figures 2(b), (c) and (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) and (d) we can see a rapid rebalancing as E and J _b (i.e. $\dot {\boldsymbol{P}}$) adjust to the recently changed χ₀.

Note also the overshoot in ${\mathcal{E}}_{J}$ just after the boundary, which can be attributed to the dampened ω₃. 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 χ₀ 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 side-effects if the model parameters λ and χ₀ 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 equation

$$\begin{equation}\dot {\boldsymbol{K}}=-\lambda \left[\boldsymbol{K}+{\chi }_{0}\boldsymbol{E}\right]+\frac{\dot {\lambda }}{\lambda }\boldsymbol{K},\end{equation} \tag{ 29 }$$

where K = J _b − χ₀ 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 [42]; it is in fact this equation we integrate, along with Maxwell's equations, to get the results shown in figures 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 χ₀, a temporal boundary in the χ₀ value is implemented simply by changing the parameter χ₀ as we integrate step by step. However, extra care need to be taken if λ changes at a temporal boundary; since its time derivative $\dot {\lambda }$ also appears. This would either be specified as part of the simulation environment defining χ₀, λ, $\dot {\lambda }$; 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 [42]. It works as intended, i.e. to create a near-constant permittivity, when (i) the desired positive permittivity shift χ₀ > 0, and when (ii) χ₀/λ₀ ⪅ 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 while still conserving energy. Alternatively, if the parameters change so that the second condition starts to fail, the three modes—two EM and one polarization—become ever more strongly coupled and eventually exhibit a complicated dynamics (and dispersion) not relevant to our presentation here.

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 [43]. Working in the time domain, we have D = ₀ E + ∑_s P _s . With s_max being the number of oscillators, we have then s_max second order equations:

$$\begin{equation}{\ddot {\boldsymbol{P}}}_{s}+\left({\lambda }_{s}{\boldsymbol{P}}_{s}\right)\dot {\enspace }+{\alpha }_{s}{\boldsymbol{P}}_{s}={\chi }_{s}\boldsymbol{E}.\end{equation} \tag{ 30 }$$

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 Lorentzian 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 sth oscillator are

$$\begin{equation}[{\boldsymbol{P}}_{s}]=0\quad \mathrm{a}\mathrm{n}\mathrm{d}\quad [2{\dot {\boldsymbol{P}}}_{s}+{\lambda }_{s}\enspace {\boldsymbol{P}}_{s}]=0.\end{equation} \tag{ 31 }$$

To see this we substitute $\boldsymbol{P}(t,\boldsymbol{x})=\theta ({t}_{\text{b}}-t)\enspace {\boldsymbol{P}}_{s}^{\text{a}}(t,\boldsymbol{x})+\theta (t-{t}_{\text{b}})\enspace {\boldsymbol{P}}_{s}^{\text{c}}(t,\boldsymbol{x})$ and likewise for λ_s , α_s and χ_s , into (30) and assume (30) holds on both sides of the boundary. This gives

$$\begin{equation*}{\delta }^{\prime }({t}_{\text{b}}-t)\enspace \left({\boldsymbol{P}}_{s}^{\text{c}}-{\boldsymbol{P}}_{s}^{\text{a}}\right)+2\enspace \delta ({t}_{\text{b}}-t)\enspace \left({\dot {\boldsymbol{P}}}_{s}^{\text{c}}-{\dot {\boldsymbol{P}}}_{s}^{\text{a}}\right)+\delta ({t}_{\text{b}}-t)\enspace \left({\lambda }_{s}^{\text{c}}{\boldsymbol{P}}_{s}^{\text{c}}-{\lambda }_{s}^{\text{a}}{\boldsymbol{P}}_{s}^{\text{a}}\right)=0.\end{equation*}$$

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 $[{\dot {\boldsymbol{P}}}_{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 ${\lambda }_{s}({\boldsymbol{P}}_{s})\dot {\enspace }$ 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 [50].

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

$$\begin{equation}\tilde {{\epsilon}}(-\omega )={{\epsilon}}_{0}+\sum\limits _{s=1}^{{s}_{\mathrm{max}}}\frac{{\chi }_{s}}{-{\omega }^{2}-{\imath}{\lambda }_{s}\omega +{\alpha }_{s}^{2}}.\end{equation} \tag{ 32 }$$

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.

We have shown that even simple time boundaries in optics cannot be described by the standard 'constant complex permittivity' 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 EM 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.

JG and PK acknowledge support provided by STFC (Cockcroft Institute, ST/G008248/1 and ST/P002056/1), and with DAJ acknowledge EPSRC (Lab in a Bubble, EP/N028694/1). DAJ also acknowledges the EU's Horizon 2020 programme (No. 871124 Laserlab-Europe), RS thanks AFOSR (FA8655-20-1-7002), and PK acknowledges recent support from EPSRC (QUANTIC, EP/T00097X/1).

The data that support the findings of this study are available upon reasonable request from the authors.

Lemma 1. Given that for $\omega \in \mathbb{R}$, ω > 0 plane waves are spatially evanescent in the propagation direction then $\text{Im}\left(\tilde {{\epsilon}}(-\omega )\right){ >}0$.

Proof. Given (1) with ω > 0. For Re(k) > 0 then the direction of propagation is positive z. If the plane waves are evanescent for positive z then Re(ikz) < 0 hence Im(k) > 0. This implies k lies in the top right quadrant of $\mathbb{C}$. Hence Im(k²) > 0.

Likewise for Re(k) < 0 then the direction of propagation is negative z. If the plane waves are evanescent for negative z then Re(ikz) < 0 hence Im(k) < 0. This implies k lies in the bottom left quadrant of $\mathbb{C}$. Hence Im(k²) > 0.

In both cases Im(k²) > 0 and since ω² > 0 then (2) implies $\text{Im}\left(\tilde {{\epsilon}}(-\omega )\right){ >}0$. □

From (3), $\text{Im}\left(\tilde {{\epsilon}}(-\omega )\right){ >}0$ implies that $\tilde {n}(-\omega )$ is either (a) in the top right quadrant of the complex plane $\left\{\text{Re}(\tilde {n}){ >}0\text{and}\text{Im}(\tilde {n}){ >}0\right\}$, or (b) in the bottom left quadrant $\left\{\text{Re}(\tilde {n}){< }0\text{and}\text{Im}(\tilde {n}){< }0\right\}$. Consequently, having $\text{Re}\left(\tilde {n}(-\omega )\right){< }0$ does not contradict our assumption that the real part of both permittivity and permeability are positive: it is still possible to have a negative index of refraction [48, 49].

Demonstration of (19)-(22). Since we chose (1) to satisfy (18) then this is (19) for E_x = exp(−iωt + ikz). Replacing k → = k then gives (19).

Replacing ω → −ω and k → −k then gives (21).

Taking the complex conjugate of (19) and using $\tilde {n}(-\overline{\omega })=\overline{\tilde {n}(\omega )}$ gives (22).

Replacing ω → −ω and k → −k in (22) then gives (21). □

Proof of (24). From (23) we have

$$\begin{align}{B}_{y}(t,z)& =\frac{{\text{e}}^{k{C}_{\text{I}}t}}{{C}_{\text{I}}-i{C}_{\text{R}}}\left({g}_{+}^{-}\enspace {\text{e}}^{-\text{i}k({C}_{\text{R}}t-z)}-{g}_{-}^{-}\enspace {\text{e}}^{-\text{i}k({C}_{\text{R}}t+z)}\right)+\frac{{e}^{k{C}_{\text{I}}t}}{{C}_{\text{I}}+i{C}_{\text{R}}}\left({g}_{-}^{+}\enspace {\text{e}}^{\text{i}k({C}_{\text{R}}t-z)}-{g}_{+}^{+}\enspace {\text{e}}^{\text{i}k({C}_{\text{R}}t+z)}\right)\\ & =\frac{\overline{N}\enspace {\text{e}}^{k{C}_{\text{I}}t}}{{c}_{0}}\left({g}_{+}^{-}\enspace {\text{e}}^{-\text{i}k({C}_{\text{R}}t-z)}-{g}_{-}^{-}\enspace {\text{e}}^{-\text{i}k({C}_{\text{R}}t+z)}\right)+\frac{N\enspace {e}^{k{C}_{\text{I}}t}}{{c}_{0}}\left({g}_{-}^{+}\enspace {\text{e}}^{\text{i}k({C}_{\text{R}}t-z)}-{g}_{+}^{+}\enspace {\text{e}}^{\text{i}k({C}_{\text{R}}t+z)}\right).\end{align} \tag{ D1 }$$

Using the fundamental TBC (9) we get for t = 0⁺ and $k\in \mathbb{R}$

$$\begin{align*}{E}_{x}({0}^{+},z)& ={\text{e}}^{\text{i}kz}({g}_{+}^{-}+{g}_{+}^{+})+{\text{e}}^{-\text{i}kz}({g}_{-}^{-}+{g}_{-}^{+})\\ {c}_{0}{B}_{y}({0}^{+},z)& ={\text{e}}^{\text{i}kz}(\overline{N}{g}_{+}^{-}-N{g}_{+}^{+})+{\text{e}}^{-\text{i}kz}(N{g}_{-}^{+}-\overline{N}{g}_{-}^{-})\end{align*}$$

and hence

$$\begin{equation*}{g}_{+}^{-}=\overline{{g}_{-}^{+}}=\frac{{E}_{0}}{2}\frac{N+1}{N+\overline{N}}\quad \text{and}\quad {g}_{-}^{-}=\overline{{g}_{+}^{+}}=\frac{{E}_{0}}{2}\frac{\overline{N}+1}{N+\overline{N}}\end{equation*}$$

Substituting into (23) gives (24). □

The concept behind our minimal model is that it should mimic the 'ideal' of a constant-like permittivity with real and imaginary components. Thus we want there to exist a dielectric polarization P that closely follows the current value of the electric field E , but allows freedom for dynamical variation about its chosen target value. Most simply, we can write

$$\begin{equation}{\partial }_{t}\boldsymbol{P}=-{\lambda}\boldsymbol{P}+{\chi }_{0}\boldsymbol{E},\end{equation} \tag{ E1 }$$

where in the steady-state limit we have the desired P = (χ₀/λ) E = Δɛ E . In this model, the dynamics are simply that P exponentially decays towards to Δɛ E with rate λ.

However, from a microscopic point of view, a charge or current is the more useful physical property. Thus to adapt the initial concept in (E1) we use the fact that J _b = ∂_t P . However, since substitution of J _b into (E1) leaves us with no dynamics, we also apply an extra time derivative to the above equation. As a result we have

$$\begin{equation}{\partial }_{t}{\boldsymbol{J}}_{b}=-{\lambda}{\boldsymbol{J}}_{b}+{\partial }_{t}\left({\chi }_{0}\boldsymbol{E}\right)-\left({\partial }_{t}{\lambda}\right)\boldsymbol{P}.\end{equation} \tag{ E2 }$$

Since this has time derivatives of fields on both sides, it is not straightforward to interpret it in a causal manner [42]. Thus we combine the J_b and E fields, both of which are subject to an applied time derivative, into a single quantity:

$$\begin{equation}\boldsymbol{K}={\boldsymbol{J}}_{b}-{\chi }_{0}\boldsymbol{E}.\end{equation} \tag{ E3 }$$

Thus we rearrange (E2), and substitute for P using a rearranged (E1), to get

$$\begin{equation}{\partial }_{t}\left({\boldsymbol{J}}_{b}-{\chi }_{0}\boldsymbol{E}\right)=-{\lambda}{\boldsymbol{J}}_{b}-\left({\partial }_{t}{\lambda}\right)\frac{1}{{\lambda}}\left[{\chi }_{0}\boldsymbol{E}-{\partial }_{t}\boldsymbol{P}\right],\end{equation} \tag{ E4 }$$

$$\begin{equation}{\partial }_{t}\boldsymbol{K}=-{\lambda}{\boldsymbol{J}}_{b}+\left({\partial }_{t}{\lambda}\right)\frac{1}{{\lambda}}\left[{\boldsymbol{J}}_{b}-{\chi }_{0}\boldsymbol{E}\right],\end{equation} \tag{ E5 }$$

$$\begin{equation}{\partial }_{t}\boldsymbol{K}=-{\lambda}\left[\boldsymbol{K}+{\chi }_{0}\boldsymbol{E}\right]+\frac{{\partial }_{t}{\lambda}}{{\lambda}}\boldsymbol{K}.\end{equation} \tag{ E6 }$$

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 χ₀. 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 .

In the harmonic case, we can Fourier transform the evolution equation; here a time-domain field A (t) becomes a frequency-domain $\tilde {A}(\omega )$. Then (E6) tells us that

$$\begin{equation}-{\imath}\omega ~{\mathbf{K}}=-{\lambda}\left(~{\boldsymbol{K}}+{\chi }_{0}~{\mathbf{E}}\right),\end{equation} \tag{ E7 }$$

$$\begin{equation}\mathrm{i}.\mathrm{e}.\quad ~{\mathbf{K}}\mathbf{=}\frac{-{\chi }_{0}{\lambda}}{-{\imath}\omega +{\lambda}}~{\boldsymbol{E}}\quad \mathrm{o}\mathrm{r}\quad {~{\boldsymbol{J}}}_{b}={\chi }_{0}\left[\frac{-{\lambda}}{{\lambda}-{\imath}\omega }+1\right]~{\boldsymbol{E}}.\end{equation} \tag{ E8 }$$

In the intended limit where λ ≫ ω, we can expand to second order so that

$$\begin{equation}{\tilde {J}}_{b}\simeq -{\imath}\omega {\chi }_{0}{{\lambda}}^{-1}\left[1-{\imath}\omega {{\lambda}}^{-1}\right]\tilde {E},\end{equation} \tag{ E9 }$$

and as a result we have

$$\begin{equation}~{\boldsymbol{P}}=\left[{\Delta}{\varepsilon }_{r}-{\imath}{\Delta}{\varepsilon }_{i}(\omega )\right]~{\mathbf{E}},\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\quad {\Delta}{\varepsilon }_{i}(\omega )={\Delta}{\varepsilon }_{r}\omega /{\lambda}.\end{equation} \tag{ E10 }$$

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.

Across our time boundary, a change in constant complex CR for a medium is just an instantaneous change, which is straightforwardly causal. Causal behaviour, however, of itself is not necessarily guaranteed to give physically reasonable predictions. Indeed, all but one of the results in this paper are causal—even the unphysical (17), where the future behaviour is by construction explicitly dependent on the past behaviour. The single exception is the narrowband result (25), because under this approximation questions of causality are moot.

Causality is often tested by applying the Kramers–Kronig relations (see e.g. [42]), but they do not apply to all situations. For example, even though the result (16) is causal, and clearly so when solving in the time domain, as an exponentially increasing function it cannot be Fourier transformed so as to allow Kramers–Kronig to be tested. Indeed, since the original function of Kramers–Kronig was to analyse, test, or correct raw data collected in the frequency domain, using them as a causality test when a time domain description is already available is redundant.

This is why the dynamical model (26) is a natural starting point for an examination of temporal boundaries; although of course more complicated models, such as the summed Lorentzians of (32), or ones involving (causal) integral kernels, can be constructed. Notably, even the highly simplified (26), designed to give results that in the appropriate limit are as close as possible to the failed constant complex CR model, is sufficient to restore physical behaviour.

In principle, we could try to expressed this time boundary situation as a Cauchy problem—i.e. namely asking what is the subsequent wave behaviour if the Cauchy (initial condition) data is E(z) = E₀ cos(k_a z). However, it is not clear how this extended field could be set up as a realistic initial condition; and, further, trying to use a straightforward convolution approach would fail because E(z, t) for t < 0 (i.e. the pre-boundary field) would be unknown.

For convolution transforms, we note that their standard use in CRs is of the time-independent form D(t) = ∫(t − t')E(t')dt', but unfortunately, this would only be valid—obviously—for time independent media. For time dependent media, such as one including a time boundary, we would need to use the two-time generalisation D(t) = ∫(t, t')E(t')dt'.

While in principle such generalisations can be written down, it is unfortunate that the time-independent convolution does not indicate in any way how a general (two-time) formulation could be arrived at, and this is complicated further since the convolution must include the time boundary itself. This is in stark contrast to dynamic response models such as (26), or the more general form given here in the appendix as (32), where insisting on (e.g.) continuity of the auxiliary fields is quite natural.

Ohmic losses are also covered by our analysis; as they can also be modelled by a complex permittivity. Since the polarisation corresponding to ohmic losses is given by P = (σ/(−iω)) E then this is an alternative dispersive CR, which means it does not contradict the statements about constant _c.