"Superlumial" phenomena in linear optics

S. M. Rochester and D. Budker
UC Berkeley Physics Dept.
Sept. 9, 2002

Introduction

Recently, both the scientific and popular (e.g., NYT and NPR) media have paid much attention to optical phenomena dubbed "slow" and "fast" light. When a pulse of weak probe light propagates through a nonlinear medium in the presence of a strong drive field, the light pulse can, under certain circumstances, be transmitted without much distortion, and the apparent group velocity can be very slow (∼1 m  (s^(-1)) in some experiments), faster than c, or negative (see recent comprehensive reviews by Matsko et. al., 2001; Boyd and Gauthier, 2002). Although the recent  studies of these phenomena have been in the context of nonlinear optics, in which at least two light fields interact with an atomic medium, these effects can be understood by considering the simpler case of linear optics.

In this tutorial, we use considerations from linear optics to show how the variation in group velocity comes about. Although the physics is quite straightforward, and has been worked out in detail both theoretically (Garrett and McCumber, 1970) and experimentally (Chu and Wong, 1982), we attempt to remove any remaining mystery from the subject by discussing it in very basic language and illustrating each step with graphics derived directly from the discussion. (See also Sprangle et. al., 2001; McDonald, 2001; Dogariu et. al. 2001).

We consider an optical medium consisting of an atomic vapor, with light whose frequency is tuned near an isolated resonance. We introduce the complex refractive index, and first describe the propagation of a monochromatic light wave through a slab of the medium. The imaginary part of the refractive index leads to light absorption, while the real part affects the light's phase. Next, we consider a smooth (Gaussian) pulse incident on the slab of our atomic medium. We decompose the pulse into Fourier components, and analyze the propagation of each of these components through the medium separately. This is justified in the case of linear medium whose optical properties are independent of the presence of light. Finally, we recombine the resulting Fourier components into a "new" pulse, and analyze its spatial and temporal distribution. Clearly, this discussion does not contain any physical or mathematical concepts beyond those found in any introductory physics course. Nevertheless, this will be sufficient to illustrate all of our conclusions.

Mathematica Setup

Turn off annoying warnings

Load packages we'll use

Set options for plotting

Set ComplexExpand to behave the way we like

Discussion of light propagation in linear media

Definition of the light field

Any plane light wave propagating in the z direction can be decomposed into a number of sinusoidal waves. We use the complex notation in which each field is described by an amplitude times the real part of

EField[
        k_, ω_] = (E^(I (k z  -  ω t)))

where

k[n_, ω_] = ((n ω)/c)

is the wavenumber, which depends on the light frequency ω and the index of refraction of the medium n(ω). c is the speed of light in vacuum. The real field is given by

ReEField[k_, ω_] = 
      ComplexExpand[Re[EField[k, ω]], k]

Table[
        Plot[ReEField[2  π, 2  π], {z, 0, 8}, 
          PlotRange -> {{0, 8}, {-1.25, 1.25}}], {t, 
          0,  .9,  .1}];

A monochromatic wave 
(double-click plot to view as a function of time)

Phase velocity — general medium

The phase velocity is the velocity of a point of constant phase of the electric field, for example, a zero crossing. To find the phase velocity we solve for z such that Re[E(z)] = 0 as a function of time:

zConstPhase = Solve[ReEField[k, ω] == 0, 
          z][[1]]

The phase velocity is given by ((d z)/(d t))

Vph[k_, ω_] = D[z /. zConstPhase, t]

In terms of the index of refraction

VphOfN[n_] = ComplexExpand[Vph[k[n, ω], ω], n]

In a vacuum, n = 1 and the phase velocity is

VphOfN[1]

Wave propagation through linear atomic media

In a linear medium near an atomic resonance of frequency ω0 and width Γ, the index of refraction is given by (see, for example, Lipson, 1995, Ch. 13.3.2; Griffiths, 1999, Ch. 8.4)

nMedium[ω_] = 
      1 - ((α c)/(ω - ω0 + I Γ/2)) (Γ/(4 ω0))

where α is the absorption coefficient on resonance. The real part of the index of refraction, responsible for phase shifts, is given by

RenMedium[ω_] = 
      Collect[ComplexExpand[
            Re[nMedium[ω]]] /. ((c α Γ)/( (((Γ^2)/4) + (((ω - ω0))^2)))) -> 
            4 ξ , ξ]

