# Classically accelerating solenoidal wave packets in two dimensions

###### Abstract

The wave function for a quantum mechanical particle in a circular box can be prepared in shape-preserving wave packets that rotate at constant angular speed around the center of the box. Some of these rotating wave packets correspond to classical trajectories undergoing uniform circular motion through the force-free interior of the box. This apparent violation of Ehrenfest’s theorem is resolved by considering the force exerted on the particle’s wave packet by the enclosing wall. Remarkably, the wave packet continues to rotate even after the wall potential is removed. We show that this force-free finite-energy rotating state actually corresponds to classical motion with constant velocity, again in agreement with Ehrenfest’s theorem. Even so, the classical angular momentum carried by the rotating states poses a conceptual challenge because it differs from its quantum mechanical angular momentum, and indeed can have the opposite sign.

###### pacs:

03.65.Ge, 03.50.-z, 42.25.FxIn 1979, Berry and Balazs reported the discovery of shape-preserving wave packets for quantum mechanical particles that translate with uniform acceleration in one dimension [1]. Taking the form of Airy functions, these wave packets appear to violate Ehrenfest’s theorem because they accelerate in the absence of any applied force. This conundrum is resolved by recognizing that an Airy wave packet is not square-integrable and thus is best interpreted as an ensemble of non-accelerating single-particle plane-wave states [1, 2]. Optical analogs to Airy wave packets have been realized in holographically-patterned laser beams, the temporal evolution of the quantum state being modeled through the spatial evolution of the light’s intensity profile [3, 4, 5]. In this case, the beam’s intensity profile translates along a parabolic path as it propagates, without otherwise distorting. The analogy between the spatial structure of a propagating light beam and the temporal evolution of a quantum mechanical wave packet reflects the homology of the paraxial wave equation with Schrödinger’s equation.

Here, we focus on an alternative class of shape-preserving wave packets in two dimensions that describe a particle undergoing uniform circular motion. The apparent radial acceleration of these self-rotating states has been described previously, although not its apparent violation of Ehrenfest’s theorem [6]. Here, we demonstrate that not all rotating states are classically accelerating in the sense identified by Berry and Balazs. We analytically demonstrate that the confined states’ radial acceleration arises from the wave packet’s interaction with the confining potential, and so is consistent with Ehrenfest’s theorem. We furthermore show that unconfined rotating states of the kind presented in [6] do not accelerate at all, nor do they even rotate. Counterintuitively, these “rotating” states actually undergo uniform translation. This surprising observation corrects a misinterpretation of the data in [6]. These states also have the surprising property that the classical angular momentum they carry differs from their quantum mechanical angular momentum, and indeed can have the opposite sign. We illustrate these properties through experimental realizations of analogous propagation-invariant laser modes projected with intermediate-plane holography [7].

The wave function, $\mathrm{\Psi}(\bm{r})$, of a nonrelativistic particle of mass $m$ moving in a two-dimensional circular box of radius $R$ can be expressed in polar coordinates, $\bm{r}=(r,\varphi )$, in terms of eigenfunctions of the force-free time-independent Schrödinger equation,

$$-\frac{{\mathrm{\hslash}}^{2}}{2m}{\nabla}^{2}\mathrm{\Psi}(\bm{r})=E\mathrm{\Psi}(\bm{r}),$$ | (1a) | ||

with polar Laplacian | |||

$${\nabla}^{2}=\frac{{\partial}^{2}}{\partial {r}^{2}}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{1}{{r}^{2}}\frac{{\partial}^{2}}{\partial {\varphi}^{2}},$$ | (1b) | ||

subject to the boundary condition | |||

$$\mathrm{\Psi}(R,\varphi )=0.$$ | (1c) |

Equation (1) is satisfied by the Bessel wave functions, $|n,\nu \u27e9$, whose spatial representation is

${\mathrm{\Psi}}_{n,\nu}(\bm{r})$ | $=\u27e8\bm{r}|n,\nu \u27e9$ | (2a) | ||

$={A}_{n,\nu}{J}_{n}\left({j}_{n,\nu}{\displaystyle \frac{r}{R}}\right){e}^{in\varphi},$ | (2b) | |||

where ${J}_{n}(x)$ is a Bessel function of the first kind of order $n$, and where ${j}_{n,\nu}$ is its $\nu $-th zero. The prefactor | ||||

