Temporal aiming

Deflecting and changing the direction of propagation of electromagnetic waves are needed in multiple applications, such as in lens–antenna systems, point-to-point communications and radars. In this realm, metamaterials have been demonstrated to be great candidates for controlling wave propagation and wave–matter interactions by offering manipulation of their electromagnetic properties at will. They have been studied mainly in the frequency domain, but their temporal manipulation has become a topic of great interest during the past few years in the design of spatiotemporally modulated artificial media. In this work, we propose an idea for changing the direction of the energy propagation of electromagnetic waves by using time-dependent metamaterials, the permittivity of which is rapidly changed from isotropic to anisotropic values, an approach that we call temporal aiming. In so doing, here, we show how the direction of the Poynting vector becomes different from that of the wavenumber. Several scenarios are analytically and numerically evaluated, such as plane waves under oblique incidence and Gaussian beams, demonstrating how proper engineering of the isotropic—anisotropic temporal function of εr(t) can lead to a redirection of waves to different spatial locations in real time. Theoretical analysis shows that beam-steering of electromagnetic waves can be accomplished by temporally changing the permittivity of metamaterials between isotropic to anisotropic values. The approach, called “temporal aiming”, has been formulated by Victor Pacheco-Peña from Newcastle University, UK and Nader Engheta from University of Pennsylvania in the US. In principle, it could open new opportunities for the flow of information around integrated photonic circuits, with flat metamaterial elements deflecting electromagnetic waves to specific targets or receivers on an optical chip as desired. Simulations performed with the software COMSOL with both plane, monochromatic waves and more complex Gaussian beams confirm the feasibility of the approach. It is proposed that the required temporal changes in the metamaterial’s relative permittivity could be achieved by the use of tunable metasurfaces or transmission lines with time-varying circuit elements.


Introduction
Achieving arbitrary control of electromagnetic wave propagation has been of great interest within the scientific community for many years. It is well known that carefully engineered spatially varying geometries and materials can be implemented to manipulate wave-matter interactions 1 . This spatial control of waves is, in fact, the mechanism behind the development of many applications we often use on a regular basis, such as lenses, sensors and radars. The field of antennas has also benefited from this spatial control of wave propagation, where, for instance, in its basic configuration, it is possible to change the direction of a transmitted wave by mechanically modifying the spatial location of the transmitter, a technique known as mechanical beam steering 2 .
The beam steering of electromagnetic waves has also benefited from the introduction of metamaterials, where it has been shown that electromagnetic radiated beams can be redirected by changing the position of the source in a metalens-antenna system, by locally designing the phase of the unit cells in the metamaterials, or by real-time tuning of the effective electromagnetic parameters of the metastructures [31][32][33][34][35][36] , among other techniques. These steering properties are important in different areas, such as in point-to-point communications and radars, where the spatial aiming of targets is required.
Metamaterials and metasurfaces have so far been studied mostly in the time-harmonic scenario, where wave propagation is controlled by engineering geometries and materials in the spatial region, i.e., spatial inhomogeneity, in which the wave is travelling. Recently, the temporal modulation of metamaterials has also gained growing attention within the scientific community [37][38][39][40][41][42][43][44] , as changing the electromagnetic properties (ε, µ) of metamaterials both in space (x, y, z) and time (t) can offer full fourdimensional spatiotemporal control of wave-matter interactions. It is important to highlight that the interaction of electromagnetic waves in a time-modulated medium has been of great interest in the scientific community for several decades, where, for instance, in the last century, it was considered a time-dependent relative permittivity ε r (t) that is rapidly changed in time from one positive value ε r1 (greater than unity) to a different greater-than-unity positive value ε r2 40,41 . With this configuration, it was demonstrated that a set of two waves is created at this temporal boundary, one of which is travelling forward (FW), and the other, backward (BW). Remarkably, an analogy between this temporal and the spatial interface between two materials with different electromagnetic parameters was demonstrated, showing that the two waves created at a temporal boundary (FW and BW) are the temporal equivalent/analogy to the transmitted and reflected waves in spatial interfaces. In this realm, temporal and spatiotemporal metamaterials have recently been proposed and applied in several exciting and intriguing applications, such as effective medium theory 45 , inverse prisms 46 , nonreciprocity 47,48 , anti-reflection temporal coatings 49 , frequency conversion 50 and time reversal 51 .
Motivated by the exciting possibilities and opportunities opened up by the spatiotemporal modulation of metamaterials in four dimensions (x, y, z, t) and the importance of the beam steering of electromagnetic waves for different applications, in this work, we introduce the concept of temporal aiming as the temporal analogue of spatial aiming. First, the fundamental physics of the proposed technique are presented, considering an oblique incident p-polarised monochromatic plane wave propagating in an unbounded medium with a time-dependent permittivity ε r (t). A temporal interface is introduced by changing the relative permittivity ε r (t) from an isotropic positive value ε r1 to an anisotropic permittivity ε r2 = {ε r2x , ε r2z } (all with positive values greater than unity) at time t = t 1 . In so doing, the wavenumber k is preserved before and after the temporal change, but the direction of the energy (Poynting vector S) is modified to a different angle compared with k. The dependence of the new direction of the Poynting vector on the incident angle before the temporal change of ε and the values of the permittivity tensor ε r2 is presented and discussed. Next, we study more complex scenarios by using Gaussian beams under normal and oblique incidence to consider the case of multiple plane waves travelling in different directions. With this set-up, it is demonstrated how the direction of the energy is modified to multiple angles when using an isotropic-toanisotropic temporal boundary. Finally, the temporal aiming technique is numerically demonstrated by using a narrowband wavepacket under oblique incidence. As will be shown, this transmitted wavepacket can be "ushered" and "re-directed" in real time to reach different receivers placed at different spatial locations by properly engineering the ε r (t) of the background medium. All the results presented here are compared with numerical simulations, demonstrating good agreement with the design and analytical calculations.

