Home Physical Sciences Space–time graded-index interfaces and related chirping
Article Open Access

Space–time graded-index interfaces and related chirping

  • Zhiyu Li ORCID logo EMAIL logo , Xikui Ma ORCID logo , Klaas De Kinder ORCID logo , Amir Bahrami ORCID logo and Christophe Caloz ORCID logo
Published/Copyright: September 5, 2025

Abstract

Space–time modulated systems have recently emerged as a powerful platform for dynamic electromagnetic processing in both space and time. Most of the related research so far has assumed abrupt parameter profiles. This paper extends the field to generalized graded-index (GRIN) interfaces, which are both more practical than ideal profiles and offer new avenues for wave manipulations. It presents an exact solution for wave propagation across arbitrary space–time modulated GRIN interfaces and describes versatile chirping effects. The solution is based on a generalization of the impulse response method from linear time-invariant to linear space–time-varying systems. The proposed framework shows that space–time GRIN systems represent a novel approach for generating a new form of chirping that is not inherently based on dispersion, with promising applications in pulse shaping and signal processing.

1 Introduction

Space–time modulation systems – media whose parameters vary dynamically in both space and time under the influence of an external traveling-wave modulation – have recently enabled a range of novel applications, including magnet-free nonreciprocity [1], frequency transitioning [2], parametric amplification [3], [4], and the breaking of fundamental bounds [5].

Space–time interfaces are the fundamental building blocks – or “meta-atoms” – of space–time modulation systems. Their understanding and control are, therefore, essential. These interfaces can be of two distinct types, as illustrated in Figure 1. The first is the step-index interfaces, where the transition width is much smaller than the wavelength, as shown in Figure 1(a). Such interfaces have been extensively studied in the literature [6], [7] and are known to induce a uniform frequency shifting of the incoming wave. The second type is the graded-index (GRIN) interfaces. These interfaces, shown in Figure 1(b), may be considered more practical, as physical modulations necessarily occur gradually at the microscopic level [8], [9], [10], and can also support a broader range of wave transformations. In particular, they induce nonuniform frequency transitions and produce wave chirping. While some studies have addressed the limiting case of pure-time GRIN interfaces [11], [12], [13], [14], general space–time GRIN interfaces [15] remain an essentially unexplored topic.

Figure 1: 
Space–time modulation interfaces, with basic structures (top) and space–time diagrams in the harmonic-wave regime (bottom). (a) Step-index modulation, where the refractive index changes abruptly from n
1 to n
2 at an interface initially located at 




z


b


0




${z}_{\text{b}}^{0}$



 and moving with a constant velocity v. (b) Graded-index (GRIN) modulation, where the refractive index traditions smoothly from n
1 to n
2 via n
G over the interval 






z


b
1


0


,


z


b
2


0






$\left[{z}_{\text{b}1}^{0},{z}_{\text{b}2}^{0}\right]$



 with width 


D
=


z


b
2


0


−


z


b
1


0




$D={z}_{\text{b}2}^{0}-{z}_{\text{b}1}^{0}$



. The subscripts G and b denote the GRIN layer and the boundary, respectively.
Figure 1:

Space–time modulation interfaces, with basic structures (top) and space–time diagrams in the harmonic-wave regime (bottom). (a) Step-index modulation, where the refractive index changes abruptly from n 1 to n 2 at an interface initially located at z b 0 and moving with a constant velocity v. (b) Graded-index (GRIN) modulation, where the refractive index traditions smoothly from n 1 to n 2 via n G over the interval z b 1 0 , z b 2 0 with width D = z b 2 0 z b 1 0 . The subscripts G and b denote the GRIN layer and the boundary, respectively.

This paper investigates the wave transformations induced by space–time GRIN interfaces, introducing a generalization of the impulse response method from linear time-invariant (LTI) to linear (space-) time-varying (LTV) systems. Unlike the Wentzel–Kramers–Brillouin (WKB) [11], [12], [13] or transfer-matrix method (TMM) [15] commonly used in prior works, this approach provides an exact solution. We derive closed-form solutions for both the analysis and synthesis problems and describe the corresponding chirping effects. All the results are validated through full-wave simulations.

For simplicity, we restrict our attention to systems that are 1 + 1D, with dimensions z and t, involving GRINs of uniform velocity, i.e., v = v z ̂ = c o n s t . , boundary and intrinsic impedance matching, i.e., η i = μ i / ϵ i = c o n s t . for i = 1, G, 2, and no dispersion, i.e., n i n i (ω) for i = 1, G, 2.

2 Generalized impulse response method

In this section, we shall derive the generalized impulse for the GRIN system in Figure 1(b), which is generically represented in Figure 2.

Figure 2: 
Generalized impulse response method applied to the space–time GRIN interface system in Figure 1(b). (a) Responses in the frequency and time domains. (b) Response of an arbitrary space–time GRIN modulation system. The top panel shows the space–time diagram with the trajectory of the impulse δ(t − t′) and its output 


δ


t
−


t


n



(

z
,


t


′



)





$\delta \left[t-{t}_{n}\left(z,{t}^{\prime }\right)\right]$



, where p
1,2 are the space–time intersection points of the impulse with the two boundaries of the GRIN layer, initially located at 




z


b
1,2


0




${z}_{\text{b}1,2}^{0}$



. The bottom panel shows the initial refractive index profile n(z, t = 0), which moves at a uniform velocity v.
Figure 2:

Generalized impulse response method applied to the space–time GRIN interface system in Figure 1(b). (a) Responses in the frequency and time domains. (b) Response of an arbitrary space–time GRIN modulation system. The top panel shows the space–time diagram with the trajectory of the impulse δ(tt′) and its output δ t t n ( z , t ) , where p 1,2 are the space–time intersection points of the impulse with the two boundaries of the GRIN layer, initially located at z b 1,2 0 . The bottom panel shows the initial refractive index profile n(z, t = 0), which moves at a uniform velocity v.

In the frequency domain, the one-dimensional Helmholtz equation may be expressed as

(1) 2 E ̃ z 2 + k 2 ( ω ) E ̃ = 0 ,

where k ( ω ) = ω μ ϵ = ω n / c is the wavenumber. The general solution to Eq. (1) can be written as

(2a) E ̃ ( z , ω ) = H ( z , ω ) E ̃ ( 0 , ω ) ,

