*Let's have a good test of Fnora's equation capabilities. Here is the first
chapter of my dissertation converted into QTML. *

The problem of solving the full vector Maxwell equations has been around for more than a century. With the exception of certain special cases, it is a problem to be avoided whenever possible. A full analytical or numerical solution can be quite difficult.

No, a general technique for fully solving Maxwell's equations is not to be found in this document. Only one step closer to this goal has been attempted. In particular, a new method for taking vector coupling into account has been developed for modeling a laser beam traveling in a guiding medium.

If one assumes that all nonlinearity can be included in the polarization density $\text{P}$ , the Maxwell wave equation [1] in Gaussian cgs units is

(1.1)

${\nabla}^{2}\text{E}-\nabla (\nabla \xb7\text{E})-\frac{\mu \epsilon}{{c}^{2}}{\partial}_{tt}\text{E}=\frac{4\pi \mu}{{c}^{2}}{\partial}_{tt}\text{P}.$

This can be written in matrix form as

(1.2)

$\left[\begin{array}{c}{\partial}_{yy}+{\partial}_{zz}\\ ----\\ -{\partial}_{y}{\partial}_{x}\\ ----\\ -{\partial}_{z}{\partial}_{x}\end{array}\begin{array}{c}|\\ |\\ |\\ |\\ |\end{array}\begin{array}{c}-{\partial}_{x}{\partial}_{y}\\ ----\\ {\partial}_{xx}+{\partial}_{zz}\\ ----\\ -{\partial}_{z}{\partial}_{y}\end{array}\begin{array}{c}|\\ |\\ |\\ |\\ |\end{array}\begin{array}{c}-{\partial}_{x}{\partial}_{z}\\ ----\\ -{\partial}_{y}{\partial}_{z}\\ ----\\ {\partial}_{xx}+{\partial}_{yy}\end{array}\right]\left[\begin{array}{c}{\text{E}}_{x}\\ --\\ {\text{E}}_{y}\\ --\\ {\text{E}}_{z}\end{array}\right]-\frac{\mu \epsilon}{{c}^{2}}{\partial}_{tt}\left[\begin{array}{c}{\text{E}}_{x}\\ --\\ {\text{E}}_{y}\\ --\\ {\text{E}}_{z}\end{array}\right]=\frac{4\pi \mu}{{c}^{2}}{\partial}_{tt}\left[\begin{array}{c}{\text{P}}_{x}\\ --\\ {\text{P}}_{y}\\ --\\ {\text{P}}_{z}\end{array}\right].$

The problem comes from the cross derivatives. First, these cross terms couple the vector components together, which prevents solving for a single component. Secondly, consider a laser beam traveling in the $z$ direction. It would simplify matters considerably if the equation of motion could be written in such a way that ${\partial}_{z}\text{E}$ is equal to a transverse operator acting on$\text{E}$ plus source terms. Unfortunately, the cross terms involving rule this out.

This problem is particularly apparent for the case of finding the modes of a dielectric waveguide. There, $\epsilon $ and are functions of and $only,\text{and}$$\text{exp}(i{\beta}_{n}z\text{\hspace{0.17em}}-\text{\hspace{0.17em}}i\omega t)$ , where ${\beta}_{n}$$n$ th mode. Replacing ${\partial}_{t}$ with $-i\omega $ and ${\partial}_{z}$ with in (1.2) results in a nice eigenvalue problem for $\omega $ . What comes next is to find a matrix representation for the various operators and send the resulting matrix to the appropriate EISPACK subroutine. But usually it is $\omega $ that is known and ${\beta}_{n}$ that needs to be found. To find ${\beta}_{n}$ by solving an inverse problem is quite inconvenient. One could set ${\beta}_{n}$ to some trial value in order to see if an eigenvalue near the known value of $\omega $ appears. Then an iteration scheme could be used to make that eigenvalue closer to the known $\omega $ . For a single-mode fiber, this is feasible but time-consuming. For a representation with numerical values of .eqni E .eqni on a mere 10x10 grid, a 300x300 matrix has to be solved for each guess at ${\beta}_{n}$ .

The usual way around these problems is to assume $\nabla (\nabla \xb7\text{E})=0$ . Since

(1.3)

$\nabla \xb7\text{E}=-\frac{4\pi (\nabla \chi \xb7\text{E}+\nabla \xb7\text{P})}{1+4\pi \chi},$

