Nonlocal effects in temporal metamaterials

2.1 Nonlocal effective medium theory

Let us consider an isotropic, generally inhomogeneous, time-modulated medium, whose EM response is described by a relative dielectric permittivity periodically modulated in time, ɛ(r, t) = ɛ(r, t + τ). In our study, we assume that the operating frequencies are much lower than any material resonance frequencies, so that temporal dispersion effects can be approximately neglected [17, 39]. From Maxwell’s equations, the field dynamics can be described by the vector wave equation for the electric induction D, namely

where L ̂ = ∇ × ∇ × is the double-curl operator, and c denotes the wavespeed in vacuum [31]. In the presence of a rapid temporal modulation (i.e., with a modulation angular frequency Ω = 2π/τ much higher than the EM carrier one ω), it is convenient to introduce the parameter η = ω/Ω. By exploiting the standard multiscale approach [6], we can develop an asymptotic analysis in the regime where η ≪ 1. It is natural to assume for the EM observables a separate dependence on the slow and fast temporal scales (t and T = t/η, respectively), and to represent them as a Maclaurin series expansion in the scale-parameter η, i.e.,

with A ̄ = ∑ m = 0 + ∞ η m A ̄ ( m ) and A ̃ = ∑ m = 0 + ∞ η m A ̃ ( m ) . Here, A represents a generic field component, the superscript (m) labels the order of each term, whereas the overline and the tilde label the averaged and rapidly varying contributions to each order, respectively. By noticing that the relative dielectric permittivity only depends on the fast coordinate T, and by substituting the multiscale series expansion in Eq. (1) (up to the second order in η), we obtain

where we have separated the fast and slow contributions, and we have assumed that the reciprocal of the relative permittivity admits the Fourier series expansion

with ε eff − 1 = a 0 , Ω ̃ = η Ω , and i denoting the imaginary unit. Equivalently, we can write Maxwell’s equations governing the dynamics of the average fields in our temporal metamaterial as

along with the constitutive relationships

which are evidently nonlocal because of the presence of field derivatives. By considering the limit η → 0, it is apparent that the term proportional to η ² in the first of Eq. (3) vanishes, and the temporal metamaterial behaves as a dielectric medium with an effective relative permittivity ɛ _eff, thereby recovering the results in previous studies [30, 31]. In the simplest case of spatial homogeneity (i.e., spatially unbounded, time-varying media [30]), the first of Eq. (6) becomes

where K = 2π/(cτ) and γ = 2 ∑ n = 1 + ∞ | a n | 2 / n 2 .

Next, we consider a monochromatic plane wave propagating in a temporal metamaterial, i.e., D ̄ = 2 R e d 0 e i k z − i ω ( k ) t . By substituting this expression in the effective Maxwell’s equations, and exploiting the constitutive relationship in Eq. (7), we derive the dispersion relationship

where

plays the role of an effective relative magnetic permeability. Note that artificial magnetism in time-modulated dielectric media has been predicted in a recent study [31], suggesting that a temporal metamaterial can behave as a resonant magneto-dielectric medium. Our results above confirm the previous findings, and establish tighter bounds on the attainable magnetic response. Specifically, for a temporally periodic metamaterial based on positive-permittivity modulation [ɛ(t) = ɛ(t + τ) > 0], it is apparent from Eq. (9) (recalling that γ ≥ 0) that only nonlocal diamagnetism with 0 < μ _eff(k) < 1 can be attained, whereas resonant, paramagnetic, and μ-negative responses are forbidden.

We highlight that, as in conventional spatially modulated metamaterials, optical magnetism stems from spatial dispersion. As previously mentioned, in spatial metamaterials, nonlocality can strongly affect the EM response, and it can produce undesirable effects [7]. Conversely, in an isotropic temporal metamaterial, the spatial dispersion up to the second order (i.e., up to η ²) is fully equivalent to optical magnetism described by the effective relative magnetic permeability μ _eff(k) in Eq. (9) (see the Methods Section 4.1 for further details).

For some quantitative assessments, we assume that the relative permittivity profile is given by