Results
Temporal aiming: isotropic-to-anisotropic change in ε r First, let us discuss the spatial scenario schematically described in Fig. 1a. In this case, we can consider a monochromatic continuous wave (CW) that is being emitted by a source located on the xz plane and tilted to an angle θ 1 . The background medium is homogeneous, isotropic and time-independent, with a relative permeability µ r (t) = µ r1 and relative permittivity ε r (t) = ε r1 (see the inset of Fig. 1a). As is well known, the wave emitted from this source will also be tilted with the wavenumber k and Poynting vector S parallel, travelling along the same direction defined by θ 1 . Now, let us place two receivers at two different spatial locations (Rx1 and Rx2), as schematically shown in Fig. 1a. As observed, this is not the best scenario if one needs to send the emitted wave to either of the two receivers since the source is not aligned to Rx1 or to Rx2. However, as described in the introduction, the most simple yet effective way to reach either Rx1 or Rx2 is to perform mechanical beam steering 2,32 . In this technique, the source is placed on a translation stage, and then spatially shifted to different locations on the xz plane to tilt the emitted wave to the correct angle θ 1 such that it reaches either Rx2 (Fig. 1b) or Rx1 (Fig. 1c). In addition to this technique, there are other effective and more sophisticated alternatives that can be used to steer electromagnetic waves, such as phased arrays (in which the phase shifts between the antenna elements can be changed in real time) and metamaterial-based antennas with tuneable properties 31,34,36 . In this context, steering electromagnetic waves is considered a key feature in applications where the spatial aiming of targets is needed, such as radars and point-to-point communications, as explained before.
The spatial aiming described in Fig. 1a-c has been discussed considering the time-harmonic scenario (frequency domain), taking into account that the relative permittivity and permeability of the background medium are constant [ε r (t) = ε r1 , µ r (t) = µ r1 ]. With this in mind, one may ask the following: would it be possible to change the direction of the energy propagation of an emitted monochromatic CW plane wave in time by considering a time-dependent permittivity ε r (t) ≠ ε r1 of the background medium? If such temporal aiming is possible, what type of function for ε r (t) may we use to achieve it? To answer these questions, in this work, we focus our attention on the rapid change in relative permittivity (approximately modelled mathematically as a "step function"), which is rapidly modified in time from an initial value (equal to or greater than unity) ε r1 to another greater-than-unity ε r2 at a time t = t 1 , considering that the fall/rise time is smaller than the period T of the incident wave. (Strictly speaking, one needs to take into account the material dispersion in these scenarios. However, if we assume that the material resonance frequencies are much larger than the frequency of operation, we can approximately assume the materials to be dispersionless).
Let us study the case schematically shown in Fig. 1d. We again assume a p-polarised monochromatic CW plane wave travelling in an unbounded medium with an incident angle θ 1 . Let us first consider that the permittivity is time-dependent ε r (t) and is isotopically changed from a positive value ε r1 to another positive value ε r2 . This scenario was studied in the last century, and it was shown that this ε r (t) can induce a temporal boundary/interface that generates a FW (temporal transmission) and a BW Fig. 1 Schematic representation and comparison of spatial and temporal aiming. a-c Sketch of the conventional mechanical beam-steering technique, where the beam emitted by a source can be directed to either receiver 1 (Rx1) or receiver 2 (Rx2) by simply placing the source on a translation stage and moving it on the xz plane. d-f Temporal aiming sketch, considering a source emitting a monochromatic wavepacket embedded in a time-dependent medium, where its relative permittivity is changed from an isotropic value ε r1 to a tensor ε r2 = [ε r2x , ε r2z ] at t = t 1 . Vectors S and k are parallel before the change in ε r for t < t 1 , d, and non-parallel for t > t 1 , e-f. The angle of the vector S can be designed to reach either e, receiver 2 or f, receiver 1, depending on the tensor ε r2 = [ε r2x , ε r2z ] (temporal reflection) wave travelling with the same angle as that of the incoming wave, i.e., the wavenumber k is preserved while the frequency is changed from f 1 to f 2 = ( ffiffiffiffiffi ffi ε r1 p / ffiffiffiffiffi ffi ε r2 p )f 1 . This type of time-dependent ε r (t) was then recently used to propose exciting applications, such as time reversal and anti-reflection temporal coatings, as explained in the introduction. Nevertheless, since k is the same before (t = t 1 − ) and after (t = t 1 + ) the change in ε from ε r1 to ε r2 and the direction of propagation for both k and energy S is not modified, it is straightforward to conclude that this type of ε r (t) is not suitable for the temporal aiming described in Fig. 1, since our aim is to redirect the energy of emitted waves to different receivers placed at different spatial locations. Now, what if we change the relative permittivity ε r from an isotropic value ε r1 to an anisotropic permittivity tensor at t = t 1 such that ε r2 = [ε r2x , ε r2z ] with ε r2x ≠ ε r2z , both of which are positive, real-valued and greater-thanunity parameters? The schematic representation of this ε(t) is shown in the inset of Fig. 1d. Note that here, we consider only the z and x components of the permittivity tensor since we have a TM polarisation with the electric field lying on the xz plane. Akbarzadeh et al. recently explored this function of ε r (t) to create what they aptly called the "inverse prisms", demonstrating that vector k is again preserved as the isotropic case while the frequency is modified, but this time with a value depending on the tensor ε r2 and the incident angle θ 1 46 , hence their coined name, i.e., "inverse prism". However, one may ask an intriguing question: what will happen to the Poynting vector S in this scenario? Would it be possible to exploit this isotropic-to-anisotropic temporal variation of ε(t) for temporal aiming, as we propose here?
To answer these questions, let us analytically evaluate this case, as shown in Fig. 1d. For t < t 1 , the background medium is isotropic and homogeneous with relative permittivity and relative permeability values ε r1 and µ r1 , respectively. With this set-up (assuming the e (iωt) time convention), the magnetic and electric fields are and c being the velocity of light in vacuum. At t = t 1 , the relative permittivity is changed to ε r2 = [ε r2x , ε r2z ], and for completeness, the permeability is changed to µ r2y . As mentioned before, this temporal boundary creates a set of FW (E 2 + , H 2 + ) and BW (E 2 − , H 2 − ) waves, with the total electric and magnetic fields defined as E 2 = E 2 + + E 2 − and H 2 = H 2 + -H 2 − , respectively. After applying the temporal boundary conditions for vectors B and D at the temporal boundary (D t1-δ = D t1+δ and B t1-δ = B t1+δ in the limit when δ→0 43 ), it is straightforward 46 to calculate the normalised amplitude of the electric field for both FW and BW waves (for the sake of completeness and easy access, we also give the complete derivation of the fields in Supplementary Materials section 1), resulting in the following expressions: with As observed, the new frequency ω 2 and the amplitude of the generated FW and BW waves clearly depend on the incident angle θ 1 and the values of µ r2y and tensor ε r2 = [ε r2x , ε r2z ], as expected. Note that Eq. (1) reduces to the one shown in ref. 46 when µ r2y = µ r1 . Moreover, in the special case where the change in permittivity/permeability is isotropic (ε r2 = ε r2x = ε r2z , µ r2 = µ r2y ), Eq. (1) is reduced to the case described in refs. 40,41 with (1) denotes the amplitude of the FW and BW waves; however, it is important to remark that each component along the x and z axes (E 2× can be demonstrated that for times t > t 1 + , the momentum k is preserved (θ 1k = θ 2k = θ 1 ), while the direction of the Poynting vector for both FW and BW waves can be calculated as , which can be reduced to the following simple expression: From the expression above, one may notice that the direction of the energy flow is now different from the direction of the phase variation [θ 2S ≠ (θ 2k = θ 1 )], with the former angle depending on the angle of the incident wave before the temporal change (θ 1 ) and the values of the relative permittivity tensor ε r2x and ε r2z . A schematic representation of this performance is shown in Fig. 1e, f, where it is shown how the direction of the energy flow (S) is different from the direction of the wavenumber (k), and that the former can be steered in time by changing the relative permittivity from isotropic to anisotropic tensorial values, reaching the receivers Rx1 or Rx2, depending on the values of ε r2 = [ε r2x , ε r2z ] and the incident angle. In the following sections, we discuss analytical and numerical calculations using this temporal change of ε r for plane waves and Gaussian beams to achieve temporal aiming with timedependent metamaterials.
Oblique incident plane wave: results obtained using isotropic-to-anisotropic ε(t) Let us now evaluate the response of the proposed temporal aiming approach using a time-dependent ε r (t) that is rapidly changed from isotropic to anisotropic values. Without loss of generality, we consider a constant permeability (µ r1 = µ r2 = 1). The isotropic relative permittivity is initially ε r1 = 10 and is modified to ε r2 at t = t 1 = 38 T, where T is the period of the monochromatic incident wave before the change in permittivity occurs. Let us first use two different values for the relative permittivity tensor, i.e., ε r2 = [ε r2x = 8, ε r2z = 12] and ε 2 = [ε r2x = 2, ε r2z = 20], noting that in both cases, ε r2z > ε r2x .
With this set-up, the analytical calculations of the angle of the Poynting vector θ 2S for t > t 1 as a function of the incident angle θ 1 using Eq. (2) are shown as blue and black circles in Fig. 2a. In this figure, it is clear how θ 2S depends on the tensor ε r2 and θ 1 , as explained in the last section. Moreover, note that for the case with ε r2 = [ε r2x = 2, ε r2z = 20], θ 2S is larger for smaller values of θ 1 than in the case with ε r2 = [ε r2x = 8, ε r2z = 12]. For instance, if θ 1 = 15°, then θ 2S will be 69.5°and 21.9°for each case, respectively. These results are expected since the amplitude of the x and z components of the electric field for t > t 1 can be increased or reduced depending on the values of ε r2x and ε r2z (see Supplementary Fig. 2 Analytical results of temporal aiming using plane waves. a Analytically derived angles of the instantaneous Poynting vector θ 2s (t > t 1 ) as a function of the incident angle θ 1 (t < t 1 ) when ε r is changed from isotropic ε r1 = 10 to anisotropic ε r2 = [ε r2x = 8, ε r2z = 12] (blue circles) and to ε r2 = [ε r2x = 2, ε r2z = 20] (black circles). b Snapshot of the H y field (colour plot) and instantaneous Poynting vector (black arrows) distributions for an incident wave with θ 1 = 45°at t < t 1 . c, d Snapshot of the H y field of the FW wave (colour plot) and instantaneous Poynting vector (black arrows) distributions at t > t 1 when ε r is changed from ε r1 = 10 to ε r2 = [ε r2x = 8, ε r2z = 12] and to ε r2 = [ε r2x = 2, ε r2z = 20], respectively, for a time t = 38.2 T > t 1 . e, The same as in panel a but when ε r is changed from ε r1 = 10 to ε r2 = [ε r2x = 12, ε r2z = 8] (blue circles) and to ε r2 = [ε r2x = 20, ε r2z = 2] (black circles). f The same as in panel b but with an incident angle θ 1 = 65°at t < t 1 . g, h The same as in panels c, d but when ε r is changed from ε r1 = 10 to ε r2 = [ε r2x = 12, ε r2z = 8] and to ε r2 = [ε r2x = 20, ε r2z = 2], respectively, for the same incident angle as in panel f, θ 1 = 65°. The time-dependent relative permittivity ε r (t) for the cases under study is shown at the top of panels c, d and g, h magnetic field (H y ) distribution before the change in ε r (t < t 1 ) is shown in Fig. 2b, along with the distribution of the instantaneous Poynting vector (black arrows). Now, at t = t 1 = 38 T, ε is changed from isotropic ε r1 = 10 to ε r2 = [ε r2x = 8, ε r2z = 12]. The analytical H y field distribution for the FW wave for a time t = 38.2 T > t 1 is shown in Fig. 2c. From these results, it can be seen how k is preserved with θ 1k = θ 2k = θ k = θ 1 = 45°while the instantaneous Poynting vector has an angle θ 2S = 56.3°, in agreement with the analytically calculated values from Fig. 2a (extracted from Eq. (2)). As shown in Fig. 2a, if one needs to further increase θ 2S , we can just use a different value of ε r2 . For instance, a snapshot of the analytical H y distribution of the FW wave using ε r2 = [ε r2x = 2, ε r2z = 20] is shown in Fig. 2d at the same time instant as in Fig. 2c. As observed, θ 2S is further increased to 81.4°for the same initial angle θ 1 = 45°. An animation showing this latter scenario can be seen in Supplementary Video 1. For completeness, let us now evaluate the case where ε r2z < ε r2x . Following the same idea as in Fig. 2a-d, ε r is changed from isotropic ε r1 = 10 to ε r2 = [ε r2x = 12, ε r2z = 8] and ε r2 = [ε r2x = 20, ε r2z = 2] at t = t 1 = 38 T. The analytically derived angles of the instantaneous Poynting vector θ 2S are shown in Fig. 2e as blue and black circles, respectively. By comparing these results with those shown in Fig. 2a, we notice that θ 2S is now larger when ε r2 = [ε r2x = 12, ε r2z = 8] than when ε r2 = [ε r2x = 20, ε r2z = 2], again because of the dependence of the x and z components of the electric field and θ 2S on the values of ε r2x and ε r2z . For instance, if we now take the same θ 1 = 15°, θ 2S will be 33.7°and 8.5°for each case. Let us now evaluate the analytically derived field distribution, and let us consider an angle θ 1 = 65°. A snapshot of the incident H y distribution for t < t 1 is shown in Fig. 2f, along with the instantaneous Poynting vector. Now, if the relative permittivity is changed to ε r2 = [ε r2x = 12, ε r2z = 8], the resulting H y distribution at t = 38.2 T > t 1 is the one shown in Fig. 2g, where we have also plotted the instantaneous Poynting vector. As observed, the angle of k is always preserved (θ 1 = 65°), but θ 2S is reduced to θ 2S = 55°, in agreement with the results shown in Fig. 2e. If we now consider ε r2 = [ε r2x = 20, ε r2z = 2], θ 2S is further reduced to θ 2S = 17.8°. An animation showing this latter scenario can be found in Supplementary Video 2. For the sake of completeness, snapshots in time of the BW waves for the cases shown in Fig. 2c, d and Fig. 2g, h are shown in Supplementary Information section 2.
These results demonstrate how the value of θ 2S will differ from the wavenumber direction θ 1k = θ 2k = θ k = θ 1 when ε r is changed from isotropic to anisotropic values. Moreover, this θ 2S can be tuned to smaller/larger angles than θ 1k = θ 2k = θ k = θ 1 by properly engineering the values of the tensor ε r2 = [ε r2x , ε r2z ], a feature that is important for the proposed temporal aiming approach.

Multiple plane waves: Gaussian beam propagation
In the previous section, we evaluated the case where the relative permittivity of the medium is modified from isotropic to anisotropic values for a monochromatic CW plane wave under oblique incidence θ i . However, what would happen if we use more complex waves, such as a Gaussian beam? To answer this question, let us study an oblique incident monochromatic p-polarised Gaussian beam (in-plane E field), as schematically shown in Fig. 3a. As is well known, a Gaussian beam can be modelled as a summation of multiple plane waves travelling in different directions (same magnitude of the wavenumber but different directions for vector k in the same figure) 52 . Moreover, it is also known that the angular aperture of the Gaussian beam will be large/small for small/large values of the beam waist diameter D. If we then use a time-dependent metamaterial for the medium through which the Gaussian beam is travelling, as in the examples from Fig. 2 (i.e., ε r from ε r1 to ε r2 = [ε r2x , ε r2z ]), one will expect to preserve k (in every direction) for t > t 1 , while the angle of the instantaneous Poynting vector θ 2S will be different for each direction of k. This is because each plane wave forming the Gaussian beam will have different θ 1k = θ 2k = θ k = θ 1 (θ 1kn , with n = a, b, c … denoting each of the multiple plane waves forming the Gaussian beam). To visualise this, let us calculate the analytical values of θ 2S as a function θ 1k , as shown in Eq. (2), considering the expression for plane waves. The results are shown as black circles in Fig. 3b when relative ε r is modified from ε r1 = 10 to ε r2 = [ε r2x = 1, ε r2z = 15]. Note that the same trend as in Fig. 2a is observed since ε r2x < ε r2z , showing a clear dependence of θ 2S on θ 1k , as expected. To better understand these results using monochromatic Gaussian beams, let us consider an incident angle θ i = 0°and a beam waist diameter D = 9λ. Here, the permittivity of the whole medium is the same as in Fig. 3b, where it is initially ε r1 = 10 and then it is changed to ε r2 = [ε r2x = 1, ε r2z = 15] at t = t 1 = 30.3 T. With this set-up, the numerical results of the power-flow distribution (i.e., the magnitude of the instantaneous power flow) on the xz plane, along with the instantaneous Poynting vector (blue arrows), for a time t = 30.2 T (t = t 1 − ) are shown in Fig.  3c. As observed, most of the energy is travelling near θ = 0°because of the large beam waist diameter (D = 9λ), as expected. The power-flow distribution and the instantaneous Poynting vector for a time instant just after the change in permittivity to an anisotropic value (t = 30.4 T, t = t 1 + ) are shown in Fig. 3d. By comparing these results with those from Fig. 3c, it can be clearly seen that the Poynting vector angles θ 2S are still close to θ = 0°, but they are indeed modified compared with θ 1kn . These results are in agreement with the analytical angles shown in Fig. 3b, where a value of θ 2S = 0°will be achieved for an angle θ 1k = 0°and can be increased to θ 2S ≈ 15°for a small angle θ 1k = 1°. For completeness, the power-flow distribution and instantaneous Poynting vector for a time t > t 1 are shown in Fig. 3e.
The results discussed in Fig. 3c-e were obtained considering a large beam waist diameter (D = 9λ). However, what if we use a smaller value of D? To evaluate this case, we use D = 2λ, and the results of the power-flow distribution and instantaneous Poynting vector for a time t = 30.2 T (before the change in ε r ) are shown in Fig. 3f. As clearly shown, the angular aperture of the Gaussian beam is increased (as expected), and the energy is now travelling along multiple θ 1kn directions. With this set-up, let us  t 1 ), respectively. f-h Numerical results of the snapshot of the power-flow distribution and instantaneous Poynting vector (blue arrows) distributions at the same times as in panels c-e and using the same time-dependent ε r (t) but considering a Gaussian beam with a beam waist diameter of D = 2λ. In all the numerical results, the incoming signal is switched off once the temporal boundary is induced at t = t 1 to appreciate the effect of using a time-dependent ε r (t) on a signal already present in the medium now change ε r to an anisotropic value, as in Fig. 3d. The results of the power-flow distribution and instantaneous Poynting vector at t = 30.4 T are shown in Fig. 3g. As observed, k is preserved, as expected, but we can now clearly see how multiple directions of θ 2S are obtained, in agreement with the description provided by the plane waves in Fig. 3b. The power-flow distribution and instantaneous Poynting vector for a time t > t 1 , as shown in Fig. 3e, are shown in Fig. 3h. From this figure, one can notice that the magnitude of the power-flow distribution is larger for angles other than 0°. This observation can be explained by looking at the red curve in Fig. 3b, which shows the amplitude of the electric field for the FW wave as a function of the incident angle θ 1k (values calculated using Eq. (1)). As observed, there is a clear dependence of E 2 + /E 1 on θ 1k , achieving an increased/reduced amplitude for large/small values of θ 1k when considering the values under study of ε r2 = [ε r2x = 1, ε r2z = 15]. For instance, E 2 + /E 1 is 0.74 when θ 1k = 0°and increases up to 6.6 when θ 1k = 90°. Similarly, for the BW wave (green circles in Fig. 3b), E 2 − /E 1 is −0.074 when θ 1k = 0°; then, it reaches its maximum of 3.4 when θ 1k = 90°. For completeness, the amplitudes of the FW and BW waves obtained using the temporal permittivity functions studied in Fig. 2 are also shown in Supplementary Fig. S3  The results discussed in Fig. 3 were calculated considering normal incident Gaussian beams at θ i = 0°. For completeness, the numerical results of the power-flow distribution and instantaneous Poynting vector using oblique incident Gaussian beams with θ i = 25°and the same beam waist diameters as in Fig. 3, i.e., D = 9λ and D = 2λ, are shown in Fig. 4a, b and Fig. 4c, d, respectively. As observed for a time t = t 1 -(i.e., when ε r1 = 10), most of the energy travels along the incident angle θ i = 25°for D = 9λ (Fig. 4a), while it spreads to more angles when D = 2λ (Fig. 4c), as expected. Now, when ε is changed to an anisotropic tensor ε r2 = [ε r2x = 1, ε r2z = 15] (Fig. 4b, d), the energy is re-directed to θ 2S , which is no longer parallel to the angle of the wavenumber k, θ 1kn , as explained before. Moreover, note that in the results shown in Fig.  4b, d, the BW waves created at the temporal boundary are ) and after the change in permittivity to ε r2 = [ε r2x = 1, ε r2z = 15] at t = 36.3 T (t > t 1 ), respectively. c, d Numerical results of the snapshot of the power-flow distribution and instantaneous Poynting vector distributions using the same set-up, times and change in ε r as in panels a, b but for a Gaussian beam with a beam waist diameter D = 2λ. As in Fig. 3, in all the numerical results, the incoming signal is switched off once the temporal boundary is induced at t = t 1 to appreciate the effect of using a time-dependent ε r (t) on a signal already present in the medium. Moreover, note that the scale bars for panels b, d were saturated from 0 to 0.8 to better appreciate the FW and BW waves produced at the temporal boundary clearer than the results shown in Fig. 3. This is because of the dependence of the amplitude of the FW and BW waves on the angle θ 1k . Since larger angles θ 1k are generated when tilting the Gaussian beam to an angle θ i = 25°(compared with θ i = 0°in Fig. 3), the amplitudes of both FW and BW waves are further increased (as detailed in Fig. 3b and Eq. (1)).