this holds very well in homogeneous media with nonlinearity small enough so that $\nabla (\nabla \xb7\text{P})$ is much smaller than$\left({\omega}^{2}/{c}^{2}\right)\text{P}$ . It also holds for inhomogeneous media as long as the variation of $\chi $ is sufficiently gradual; i.e., the weakly guided approximation. This encompasses many cases of physical interest.

The resulting equation to be solved is

(1.4)

${\partial}_{zz}\text{E}+{\nabla}_{T}^{2}\text{E}-\frac{\mu \epsilon}{{c}^{2}}{\partial}_{tt}\text{E}=\frac{4\pi \mu}{{c}^{2}}{\partial}_{tt}\text{P}.$

There is no explicit coupling between the components. It is now possible to solve for just one component of the electric field if there is no significant vector coupling in the constitutive relation for $\text{P}$ . There are two basic methods of solving (1.4) for a laser beam traveling in the $z$ direction. The first is to write $\text{E}$ in terms of a carrier and envelope, $\text{E}=\stackrel{\text{~}}{\text{E}}(x,y,z,t)\text{exp}(ikz-i\omega t)$ where $k=\sqrt{{\epsilon}_{0}{\mu}_{0}}\omega /c$ and where${\epsilon}_{0}$ and ${\mu}_{0}$ are the constant background parts of the permittivity and permeability. Then one assumes that the envelope varies slowly compared to the oscillations of the carrier so that ${\partial}_{tt}\stackrel{\text{~}}{\text{E}}$ , ${\partial}_{tt}\stackrel{\text{~}}{\text{P}}$ , ${\partial}_{t}\stackrel{\text{~}}{\text{P}}$ , and${\partial}_{zz}\stackrel{\text{~}}{\text{E}}$ can be dropped. The last deletion is the paraxial approximation, which assumes that all rays are nearly parallel to the $z$ axis. For constant $\epsilon $ the remaining time derivative of $\text{E}$ can be eliminated by changing to a comoving frame. The result is a separate equation for each time slice, the only time coupling being through $\text{P}$ . For many situations of interest it is reasonable to entirely ignore the time variations of the envelope. The resulting equations can be solved by numerous techniques such as finite-difference or transform techniques. For an example, see Adams, Cantrell, Lee [2] or Crenshaw and Cantrell [3].

The second method is to incorporate all nonlinearity into a local$\epsilon $ . Then the time independent equation for one component takes the form

(1.5)

$({\partial}_{zz}+T)\psi =0,$

where

(1.6)

$T={\nabla}_{T}^{2}+\frac{{\omega}^{2}\mu \epsilon}{{c}^{2}}.$

This can be factored as $\left(({\partial}_{z}+i\sqrt{T})({\partial}_{z}-i\sqrt{T})+i[{\partial}_{z},\sqrt{T}]\right)\psi =0$ . If $[{\partial}_{z},\sqrt{T}]$ is small enough to be ignored, we have

(1.7)

${\partial}_{z}\psi =\pm i\sqrt{T}\psi $

This can be formally solved for the forward-traveling wave as

(1.8)

${\psi}_{i+1}=\text{exp}\left[i\underset{0}{\overset{\Delta z}{\int}}\sqrt{T}dz\right]{\psi}_{i}\approx \text{exp}\left[i\Delta z\sqrt{T}\right]{\psi}_{i}.$

This is the idea used for one version of the Beam Propagation Method developed by Feit and Fleck [4,5]. The exponentiation of the square root of the transverse operator is accomplished by using a combination of Fourier-transform techniques and operator splitting. The possibility of back-scattered waves is taken into account by the fact that $\text{exp}(i\Delta zQ)$ is not unitary. Energy can get lost from step to step.

The above is a powerful method for getting rid of the paraxial approximation. But the scalar approximation remains.

For many problems of physical (and economic) interest, the scalar approximation is questionable. Modern communication technology is moving rapidly toward single-mode optical fibers as the transmission medium of choice. For a single-mode fiber $\chi $ varies significantly over the course of a wavelength. The same is true for the semiconductor lasers used to generate the signals. An entire optical counterpart to electronics is in development. The speed and bandwidth of optical devices outstrip those of their electronic counterparts [6]. Any improved model for such devices could be quite useful. As Dr. Lars Thylen [needs accent] recently stated, "A simple BPM type formalism for solving the vector wave equation is a very interesting research topic [7]."