where

(2b) H ( z , ω ) = e i k ( ω ) z

represents the frequency response of the system.

Applying the temporal inverse Fourier transform to Eq. (2), as suggested in Figure 2(a), we obtain the corresponding time-domain relation

(3) E ( z , t ) = h ( z , t , t ) E ( 0 , t ) d t ,

where h(z, t, t′) is the impulse response [16], representing the system’s response to the impulse E(0, t) = δ(tt′).[1]

In LTI systems, the impulse response depends solely on the time difference between the output and input signals, following the time-shift invariance property h(z, t, t′) = h(z, tt′) [18]. Xiao et al. extended the impulse response method from LTI to purely time-varying systems [19], where the lack of time invariance – or nonstationarity – leads to an impulse response that independently depends on t and t′, i.e.,

(4) h ( z , t , t ) h ( z , t t ) .

In space–time varying systems, h(z, t, t′) also involves the modulation velocity v, which further increases the complexity of the analysis. We now introduce the generalization of the impulse response method to space–time GRIN systems, explicitly incorporating v, with the aid of Figure 2(b), where a space–time diagram illustrates the impulse trajectory.

As shown in the top panel of Figure 2(b), the input impulse experiences a propagation delay as it propagates through the system, resulting in the impulse response

(5) h ( z , t , t ) = δ t t n ( z , t ) ,

where t n (z, t′) is the arrival time of the impulse at the position z assuming an input time of t′ in the GRIN system with space–time varying refractive index n [brown line in the top panel of Figure 2(b)].[2] The function t n (z, t′) is determined by the wave trajectory equation

(6) c d t n d z = n ( z , t n ) ,

subject to the boundary conditions in different regions

(7) t n ( 0 ) = t and t n ( z b 1,2 ) = t b 1,2 ,

where z b1,2 and t b1,2 are the coordinates of the intersection points p 1,2 of the impulse with the first and second modulation interfaces, respectively [Figure 2(b)].

To derive the electromagnetic field in each region of the space–time GRIN system, we will apply the following four steps to each of the three regions in Figure 2(b): (i) determine the trajectory equation and corresponding boundary condition for the concerned region; (ii) solve the trajectory equation for the arrival time t n (z, t′); (iii) substitute the resulting t n (z) into Eq. (5) to obtain the impulse response h(z, t, t′); and (iv) substitute that response into Eq. (3) to evaluate the output field E(z, t).

3 Field solutions

In this section, we shall derive the general field solutions for the GRIN system in Figure 2(b) using the generalized impulse response derived in Section 2.

The refractive index profile along the impulse trajectory is given by

(8) n ( z , t ) = n 1 , 0 < z < z b 1 0 + v t , n G ( z v t ) , z b 1 0 + v t < z < z b 2 0 + v t , n 2 , z > z b 2 0 + v t ,

where n 1 and n 2 are the (constant) refractive indices of media 1 and 2, respectively, and

(9) z b 2 0 = z b 1 0 + D .

The system can then be divided into the three corresponding regions, which we will address one by one.

3.1 First-medium (n 1) region

Substituting n(z, t n ) = n 1 into Eq. (6), we obtain the wave trajectory equation

(10a) c d t n d z = n 1 ,

where the boundary condition is

(10b) t n ( 0 ) = t .

Solving Eq. (10) for t n (z, t′) yields then the impulse trajectory function

(11) t n ( z , t ) = t + n 1 c z ,

whose insertion into Eq. (5) provides the impulse response

(12) h ( z , t , t ) = δ t t n 1 z c .

Note that the argument of the impulse function in this relation represents a traveling-wave, due to the uniform nature of the propagation medium. Finally, substituting Eq. (12) into (3) and solving for the output field E(z, t), we get the field,

(13) E 1 ( z , t ) = E 0 , t n 1 z c ,

where E(⋅, ⋅) is an arbitrary field function (e.g., harmonic plane wave or Gaussian pulse) of the space (first entry) and time (second entry) variables.[3]

3.2 GRIN (n G) region

Substituting now n(z, t n ) = n G(zvt n ) into Eq. (6), we obtain the wave trajectory equation in the GRIN region,

(14) c d t n d z = n G ( z v t n ) .

The corresponding space–time boundary condition with medium 1 corresponds to the intersection point p 1 [see Figure 2(b)], whose coordinates are related as

(15a) t b 1 = t + n 1 c z b 1

and

(15b) z b 1 = z b 1 0 + v t b 1 .

Solving Eq. (15) for z b1 and t b1, we find the coordinates of p 1 to be

(16) z b 1 = v t + z b 1 0 1 n 1 v / c and t b 1 = t + n 1 z b 1 0 / c 1 n 1 v / c .

To solve Eq. (14) with the related boundary condition, we let ξ = zvt n , which implies d ξ d z = 1 v d t n d z , or

(17) d t n d z = 1 v 1 d ξ d z .

Substituting now Eq. (17) into (14), and separating the ξ and z terms, leads to the differential equation

(18) d ξ 1 n G ( ξ ) v / c = d z ,

which integrates to

(19) 1 1 n G ( ξ ) v / c d ξ = z + C G ,

where C G is an integration constant. For generalization to arbitrary GRIN profiles, we define the left-hand side term of this relation as the function

(20) F ( ξ ) = 1 1 n G ( ξ ) v / c d ξ ,

which allows to express Eq. (19) in the compact form

(21) F ( ξ ) = z + C G .

To determine the integration constant C G in this relation, we first apply the boundary condition t n (z b1) = t b1 (point p 1), which yields

(22) F ( z b 1 v t b 1 ) = z b 1 + C G .

Next substituting Eq. (16) into this relation and solve for C G, we obtain

(23) C G ( t ) = v 1 n 1 v / c t z b 1 0 1 n 1 v / c + F z b 1 0 .

Finally, substituting Eq. (23) and the relation ξ = zvt n into Eq. (21), and solving for t n , we find

(24) t n ( z , t ) = 1 v F 1 z + C G ( t ) + z v ,

where F −1(⋅) is the inverse function of F(⋅).

Substituting Eq. (24) into (5), we obtain now the impulse response

(25) h ( z , t , t ) = δ t + 1 v F 1 z + C G ( t ) z v ,