$${A}_{n,\nu}={\left[{\pi}^{1/2}R{J}_{n+1}({j}_{n,\nu})\right]}^{-1}$$ | (2c) |

ensures that the corresponding probability density

$${\rho}_{n,\nu}(\bm{r})={\left|{\mathrm{\Psi}}_{n,\nu}(\bm{r})\right|}^{2},$$ | (3) |

is properly normalized. The Bessel states then are orthonormal:

$$\u27e8{n}^{\prime},{\nu}^{\prime}|n,\nu \u27e9={\delta}_{n,{n}^{\prime}}{\delta}_{\nu ,{\nu}^{\prime}}.$$ | (4) |

Their eigenenergies,

$${E}_{n,\nu}=\frac{{\mathrm{\hslash}}^{2}{j}_{n,\nu}^{2}}{2m{R}^{2}},$$ | (5) |

depend on both the azimuthal quantum number, $n$, and the radial quantum number, $\nu $. The associated frequency, ${\omega}_{n,\nu}={E}_{n,\nu}/\mathrm{\hslash}$, establishes the eigenstates’ time evolution,

$${\mathrm{\Psi}}_{n,\nu}(\bm{r},t)=\u27e8\bm{r}|n,\nu \u27e9{e}^{-i{\omega}_{n,\nu}t}.$$ | (6) |

The states described by Eqs. (2) and (6) are analogous to optical Bessel beams [8, 9], with the time in Eq. (6) serving as an analog to the light wave’s axial coordinate. Like their optical counterparts, Bessel wavefunctions carry angular momentum with expectation value

$$\u27e8{L}_{z}\u27e9=-i\mathrm{\hslash}\int {\mathrm{\Psi}}_{n,\nu}^{\ast}(\bm{r})\frac{\partial}{\partial \varphi}{\mathrm{\Psi}}_{n,\nu}(\bm{r}){d}^{2}r=n\mathrm{\hslash}$$ | (7) |

that depends on $n$, but not on $\nu $. For optical Bessel beams, this orbital angular momentum is a classical property of the electromagnetic field [10, 11] that also is a quantum mechanical property of the individual photons [12]. For the particle in a circular box, it is strictly a quantum mechanical property. The particle’s probability density, ${\rho}_{n,\nu}(\bm{r})$, is independent of time, which means that the particle is stationary in the classical sense and therefore carries no classical angular momentum. Similar discrepancies between the classical and quantum mechanical angular momentum have been noted for Landau states in free electron beams [13].

Although individual Bessel eigenmodes are time-invariant, some of their superpositions have probability densities that rotate at a uniform rate without otherwise distorting [14, 15, 16, 6, 17]. Some of these rotating wave packets constitute accelerating states in the sense that the expectation value of the particle’s position traces out an accelerating trajectory [6]. These are not simply related to two-dimensional Airy states or to related Matthieu and Weber states [18] or to their generalizations [19, 20, 21], and thus constitute a distinct class of accelerating states in two dimensions.

Minimal examples of rotating wave packets can be constructed by superposing two Bessel states:

$$\mathrm{\Psi}(\bm{r},t)=\frac{1}{\sqrt{2}}\left[{\mathrm{\Psi}}_{n,\nu}(\bm{r},t)+{\mathrm{\Psi}}_{{n}^{\prime},{\nu}^{\prime}}(\bm{r},t)\right].$$ | (8) |

Similar superpositions describe helicon waves in plasmas [14], which also are known as whistlers. For clarity, we arrange indices so that $\mathrm{\Delta}n={n}^{\prime}-n>0$. This superposition’s probability density

$$\begin{array}{c}\rho (\bm{r},t)=\frac{1}{2}\left[{\rho}_{n,\nu}(\bm{r})+{\rho}_{{n}^{\prime},{\nu}^{\prime}}(\bm{r})\right]+\hfill \\ \hfill {\left[{\rho}_{n,\nu}(\bm{r}){\rho}_{{n}^{\prime},{\nu}^{\prime}}(\bm{r})\right]}^{1/2}\mathrm{cos}\left(\mathrm{\Delta}n[\varphi -\mathrm{\Omega}t]\right),\end{array}$$ | (9) |

rotates around the origin with an angular frequency

