Reflection and transmission of electromagnetic pulses at a planar dielectric interface: Theory and quantum lattice simulations

Ram, Abhay K.; Vahala, George; Vahala, Linda; Soe, Min

doi:10.1063/5.0067204

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 $\sqrt{n_{2} / n_{1}}$ times that for a plane wave; the incident and transmitted pulses propagate in dielectric media with refractive indices n₁ and n₂, 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.

I. INTRODUCTION

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 Oppenheimer¹ that it is possible to cast Maxwell equations for electromagnetic fields in vacuum into a form that is similar to the Dirac equation.² 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;³ 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.⁴ For a simulation code that will complement the theory, we cast Maxwell equations in a matrix representation that is akin to the Dirac equation.³ 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,

ϵ (z) = \{\begin{cases} ϵ_{1} & f o r z < 0 (r e g i o n 1), \\ ϵ_{2} & f o r z > 0 (r e g i o n 2) . \end{cases})

(1)

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:

E_{I} = E_{0} e^{i (k_{1} z - ω t)} \hat{x}, B_{I} = \frac{1}{c_{1}} E_{0} e^{i (k_{1} z - ω t)} \hat{y},

(2)

E_{R} = E_{r} e^{i (- k_{1} z - ω t)} \hat{x}, B_{R} = - \frac{1}{c_{1}} E_{r} e^{i (- k_{1} z - ω t)} \hat{y},

(3)

E_{T} = E_{t} e^{i (k_{2} z - ω t)} \hat{x}, B_{T} = \frac{1}{c_{2}} E_{t} e^{i (k_{2} z - ω t)} \hat{y},

(4)

where E₀, E_r, and E_t are the complex field amplitudes and ω is the angular frequency of the wave. Within each dielectric region, the speed of light is $v_{i} = 1 / \sqrt{μ_{o} ϵ_{i}}$ ⁠, the refractive index is $n_{i} = \sqrt{ϵ_{i} / ϵ_{0}} = c / v_{i}$ ⁠, k_i = ωn_i/c, μ₀ is the vacuum permeability, ϵ₀ 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 relations⁴

E_{r} = - (\frac{v_{1} - v_{2}}{v_{1} + v_{2}}) E_{0} = (\frac{n_{1} - n_{2}}{n_{1} + n_{2}}) E_{0},

(5)

E_{t} = (\frac{2 v_{2}}{v_{1} + v_{2}}) E_{0} = (\frac{2 n_{1}}{n_{1} + n_{2}}) E_{0} .

(6)

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 $\sqrt{n_{2} / n_{1}}$ ⁠. 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

E_{I} (z) = E_{0} α e^{- {(z + z_{0})}^{2} / (2 σ^{2})} \hat{x},

(7)

where α is a normalization constant, z₀ > 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

\int_{- \infty}^{\infty} d z E_{I} (z) \cdot E_{I}^{*} (z) = {|E_{0}|}^{2}

(8)

leads to

α = \frac{1}{π^{1 / 4} \sqrt{σ}} .

(9)

In (8), * denotes the complex conjugate of the field. The time evolution of the normalized Gaussian pulse in medium 1 is

E_{I} (z, t) = \frac{E_{0}}{π^{1 / 4} \sqrt{σ}} e^{- {(z - v_{1} t + z_{0})}^{2} / (2 σ^{2})} \hat{x} f o r z < 0 .

(10)

Upon taking the Fourier transform of (10), the plane wave representation of the incident Gaussian pulse is

\begin{align} E_{I k} (k, z, t) & = \frac{E_{0}}{π^{1 / 4}} \sqrt{2 π σ} e^{- σ^{2} k^{2} / 2} e^{i k (z + z_{0} - v_{1} t)} \hat{x} \\ \equiv E_{I k} (k, z, t) \hat{x}, \end{align}

(11)

with

E_{I} (z, t) = \frac{1}{2 π} \int_{- \infty}^{\infty} d k E_{I k} (k, z, t) .

(12)

The magnetic field associated with each plane wave is

B_{I k} (k, z, t) = \frac{1}{v_{1}} E_{I k} (k, z, t) \hat{y} .

(13)

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,

E_{R} (z, t) = \frac{E_{r}}{π^{1 / 4} \sqrt{σ_{1}}} e^{- {(z + z_{1} + v_{1} t)}^{2} / (2 σ_{1}^{2})} \hat{x} f o r z < 0,

(14)

E_{T} (z, t) = \frac{E_{t}}{π^{1 / 4} \sqrt{σ_{2}}} e^{- {(z - z_{2} - v_{2} t)}^{2} / (2 σ_{2}^{2})} \hat{x} f o r z > 0,

(15)

where z₁ and z₂ are the constants that shift the Gaussian pulses away from z = 0 and σ₁ and σ₂ 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

\begin{align} E_{R k} (k, z, t) & = \frac{E_{r}}{π^{1 / 4}} \sqrt{2 π σ_{1}} e^{- σ_{1}^{2} k_{1}^{2} / 2} e^{i k_{1} (z + z_{1} + v_{1} t)} \hat{x} \\ \equiv E_{R k} (k, z, t) \hat{x}, \end{align}

(16)