where C G(t′) is given in Eq. (23) and F(⋅) is defined in Eq. (20).

Finally, substituting Eq. (25) into (3), and solving for the output field E(z, t), yields (see Appendix A)

(26) E G ( z , t ) = 1 n 1 v / c 1 n G ( z v t ) v / c × E 0 , 1 n 1 v / c v z b 1 0 + v t z 1 1 n G ( z v t ) v / c d z n 1 z c + z z b 1 0 v ,

where E(⋅, ⋅) is the same field function of space and time as in Eq. (13).

3.3 Second-medium (n 2) region

Substituting now n(z, t n ) = n 2 into Eq. (6), we obtain

(27) c d t n d z = n 2 ,

where the space–time boundary condition with the GRIN medium corresponding to the intersection point p 2, whose coordinates are related as

(28a) t b 2 = 1 v F 1 z b 2 + C G ( t ) + z b 2 v

and

(28b) z b 2 = z b 2 0 + v t b 2 ,

and found by solving Eq. (28) for z b2 and t b2 as

(29a) z b 2 = F z b 2 0 C G ( t )

and

(29b) t b 2 = F z b 2 0 C G ( t ) z b 2 0 v .

Solving next Eq. (27) for t n (z, t′), we obtain the impulse trajectory function

(30a) t n ( z , t ) = n 2 c z + C 2 ,

where C 2 is a new integration constant, which is obtained by applying the boundary condition t n (z b2) = t b2 [Eq. (29), point p 2] to Eq. (30a) and solving the resulting expression for C 2 as

(30b) C 2 ( t ) = 1 n 2 v / c 1 n 1 v / c t + 1 n 2 v / c v F z b 1 0 + D F z b 1 0 n 2 n 1 1 n 1 v / c z b 1 0 c D v .

Finally, substituting Eq. (30) into (5), we get the impulse response

(31) h ( z , t , t ) = δ t n 2 c z C 2 ( t ) ,

where C 2(t′) is given in Eq. (30b).

Substituting Eq. (31) into (3) and solving for the output field E(z, t) (see Appendix B), we obtain the field[4]

(32) E 2 ( z , t ) = 1 n 1 v / c 1 n 2 v / c × E 0 , 1 n 1 v / c 1 n 2 v / c t n 2 c z + n 2 n 1 1 n 1 v / c z b 1 0 c + D v 1 n 2 v / c v z b 1 0 + v t z b 1 0 + D + v t 1 1 n G ( z v t ) v / c d z .

Equations (13), (26), and (32) represent the key results of this paper. They accommodate arbitrary field waveforms E(0, t) and arbitrary GRIN profiles n G(zvt).

4 Chirping physics

Due to the space- and time-varying properties of the GRIN medium, the wave behavior in our system is fairly complex, as evident in Eq. (26). In this section, we show that this system produces a new type of wave chirping and compare it with other chirping mechanisms.

4.1 Chirping in the GRIN layer

We consider a time-harmonic incident field, which reads at z = 0,

(33) E ( 0 , t ) = e i ω 0 t ,

where ω 0 is assumed to be constant. Other types of fields can be treated similarly. Substituting Eq. (33) into (26) – specifically, inserting the content of the second slot of E(⋅, ⋅) in Eq. (26) into the relation ϕ G = −ω 0 t(z) – we obtain the wave phase in the GRIN layer as

(34) ϕ G = ω 0 1 n 1 v / c v z b 1 0 + v t z 1 1 n G ( z v t ) v / c d z + ω 0 n 1 z c z z b 1 0 v .

The related instantaneous frequency at a given position in the GRIN layer is obtained by differentiating the phase in Eq. (34) with respect to time, yielding

(35) ω G = ω 0 1 n 1 v / c v t z b 1 0 + v t z 1 1 n G ( z v t ) v / c d z = ω 0 1 n 1 v / c 1 n G ( z v t ) v / c ,

where we have used the Leibniz integral rule. Equation (35) reveals that the wave frequency in the GRIN layer varies with time, through the function n G(zvt) in the denominator, indicating a space–time chirping effect. This effect is fundamentally different from the group-velocity dispersion (GVD) chirping effect occurring in dispersive media [18], since no dispersion is present. It will be explained shortly.

Furthermore, the chirp parameter α, whose sign determines whether the field is up-chirping (α > 0) or down-chirping (α < 0), is obtained by time-differentiating the instantaneous frequency in Eq. (35), which gives

(36) α = ω G t = ω 0 ( 1 n 1 v / c ) c [ 1 n G ( z , t ) v / c ] 2 n G ( z , t ) z .

This result indicates that the sign of the chirp, α ≷ 0, depends on

(37) ( 1 n 1 v / c ) n G ( z , t ) / z 0 .

According to Eq. (37), the up- or down-chirping behavior of the system is governed by the velocity regime, which may be subluminal [v < c/max(n 1, n 2)] or superluminal [v > c/min(n 1, n 2)] [21]. To understand the chirping mechanism within the GRIN layer, we now focus on the subluminal case – the superluminal case can be analyzed analogously. In the subluminal regime, where v < c/n 1,[5] Eq. (37) simplifies to

(38) n G ( z , t ) / z 0 ,

indicating that a decreasing refractive index slope leads to an increase in the instantaneous frequency, i.e., up-chirping, and vice versa, i.e., down-chirping.

For simplicity, but without loss of generality, we consider the simplest GRIN profile, the linear profile.

(39) n G ( z v t ) = n 1 + n 2 n 1 D z v t z b 1 0 .

Figure 3 represents the corresponding trajectories for the wave crests incident at z = 0 at the times t i1,2,3,4,5, as derived from the analytical solutions in the different regions, given by Eqs. (13), (26), and (32). Figure 3(a) considers the case where the refractive index increases within the GRIN layer, i.e., n 2 > n 1, corresponding to the positive spatial gradient ∂n G(z, t)/∂z > 0, which leads to a down-chirping according to the condition in Eq. (38). For a fixed observation position z o within the GRIN region, the arrival times of the five incident crests are denoted as t o1,2,3,4,5. Due to the space–time variation of n G, each crest experiences a different local refractive index at the observation point, as indicated by the color-coded dots in the figure. The latest crest (purple dot) propagates fastest, as it has just entered the GRIN region and hence encounters the lowest refractive index. In contrast, the earliest crest (red dot) travels slowest, as it is near the exit of the GRIN region and hence sees the highest refractive index. As a result, the later, faster crest gradually catches up with the earlier, slower one, leading to an increasing temporal separation between adjacent crests. This corresponds to a decreasing frequency over time, i.e., down-chirping, consistent with the down-chirping condition in Eq. (38). A similar mechanism applies in the case of a decreasing refractive index, i.e., n 2 < n 1, where the crests compress in time, resulting in up-chirping [Eq. (38)], as illustrated in Figure 3(b).