Temporal aiming narrowband Gaussian wavepacket
In the previous section, we analysed the temporal aiming approach using time-dependent metamaterials with an isotropic-to-anisotropic change in the relative ε r by considering a monochromatic plane wave under oblique incidence and Gaussian beams with different beam waist diameters. In this section, we discuss how temporal aiming can be achieved when the source generates a narrowband Gaussian wavepacket. A schematic representation of this scenario is shown in Fig. 5. Let us consider an oblique incident (θ 1 ) Gaussian wavepacket propagating in a medium with a timedependent ε r [ε r (t), µ r = 1]. Moreover, consider that we have three receivers (Rx1, Rx2 and Rx3) placed at different spatial locations. Without loss of generality, Rx1 is directly aligned with the incoming oblique incident wavepacket, while Rx2 and Rx3 are placed at different locations (see Fig. 5). Now, if we want to send the wavepacket to Rx1, we need only to keep ε(t) = ε 1 = constant since both the source and Rx1 are aligned. However, if the wavepacket is already propagating in the medium, would it be possible to redirect it to reach either Rx2 or Rx3? In the previous section, it was shown how the Poynting vector S is modified to an angle different from the wavenumber k once the relative permittivity is changed from an isotropic value to an anisotropic tensor. Based on this, a way to deflect the propagating wavepacket and reach Rx2, for instance, would be to engineer a time-dependent ε(t), as shown in the inset of Fig. 5 for this receiver. In this case, ε r can be changed from ε r1 to ε r2 = [ε r2x ≠ ε r2z ] at t = t 1 and kept at this value for a certain time duration Δt = t 2 -t 1 . During this time interval, the wavepacket preserves k (as explained in the previous sections), but the direction of the energy flow (S) changes. Hence, the interval Δt should be selected such that the wavepacket moves along the xz plane until it is aligned with the receiver Rx2 at time t = t 2 − . Once this step is achieved, then ε r can be returned from ε r2 = [ε r2x ≠ ε r2z ] to its original value ε r1 at t = t 2 to allow the wavepacket to travel again with the initial angle θ 1 , reaching Rx2. Similarly, this process can be performed to redirect the pulse to Rx3, and the same values of ε r1 and ε r2 = [ε r2x ≠ ε r2z ] can be applied as in the previous case. The only difference would be that the time interval Δt z x Rx 1 Fig. 5 Sketch of the idea of temporal aiming with narrowband Gaussian beams under oblique incidence. The time-dependent ε r functions used to redirect the wavepacket to each receiver (Rx1, Rx2 and Rx3) are shown in the insets should now be modified accordingly to align the wavepacket to this receiver. An example of the temporal aiming described in Fig. 5 is presented in Fig. 6 using an oblique incident narrowband wavepacket with θ 1 = 25°. The specific time-dependent ε r (t) of the medium for re-directing the wavepacket to Rx1, Rx2 and Rx3 is shown in Fig. 6a-c, respectively. For Rx1 (Fig. 6a), the permittivity ε r is constant at ε r1 = 10, with no change because the source is aligned to this receiver, as explained before. The numerical results of the H y field distribution at different times for this case are shown in Fig. 6d-f, where it can be seen how the pulse is directly sent to Rx1. Now, to reach Rx2, the relative ε r is changed from isotropic ε r1 = 10 to anisotropic ε r2 = [ε r2x = 1, ε r2z = 15] (the same values as in Fig. 4) at t = t 1 = 30.3 T and kept at this value until it is returned to ε r1 = 10 at t = t 2 = 33.5 T (Δt = 3.2 T). The numerical results of the H y field distribution for a time t < t 1 , t 1 < t < t 2 and t > t 2 for this case are shown in Fig. 6g-i, respectively. As observed, when t 1 < t < t 2 (Fig. 6h), the wavepacket propagates with an angle defined by the Poynting vector (θ 2S ≈ 82°, in agreement with the analytical values from Eq. (2), which predicts θ 2S = 81.4°). Finally, once ε is returned to the isotropic state with ε r1 = 10, the  Fig. 6 Results of temporal aiming with narrowband Gaussian wavepacket under oblique incidence. a-c Time-dependent ε r to redirect an oblique incident (θ 1 = 25°) narrowband wavepacket towards each receiver Rx1, Rx2 and Rx3, respectively. d-f Snapshot of the H y field distribution at times t < t 1 , t 1 < t < t 2 and t > t 2 , respectively, considering the time-dependent ε r shown in panel a. g-i Snapshot of the H y field distribution at times t < t 1 , t 1 < t < t 2 and t > t 2 , respectively, considering the time-dependent ε shown in panel b. j-l Snapshot of the H y field distribution at times t < t 1 , t 1 < t < t 2 and t > t 2 , respectively, considering the time-dependent ε r shown in panel c. The times at which panels d-l have been obtained are shown as circles in panels a-c to guide the eye wavepacket propagates with the same incident angle (θ 1 = 25°) and is able to reach Rx2 (Fig. 6i). The same process is then applied to the wavepacket to reach Rx3, but now the parameter Δt is increased to Δt = 6.3 T. The ε r function for this receiver and the H y field distribution at different times are shown in Fig. 6c, j-l, demonstrating how the wavepacket can reach Rx3 using this timedependent function of ε r . An animation showing the results of the temporal aiming described in Fig. 6 can be found in Supplementary Video 4. Finally, it is important to note that we change here the relative permittivity of the whole medium through which the wave is travelling. This temporal aiming may be achieved using 2D transmission lines loaded with time-varying circuit elements or with tuneable metasurfaces 36,53 .