Thus, this work is centered around designing a new numerical method for solving the vector wave equation for guided-wave optical devices. As will be shown in the next chapter, the use of the Hertz vectors, with appropriate gauge conditions, allows the complete elimination of the cross terms with respect to $z$ for the case of a transversely varying linear dielectric. A cross derivative with respect to $z$ shows up peripherally for the case of nonlinear source terms from the conversion from the Hertz vector potential to the physical field.

Chapter 3 outlines the numerical method used. Since the Hertz vector equations are more complex than the scalar wave equation, the trick of exponentiating the square root of the transverse operator is not used. Instead, finite differences are used. Probably the most interesting part of this chapter is a fast method for inverting a small class of pentadiagonal matrices. This allows the practical use of five-point finite differences to estimate ${\partial}_{xx}$ and ${\partial}_{yy}$ by an implicit method. Chapter 4 shows the results of some tests of the numerical method against a couple of analytic solutions. Chapter 5 shows the results of applying the resulting code to the problem of self-focusing, and Chapter 6 summarizes the overall results. The computer source program can be found in the appendix.

Since finding the modal properties of optical fibers is interesting from an engineering standpoint, but is perhaps less interesting from the point of view of pure physics, this new algorithm has been applied to the study of self-focusing.

One of the simplest models in nonlinear optics is an index of refraction or dielectric permittivity which is linearly dependent on the intensity of light.

(1.9)

$n={n}_{0}+{n}_{2}{\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}},$

or

(1.10)

$\epsilon ={\epsilon}_{0}+{\epsilon}_{2}{\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}}.$

If ${n}_{2}{\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}}$ is much less than ${n}_{0}$

the two models are effectively equivalent. In actual physical systems modeled by the above equations, it may actually be the intensity averaged over some period of time for which $n$ depends. As long as this time is short compared to the time variation of the envelope, the above model can be used allowing the use of a time-independent numerical model. Another way of stating (1.10) is

(1.11)

$\stackrel{\text{~}}{\text{P}}={\chi}^{\left(3\right)}({\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}})\stackrel{\text{~}}{\text{E}},$

where

(1.12)

${\chi}^{\left(3\right)}=\frac{{\epsilon}_{2}}{4\pi}\approx \frac{2{n}_{0}{n}_{2}}{4\pi \mu}.$

The last term assumes that $\mu $ is constant. In the literature it is not uncommon for (1.9) to be the initially stated model, while it is actually (1.10) that actually goes into the equation to be solved [7,8]. In this work it is (1.11) that shall be used although parameters will often be stated as ${n}_{2}$ to relate to other works. The assumption to keep in mind is that ${n}_{2}$ is not used directly, but the ${\chi}^{\left(3\right)}$ which is equivalent in the small field limit.

If ${\chi}^{\left(3\right)}$ has a positive real part, there is an effective increase in $\epsilon $ and therefore $n$ in areas where the intensity is high. This increase in the index of refraction acts like a convex lens which focuses the beam. As the beam focuses, the intensity increases, causing an increase in the lens effect for positive feedback.

Self-focusing was first predicted by Askar'yan in 1961 [10] and first observed by Hercher in 1964 [11]. The observation was made that a high-power Q-switched laser beam traveling through glass produced long straight filaments of damage a few microns across and up to over two centimeters long. Chiao, Garmire and Townes explained these filaments as self-trapped light [9]. That is, by letting $\epsilon ={\epsilon}_{0}+{\epsilon}_{2}|\text{E}{|}^{2}$ (scalar approximation), they obtained transverse profiles in which the induced change in $\epsilon $ compensated exactly for diffraction. The light forms its own waveguide. These profiles were found directly, not through numerical propagation.

These filaments present several problems. First, they put a limit on how much power can be sent through otherwise transparent optical elements. For each local maximum in a beam of high enough power a damage filament can be produced. Since this is such a low-order effect, one can expect it to show up through many different physical mechanisms, much as harmonic motion is ubiquitous in nature. Self-focusing by many different mechanisms has been observed in solids, liquids and gases [8]. It is a limiting factor in the design of high-power lasers [8]. Filamentation has been observed in plasmas and is a major impediment to inertial-confinement fusion [12,13]. Filamentation has even been observed using a mere 5 mW helium-neon laser beam in Chinese green tea [14].