Figure 3: 
Space–time diagrams of a linearly varying GRIN medium [Eq. (39)] in the subluminal velocity regime (v = 0.25c) used for the chirping analysis, with (a) a positive n
G slope, where n
1 = 1, n
2 = 2 and D = 1.2λ
0 with λ
0 = cT
0 being the free-space wavelength, and (b) a negative n
G slope, where n
1 = 2, n
2 = 1 and D = 1.2λ
0. The bottom panels show the corresponding initial refractive index profiles, n(z, t = 0).
Figure 3:

Space–time diagrams of a linearly varying GRIN medium [Eq. (39)] in the subluminal velocity regime (v = 0.25c) used for the chirping analysis, with (a) a positive n G slope, where n 1 = 1, n 2 = 2 and D = 1.2λ 0 with λ 0 = cT 0 being the free-space wavelength, and (b) a negative n G slope, where n 1 = 2, n 2 = 1 and D = 1.2λ 0. The bottom panels show the corresponding initial refractive index profiles, n(z, t = 0).

Note that the wave exiting the GRIN layer is no longer chirped. This because, at the exit points, all the crest have experienced the entire GRIN profile and find, therefore, themselves resynchronized to the velocity of the second medium.

The relation (35) offers a practical foundation for designing systems with a prescribed chirp profile, where the frequency varies according to a desired function f(t),

(40) ω G = f ( t ) .

To realize such a frequency evolution, one can tailor the refractive index profile of the GRIN layer, by substituting Eq. (35) into (40), and solving for n G, yielding

(41) n G = c v 1 ω 0 ( 1 n 1 v / c ) f ( t ) ,

which provides a closed-form expression for engineering a GRIN profile generating the desired chirp. As an example, for a linear chirp, where

(42) f ( t ) = a + b ( z v t )

with a and b being constants, the refractive index profile becomes

(43) n G = c v 1 ω 0 ( 1 n 1 v / c ) a + b ( z v t ) .

4.2 Comparison with other chirping mechanisms

In this section, we review the main chirping mechanisms and compare them with the GRIN-based mechanism described in the previous section.

Dispersive systems are the most common conventional means of generating chirping. In these systems, group-velocity dispersion (GVD) [18] causes a frequency-dependent group delay, resulting in a quadratic (or higher-order) spectral phase and corresponding variations in the instantaneous frequency across the pulse. Since this process alters only the spectral phase while leaving the spectral amplitude unchanged, the temporal waveform is scaled without spectral alteration, which limits the achievable chirp range for a given input bandwidth. Chirping can also arise in nonlinear systems [22]. A common example is self-phase modulation (SPM) [18], in which a high-intensity pulse propagating through a Kerr medium experiences an intensity-dependent refractive index change, imprinting a corresponding nonlinear phase on the pulse. The resulting time-varying phase produces instantaneous frequency shifts that evolve along the propagation distance. Because the effect scales nonlinearly with the optical intensity and interaction length, SPM-induced chirping offers limited independent control over the chirp profile and is less straightforward to tune than dispersive methods.

Recently, other chirping mechanisms have been explored in LTV systems. In these structures, chirping arises from various effects at a time-varying interface, which generate new frequency components and thereby alters the magnitude spectrum – in both shape and bandwidth – within the linear regime. For instance, Shlivinski and Hadad demonstrated transient chirped-like radiation from a lossy, dispersive time-varying slab into air, attributed to the conservation of longitudinal wavenumber at the temporal interface [23]. Another example involves accelerated space–time step-index interfaces [24], [25], [26], where nonuniform modulation velocities induce Doppler-based time-dependent frequency transitions and associated chirping effects.

The space–time GRIN-based chirping system presented in this paper represents a new class of LTV chirping mechanisms. In these systems, chirping arises from the space–time variation of medium properties within the modulation slab, which induces local wave velocity differences that modify the instantaneous frequency. Compared to the other two LTV structures, GRIN-based systems offer greater design flexibility for chirping and do not require complicated nonuniform modulation velocities.

5 Illustrative results

Figure 4 plots the electric field magnitudes across two space–time GRIN interfaces computed by Eqs. (13), (26), and (32). Figure 4(a) corresponds to a hyperbolic tangent GRIN interface profile. It may be observed that the field experiences a gradual down-chirping in the GRIN region, as expected from α < 0 [Eq. (38)], before reaching a steady frequency in the second medium. Figure 4(b) corresponds to a sinusoidal GRIN interface profile. In this case, the field undergoes nonmonotonic, twisted chirping within the GRIN layer with varying alternating chirping sign [Eq. (38)], and eventually recovers its original frequency after exiting the modulated region. In both cases, the closed-form field distributions in the figure have been validated against full-wave finite-difference time-domain (FDTD) simulations [15], [24] (see Appendix C).

Figure 4: 
Electric field magnitude |E

x
| across space–time GRIN interfaces, computed by Eqs. (13), (26), and (32), for the input pulse 


E

(

0
,
c
t

)

=


e


−



(

t
−
4


T


0



)



2


/
2


T


0


2






e


−
i


ω


0


t




$E\left(0,ct\right)={\mathrm{e}}^{-{\left(t-4{T}_{0}\right)}^{2}/2{T}_{0}^{2}}{\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}$



, interface velocity v = 0.2c, and different GRIN profiles, (a) a hyperbolic tangent profile, 




n


G



(

z
−
v
t

)

=


n


1


+

[


(



n


2


−


n


1



)

/
2

]


[

1
+
10
⁡
tanh


z
−
v
t
−


z


b
1


0


−
D
/
2


/
D

]