with Δ parameterizing the modulation depth, and ϕ being an arbitrary phase. To validate and calibrate the predictions of our proposed nonlocal EMT model, we compare them with the results from a rigorous full-wave theory (see the Methods Section 4.1 for further details), as well as with the conventional (local) EMT [16]. Figure 1 compares these three predictions for the dispersion relationship ω(k), for ε ̄ = 5 , and three different values of the modulation depth Δ. Here and henceforth, consistently with our assumption of a negligible temporal dispersion, parameters are chosen so that the relative permittivity in time is always greater than one. Note that the effective parameters ɛ _eff and μ _eff(k) do not depend on ϕ, since this latter merely corresponds to a time shift. As expected, our nonlocal EMT is in fair agreement with the full-wave results and the conventional (local) EMT in the quasi-homogenized (weak-dispersion) regime k/K ≪ 1, and we observe that it works well within the regions ω(k)/Ω < 0.45, 0.47, 0.49 for Δ = 0.5, 0.65, 0.8, respectively, as shown in Figure 1. Note that ω(k)/Ω corresponds to the scale parameter η that rules the multiscale homogenization process (see Section 2.1). By increasing the modulation depth Δ, we expect the nonlocal effects to become more pronounced. This is evident in the three panels of Figure 1, as the wavenumber region where the temporal metamaterial exhibits a homogeneous behavior [i.e., where ω(k) is approximately linear] progressively shrinks as Δ increases.

Figure 1:

Comparison between the predictions from the full-wave theory and both local an nonlocal EMT models for a temporal metamaterial with relative permittivity as in Eq. (10), with ε ̄ = 5 .

(A)–(C) Normalized angular frequency ω/Ω as function of the normalized wavenumber k/K, for Δ = 0.5, 0.65, and 0.8, respectively.The light-blue shaded area indicates the region where the nonlocal EMT works well (≲10% error with respect to full-wave theory).

Figure 2(A) shows the effective relative permittivity ɛ _eff as a function of Δ, whereas Figure 2(B) shows the effective relative permeability μ _eff(k) as a function of k/K, for different values of Δ; these effective parameters are only shown within the region where the nonlocal EMT is in good agreement with the full-wave theory (corresponding to the light-blue shaded area in Figure 1). We observe that, in the case of a deep temporal modulation (e.g., Δ = 0.8), the metamaterial exhibits a significant diamagnetic response (e.g., μ _eff ≃ 0.91 for k/K = 0.73).

Figure 2:

Effective parameters for the temporal metamaterial considered in Figure 1.

(A) Effective relative permittivity ϵ _eff as a function of Δ; (B) effective relative permeability μ _eff(k) as a function of k/K, for different values of the modulation depth Δ. Note that μ _eff(k) is only shown within the region where the nonlocal EMT holds (≲10% error with respect to full-wave theory).

2.2 Nonlocal temporal boundary

In analogy with the spatial scenario, the results derived in Section 2.1 play the role of the “bulk” response of an infinite-extent metamaterial. Since a series of recent studies on spatially truncated multilayered dielectric metamaterials (see, e.g., [40–42]) have indicated the possible onset of intriguing “boundary” effects in critical parameter regimes, it appears suggestive to explore similar effects in our temporal scenario here. To this aim, we derive the conditions at the boundary of a nonlocal temporal metamaterial. More precisely, we consider a spatially homogeneous, unbounded, temporal metamaterial exhibiting a temporal boundary at a given time instant t = t ₀, at which the relative permittivity abruptly changes from a constant value ɛ ₁ to a time-varying periodic function ɛ(t), viz.,

Obviously, since the medium response to modulation cannot be infinitely fast, the assumption of abrupt temporal transitions is highly idealized, and a finite rise/fall time should be considered. However, for fall/rise times much shorter than the period of the incident wave, our approach remains valid. Within this framework, it is also worth pointing out that, in our numerical simulations (see the Methods Section 4.3), we actually assume smooth transitions with sufficiently short rise/fall times in order to favor convergence. As for the canonical temporal boundary [43], where the dielectric permittivity exhibits a temporal transition between two constant values, the microscopic electric induction D and magnetic induction B remain continuous across the boundary. Here, we consider plane waves, D = 2 Re[d(k, t)e^ik⋅r ] and B = 2 Re[b(k, t)e^ik⋅r ], experiencing the temporal boundary described by Eq. (11). By enforcing the standard boundary conditions [43], and taking into account the second of Eq. (3) and the microscopic Maxwell’s curl equation ∇ × B = μ 0 ∂ D ∂ t , we obtain