The second problem is in the theory. With the scalar and paraxial approximations, the power and therefore the change in $\epsilon $ is inversely proportional to the square of the width of the beam, but so is diffraction [8]. Thus, the positive feedback mentioned earlier goes on indefinitely. A beam which starts to self-focus will continue to focus all the way down to a singularity. A beam which is below threshold will continue to diffract. So the self-trapped solutions should rarely show up in experiment. For a real physical system other nonlinear effects come into play before the field can focus to a point. Actually, there is a middle ground where a beam can start to focus and then diffract when only near the peak of, say, a Gaussian, does the change in $\epsilon $ counteract the diffraction term. Outside this area the beam diffracts, providing an "absorbing boundary" for the inner area. But this mechanism produces a single maximum and then a die-off of on-axis power. It does not produce self-trapped filament solutions. This effect can be seen in computer studies by Dawes and Marburger [15]. In the same paper, Dawes and Marburger did get filament-like solutions by making the change in

$\epsilon $ saturable. That is, they let $\epsilon ={\epsilon}_{0}+{\epsilon}_{2}|\text{E}{|}^{2}/(1+(\text{E}/{\text{E}}_{s}{)}^{2})$ . This allows $\epsilon $ to approach an asymptotic value for high intensities which flattens out the inner part of the "lens," stopping further focusing.

Another possibility is that filamentation is not caused by self-trapping at all. In 1967 Lugovoi and Prokhorov proposed that the filaments are not produced by self-trapped waves but are an artifact of the time variation of the Q-switched laser pulses used in experiment [16]. The idea is that at different points in time of the pulse, the peak power is different though the transverse shape is the same. If the material relaxes quickly compared to the period of the pulse, each time slice can be considered separately. Since each slice has a different peak intensity, it will focus with a different focal length. The resulting string of self-focal points produces the damage filament. In 1968 Zverev, Maldutis and Pashkov provided some experimental confirmation of this explanation [17]. By using successively higher power pulses with a single mode laser, they were able to produce successively longer filaments. However, these filaments grew only in the direction towards the laser, which corresponds to the closer focus from the greater intensity peak.

Even though later writers have added to the verification of the moving focus theory [18-21], does this mean that self-trapping never occurs? The ability to produce induced waveguiding could possibly be useful somewhere. And the catastrophic focus solution is based on two approximations to Maxwell's equations which break down before the field could focus to a point.

In 1987 Feit and Fleck applied their non-paraxial version of the Beam Propagation Method to the problem of self-focusing [5]. (They actually use (1.9) in their model as opposed to (1.10).) With their non-paraxial method a limit gets placed on the focusing in that energy gets backscattered when the beam tries to focus smaller than a wavelength. With this model multiple foci are produced although the power dies off eventually due to backscattering. Such a succession of foci could produce filaments although for the parameters used the die-off occurs over much less than a centimeter.

Backscattering is obtained in their model through the non-unitary nature of $\text{exp}(i\sqrt{T}\Delta z)$ . Using their approximations, this is approximately

(1.13)

$\text{exp}\left(i\Delta z\left[\frac{{\nabla}_{T}^{2}}{({\nabla}_{T}^{2}+{k}^{2}{)}^{1/2}+k}+\chi \right]\right),\chi =k\left(\frac{n}{{n}_{0}}-1\right).$

To second order in $\Delta z$ this can be factored as

(1.14)

$\text{exp}\left(\frac{i\Delta z}{2}\left[\frac{{\nabla}_{T}^{2}}{({\nabla}_{T}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{k}^{2}{)}^{1/2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}k}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\chi \right]\right)\text{exp}\left(i\Delta z\chi \right)\text{exp}\left(\frac{i\Delta z}{2}\left[\frac{{\nabla}_{T}^{2}}{({\nabla}_{T}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{k}^{2}{)}^{1/2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}k}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\chi \right]\right).$

The square roots in the denominators of the outer exponentials become imaginary for transverse Fourier components of frequencies higher than $k$ . This produces energy loss when the field focuses to a sharp peak. One strike against this model is that the observed filaments are several wavelengths across, with the width depending on the medium.