\begin{align} E_{T k} (k, z, t) & = \frac{E_{t}}{π^{1 / 4}} \sqrt{2 π σ_{2}} e^{- σ_{2}^{2} k_{2}^{2} / 2} e^{i k_{2} (z - z_{2} - v_{2} t)} \hat{x} \\ \equiv E_{T k} (k, z, t) \hat{x}, \end{align}

(17)

where k₁ and k₂ are the Fourier space variables for the reflected and transmitted waves, respectively. The corresponding magnetic fields are

B_{R k} (k, z, t) = - \frac{1}{v_{1}} E_{R k} (k, z, t) \hat{y},

(18)

B_{T k} (k, z, t) = \frac{1}{v_{2}} E_{T k} (k, z, t) \hat{y} .

(19)

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,

ω \equiv k v_{1} = - k_{1} v_{1} = k_{2} v_{2} .

(20)

Following the discussion in Sec. I, the Fresnel jump conditions yield

\frac{E_{r}}{π^{1 / 4}} \sqrt{2 π σ_{1}} e^{- σ_{1}^{2} k_{1}^{2} / 2} e^{i k_{1} z_{1}} = \frac{n_{1} - n_{2}}{n_{1} + n_{2}} \frac{E_{0}}{π^{1 / 4}} \sqrt{2 π σ} e^{- σ^{2} k^{2} / 2} e^{i k z_{0}},

(21)

\frac{E_{t}}{π^{1 / 4}} \sqrt{2 π σ_{2}} e^{- σ_{2}^{2} k_{2}^{2} / 2} e^{- i k_{2} z_{2}} = \frac{2 n_{1}}{n_{1} + n_{2}} \frac{E_{0}}{π^{1 / 4}} \sqrt{2 π σ} e^{- σ^{2} k^{2} / 2} e^{i k z_{0}} .

(22)

From Eqs. (21) and (22),

E_{r} = \frac{n_{1} - n_{2}}{n_{1} + n_{2}} E_{0} \sqrt{\frac{σ}{σ_{1}}} e^{- \frac{1}{2} (σ^{2} k^{2} - σ_{1}^{2} k_{1}^{2})} e^{i (k z_{0} - k_{1} z_{1})},

(23)

E_{t} = \frac{2 n_{1}}{n_{1} + n_{2}} E_{0} \sqrt{\frac{σ}{σ_{2}}} e^{- \frac{1}{2} (σ^{2} k^{2} - σ_{2}^{2} k_{2}^{2})} e^{i (k z_{0} + k_{2} z_{2})} .

(24)

In defining the wave pulses (10), (14), and (15), it was assumed that the amplitudes E₀, E_r, and E_t are constants independent of space and time. Consequently, in (23) and (24), E_r and E_t cannot depend on the Fourier variables k, k₁, and k₂. Hence,

σ^{2} k^{2} - σ_{1}^{2} k_{1}^{2} = 0,

(25)

σ^{2} k^{2} - σ_{2}^{2} k_{2}^{2} = 0,

(26)

k z_{0} - k_{1} z_{1} = 0,

(27)

k z_{0} + k_{2} z_{2} = 0 .

(28)

From (20), k₁ = −k, and $k_{2} = k (v_{1} / v_{2})$ ⁠. Making use of these equalities in (25) and (26), respectively, leads to

σ_{1} = |\frac{k}{k_{1}}| σ = σ,

(29)

σ_{2} = |\frac{k}{k_{2}}| σ = \frac{v_{2}}{v_{1}} σ = \frac{n_{1}}{n_{2}} σ .

(30)

The width of the reflected pulse is the same as that of the incident pulse. For n₁ > n₂, the width of the transmitted pulse is broader than that of the incident pulse, while for n₂ > n₁, 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/v₁. The reflected and transmitted pulses are formed during the time interval T. The effective width of the reflected pulse is T/v₁, which is the same as that of the incident pulse. However, the effective width of the transmitted pulse is T/v₂—in agreement with (30).

Equations (27) and (28) give a relationship between the shifts,

z_{1} = - z_{0} a n d z_{2} = - \frac{v_{2}}{v_{1}} z_{0} = - \frac{n_{1}}{n_{2}} z_{0} .

(31)

Finally, Eqs. (23) and (24) lead to the following conclusion:

E_{r} = \frac{n_{1} - n_{2}}{n_{1} + n_{2}} E_{0},

(32)

E_{t} = \frac{2 n_{1}}{n_{1} + n_{2}} \sqrt{\frac{n_{2}}{n_{1}}} E_{0} .

(33)

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 n₂ > n₁ and smaller for n₂ < n₁ by the square-root of the ratio n₂/n₁. Note that (33) is unchanged when n₁ and n₂ 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

\nabla \cdot \{ϵ (r) E (r, t)\} = 0, \nabla \cdot B (r, t) = 0,

(34)

\nabla \times E (r, t) = - \frac{\partial B (r, t)}{\partial t}, \nabla \times B (r, t) = μ_{0} ϵ (r) \frac{\partial E (r, t)}{\partial t} .

(35)

Here, we have assumed that there are no free charges or currents in the materials. The Riemann–Silberstein–Weber (RSW) vectors^5,6 for a dielectric medium are defined as³

F^{\pm} (r, t) = \frac{1}{\sqrt{2}} [\sqrt{ϵ (r)} E (r, t) \pm \frac{i}{\sqrt{μ_{0}}} B (r, t)] .

(36)