${n}_{\mathrm{G}}\left(z-vt\right)={n}_{1}+\left[\left({n}_{2}-{n}_{1}\right)/2\right]\left[1+10\mathrm{tanh}\left(z-vt-{z}_{\text{b}1}^{0}-D/2\right)/D\right]$



, with n
1 = 1.5, n
2 = 3 and D = 2λ
0, and (b) a sine profile, 




n


G



(

z
−
v
t

)

=


n


1


+

[

0.8
+


n


2


−


n


1



]

sin


2
π


z
−
v
t
−


z


b
1


0




/
D




${n}_{\mathrm{G}}\left(z-vt\right)={n}_{1}+\left[0.8+{n}_{2}-{n}_{1}\right]\mathrm{sin}\left[2\pi \left(z-vt-{z}_{\text{b}1}^{0}\right)/D\right]$



, with n
1 = n
2 = 2 and D = 2λ
0. The top panels show the space–time diagrams of the normalized electric field magnitude 



[

|




E

̄



x



(

z
,
c
t

)

|
=
|


E


x



(

z
,
c
t

)

|
/
max

(

|


E


x



(

0
,
c
t

)

|

)


]



$\left[\vert {\bar{E}}_{x}\left(z,ct\right)\vert =\vert {E}_{x}\left(z,ct\right)\vert /\mathrm{max}\left(\vert {E}_{x}\left(0,ct\right)\vert \right)\right]$



 under modulated Gaussian pulse excitation, where the white solid lines mark the two boundaries of the GRIN region. The middle panels show the refractive index profiles n(z, t) [Eq. (8)] at t = 0. The bottom panels show the normalized spectrograms, 


|




E

̄



x



(

t
,
ω

)



|


2




$\vert {\bar{E}}_{x}\left(t,\omega \right){\vert }^{2}$



, with the input pulse being replaced by the quasi-continuous wave 


E

(

0
,
t

)

=


e


−
i


ω


0


t


rect

(

t
/
τ

)



$E\left(0,t\right)={\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}\text{rect}\left(t/\tau \right)$



 (with τ = 30T
0), for easier visualization, across the GRIN layer at z = 4.5λ
0 between points Q
1 and Q
2. The dashed black line corresponds to the instantaneous frequency ω
G(t) given by Eq. (35).
Figure 4:

Electric field magnitude |E x | across space–time GRIN interfaces, computed by Eqs. (13), (26), and (32), for the input pulse E ( 0 , c t ) = e ( t 4 T 0 ) 2 / 2 T 0 2 e i ω 0 t , interface velocity v = 0.2c, and different GRIN profiles, (a) a hyperbolic tangent profile, n G ( z v t ) = n 1 + [ ( n 2 n 1 ) / 2 ] [ 1 + 10 tanh z v t z b 1 0 D / 2 / D ] , with n 1 = 1.5, n 2 = 3 and D = 2λ 0, and (b) a sine profile, n G ( z v t ) = n 1 + [ 0.8 + n 2 n 1 ] sin 2 π z v t z b 1 0 / D , with n 1 = n 2 = 2 and D = 2λ 0. The top panels show the space–time diagrams of the normalized electric field magnitude [ | E ̄ x ( z , c t ) | = | E x ( z , c t ) | / max ( | E x ( 0 , c t ) | ) ] under modulated Gaussian pulse excitation, where the white solid lines mark the two boundaries of the GRIN region. The middle panels show the refractive index profiles n(z, t) [Eq. (8)] at t = 0. The bottom panels show the normalized spectrograms, | E ̄ x ( t , ω ) | 2 , with the input pulse being replaced by the quasi-continuous wave E ( 0 , t ) = e i ω 0 t rect ( t / τ ) (with τ = 30T 0), for easier visualization, across the GRIN layer at z = 4.5λ 0 between points Q 1 and Q 2. The dashed black line corresponds to the instantaneous frequency ω G(t) given by Eq. (35).

Figure 5 presents two linear-chirping GRIN designs using Eq. (43). In Figure 5(a), the design is performed in a co-moving (up-chirping) subluminal regime with v = 0.3c, while in Figure 5(b), it corresponds to a contra-moving (down-chirping) superluminal regime with v = −0.85c. The spectrograms in the bottom panels, corresponding to the GRIN profiles in the top panels, precisely match the theoretical predictions (see Appendix C).

Figure 5: 
Design of a linear-chirping [Eq. (42)] GRIN system, corresponding to the refractive index profiles obtained from Eq. (43) and shown in the top panels. The input pulse is the same quasi-continuous wave as in the bottom panels of Figure 4. (a) Up-chirping system with n
1 = 2, n
2 = 1.58, v = 0.3c (subluminal regime), D = 2λ
0, a = 1 and b = −0.4. (b) Down-chirping system with n
1 = 1.2, n
2 = 2.22, v = −0.85c (superluminal regime), D = 5λ
0, a = 1 and b = −0.2. The top and bottom panels show the GRIN profiles at t = 0 and the corresponding spectrograms, 


|




E

̄



x



(

t
,
ω

)



|


2




$\vert {\bar{E}}_{x}\left(t,\omega \right){\vert }^{2}$



, respectively. The dashed black lines represent the target linear-chirp function f(t) [Eq. (42)].
Figure 5:

Design of a linear-chirping [Eq. (42)] GRIN system, corresponding to the refractive index profiles obtained from Eq. (43) and shown in the top panels. The input pulse is the same quasi-continuous wave as in the bottom panels of Figure 4. (a) Up-chirping system with n 1 = 2, n 2 = 1.58, v = 0.3c (subluminal regime), D = 2λ 0, a = 1 and b = −0.4. (b) Down-chirping system with n 1 = 1.2, n 2 = 2.22, v = −0.85c (superluminal regime), D = 5λ 0, a = 1 and b = −0.2. The top and bottom panels show the GRIN profiles at t = 0 and the corresponding spectrograms, | E ̄ x ( t , ω ) | 2 , respectively. The dashed black lines represent the target linear-chirp function f(t) [Eq. (42)].

6 Conclusion and discussion

We have presented an exact electromagnetic solution to the problem of wave propagation across arbitrary space–time modulated GRIN interfaces and a detailed description of the related chirping effects. This solution offers a novel approach to chirp generation, which can be realized using artificial transmission lines at microwave frequencies [27], [28], [29] and acoustic or optical wave-based modulation techniques at optical frequencies [8], [9], [10], [30].