$$\mathrm{\Omega}=\frac{\mathrm{\hslash}}{2m{R}^{2}}\frac{{j}_{{n}^{\prime},{\nu}^{\prime}}^{2}-{j}_{n,\nu}^{2}}{\mathrm{\Delta}n}.$$ | (10) |

Aside from this rotation, the state neither broadens nor otherwise distorts. The resulting periodic recurrence differs from the breathing modes identified in generalized Airy states [5]. Instead, it closely resembles the discrete propagation invariance of rotating optical modes [22], particularly solenoidal beams [16, 6]. For this reason, we refer to the rotating wave functions described by Eq. (8) as solenoidal states. Figure 1 shows the time evolution of two illustrative examples.

Most solenoidal wave packets are not accelerating states in the sense identified by Berry and Balazs. Those with $\mathrm{\Delta}n>1$, such as the example in Fig. 1(a), are symmetric about the origin; the expectation value of the particle’s position coincides with the center of the box. Classically, therefore, such states resemble their constituent Bessel states in that the particle remains motionless at the origin even as its probability density rotates.

Solenoidal states with $\mathrm{\Delta}n=1$ are asymmetric, as can be seen in Fig. 1(b). The expectation value for the particle’s position,

$$\u27e8\bm{r}(t)\u27e9=\alpha R\left[\mathrm{cos}(\mathrm{\Omega}t)\widehat{x}+\mathrm{sin}(\mathrm{\Omega}t)\widehat{y}\right],$$ | (11) |

undergoes uniform circulation motion at angular frequency $\mathrm{\Omega}$ and radius $\u27e8r\u27e9=\alpha R$, where

$\alpha $ | $={\displaystyle \frac{{\int}_{0}^{1}{J}_{n+1}({j}_{n+1,{\nu}^{\prime}}x){J}_{n}({j}_{n,\nu}x){x}^{2}\mathit{d}x}{{J}_{n+2}({j}_{n+1,{\nu}^{\prime}}){J}_{n+1}({j}_{n,\nu})}}$ | (12a) | ||

$={\displaystyle \frac{2{j}_{n+1,{\nu}^{\prime}}{j}_{n,\nu}}{{({j}_{n+1,{\nu}^{\prime}}^{2}-{j}_{n,\nu}^{2})}^{2}}}.$ | (12b) |

This seems remarkable because no force acts on the particle within the box. Such force-free acceleration appears to violate Ehrenfest’s theorem,

$$\frac{{d}^{2}\u27e8\bm{r}(t)\u27e9}{d{t}^{2}}=\frac{1}{m}\u27e8\bm{F}(\bm{r}(t))\u27e9,$$ | (13) |

which relates the expectation value of the particle’s acceleration to the expectation value of the force acting on the particle. In the present case, $\bm{F}(\u27e8\bm{r}(t)\u27e9)=0$ in the force-free region within the box, but

$$\frac{{d}^{2}\u27e8\bm{r}(t)\u27e9}{d{t}^{2}}=-\alpha R{\mathrm{\Omega}}^{2}\widehat{r}.$$ | (14) |

The apparent discrepancy can be explained because $\bm{F}(\u27e8\bm{r}(t)\u27e9)\ne \u27e8\bm{F}(\bm{r}(t))\u27e9$.

The integrability of the solenoidal wave packets comes at the cost of applying the boundary condition from Eq. (1c) at $r=R$. The confined particle exerts a pressure on the wall,

$P$ | $=-{\displaystyle \frac{1}{2\pi R}}{\displaystyle \frac{dE}{dR}}$ | (15a) | ||

$={\displaystyle \frac{{\mathrm{\hslash}}^{2}}{4\pi m{R}^{4}}}\left({j}_{n+1,{\nu}^{\prime}}^{2}+{j}_{n,\nu}^{2}\right).$ | (15b) |

By Newton’s third law, the wall exerts a complementary force on the wave packet that is directed radially inward. When averaged over angles, the net force acting on the particle located at $\u27e8\bm{r}\u27e9$ is

$$\u27e8\bm{F}(\bm{r})\u27e9=-\beta RP\widehat{r},$$ | (16a) | |||

where | ||||