After some straightforward algebraic manipulations, and making use of Maxwell equations, we find

\nabla \cdot F^{\pm} (r, t) = \frac{1}{2 v (r)} \nabla v (r) \cdot [F^{\pm} (r, t) + F^{\mp} (r, t)],

(37)

i \frac{\partial F^{\pm} (r, t)}{\partial t} = \pm v (r) \nabla \times F^{\pm} \pm \frac{1}{2} \nabla v (r) \times [F^{\pm} (r, t) + F^{\mp} (r, t)],

(38)

where $v (r) = 1 / \sqrt{μ_{0} ϵ (r)}$ has the dimensions of speed. In this representation, Eq. (38) is the evolution equation, with (37) being the constraint.

If $ϵ (r) = ϵ_{0}$ ⁠, 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.⁷

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:³

Ψ^{+} (r, t) = [\begin{matrix} - F_{x}^{+} + i F_{y}^{+} \\ F_{z}^{+} \\ F_{z}^{+} \\ F_{x}^{+} + i F_{y}^{+} \end{matrix}], Ψ^{-} (r, t) = [\begin{matrix} - F_{x}^{-} - i F_{y}^{-} \\ F_{z}^{-} \\ F_{z}^{-} \\ F_{x}^{-} - i F_{y}^{-} \end{matrix}],

(39)

where $F^{\pm} = F_{x}^{\pm} \hat{x} + F_{y}^{\pm} \hat{y} + F_{z}^{\pm} \hat{z}$ ⁠. The evolution equations for Ψ⁺ and Ψ⁻ that follow from (37) and (38) are of the form³

\frac{\partial}{\partial t} [\begin{matrix} \begin{matrix} I_{4} & 0_{4} \\ 0_{4} & I_{4} \end{matrix} \end{matrix}] [\begin{matrix} Ψ^{+} \\ Ψ^{-} \end{matrix}] = - v (r) [\begin{matrix} M \cdot \nabla + Σ \cdot \frac{1}{2} \nabla \{\ln v (r)\} & - i M_{z} [Σ \cdot \frac{1}{2} \nabla \{\ln v (r)\}] α_{y} \\ - i M_{z} [Σ^{*} \cdot \frac{1}{2} \nabla \{\ln v (r)\}] α_{y} & M^{*} \cdot \nabla + Σ^{*} \cdot \frac{1}{2} \nabla \{\ln v (r)\} \end{matrix}] [\begin{matrix} Ψ^{+} \\ Ψ^{-} \end{matrix}] .

(40)

In this equation, I₄ is a 4 × 4 identity matrix, 0₄ is a 4 × 4 null matrix, and

Σ = [\begin{matrix} \begin{matrix} σ & 0_{2} \\ 0_{2} & σ \end{matrix} \end{matrix}], α = [\begin{matrix} \begin{matrix} 0_{2} & σ \\ σ & 0_{2} \end{matrix} \end{matrix}],

(41)

and $σ = σ_{x} \hat{x} + σ_{y} \hat{y} + σ_{z} \hat{z}$ are the Pauli matrices, where

σ_{x} = [\begin{matrix} \begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix} \end{matrix}], σ_{y} = [\begin{matrix} \begin{matrix} 0 & - i \\ i & 0 \end{matrix} \end{matrix}], σ_{z} = [\begin{matrix} \begin{matrix} 1 & 0 \\ 0 & - 1 \end{matrix} \end{matrix}] .

(42)

In the expressions for Σ and σ, 0₂ is a 2 × 2 null matrix. The three Cartesian components of the matrix M are

\begin{gathered} M_{x} = [\begin{matrix} \begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{matrix} \end{matrix}], M_{y} = [\begin{matrix} \begin{matrix} 0 & 0 & - i & 0 \\ 0 & 0 & 0 & - i \\ i & 0 & 0 & 0 \\ 0 & i & 0 & 0 \end{matrix} \end{matrix}], \\ M_{z} = [\begin{matrix} \begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & - 1 \end{matrix} \end{matrix}] . \end{gathered}

(43)

By the taking the difference of the first and fourth (fifth and eighth) rows of (40) we obtain the evolution equation for $F_{x}^{+}$ (⁠ $F_{x}^{-}$ ⁠) in (38); the sum of the first and fourth (fifth and eighth) rows gives the evolution of $F_{y}^{+}$ (⁠ $F_{y}^{-}$ ⁠). The sum of the second and third (sixth and seventh) rows leads to the evolution equation for $F_{z}^{+}$ (⁠ $F_{z}^{-}$ ⁠) 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

As in Sec. II, we consider one-dimensional (1D) propagation of electromagnetic waves along the z-direction in a medium with permittivity that is a function of z only. Then, Eq. (40) reduces to

\begin{align} \frac{\partial}{\partial t} [\begin{matrix} Ψ^{+} \\ Ψ^{-} \end{matrix}] & = - v (z) \frac{\partial}{\partial z} [\begin{matrix} M_{z} & 0_{4} \\ 0_{4} & M_{z}^{*} \end{matrix}] [\begin{matrix} Ψ^{+} \\ Ψ^{-} \end{matrix}] \\ - \frac{1}{2} \frac{\partial v (z)}{\partial z} [\begin{matrix} Σ_{z} & - i M_{z} Σ_{z} α_{y} \\ - i M_{z} Σ_{z}^{*} α_{y} & Σ_{z}^{*} \end{matrix}] [\begin{matrix} Ψ^{+} \\ Ψ^{-} \end{matrix}] . \end{align}