In the microwave regime, such interfaces can be realized in the subluminal regime by injecting pump pulses into transmission lines loaded with nonlinear capacitive and inductive elements, such as varactors and ferrite cores [28]. The pump signal applied at the terminal propagates along the line, inducing a dynamic, intensity-dependent refractive index variation that forms a moving GRIN interface between regions with different electromagnetic properties. A linear probe signal interacting with this modulation interface undergoes the designed chirping effect. While achieving superluminal modulation velocities is impossible with such a pump-probe platform, it could be potentially realized using switched transmission lines composed of subwavelength units spaced by Δz, each loaded with a sequence of elements having different parameters controlled by switches to create an effective spatial gradient, n G(z). An external controller, such as a field-programmable gate array (FPGA), actuates these switches with a time interval Δt, enabling sequential time delays [27]. By adjusting the ratio Δzt, this spatial gradient can propagate at an effectively unlimited velocity v = Δzt ∈ (0, ∞), forming an effective moving GRIN interface, n G(zvt).

In the optical regime, similar pump–probe setup can be achieved by launching a strong pump obliquely at an angle θ relative to the probe, producing a modulation velocity v = c/sinθ ∈ (0, ∞) that spans both subluminal and superluminal regimes [30]. Dynamic permittivity modulation can be realized via surface or bulk acoustic waves in piezoelectric crystals, ultrafast laser pulses in semiconductor slabs, or epsilon-near-zero (ENZ) materials for higher refractive index contrast [10]. Dynamic control of permeability can be achieved through the magneto-optical response of gyromagnetic materials. Experimental challenges in realizing arbitrary GRIN profiles may be addressed by employing arbitrary waveform generators, enabling synthesis of the desired modulation profiles with high precision. This paper advances the modeling of space–time dispersive systems and introduces a new paradigm for linear pulse shaping. The approach can also be integrated with various platforms for dispersion compensation, enabling enhanced control over wave propagation in dynamic media.


Corresponding author: Zhiyu Li, Xi’an Jiaotong University, Xi’an, China, E-mail: 

Award Identifier / Grant number: G0B0623N

  1. Research funding: This work was supported by Grant #G0B0623N from the Fonds voor Wetenschappelijk Onderzoek (FWO).

  2. Author contributions: All authors have accepted responsibility for the entire content of this manuscript and consented to its submission to the journal, reviewed all the results, and approved the final version of the manuscript. ZL performed the bulk of the work. XM provided supervision and guidance. KDK contributed to the development of Section 4. AB assisted in simplifying the theoretical formulation in Section 3. CC supervised the research and contributed to the manuscript preparation.

  3. Conflict of interest: Authors state no conflict of interest.

  4. Data availability: The datasets generated during and/or analyzed during the current study are available from the corresponding author on a reasonable request.

Appendix A Derivation of Eq. (26)

In this appendix, we provide a detailed derivation of Eq. (26) given in Section 3.2. Substituting Eq. (25) into (3), we get

(44) E G ( z , t ) = δ t + 1 v F 1 z + C G ( t ) z v E ( 0 , t ) d t .

To simplify the notation, we define the argument of the delta function in Eq. (44) as

(45) g G ( t ) = t + 1 v F 1 z + C G ( t ) z v .

Moreover, we use the following property of the delta function [31]

(46) δ g G ( t ) f ( 0 , t ) d t = f ( 0 , t G ) | g G t G | ,

where g G t G = 0 and g G ( ) = d g G ( ) / d t . Therefore, we first solve for the root t G of g G(t′) = 0, i.e.,

(47) t + 1 v F 1 z + C G t G z v = 0 ,

which gives

(48) C G t G = F ( ξ ) z .

Substituting Eq. (23) into the left-hand side of Eq. (48) and solving for t G , we find

(49) t G = 1 n 1 v / c v F ( ξ ) F z b 1 0 n 1 z c + z z b 1 0 v .

Applying Eq. (20) and performing the change of variables z′ = ξ + vt, the difference in the square brackets of Eq. (49) can be written as

(50) F ( ξ ) F z b 1 0 = z b 1 0 + v t z 1 1 n G ( z v t ) v / c d z .

Substituting Eq. (50) into (49) yields the final expression

(51) t G = 1 n 1 v / c v z b 1 0 + v t z 1 1 n G ( z v t ) v / c d z n 1 z c + z z b 1 0 v .

Next, we determine the derivative of g G(t′) [Eq. (45)] with respect to t′ by applying the chain rule and introducing the substitution u = z + C G(t′), yielding

(52) d g G ( t ) d t = 1 v d F 1 ( u ) d u d C G ( t ) d t .

Let w = F −1(u), i.e., u = F(w). Then the derivative dF −1(u)/du on the right-hand side of Eq. (52) can be expressed as

(53) d F 1 ( u ) d u = d w d u = 1 d u d w = 1 F ( w ) = 1 F [ F 1 ( u ) ] ,

where F′(⋅) denotes the derivative of F with respect to its argument. Using Eq. (23), the derivative dC G(t′)/dt′ on the right-hand side of Eq. (52) is given by

(54) d C G ( t ) d t = v 1 n 1 v / c .

Substituting Eqs. (53) and (54) into (52), we obtain

(55) d g G ( t ) d t = 1 1 n 1 v / c 1 F [ F 1 ( u ) ] .

At the root t G , using Eq. (47), we find

(56) F 1 u t G = z v t = ξ .

Then, applying the definition in Eq. (20), we get

(57) F F 1 u t G = F ( ξ ) = 1 1 n G ( ξ ) v / c .

Substituting Eq. (57) into (55) and using the relation ξ = zvt, we find the derivative of g G(t′) evaluated at the root

(58) d g G ( t ) d t t G = 1 n G ( z v t ) v / c 1 n 1 v / c .

Finally, substituting Eqs. (45), (51), and (58), along with the relation f(0, t′) = E(0, t′), into the delta function identity [Eq. (46)], we evaluate the integral in Eq. (44) and find the result in Eq. (26).

Appendix B Derivation of Eq. (32)

In this appendix, we provide a detailed derivation of Eq. (32) given in Section 3.3. Substituting Eq. (31) into (3), we get

