Generalized temporal transfer matrix method: a systematic approach to solving electromagnetic wave scattering in temporally stratified structures

3.1 TTMM formalism

Before establishing the TTMM formalism, let us define the temporal system (Figure 3d). It consists of an unbounded homogenous medium, whose permittivity or permeability undergoes abrupt changes n times at t ( i ) ( i = 1,2,3 … n ) . It can be regarded as a ‘temporal stratified structure’ consisting of n − 1 temporal layers, whose durations are τ ( i ) = t ( i ) − t ( i − 1 ) ( i = 2,3,4 … n ) . Moreover, an EM wave is incident before t ( 1 ) .

Several assumptions are made about the system: (a) The anisotropy of each temporal layer shares the same principal axes ( X Y Z ), and (b) in the first temporal layer (i.e., when the wave is incident), the material is isotropic, which guarantees that k , D , and B are perpendicular to each other. These two assumptions are the foundations of the coordinate transformation. In S1, the material property in the i -th layer is expressed as

where ϵ ′ ( i ) and μ ′ ( i ) denote relative permittivity and permeability (here, and in all of the following discussions, the subscript ‘r’ is dropped, for simplicity). Notice that the explicit expressions for ϵ and μ depend on the choice of coordinate system. Therefore, in S3, we have

where A is the rotation matrix between S1 and S3:

Notice that ϵ ( i ) and μ ( i ) have off-diagonal terms. For the remainder of this section, all the physical quantities and expressions derived are considered in the framework of the S3 coordinate system, unless specified otherwise. Based on Maxwell’s equations, the D field should satisfy the relation:

where k = ( k x , k y ) is the wave vector (a known, unchanged value in homogenous media), and ω is the radian frequency.

Similar to the discussions in [24, 40], one can express the D and B field in the i -th layer as linear combinations of the four modes:

where { A 1 ( i ) , A 2 ( i ) , A 3 ( i ) , A 4 ( n ) } are a set of expansion coefficients of the modes. The terms d σ ( i ) and b σ ( i ) denote the unit field vectors of D and B inside the i -th layer. According to Maxwell’s equations, the following relation is satisfied:

By applying the boundary conditions at the i -th temporal boundary:

we have

where P n is the temporal propagation matrix defined as (Eq. S5 in [24]):

and D i is defined as (similar to Eq. S4 in [24] and Eq. (14) in [40]):

Importantly, at this point in the development, it is not clear what role the incidence angle would play in a generalized formulation. In the following section, however, we will see that when the wave is obliquely incident (i.e., S1 and S3 do not overlap), the relation between k and ω σ ( i ) will be different. Therefore, the explicit form of the dispersion relation (Eq. (5)) and the expression for the matrix D n (Eq. (10)) will also change accordingly.

3.2 Explicit form of the transfer matrices

First, let us investigate the dispersion relation, as determined by Eq. (5). For each temporal layer, we replace ϵ and μ with ϵ ( i ) and μ ( i ) in Eq. (5), then replace D with D ( i ) as defined in Eq. (6.1). It can be shown that there are four non-trivial solutions of the dispersion relation: k ∼ ω σ ( i )   ( σ = 1,2,3,4 ) . After carrying out the appropriate mathematical analysis (see Supplemental Document 1A), we arrive at

where c 1 = cos ( θ ) and s 1 = sin ( θ ) . Also, we introduce the following definitions:

These quantities can be viewed as the equivalent refractive indices seen by x- and y-polarized waves inside the i -th temporal layer.

It then follows that the corresponding field vectors of D are

Clearly, these four solutions correspond to the two orthogonal modes (i.e., y- and x-polarized waves) propagating in opposite directions. Next, we use Eqs. (6.2), (10) and (11.4) to obtain:

which can be rewritten as:

At this point we define the following identities:

which can be viewed as the equivalent impedances seen by the x- and y-polarized waves inside the i -th temporal layer (normalized by the free space impedance Z 0 ). It is obvious that, if θ = 0 (i.e., normal incidence), then Eq. (12.3) will reduce to the definition of impedance in the normal sense, which is the same as in Eq. (S13) of [24] if i = 2 . It is worth mentioning that Eq. (12.1) and (12.2) hold only if the material is anisotropic. If the material is isotropic, the D matrix can take different forms depending on the angle ϕ . This is elaborated on in the Supplemental Document 1B.

With the explicit form of matrices D and P given, one can calculate the total transfer matrix of the temporal structure illustrated in Figure 3d:

From the discussion above, we know that the Q matrix is solely determined by ϵ , μ , and the durations of each temporal layer, which is referred to as the (temporal) profile of this temporal structure.