(44)

If we define the four elements of Ψ⁺ and Ψ⁻ as

Ψ^{+} (z, t) = [\begin{matrix} ψ_{0} (z, t) \\ ψ_{1} (z, t) \\ ψ_{2} (z, t) \\ ψ_{3} (z, t) \end{matrix}], Ψ^{-} (z, t) = [\begin{matrix} ψ_{4} (z, t) \\ ψ_{5} (z, t) \\ ψ_{6} (z, t) \\ ψ_{7} (z, t) \end{matrix}],

(45)

then (44) leads to

\begin{align} \frac{\partial}{\partial t} [\begin{matrix} \begin{matrix} ψ_{0} \\ ψ_{1} \\ ψ_{2} \\ ψ_{3} \end{matrix} \end{matrix}] & = - v (z) \frac{\partial}{\partial z} [\begin{matrix} \begin{matrix} ψ_{0} \\ ψ_{1} \\ - ψ_{2} \\ - ψ_{3} \end{matrix} \end{matrix}] - \frac{1}{2} \frac{\partial v (z)}{\partial z} [\begin{matrix} \begin{matrix} ψ_{0} - ψ_{7} \\ - ψ_{1} - ψ_{6} \\ ψ_{2} + ψ_{5} \\ - ψ_{3} + ψ_{4} \end{matrix} \end{matrix}], \end{align}

(46)

\begin{align} \frac{\partial}{\partial t} [\begin{matrix} ψ_{4} \\ ψ_{5} \\ ψ_{6} \\ ψ_{7} \end{matrix}] & = - v (z) \frac{\partial}{\partial z} [\begin{matrix} \begin{matrix} ψ_{4} \\ ψ_{5} \\ - ψ_{6} \\ - ψ_{7} \end{matrix} \end{matrix}] - \frac{1}{2} \frac{\partial v (z)}{\partial z} [\begin{matrix} \begin{matrix} ψ_{4} - ψ_{3} \\ - ψ_{5} - ψ_{2} \\ ψ_{6} + ψ_{1} \\ - ψ_{7} + ψ_{0} \end{matrix} \end{matrix}] . \end{align}

(47)

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.⁸ 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 condensates¹⁰ 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

\begin{aligned} q_{0} (z, 0) & = \frac{1}{2} ψ_{0} (z, 0), & q_{1} (z, 0) & = \frac{1}{2} ψ_{1} (z, 0), & q_{2} (z, 0) & = q_{0} (z, 0), & q_{3} (z, 0) & = q_{1} (z, 0), \\ q_{4} (z, 0) & = \frac{1}{2} ψ_{2} (z, 0), & q_{5} (z, 0) & = \frac{1}{2} ψ_{3} (z, 0), & q_{6} (z, 0) & = q_{4} (z, 0), & q_{7} (z, 0) & = q_{5} (z, 0), \\ q_{8} (z, 0) & = \frac{1}{2} ψ_{4} (z, 0), & q_{9} (z, 0) & = \frac{1}{2} ψ_{5} (z, 0), & q_{10} (z, 0) & = q_{8} (z, 0), & q_{11} (z, 0) & = q_{9} (z, 0), \\ q_{12} (z, 0) & = \frac{1}{2} ψ_{6} (z, 0), & q_{13} (z, 0) & = \frac{1}{2} ψ_{7} (z, 0), & q_{14} (z, 0) & = q_{12} (z, 0), & q_{15} (z, 0) & = q_{13} (z, 0) . \end{aligned}

(48)

The choice expressed in (48) requires the collision matrix at every spatial lattice point to couple each pair of the qubits: $(q_{0}, q_{2})$ ⁠, $(q_{1}, q_{3})$ ⁠, $(q_{4}, q_{6})$ ⁠, $(q_{5}, q_{7})$ ⁠, $(q_{8}, q_{10})$ ⁠, $(q_{9}, q_{11})$ ⁠, $(q_{12}, q_{14})$ ⁠, and $(q_{13}, q_{15})$ ⁠. An appropriate form of the unitary collision matrix that couples the qubits is

C (θ) = [\begin{matrix} V_{4} (θ) & 0_{4} & 0_{4} & 0_{4} \\ 0_{4} & V_{4}^{T} (θ) & 0_{4} & 0_{4} \\ 0_{4} & 0_{4} & V_{4} (θ) & 0_{4} \\ 0_{4} & 0_{4} & 0_{4} & V_{4}^{T} (θ) \end{matrix}],

(49)

where

V_{4} (θ) = \cos (θ) [\begin{matrix} \begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix} \end{matrix}] + \sin (θ) [\begin{matrix} \begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ - 1 & 0 & 0 & 0 \\ 0 & - 1 & 0 & 1 \end{matrix} \end{matrix}],

(50)

$V_{4}^{T} (θ)$ is the transpose of $V_{4} (θ)$ ⁠, 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 $S_{\pm ϵ}^{(1), (2)}$ acting on the 16-qubit representation lead to

