22 Wave Equations II: Quantum Packets and Electromagnetic
22.1 Quantum Wave Packets
Problem An experiment places an electron with a definite momentum and position in a 1D region of space the size of an atom. It is confined to that region by some kind of attractive potential. Your problem is to determine the resultant electron behavior in time and space.
22.2 Time-Dependent Schrödinger Equation (Theory)
Because the region of confinement is the size of an atom, we must solve this problem quantum mechanically. Because the particle has both a definite momentum and position, it is best described as a wave packet, which implies that we must now solve the time-dependent Schrödinger equation for both the spatial and time dependence of the wave packet. Accordingly, this is a different problem from the bound state one of a particle confined to a box, considered in Chapters 7 and 9, where we solved the eigenvalue problem for stationary states of the time-independent Schrödinger equation.
We model an electron initially localized in space at x = 5 with momentum k0 ( = 1 in our units) by a wave function that is a wave packet consisting of a Gaussian multiplying a plane wave:
Figure 22.1 The position as a function of time of a localized electron confined to a square well (computed with the code SqWell.py available in the instructor’s site). The electron is initially on the left with a Gaussian wave packet. In time, the wave packet spreads out and collides with the walls.
Figure 22.2 The probability density as a function of time for an electron confined to a 1D harmonic oscillator potential well. (a) A conventional surface plot, (b) a color visualization.
To solve the problem, we must determine the wave function for all later times. If (22.1) was an eigenstate of the Hamiltonian, its exp(−iωt) time dependence can be factored out of the Schrödinger equation (as is usually carried out in textbooks). However, for this ψ, and so we must solve the full time-dependent Schrödinger equation. To show you where we are going, the resulting wave packet behavior is shown in Figures 22.1 and 22.2.
The time and space evolution of a quantum particle is described by the 1D time-dependent Schrödinger equation,
(22.2)
(22.3)
where we have set 2m = 1 to keep the equations simple. Because the initial wave function is complex (in order to have a definite momentum associated with it), the wave function will be complex for all times. Accordingly, we decompose the wave function into its real and imaginary parts:
(22.4)
(22.5)
(22.6)
(22.7)
where V(x) is the potential acting on the particle.
22.2.1 Finite-Difference Algorithm
The time-dependent Schrödinger equation can be solved with both implicit (large-matrix) and explicit (leapfrog) methods. The extra challenge with the Schrödinger equation is to establish that the integral of the probability density remains constant (conserved) to a high level of precision for all time. For our project, we use an explicit method that improves the numerical conservation of probability by solving for the real and imaginary parts of the wave function at slightly different or “staggered” times (Askar and Cakmak, 1977; Visscher, 1991; Maestri et al., 2000). Explicitly, the real part R is determined at times 0, Δt, …, and the imaginary part I at 1/2Δt, 3/2Δt, …. The algorithm is based on (what else?) the Taylor expansions of R and I:
(22.8)
where α = Δt/2(Δx)2. In the discrete form with , we have
Although probability is not conserved exactly with this algorithm, the error is two orders higher than that in the wave function, and this is usually quite satisfactory. If it is not, then we need to use smaller steps. While this definition of ρ may seem strange, it reduces to the usual one for Δt → 0 and so can be viewed as part of the art of numerical analysis. We will ask you to investigate just how well probability is conserved. We refer the reader to Koonin (1986) and Visscher (1991) for details on the stability of the algorithm.
22.2.2 Wave Packet Implementation, Animation
In Listing 22.1, you will find the program HarmosAnimate.py that solves for the motion of the wave packet (22.1) inside a harmonic oscillator potential. The program Slit.py on the instructor’s site solves for the motion of a Gaussian wave packet as it passes through a slit (Figure 22.3). You should solve for a wave packet confined to the square well:
Define arrays psr[751,2] and psi[751,2] for the real and imaginary parts of ψ, and Rho[751] for the probability. The first subscript refers to the x position on the grid, and the second to the present and future times.
Use the values σ0 = 0.5, Δx = 0.02, k0 = 17π, and Δt = 1/2Δx2.
Figure 22.3 The probability density as a function of position and time for an electron incident upon and passing through a slit. Significant reflection is seen to occur.
Use (22.1) for the initial wave packet to define psr[j,1] for all j at t = 0 and to define psi[j,1] at t = 1/2Δt.
Set Rho[1] = Rho[751] = 0.0 because the wave function must vanish at the infinitely high well walls.
Increment time by 1/2Δt. Use (22.9) to compute psr[j,2] in terms of psr[j,1], and (22.10) to compute psi[j,2] in terms of psi[j,1].
Repeat the steps through all of space, that is, for i = 2−750.
Throughout all of space, replace the present wave packet (second index equal to 1) by the future wave packet (second index 2).
After you are sure that the program is running properly, repeat the time-stepping for ∼ 5000 steps.
Listing 22.1HarmosAnimate.py solves the time-dependent Schrödinger equation for a particle described by a Gaussian wave packet moving within a harmonic oscillator potential.
Animation: Output the probability density after every 200 steps for use in animation.
Make a surface plot of probability vs. position vs. time. This should look like Figure 22.1 or 22.2.
Make an animation showing the wave function as a function of time.
Check how well the probability is conserved for early and late times by determining the integral of the probability over all of space, , and seeing by how much it changes in time (its specific value doesn’t matter because that is just normalization).
What might be a good explanation of why collisions with the walls cause the wave packet to broaden and break up? (Hint: The collisions do not appear so disruptive when a Gaussian wave packet is confined within a harmonic oscillator potential well.)
Figure 22.4 The probability density as a function of x and y at three different times for an electron confined to a 2D parabolic tube (infinite in the y direction). The electron is initially placed in a Gaussian wave packet in both the x and y directions, and it is to be noted how there is spreading of the wave packet in the y direction, but not in the x direction.
22.2.3 Wave Packets in Other Wells (Exploration)
1D Well Now confine the electron to a harmonic oscillator potential:
(22.12)
Take the momentum k0 = 3π, the space step Δx = 0.02, and the time step . Note that the wave packet broadens yet returns to its initial shape!
2D Well Now confine the electron to a 2D parabolic tube (Figure 22.4):
(22.13)
The extra degree of freedom means that we must solve the 2D PDE:
(22.14)
Assume that the electron’s initial localization is described by the 2D Gaussian wave packet:
(22.15)
Note that you can solve the 2D equation by extending the method we just used in 1D or you can look at the next section where we develop a special algorithm.
22.3 Algorithm for the 2D Schrödinger Equation
One way to develop an algorithm for solving the time-dependent Schrödinger equation in 2D is to extend the 1D algorithm to another dimension. Rather than that, we apply quantum theory directly to obtain a more powerful algorithm (Maestri et al., 2000). First, we note that Equation 22.14 can be integrated in a formal sense (Landau and Lifshitz, 1976) to obtain the operator solution:
where is an operator that translates a wave function by an amount of time t and is the Hamiltonian operator. From this formal solution, we deduce that a wave packet can be translated ahead by time Δt via
(22.18)
where the superscripts denote time t = nΔt and the subscripts denote the two spatial variables x = iΔx and y = jΔy. Likewise, the inverse of the time evolution operator moves the solution back one time step:
While it would be nice to have an algorithm based on a direct application of (22.19), the references show that the resulting algorithm would not be stable. That being so, we base our algorithm on an indirect application (Askar and Cakmak, 1977), namely, the relation between the difference in ψn+1 and ψn-1:
where the difference in sign of the exponents is to be noted. The algorithm derives from combining the 0(Δx2) expression for the second derivative obtained from the Taylor expansion,
(22.21)
with the corresponding-order expansion of the evolution equation (22.20). Substituting the resulting expression for the second derivative into the 2D time-dependent Schrödinger equation results in1)
(22.22)
where α = Δt/2(Δx)2. We convert this complex equations to coupled real equations by substituting them into the wave function ψ = R + iI,
(22.23)
(22.24)
This is the algorithm we use to integrate the 2D Schrödinger equation. To determine the probability, we use the same expression (22.11) used in 1D.
22.3.1 Exploration: Bound and Diffracted 2D Packet
Determine the motion of a 2D Gaussian wave packet within a 2D harmonic oscillator potential:
(22.25)
Center the initial wave packet at (x, y) = (3.0, −3) and give it momentum (k0x, k0y) = (3.0, 1.5).
Young’s single-slit experiment has a wave passing through a small slit with the transmitted wave showing interference effects. In quantum mechanics, where we represent a particle by a wave packet, this means that an interference pattern should be formed when a particle passes through a small slit. Pass a Gaussian wave packet of width 3 through a slit of width 5 (Figure 22.3), and look for the resultant quantum interference.
We have just seen how to represent a quantum particle as a wave packet and how to compute the interaction of that wave packet/particle with an external potential. Although external potentials do exist in nature, realistic scattering often involves the interaction of one particle with another, which in turn would be represented as the interaction of a wave packet with a different wave packet. We have already done the hard work needed to compute wave packet–wave packet scattering in our implementation of wave packet–potential scattering even in 2D. We now need only to generalize it a bit.
Two interacting particles are described by the time-dependent Schrödinger equation in the coordinates of the two particles
where, for simplicity, we assume a one-dimensional space and again set = 1. Here mi and xi are the mass and position of particle i = 1, 2. Knowledge of the two-particle wavefunction ψ(x1, x2, t) at time t permits the calculation of the probability density of particle 1 being at x1 and particle 2 being at x2:
(22.27)
The fact that particles 1 and 2 must be located somewhere in space leads to the normalization constraint on the wavefunction
The description of a particle within a multiparticle system by a single-particle wavefunction is an approximation unless the system is uncorrelated, in which case the total wavefunction can be written in product form. However, it is possible to deduce meaningful one-particle probabilities (also denoted by ρ) from the two-particle density by integrating over the other particle:
Of course, the true solution of the Schrödinger equation is ψ(x1,x2, t), but we find it hard to unravel the physics in a three-variable complex function, and so will usually view ρ(x1, t) and ρ(x2, t) as two separate wave packets colliding.
If particles 1 and 2 are identical, then their total wavefunction must be symmetric (s) or antisymmetric (a) under interchange of the particles. We impose this condition on our numerical solution ψ(x1, x2), by forming the combinations
The cross term in (22.31) places an additional correlation into the wave packets.
22.4.1 Algorithm
The algorithm for the solution of the two-particle Schrödinger equation (22.26) in one dimension x, is similar to the one we outlined for the solution for a single particle in the two dimensions x and y, only now with x → x1 and y → x2. As usual, we assume discrete space and time
(22.32)
We employ the improved algorithm for the time derivative introduced by Askar and Cakmak (1977), which uses a central difference algorithm applied to the formal solution (22.16)
where we have assumed Δx1 = Δx2 = Δx and formed the ratio λ = Δt/Δx2. Again we will take advantage of the extra degree of freedom provided by the complexity of the wavefunction
(22.35)
and staggered times to preserve probability better (Visscher, 1991). The algorithm (22.34) then separates into the pair of coupled equations
(22.36)
(22.37)
22.4.2 Implementation
We assume that the particle–particle potential is central and depends only on the relative distance between particles 1 and 2 (the method can handle any x1 and x2 functional dependences if needed). For example, we have used both a “soft” potential with a Gaussian dependence and a “hard” potential with a square-well dependence:
(22.38)
where α is the range parameter and V0 is the depth parameter.
Because we are solving a PDE, we must specify both initial and boundary conditions. For scattering, we assume that particle 1 is initially at with momentum k1, and that particle 2 is initially far away at with momentum k2. Since the particles are initially too far apart to be interacting, we assume that the initial wave packet is a product of independent wave packets for particles 1 and 2
Because of the Gaussian factors here, ψ is not an eigenstate of the momentum operator for either particle 1 or 2, but instead contains a spread of momenta about the mean, initial momenta k1 and k2. If the wave packet is made very broad (σ → ∞), we would obtain momentum eigenstates, but then would have effectively eliminated the wave packets. Note that while the Schrödinger equation may separate into one equation in the relative coordinate x = x1 − x2 and another in the center-of-mass coordinate X = (m1x1 + m2x2)/(m1 + m2), the initial condition (22.39), or more general ones, cannot be written as a product of separate functions of x and X, and accordingly, a solution of the partial differential equation in two variables is required (Landau, 1996).
We start the staggered-time algorithm with the real part the wavefunction (22.39) at t = 0 and the imaginary part at t = Δt/2. The initial imaginary part follows by assuming that Δt/2 is small enough and σ is large enough for the initial time dependence of the wave packet to be that of the plane wave parts:
(22.40)
In an actual scattering experiment, a projectile starts at infinity and the scattered particles are observed also at infinity. We model that by solving our partial differential equation within a box of side L that in ideal world would be much larger than both the range of the potential and the width of the initial wave packet. This leads to the boundary conditions
(22.41)
The largeness of the box minimizes the effects of the boundary conditions during the collision of the wave packets, although at large times there will be artificial collisions with the box that do not correspond to actual experimental conditions.
Some typical parameters we used are , , , and . The original C code is available on the Web (Maestri et al., 2000). Note that our space step size is 1/1400th of the size of the box L, and l/70th of the size of the wave packet. Our time step is 1/20 000th of the total time T, and l/2000th of a typical time for the wave packet. In all cases, the potential and wave packet parameters are chosen to be similar to those used in the one-particle studies by Goldberg et al. (1967). The time and space step sizes were determined by trial and error until values were found that provided stability and precision. Importantly, ripples during interactions found in earlier studies essentially disappear when (more accurate) small values of Δx are employed. In general, stability is obtained by making Δt small enough, with simultaneous changes in Δt and Δx made to keep λ = Δt/Δx2 constant. Total probability, as determined by a double Simpson’s rule integration of (22.28), is typically conserved to 13 decimal places, impressively close to machine precision. In contrast, the mean energy, for which we do not use an optimized algorithm, is conserved only to 3 places.
22.4.3 Results and Visualization
We solve our problem in the center-of-momentum system by taking k2 = −k1 (particle 1 moving to larger x values and particle 2 to smaller x). Because the results are time dependent, we make movies of them, and since the movies show the physics much better than the still images, we recommend that the reader look at them (Maestri et al., 2000). We first tested the procedure by emulating the one-particle collisions with barriers and wells studied by Goldberg et al. (1967) and presented by Schiff. We made particle 2 ten times heavier than particle 1, which means that particle 2’s initial wave packet moves at 1/10th the speed of particle 1’s, and so is similar to a barrier. On the left of Fig. 22.5, we see six frames from an animation of the two-particle density ρ(x1, x2, t) as a simultaneous function of the particle positions x1 and x2. On the right of Fig. 22.5 we show, for this same collision, the single-particle densities ρ1(x1, t) and ρ2(x2, t) extracted from ρ(x1, x2, t) by integrating out the dependence on the other particle via (22.29). Because the mean kinetic energy equals twice the maximum height of the potential barrier, we expect complete penetration of the packets, and indeed, at time 18 we see on the right that the wave packets have large overlap, with the repulsive interaction “squeezing” particle 2 (it gets narrower and taller). During time 22–40, we see part of wave packet 1 reflecting off wave packet 2 and then moving back to smaller x (the left). From times 26–55, we also see that a major part of wave packet 1 gets “trapped” inside of wave packet 2 and then leaks out rather slowly.
When looking at the two-particle density ρ(x1,x2, t) on the left of Fig. 22.5, we see that for times 1-26, the x2 position of the peak of changes very little with time, which is to be expected since particle 2 is heavy. In contrast, the x1 dependence in ρ(x1, x2, t) gets broader with time, develops into two peaks at time 26, separates into two distinct parts by time 36, and then, at time 86 after reflecting off the walls, returns to particle 2’s position. We also notice in both these figures that at time 40 and thereafter, particle 2 (our “barrier”) fissions into reflected and transmitted waves. As this comparison of the visualizations on the right and left of Figures 22.5 demonstrates, it seems easier to understand the physics by superimposing two single-particle densities (thereby discarding information on correlations) than by examining the two-particle density.
Figure 22.5 (a) Six frames from an animation of the two-particle density ρ(x1, x2, t)as a function of the position of particle 1 with mass m and of the position of particle 2 with mass 10m. (b) This same collision as seen with the single-particle densities ρ(x1, t) and ρ(x2, t). The numbers in the left-hand corners are the times in units of 100Δt. Note that each plot ends at the walls of the containing box, and that particle 1 “bounces off” a wall between times 36 and 86.
Figure 22.6 A time sequence of two single-particle wave packets scattering from each other. The particles have equal mass, a mean kinetic energy equal to a quarter of the well’s depth, and the wavefunction has been symmetrized.
In Fig. 22.6, we see nine frames from the movie of an attractive m–m collision in which the mean energy equals one-quarter of the well depth. The initial packets speed up as they approach each other, and at time 60, the centers have already passed through each other. After that, a transmitted and reflected wave for each packet is seen to develop (times 66–78). Although this may be just an artifact of having two particles of equal mass, from times 110–180, we see that each packet appears to capture or “pick up” a part of the other packet and move off with it.
Note in Fig. 22.6 that at time 180, the wave packets are seen to be interacting with the wall (the edges of the frames), as indicated by the interference ripples between incident and reflected waves. Also note that at time 46 and thereafter two additional small wave packets are seen to be traveling in opposite directions to the larger initial wave packets. These small packets are numerical artifacts and arise because outgoing waves satisfy the same differential equations as incoming waves.
22.5 E&M Waves via Finite-Difference Time Domain
Problem You are given a region in space in which the E and H fields are known to have a sinusoidal spatial variation
(22.42)
(22.43)
with all other components vanishing. Determine the fields for all z values at all subsequent times.
22.6 Maxwell’s Equations
The description of EM waves via Maxwell’s equations is given in many textbooks. For propagation in just one dimension (z) and for free space with no sinks or sources, four coupled PDEs result:
Figure 22.7 A single EM pulse traveling along the z-axis. The coupled E and H pulses are indicated by solid and dashed curves, respectively, and the pulses at different z values correspond to different times.
As indicated in Figure 22.7, we have chosen the electric field E(z, t) to oscillate (be polarized) in the x direction and the magnetic field H(z, t) to be polarized in the y direction. As indicated by the bold arrow in Figure 22.7, the direction of power flow for the assumed transverse electromagnetic (TEM) wave is given by the right-hand rule applied to E × H. Note that although we have set the initial conditions such that the EM wave is traveling in only one dimension (z), its electric field oscillates in a perpendicular direction (x), and its magnetic field oscillates in yet a third direction (y); so while some may call this a 1D wave, the vector nature of the fields means that the wave occupies all three dimensions.
22.7 FDTD Algorithm
We need to solve the two coupled PDEs (22.46) and (22.47) appropriate for our problem. As is usual for PDEs, we approximate the derivatives via the central-difference approximation, here in both time and space. For example,
(22.48)
(22.49)
We next substitute the approximations into Maxwell’s equations and rearrange the equations into the form of an algorithm that advances the solution through time. Because only first derivatives occur in Maxwell’s equations, the equations are simple, although the electric and magnetic fields are intermixed.
As introduced by Yee (Yee, 1966), we set up a space–time lattice (Figure 22.8) in which there are half-integer time steps as well as half-integer space steps. The magnetic field will be determined at integer time sites and half-integer space sites (open circles), and the electric field will be determined at half-integer time sites and integer space sites (filled circles). While this is an extra level of complication, the transposed lattices do lead to an accurate and robust algorithm. Because the fields already have subscripts indicating their vector nature, we indicate the lattice position as superscripts, for example,
Figure 22.8 The algorithm for using the known values of Ex and Hy at three earlier times and three different space positions to obtain the solution at the present time. Note that the values of Ex are determined on the lattice of filled circles, corresponding to integer space indices and half-integer time indices. In contrast, the values of Hy are determined on the lattice of open circles, corresponding to half-integer space indices and integer time indices.
To repeat, this formulation solves for the electric field at integer space steps (k) but half-integer time steps (n), while the magnetic field is solved for at half-integer space steps but integer time steps.
We convert these equations into two simultaneous algorithms by solving for Ex at time n + 1/2, and Hy at time n:
The algorithms must be applied simultaneously because the space variation of Hy determines the time derivative of Ex, while the space variation of Ex determines the time derivative of Hy (Figure 22.8). These algorithms are more involved than our usual time-stepping ones in that the electric fields (filled circles in Figure 22.8) at future times t = n + 1/2 are determined from the electric fields at one time step past t = n − 1/2, and the magnetic fields at half a time step past t = n. Likewise, the magnetic fields (open circles in Figure 22.8) at future times t = n + 1 are determined from the magnetic fields at one time step past t = n, and the electric field at half a time step past t = n + 1/2. In other words, it is as if we have two interleaved lattices, with the electric fields determined for half-integer times on lattice 1 and the magnetic fields at integer times on lattice 2.
Although these half-integer times appear to be the norm for FDTD methods (Taflove and Hagness , 1989; Sullivan, 2000), it may be easier for some readers to understand the algorithm by doubling the index values and referring to even and odd times:
(22.53)
(22.54)
This makes it clear that E is determined for even space indices and odd times, while H is determined for odd space indices and even times.
We simplify the algorithm and make its stability analysis simpler by renormalizing the electric fields to have the same dimensions as the magnetic fields,
Here, c is the speed of light in vacuum and β is the ratio of the speed of light to grid velocity Δz/Δt.
The space step Δz and the time step Δt must be chosen so that the algorithms are stable. The scales of the space and time dimensions are set by the wavelength and frequency, respectively, of the propagating wave. As a minimum, we want at least 10 grid points to fall within a wavelength:
(22.59)
The time step is then determined by the Courant stability condition (Taflove and Hagness , 1989; Sullivan, 2000) to be
As we have seen before, (22.60) implies that making the time step smaller improves precision and maintains stability, but making the space step smaller must be accompanied by a simultaneous decrease in the time step in order to maintain stability (you should check this).
Listing 22.2FDTD.py solves Maxwell’s equations via FDTD time stepping (finite-difference time domain) for linearly polarized wave propagation in the z direction in free space.
Figure 22.9 The E field (light) and the H field (dark) at the initial time (a) and at a later time (b). Periodic boundary conditions are used at the ends of the spatial region, which means that the large z wave continues into the z = 0 wave.
22.7.1 Implementation
In Listing 22.2, we provide a simple implementation of the FDTD algorithm for a z lattice of 200 sites and in Figure 22.9 we show some results. The initial conditions correspond to a sinusoidal variation of the E and H fields for all z values in for 0 ≤ z ≤ 200:
(22.61)
The algorithm then steps out in time for as long as the user desires. The discrete form of Maxwell equations used are:
where 1 ≤ k ≤ 200, and beta is a constant. The second index takes the values 0 and 1, with 0 being the old time and 1 the new. At the end of each iteration, the new field throughout all of space becomes the old one, and a new one is computed. With this algorithm, the spatial endpoints k=0 and k=xmax-1 remain undefined. We define them by assuming periodic boundary conditions:
22.7.2 Assessment
Impose boundary conditions such that all fields vanish on the boundaries. Compare the solutions so obtained to those without explicit conditions for times less than and greater than those at which the pulses hit the walls.
Examine the stability of the solution for different values of Δz and Δt and thereby test the Courant condition (22.60).
Extend the algorithm to include the effect of entering, propagating through, and exiting a dielectric material placed within the z integration region.
Ensure that you see both transmission and reflection at the boundaries.
Investigate the effect of varying the dielectric’s index of refraction.
The direction of propagation of the pulse is in the direction of E × H, which depends on the relative phase between the E and H fields. (With no initial H field, we obtain pulses both to the right and the left.)
Modify the program so that there is an initial H pulse as well as an initial E pulse, both with a Gaussian times a sinusoidal shape.
Verify that the direction of propagation changes if the E and H fields have relative phases of 0 or π.
Investigate the resonator modes of a wave guide by picking the initial conditions corresponding to plane waves with nodes at the boundaries.
Investigate standing waves with wavelengths longer than the size of the integration region.
Simulate unbounded propagation by building in periodic boundary conditions into the algorithm.
Place a medium with periodic permittivity in the integration volume. This should act as a frequency-dependent filter, which does not propagate certain frequencies at all.
22.7.3 Extension: Circularly Polarized Waves
We now extend our treatment to EM waves in which the E and H fields, while still transverse and propagating in the z direction, are not restricted to linear polarizations along just one axis. Accordingly, we add to (22.46) and (22.47):
(22.62)
(22.63)
When discretized in the same way as (22.51) and (22.52), we obtain
(22.64)
(22.65)
To produce a circularly polarized traveling wave, we set the initial conditions
(22.66)
(22.67)
We take the phases to be so that their difference which leads to circular polarization. We include the initial conditions in the same manner as we did the Gaussian pulse, only now with these cosine functions.
Listing 22.3 gives our implementation EMcirc.py for waves with transverse two-component E and H fields. Some results of the simulation are shown in Figure 22.10, where you will note the difference in phase between E and H.
Figure 22.10E and H fields at t = 100 for a circularly polarized wave in free space.
Listing 22.3CircPolartzn.py solves Maxwell’s equations via FDTD time-stepping for circularly polarized wave propagation in the z direction in free space.
22.8 Application:Wave Plates
Problem Develop a numerical model for a wave plate that convert a linearly polarized EM wave into a circularly polarized one.
As can be seen in Figure 22.11a wave plate is an optical device that alters the polarization of light traveling through it by shifting the relative phase of the components of the polarization vector. A quarter-wave plate introduces a relative phase of λ/4, where λ is the wavelength of the light. Physically, a wave plate is often a birefringent crystal in which the different propagation velocities of waves in two orthogonal directions leads to the phase change. The amount of phase change is adjusted by varying the thickness of the plate.
Figure 22.11 One frame from the program quarterwave.py (on Instructor’s site) showing a linearly polarized EM wave entering a quarter-wave plate from the left and leaving as a circularly polarized wave on the right (the arrow on the left oscillates back and forth at 45° while the one on the right rotates).
To attack this problem, we apply our FDTD method of solving Maxwell equations. We start with a linear polarized wave with both Ex and Ey components propagating along the z direction. The wave enters the plate and emerges from it still traveling in the z direction, but now with the relative phase of these fields shifted. Of course, because this is an EM wave, there will also be coupled magnetic field components present, in this case Hx and Hy, and they too will need be computed.
Theory Maxwell equations for a wave propagating along the z-axis are:
(22.68)
(22.69)
We take as initial conditions a wave incident from the left along the z-axis, linearly polarized (electric field direction of 45°), with corresponding, and perpendicular, H components:
(22.70)
(22.71)
Because only the relative phases matter, we simplify the calculation by assuming that the Ey and Hx components do not have their phases changed, but that the Ex and Hy components do (in this case by λ/4 when they leave the plate). Of course, after leaving the plate and traveling in free space, there are no further changes in the relative phase.
22.9 Algorithm
As in Section 22.7 and Figure 22.8, we follow the FDTD approach of using known values of Ex and Hy at three earlier times and three different space positions to obtain the solution at the present time. With the renormalized electric fields as in (22.55), this leads to the beautifully symmetric equations:
After each time step, impose a gradual increment of the phase so that the total phase change will be one-quarter of a wavelength. Our program for this, quarterplat.py, is on the instructor’s page.
Verify that the plate converts an initially linearly polarized wave into a circularly polarized one.
Verify that the plate converts an initially circularly polarized wave into a linearly polarized one.
What happens if you put two plates together? Three? Four? (Verify!)
= 1 in our units) by a wave function that is a wave packet consisting of a Gaussian multiplying a plane wave:
for this ψ, and so we must solve the full time-dependent Schrödinger equation. To show you where we are going, the resulting wave packet behavior is shown in Figures 22.1 and 22.2.
remains constant (conserved) to a high level of precision for all time. For our project, we use an explicit method that improves the numerical conservation of probability by solving for the real and imaginary parts of the wave function at slightly different or “staggered” times (Askar and Cakmak, 1977; Visscher, 1991; Maestri et al., 2000). Explicitly, the real part R is determined at times 0, Δt, …, and the imaginary part I at 1/2Δt, 3/2Δt, …. The algorithm is based on (what else?) the Taylor expansions of R and I:
, we have
, and seeing by how much it changes in time (its specific value doesn’t matter because that is just normalization).
. Note that the wave packet broadens yet returns to its initial shape!
is an operator that translates a wave function by an amount of time t and 
is the Hamiltonian operator. From this formal solution, we deduce that a wave packet can be translated ahead by time Δt via
with momentum k1, and that particle 2 is initially far away at 
with momentum k2. Since the particles are initially too far apart to be interacting, we assume that the initial wave packet is a product of independent wave packets for particles 1 and 2
for either particle 1 or 2, but instead contains a spread of momenta about the mean, initial momenta k1 and k2. If the wave packet is made very broad (σ → ∞), we would obtain momentum eigenstates, but then would have effectively eliminated the wave packets. Note that while the Schrödinger equation may separate into one equation in the relative coordinate x = x1 − x2 and another in the center-of-mass coordinate X = (m1x1 + m2x2)/(m1 + m2), the initial condition 
, 
, 
, and 
. The original C code is available on the Web (Maestri et al., 2000). Note that our space step size is 1/1400th of the size of the box L, and l/70th of the size of the wave packet. Our time step is 1/20 000th of the total time T, and l/2000th of a typical time for the wave packet. In all cases, the potential and wave packet parameters are chosen to be similar to those used in the one-particle studies by Goldberg et al. (1967). The time and space step sizes were determined by trial and error until values were found that provided stability and precision. Importantly, ripples during interactions found in earlier studies essentially disappear when (more accurate) small values of Δx are employed. In general, stability is obtained by making Δt small enough, with simultaneous changes in Δt and Δx made to keep λ = Δt/Δx2 constant. Total probability, as determined by a double Simpson’s rule integration of 
so that their difference 
which leads to circular polarization. We include the initial conditions in the same manner as we did the Gaussian pulse, only now with these cosine functions.