$\beta $ | $=\underset{\u03f5\to 0}{lim}{\displaystyle \frac{{{\int}_{0}^{2\pi}\rho (\bm{r},t)\mathrm{cos}(\varphi -\mathrm{\Omega}t)\mathit{d}\varphi |}_{r=R-\u03f5}}{{{\int}_{0}^{2\pi}\rho (\bm{r},t)\mathit{d}\varphi |}_{r=R-\u03f5}}}$ | (16b) | ||

$={\displaystyle \frac{2{j}_{n+1,{\nu}^{\prime}}{j}_{n,\nu}}{{j}_{n+1,{\nu}^{\prime}}^{2}+{j}_{n,\nu}^{2}}}.$ | (16c) |

Comparison with Eq. (14) confirms that Ehrenfest’s theorem is satisfied: the force responsible for the particle’s classical circular motion is exerted by the wave function’s interaction with the bounding wall.

In principle, accelerating solenoidal states can be prepared as superpositions of Bessel beams without confining boundary conditions. Such unconfined states rotate without a central force, and so appear to violate Ehrenfest’s theorem. Because Bessel states are not square-integrable, however, they are best interpreted as superpositions of plane-wave states [1]. The rotation of solenoidal wave packets then represents the ensemble-average of multiple particles’ non-accelerating trajectories.

The possibility of Ehrenfest violations arises again for normalizable approximations to unconfined solenoidal wave packets that are prepared by truncating $\mathrm{\Psi}(\bm{r},t)$ at $r=R$. These truncated wave packets still appear to rotate, and they evolve under truly force-free conditions. To illustrate this phenomenon, Fig. 2 presents optical analogs of truncated solenoidal states. We prepare these solenoidal beams of light by using intermediate-plane holograms [7] to convert the linearly polarized Gaussian TEM${}_{00}$ beam from a solid state laser (Coherent Verdi 5W) into a superposition of helical Bessel modes of the form described by Eq. (8). The mode-converting holograms are imprinted on the beam with a phase-only liquid crystal spatial light modulator (Hamamatsu X10468-16).

Figure 2(a) and Fig. 2(b) show volumetric reconstructions of two solenoidal laser beams with $\mathrm{\Delta}n=1$. The data for these reconstructions were obtained by translating a video camera (NEC TI-324AII) along the optical axis and combining the resulting stack of images. Each reconstruction shows three complete cycles of shape-preserving propagation. The solenoid in Fig. 2(a) has a right-handed twist while that in Fig. 2(b) is left-handed. Each volumetric reconstruction is paired with three transverse slices from the planes labeled (i), (ii) and (iii) in the renderings. These slices show three stages in the intensity distributions’ rotation about the optical axis at $120\mathrm{\xb0}$ intervals. Additional planes in the renderings correspond to each of the three complete rotations captured over the course of $1\mathrm{m}$ of propagation. The full non-diffracting range of these beams extends beyond $2\mathrm{m}$.

The two-state superpositions discussed so far are not the only accelerating wave packets. Figure 2(c) shows a four-state superposition designed to better localize the wave packet as it spirals around the origin. In this case, there can be no doubt that the intensity maximum precesses steadily around the optical axis. Incorporating higher-order terms in the superposition can further improve localization [16] at the cost of reducing diffraction efficiency.

Figure 3(a) shows a representative simulation of ${\bm{r}}^{\ast}(t)$ superimposed on a snapshot of the initial distribution, $\rho (\bm{r},0)$. Figure 3(b) shows the corresponding experimental measurement of ${\bm{r}}^{\ast}(z)$. The angular position of the peak, ${\theta}^{\ast}(t)$, advances uniformly, as can be seen in Fig. 3(c). For clarity, we have scaled the simulation time to best superimpose the temporal evolution of the simulation data on the spatial evolution of the experimental data, ${\theta}^{\ast}(z)$. The shaded experimental error range was estimated by measuring variations in the apparent position of a non-rotating $n=0$ Bessel state using the same protocol as we use for the rotating superposition states. These variations arise from mechanical vibrations in the translation stage used to move the camera and from the influence of Gaussian noise in the measured images. Discrepancies between the measured and simulated results also can be attributed to speckle and other imperfections in the projected beam’s mode structure.