\begin{gathered} S_{\pm ϵ}^{(1)} [\begin{matrix} \begin{matrix} q_{0} (z, t) \\ q_{1} (z, t) \\ q_{2} (z, t) \\ q_{3} (z, t) \\ q_{4} (z, t) \\ q_{5} (z, t) \\ q_{6} (z, t) \\ q_{7} (z, t) \\ q_{8} (z, t) \\ q_{9} (z, t) \\ q_{10} (z, t) \\ q_{11} (z, t) \\ q_{12} (z, t) \\ q_{13} (z, t) \\ q_{14} (z, t) \\ q_{15} (z, t) \end{matrix} \end{matrix}] = [\begin{matrix} \begin{matrix} q_{0} (z \pm ϵ, t) \\ q_{1} (z \pm ϵ, t) \\ q_{2} (z, t) \\ q_{3} (z, t) \\ q_{4} (z \pm ϵ, t) \\ q_{5} (z \pm ϵ, t) \\ q_{6} (z, t) \\ q_{7} (z, t) \\ q_{8} (z \pm ϵ, t) \\ q_{9} (z \pm ϵ, t) \\ q_{10} (z, t) \\ q_{11} (z, t) \\ q_{12} (z \pm ϵ, t) \\ q_{13} (z \pm ϵ, t) \\ q_{14} (z, t) \\ q_{15} (z, t) \end{matrix} \end{matrix}], \\ S_{\pm ϵ}^{(2)} [\begin{matrix} \begin{matrix} q_{0} (z, t) \\ q_{1} (z, t) \\ q_{2} (z, t) \\ q_{3} (z, t) \\ q_{4} (z, t) \\ q_{5} (z, t) \\ q_{6} (z, t) \\ q_{7} (z, t) \\ q_{8} (z, t) \\ q_{9} (z, t) \\ q_{10} (z, t) \\ q_{11} (z, t) \\ q_{12} (z, t) \\ q_{13} (z, t) \\ q_{14} (z, t) \\ q_{15} (z, t) \end{matrix} \end{matrix}] = [\begin{matrix} \begin{matrix} q_{0} (z, t) \\ q_{1} (z, t) \\ q_{2} (z \pm ϵ, t) \\ q_{3} (z \pm ϵ, t) \\ q_{4} (z, t) \\ q_{5} (z, t) \\ q_{6} (z \pm ϵ, t) \\ q_{7} (z \pm ϵ, t) \\ q_{8} (z, t) \\ q_{9} (z, t) \\ q_{10} (z \pm ϵ, t) \\ q_{11} (z \pm ϵ, t) \\ q_{12} (z, t) \\ q_{13} (z, t) \\ q_{14} (z \pm ϵ, t) \\ q_{15} (z \pm ϵ, t) \end{matrix} \end{matrix}], \end{gathered}

(51)

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⁽¹⁾ and P⁽²⁾, referred to as potential collision operators, are

P^{(1)} (γ) = [\begin{matrix} \begin{matrix} Φ_{4} & 0_{4} & 0_{4} & 0_{4} \\ 0_{4} & Φ_{4} & 0_{4} & 0_{4} \\ 0_{4} & 0_{4} & Φ_{4} & 0_{4} \\ 0_{4} & 0_{4} & 0_{4} & Φ_{4} \end{matrix} \end{matrix}],

(52)

P^{(2)} (γ) = [\begin{matrix} \begin{matrix} {\bar{Φ}}_{8}^{(1)} (γ) & {\bar{Φ}}_{8}^{(2)} (γ) \\ {\bar{Φ}}_{8}^{(2)} (γ) & {\bar{Φ}}_{8}^{(1)} (γ) \end{matrix} \end{matrix}],

(53)

where

Φ_{4} = \cos (γ) I_{4} + \sin (γ) [\begin{matrix} \begin{matrix} 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & - 1 \\ - 1 & 0 & 0 & 0 \\ 0 & - 1 & 0 & 0 \end{matrix} \end{matrix}],

(54)

\begin{gathered} {\bar{Φ}}^{(1)} = \cos (γ) I_{8}, \\ {\bar{Φ}}^{(2)} = \sin (γ) [\begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & - 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & - 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{matrix} \end{matrix}], \end{gathered}

(55)

where γ is a mixing angle.

From the collide-stream operators, we can construct the following unitary operators:

\begin{gathered} U = S_{- ϵ}^{(1)} C (θ) S_{ϵ}^{(1)} C^{†} (θ) S_{ϵ}^{(2)} C (θ) S_{- ϵ}^{(2)} C^{†} (θ), \\ U^{†} = S_{ϵ}^{(1)} C^{†} (θ) S_{- ϵ}^{(1)} C (θ) S_{- ϵ}^{(2)} C^{†} (θ) S_{ϵ}^{(2)} C (θ), \end{gathered}

(56)

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,

Q (z, t + δ t) = P^{(2)} P^{(1)} U^{†} U Q (z, t),

(57)

where $Q^{T} (z, t) = [q_{0} (z, t), q_{1} (z, t), \dots, q_{15} (z, t)]$ ⁠.

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 γ ∝ ϵ². This ordering is akin to diffusion ordering. Subsequently, we expand the evolution equation (57) out to order ϵ² using Mathematica. There is an additional constraint that guides our choice of the mixing angles. At order ϵ², 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