where we have defined the parameter ξ to be

ξ[ω_] = ((c α Γ)/( 4 (((Γ^2)/4) + (((ω - ω0))^2))))

The following function plots the real and imaginary parts of the index of refraction vs. frequency:

The phase velocity in the medium is given by

VphMedium[ω_] = 
      Vph[k[RenMedium[ω], ω], ω] // 
        ComplexExpand

Since ξ is always positive, we see that the phase velocity is c on resonance, less than  c below resonance, and greater than c above resonance. To illustrate this, we can plot the light electric field as it traverses a medium. The wavenumber in the medium is given by

kMedium[ω_] = 
      k[nMedium[ω], ω] // Simplify

and the electric field is given by

EFieldMedium[z_, t_, ω_] = 
      EField[kMedium[ω], ω] // Simplify

The light intensity is proportional to Abs[(E^2)], where Abs stands for absolute value. We can verify that the intensity has the correct dependence on the absorption coefficient:

Abs[EFieldMedium[z, t, ω0]^2] // ComplexExpand // 
      PowerExpand

The following function produces an animation of a monochromatic wave traversing a medium.

Here is a plot of the electric field for a light wave tuned above resonance. The phase of the wave at the output is advanced relative to what it would have been, had there been no medium (dotted line). The amplitude of the wave is reduced because of absorption (imaginary part of the index of refraction). Try changing the parameters and re-evaluating. Tuning the light frequency closer to resonance reduces the phase shift and increases absorption. Tuning below resonance (i.e. making ω < ω0), one sees that the phase of the wave at the output is retarded relative to what it would have been, had there been no medium (dotted line) because the phase velocity is <c in this case.

PlotMonochromatic[{0, 3, 5, 8}, {0. , 
        0.95,  .05}, {c -> 1, α -> 3. , ω -> 
          2.1  π, Γ -> 1, ω0 -> 
          2  π}]

A monochromatic wave propagating through a linear medium.
(double-click plot to view as a function of time)

Constructing a light pulse

Introduction

To construct a Gaussian pulse of width Γp, we superpose waves with amplitudes given by

Envelope[Δ_] = (1/(Γp (Sqrt[ 2  π]))) (E^(-((Δ^2)/(2  Γ(p^2)))))

where Δ is the detuning from the central frequency ω, and the first factor provides overall normalization. The more waves of different frequencies we add, the more the resultant envelope looks like a pulse, as we will now show.

Discrete frequencies

Two Frequencies

For illustration, let us look at waves traveling in vacuum, in which n = 1. A "pulse" with two frequencies is given by

TwoFreq = Sum[
            Envelope[Δ] 
              EField[k[
                  1, ω + Δ], ω + Δ], {Δ, -3 π/2, 3  π/2, 3  π}] // 
          ComplexExpand // FullSimplify

The envelope is given by the absolute value of the field

TwoFreqEnvelope = Abs[TwoFreq] // ComplexExpand // Simplify // 
        PowerExpand

Re[TwoFreq /. c -> 1] // ComplexExpand

Here is the plot. The sinusoidal envelope is an indication of the pulse-like behavior that will come when we add waves of more frequencies.

Table[
        Plot[Evaluate[
            Re[TwoFreq /. {c -> 1, ω -> 
                    15  π, Γp -> π}]], {z, -8/(2  π), 8/(2  π)}, 
          PlotRange -> {{-8/(2  π), 
                8/(2  π)}, {- .1,  .1}}, 
          PlotPoints -> 70], {t, 0, 4/3 - 1/9, 1/9}];

Superposition of two monochromatic waves 
(double-click plot to view as a function of time)

Three frequencies

Superposing waves of three frequencies, we get

ThreeFreq = Sum[
            Envelope[Δ] 
              EField[k[
                  1, ω + Δ], ω + Δ], {Δ, -π, π, π}] // 
          ComplexExpand // FullSimplify

Abs[ThreeFreq] // ComplexExpand // Simplify // 
      PowerExpand

Re[ThreeFreq] // ComplexExpand // FullSimplify // 
      PowerExpand

Table[
        Plot[Evaluate[
            Re[ThreeFreq /. {c -> 1, ω -> 
                    20  π, Γp -> π}]], {z, -12/(2  π), 12/(2  π)}, 
          PlotRange -> {{-12/(2  π), 
                12/(2  π)}, {- .3,  .3}}, 
          PlotPoints -> 70], {t, 0, 2 - 1/9, 1/9}];

