There is considerable interest in the application of quantum information science to advance computations in plasma physics. A particular point of curiosity is whether it is possible to take advantage of quantum computers to speed up numerical simulations relative to conventional computers. Many of the topics in fusion plasma physics are classical in nature. In order to implement them on quantum computers, it will require couching a classical problem in the language of quantum mechanics. Electromagnetic waves are routinely used in fusion experiments to heat a plasma or to generate currents in the plasma. The propagation of electromagnetic waves is described by Maxwell equations with an appropriate description of the plasma as a dielectric medium. Before advancing to the tensor dielectric of a magnetized plasma, this paper considers electromagnetic wave propagation in a one-dimensional inhomogeneous scalar dielectric. The classic theory of scattering of plane electromagnetic waves at a planar interface, separating two different dielectric media, leads to Fresnel equations for reflection and transmission coefficients. In contrast to plane waves, this paper is on the reflection and transmission of a spatially confined electromagnetic pulse. Following an analytical formulation for the scattering of a Gaussian pulse, it is deduced that the maximum transmission coefficient for a pulse is times that for a plane wave; the incident and transmitted pulses propagate in dielectric media with refractive indices n1 and n2, respectively. The analytical theory is complemented by numerical simulations using a quantum lattice algorithm for Maxwell equations. The algorithm, based on the Riemann–Silberstein–Weber representation of the electromagnetic fields and expressed in terms of qubits, is an interleaved sequence of entangling operators at each lattice site and unitary streaming operators, which transmit information from one site to an adjacent lattice site. Besides substantiating results from the theory for Gaussian pulses, numerical simulations show their validity for non-Gaussian pulses. Apart from their time-asymptotic forms, the simulations display an interplay between the incident, reflected, and transmitted pulses in the vicinity of the transition region between two dielectric media.
The anticipated availability of quantum computers in the near future, and the associated speed up in computations, has provided an impetus for exploring application of quantum information science to plasma physics. Our motivation in this paper is to study the propagation of electromagnetic waves in inhomogeneous dielectric media within the realm of quantum information science. The propagation of waves is described by Maxwell equations in which information about a dielectric medium is expressed through its permittivity. In ordinary materials, the permittivity is usually a scalar, but in a magnetized plasma, it is a tensor. This entire description of waves is classical in the sense that no quantum effects come into play. However, it was recognized early on, in 1931, by Oppenheimer1 that it is possible to cast Maxwell equations for electromagnetic fields in vacuum into a form that is similar to the Dirac equation.2 More recently, this formalism has been extended to waves propagating in a homogeneous dielectric medium and, subsequently, to waves propagating in a spatially inhomogeneous dielectric medium;3 in both cases, the permittivity of the medium is a scalar.
In this paper, we study the propagation of electromagnetic waves in spatially inhomogeneous scalar dielectrics. This is a prelude to future studies on wave propagation in magnetized plasmas. We construct a theoretical model for the reflection and transmission of a spatially confined electromagnetic pulse incident on a surface separating two different dielectric media. The plane wave version of this model is a standard example in textbooks on electromagnetism.4 For a simulation code that will complement the theory, we cast Maxwell equations in a matrix representation that is akin to the Dirac equation.3 We formulate a quantum lattice algorithm (QLA) based on this representation. The results from simulations using this algorithm motivated the analytical formalism and are found to be in excellent agreement with those given by the theory. A QLA is an interleaved sequence of unitary collision and streaming operators that can be modeled by qubit gates. Such an algorithm is ideally parallelized on traditional computers and can be implemented on a quantum computer. Thus, a QLA can be tested before quantum computers become readily available.
We outline the plane wave model of reflection and transmission at an interface separating two different dielectrics. This will set up the theory for propagation of an electromagnetic pulse. In the Cartesian coordinate system, consider a plane interface at z = 0, which divides space into two non-magnetic regions with different scalar dielectric permittivities,
An incident plane wave in region 1, propagating toward the interface along the z-direction—normal to the interface—will lead to reflected and transmitted plane waves. The electric and magnetic fields of the three waves, respectively, are of the following form:
where E0, Er, and Et are the complex field amplitudes and ω is the angular frequency of the wave. Within each dielectric region, the speed of light is , the refractive index is , ki = ωni/c, μ0 is the vacuum permeability, ϵ0 is the vacuum permittivity, c is the speed of light in vacuum, and i = 1, 2. The boundary conditions that follow from Ampere and Faraday equations require continuity of the tangential electric and magnetic fields at z = 0, resulting in the Fresnel relations4
In this paper, we develop a model for the reflection and transmission of a Gaussian pulse at the interface given in (1). The incident Gaussian pulse is propagating along the normal to the interface. We show that the amplitude of the transmitted pulse is different from that obtained for plane waves—it is larger by a factor of . The analytical model is complemented by QLA simulations for electromagnetic wave propagation in a dielectric medium. While the analytical theory is based on satisfying electromagnetic boundary conditions at the discontinuous interface, the QLA simulations propagate the pulse through a continuous, monotonic, narrow interface representing the discontinuity. The two approaches give identical results. Moreover, the simulations display the interplay between the incident, reflected, and transmitted pulses in the vicinity of the interface. This is not possible within the realm of our analytical theory.
The rest of the paper is organized as follows: The analytical model is developed in Sec. II followed by a matrix formulation of the Maxwell equation for an inhomogeneous scalar dielectric in Sec. III. Based on this formulation, the QLA is set up in Sec. IV. In Sec. V, we display results from QLA simulations and compare them with the analytical model.
II. REFLECTION AND TRANSMISSION OF A GAUSSIAN PULSE
In what follows, we assume that the Gaussian pulse has a compact support. However, for mathematical convenience, any integrals involving the pulse will be extended from −∞ to ∞. In addition, the incident, reflected, and transmitted pulses will be considered as separate entities, thereby avoiding those times when the pulses overlap near the interface.
A. The incident Gaussian pulse
At time t = 0, we assume that the incident pulse has the form
where α is a normalization constant, z0 > 0 is the location of the peak of the pulse, σ is its effective width, and the pulse is assumed to be localized in the region z < 0. The requirement
In (8), * denotes the complex conjugate of the field. The time evolution of the normalized Gaussian pulse in medium 1 is
Upon taking the Fourier transform of (10), the plane wave representation of the incident Gaussian pulse is
The magnetic field associated with each plane wave is
B. The reflected and transmitted Gaussian pulses
The reflected (medium 1) and the transmitted (medium 2) pulses are taken to be of the following forms, respectively,
where z1 and z2 are the constants that shift the Gaussian pulses away from z = 0 and σ1 and σ2 are the effective widths of each of the pulses. Analogous to the incident pulse, the plane wave representation of the reflected and transmitted pulses is
where k1 and k2 are the Fourier space variables for the reflected and transmitted waves, respectively. The corresponding magnetic fields are
C. Boundary conditions: Amplitudes and widths of the reflected and transmitted pulses
The boundary conditions—continuity of the tangential electric and magnetic fields—imposed at z = 0 have to be satisfied for all times. Thus,
Following the discussion in Sec. I, the Fresnel jump conditions yield
In defining the wave pulses (10), (14), and (15), it was assumed that the amplitudes E0, Er, and Et are constants independent of space and time. Consequently, in (23) and (24), Er and Et cannot depend on the Fourier variables k, k1, and k2. Hence,
The width of the reflected pulse is the same as that of the incident pulse. For n1 > n2, the width of the transmitted pulse is broader than that of the incident pulse, while for n2 > n1, the transmitted pulse is narrower. There is an intuitive explanation of this result. Let T be the time interval between the leading edge and the trailing edge of the incident pulse arriving at z = 0. The “effective width” of the incident pulse is T/v1. The reflected and transmitted pulses are formed during the time interval T. The effective width of the reflected pulse is T/v1, which is the same as that of the incident pulse. However, the effective width of the transmitted pulse is T/v2—in agreement with (30).
The ratio of the amplitude of the reflected pulse to that of the incident pulse is the same as in (5) for plane wave scattering. However, the ratio of the amplitude of the transmitted pulse to that of the initial pulse is different from (6). It is larger for n2 > n1 and smaller for n2 < n1 by the square-root of the ratio n2/n1. Note that (33) is unchanged when n1 and n2 are interchanged.
III. MAXWELL EQUATIONS: REPRESENTATION IN TERMS OF RIEMANN–SILBERSTEIN–WEBER VECTORS
For non-magnetic materials with a spatially dependent permittivity ϵ(r), Maxwell equations are of the form
After some straightforward algebraic manipulations, and making use of Maxwell equations, we find
If , then, from (37), we note that the equations for F+ and F− decouple and, from (38), the divergences of F+ and of F− are zero. It is easy to see the relationship between Eqs. (37) and (34) and between Eqs. (38) and (35). The advantage of using the RSW vectors is that, in vacuum, the two “polarizations” F+ and F− propagate independently.7
The evolution equations (38) for F+ and F− can each be cast in a 3 × 3 matrix form. In order to include (37) in a unified matrix representation, we need two four-vectors Ψ+ and Ψ− defined as follows:3
In this equation, I4 is a 4 × 4 identity matrix, 04 is a 4 × 4 null matrix, and
and are the Pauli matrices, where
In the expressions for Σ and σ, 02 is a 2 × 2 null matrix. The three Cartesian components of the matrix M are
By the taking the difference of the first and fourth (fifth and eighth) rows of (40) we obtain the evolution equation for () in (38); the sum of the first and fourth (fifth and eighth) rows gives the evolution of (). The sum of the second and third (sixth and seventh) rows leads to the evolution equation for () in (38). The difference between the second and third (sixth and seventh) rows leads to the divergence equation (37). Thus, the compact form (40) accounts for all the four Maxwell equations.
Equation (40) can be separated into three equations corresponding to the principal directions of the Cartesian coordinate system. For a medium in which the permittivity varies along one particular direction, the evolution equation (40) is simplified. However, the three principal directions do not lead to the same evolution equation since the three Pauli matrices are different.8,9
A. Propagation and inhomogeneity along the z-direction
If we define the four elements of Ψ+ and Ψ− as
then (44) leads to
In Eqs. (46) and (47), it is worth noting that the time derivative of each component ψi (i = 0,...,7) is related to the spatial derivative of the same component. Instead, if we had assumed spatial variation in the x or y directions, the time derivative of one component ψi would be related to the spatial derivative of a different component ψj (j ≠ i).8,9 This is a consequence of σz being a diagonal matrix while σx and σy have only off-diagonal, non-zero, elements.
In the sections to follow, we will assume that the speed of light in the medium v(z) is normalized to c. Even though we will continue to use z and t as the space–time variables, it will be implicit that the dimensions of z and t are the same. The time will be related to the discrete temporal step used in the evolution equations for the fields. In this system of units, v(z) = 1/n(z), with n(z) being the “spatially” varying index of refraction.
IV. QUANTUM LATTICE ALGORITHM FOR MAXWELL EQUATIONS
In setting up the quantum lattice algorithm (QLA) for wave propagation in the z-direction, the four-spinor representations in Eqs. (46) and (47) are not in a suitable form. For wave propagation in the x-direction, each spinor ψi was classified as a qubit, and the spatial derivative entangled the qubits.8 In order to include entanglement, we have to increase the dimensionality of the spinor representation from 8 to 16. This follows an earlier precedence where the Gross–Pitaevskii equation for Bose–Einstein condensates10 was formulated at the mesoscopic level by twice as many qubits as field components. At time t = 0, when the electromagnetic fields of an initial incoming pulse are set up, we initialize the qubits in the 16 spinor representation to be
The choice expressed in (48) requires the collision matrix at every spatial lattice point to couple each pair of the qubits: , , , , , , , and . An appropriate form of the unitary collision matrix that couples the qubits is
is the transpose of , and the mixing angle θ will be given later in this section.
There are two different streaming operators that translate qubits from one lattice site to a neighboring site. Each streaming operator is unitary and diagonal and operates only on one element of each pair of qubits discussed above. The streaming operators acting on the 16-qubit representation lead to
where ϵ is the step-size to the adjacent lattice site. We have also used ϵ for the permittivity of a medium. In the rest of the narrative, ϵ is the step-size, and any dielectric medium will be described by its index of refraction.
In order to include spatial inhomogeneity in the refractive index, we need two collision operators, which provide the 16-qubit coupling, similar in vein to the eight-spinor coupling in the second term on the right hand side of Eqs. (46) and (47). These operators, P(1) and P(2), referred to as potential collision operators, are
where γ is a mixing angle.
From the collide-stream operators, we can construct the following unitary operators:
where † is the complex conjugate transpose of the matrix. The time evolution of the 16 qubits can be expressed in terms of these unitary operators,
The quantum lattice algorithm is complete once the mixing angles θ and γ are defined. For this, we assume that ϵ is a perturbation parameter and make the ansatz that θ ∝ ϵ and γ ∝ ϵ2. This ordering is akin to diffusion ordering. Subsequently, we expand the evolution equation (57) out to order ϵ2 using Mathematica. There is an additional constraint that guides our choice of the mixing angles. At order ϵ2, with appropriate combinations of the 16 qubits, we should retrieve, in the continuum limit, the eight-spinor form of Maxwell equations given in (46) and (47). We find that
where ′ denotes a derivative with respect to the argument. Note that, in our normalized set of units, n(z) = 1/v(z). In the continuum limit, and to order ϵ2, (57) leads to the following mesoscopic evolution equation for the 16 qubits:
V. QLA SIMULATIONS FOR TIME EVOLUTION OF AN ELECTROMAGNETIC PULSE
The analytical expressions for the amplitude and width of the reflected and transmitted pulses in Sec. II are obtained by imposing electromagnetic boundary conditions at the interface. In essence, the scattering is treated as a boundary value problem, and the results for the relative amplitudes are in the time-asymptotic limit. In contrast, the QLA simulations treat the scattering as an initial value problem. At time t = 0, a well-defined pulse is initiated in the vicinity of the boundary of a simulation domain. The pulse propagates toward the dielectric interface following the QLA discussed in Sec. IV. The interface is set up to be a smooth, monotonic, narrow transition layer. The time evolution of the pulse, as it propagates through the transition layer, is governed by the same QLA—no boundary conditions are imposed anywhere inside the simulation domain. In the “time-asymptotic” limit, when the reflected and transmitted pulses are well separated, we measure their amplitudes and compare with theory. At the edge of the simulation domain, we impose periodic boundary conditions. However, we stop our simulations well before the pulses reach the domain boundaries. In the QLA computations discussed below, we have set ϵ = 0.2.
For the QLA simulations, we assume a dielectric slab with n2 = 1.5 surrounded by vacuum with n1 = 1. Figure 1(a) shows the index of refraction along the z-axis, while Fig. 1(b) is a magnified view in the vicinity of the transition layer. In the discussion to follow, we refer to the side of the dielectric slab facing the incoming pulse as the front-end of the slab, while the opposite side is the back-end.
We will present results from QLA simulations for two different profiles of the initial electromagnetic pulse. In the normalized system of units, the first profile for the x-component of the electric field is a hyperbolic secant function,
while the second profile is an exponential cusp of the form
In vacuum, the profile for By is the same as that for Ex. In the simulations, we will set E0 = 0.01.
We carried out a set of simulations for a Gaussian pulse, corresponding to the theory in Sec. II. The results, as it turns out, are the same as those for the hyperbolic secant pulse. Consequently, we will concentrate our discussion on the hyperbolic secant pulse.
A. Time evolution of a hyperbolic secant pulse
For an initial profile of the form in (61), the pulse at time t = 6500 is plotted in Fig. 2(a) as it approaches the front-end of the dielectric slab. The Ex and By profiles are on top of each other in vacuum. At time t = 12 000, the initial pulse does not exist as it has completely crossed the front-end of the slab. Only the reflected and transmitted pulses exist with their respective Ex (blue) and By (red), distinctly visible in Fig. 2(b). The reflected pulse is propagating along the −z-direction, while the transmitted pulse is inside the dielectric slab and propagating along the z-direction. From Faraday’s law in (3), we expect the reflected Ex and By to be out of phase. The simulation results in Fig. 2(b) are in agreement. Furthermore, from (32), since n1 < n2, it is the reflected electric field that will flip sign with respect to the electric field of the incoming pulse. For the transmitted pulse, from (4), we note that By/Ex = n2 = 1.5 and the two field components are in phase. The simulation results are in accordance with these theoretical expectations. From (32) and (33), we see that the maximum amplitudes of the reflected and transmitted pulses are and , respectively. Subscript 1 indicates the first encounter with the vacuum–dielectric interface.
As the simulation continues to advance in time, the transmitted pulse reaches the back-end of the dielectric slab. A part of this pulse will get transmitted out into vacuum, and a part of it will get reflected and remain inside the dielectric slab. Figure 3(a) shows this stage. The pulse transmitted through the back-end of the slab will have Ex = By from Faraday’s law. From (33), the transmitted electric field has the same phase as the wave incident on the back-end. Moreover, the ratio of the electric field of the transmitted pulse to the incident pulse is Et2/Et1 = 0.9798. For the pulse reflected from the back-end of the slab, (32) gives Er2/Et1 = 0.2 with the two electric fields being in phase. Thus, the reflected and incident electric fields have the same phase after scattering from the back-end of the dielectric slab. In order to satisfy Faraday’s law, for the reflected pulse, the magnetic field has to be out of phase with the electric field. The results in Fig. 3(b) are in accordance with these properties of the pulses.
For longer times, the pulse reflected from the back-end reaches the front-end dielectric boundary and undergoes a reflection and transmission. Figure 3(b) shows the various pulses that are propagating in vacuum and in the dielectric slab. It is interesting to note that the electric and magnetic fields of the two pulses propagating in vacuum to the left of the dielectric slab are out of phase with respect to each other. The phase difference can be explained using the same arguments, based on Faraday’s law and (33), as discussed above. Furthermore, for the wave transmitted through the front-end, Et3/Er2 = 0.9798.
From our theory, we know that the width of the reflected pulse is the same as the width of the incident pulse. However, the width of the transmitted pulse either narrows or broadens depending on whether the incident pulse is propagating from a medium of a lower refractive index to a medium of a higher refractive index or vice versa. This is visually discernible in Figs. 2(b), 3(a), and 3(b). An important conclusion from the analytical analysis is that the maximum amplitude of the transmitted pulse is different from that obtained for a single plane wave model. The maximum amplitude of the transmitted pulse changes by a factor of for propagation from a medium with a refractive index n1 to a medium with a refractive index n2. These results are borne out as shown in Table I. The peak amplitudes of the pulses at various stages of propagation are in excellent agreement with theory.
|Time .||z .||Pulse .||Bmax .||Emax .||Et/Ei .||Er/Ei .||Bmax/Emax .|
|Time .||z .||Pulse .||Bmax .||Emax .||Et/Ei .||Er/Ei .||Bmax/Emax .|
B. Poynting flux
The instantaneous Poynting flux is defined as
where is the outward pointing normal at the vacuum–dielectric interface. In Fig. 4, we plot S(t) as a function of time. The QLA conserves the Poynting flux reasonably well over the entire time of the simulation, except for certain gaps in time. In the first gap around t = 9000, the initial pulse is in the vicinity of the front-end boundary of the dielectric where the incident and reflected pulses overlap and are not clearly separated. The second gap is around the time when the initial transmitted pulse reaches the back-end boundary of the dielectric. It is the same with the other gaps when the pulse inside the dielectric reaches one boundary or the other. In these cases, it is not straightforward to explicitly isolate the incident and reflected pulses and easily assign the outward pointing normal to each pulse as required in (63). In the Appendix, we discuss one possible method for working around this dilemma. Regardless, it should be noted that S(t) is well conserved away from these isolated instances of time.
C. Time evolution of an exponential cusp pulse
In this section, we consider the propagation of an exponential cusp profile (62) through a similar dielectric slab shown in Fig. 1(a). Since this pulse has a longer tail than the hyperbolic secant pulse, we increase the spatial domain for the simulations as well as the spatial extent of the dielectric slab: 12 000 ≲ z ≲ 15 000. The initial profile for the fields is shown in Fig. 5(a). The subsequent propagation and splitting of the initial pulse into reflected and transmitted pulses are shown in Figs. 5(b), 6(a), and 6(b). The similarity with the propagation of a hyperbolic secant pulse, discussed in Subsection V A, is quite obvious. Consequently, all analyses of the evolution and splitting of pulses is essentially the same for the two profiles.
The numerical results shown in Table II not only agree with the analytical calculations but also bear remarkable resemblance to those in Table I. The theoretical formulation was for an electromagnetic Gaussian pulse. However, the QLA simulations for three different pulse shapes, including a Gaussian, have been in remarkable agreement with the theory. Consequently, we believe that the physics of scattering by a dielectric interface is common for different pulse shapes that are initially spatially confined.
|Time .||z .||Pulse .||Bmax .||Emax .||Et/Ei .||Er/Ei .||Bmax/Emax .|
|Time .||z .||Pulse .||Bmax .||Emax .||Et/Ei .||Er/Ei .||Bmax/Emax .|
The instantaneous Poynting flux S(t) is plotted as a function of time in Fig. 7. This plot is similar to that for the hyperbolic secant profile shown in Fig. 4. The explanation put forth for the hyperbolic secant profile applies here as well.
VI. SUMMARY AND CONCLUSION
We have shown, analytically and computationally, that the scattering of a spatially confined electromagnetic pulse by an interface, separating two disparate dielectric media, is different from the scattering of a plane wave. In particular, the transmission coefficient is modified by a factor for a pulse traveling from a medium with a refractive index n1 to a medium with a refractive index n2. The analytical model is based on a Fourier expansion of a Gaussian pulse and the matching of electromagnetic boundary conditions at the interface separating the two media. The computational results are obtained from a code that has several layers of formalism embedded into it. We start off by expressing Maxwell equations in a matrix form using the Reimann–Silberstein–Weber vectors.3 This is an eight-spinor representation of Maxwell equations, and this has similarities to the Dirac equation for a massless particle. It also forms a basis for the quantum lattice algorithm that solves Maxwell equations. Since each spinor can be cast as a qubit, the initial expectation is that we need an eight-qubit algorithm. However, for wave propagation along the z-direction, the Pauli spin matrix σz is diagonal and does not entangle the qubits. As a result, we developed a 16-qubit algorithm that allows for entanglement of qubits. The subsequent QLA is a series of streaming and collision operators that advance the 16 qubits from one lattice site to another and entangle them at each site. The QLA recovers the full set of four Maxwell equations when expanded to second order in ϵ—the separation between adjacent lattice sites.
The 16-qubit representation for wave propagation in the z-direction can be contrasted with the eight-qubit representation of Maxwell equations in the x-direction.8 The differences arise from the Pauli matrices σx and σz—the former has non-zero values for the off-diagonal elements while the latter has non-zero values for the diagonal elements. This leads to very different collision matrices. Since both representations are for the same set of Maxwell equations, we should not expect any differences in the results for wave propagation in dielectric media. Indeed, this is the case.
The maximum amplitudes of the scattered waves obtained from QLA simulations are in excellent agreement with those given by the analytical model. The QLA simulations for two different initial pulses lead to the same results for the maximum amplitudes of the reflected and transmitted pulses. We did not display any results for a Gaussian pulse of the type used in the theoretical model since the QLA simulations yielded the same ratios for the amplitudes, as shown in Table I. Hence, the theoretical results are applicable to the scattering of different pulse shapes of finite spatial width.
Finally, we like to point out that the simulation results do not change as the parameter ϵ is varied. We used ϵ = 0.2 for the results displayed in this paper. Even for values of ϵ approaching unity, we still recovered the factor associated with the transmission coefficient.
This research was supported by the Department of Energy under Grant Nos. DE-SC0021647, DE-FG02-91ER-54109, DE-SC0021651, DE-SC0021857, and DE-SC0021653.
Conflict of Interest
The authors have no conflicts to disclose.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: DETERMINATION OF THE POYNTING FLUX WHEN THE PULSES OVERLAP
The dips in the instantaneous Poynting flux in Figs. 4 and 7 occur when the incident and reflected pulses overlap in the vacuum–dielectric transition region. In such instances, the assignment of an outward pointing normal for the appropriate fields is tricky. In this section, we resolve this issue for a Gaussian pulse.
A new set of QLA simulations is performed in which the initial pulse has a Gaussian profile. The pulse propagates from the vacuum toward a dielectric medium with a refractive index n2 = 2. The vacuum dielectric boundary layer shown in Fig. 8 has a transition layer that is 12 lattice units wide. Figures 9(a) and 9(b) display, in blue, the electric field Ex and the magnetic field By, respectively, of the pulse at t = 3500. At this time, the incident and reflected pulses overlap in the vicinity of the transition layer. We subsequently carry out QLA simulations in vacuum and superimpose the associated Ex and By fields, in red, in Figs. 9(a) and 9(b), respectively.
In order to correctly determine the reflected part of the pulse and associate it properly with the outward pointing normal, we subtract, in Figs. 9(a) and 9(b), the fields in blue from the vacuum fields in red for z < 3000. We know that for z > 3000, there is only the transmitted pulse, while for z < 3000, the initial and the reflected pulses coexist. The subtraction separates out the reflected pulse from the incident pulse; the incident pulse propagates along the z-direction while the reflected pulse propagates along the −z-direction. This procedure helps remove any ambiguity and allows for correct evaluation of S(t). It has to be carried out for all the time steps during which the incident and reflected pulses overlap. In Fig. 10, we compare the results from this procedure with the one used to evaluate S(t) in Figs. 4 and 7. Figure 10 illustrates that the modified evaluation of S(t) conserves the flow of energy when we properly account for the reflected part of the pulse in the overlap region. The dips in Figs. 4 and 7 are just due to the approximate way in which we try to separate the incident and reflected pulses.
Earlier in this paper, it was mentioned that, through QLA simulations, we can visualize the interplay between incident, reflected, and transmitted pulses in the neighborhood of the transition region between two dielectric media. The blue curves in Figs. 9(a) and 9(b) illustrate that point. This would be difficult to realize with the theoretical model since it does not solve an initial value problem.