At this point we have completed our introduction to the GTTMM. In the next two sections, we will apply this tool to four different temporal systems, which are schematically represented in Figure 4. All of these systems can be regarded as temporal multilayer structures, as have been described in Figure 3d. Moreover, the desired physical quantities, for example, the permittivity tensor of the anisotropic antireflection temporal coating or the transmission coefficients in the case of polarization conversion, could be calculated directly from the Q matrix. This procedure will be elaborated on in the following sections.

Figure 4:

Schematics of several temporal systems.

The different temporal regions are represented by colored boxes. The temporal profile in the case of (a) an arbitrary temporal structure, (b) an anisotropic antireflection temporal coating, (c) polarization conversion, and (d) redirection of energy propagation. Notice that anisotropic temporal regions are denoted by checkerboard patterns.

4.1 Arbitrary multi-layer temporal structure

First, in order to demonstrate the validity of our GTTMM algorithm, we consider a 6-layer structure, which has a random profile. The information on the composition of this structure is displayed in Table 1. Notice that the first and last layers are isotropic, and they are semi-infinite in time. Now, let us study the interactions between the electromagnetic waves and this structure. As mentioned before, two independent parameters, θ and ϕ , are required to determine the orientation of the incident wave. Then we define two orthogonal polarization modes as s and p , as depicted in Figure 5a. The orientation of the latter polarization state can be rotated counter-clockwise from the former one by 90 ° along k . To this end, any polarization state of the incident, reflected, and transmitted waves can be written as a linear combination of the two modes. Hence, one can naturally define the S parameters of a given temporal structure as t i j and r i j   ( i , j = s , p ) . Please refer to Supplemental Document 3 for the detailed calculation methodology for the S parameters.

Figure 5 shows the spectral response of a p -polarized wave ‘propagating’ through this temporal structure. The dots and solid lines represent the results obtained by the FDTD simulations and the GTTMM theory, respectively. It can be clearly seen that, for all four S parameters (i.e.,  t p p , t p s , r p p , r p s ), the numerical and theoretical results match very well over a wide frequency range. Considering the random nature of the temporal structure, this analysis confirms the validity and robustness of the GTTMM. In fact, if the material becomes isotropic, then this case would reduce to Example 1 in [38].

4.2 Polarization conversion

Polarization conversion is an important property with many practical applications in optics and electromagnetics. Traditionally, one could convert polarization of an EM wave by utilizing the interface between an isotropic and an anisotropic material [43]. In our previous work [24], we extended this idea to the temporal domain, and achieved complete polarization conversions in real time, as schematically demonstrated in Figure 4c. In that work, however, we only considered the normal incidence case ( θ = 0 ). What if the wave is obliquely incident ( θ ≠ 0 )? Specifically, we study the temporal system which has a material profile illustrated in Figure 6a. It comprises an anisotropic temporal layer, which has a duration Δ t ( 2 ) and permittivity ϵ ( 2 ) = diag ( ϵ X ( 2 ) , ϵ Y ( 2 ) , ϵ Z ( 2 ) ) , embedded in an isotropic medium with a permittivity ϵ ( 1 ) = 1 . The Q matrix of the system can be represented as

Figure 6:

Illustrations of the polarization conversion effect. (a) Schematic of the temporal profile of the system (b)–(d) the transmission spectra.

After some mathematical manipulations, we arrive at the following result:

for which

where n j ′ ( 2 ) and Z j ′ ( 2 ) have been defined in Eq. (11.3) and (12.3), respectively. Based on these equations, it follows that the property of the anisotropic layer ( Δ t ( 2 ) , ϵ ( 2 ) ) will solely determine the response of this temporal system to EM waves. Notice that if θ = 0 , Eqs. (15) and (16) will reduce to the known results, i.e., Eq. (2) in [24]. In order to better illustrate the polarization conversion effect, let us consider the special case of Complete Polarization Conversion (CPC):

CPC means that the incident wave is completely converted to the orthogonal polarization for certain values of θ , ϕ and frequency, which are denoted as θ 0 , ϕ 0 , and f 0 . Assuming that the incident wave is p -polarized, the CPC condition implies that t p p ( θ 0 , ϕ 0 , f 0 ) = 0 . Then, we choose values for Δ t ( 2 ) and ϵ ( 2 ) that satisfy the following criteria: Δ t ( 2 ) = 1.11 / f 0 , ϵ ( 2 ) = diag ( 2.5 , 12.5 , 1.25 ) . Since t p p is a complicated function of θ ,

ϕ and f , it is convenient to visualize the theoretical results from different perspectives. In Figure 6b, t p p is plotted versus f and θ with ϕ = ϕ 0 . While in Figure 6c, t p p is plotted versus f and ϕ with θ = θ 0 . Then in Figure 6d, t p p and t p s are plotted versus f with θ = θ 0 and ϕ = ϕ 0 . FDTD simulations were also performed as a validation, and they agree well with the analytical results. From these figures, one can clearly observe that the desired CPC phenomenon is achieved. From this example, we find that, by generalizing to oblique incidence, there is greater flexibility in real-time polarization control compared to what has been achieved in [24].