(59) E 2 ( z , t ) = δ t n 2 c z C 2 ( t ) E ( 0 , t ) d t .

We define the argument of the delta function in Eq. (59) as

(60) g 2 ( t ) = t n 2 c z C 2 ( t ) ,

where C 2(t′) is given in Eq. (30b).

To evaluate the convolution in Eq. (59), we use again the property [31]

(61) δ g 2 ( t ) f ( 0 , t ) d t = f ( 0 , t 2 ) | g 2 t 2 | ,

where g 2 t 2 = 0 and g 2 ( ) = d g 2 ( ) / d t . Substituting Eq. (30b) into (59) and solving g 2(t′) = 0 for t 2 , we obtain

(62) t 2 = 1 n 1 v / c 1 n 2 v / c t n 2 c z + n 2 n 1 1 n 1 v / c z b 1 0 c + D v 1 n 2 v / c v F z b 1 0 + D F z b 1 0 .

Applying Eq. (20) and performing the substitution z′ = ξ + vt, the difference in the square brackets of Eq. (62) can be written as

(63) F z b 1 0 + D F z b 1 0 = z b 1 0 + v t z b 1 0 + D + v t 1 1 n G ( z v t ) v / c d z .

Substituting Eq. (63) into (62) yields the final expression

(64) t 2 = 1 n 1 v / c 1 n 2 v / c t n 2 c z + n 2 n 1 1 n 1 v / c z b 1 0 c + D v 1 n 2 v / c v z b 1 0 + v t z b 1 0 + D + v t 1 1 n G ( z v t ) v / c d z .

The derivative of g 2(t′) [Eq. (60)] with respect to t′, evaluated at the root, is then found as

(65) d g 2 ( t ) d t t 2 = d g 2 ( t ) d t = 1 n 2 v / c 1 n 1 v / c .

Finally, substituting Eqs. (60), (64), and (65), along with the relation f(0, t′) = E(0, t′), into the delta function identity [Eq. (61)], we evaluate the integral in Eq. (59) and find the result in Eq. (32).

Appendix C FDTD validation

Figure 6 compares the results obtained by the theory in Figure 4 and by the FDTD method presented in [15], [24]. Specifically, it plots the field (top panels) and spectrogram (bottom panels) differences

(66a) δ E = | E x E FDTD | max ( E x ) × 100 %

and

(66b) δ S = | E x | 2 | E FDTD | 2 max | E x | 2 × 100 % ,

between the theoretical and simulation results, where E x corresponds to the theoretical result and E FDTD corresponds to the simulation result. As shown in Figure 6, the theoretical results closely match the simulation results (less than 5 % difference). This qualitatively validates the theory, while the fact that the theoretical results are exact (no approximation) suggests that the difference is due to the simulation (meshing) approximation errors, as we could verify by decreasing the mesh size up to the memory capability of our computer.

Figure 6: 
Difference between the exact theoretical results in Figure 4 and full-wave FDTD results, with difference attributed to simulation (meshing) approximation errors.
Figure 6:

Difference between the exact theoretical results in Figure 4 and full-wave FDTD results, with difference attributed to simulation (meshing) approximation errors.

Figure 7 shows the corresponding difference for the designed linear-chirping GRIN interfaces shown in Figure 5, leading to the same conclusion.

Figure 7: 
Same as in Figure 6 but for Figure 5.
Figure 7:

Same as in Figure 6 but for Figure 5.

References

[1] N. A. Estep, D. L. Sounas, J. Soric, and A. Alù, “Magnetic-free non-reciprocity and isolation based on parametrically modulated coupled-resonator loops,” Nat. Phys., vol. 10, no. 12, pp. 923–927, 2014. https://doi.org/10.1038/nphys3134.Search in Google Scholar

[2] S. Taravati, “Aperiodic space-time modulation for pure frequency mixing,” Phys. Rev. B, vol. 97, no. 11, 2018, Art. no. 115131. https://doi.org/10.1103/physrevb.97.115131.Search in Google Scholar

[3] P. K. Tien, “Parametric amplification and frequency mixing in propagating circuits,” J. Appl. Phys., vol. 29, no. 9, pp. 1347–1357, 1958. https://doi.org/10.1063/1.1723440.Search in Google Scholar

[4] J. B. Pendry, E. Galiffi, and P. A. Huidobro, “Gain in time-dependent media – a new mechanism,” J. Opt. Soc. Am., vol. 38, no. 11, pp. 3360–3366, 2021. https://doi.org/10.1364/josab.427682.Search in Google Scholar

[5] A. Shlivinski and Y. Hadad, “Beyond the Bode-Fano bound: wideband impedance matching for short pulses using temporal switching of transmission-line parameters,” Phys. Rev. Lett., vol. 121, no. 20, 2018, Art. no. 204301. https://doi.org/10.1103/physrevlett.121.204301.Search in Google Scholar

[6] F. Biancalana, A. Amann, A. V. Uskov, and E. P. O’Reilly, “Dynamics of light propagation in spatiotemporal dielectric structures,” Phys. Rev. E, vol. 75, no. 4, 2007, Art. no. 046607. https://doi.org/10.1103/physreve.75.046607.Search in Google Scholar PubMed

[7] Z.-L. Deck-Léger, A. Akbarzadeh, and C. Caloz, “Wave deflection and shifted refocusing in a medium modulated by a superluminal rectangular pulse,” Phys. Rev. B, vol. 97, no. 10, 2018, Art. no. 104305. https://doi.org/10.1103/physrevb.97.104305.Search in Google Scholar

[8] M. A. Gaafar, T. Baba, M. Eich, and A. Y. Petrov, “Front-induced transitions,” Nat. Photonics, vol. 13, no. 11, pp. 737–748, 2019. https://doi.org/10.1038/s41566-019-0511-6.Search in Google Scholar

[9] J. Sloan, N. Rivera, J. D. Joannopoulos, and M. Soljačić, “Controlling two-photon emission from superluminal and accelerating index perturbations,” Nat. Phys., vol. 18, no. 1, pp. 67–74, 2022. https://doi.org/10.1038/s41567-021-01428-4.Search in Google Scholar

[10] A. Ball, et al.., “A space-time knife-edge in epsilon-near-zero films for ultrafast pulse characterization,” Laser Photonics Rev., vol. 19, no. 5, 2025, Art. no. 2401462. https://doi.org/10.1002/lpor.202401462.Search in Google Scholar