where d ± = d ̄ k , t 0 ± and b ± = b ̄ k , t 0 ± , with

From Eq. (12), it is evident that, in the limit k/K → 0, the inductions d and b are continuous, and our approach reproduces the known results for the canonical temporal boundary [43]. Generally, the homogenization process can hide the asymmetry of the system. In the spatial domain, the EM chirality is washed out in the conventional homogenized response of a composite metamaterial [6, 42]. Here, the system exhibits a time-reversal symmetry breaking due to the permittivity modulation in time and, while this effect is lost in the effective “bulk” response, it is restored in the boundary conditions expressed by Eq. (12). This is somehow analogous to what is observed in spatial multilayered metamaterials, which exhibit chiral boundary effects attributable to the parity symmetry breaking [42]. The abrupt switching of the permittivity breaks the time-reversal symmetry. On the other hand, by comparison with the conventional boundary conditions (i.e., continuity of d and b at the temporal boundary), we highlight that, in Eq. (12), the term proportional to β ₀ breaks explicitly the time-reversal symmetry. The nonlocality preserves the time-reversal asymmetry of the “microscopic” permittivity modulation, and this peculiar wave-matter interaction endows us with novel degrees of freedom for manipulating the wave propagation.

As an example, we consider once again the temporal metamaterial described by Eq. (10), with ε ̄ = 5 . Figure 3 shows the nonlocal effective parameters α ₀ and β ₀ as a function of the phase parameter ϕ, for Δ = 0.5, 0.65, 0.8. It is evident that the EM boundary response strongly depends on ϕ, which thus constitutes an additional “knob” for tailoring the wave propagation. As a general result, it is worth noting that α ₀ = 0 for ϕ = π/2, 3π/2, whereas β ₀ = 0 for ϕ = 0, π. By exploiting Eq. (13), one can prove that the parameters α ₀ and β ₀ vanish when the dielectric function exhibits even [ɛ(t − t ₀) = ɛ(t ₀ − t)] or odd [ɛ(t − t ₀) = −ɛ(t ₀ − t)] parity, respectively. Similar and related effects have been investigated in multilayered dielectric metamaterials [42], and topological photonic crystals [44]. To better understand the impact of the novel nonlocal terms in the boundary conditions expressed by Eq. (12), Figure 4 shows the two nonlocal contributions

as a function of k/K and ϕ, for Δ = 0.5 [panels (A) and (B)], Δ = 0.65 [panels (C) and (D)], and Δ = 0.8 [panels (E) and (F)]; as in previous examples, the parameters are only shown within the region where the nonlocal EMT holds (i.e., the deviation is less than 10% with respect to full-wave theory). Recalling the results in Figure 1, the larger the modulation depth Δ, the smaller the wavenumber region where the nonlocal EMT is in agreement with the full-wave theory. The results in Figure 4 indicate that the impact of nonlocality at the temporal boundary increases for deeper modulations. In particular, for Δ = 0.8, α(k) and β(k) reach the values 0.21 and 0.26, respectively, for k/K ≃ 0.74.

Figure 3:

Nonlocal temporal boundary. (A) and (B) Nonlocal effective boundary parameters α ₀ and β ₀, respectively, as a function of the modulation phase ϕ. Here, t ₀ = 0, and a metamaterial with the permittivity profile in Eq. (10) is considered with ε ̄ = 5 .

Figure 4:

Nonlocal contributions in Eq. (14) as a function of k/K and ϕ, for Δ = 0.5 ((A) and (B)), Δ = 0.65 ((C) and (D)), and Δ = 0.8 ((E) and (F)). Here, t ₀ = 0 and a temporal metamaterial with the permittivity profile in Eq. (10) is considered, with ε ̄ = 5 . The parameters are only shown within the region where the nonlocal EMT holds (≲10% error with respect to full-wave theory).

2.3 Reflection and transmission at a nonlocal temporal boundary