4.3 Anisotropic antireflection temporal coating (AATC)

V. Pacheco-Pena et al. recently introduced the intriguing concept of an antireflection temporal coating (ATC) in the time domain [15]. We know from the theory of temporal boundaries that there is a reflection when the material permittivity undergoes a sudden change (e.g., from ϵ ( 1 ) to ϵ ( 1 ) ). However, one can reduce the reflection to zero by ‘inserting’ an intermediate temporal region whose permittivity and duration satisfy the following relations:

where n = 1,2,3 … and ω 0 is the radian frequency of the wave when it is inside the intermediate layer.

In [15], however, the authors only consider the case where the material is isotropic. One may be naturally curious as to what conditions the material properties should satisfy in order to minimize reflection if anisotropic materials are involved. In other words, is it possible to derive a similar expression to Eq. (18) using the GTTMM?

To be more specific, we consider the following temporal structure comprising of three temporal regions, whose permittivities are ϵ ( 1 ) , ϵ ( 2 ) = ( ϵ X ( 2 ) , ϵ Y ( 2 ) , ϵ Z ( 2 ) ) , and  ϵ ( 3 ) (see Figure 6a). The second anisotropic layer acts as an anti-reflection temporal coating, which we call an anisotropic antireflection temporal coating (AATC). Our goal is to calculate the permittivity and duration of the AATC, so that there is no net reflection for a certain frequency f 0 , angle θ 0 and ϕ 0 ((i.e.,  r p p ( θ 0 . ϕ 0 , f 0 ) = r p s ( θ 0 . ϕ 0 , f 0 ) = 0 )), assuming that the incident wave is p - polarized. After some mathematical manipulations (see Supplemental Document 4), we arrive at:

where the definition of ω 1 ( 2 ) and ω 3 ( 2 ) have been given in Eq. (11.2). These criteria are easy to interpret: Eq. (19.1) is the impedance matching condition, while Eq. (19.2) represents the π / 2 phase accumulation condition. Interestingly, ϕ does not appear in Eq. (19). This means that, as long as the incident angle is θ , there will be no reflection regardless of the orientation of the D and B fields. Moreover, if the wave is normally incident (i.e.,  c 1 = 1 , s 1 = 0 ), then Eq. (19) reduces to Eq. (18), assuming ϵ X ( 2 ) = ϵ Y ( 2 ) .

Next, for a practical illustration, we arbitrarily choose a combination of parameters that satisfy Eq. (19): ε ( 1 ) = 1 , ε ( 2 ) = ( 1.5 , 2,3 ) , ε ( 3 ) = 4 , Δ t ( 2 ) = 2 / 4   T 0 , and θ 0 = π / 4 . The temporal profile is depicted by the dotted lines in Figure 7a.

Figure 7:

Illustration of an anisotropic anti-reflection temporal coating (AATC).

(a) The temporal profiles with and without the AATC. (b) Transmission and reflection spectra with and without the AATC, where the GTTMM and FDTD results are represented by solid lines and dots, respectively. (c) and (e) FDTD simulation results of E p ( z , t ) with and without the AATC, at central frequency f 0 . (f) A magnified plot of a portion of (c), which is indicated by the black rectangle. (d) FDTD simulation results of E p ( z = z 1 , t ) with and without the AATC, where the position of z 1 is indicated by the red dashed lines in (c) and (e).

Similar to the approach adopted in Section 4A, we have plotted the spectral response of a p wave propagating in this temporal structure (see Figure 7b). From this figure, we can observe that the reflection ( r p p ) drops to zero at the central frequency f 0 = 1 / T 0 . Notice that the cross-polarization S parameters are not presented because r p s is zero across the entire frequency range, if Eq. (19) is satisfied (see a detailed explanation of this phenomenon in Supplemental Document 4). Moreover, we show the real time electric field E p ( z , t ) in Figure 7c, where the incident wave is a narrow band Gaussian pulse with a center frequency f 0 . From this figure, we find that there is no reflection, after the wave meets two temporal boundaries.

As a comparison, we consider another case where ϵ changes suddenly from 1 to four at an instant in time (i.e., the case without an AATC, see the black solid line in Figure 7a). The corresponding S parameters and real time electric fields are also plotted in Figure 7b and e. It is important to note that r p p ≠ 0 at f = f 0 in this case, which is expected. Besides, in order to better understand the effect of the AATC, we show the electric field distribution at a certain position z 1 for both with and without AATC cases, in the same plot (Figure 7d). It is seen that the existence of the AATC effectively reduces reflected waves.

This example represents a natural but highly nontrivial generalization of the work reported in [15]. By comparing with the work in [15], our presented results cannot be easily understood from a space-time symmetry perspective. Rather, a rigorous GTTMM analysis is required to derive the conditions for reflection cancellation.