[11] L. Felsen and G. Whitman, “Wave propagation in time-varying media,” IEEE Trans. Antennas Propag., vol. 18, no. 2, pp. 242–253, 1970. https://doi.org/10.1109/tap.1970.1139657.Search in Google Scholar

[12] R. Fante, “Transmission of electromagnetic waves into time-varying media,” IEEE Trans. Antennas Propag., vol. 19, no. 3, pp. 417–424, 1971. https://doi.org/10.1109/tap.1971.1139931.Search in Google Scholar

[13] Y. Hadad and A. Shlivinski, “Soft temporal switching of transmission line parameters: wave-field, energy balance, and applications,” IEEE Trans. Antennas Propag., vol. 68, no. 3, pp. 1643–1654, 2020. https://doi.org/10.1109/tap.2020.2967302.Search in Google Scholar

[14] Z. Hayran, J. B. Khurgin, and F. Monticone, “ℏω versus ℏk: dispersion and energy constraints on time-varying photonic materials and time crystals,” Opt. Mater. Express, vol. 12, no. 10, pp. 3904–3917, 2022. https://doi.org/10.1364/ome.471672.Search in Google Scholar

[15] Z.-L. Deck-Léger, A. Bahrami, Z. Li, and C. Caloz, “Generalized FDTD scheme for the simulation of electromagnetic scattering in moving structures,” Opt. Express, vol. 31, no. 14, pp. 23 214–23 228, 2023. https://doi.org/10.1364/oe.493099.Search in Google Scholar

[16] S. J. Orfanidis, Electromagnetic Waves and Antennas, New Brunswick, NJ, USA, Rutgers University, 2002.Search in Google Scholar

[17] L. Novotny and B. Hecht, Principles of Nano-Optics, 2nd ed. Cambridge, UK, Cambridge University Press, 2012.10.1017/CBO9780511794193Search in Google Scholar

[18] B. E. Saleh and M. C. Teich, Fundamentals of Photonics, 3rd ed. Hoboken, NJ, USA, John Wiley & Sons, 2019.Search in Google Scholar

[19] Y. Xiao, G. P. Agrawal, and D. N. Maywar, “Spectral and temporal changes of optical pulses propagating through time-varying linear media,” Opt. Lett., vol. 36, no. 4, pp. 505–507, 2011. https://doi.org/10.1364/ol.36.000505.Search in Google Scholar

[20] L. A. Ostrovskiĭ, “Some “moving boundaries paradoxes” in electrodynamics,” Sov. Phys. Usp., vol. 18, no. 6, p. 452, 1975. https://doi.org/10.1070/pu1975v018n06abeh001967.Search in Google Scholar

[21] C. Caloz and Z.-L. Deck-Léger, “Spacetime metamaterials – part II: theory and applications,” IEEE Trans. Antennas Propag., vol. 68, no. 3, pp. 1583–1598, 2019.10.1109/TAP.2019.2944216Search in Google Scholar

[22] G. P. Agrawal, Nonlinear Fiber Optics, 6th ed. San Diego, CA, USA, Academic Press, 2019.Search in Google Scholar

[23] A. Shlivinski and Y. Hadad, “Universal transient radiation dynamics by abrupt and soft temporal transitions in optical waveguides,” Nanophotonics, vol. 14, no. 6, pp. 785–793, 2025. https://doi.org/10.1515/nanoph-2024-0525.Search in Google Scholar PubMed PubMed Central

[24] A. Bahrami, Z.-L. Deck-Léger, Z. Li, and C. Caloz, “A generalized FDTD scheme for moving electromagnetic structures with arbitrary space-time configurations,” IEEE Trans. Antennas Propag., vol. 72, no. 2, pp. 1721–1734, 2024.10.1109/TAP.2023.3332491Search in Google Scholar

[25] K. De Kinder and C. Caloz, “Electromagnetic scattering at an arbitrarily accelerated interface,” in 2024 Eighteenth International Congress on Artificial Materials for Novel Wave Phenomena (Metamaterials), 2024, pp. 1–3.10.1109/Metamaterials62190.2024.10703290Search in Google Scholar

[26] K. De Kinder, A. Bahrami, and C. Caloz, “Scattering and chirping at accelerated interfaces,” arXiv preprint arXiv:2506.19575, 2025. https://arxiv.org/abs/2506.19575.Search in Google Scholar

[27] J. Ran, et al.., “Realizing tunable inverse and normal Doppler shifts in reconfigurable RF metamaterials,” Sci. Rep., vol. 5, no. 1, 2015, Art. no. 11659. https://doi.org/10.1038/srep11659.Search in Google Scholar PubMed PubMed Central

[28] S. Taravati and C. Caloz, “Mixer-duplexer-antenna leaky-wave system based on periodic space-time modulation,” IEEE Trans. Antennas Propag., vol. 65, no. 2, pp. 442–452, 2016. https://doi.org/10.1109/tap.2016.2632735.Search in Google Scholar

[29] H. Moussa, G. Xu, S. Yin, E. Galiffi, Y. Ra’di, and A. Alù, “Observation of temporal reflection and broadband frequency translation at photonic time interfaces,” Nat. Phys., vol. 19, no. 6, pp. 863–868, 2023. https://doi.org/10.1038/s41567-023-01975-y.Search in Google Scholar

[30] C. Caloz, Z.-L. Deck-Léger, A. Bahrami, O. C. Vicente, and Z. Li, “Generalized space-time engineered modulation (GSTEM) metamaterials: a global and extended perspective,” IEEE Antennas Propag. Mag., vol. 65, no. 4, pp. 50–60, 2023. https://doi.org/10.1109/map.2022.3216773.Search in Google Scholar

[31] B. R. Kusse and E. A. Westwig, Mathematical Physics: Applied Mathematics for Scientists and Engineers, 2nd ed. Weinheim, Germany, Wiley VCH, 2006.10.1002/9783527618132Search in Google Scholar

Received: 2025-07-07
Accepted: 2025-08-12
Published Online: 2025-09-05

© 2025 the author(s), published by De Gruyter, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 6.1.2026 from https://www.degruyterbrill.com/document/doi/10.1515/nanoph-2025-0308/html
Scroll to top button