Superposition of three monochromatic waves 
(double-click plot to view as a function of time)

Four frequencies

With four frequencies:

FourFreq = 
      Sum[Envelope[Δ] 
          EField[k[
              1, ω + Δ], ω + Δ], {Δ, -3 π/2, 3  π/2, π}]

Table[
        Plot[Evaluate[
            Re[FourFreq /. {c -> 1, ω -> 
                    20  π, Γp -> π}]], {z, -12/(2  π), 12/(2  π)}, 
          PlotRange -> {{-12/(2  π), 
                12/(2  π)}, {- .5,  .5}}, 
          PlotPoints -> 70], {t, 0, 6/3 - 1/9, 1/9}];

Superposition of four monochromatic waves 
(double-click plot to view as a function of time)

Note that we see periodic "pulse trains" because the frequencies that we use have small least common multiples. We could vary the frequencies to remove this periodic behavior, at the expense of destroying the smooth "pulse-like" appearance. Or, by adding more and more frequencies, we can obtain nonperiodic pulses with smooth shape, as shown in the next section.

Continuous frequency - Gaussian pulse

A real pulse contains every frequency within its bandwidth. Integration over all frequencies yields the formula

GaussianPulse[z_, 
        t_] = (Integrate[
            Envelope[Δ] 
              EField[k[
                  1, ω0 + Δ], ω0 + Δ], {Δ, -∞, ∞}, 
            Assumptions -> Γp > 0] // 
          FullSimplify)

The envelope is, not surprisingly, a Gaussian:

Abs[GaussianPulse[z, t]] // ComplexExpand // Simplify // 
      PowerExpand

Table[
        Plot[Evaluate[
            Re[GaussianPulse[z, 
                  t] /. {c -> 
                    1, Γp -> π, ω0 -> 
                    20  π}]], {z, -12/(2  π), 
            12/(2  π)}, 
          PlotRange -> {{-12/(2  π), 
                12/(2  π)}, {-1.5, 1.5}}, 
          PlotPoints -> 70], {t, 0, π, π/16}];

A Gaussian pulse
(double-click plot to view as a function of time)

Group velocity — general transparent medium

Consider two waves with slightly different frequencies

Table[(Subscript[ω, q]) = ωc + (((-1))^q) δ/2, {q, 1, 
          2}] // TableForm

Since the difference δ between the freqencies (Subscript[ω, 1]) and (Subscript[ω, 2]) is small, we can approximate the frequency dependence of k by a series about the central frequency ωc:

kSeries[ω_] = 
      Series[k[ω], {ω, ωc, 1}]

Here k'[ωc]=((d k)/(d ω))(Subscript[|, ωc]). The wavenumbers of the two waves are given in this approximation by