θ = \frac{1}{4 n (z)} ϵ, γ = \frac{1}{2} \frac{n^{'} (z)}{n^{2} (z)} ϵ^{2},

(58)

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 ϵ², (57) leads to the following mesoscopic evolution equation for the 16 qubits:

\frac{\partial}{\partial t} [\begin{matrix} \begin{matrix} q_{0} \\ q_{1} \\ q_{2} \\ q_{3} \\ q_{4} \\ q_{5} \\ q_{6} \\ q_{7} \\ q_{8} \\ q_{9} \\ q_{10} \\ q_{11} \\ q_{12} \\ q_{13} \\ q_{14} \\ q_{15} \end{matrix} \end{matrix}] = - \frac{1}{n (z)} \frac{\partial}{\partial z} [\begin{matrix} \begin{matrix} q_{2} \\ q_{3} \\ q_{0} \\ q_{1} \\ - q_{6} \\ - q_{7} \\ - q_{4} \\ - q_{5} \\ q_{6} \\ q_{11} \\ q_{8} \\ q_{9} \\ - q_{14} \\ - q_{11} \\ - q_{12} \\ - q_{13} \end{matrix} \end{matrix}] - \frac{1}{2} \frac{n^{'} (z)}{n^{2} (z)} [\begin{matrix} \begin{matrix} - q_{2} + q_{15} \\ q_{3} + q_{14} \\ - q_{0} + q_{13} \\ q_{1} + q_{12} \\ - q_{6} - q_{11} \\ q_{7} - q_{10} \\ - q_{4} - q_{9} \\ q_{5} - q_{8} \\ - q_{10} + q_{7} \\ q_{11} + q_{6} \\ - q_{8} + q_{5} \\ q_{9} + q_{4} \\ - q_{14} - q_{3} \\ q_{15} - q_{2} \\ - q_{12} - q_{1} \\ q_{13} - q_{0} \end{matrix} \end{matrix}],

(59)

where all the q_i’s $(i = 0, 1, \dots, 15)$ are functions of z and t. It is straightforward to show that (59) leads to the eight-spinor evolution equations (46) and (47) with the following substitutions:

\begin{gathered} ψ_{0} = q_{0} + q_{2}, ψ_{1} = q_{1} + q_{3}, ψ_{2} = q_{4} + q_{6}, \\ ψ_{3} = q_{5} + q_{7}, ψ_{4} = q_{8} + q_{10}, ψ_{5} = q_{9} + q_{11}, \\ ψ_{6} = q_{12} + q_{14}, ψ_{7} = q_{13} + q_{15} . \end{gathered}

(60)

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 n₂ = 1.5 surrounded by vacuum with n₁ = 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.

FIG. 1.

View large Download slide

Permittivity as a function of z with n₁ = 1 and n₂ = 1.5: (a) a dielectric slab embedded in vacuum; (b) a magnified view of the transition layer.

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,

E_{x} (t = 0) = E_{0} sech \{0.3 (\frac{z - z_{0}}{24})\},

(61)

while the second profile is an exponential cusp of the form

E_{x} (t = 0) = E_{0} \exp \{- 0.3 |\frac{z - z_{0}}{70}|\} .

(62)

In vacuum, the profile for B_y is the same as that for E_x. In the simulations, we will set E₀ = 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 E_x and B_y 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 E_x (blue) and B_y (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 E_x and B_y to be out of phase. The simulation results in Fig. 2(b) are in agreement. Furthermore, from (32), since n₁ < n₂, 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 B_y/E_x = n₂ = 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 $|) E_{r 1} / E_{0} |) = 0.2$ and $|) E_{t 1} / E_{0} |) = 0.9798$ ⁠, respectively. Subscript 1 indicates the first encounter with the vacuum–dielectric interface.

FIG. 2.

View large Download slide

Electric and magnetic field components E_x (blue) and B_y (red), respectively, for (a) the incident pulse at t = 6500 and (b) the reflected and transmitted pulses at t = 12 000. For the incident pulse, E_x = B_y. The boundaries of the dielectric slab are denoted by the vertical dashed lines.

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 E_x = B_y 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 E_t2/E_t1 = 0.9798. For the pulse reflected from the back-end of the slab, (32) gives E_r2/E_t1 = 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.

FIG. 3.

View large Download slide

(a) Pulse profiles at t = 19 000 after scattering from the back-end of the dielectric slab. (b) Different pulse profiles at t = 25 000 after the second scattering event at the front-end of the dielectric slab.

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, E_t3/E_r2 = 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 $\sqrt{n_{2} / n_{1}}$ for propagation from a medium with a refractive index n₁ to a medium with a refractive index n₂. 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.

TABLE I.

Results from QLA simulations of an initial hyperbolic secant pulse. The first column is the simulation time, which correspond to Figs. 2(b), 3(a), and 3(b). The second column is the z location of the peak of a pulse in the simulation domain. The third column gives more information about the pulse: R and T stand for the reflected and transmitted pulses, respectively; 1 and 2 indicate the front-end and the back-end, respectively, of the dielectric where the pulse is generated; V and D stand, respectively, for vacuum and the dielectric slab where the pulse is propagating. B_max and E_max are the maximum values of B_y and E_x, respectively, for a pulse. E_t/E_i and E_r/E_i are, respectively, the transmitted and reflected E_y at the peak of a pulse normalized to the maximum E_y of the incoming pulse at a given vacuum–dielectric interface. From the analytical model, with $\sqrt{n_{2} / n_{1}}$ included in the expression for the transmitted wave amplitude (33), we get E_t/E_i = 0.9798 for all the pulses. In addition, from (32), $|) E_{r} / E_{i} |) = 0.2$ ⁠. The results in columns 6 and 7 are in excellent agreement with these values.