Let us consider a plane wave propagating in a spatially homogeneous unbounded medium with a relative permittivity as described by Eq. (11). Without loss of generality, we assume propagation along the positive z direction, and an x-polarized electric induction, i.e. D ( t , z ) = 2 Re [ d ( k , t ) e i k z ] e ̂ x . Following the nonlocal EMT model developed in Section 2.1, the average EM fields are governed by Eq. (5), where the effective wave-matter coupling is given by

for t < t ₀, and by Eq. (7) for t > t ₀. As a consequence, by neglecting the fast components, the electric induction d(k, t) is merely given by

where ω 1 ( k ) = c k / ε 1 , ω ( k ) = c k / ε eff μ eff ( k ) , and d _in, d _r, d _t are the incident, reflected (backward) and transmitted (forward) amplitudes, respectively. By means of the third of Eqs. (5) and (7), we can derive the magnetic field associated to the electric induction in Eq. (16). Then, by enforcing the temporal boundary conditions in Eq. (12), after some straightforward algebra, we obtain the temporal transmission t d = d t d in e − i [ ω ( k ) − ω 1 ( k ) ] t 0 and reflection r d = d r d in e i [ ω ( k ) + ω 1 ( k ) ] t 0 coefficients for the electric induction

where n 1 = ε 1 and n eff ( k ) = ε eff μ eff ( k ) . In the weak-dispersion regime (k/K ≪ 1), the magnetic effect is negligible [i.e., ω ( k ) ≃ c k / ε eff ] and the temporal transmission and reflection coefficients can be approximated as first-order Maclaurin series in k, viz.,

From Eq. (18), by taking into account Eqs. (15) and (7), and maintaining terms up to the first order in k/K, we obtain the corresponding expressions for the electric field

where the k-dependent terms account for the (weak) dispersion. By neglecting these nonlocal terms, we recover the well-known analytical expressions derived in previous studies on conventional (local) temporal boundaries [43].

2.4 Harnessing the temporal-boundary nonlocality

As for the spatial case [10], it is insightful to explore to what extent nonlocality in temporal metamaterials can be harnessed to attain some elementary pulse-shaping or analog-computing functionalities. By inspecting Eqs. (18) and (19), it is apparent that, in the weak-dispersion regime, backward and forward fields at a nonlocal temporal boundary are dominated by a linear combination of a local (constant) term and a nonlocal correction that is proportional to k, and hence corresponds to a first derivative. Remarkably, in the temporal reflection coefficients r _d(k) and r _e(k), this nonlocal correction can be made dominant by enforcing the impedance-matching condition ɛ _eff = ɛ ₁, so that the local terms vanish. Thus, the backward waveform is proportional to the first derivative of the incident one.

Via Fourier transform, a generic wavepacket can be synthesized as the superposition of the time-harmonic modes in Eq. (16), i.e., D ( z , t ) = ∫ − ∞ + ∞ d ( k , t ) e i k z d k . Then, by recalling Eq. (17), we can calculate the forward and backward wavepacket signals. As a representative example, we consider a Gaussian pulse interacting with the time-varying medium described by Eq. (11), where ɛ is given by Eq. (10), with ϕ = 3π/2. For t < t ₀, we assume that the incident pulse propagates in a dielectric medium (ɛ ₁ = 3) with the following profile

with v 1 = c / ε 1 , σ _t = 4τ, t ₀ = 10τ, and D ₀ being a constant. After the time instant t = t ₀, the pulse is partially reflected by the temporal boundary, as illustrated in Figure 5 (where the dielectric modulation depth is Δ = 0.8). It is evident that, for impedance-mismatched configurations, the backward pulses exhibit the same shape as the incident one [see panels (A) and (C), where ε ̄ = 4,6 , corresponding to ɛ _eff = 2.4, 3.6, respectively]. On the contrary, when the impedance is matched [see panel (B), where ε ̄ = 5 , corresponding to ɛ _eff = 3], the backward pulse reproduces the spatial first derivative of the incident one, as predicted by the second of Eq. (18). For completeness, Figure 6 shows the space-time maps for the three cases considered: ɛ _eff = 2.4 [panels (A)–(C)], ɛ _eff = 3 [panels (D)–(F)], and ɛ _eff = 3.6 [panels (G)–(I)].