Discussion
In conclusion, we have discussed time-dependent metamaterials by inducing temporal boundaries using a rapid change in permittivity from isotropic to anisotropic values. The physics behind this temporal variation of the electromagnetic properties of the medium has been presented, highlighting that the wavenumber and instantaneous Poynting vector of the electromagnetic waves already propagating in this medium exhibit different directions determined by the change in the permittivity tensor. This performance has been exploited to achieve temporal aiming, where an electromagnetic wavepacket can be ushered and re-directed to desired angles by engineering the isotropic-anisotropic temporal variation of the relative permittivity. Different examples have been evaluated both numerically and analytically, such as plane waves under oblique incidence and monochromatic normal and oblique incident Gaussian beams. Moreover, our temporal aiming approach was also evaluated considering the case of an oblique incident narrowband Gaussian wavepacket, showing that it can be possible to deflect and redirect a wavepacket to reach receivers placed at different spatial locations by using isotropic-anisotropicisotropic temporal metamaterials. The ideas presented here may find applications in integrated photonics scenarios where it may be required to redirect and send waves to specific targets/receivers on photonic chips in real time, and may open up new avenues in manipulating waves by ushering and guiding wavepackets at will.

Materials and methods
All numerical simulations were performed using the time-domain solver of the commercial software COMSOL Multiphysics ® . For all the simulations, a rectangular box of dimensions 105λ × 62.5λ was implemented. The incident field was applied from the top boundary of the simulation box via a scattering boundary condition with an out-ofplane magnetic field. The complete non-paraxial Gaussian beam expression was introduced using the angular spectrum technique for plane waves. In this method, the Gaussian beam is calculated as a summation (integral in our case) of multiple plane waves propagating with the same magnitude of wavenumber k but in different directions 52 . Scattering boundary conditions were also implemented on the bottom, left and right boundaries of the simulation box to avoid undesirable reflections. Finally, a triangular mesh was implemented with minimum and maximum sizes of 1.5 × 10 −8 λ and 0.1 λ, respectively, to ensure accurate results. The rapid changes in ε in all the studies were modelled by implementing rectangular analytical functions with smooth transitions using two continuous derivatives to ensure convergence in the calculations.