Table[(Subscript[k, q]) = kSeries[(Subscript[ω, q])] // Normal, {q, 1, 2}] // 
      TableForm

Table[Plot[Evaluate[{Re[EField[(Subscript[k, 1]), (Subscript[ω, 1])]], 
                    Re[EField[(Subscript[k, 2]), (Subscript[ω, 2])]]}/.{ωc -> 6,k[ωc] -> 6, k'[ωc]->1,δ -> 1}],{z, 0, 2  π},PlotStyle -> {RGBColor[0, 0, 0], 
                  RGBColor[1, 0, 0]},PlotRange -> {{0, 2  π}, {-2.5, 2.5}},PlotPoints -> 30],{t, 0, 2  π (17/20), 2  π/30}];

Two monochromatic waves with different frequencies and wavenumbers
(double-click plot to view as a function of time)

We see that the two waves are "in phase" at some regions of space and time, and are out of phase in others. Now let us look at the field corresponding to the sum of the two waves given by

EField[(Subscript[k, 1]), (Subscript[ω, 1])] + EField[(Subscript[k, 2]), (Subscript[ω, 2])]

Simplifying:

SumField = % // FullSimplify

The real field is given by

Re[SumField] // ComplexExpand // PowerExpand

and the envelope is given by

SumFieldEnv = Abs[SumField] // ComplexExpand // 
        PowerExpand

We assume here that k is real, i.e. that the medium is transparent.

Table[Plot[Evaluate[{Re[SumField], SumFieldEnv, -SumFieldEnv}/.{ωc -> 6,k[ωc] -> 4, k'[ωc]->1,δ -> 1. }],{z, 0, 8  π},PlotRange -> {{0, 8  π}, {-2.5, 2.5}},PlotPoints -> 50,PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0], 
                  RGBColor[1, 0, 0]}],{t, 0, 4  π (19/20), 4  π/20}];

Two superposed monochromatic waves in a medium, plotted with their envelope. Note that since the group velocity does not equal the phase velocity, a phase difference between the light field and the envelope is acumulated in time.
(double-click plot to view as a function of time)

To determine the velocity of the pulse, we consider the motion of a given point on the envelope, for example, a zero crossing: z[t] such that the |E[t]| = 0:

zConstantEnv = Solve[SumFieldEnv == 0, 
          z][[1, 1]]

The derivative of z with respect to time is the group velocity

Vg[k_] = D[z /. zConstantEnv, t]

as a function of ω, k is given by

kOfω[ω_] = k[n[ω], ω]

so the group velocity is given in terms of the index of refraction by

VgOfn[n_] = Vg[kOfω] // Simplify

This derivation is for the case of a transparent medium. As is pointed out in most optics texts, one has to be careful when applying the concept of group velocity to a medium with strong absorption, because pulse reshaping will occur as the pulse travels through the medium. We will look at this more closely next.

Pulse propagation in linear atomic media

Nontransparent media

We would like to determine when our formula for group velocity in transparent media will apply to nontransparent media. For a Gaussian pulse, the group velocity refers to the velocity of the peak of the pulse. Although inside an absorptive medium the pulse shape is distorted spatially, we can still measure the group velocity by comparing the times at which the pulse reaches a maximum at different points inside the medium. Our formula for the group velocity should hold as long as there is not too much distortion of the pulse shape, which would render the concept of a "group" velocity (the velocity at which the entire pulse travels) meaningless. Distortion is caused when resonant spectral components are absorbed much more strongly than off-resonant ones.

Absorption is governed by the imaginary part of k:

ImkMedium[ω_] = I Im[kMedium[ω]] // 
          ComplexExpand // FullSimplify

Let's consider a Gaussian input pulse. The spectrum of the input intensity is given by

Abs[Envelope[Δ]]^2 // ComplexExpand

The spectrum of the output intensity is given by

Transmission[Δ_] = Abs[Envelope[Δ] EField[ImkMedium[Δ + ω0], ω] /. Δ + ω0 -> ω0]^2 // ComplexExpand // PowerExpand

We assume here that ω ≃ ω0.

A Gaussian spectrum has a maxium only in the center. If we require this to be true for the output spectrum, we can obtain a condition for avoiding pulse reshaping. Setting the derivative with respect to Δ equal to zero, we can find the extrema:

Extrema = 
      Solve[D[Transmission[Δ], Δ] == 
          0, Δ]

The first solution is zero, and the second and third are imaginary. To find the condition for excluding the last two (i.e., requiring them to be imaginary), we solve for Δ = 0 in the fourth solution:

Solve[Δ == 0 /. 
        Extrema[[4]], α]

So we need \(TraditionalForm\`α\ z << \(Γ\^2\)/Γp\^2\) to avoid significant pulse reshaping. If this condition is satisfied, the formula we derived for group velocity in a transparent medium will remain approximately correct.

Group velocity

We now apply the expression for group velocity that we derived for transparent media to a linear atomic medium near resonance. Substituting our formula for the real part of the index of refraction in a atomic medium into the expression for group velocity, we obtain for the group velocity in a linear atomic medium near resonance

VgOfn[RenMedium]

Using our definition of ξ and assuming that we are directly on resonance gives

VgMedium = % /. ξ -> ξ[ωc] /. ωc -> ω0

We can see that on resonance, when α c < Γ the group velocity is greater than c and when α c > Γ the group velocity is negative. While at first glance, these results may seem unphysical (especially the second, which implies that the pulse moves backwards in the medium!) they actually have simple, consistent interpretations. We illustrate this below.

Setup for plotting

There is a fast-oscillating phase (e^(-i ω0 (t - z))) in the electric field that we don't need to include if we are plotting the envelope. We'll divide it out of the expression for the E-field:

Map[Factor, 
      ExpandAll[
        EFieldMedium[z, 
            t, ω0 + Δ] (E^(I ω0 (t - z/c)))], {2}]

Assuming that Γ << ω0, we can set the term containing Γ/ω0 to zero:

EFieldEnv[z_, 
        t_, Δ_] = % /. Γ/ω0 -> 
          0

This function plots a pulse traversing a medium

PlotField[FieldFunction_, {z1_, z2_}, 
        EnvelopePlotOptions_] := 
Plot[
        Abs[FieldFunction], {z, z1, z2}, 
        Evaluate[
          Join[EnvelopePlotOptions, {DisplayFunction -> Identity, 
              PlotDivision -> 2, PlotPoints -> 10}]]]

NPlotField[FieldFunction_, {z1_, z2_}, IntRange_, 
        EnvelopePlotOptions_] := 
Plot[
        Abs[NIntegrate[
            FieldFunction, {Δ, -IntRange, IntRange}, 
            PrecisionGoal -> 3, WorkingPrecision -> 12, 
            MaxRecursion -> 12, MinRecursion -> 2]], {z, z1, z2}, 
        Evaluate[
          Join[EnvelopePlotOptions, {DisplayFunction -> Identity, 
              PlotDivision -> 2, PlotPoints -> 10}]]]

Plots of pulse propagation

Superluminal pulse

We start with the case in which α c < Γ. Here α c/Γ = 0.5, and the group velocity is

VgMedium /. {α -> 1/(2  c), Γ -> 
          1}

For all the plots below, we set the light frequency to resonance. We see that even though the pulse is partly absorbed (i.e. its amplitude at the output is smaller than that at the input), the output pulse has approximately the same shape as the input. The peak of the output pulse is advanced relative to where the peak would have been hand there been no medium (dashed line), i.e., the pulse appears to travel faster than c inside the medium. No photons travel faster than c here, rather, the trailing edge of the pulse is absorbed more strongly than the leading edge creating the illusion of superluminal propagation.

PlotPacket[{0. , 25. , 35. , 60. }, {-30. , 90. , 
        3. }, {c -> 
          1, α ->  .5, Γp ->  .1, Γ -> 1. }]

A Gaussian pulse propagating through a linear medium.
(double-click plot to view as a function of time)

Negative group velocity

When α c < Γ, the group velocity is negative:

VgMedium /. {α -> 6/c, Γ -> 
          2}

This implies that the peak of the output pulse will exit the medium before the peak of the input pulse enters it:

PlotPacket[{0. , 25. , 35. , 60. }, {-30. , 90. , 
        3. }, {c -> 1, α -> 
          6. , Γp ->  .1, Γ -> 2. , 
        Region3Mag -> 8}]

A Gaussian pulse propagating through a linear medium.
(double-click plot to view as a function of time)

The absorption is now so large that we have to scale the amplitude of the output pulse  by (10^8) in order to see what is going on. Inside the medium the scaling is ramped up exponentially so that the plot remains smooth.

By adjusting the parameters, we can actually see the pulse moving backwards in the medium, although the (10^34)scaling required to observe the output makes this situation rather unphysical. We can also see that the pulse is compressed inside the medium.

VgMedium /. {α -> 16/c, Γ -> 
          4}

PlotPacket[{0. , 25. , 35. , 60. }, {-30. , 90. , 
        3. }, {c -> 1, α -> 
          16. , Γp ->  .1, Γ -> 4. , 
        Region3Mag -> 34}]

A Gaussian pulse propagating through a linear medium.
(double-click plot to view as a function of time)

Actually, the spatial pulse-like shape in the medium is an artifact of our scaling. The electric field at any given time is, in fact, a monotonically decreasing function of z. However, if we look at the value of the field at a given point as a function of time, we can observe pulse-like behavior even inside the medium.

Negative group velocities can be observed in experiments in linear optics. A situation much closer to the common experimental setup than the above case is one in which the pulse advance is much smaller than the total pulse width. Here the pulse is much wider than in the previous examples, but we still see that the output pulse peaks at the output before the input pulse peaks at the input.

PlotPacket[{0. , 25. , 35. , 60. }, {-3.  (10^10), 3.  (10^10), 
        3.  (10^9)}, {c -> 1, α -> 
          1. , Γp -> (10^(-10)), Γ -> 
          (10^(-9)), Region3Mag -> 2}]

A Gaussian pulse propagating through a linear medium.
(double-click plot to view as a function of time)

The pulse shape is not restricted to a simple Gaussian, as long as the pulse shape is analytic,  the entire pulse spectrum is within the region of linear anomalous dispersion, and \(TraditionalForm\`α\ z << \(Γ\^2\)/Γp\^2\), as derived above. Although complex pulse shapes can pass through a medium with group velocity greater than c, information does not travel faster than the speed of light, since the requirement that the pulse shape is analytic means that all the information in the pulse is contained in the front "tail" of the pulse.  The fact that the peak of the pulse moves faster than c does not mean that the information in the pulse does.

We illustrate the propagation of a non-Gaussian pulse with a double-peaked pulse (the sum of two Gaussians).

Setup for plotting

Plots

Here we plot a double-peaked pulse propagating through a linear medium. The entire pulse propagates at the same group velocity, and exits the medium advanced relative to where it would have been, had there been no medium.

Plot2Packets[{0. , 25. , 35. , 60. }, {-30. , 90. , 
        3. }, {c -> 
          1, α ->  .5, Γp ->  .1, Γ -> 1. , ts -> 25. , 
        Region3Mag -> 1}]

A double-peaked pulse propagating through a linear medium.
(double-click plot to view as a function of time)

Conclusion

We hope that this tutorial has demonstrated that:

All the "unusual" light propagation effects in this system follow directly from the form of the complex refractive index and elementary Fourier analysis.

Even though it may appear as if the peak of the light pulse is moving faster than c (under certain conditions), nothing is actually traveling that fast. The "superluminal" effects arise because the back of the pulse is absorbed more than the front, shifting the apparent peak of the pulse forward.

The usual formula for the group velocity 
	({{(Subscript[v, g]), =, Re((dω)/(d k))}, { , =, \(c\/\(Re\ n\  + \ ω\ Re\ \(d\ n\)\/\(d\ ω\)\)\)}}),
where ω is the light frequency, k is the complex wavenumber, and n is the complex refractive index, remains valid even in the cases where such velocity becomes infinite or negative, with all the seemingly strange consequencies to which this leads.

We would like to acknowldge helpful discussions with M. Zolotorev. Creation of this tutorial was supported by a UC Berkeley Classroom Technologies Grant (2002). The material is based in part upon work supported by the National Science Foundation (NSF) under Grant No. PHY-9733479 and by the Office of Naval Research (ONR) under Grant No. N00014-97-1-0214. Any opinions, findings, conclusions or recomendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF or the ONR.

R.W. Boyd and D.J. Gauthier, in Progress in Optics, edited by E. Wolf (Elsevier, Amsterdam, 2002), Vol. 43.

R.Y. Chiao, in Amazing light: a volume dedicated to Charles Hard Townes on his 80th birthday, edited by R.Y. Chiao (Springer, New York, 1996).

S. Chu and S. Wong, Physical Review Letters 48, 738 (1982).

A. Dogariu, A. Kuzmich, and L.J. Wang, Physical Review A (Atomic, Molecular, and Optical Physics) 63, 053806/1 (2001).

C.G.B. Garrett and D.E. McCumber, Physical Review A (General Physics) 1, 305 (1970).

D. J. Griffiths, Introduction to electrodynamics (Prentice Hall, Upper Saddle River, N.J., 1999).

A. Kuzmich, A. Dogariu, L.J. Wang, et al., Physical Review Letters 86, 3925 (2001).

L.D. Landau, E.M. Lifshits, and L.P. Pitaevskii, Electrodynamics of continuous media (Butterworth-Heinemann, Oxford England, 1995).

S.G. Lipson, H. Lipson, and D.S. Tannhauser, Optical physics (Cambridge University Press, Cambridge 1995).

A. B. Matsko, O. A. Kocharovskaya, Y. Rostovtsev, et al., in Advances in Atomic, Molecular, and Optical Physics, edited by B. Bederson (Academic Press, New York, 2001), Vol. 46, p. 191.

K.T. McDonald, American Journal of Physics 69, 607 (2001).

P. Sprangle, J.R. Penano, and B. Hafizi, Physical Review E (Statistical, Nonlinear, and Soft Matter Physics) 64, 026504/1 (2001).