Figure 5:

Spatial profiles of reflected (backward) pulses D _r (z, t ₁) (t ₁ ≃ 20τ) in the configuration described by Eqs. (11) and (10), with ɛ ₁ = 3, t ₀ = 10τ, ϕ = 3π/2, Δ = 0.8, and ε ̄ = 4,5,6 shown in panels (A)–(C), respectively. The blue curve represents the incident pulse profile D _in(z, t = 0), and the backward pulse D _r (z, t ₁) (orange-dashed) is superposed to the first spatial derivative of D _in(z, t = 0) (red).

Figure 6:

As in Figure 5, but space-time maps of the electric induction D and the corresponding forward (D _t ) and backward (D _r ) components for ɛ _eff = 2.4 ((A)–(C)), ɛ _eff = 3 ((D)–(F)), and ɛ _eff = 3.6 ((G)–(I)).

Finally, we validate these predictions from the nonlocal EMT model against full-wave numerical simulations (see the Methods Section 4.3 for details). Here, we consider the same EM configuration analyzed above with ε ̄ = 5 , and Δ = 0.8, so that the impedance matching condition ɛ _eff = ɛ ₁ = 3 is satisfied.

Figure 7 shows the results on three cases where ϕ = 3π/2 ((A)–(C)), ϕ = 0.67π ((D)–(F)), and ϕ = π ((G)–(I)). Specifically, we compare the electric induction and electric field distributions predicted by the conventional (local), nonlocal EMT, and full-wave simulations for both cases. Panels (B), (E), and (H) show the backward pulse profiles D _r/D ₀ at z = −30cτ for the configurations in panels (A), (D) and (G), respectively. At the impedance matching condition, the local EMT predicts zero temporal reflection, and it does not properly describe the backward wave dynamics. For ϕ = 3π/2 and ϕ = 0.67π, our proposed nonlocal EMT is in very good agreement with the full-wave simulations. For ϕ = π, the local and nonlocal EMT predict zero temporal reflection, whereas full-wave simulations yield a very small backward reflection signal (about an order of magnitude weaker than the previous cases). In this latter case, the parameter β ₀ ruling the nonlocal effect vanishes, and the temporal reflection is negligible [see Figure 3 along with the second of Eq. (18)]. Also shown in panels (C), (F) and (I) are the corresponding profiles of the normalized electric fields ɛ ₀ E _r/D ₀. In panels (C), (F) and (I), we observe that the full-wave predictions exhibit fast modulations due to the temporal modulation of the permittivity, whereas the nonlocal EMT prediction obtained from Eq. (19) is only representative of the slow component. Clearly, Eq. (19) describe the temporal transmission and reflection coefficients of the average component of the electric field, and all fast scales are not considered.

Figure 7:

Comparison between the backward wavepacket profiles predicted by full-wave simulations, local and non-local EMTs. The considered temporal dielectric profiles ensuring the impedance-matching conditions [i.e., the configuration described by Eqs. (11) and (10), with ɛ ₁ = 3, t ₀ = 10τ, Δ = 0.8, and ε ̄ = 5 ] is plotted with ϕ = 3π/2 (panel (A)), ϕ = 0.67π (panel (D)), and ϕ = π (panel (G)). Panels (B), (E), and (H) show the corresponding normalized electric inductions D _r/D ₀ for the backward pulses. Panels (C), (F), and (I) show the corresponding normalized electric fields ɛ ₀ E _r/D ₀. All temporal profiles are evaluated at z = −30cτ.

Overall, the above results confirm the validity of the nonlocal EMT predictions to harness the temporal-boundary nonlocality. In the configurations with ϕ = 3π/2 and ϕ = 0.67π, the nonlocal contributions, appearing in Eqs. (18) and (19), are significant (β ₀ = 0.31, −0.35 for ϕ = 3π/2, 0.67π, respectively), so that these results also confirm the possibility to attain the spatial first derivative of the incident wavepacket. Remarkably, by tailoring the modulation phase ϕ, the amplitude and the symmetry of the backward pulse can be suitably tuned. Comparing the field profiles in panels ((B) and (C)) with those in panels ((C) and (F)), it is evident that the pulse reverses its temporal symmetry profile by switching ϕ from 3π/2 to 0.67π. Quite interestingly, for ϕ = 0.67π, the relative-permittivity function ɛ _b (t) is continuous, thereby indicating that the predicted nonlocal effects do not depend critically on abrupt temporal changes.