This apparent contradiction of the Ehrenfest theorem is resolved by tracking the mean particle position, $\u27e8\bm{r}(t)\u27e9$, which also is plotted in Fig. 3. Results for $\u27e8x(t)\u27e9$ and $\u27e8x(z)\u27e9$ are plotted in Fig. 3(d), with the simulation time again scaled as in Fig. 3(c). In both simulation and experiment, the classical trajectory of the unconfined rotating wave packet actually translates steadily away from the beam’s axis with an impact parameter, $b$ set by the average position at $t=0$. The experimental center of brightness closely follows a linear trajectory along $b=0.11\pm 0.01\mathrm{mm}$ with a standard deviation of $\mathrm{\Delta}x=0.02\mathrm{mm}$. The experimental trajectory’s impact parameter is shifted by $\mathrm{\Delta}b=0.05\mathrm{mm}$ relative to the simulation value due to imperfections in the projected mode’s shape.

The unconfined solenoidal mode’s linear translation
arises because the truncated
wave packet diffracts beyond $r=R$ by precisely
the amount needed to conserve momentum
under force-free conditions
^{1}^{1}
1
See Supplemental Material at [URL] for
an animated visualization that shows
how diffraction influences the structure
of solenoidal states. This video presents
experimental data for the spatial development
of an optical solenoid, and
a simulation of the time evolution of the
corresponding quantum mechanical wave packet..
The nature of the state’s time
evolution is masked because diffraction has
little apparent influence on the wave packet’s
structure at early times, particularly near
the center of the system.
Remarkably, this means that freely propagating
solenoidal laser modes neither
accelerate [6] nor even rotate [14]; they are
translating states.

Although the confined particle’s classical acceleration can be accounted for by the influence of boundary conditions, its rate of circulation is less straightforward to interpret. The classical angular momentum carried by an accelerating solenoidal wave packet is

${L}_{z}^{(c)}$ | $=m{\alpha}^{2}{R}^{2}\mathrm{\Omega}$ | (17a) | ||

$=2{\displaystyle \frac{{j}_{n+1,{\nu}^{\prime}}^{2}{j}_{n,\nu}^{2}}{{({j}_{n+1,{\nu}^{\prime}}^{2}-{j}_{n,\nu}^{2})}^{3}}}\mathrm{\hslash},$ | (17b) |

which differs from the state’s quantum mechanical angular momentum,

$$\u27e8{L}_{z}\u27e9=\left(n+\frac{1}{2}\right)\mathrm{\hslash}.$$ | (18) |

Indeed, ${L}_{z}^{(c)}$ and $\u27e8{L}_{z}\u27e9$ can have opposite signs depending on the choice of radial quantum numbers $\nu $ and ${\nu}^{\prime}$. The solenoidal states represented in Figs. 2(a) and 2(b), for example, have opposite classical angular momentum even though they carry the same quantum mechanical angular momentum. Such discrepancies are surprising because they are not generically observed in quantum mechanical systems carrying angular momentum [24, 25].

Unconfined solenoidal states also carry classical angular momentum

$${L}_{z}^{(c)}=m\left[\u27e8\bm{r}(0)\u27e9\times \frac{d\u27e8\bm{r}(t)\u27e9}{dt}\right]\cdot \widehat{z}$$ | (19) |

that is equal to the confined value from Eq. (17) and generally differs from the quantum-mechanical value, Eq. (18). Analogous discrepancies have been been reported in holographically-patterned electron beams [13].

We have demonstrated that the radial acceleration associated with confined solenoidal states does not carry over to the unconfined states that are experimentally realizable with light beams [14, 16, 6]. Instead, the small amount of diffraction that occurs for these states is precisely what is required to transform their acceleration into uniform translation. This does not account for the observed discrepancy between classical and quantum mechanical angular momentum, which therefore remains an outstanding challenge.

This work was supported primarily by the National Science Foundation under award number DMR-1305875 and in part by the MRSEC program of the NSF through award number DMR-1420073. The authors are grateful for helpful conversations with David Pine, Paul Chaikin and Aaron Yevick.

## References