Time $(\times 1 0^{3})$	z	Pulse	B_max $(\times 1 0^{- 3})$	E_max $(\times 1 0^{- 3})$	E_t/E_i	E_r/E_i	B_max/E_max
12	6 615	R1V	1.9505	−1.9505		−0.1951	−1.0
	7 926	T1D	14.6973	9.7982	0.9799		1.5
19	4 523	R1V	1.9505	−1.9505			−1.0
	7 876	R2D	−2.8556	1.9037		0.1943	−1.5
	9 682	T2V	9.6094	9.6094	0.9807		1.0
25	2 730	R1V	1.9504	−1.9504			−1.0
	6 319	T1V	−1.8668	1.8668	0.9806		−1.0
	8 124	R1D	0.5563	0.3709		0.1948	1.4999
	11 475	T2V	9.6088	9.6088			1.0

Time $(\times 1 0^{3})$	z	Pulse	B_max $(\times 1 0^{- 3})$	E_max $(\times 1 0^{- 3})$	E_t/E_i	E_r/E_i	B_max/E_max
12	6 615	R1V	1.9505	−1.9505		−0.1951	−1.0
	7 926	T1D	14.6973	9.7982	0.9799		1.5
19	4 523	R1V	1.9505	−1.9505			−1.0
	7 876	R2D	−2.8556	1.9037		0.1943	−1.5
	9 682	T2V	9.6094	9.6094	0.9807		1.0
25	2 730	R1V	1.9504	−1.9504			−1.0
	6 319	T1V	−1.8668	1.8668	0.9806		−1.0
	8 124	R1D	0.5563	0.3709		0.1948	1.4999
	11 475	T2V	9.6088	9.6088			1.0

FIG. 4.

View large Download slide

Instantaneous Poynting flux S(t) as a function of time. The presence of the dips is discussed in the text.

B. Poynting flux

The instantaneous Poynting flux is defined as

S (t) = \int_{0}^{L} d z E (z, t) \times B (z, t) \cdot \hat{n},

(63)

where $\hat{n}$ 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.

FIG. 5.

View large Download slide

Spatial span of the dielectric slab is 12 000 ≤ z ≤ 15 000. The plots show E_x (blue) and B_y (red): (a) profile of the initial exponential cusp pulse in vacuum where E_x = B_y; (b) reflected and transmitted pulse profiles at t = 20 000 following the passage of the initial pulse through the front-end of the dielectric.

FIG. 6.

View large Download slide

(a) Different pulse profiles at t = 38 000 following scattering from the back-end of the dielectric. (b) The pulse profiles at t = 48 500 after second scattering from the front-end of the dielectric.

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.

TABLE II.

Results from QLA simulations of an initial exponential cusp pulse. The notation used is the same as in Table I.

Time $(\times 1 0^{3})$	z	Pulse	B_max $(\times 1 0^{- 3})$	E_max $(\times 1 0^{- 3})$	E_t/E_i	E_r/E_i	B_max/E_max
20	10 028	R1V	1.8808	−1.8818		−0.1917	0.9995
	13 317	T1D	14.2499	9.5	0.9677		1.5
38	4 647	R1V	1.8773	−1.8777			0.9998
	13 087	R2D	−2.737	1.8241		0.192	−1.5005
	17 858	T2V	9.2389	9.2374	0.9724		1.0001
48.5	1 510	R1V	1.8756	−1.8755			1.0
	10 490	T1V	−1.7857	1.7855	0.9788		1.0001
	13 015	R1D	0.5337	0.3558		0.195	1.5
	20 996	T2V	9.2176	9.2176			1.0

Time $(\times 1 0^{3})$	z	Pulse	B_max $(\times 1 0^{- 3})$	E_max $(\times 1 0^{- 3})$	E_t/E_i	E_r/E_i	B_max/E_max
20	10 028	R1V	1.8808	−1.8818		−0.1917	0.9995
	13 317	T1D	14.2499	9.5	0.9677		1.5
38	4 647	R1V	1.8773	−1.8777			0.9998
	13 087	R2D	−2.737	1.8241		0.192	−1.5005
	17 858	T2V	9.2389	9.2374	0.9724		1.0001
48.5	1 510	R1V	1.8756	−1.8755			1.0
	10 490	T1V	−1.7857	1.7855	0.9788		1.0001
	13 015	R1D	0.5337	0.3558		0.195	1.5
	20 996	T2V	9.2176	9.2176			1.0

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.

FIG. 7.

View large Download slide

Instantaneous Poynting flux S(t) as a function of time. The presence of the dips is discussed in the text.

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 $\sqrt{n_{2} / n_{1}}$ for a pulse traveling from a medium with a refractive index n₁ to a medium with a refractive index n₂. 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.³ 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.⁸ 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 $\sqrt{n_{2} / n_{1}}$ associated with the transmission coefficient.

ACKNOWLEDGMENTS

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.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