In principle, by exploiting more complex modulation schemes with additional degrees of freedom (e.g., temporal dielectric structures with multiple harmonics), it could be possible to tailor the parameters so as to perform higher-derivative orders (and their linear combinations).

4.1 Nonlocal magnetism

Recalling that the effective Maxwell’s equations are invariant with respect to the Serdyukov–Fedorov transformation [50] D ′ = D ̄ − ∇ × Q and H ′ = B ̄ / μ 0 − ∂ t Q , after setting Q = ∂ _t M, we obtain the equivalent effective constitutive relationships

where

E ′ = E ̄ , and B ′ = B ̄ . In Eq. (21b), the vector M plays the role of an effective magnetic polarization. Therefore, in an isotropic and spatially homogeneous temporal metamaterial, the spatial dispersion up to the second-order (i.e., up to η ²) is fully equivalent to optical magnetism with the magnetic polarization given by Eq. (22). Considering the propagation of a monochromatic plane wave, we obtain that the EM fields experience a nonlocal magnetic response described by the effective relative magnetic permeability given in Eq. (9).

4.2 Rigorous dispersion relationship in time-periodic varying media

Following the rigorous approach in [51], we focus on the wave equation describing the electric field dynamics in a time-periodic varying medium. The propagation of a plane wave E ( z , t ) = 2 R e e ( t ) e i k z e ̂ x is described by the equation

where ɛ(t) = ɛ(t + τ). Since the permittivity is periodic in time, this equation admits Bloch-type modes e ( t ) = e ̃ ( ω , t ) e − i ω t , where e ̃ ( ω , t ) is a periodic function of period τ. By expanding in Fourier series the observables, Eq. (23) becomes

where c ̃ n , e ̃ n are the Fourier coefficients for the relative permittivity and electric field, respectively, and δ _m,n is the standard Kronecker-delta tensor (m, n = 0, ±1, ±2, … ). Eq. (24) is a set of linear equations that exhibits a nontrivial solution only if the associated determinant vanishes. By suitably truncating the Fourier series expansions, we numerically obtain the rigorous dispersion relationship ω(k) for waves propagating in a periodic temporal medium.

4.3 Full-wave simulations

We consider an arbitrary wavepacket propagating in a spatially homogeneous unbounded, time-varying metamaterial, with the temporal boundary described by Eq. (11). The wavepacket electric induction can be synthesized via Fourier transform as

Then, for each value of the wavenumber k, we define the auxiliary functions

where d k , t and b k , t are the plane-wave spectra of the electric and magnetic inductions, respectively, at a given time t, with Z ₀ and D ₀ denoting the vacuum intrinsic impedance and a dimensional normalization constant, respectively. From Maxwell’s curl equations, we derive a pair of coupled ordinary differential equations, namely

with initial conditions

where d _in(k, 0) is the plane wave spectrum of the incident electric induction field at t = 0. Next, we solve numerically Eq. (27) by means of the NDSolve routine available in Mathematica™ [52]. This routine provides the numerical solution of generic systems of ordinary differential equations, via a broad arsenal of methods (including Runge–Kutta, predictor-corrector, implicit backward differentiation) that can be tailored adaptively to the specific scenario of interest and, in principle, it can automatically handle discontinuities in the equations [52]. In our implementation, we utilize default settings and parameters. Moreover, in order to favor numerical convergence, we implement the abrupt permittivity changes via an analytical, smooth unit-step function U s t = tanh t / T s + 1 / 2 , where T _s = τ/100.

Once a numerical solution is available for Eq. (27), the electric induction is synthesized via Eq. (25) (where d = D ₀ u ₁), whereas the corresponding electric field can be readily obtained via division by ε b t . In our numerical implementation, this synthesis is implemented via fast-Fourier-transform by means of the Four er routine available in Mathematica™ [52].