- [1] M. V. Berry and N. L. Balazs. “Nonspreading wave packets.” Am. J. Phys. 47, 264–267 (1979).
- [2] L. E. Ballentine, Y. Yang and J. P. Zibin. “Inadequacy of Ehrenfest’s theorem to characterize the classical regime.” Phys. Rev. A 50, 2854 (1994).
- [3] G. A. Siviloglou, J. Broky, A. Dogariu and D. N. Christodoulides. “Observation of accelerating Airy beams.” Phys. Rev. Lett. 99, 213901 (2007).
- [4] G. A. Siviloglou and D. N. Christodoulides. “Accelerating finite energy Airy beams.” Opt. Lett. 32, 979–981 (2007).
- [5] I. Kaminer, R. Bekenstein, J. Nemirovsky and M. Segev. “Nondiffracting accelerating wave packets of Maxwell’s equations.” Phys. Rev. Lett. 108, 163901 (2012).
- [6] C. Vetter, T. Eichelkraut, M. Ornigotti and A. Szameit. “Generalized radially self-accelerating helicon beams.” Phys. Rev. Lett. 113, 183901 (2014).
- [7] A. Mondal, Yevick, Aaron, L. C. Blackburn, N. Kanellakopoulos and D. G. Grier. “Projecting non-diffracting waves with intermediate-plane holography.” Opt. Express 26, 3926–3931 (2018).
- [8] J. Durnin. “Exact-solutions for nondiffracting beams. 1. The scalar theory.” J. Opt. Soc. Am. A 4, 651–654 (1987).
- [9] J. Durnin, J. J. Miceli, Jr and J. H. Eberly. “Diffraction-free beams.” Phys. Rev. Lett. 58, 1499–1501 (1987).
- [10] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw and J. P. Woerdman. ‘‘Orbital angular-momentum of light and the transformation of Laguerre-Gaussian laser modes.” Phys. Rev. A 45, 8185–8189 (1992).
- [11] L. Allen, M. J. Padgett and M. Babiker. “The orbital angular momentum of light.” Prog. Opt. 39, 291–372 (1999).
- [12] J. Leach, M. J. Padgett, S. M. Barnett, S. Franke-Arnold and J. Courtial. “Measuring the orbital angular momentum of a single photon.” Phys. Rev. Lett. 88, 257901 (2002).
- [13] P. Schattschneider, T. Schachinger, M. Stöger-Pollach, S. Löffler, A. Steiger-Thirsfeld, K. Y. Bliokh and F. Nori. “Imaging the dynamics of free-electron Landau states.” Nat. Comm. 5, 1–6 (2014).
- [14] C. Paterson and R. Smith. “Helicon waves: propagation-invariant waves in a rotating coordinate system.” Opt. Commun. 124, 131–140 (1996).
- [15] R. Vasilyeu, A. Dudley, N. Khilo and A. Forbes. “Generating superpositions of higher-order Bessel beams.” Opt. Express 17, 23389–23395 (2009).
- [16] S.-H. Lee, Y. Roichman and D. G. Grier. “Optical solenoid beams.” Opt. Express 18, 6988–6993 (2010).
- [17] C. Schulze, F. S. Roux, A. Dudley, R. Rop, M. Duparré and A. Forbes. “Accelerated rotation with orbital angular momentum modes.” Phys. Rev. A 91, 043821 (2015).
- [18] P. Zhang, Y. Hu, T. Li, D. Cannan, X. Yin, R. Morandotti, Z. Chen and X. Zhang. “Nonparaxial Matthieu and Weber accelerating beams.” Phys. Rev. Lett. 109, 193901 (2012).
- [19] M. A. Bandres, M. A. Alonso, I. Kaminer and M. Segev. “Three-dimensional accelerating electromagnetic waves.” Opt. Express 21, 13917–13929 (2013).
- [20] I. D. Chremmos and N. K. Efremidis. “Nonparaxial accelerating Bessel-like beams.” Phys. Rev. A 88, 063816 (2013).
- [21] Q. Zhao, L. Gong and Y.-M. Li. “Shaping diffraction-free Lommel beams with digital binary amplitude masks.” Appl. Opt. 54, 7553–7558 (2015).
- [22] J. Tervo and J. Turunen. “Rotating scale-invariant electromagnetic fields.” Opt. Express 9, 9–15 (2001).
- [23]
- [24] W. Firth and D. Skryabin. “Optical solitons carrying orbital angular momentum.” Phys. Rev. Lett. 79, 2450–2453 (1997).
- [25] D. Skryabin and W. Firth. “Dynamics of self-trapped beams with phase dislocation in saturable Kerr and quadratic nonlinear media.” Phys. Rev. E 58, 3916–3930 (1998).