DATA AVAILABILITY

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 n₂ = 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 E_x and the magnetic field B_y, 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 E_x and B_y fields, in red, in Figs. 9(a) and 9(b), respectively.

FIG. 8.

$FIG. 8. Transition layer, from vacuum to a dielectric with a refractive index n2 = 2, is 12 lattice units wide.$

View large Download slide

Transition layer, from vacuum to a dielectric with a refractive index n₂ = 2, is 12 lattice units wide.

FIG. 9.

View large Download slide

(a) E_x field component at t = 3500 (blue). The reversal in E_x is due to reflection at the interface. Superimposed in red is E_x of the initial Gaussian pulse in vacuum. (b) The B_y field component at t = 3500 (blue). Superimposed in red is B_y of the initial Gaussian pulse in vacuum.

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.

FIG. 10.

View large Download slide

In red is the instantaneous Poynting flux as a function of time evaluated using the subtraction technique. In blue is the result obtained using the same procedure that led to gaps in Figs. 4 and 7.

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.

REFERENCES

1.

J. R.

Oppenheimer

,

Phys. Rev.

38

,

725

–

746

(

1931

).

https://doi.org/10.1103/physrev.38.725

Crossref

2.

R. H.

Good

, Jr.,

Phys. Rev.

105

,

1914

–

1919

(

1957

).

https://doi.org/10.1103/physrev.105.1914

Crossref

3.

S. A.

Khan

,

Phys. Scr.

71

,

440

–

442

(

2005

); arXiv:0205083v1 (2002).

https://doi.org/10.1238/physica.regular.071a00440

Crossref

4.

D. J.

Griffiths

,

Introduction to Electrodynamics

, 4th ed. (

Pearson

,

2013

), Chap. 7.

Google Scholar

5.

M. K.-H.

Kiesslinga

and

A. S.

Tahvildar-Zadehb

,

J. Math. Phys.

59

,

112302-1

–

112302-33

(

2018

).

https://doi.org/10.1063/1.5021066

Crossref

6.

C. T.

Sebens

,

Found. Phys.

49

,

365

–

389

(

2019

).

https://doi.org/10.1007/s10701-019-00253-3

Crossref

7.

I.

Bialynicki-Birula

and

Z.

Bialynicki-Birula

,

J. Phys. A: Math. Theor.

46

,

053001-1

–

053001-32

(

2013

).

https://doi.org/10.1088/1751-8113/46/5/053001

Crossref

8.

G.

Vahala

,

L.

Vahala

,

M.

Soe

, and

A. K.

Ram

,

J. Plasma Phys.

86

,

905860518-1

–

905860518-17

(

2020

).

https://doi.org/10.1017/s0022377820001166

Crossref

9.

G.

Vahala

,

L.

Vahala

,

M.

Soe

, and

A. K.

Ram

,

Radiat. Eff. Defects Solids

176

,

49

–

63

(

2021

).

https://doi.org/10.1080/10420150.2021.1891058

Crossref

10.

G.

Vahala

,

M.

Soe

, and

L.

Vahala

,

Radiat. Eff. Defects Solids

175

,

113

–

119

(

2020

).

https://doi.org/10.1080/10420150.2020.1718136

Crossref

2021

Author(s)

Reflection and transmission of electromagnetic pulses at a planar dielectric interface: Theory and quantum lattice simulations

I. INTRODUCTION

II. REFLECTION AND TRANSMISSION OF A GAUSSIAN PULSE

A. The incident Gaussian pulse

B. The reflected and transmitted Gaussian pulses

C. Boundary conditions: Amplitudes and widths of the reflected and transmitted pulses

III. MAXWELL EQUATIONS: REPRESENTATION IN TERMS OF RIEMANN–SILBERSTEIN–WEBER VECTORS

A. Propagation and inhomogeneity along the z-direction

IV. QUANTUM LATTICE ALGORITHM FOR MAXWELL EQUATIONS

V. QLA SIMULATIONS FOR TIME EVOLUTION OF AN ELECTROMAGNETIC PULSE

A. Time evolution of a hyperbolic secant pulse

B. Poynting flux

C. Time evolution of an exponential cusp pulse

VI. SUMMARY AND CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

APPENDIX: DETERMINATION OF THE POYNTING FLUX WHEN THE PULSES OVERLAP

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Reflection and transmission of electromagnetic pulses at a planar dielectric interface: Theory and quantum lattice simulations

I. INTRODUCTION

II. REFLECTION AND TRANSMISSION OF A GAUSSIAN PULSE

A. The incident Gaussian pulse

B. The reflected and transmitted Gaussian pulses

C. Boundary conditions: Amplitudes and widths of the reflected and transmitted pulses

III. MAXWELL EQUATIONS: REPRESENTATION IN TERMS OF RIEMANN–SILBERSTEIN–WEBER VECTORS

A. Propagation and inhomogeneity along the z-direction

IV. QUANTUM LATTICE ALGORITHM FOR MAXWELL EQUATIONS

V. QLA SIMULATIONS FOR TIME EVOLUTION OF AN ELECTROMAGNETIC PULSE

A. Time evolution of a hyperbolic secant pulse

B. Poynting flux

C. Time evolution of an exponential cusp pulse

VI. SUMMARY AND CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

APPENDIX: DETERMINATION OF THE POYNTING FLUX WHEN THE PULSES OVERLAP

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only