4.4 Redirection of energy flow

One of the interesting features of temporal boundaries is that the D and B fields are conserved before and after the boundary. This means that the E and H fields will follow the same changes with time as the anisotropy of the material. Hence, the energy flux of the wave (i.e., the Poynting vector), defined as S = E × H , could also change with time. V. Pacheco-Pena et al. investigated this phenomenon in [29] and called it ‘temporal aiming’. The authors of that work, however, only considered a special case when ϕ = 0 . One may naturally ask the question: what happens if ϕ ≠ 0 ? To answer this question, let us first define the temporal system as shown in Figure 8a. The material is isotropic ( ϵ ( 1 ) = μ ( 1 ) = 1 ) in the beginning, and then suddenly switches to an anisotropic form where ϵ ( 2 ) = diag ( ϵ X ( 2 ) , ϵ Y ( 2 ) , ϵ Z ( 2 ) ) ) at t = t ( 1 ) . Here we assume that the material is non-magnetic for simplicity. The orientation of the incident wave is described by the angle θ and ϕ , as depicted in Figure 3b. The Poynting vectors of the incident and transmitted waves (i.e., before and after the temporal boundary) are S ( 1 ) and S ( 2 ) , respectively, and their angle is Δ θ S = < S ( 1 ) , S ( 2 ) > .

Figure 8:

Demonstrations of the redirection of energy propagation.

(a) The temporal profile of the system. (b) Δ θ S  as a function   of   θ for ϕ = 0 (blue) and ϕ = 60 ° (magenta), as dictated by Eq. (22). (c)–(f) Normalized E X field (color) and Poynting vector (black arrow) of the incident wave at = t ( 1 ) − , and of the transmitted wave at t = t ( 1 ) + , for the two cases. The spatial positions of the plots are chosen to contain the center of the incident Gaussian beam at the given times. In (e) and (f), the physical quantities ( E and S fields) are plotted on a plane parallel to both S ( 1 ) and S ( 2 ) in order to depict the redirection angle Δ θ S between them.

Next, let us calculate the value of Δ θ S . Notice that the medium is isotropic before t ( 1 ) , such that S ( 1 ) is parallel to k and Δ θ S = < k , S ( 2 ) > . Recalling the coordinate transformation technique introduced in Section 3A (Figure 3c), we have Δ S = < z , S ( 2 ) > . This means that, if we express S ( 2 ) in x − y − z coordinates, S ( 2 ) = ( S x ( 2 ) , S y ( 2 ) , S z ( 2 ) ) , then Δ θ S can be easily calculated using the following expression:

where S ( 2 ) = E ( 2 ) + × H ( 2 ) + . Notice that the ‘+’ sign in the superscript denotes the transmitted wave. The transmission coefficients can be computed using the Q matrix of the system:

After simplification (see Supplemental Document 3), we have

Clearly, Δ θ S is a function of both the material anisotropy ( ϵ ( 2 ) ), and the orientation of the incident wave ( θ , ϕ ). It is interesting to observe that if ϕ = 0 , then Eq. (22) will reduce to Δ θ S = tan − 1 ( tan ( θ ) ⋅ ϵ X ( 2 ) / ϵ Z ( 2 ) ) , which is in agreement with Eq. (2) in [29].

Now, we consider a specific material profile: ϵ ( 2 ) = diag ( 2,5,20 ) , which is schematically illustrated in Figure 8a. Using Eq. (22), Δ θ S is plotted against θ for ϕ = 0 and ϕ = 60 ° in Figure 8b. Following this figure, we study two cases: (1) θ = 45 ° , ϕ = 0 , and (2) θ = 45 ° , ϕ = 60 ° , which correspond to Δ θ S = 39 ° and 27 ° , respectively. In Figure 8c–f, we present overlapping plots of the electric fields and the instantaneous Poynting vectors for the two cases just before and after the temporal boundary (at t = t ( 1 ) + and t = t ( 1 ) − ). We observe that the normal of the phase front ( k ) and Poynting vector ( S ) are along the same direction at t = t ( 1 ) − , while they have a nonzero angle at t = t ( 1 ) + . Using the simulated data from lumerical FDTD, we confirm that this angle matches very well the theoretical predications in Figure 8b. Moreover, it is interesting to note that in case 1, S ( 1 ) and S ( 2 ) lie in the XZ plane, therefore, Figure 8c–d are rendered in 2D. For case 2, however, S ( 1 ) and S ( 2 ) are not along any specific axial direction and, therefore, the results in Figure 8e and f are displayed on a plane in 3D space. From the results, one can clearly see that  ϕ ≠ 0 represents a highly nontrivial scenario, which includes the special case of ϕ = 0 considered in [29].