The other approximation that breaks down as a beam focuses toward a point is the scalar approximation. For $\stackrel{\text{~}}{\text{P}}={\chi}^{\left(3\right)}({\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}})\stackrel{\text{~}}{\text{E}}$ , we have

(1.15)

$\nabla \prime \xb7\stackrel{\text{~}}{\text{P}}={\chi}^{\left(3\right)}\left[\stackrel{\text{~}}{\text{E}}\xb7\nabla \prime ({\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}})+({\stackrel{\text{~}}{\text{E}}}^{*}\xb7\stackrel{\text{~}}{\text{E}})\nabla \prime \xb7\stackrel{\text{~}}{\text{E}}\right].$

(The prime is a reminder that ${\partial}_{z}$ needs to be replaced by

${\partial}_{z}+ik$ since we are looking at the envelope. A look at (1.3) reveals that for $\nabla \prime \xb7\stackrel{\text{~}}{\text{E}}$ to be zero,

$\nabla \prime \xb7\stackrel{\text{~}}{\text{P}}$ must also be zero which is generally inconsistent with the above relation. For sufficiently tight self-focusing, the scalar approximation breaks down in a big way. In Chapter 5 I shall show what happens when the scalar approximation is not made.

[1] J. D. Jackson,
*Classical Electrodynamics*
(New York: John Wiley and Sons, 1975).

[2] D.R. Adams, C. D. Cantrell, and W. M. Lee, pp. 456-463 in
*Proceedings of the International Conference on Lasers '83,*
R. C. Powell, ed. (STS Press, 1983).

[3] M. E. Crenshaw and C. D. Cantrell, Physical Review A
**39,**
126-148 (1989).

[4] M. D. Feit and J. A. Fleck, Jr., Applied Optics
**17,**
3990-3998 (1978).

[5] M. D. Feit and J. A. Fleck, Jr., Journal of the Optical
Society of America B
**5,**
633-640 (1988).

[6] J. E. Midwinter, IEE Proceedings Part J
**132,**
371-383 (1985).

[7] L. Thyl\*'en, Paper SB1-1, Numerical Simulation and Analysis in Guided-Wave Optics and Optoelectronics (Optical Society of America, 1989 Technical Digest Series, Volume 3), pp. 20-23 (1989).

[8] Y. R. Shen,
*The Principles of Nonlinear Optics *
(New York: John Wiley and Sons, 1984).

[9] R. Y. Chiao, E. Garmire, and C. H. Townes, Physical Review
Letters
**13,**
479-482 (1964) [Erratum,
**14,**
1056 (1965)].

[10] G. A. Askar'yan, Soviet Physics JETP
**15, **
1088-1090, 1161-1162 (1962).

[11] M. Hercher, Journal of the Optical Society of America
**15, **
563 (1964).

[12] C. Joshi, C. E. Clayton, K. Marsh, Y. Sakawa, and R. L. Savager, Jr.,
Optics Communications
**70,**
44-49 (1989).

[13] S. E. Coe, T. Afshar-rad, and O. Willi, Optics Communications
**73,**
299-303 (1989).

[14] H. J. Zhang, J. H. Dai, P. Y. Wang and L. A. Wu,
Optics Letters
**14,**
695-696 (1989).

[15] E. L. Dawes and J. H. Marburger,
Physical Review
**179,**
862-868 (1969).

[16] V. N. Lugovoi and A. M. Prokhorov, JETP Letters
**7,**
117-119 (1968).

[17] G. M. Zverev, E. K. Maldutis, and V. A. Pashkov,
JETP Letters
**9,**
61-63 (1969).

[18] M. M. T. Loy and Y. R. Shen, Physical Review Letters
**22,**
994-997 (1969).

[19] M. M. T. Loy and Y. R. Shen, Applied Physics Letters
**19,**
285-287 (1971).

[20] C. R. Giuliano and J. H. Marburger, Physical Review Letters
**27,**
905-908 (1971).

[21] M. M. T. Loy and Y. R. Shen, IEEE Journal of Quantum
Electronics
**QE-9,**
409-422 (1973).

**Tags:**
physics
nonlinear optics
QTML equations

1 COMMENTS

* Stephen J. Douglass on Apr 14, 2023 8:16 PM*

I understand very little of this but I'm loving the QTML. You'll have to dumb down the physics a good bit and explain it to me this weekend.

You must be logged in to comment