WAVE OPTICS AND GAUSSIAN BEAMS

This is a continuation from the previous tutorial - RAY OPTICS AND RAY MATRICES

 

A more accurate treatment of optical beams and laser resonators must take into account diffraction and the wave nature of light. Practical laser beams are almost always well enough collimated even under worst conditions, however, that we can describe their diffraction properties using a scalar wave theory, and working in the paraxial wave approximation.

In this tutorial, therefore, we introduce the paraxial wave analysis and the equivalent Huygens-Fresnel integral approach for optical beams in free space. We also introduce the lowest and higher-order gaussian mode solutions of these equations as a widely useful set of "normal modes of free space."

The Hermite-gaussian or Laguerre-gaussian modes which we introduce in this chapter are exact and yet mathematically convenient solutions to the paraxial wave equation in free space. They also provide very close (though not quite exact) approximations for the transverse eigenmodes of stable laser resonators with finite diameter mirrors.

Gaussian beams are therefore very widely used in analyzing laser beams and related optical systems. Our approach in this chapter is to focus primarily on the mathematical derivation of these modes, whereas in the following chapter we summarize most of the important practical properties of gaussian beams in considerable detail.

 

1. THE PARAXIAL WAVE EQUATION

One fundamental way of analyzing free-space wave propagation, using a differential approach, is through the paraxial wave equation, which we can derive once again here in the following fashion.

 

Derivation of the Paraxial Wave Equation

Electromagnetic fields in free space (or in any uniform and isotropic medium) are governed in general by the scalar wave equation

\[\tag{1}[\nabla^2+k^2]\quad\tilde{E}(x,y,z)=0,\]

 

 

 

where \(\tilde{E}(x, y, z)\) is the phasor amplitude of a field distribution that is sinusoidal in time. We will be concerned in this section with optical beams propagating primarily along the \(z\) direction, so that the primary spatial dependence of \(\tilde{E}(x, y, z)\) will be an exp\((-jkz)\) variation. This exp\((-jkz)\) variation has a spatial period of one wavelength \(\lambda\) in the \(z\) direction. 

In addition, for any beam of practical interest the amplitude and phase of the beam will generally have some transverse variation in \(x\) and \(y\) which specifies the beam's transverse profile, as shown in Figure 1; and this transverse amplitude and phase profile will change slowly with distance z due to diffraction and propagation effects.

Both the transverse variations across any plane \(z\), however, and especially the variation in beam profile with distance along the \(z\) axis, will usually be slow compared to the plane-wave exp\((-jkz)\) variation in the \(z\) direction for a reasonably well-collimated beam.

It is then convenient to extract the primary exp\((-jkz)\) propagation factor out of \(\tilde{E}(x, y, z)\), by writing each relevant vector component of the field (such as \(E_x\) or \(E_y)\) in the form

\[\tag{2}\tilde{E}(x,y,z)\equiv\tilde{u}(x,y,z)e^{-jkz},\]

where \(u\) is a complex scalar wave amplitude which describes the transverse profile of the beam. Substituting this into the wave equation 1 then yields, in rectangular coordinates, the reduced equation

\[\tag{3}\frac{\partial^2\tilde{u}}{\partial x^2}+\frac{\partial^2\tilde{u}}{\partial y^2}+\frac{\partial^2\tilde{u}}{\partial z^2}-2jk\frac{\partial\tilde{u}}{\partial z}=0\]

Now, we emphasize once again that with the exp\((-jkz)\) dependence factored out, the remaining z dependence of the wave amplitude \(\tilde{u}(x, y, z)\), is caused basically

by diffraction effects, and this \(z\) dependence will in general be slow compared not only to one optical wavelength, as in exp\((-jkz)\), but also to the transverse variations due to the finite width of the beam. This slowly varying dependence of \(\tilde{u}(x, y,z)\) on \(z\) can be expressed mathematically by the paraxial approximation

\[\tag{4}\bigg|\frac{\partial^2\tilde{u}}{\partial z^2}\bigg|\ll\bigg|2k\frac{\partial\tilde{u}}{\partial z}\bigg|\quad\text{or}\quad\bigg|\frac{\partial^2\tilde{u}}{\partial x^2}\bigg|\quad\text{or}\quad\bigg|\frac{\partial^2\tilde{u}}{\partial y^2}\bigg|.\]

By dropping the second partial derivative in \(z\), we thus reduce the exact wave equation 3. to the paraxial wave equation

\[\tag{5}\frac{\partial^2\tilde{u}}{\partial x^2}+\frac{\partial^2\tilde{u}}{\partial y^2}-2jk\frac{\partial\tilde{u}}{\partial z}=0.\]

More generally we may write this paraxial wave equation as

\[\tag{6}\nabla^2_t\tilde{u}(s,z)-2jk\frac{\partial\tilde{u}(s,z)}{\partial z}=0,\]

where \(s\) refers to the transverse coordinates \(s\equiv(x,y)\;\text{or}\;s\equiv(r,\theta)\), depending on what coordinate system (rectangular or cylindrical) we elect to use, and \(\nabla^2_t\) means the laplacian operator operating on these coordinates in the transverse plane.

This equation will be the primary governing equation for all the analysis of this and the future tutorials.

 

Paraxial Wave Propagation: Finite Difference Approach

The paraxial wave equation can also be turned around and written in the form

\[\tag{7}\frac{\partial\tilde{u}(s,z)}{\partial z}=-\frac{j}{2k}\nabla^2_t\tilde{u}(s,z).\]

This equation can then be integrated forward in the \(z\) direction in order to compute the forward propagation and diffraction spreading of an arbitrary paraxial optical beam.

That is, we can employ any suitable numerical differentiation and integration algorithms, first to evaluate the transverse derivative \(\nabla^2_t\tilde{u}\) at a given plane \(z\), and then to step forward to a new plane \(z+\Delta z\).

We can thus accomplish numerical forward propagation of an arbitrary optical wavefront, making sure to use adequate numbers of sampling points in both the transverse and longitudinal directions.

This numerical approach, sometimes referred to as the "finite difference approach," has been applied to practical beam propagation problems by several workers. For almost any free-space beam propagation problem that we may consider, however, the integral formulation that we will consider in the next section is probably a better choice for numerical calculations, because of the much greater computational efficiency of fast Fourier transforms that can be employed.

 

Validity of the Paraxial Approximation

The paraxial wave equation in either of the above forms is fully adequate for describing nearly all optical resonator and beam propagation problems that arise with real lasers. As perhaps the simplest but most effective way to confirm

 

 

this, and also to illustrate the physical limitations of the paraxial approximation, we can consider the following argument. 

Any optical beam can always be viewed as being made up of a superposition of plane wave components traveling at various angles to the \(z\) axis. (We will discuss this kind of plane wave expansion more rigorously in a later section.) Consider, for example, a plane wave component \(\tilde{E}(x, z)\) traveling at angle \(\theta\) to the \(z\) axis in the \(x,z\) plane, as shown in Figure 2. (For simplicity we consider only one transverse coordinate.) The axial and transverse variations of this plane wave component are then given by

\[\tag{8}\tilde{E}(x,z)=\text{exp}[-jkx\sin\theta-jkz\cos\theta]=\tilde{u}(x,z)e^{-jkz}.\]

The exact form for the reduced wave amplitude \(\tilde{u}(x,y, z)\), and its approximate form within the paraxial approximation, then become 

\[\tag{9}\tilde{u}(x,z)=\text{exp}[-jkx\sin\theta-jkz(1-\cos\theta)]\approx\text{exp}\left[-jk\theta x+jk\frac{\theta^2z}{2}\right].\]

The normalized first and second derivatives of \(\tilde{u}(x,z)\) in the transverse direction then take on the values 

\[\tag{10}\begin{array}&-j\frac{2k}{\tilde{u}}\frac{\partial\tilde{u}}{\partial z}=+2k^2(1-\cos\theta)\approx k^2\theta^2\\\frac{1}{\tilde{u}}\frac{\partial^2\tilde{u}}{\partial x^2}=-k^2\sin^2\theta\approx-k^2\theta^2\end{array}\]

However, the second derivative in the \(z\) direction takes on the form 

\[\tag{11}\frac{1}{\tilde{u}}\frac{\partial^2\tilde{u}}{\partial z^2}=-k^2(1-\cos\theta)^2\approx -\frac{k^2\theta^4}{4}.\]

This particular derivative is smaller than either of the preceding terms by the ratio \(\theta^2/4\) (with \(\theta\) measured in radians)—a ratio which will be \(\ll1\) so long as \(\theta\) is \(\leq1/2\) radian.

We can conclude that so long as all (or at least most) of the plane wave components making up any optical beam are traveling at angles \(\theta\leq 0.5\) rad, or \(\theta\) less than about \(30^\circ\) the \(\partial^2\tilde{u}/\theta z^2\) terms will be at least an order of magnitude smaller than either of the other two terms, in agreement with the basic paraxial approximation.

Paraxial optical beams can thus be focused or can diverge at cone angles up to \(\approx 30^\circ\) before significant corrections to the paraxial wave approximation become necessary.

 

 

2. HUYGENS' INTEGRAL

Another equally valid and effective way of analyzing paraxial wave propagation, but now using an integral approach, is to employ Huygens, principle, expressed in the Fresnel approximation. We can derive this alternative approach to paraxial beam propagation as follows. 

 

Spherical Waves, and the Fresnel Approximation

Let us first note that one very general solution to the exact wave equation, which corresponds physically to a uniform spherical wave diverging from a source point \(r_0\) (Figure 3.) may be written in the form

\[\tag{12}\tilde{E}(\boldsymbol r;\boldsymbol r_0)=\frac{\text{exp}[-jk\rho(\boldsymbol r,\boldsymbol r_0)]}{\rho(\boldsymbol r,\boldsymbol r_0)},\]

where \(\tilde{E}(\boldsymbol r;\boldsymbol r_0)\) means the field at point \(\boldsymbol r\) due to a source at point \(\boldsymbol r_0\), and where the distance \(\rho(\boldsymbol r,\boldsymbol r_0)\) from the source point \(s_0\), \(z_0\) to the observation point \(s,z\) is given by

\[\tag{13}\rho(\boldsymbol r,\boldsymbol r_0\equiv\sqrt{(x-x_0)^2+(y-y_0)^2+(z-z_0)^2}.\]

We emphasize that a spherical wavefunction in this form is an exact solution to the full scalar wave equation of the previous section.

Consider, however, a situation in which the source point \(x_0,y_0,z_0\) for this wave is located somewhere not too far off the \(z\) axis; and suppose we only wish

 

 

 

 

 

to write the resulting field distribution \(\tilde{E}(x,y,z)\;\text{or}\;\tilde{u}(x,y,z)\) of this spherical wave on some transverse plane \(x,y\) that is farther along the \(z\) axis, for values of \(x\) and y that are also not too far off the axis, as in Figure 4. The Fresnel approximation to diffraction theory says that if we expand the distance \(\rho(\boldsymbol r,\boldsymbol r_0)\) in Equation 13. in a power series in the form 

\[\tag{14}\rho(\boldsymbol r,\boldsymbol r_0)=z-z_0+\frac{(x-x_0)^2+(y-y_0)^2}{2(z-z_0)}+\cdots,\]

then we can drop all terms higher than quadratic in this expression, at least in writing the phase shift factor \(\text{exp}[-jk\rho(\boldsymbol r,\boldsymbol r_o)]\). (We will examine the validity of this assumption in more detail a bit further on.) In the 1/p denominator of Equation 12., on the other hand, we will drop even the quadratic terms, and replace \(\rho(\boldsymbol r,\boldsymbol r_0)\) by simply \(z-z_0\).

The spherical wave of Equation 12. is then converted, in this Fresnel approximation, into what we might call a "paraxial-spherical wave" given by

\[\tag{15}\tilde{E}(x,y,z)\approx\frac{1}{z-z_0}\text{exp}\left[-jk(z-z_0)-jk\frac{(x-x_0)^2+(y-y_0)^2}{2(z-z_0)}\right],\]

 

 

 or 

\[\tag{16}\tilde{u}(x,y,z)=\frac{1}{z-z_0}\text{exp}\left[-jk\frac{(x-x_0)^2+(y-y_0)^2}{2(z-z_0)}\right].\]

These expressions approximate the spherical wave by a quadratic phase variation as observed on a transverse plane \(x,y\) located at a distance \(z-z_0\) away from the source point along the \(z\) axis, as shown in Figure 4. The "paraxial-spherical wave" given by this Fresnel approximation is, as the reader can verify, an exact; analytical solution to the paraxial wave equation rather than to the exact wave equation.

 

Huygens' Integral 

We can next connect these spherical-wave and paraxial-spherical-wave ideas to Huygens' integral. 

Huygens' integral originated as an intuitive physical principle, which was only later put into more rigorous mathematical terms. This principle says in physical terms that if we are given an incident field distribution \(\tilde{E}(x_0,y_0,z_0)\) over some closed surface \(S_0\), we may regard the field at each point on that surface as the source for a uniform spherical wave or "Huygens' wavelet" which radiates from that point on the surface, as illustrated in Figure 5. The total field at any other point \(s,z\) inside, or beyond, the surface \(S_0\) can then be calculated by summing the fields of all these Huygens' wavelets coming from all the points on the surface \(S_0\).

Huygens' intuitive ideas concerning this principle were put into more formal mathematical form, first by Fresnel and Kirchoff, and later by Rayleigh and Sommerfeld. The general idea is that each of the Huygens' wavelets should be viewed as a spherical wave with a form like Equation 12., leading to Huygens' integral equation in the form

\[\tag{17}\tilde{E}(s,z)=\frac{j}{\lambda}\int\int_{s_0}\tilde{E_0}(s_0,z_0)\frac{\text{exp}[-jk\rho(\boldsymbol r,\boldsymbol r_0)]}{\rho(\boldsymbol r,\boldsymbol r_0)}\cos\theta(\boldsymbol r,\boldsymbol r_0)d\boldsymbol S_0,\]

where \(\rho(\boldsymbol r,\boldsymbol r_0)\) is the distance between source and observation points as defined earlier.

In this formulation \(d\boldsymbol S_0\) is an incremental element of surface area at point \(s_0, z_0\) on the surface \(S_0\), and the factor \(\cos\theta(\boldsymbol r, \boldsymbol r_0)\) is an "obliquity factor" which depends on the angle \(\theta(\boldsymbol r,\boldsymbol r_0)\) between the line element \(\rho(\boldsymbol r,\boldsymbol r_0)\) and the normal to the surface element \(d\bf S_0\).

 

 

 

Slightly different forms for the obliquity factor in Equation 17. are predicted by the Kirchoff-Fresnel or the Rayleigh-Sommerfeld approaches to diffraction theory, although in either situation this factor goes to unity if the angle \(\theta\) is limited to small values.

The \(j/\lambda\) factor in front of Huygens' integral is a normalization factor which comes out of a more detailed approach to the theory, and which is necessary in order to get the correct near-field and far-field dependence from Huygens' integral. 

 

Fresnel Approximation to Huygens' Integral 

Suppose now that we are given the input field distribution or beam profile \(\tilde{E_0}(x_0,y_0,z_0)\) of a paraxial optical beam across an input transverse plane at \(z=z_0\), as in Figure 6., and we wish to calculate the output beam profile \(\tilde{E_0}(x,y,z)\) across another plane at distance \(z=z_0+L\), accurate to within the paraxial degree of approximation.

We can then make this calculation by using Huygens' principle in the form just given, except that we can replace the exact spherical wave \(\text{exp}[-jk\rho]/\rho\) from Equation 12. with a "paraxial-spherical wave" in the form derived in Equation 15. 

If we do this, we obtain Huygens' integral in the Fresnel approximation, as given by

\[\tag{18}\tilde{E}(x,y,z)\approx\frac{je^{jk(z-z_0)}}{(z-z_0)\lambda}\int\int\tilde{E_0}(x_0,y_0,z_0)\;\text{exp}\left[-jk\frac{(x-x_0)^2+(y-y_0)^2}{2(z-z_0)}\right]dx_0dy_0.\]

We have assumed here that the distance \(z-z_0\) between input and output planes in Figure 6. is large enough so that the angle of the line element connecting source and observation points is always \(\leq 0.5\) rad.

Hence we can approximate the obliquity factor \(\cos\theta\) in Equation 17 by unity, as well as using the Fresnel approximation in the exponent.

This formulation of Huygens' integral can equally well be written for the reduced wavefunction \(\tilde{u}(x,y,z)\) in the form

\[\tag{19}\tilde{u}(x,y,z)=\frac{j}{L\lambda}\int\int\tilde{u}_0(x_0,y_0,z_0)\;\text{exp}\left[-jk\frac{(x-x_0)^2+(y-y_0)}{2L}\right]dx_0dy_0,\]

where \(L\equiv z-z_0\) is again the distance from input to output plane, and the integrations are over the transverse input plane located at \(z=z_0\). We will most often use Huygens' integral in this second form, with the on-axis or plane-wave phase shift exp\((-jkL)\) or exp\([-jk(z-z_0)]\) omitted, since it is the transverse variation of \(\tilde{u}(x,y)\) that is usually of primary interest. The plane-wave phase factor must be brought back into the calculation if resonant frequencies or total phase shifts are to be calculated.

 

Validity of the Fresnel Approximation

The reader can verify that the "paraxial spherical wave" given in Equations 15. or 16. satisfies exactly the paraxial wave equation, just as the exact spherical wave in Equation 12. satisfies the exact wave equation.

Any field expression for \(\tilde{u}(x,y,z)\) given by Huygens' integral in the Fresnel approximation, as given in Equation 19 with \(L=z-z_0\), will therefore satisfy the paraxial wave equation exactly, and vice versa. Huygens' integral and the paraxial wave equation represent exactly the same mathematical (and physical) approximations. 

There are, however, some subtleties in the physical interpretation of the Fresnel approximation which it is useful to understand. Suppose we wish to calculate the forward propagation of a paraxial beam through a distance L, from an input wave \(\tilde{u}(x_1,y_1)\) at plane \(z=z_1\) to an output wave \(\tilde{u}_2(x_2,y_2)\) at plane \(z=z_2=z_1+L\). Suppose also that this beam is confined to a width \(\approx 2a\), i.e., \(\tilde{u}_1\) and \(\tilde{u}_2\) have negligible values outside of \(-a\leq x, y\leq a\).

Now the Fresnel approximation of Equations 14 to 19 assumes that the quartic and higher-order terms in the expansion of the exponent \(e^{-jk\rho}\) can be dropped, because their contribution to the complex exponent will be small compared to, say, \(\pi/2\). If we consider either of the transverse coordinates \(x\) or \(y\), this condition on the next higher-order terms seems to require that

\[\tag{20}\bigg|\frac{k(x_2-x_1)^4}{4L^3}\bigg|\approx\bigg|\frac{2\pi}{\lambda}\frac{(2a)^4}{4L^3}\bigg|\leq\frac{\pi}{2},\]

since by assumption \(|x_2-x_1|\leq 2a\). This condition can be rewritten in the alternative form

\[\tag{21}\frac{L}{2a}\geq\left(\frac{2a}{\lambda}\right)^{1/3}.\]

In any normal situation, where the beam width is large compared to a wavelength, or \(2a\gg\lambda\), this condition apparently says that Huygens' integral in the Fresnel approximation can only be applied to calculate forward beam propagation over lengths that are significantly greater than the beam diameter, or \(L\gg2a\).

There seems to be an inconsistency in this limitation, however, in that the paraxial wave equation 7 of the previous section can obviously be applied over arbitrarily short forward steps in \(z\), and the paraxial wave equation and the Huygens-Fresnel integral are supposed to be mathematically equivalent.

Is the condition in Equation 21 really a limitation on the use of Huygens' integral in the Fresnel approximation?

This question can be answered (in the negative) as follows. Suppose the optical wavefront \(\tilde{u}(x_1,y_1,z_1)\) at the input plane \(z_1\) is a freely propagating paraxial beam, without any sharp discontinuities or aperturing at or near that particular plane.

Then, the effective sources for this beam are really located at source points \(x_0,y_0,z_0\) that are located far behind (or, for a converging beam, far ahead) of either of the planes \(z_1\) or \(z_2\).

We can then clearly apply the Huygens-Fresnel integral to calculate the wave propagation from plane \(z_0\) to plane \(z_1\), and we can equally well apply this integral to calculate the propagation from plane \(z_0\) to plane \(z_2\).

But the inherent cascading properties of the Huygens-Fresnel integral (which the reader may want to verify; see the Problems) then imply that the Huygens-Fresnel integral in exactly the same form must also be valid over the distance \(z_2-z_1=L\), even if this distance is too small to satisfy Equation 21, and in fact even if \(L \rightarrow 0\).

The essential physical point here is that the paraxial or Fresnel approximation is a physical property of the optical beam, not a mathematical property of the Huygens-Frdsnel formulation.

The paraxial wave equation 7 and also the Huygens-Fresnel integral 19 can be applied over arbitrarily short distances \(L\) if the optical beam itself is truly paraxial.

 

Sharp-Edged Aperture Effects

We can give a more physical interpretation to this assertion about paraxial beams—which may seem somewhat confusing at first—by the following illustration. 

Suppose for example that an input wave \(\tilde{u}(x_1,y_1)\) does have some sharp discontinuity in its wavefunction at the input plane z\, either in amplitude or phase, such as would be caused by a hard-edged aperture or by a discontinuous phase step in the input plane \(z_1\), say, at radius \(x=a\).

We will see later that such discontinuities appear to act in effect as sources for quasi spherical diffraction waves or "edge waves" which appear to radiate from the discontinuous edges of the aperture. The criterion of Equation 21. must then be applied, not only to the kernel of Huygens' integral, but also to these "edge waves" as seen at any plane a distance \(L\) away from the aperture plane. Huygens' integral in the Fresnel approximation can thus not be used for distances closer than \(L/2a\approx(2a/\lambda)^{1/3}\) beyond such an aperture. 

But in fact, neither Huygens' integral in the Fresnel approximation nor the paraxial wave equation can be applied to this sort of sharp-edged diffraction situation over distances L shorter than given by the preceding condition.

The Huygens-Fresnel integral cannot be used because this would violate the Fresnel approximation for those diffracted wavelets which appear to be scattered from the edges of the aperture.

The paraxial wave equation cannot be applied accurately in this region, because the rapid changes in \(\tilde{u}(x,z)\) with both x and \(z\) near the sharp aperture edges violate the approximations inherent in the paraxial wave equation.

The fundamental point, then, is that paraxial methods, however they may be formulated, can only be applied to paraxial beams, and a beam diffracted by a sharp-edged aperture does not again become paraxial until we move far enough past the aperture to satisfy Equation 21.

 

Huygens' Integral in One Dimension

Huygens' integral may be rewritten in the general form

\[\tag{22}\tilde{u}(s,z)=\int\int \tilde{K}(\boldsymbol r,\boldsymbol r_0)\tilde{u}_0(\boldsymbol s_0,z_0)d\boldsymbol s_0,\]

where \(\boldsymbol r_0\) and \(\boldsymbol r\) indicate points on the input and output transverse planes; \(\tilde{K}(\boldsymbol r,\boldsymbol r_0)\) is shorthand for the Huygens kernel of Equation 19.; and \(d\boldsymbol s_0\) indicates an element of area on the input plane. In rectangular coordinates the Huygens-Fresnel kernel separates into a product kernel in the form

\[\tag{23}\tilde{K}(\boldsymbol r,\boldsymbol r_0)=\tilde{K}(x-x_0)\times\tilde{K}_1(y-y_0),\]

where the one-dimensional kernel \(\tilde{K}_1\) has the form

\[\tag{24}\tilde{K}_1(x-x_0)=\sqrt{\frac{j}{L\lambda}}\text{exp}\left[-j\frac{\pi(x-x_0)^2}{L\lambda}\right].\]

If the wavefunction \(\tilde{u}(\boldsymbol s, z)\) is also separable in \(x,y\) coordinates, the entire integral can be separated into two one-dimensional integrals of the form

\[\tag{25}\tilde{u}(x,z)=\sqrt{\frac{j}{L\lambda}}\int\tilde{u}_0(x_0,z_0)\;\text{exp}\left[-j\frac{\pi(x-x_0)^2}{L\lambda}\right]dx_0\]

and the same for \(\tilde{u}(y,z)\) and \(\tilde{u}_0(y_0,z_0)\). We will frequently write such integrals in only one transverse dimension, in order to simplify the mathematical expressions

Note that the \(j/L\lambda\) factor in front of the three-dimensional Huygens' integral 19. reduces to \(\sqrt{j/L\lambda}\) if there is only one transverse dimension. If only one transverse coordinate is included, the Huygens' wavelet is in essence a cylindrical wave rather than a spherical wave.

Hence its amplitude decreases with distance as \(1/\sqrt\rho\) or \(1/\sqrt L\), rather than as \(1/\rho\).

The phase shift of \(\pi/2\) due to the \(j\) factor in two dimensions also reduces to \(\pi/4\) in one dimension.

This phase factor represents in essence an initial phase shift of the Huygens' wavelet compared to the actual field value at the input point. We will see later that this corresponds to a well-known "Guoy phase shift" associated with any wave passing through a focus, or coming from a small enough source point.

 

 

 

 

3. GAUSSIAN SPHERICAL WAVES 

Our next important step is to derive the analytical form for a gaussian spherical wave, or a so-called "gaussian beam" in free space; and then to show that this gaussian beam is a very useful exact solution to the paraxial wave equation, or to Huygens' integral in the Fresnel approximation. 

 

Paraxial Spherical Waves 

Consider a uniform spherical wave diverging from a source point located at \(x_0,y_0,z_0\) and observed at an observation point \(x,y,z\), as in Figure 7. If the axial distance \(z-z_0\) between the source and observation points is sufficiently large compared to the transverse coordinates \(x_0\), \(y_0\) and \(x,y\), then the field distribution produced by this wave at point \(x,y\) on the plane located at distance \(z\) can be written, using the paraxial approximation, in the form

\[\tag{26}\tilde{u}(x,y,z)=\frac{1}{z-z_0}\;\text{exp}\left[-jk\frac{(x-x_0)^2+(y-y_0)^2}{2(z-z_0)}\right]=\frac{1}{R(z)}\;\text{exp}\left[-jk\frac{(x-x_0)^2+(y-y_0)^2}{2R(z)}\right],\]

where \(R(z)=z-z_0\) gives the radius of curvature of the spherical wave at plane \(z\). The phase variation exp\([-j\phi(x,y,z)]\) across a transverse plane at fixed \(z\) for such a paraxial spherical wave with radius of curvature \(R(z)\) thus has the quadratic form

\[\tag{27}\phi(x,y,z)\equiv k\frac{(x-x_0)^2+(y-y_0)^2}{2(z-z_0)}=\frac{\pi}{\lambda}\frac{(x-x_0)^2+(y-y_0)^2}{R(z)}.\]

The radius of curvature \(R(z)\) of the wave at plane \(z\) can be written in a more general form as

\[\tag{28}R(z)=R_0+z-z_0,\]

with \(R_0\) being the value at the earlier plane \(z_0\) (and with \(R_0=0\) if the earlier plane is the source plane, as in the present situation). As such a spherical wave propagates forward to any other plane \(z\), the radius of curvature of the wave thus increases linearly with distance.

Note that with our sign convention, a value of \(R>0\) indicates a diverging or expanding wave, whereas \(R<0\) indicates a converging wave moving inward toward the source point.

Note also that if the wave is propagating in some dielectric medium, \(k\) and \(\lambda\) are the values in that medium, not in free space. This quadratic phase variation of course represents only a paraxial or Fresnel approximation to the true surface of a sphere, so that this form will have a sizable phase error if we move far enough out from the optic axis.

 

Introducing Complex Source Point Coordinates

A paraxial spherical wave in the form given in Equation 26. cannot by itself be a very useful analytical form for a real physical beam, however, because the amplitude of the spherical wave does not fall off with transverse distance from the axis. Such a wave instead extends out to infinity in the transverse direction, and carries infinite energy and power across the transverse plane (as well as having large deviation from a true sphere far off the axis). 

A very simple way to overcome these difficulties can be developed, however, as follows. Let us first note that the spherical wave expressions in Equations 26.-28. satisfy the paraxial wave equation, or the Huygens-Fresnel integral, exactly for any arbitrary choice of the source point coordinates \(x_0,y_0,z_0\).

That is, these coordinates are simply constant parameters, which cancel out identically when the spherical wave expression is put into the paraxial wave equation or the Huygens-Fresnel integral. What then will happen if we explore the possibility of employing complex values for these source point coordinates?

In particular, suppose that for simplicity we set \(x_0\) and \(y_0\) to zero, but that we convert the axial location \(z_0\) of the source point into a complex number, by subtracting from it an arbitrary complex quantity which we will call \(\tilde q_0\). That is, we replace the purely real value \(z_0\) in the spherical wave expression by the

 

 

 

complex value \(z_0-\tilde{q}_0\). This amounts to the same thing as replacing the radius of curvature \(R(z)=R_0+z-z_0\) by the complex quantity \(\tilde{q}(z)=z-(z_0-\tilde{q}_0)=\tilde{q}_0+z-z_0\). [We introduce a new notation \(\tilde{q}(z)\) in place of \(R(z)\) because this quantity will turn out to be a complex generalization of the purely real sphericalwave radius of curvature \(R(z)]\).

Note that the value of this new "complex radius" at the source plane \(z=z_0\) is just \(\tilde{q}(z_0)=\tilde{q}_0\).

The spherical wave diverging from this complex source point is now written, following the same form as Equation 26., as

\[\tag{29}\begin{array}\tilde{u}(x,y,z)&=\frac{1}{z-z_0+\tilde{q}_0}\text{exp}\left[-jk\frac{x^2+y^2}{2(z-z_0+\tilde{q}_0)}\right]\\&=\frac{1}{\tilde{q}(z)}\text{exp}\left[-jk\frac{x^2+y^2}{2\tilde{q}(z)}\right]\end{array}\],

where the complex radius \(\tilde{q}(z)\) is given by 

\[\tag{30}\tilde{q}(z)=\tilde{q}_0+z-z_0\]

in direct analogy to the expression 28. for \(R(z)\).

But if \(\tilde{q}(z)\) is complex, then we can separate the exponent in Equation 29. into real and imaginary parts, by first separating the quantity \(1/\tilde{q}(z)\) into real and imaginary parts in the form 

\[\tag{31}\frac{1}{\tilde{q}(z)}\equiv\frac{1}{q_r(z)}-j\frac{1}{q_i(z)},\]

and then writing the spherical wave expression in the form 

\[\tag{32}\tilde{u}(x,y,z)=\frac{1}{\tilde{q}(z)}\text{exp}\left[-jk\frac{x^2+y^2}{2q_r(z)}-k\frac{x^2+y^2}{2q_i(z)}\right].\]

(Note that \(q_r\) and \(q_i\) as defined here are not in general the real and imaginary parts of \(\tilde{q}(z)\), but rather are defined by Equation 31.) 

The resulting exponent for this complex-source-point beam now has both an imaginary quadratic transverse variation, corresponding to a quadratic phase front, or a spherical wave with a real radius of curvature; and also a purely real quadratic transverse variation, which gives a gaussian transverse amplitude profile or amplitude variation, with a transverse fall-off determined by the imaginary part of \(1/\tilde{q}\).

Both of these variations are contained in the complex radius of curvature \(\tilde{q}\) through Equation 31.

 

Gaussian-Spherical Beams

At this point let us convert this result into the standard notation which is very widely used in the laser field, by rewriting it in the form 

\[\tag{33}\begin{array}\tilde{u}(x,y,z)&=\frac{1}{\tilde{q}(z)}\text{exp}\left[-jk\frac{x^2+y^2}{2\tilde{q}(z)}\right]\\&\equiv\frac{1}{\tilde{q}(z)}\text{exp}\left[-jk\frac{x^2+y^2}{2R(z)}-\frac{x^2+y^2}{w^2(z)}\right]\end{array}\]

where \(R(z)\) is now the radius of curvature and \(w(z)\) the so-called "gaussian spot size" of this solution. Equation 33. gives the lowest-order spherical-gaussian beam solution in free space.

We can view this solution if we like as being simply a paraxial-spherical wave diverging from a complex source point, which is located a complex distance \(z-z_0+\tilde{q}_0\) rather than a real distance \(z-z_0\) behind the observation plane \(z\). 

It is very important to note here that, although we use the same notation for \(R(z)\) as before, this radius ofcurvature can no longer be calculated by the formula \(R(z)=R(z_0)+z-z_0\) from Equation 28. Rather, the radius of curvature \(R(z)\) and the spot size \(w(z)\) of the wave at any plane \(z\) are derived from the complex radius \(\tilde{q}(z)\) by the definition that

\[\tag{34}\frac{1}{\tilde{q}(z)}\equiv\frac{1}{R(z)}-j\frac{\lambda}{\pi w^2(z)}\]

The variation of the complex radius of curvature \(\tilde{q}(z)\) with distance is then determined by the formula

\[\tag{35}\tilde{q}(z)=\tilde{q}_0+z-z_0.\]

The fundamental propagation law for all such gaussian beams in free space is entirely contained in this simple relation for \(\tilde{q}(z)\).

 

Summary

The primary result of this section, then, is that replacing the purely real radius of curvature \(R(z)\) for a spherical wave coming from a real source point \(z_0\) with a complex radius of curvature \(\tilde{q}(z)\), or a complex source point \((z_0-\tilde{q}_0)\), converts the paraxial-spherical wave solution given in the opening paragraphs of this section into the gaussian-spherical wave solution given by Equations 29.-35.

This gaussian-spherical wave solution is still an exact mathematical solution to either the paraxial wave equation or the Huygens-Fresnel integral. Now, however, it has in addition the physically desirable properties that its amplitude falls off smoothly and rapidly with distances from the \(z\) axis; it carries finite total power across the beam cross section; and it also remains complex gaussian in profile at all later planes \(z\).

This gaussian-spherical solution, together with various higher-order Hermitegaussian or Laguerre-gaussian extensions, will prove to be extraordinarily useful in the analysis of optical resonators and beams.

This gaussian-spherical wave solution and its higher-order extensions for paraxial beams in free space can in fact be derived in at least four different ways, by using:

\((\text i)\) The complex source point derivation used in the present section; or 

\((\text{ii})\) A differential equation approach based on the paraxial wave equation, which we will use in the following two sections; or 

\((\text{iii})\) Direct substitution of the spherical-gaussian solution into the Huygens-Fresnel integral; or

\((\text{iv})\) A plane-wave expansion approach, which we will introduce in several later sections.

The complex-source-point approach we have used here is, despite its simplicity, both subtle and entirely rigorous. We will explore the physical properties of this gaussian mode solution, and its higher-order extensions, in great detail in the following tutorial. 

 

4. HIGHER-ORDER GAUSSIAN MODES

The gaussian beam expression we introduced in the previous section is really only the lowest-order solution in an infinite family of higher-order solutions to the same Huygens' integral, or the same free-space paraxial wave equation. The higher-order solutions to the same equations can take the form either of Hermitegaussian functions in rectangular coordinates, or of Laguerre-gaussian functions in cylindrical coordinates.

These higher-order gaussian modes are of considerable importance both in practical lasers and in optical beam analyses. Hence we reproduce their mathematical derivation in some detail in this section.

 

Lowest-Order Mode: The Differential Approach

A slightly more formal way of deriving the expressions we have just given for a lowest-order gaussian beam is to assume a trial solution to the paraxial wave equation of the form 

\[\tag{36}\tilde{u}(x,y,z)=A(z)\times\text{exp}\left[-jk\frac{x^2+y^2}{2\tilde{q}(z)}\right],\]

where \(A(z)\) and \(\tilde{q}(z)\) are initially assumed to be unknown. If we substitute this trial solution into the paraxial wave equation, we obtain 

\[\tag{37}\left[\left(\frac{k}{2}\right)^2\left(\frac{d\tilde{q}}{dz}-1\right)(x^2+y^2)-\frac{2jk}{\tilde{q}}\left(\frac{\tilde{q}}{A}\frac{dA}{dz}+1\right)\right]A(z)=0.\]

Now, the only way in which this equation can be satisfied for all \(x\) and \(y\) is to set both of the differential expressions inside the large square brackets to zero. We then have the two differential equations 

\[\tag{38}\frac{d\tilde{q}}{dz}=1\quad\text{and}\quad\frac{dA(z)}{dz}=-\frac{A(z)}{\tilde{q}(z)}\]

whose solutions are given by 

\[\tag{39}\tilde{q}(z)=\tilde{q}_0+z-z_0\quad\text{and}\quad\frac{A(z)}{A_0}=\frac{\tilde{q}_0}{\tilde{q}(z)}.\]

These correspond exactly to the complex-source-point expressions derived in the preceding section. 

 

Higher-Order Solutions in Rectangular Coordinates

We will now show how this same trial solution approach can be extended to find higher-order Hermite-gaussian eigensolutions \(\tilde{u}_{nm}(x,y,z)\) or Laguerregaussian eigensolutions \(\tilde{u}_{pm}(r,\theta,z)\) to the paraxial wave equation. 

We will derive first the higher-order Hermite-gaussian solutions to the wave equation in rectangular coordinates. In rectangular coordinates the elementary solutions can be separated into products of identical solutions in the \(x\) and \(y\) directions, i.e.,

\[\tag{40}\tilde{u}_{nm}(x,y,z)=\tilde{u}_n(x,z)\times\tilde{u}_m(y,z),\]

where \(\tilde{u}_n(x,z)\) and \(\tilde{u}_m(y,z)\) have the same mathematical form. We can therefore find the solutions in only one rectangular coordinate and bring in the other coordinate by analogy. 

The paraxial wave equation in one transverse coordinate reduces to

\[\tag{41}\frac{\partial^2\tilde{u}_n(x,z)}{\partial x^2}-2jk\frac{\partial\tilde{u}_n(x,z)}{\partial z}=0.\]

As one way of looking for higher-order solutions, let us now write a more general trial solution for the wave amplitude \(\tilde{u}(x,z)\) in the form

\[\tag{42}\tilde{u}_n(x,z)=A(\tilde{q}(z))\times h_n\left(\frac{x}{\tilde{p}(z)}\right)\times\text{exp}\left[-jk\frac{x^2}{2\tilde{q}(z)}\right],\]

where \(\tilde{q}=\tilde{q}(z)\) is the same as in the preceding; \(A\tilde({q})\) and \(h_n(x/\tilde{p})\) are initially unknown functions; and \(\tilde{p}=\tilde{p}(z)\) is a distance-dependent scaling factor in the argument of \(h_n\).

Substituting this form into the paraxial wave equation, and assuming that \(\tilde{q}(z)\) will continue to obey the propagation rule \(d\tilde{q}/dz=1\), then converts the paraxial wave equation into a differential relation for \(h_n(x/\tilde{p})\), namely

\[\tag{43}H^"_n-2jk\left[\frac{\tilde{p}}{\tilde{q}}-\tilde{p}'\right]xh'_n-\frac{jk\tilde{p}^2}{\tilde{q}}\left[1+\frac{2\tilde{q}}{A}\frac{dA}{d\tilde{q}}\right]h_n=0,\]

where \(h'_n\) and \(h^"_n\) mean the first and second derivatives of hn with respect to its total argument, e.g., \(h_n(y)=dh_n(y)/dy\), and \(\tilde{p}'\equiv d\tilde{p}(z)/dz\).

But, Equation 43. is very similar to the standard differential equation for the Hermite polynomials \(H_n(x/\tilde{p})\), which has the form 

\[\tag{44}H^"_n-2(x/\tilde{p})H'_n+2nH_n=0.\]

The two equations for \(h_n\) and \(H_n\) will in fact become the same if we can find solutions for \(\tilde{p}(z)\) and \(A\tilde{(q)}\) which satisfy simultaneously the two conditions

\[\tag{45}2jk\left[\frac{\tilde{p}}{\tilde{q}}-\frac{d\tilde{p}}{dz}\right]=\frac{2}{\tilde{p}}\quad\text{or}\quad\frac{d\tilde{p}}{dz}=\frac{\tilde{p}}{\tilde{q}}+\frac{j}{k\tilde{p}},\]

and

\[\tag{46}\frac{-jk\tilde{p}^2}{\tilde{q}}\left[1+\frac{2q}{A}\frac{dA}{d\tilde{q}}\right]=2n\quad\text{or}\quad\frac{2q}{A}\frac{dA}{d\tilde{q}}=\frac{2jnk\tilde{p}^2}{\tilde{q}}-1.\]

Now, there are at least two, and probably many different ways in which Equations 45. and 46. can be solved, with each solution leading to a different family of higher-order Hermite-gaussian solutions. We will describe in this section one family of such solutions, which we will refer to as the "standard" Hermite-gaussian solutions. In the following section we will describe another alternative set of solutions which we will refer to as the "elegant" solutions.

 

The "Standard" Hermite Polynomial Solutions

The set of Hermite-gaussian solutions that we will derive in this section are by far the most widely used set of such solutions, as well as being the closest to simple physical solutions in ordinary stable lasers.

However, this set is also perhaps the most complicated and inelegant approach from a mathematical viewpoint. This standard approach to Hermite polynomial solutions is obtained by assuming that the scale factor \(\tilde{p}(z)\) in the function \(h_n(x/\tilde{p})\) will be purely real, and in fact will be related to the gaussian spot size \(w(z)\) in the form 

\[\tag{47}\frac{1}{\tilde{p}(z)}\equiv\frac{\sqrt 2}{w(z)}.\]

As motivation for this approach we can note that if this is valid then the higherorder solutions of Equation 42., namely 

\[\tag{48}\tilde{u}_n(x,z)=h_n\left(\frac{\sqrt2x}{w(z)}\right)\text{exp}\left[\frac{-jkx^2}{2R(z)}-\frac{x^2}{w^2(z)}\right],\]

will have the same normalized shape at every transverse plane \(z\). That is, these functions will change in transverse scale like \(w(z)\), and will acquire spherical curvature \(R(z)\), but their amplitude profiles will remain unchanged in shape at any plane \(z\).

We must first verify that this form for \(\tilde{p}(z)\) will satisfy the differential equation 45. Since the spot size \(w(z)\) is related to \(\tilde{q}(z)\) by

\[\tag{49}\frac{1}{\tilde{q}(z)}=\frac{1}{R(z)}-j\frac{\lambda}{\pi w^2(z)}\]

we can use this to obtain

\[\tag{50}\frac{1}{w^2(z)}=\frac{jk}{4}\left[\frac{1}{\tilde{q}(z)}-\frac{1}{\tilde{q}^*(z)}\right]=\frac{jk}{4}\frac{\tilde{q}(z)-\tilde{q}(z)}{\tilde{q}(z)\tilde{q}^*(z)}.\]

We can also note that the formulas for \(\tilde{q}(z)\) imply the useful relations that 

\[\tag{51}d\tilde{q}^*/dz=d\tilde{q}/dz=1\quad\text{and}\;\text{hence}\quad\tilde{q}^*-\tilde{q}=\tilde{q}_0^*-\tilde{q}_0^*.\]

The reader can then verify that the form for \(\tilde{p}(z)\) given in Equation 47. does indeed satisfy Equation 45.

The simplest way to satisfy the equation for \(A\tilde{(q)}\) is then perhaps to use the definition of \(1/\tilde{q}\), and the fact that \(d\tilde{q}^*=d\tilde{q}\) to rewrite Equation 46 as

\[\tag{52}\frac{dA}{A}=\frac{1}{2}\frac{d\tilde{q}}{\tilde{q}}+\frac{n}{2}\left(\frac{d\tilde{q}^*}{\tilde{q}^*}-\frac{d\tilde{q}}{\tilde{q}}\right).\]

Integrating this equation then yields 

\[\tag{53}A(\tilde{q})=A_0\times\left(\frac{\tilde{q}_0}{\tilde{q}(z)}\right)^{1/2}\times\left(\frac{\tilde{q}_0}{\tilde{q}_0}\frac{\tilde{q^*}(z)}{\tilde{q}(z)}\right)^{1/2}.\]

A complete set of properly normalized higher-order Hermite-gaussian mode functions for a beam propagating in free-space are thus given, in one transverse dimension, by

\[\tag{54}\tilde{u}_n(x,z)=\left(\frac{2}{\pi}\right)^{1/4}\left(\frac{1}{2^nn!w_0}\right)^{(1/2)}\left(\frac{\tilde{q}_0}{\tilde{q}(z)}\right)^{1/2}\left[\frac{\tilde{q}_0}{\tilde{q}_0^*}\frac{\tilde{q}^*(z)}{\tilde{q}(z)}\right]^{n/2}\times H_n\left(\frac{\sqrt 2x}{w(z)}\right)\text{exp}\left[-j\frac{kx^2}{2\tilde{q}(z)}\right],\]

where the \(H_n's\) are the Hermite polynomials of order \(n\), and \(\tilde{q}(z)\) and \(w(z)\) are exactly the same as for the lowest-order gaussian mode.

 

The Guoy Phase Shift

The most compact and efficient way of writing the higher-order Hermite-gaussian eigenmodes is as in Equation 54., using the ratios of \(\tilde{q}{z)\) and \(\tilde{q}^*(z)\) raised to appropriate powers.

This form can also be converted, however, to a more commonly used form involving the real spot size \(w(z)\) and a phase angle \(\psi(z)\), as follows.

Let us associate a magnitude and especially a phase angle \(\psi(z)\) with the complex \(\tilde{q}\) parameter at any plane \(z\) by writing

\[\tag{55}\frac{j}{\tilde{q}}=\frac{\lambda}{\pi w^2}\left[1+j\frac{\pi w^2}{R\lambda}\right]\equiv\frac{\text{exp}[j\psi(z)]}{|\tilde{q}|},\]

so that the phase angle \(\psi=\psi(z)\) is given at any plane \(z\) by

\[\tag{56}\tan\psi(z)\equiv\frac{\pi w^2(z)}{R(z)\lambda}.\]

We have included the factor of j in Equation 55. because it will be convenient later on to have \(\psi(z)=0\) at the "waist" of a gaussian beam, where the spot size w is finite but the radius of curvature \(R\) becomes infinite.

If we use this definition, we can then show (after some algebra) that the first part of the \(\tilde{q}_0/w_0\tilde{q}(z)\) normalization factor in Equation 54. can be written as

\[\tag{57}\frac{1}{w_0}\frac{\tilde{q}_0}{\tilde{q}(z)}=\frac{\text{exp}[j\psi(z)-\psi_0]}{w(z)},\]

where \(\psi_0\equiv\psi(z_0)\)  is the initial value of \(\psi(z) at \(z=z_0\).

The lowest-order gaussian-spherical wave (Equation 33.) may then be written in the alternative form 

\[\tag{58}\tilde{u}_0(x,z)=\left(\frac{2}{\pi}\right)^{1/4}\sqrt{\frac{\text{exp}j[\psi(z)-\psi_0]}{w(z)}}\text{exp}\left[-j\frac{kx^2}{2\tilde{q}(z)}\right].\]

In other words, the factor \((1/w_0)\times(\tilde{q}_0/\tilde{q}(z))\) contains the necessary \(1/w(z)\) normalization factor in front of the gaussian beam expression, along with an added phase shift term given by \(\psi(z)-\psi_0\).

For all higher-order Hermite-gaussian modes we must also include the additional factors given by the term

\[\tag{59}\left[\frac{\tilde{q}_0}{\tilde{q}^*_0}\frac{\tilde{q}^*(z)}{\tilde{q}(z)}\right]^{n/2}\equiv\text{exp}[jn[\psi(z)-\psi_0]]\]

appearing on the right-hand side of Equation 54. This factor gives rise to a pure phase shift. With this factor included, the higher-order Hermite-gaussian mode functions of Equation 54. can be written in the alternative form

\[\tag{60}\tilde{u}_n(x)=\left(\frac{2}{\pi}\right)^{1/4}\sqrt{\frac{\text{exp}[-j(2n+1)(\psi(z)-\psi_0)]}{2^nn!w(z)}}\times H_n\left(\frac{\sqrt 2x}{w(z)}\right)\text{exp}\left[-j\frac{kx^2}{2R(z)}-\frac{x^2}{w^2(z)}\right].\]

We will discuss the physical significance of the so-called "Guoy phase shift term" \(\psi(z)\) in the following tutorial.

 

Hermite-Gaussian Mode Expansions

The Hermite-gaussian functions \(\tilde{u}_n(x,z)\) we have derived here provide a complete basis set of orthogonal functions characterized by a single complex parameter, the complex \(\tilde{q}_0\) parameter at any arbitrary reference plane \(z_0\). (We will discuss the physical significance of this parameter in the following tutorial.) These functions obey the orthonormality condition

\[\tag{61}\int^\infty_{-\infty}u^*_n(x,z)\tilde{u}_m(x,z)dx=\delta_{nm},\]

independent of either \(z\) or of \(\tilde{q}_0\). They can thus be used as a basis set to expand any arbitrary paraxial optical beam \(\tilde{E}(x,y,z)\) in the form

\[\tag{62}\tilde{E}(x,y,z)=\sum_n\sum_m c_{nm}\tilde{u}_n(x,z)\tilde{u}_m|(y,z)e^{-jkz}.\]

If we multiply both sides of this by \(u^*_n(x,z)u_m^*(y,z)\) and integrate across the full cross section, we can find that the expansion coefficients \(c_{nm}\) are given by

\[\tag{63}c_{nm}=\int^\infty_{-\infty}\int_{-\infty}^\infty\tilde{E}(x,y,z)u^*_n(x,z)u^*_m(y,z)dx\;dy.\]

The coefficients \(c_{nm}\) will depend upon the arbitrary choice of \(\tilde{q}_0\) at \(z_0\). There is thus in general no unique or necessary way of choosing the waist size \(w_0\) or waist

location \(z_0\) for the basis set to expand a given beam pattern \(\tilde{E}(x,y)\). We may attempt to choose these parameters to given an expansion which best fits some other physical constraint in the problem, or which gives an expansion that best fits the actual field \(\tilde{E}(x,y,z)\) with the smallest number of terms.

 

Astigmatic Mode Functions

As we noted earlier, the Hermite-gaussian function \(\tilde{u}_n(x,z)\) in one dimension corresponds essentially to a cylindrical wave, and hence has a normalization factor of \(1/\tilde{q}^{1/2}(z)\) or \(1/w^{1/2}(z)\) rather than \(1/\tilde{q}(z)\) or \(1/w(z)\).

The overall normalization factor for the field function \(\tilde{u}_{nm}(x,y,z)\) is then the product of the individual normalization functions for \(\tilde{u}_n(x,z)\) and \(\tilde{u}_m(y,z)\).

There is no fundamental reason in fact why the complex beam parameters \(\tilde{q}_{0x}\), \(w_{0x}\) and \(z_{ox}\) associated with the functions \(\tilde{u}_n{x,z)\) in the \(x\) transverse coordinate cannot have entirely different values from the corresponding parameters \(\tilde{q}_{0y}\), \(w_{0y}\) and \(z_{0y}\) associated with the functions \(\tilde{u}_m(y,z)\) in the \(y\) direction.

The gaussian beam solutions can thus be converted, where it seems necessary or useful, into a set of somewhat more general astigmatic gaussian beam modes by replacing \(\tilde{q}(z)\), and the related quantities \(w(z)\) and \(\psi(z)\), by separate values for the \(x\) and \(y\) coordinates wherever these quantities appear in the Hermite-gaussian solutions.

Note in particular that only half of the fundamental Guoy phase comes from each transverse coordinate. The Guoy phase shift for a uniform cylindrical wave through a focus is only \(90^\circ\) rather than \(180^\circ\), and the phase factor in Huygens' integral in one transverse dimension is \(\sqrt j\) rather than \(j\).

 

Cylindrical Coordinates: Laguerre-Gaussian Modes

An alternative but equally valid family of solutions to the paraxial wave equation can be written in cylindrical rather than rectangular coordinates. These solutions are in general the Laguerre-gaussian solutions of the form

\[\tag{64}\begin{array}&\tilde{u}_{pm}(r,\theta,z)=\sqrt{\frac{2p!}{(1+8_{0m})\pi(m+p)!}}\frac{\text{exp}j(2p+m+1)(\psi(z)-\psi_0)}{w(z)}\\\times\left(\frac{\sqrt 2r}{w(z)}\right)^m\;L^m_p\left(\frac{2r^2}{w(z)^2}\right)\;\text{exp}\left[-jk\frac{r^2}{2\tilde{q(z)}}+im\theta\right]\end{array}\] 

In these solutions the integer \(p\geq 0\) is the radial index and the integer \(m\) is the azimuthal mode index; the \(L^l_p\) functions are the generalized Laguerre polynomials; and all the other quantities \(\tilde{q}\), \(w\) and \(\psi\) are exactly the same as in the Hermite-gaussian situation. 

These solutions are written using the "standard" transverse scaling \(r/\tilde{p}(z)=\sqrt 2r/w(z)\) that we used for the Hermite-gaussian solutions in Equation 47. of this section.

An alternative set of complex Laguerre-gaussian cylindrical solutions using the complex scaling \(r/\tilde{p}(z)=\sqrt{jkr^2/2\tilde{q}(z)}\) which we will introduce in the following section could equally well be developed.

In either situation these modes have cylindrical symmetry, with modes having circles of constant intensity in the radial direction and an \(e^{im\theta}\) variation in the azimuthal direction.

Alternatively linear combinations of the \(\pm m\) terms can be formed to give \(\cos\;m\theta\) and/or \(\sin\;m\theta\) variations, leading to \(2m\) nodal lines running radial outward from the mode axis. Laguerre-gaussian exhibit the same Guoy phase shift as the rectangular modes.

The Laguerre-gaussian solutions provide an equally general but alternative basis set to the Hermite-gaussian solutions for expanding an arbitrary optical beam \(\tilde{u}(r,\theta,z)\) in free space (provided we are knowledgeable about generalized Laguerre polynomials).

Since both the Hermite-gaussian and Laguerre-gaussian functions form complete sets, we must be able to expand any Hermite solution in terms of the Laguerre functions and vice versa.

The Laguerre-gaussian functions will perhaps be more convenient for problems having a large amount of cylindrical symmetry, and will probably not provide the most convenient set for expanding any real optical beam having substantial astigmatism between \(x\) and \(y\) axes.

In real lasers the Brewster windows and any other tilted surfaces or distorted elements usually provide a small but inherent rectangular symmetry to the laser cavity. Real lasers, therefore, overwhelmingly elect to oscillate in near-Hermite-gaussian rather than near- Laguerre-gaussian modes. Experiments with very carefully aligned gas lasers having internal mirrors and no Brewster windows have, however, clearly demonstrated oscillations in Laguerre-gaussian modes with higher-order radial and azimuthal symmetry.

 

 

5. COMPLEX-ARGUMENT GAUSSIAN MODES

The "standard" Hermite-gaussian solutions developed in the preceding section are the most commonly used form for the Hermite-gaussian eigenmodes, and the form most often given in the laser literature.

They represent, among other things, the set of Hermite-gaussian modes that are the closest approximation to the actual higher-order modes of finite-mirror stable resonators. These modes are perhaps somewhat inelegant mathematically, however, in that we have a complex scaling factor in the gaussian function but only a real scaling factor in the Hermite polynomials.

As a result, Equation 54. is somewhat messy, since the inelegant combinations of \(\tilde{q}\) and \(\tilde{q}^*\) values must be carried along in all the normalization factors.

There are, however, many other alternative choices for \(\tilde{p}(z)\) and \(A(\tilde{q})\) that will satisfy the differential equations 45. and 46. for these quantities derived in the preceding section.

To illustrate this we will develop one of these alternative families of solutions in this section. The motivation behind this particular solution is to use the same complex scaling factor, that is, the quantity \(\sqrt{jkx^2/2\tilde{q}}\), as the argument both in the gaussian exponent and in the Hermite polynomial functions.

 

The "Elegant" Hermite Polynomial Solutions

To do this we will define the complex scale factor \(\tilde{p}\) in the functions \(h_n(x/p)\) of the previous section not by \(1/\tilde{p}=\sqrt 2/w\) as in Equation 47., but by 

\[\tag{65}\frac{1}{\tilde{p}(z)}\equiv\frac{jk}{2\tilde{q}(z)}.\]

The reader can verify that this choice will also satisfy the differential equation 45., and that the differential condition 46. for A\((\tilde{q})\) then takes on the significantly simpler form

\[\tag{66}\frac{\tilde{q}}{a}\frac{dA}{d\tilde{q}}=-\frac{n+1}{2}.\]

If we solve this, the Hermite-gaussian eigenfunctions then take on the alternative form 

\[\tag{67}\hat{u}_n(x,z)=\hat{u_0}\left[\frac{\tilde{q}_0}{\tilde{q}(z)}\right]^{n+1/2}H_n\left(\sqrt{\frac{jkx^2}{2\tilde{q}(z)}}\right)\text{exp}\left[-j\frac{kx^2}{2\tilde{q}(z)}\right].\]

 

 

 

This alternative solution puts the same complex argument \(\sqrt{jkx^2/2\tilde{q}(z)}\) into both the Hermite polynomial and the gaussian exponent, and is much simpler than the "standard" solutions given in Equation 54. 

This alternative set of Hermite-gaussian solutions represents an equally valid complete set of solutions to the paraxial wave equation in free space. This alternative set is more compact and elegant, with some interesting analytical properties, but it is also perhaps less directly useful physically.

For example, whereas these alternative functions \(\tilde{u}_n(x,z)\) still form a mathematically complete set, they are no longer orthogonal to each other in the usual sense. Rather the alternative functions \(\tilde{u}_n\) are biorthogonal to a set of adjoint functions \(v_n\) given by

\[\tag{68}\hat{v}_n(x,z)=H_n\left(\sqrt{\frac{-jk}{2\tilde{q}^*}}x\right),\]

(no gaussian factor) with the orthogonality relation now being

\[\tag{69}\int^\infty_{-\infty}\hat{u}_n^*(x,z)dz=c_n\delta_{nm},\] 

where \(c_n\) is an appropriate normalization constant.

 

Properties of the "Elegant"Solutions

The lowest-order or \(n=0\) and \(n=1\) members of the "standard" and the alternative or "elegant" sets of Hermite-gaussian functions given in Equations 54. and 67. are indistinguishable from each other, since they consist only of the gaussian exponential, or of this exponential multiplied by \(x\).

There are, however, significant differences between the higher-order modes in the two sets. 

The next higher-order function \(n=2\), for example, uses the Hermite polynomial \(H^2(x)=4x^2-2\), so that the "standard" solution \(\tilde{u}_2(x)\) has the form 

\[\tag{70}\tilde{u}^2(x,z)=\text{const}\times\left[\frac{4x^2}{w^2}-1\right]\text{exp}\left[-j\frac{kx^2}{2\tilde{q}}\right],\]

 

 

 

whereas the "elegant" solution \(\hat{u}_2{x)\) becomes 

\[\tag{71}\begin{array}\&\hat{u}_2(x)=\text{const}\times\left[\frac{jkx^2}{\tilde{q}}-1\right]\text{exp}\left[-j\frac{kx^2}{2\tilde{q}}\right]\\=\text{const}\times\left[\frac{2(1+ja)x^2}{w^2}-1\right]\text{exp}\left[-j\frac{kx^2}{2\tilde{q}}\right]\end{array}\]

where \(a=\pi w^2/R\lambda\). The two functions are different, though somewhat similar, even for \(a=0\); and for nonzero values of a the complex argument in the polynomial causes the usual nulls in the function to be filled in, as illustrated in Figure 10., and also produces an additional phase variation across the beam which is not purely spherical in form, as Figure 10. also shows.

For \(n\geq 2\) the amplitude and phase patterns of the higher-order alternative modes also change in shape with propagation distance \(z\) (unlike the "standard" modes) because of the change in \(\tilde{q}(z)\) with distance.

Either family of Hermite-gaussian solutions, 54. or 67., is equally valid as a general basis set for analytic expansions of arbitrary optical beams. Stable laser resonators with spherical mirrors and neligible beam aperturing will generally have real eigenmodes that are much closer to the "standard" or \(\sqrt 2x/w\) family of Hermite-gaussian modes, and hence this form for the Hermite-gaussian solutions is much more widely used in the laser literature.

More general complex paraxial systems including soft gaussian apertures, such as we will discuss later, do lead to Hermite modes with complex arguments, like the "elegant" or \(\sqrt{jkx^2/2\tilde{q}}\) family, and this type of solution is now being more extensively considered.

In all situations we can develop astigmatic mode solutions with different fundamental parameters in the \(x\) and \(y\) coordinates if this seems useful.

 

 

6.  GAUSSIAN BEAM PROPAGATION IN DUCTS 

Gaussian beams in free space always remain gaussian, but do of course spread outward due to diffraction effects as they propagate. In a medium with a quadratic transverse variation of index of refraction, however, it becomes possible to trap and propagate a particular confined gaussian beam which neither spreads nor contracts with distance.

We have already discussed the ray-trapping properties of graded-index optical waveguides or ducts. Such ducts are also of substantial practical and analytical interest in gaussian beam optics as well as in ray optics.

The two topics are, in fact, essentially identical in concept and in results. 

 

Gaussian Beam Propagation in Ducts

Suppose again that the index of refraction \(n(r)\) in a duct has a radial (or transverse) variation given by 

\[\tag{72}n(r)=n_0-\frac{1}{2}n_2r^2.\]

The wave equation 1. in a medium with a quadratic transverse variation such as this must then be expanded to the form 

\[\tag{73}\left[\nabla^2+\omega^2\mu\epsilon[1-n_2(x^2+y^2)]\right]\tilde{E}(x,y,z)=0.\]

Converting this to the paraxial approximation in the same form as used in deriving Equation 5. then gives 

\[\tag{74}\nabla^2_{xy}-k^2n_2(x^2+y^2)-2jk\frac{\partial}{\partial z}\tilde{u}(x,y,z)=0.\]

 

 

 

The reader can verify that a stable solution to this equation is given by the confined or trapped gaussian beam 

\[\tag{75}\tilde{u}(x,y,z)=\tilde{u}_0\text{exp}\left[-\frac{x^2+y^2}{w^2_1}+j\frac{\lambda z}{w^2_1}\right],\]

where the spot size \(w_1\) is given by 

\[\tag{76}w^2_1=\frac{\lambda}{\pi\sqrt{n_2}}.\]

This solution thus represents a stable trapped gaussian eigenmode of fixed diameter in the waveguide or duct. Note that the quadratic exponent for this wave is purely real, i.e., the wavefronts in this guided beam are exactly plane waves. Higher-order modes with Hermite-gaussian or Laguerre-gaussian form can also propagate in the same duct. 

 

Physical Interpretation

One way of understanding this confined mode is the following. In a duct with an index variation like Equation 72., each axial segment of length \(\Delta  z\) is like a thin lens with focal length \(f=1/n_2\Delta z\), as illustrated in Figure 11.

For a gaussian beam with spot size \(w_1\) as given in the preceding, the effects of diffraction spreading in each unit length are just canceled by this focusing effect, so that the beam size remains constant.

Note again that the steady-state gaussian spot size \(w_1\) varies inversely as the strength \(n_2\) of the transverse index variation—the stronger the focusing the smaller the steady-state beam profile that will be propagated in the duct.

This gaussian eigenmode also acquires a small added phase shift per unit length, over and above the \(e^{jkz}\) factor, that is expressed by the \(+j\lambda z/\pi w^2_1\) term.

This indicates that the \(z\)-directed propagation constant in the duct, call it \(k_d\), is not the on-axis value of \(k\) in the medium but rather a guided-wave value given by

\[\tag{77}k_d=k-\lambda/\pi w^2_1.\]

The added phase shift \(\psi(z)\) associated with the \(-\lambda z/\pi w^2_1\) per unit length in the duct is just given by \(d\psi(z)/dz=\lambda\pi w^2_1=1/z_R\) where \(z_R\) would be the Rayleigh range for the guided beam without the ducting effects.

This is in fact exactly the same as the derivative \(d\psi(z)/dz\) exactly at the waist for the Guoy effect which we will discuss in the following tutorial.

 

 

 

Beam Oscillations or Beam Scalloping in Ducts 

Suppose we introduce into such a duct a gaussian beam which does not match the gaussian eigenmode of the duct, either in spot size or in wavefront curvature. Suppose the input beam is initially smaller than the steady-state spot size \(w_1\).

The diffraction spreading for this smaller beam will then be larger than for the steady-state spot size, whereas the refocusing produced by each unit length of the duct will be the same.

The spot size will therefore begin to grow, and the gaussian beam will begin to spread with distance. 

As soon as its spot size becomes larger than the steady-state value, however, the opposite condition will prevail, and the beam will be refocused again. An input beam with a spot size \(w(z)\) either larger or smaller than the steady-state value wi will thus oscillate or scallop periodically inward and outward about the steady-state value in sausage-like fashion as it propagates down the guide, as shown in Figure 12.

This gaussian beam behavior is very much analogous to the oscillatory behavior of rays trapped in a stable duct, and to first order the beam oscillations will have the same oscillation period as the stable ray solutions we derived earlier.

As an equally valid alternative explanation, we can say that the mismatched input beam excites in the duct a mixture of the lowest and higher-order eigenmodes of the duct.

These eigenmodes then propagate with slightly different phase velocities, since the added phase shift term \(\lambda z/\pi w^2_1\) is different for the lowest-order and for the higher-order eigenmodes.

The periodic scalloping behavior then represents a beating phenomena in the axial direction as the higher-order modes propagate and interfere with different phases at difference distances along the duct.

 

Ducts in Real Systems

The paraxial wave solutions in a duct will be exactly gaussian only if the index variation about the axis is exactly quadratic—and a purely quadratic decrease of the index with radius cannot continue indefinitely.

The Hermitegaussian solutions for propagation in a duct are thus generally approximations to the real modes, although they will be quite good approximations if the index variation is in fact approximately quadratic out to at least a few spot sizes \(w_1\) of the trapped gaussian beam (as is quite often the situation).

If this approximation fails, or if the index variation is other than quadratic, then other mathematical forms for the trapped modes must be sought; and there is a very large published literature on such trapped modes, especially of course for the modes in simple optical fibers which have step variations rather than continuous variations in index of refraction.

Transverse index variations leading to duct-like effects occur naturally on a random basis in many situations—for example, in laser rods, in inhomogeneous optical materials, in optical (or radio-wave) propagation through the atmosphere, or in sound-wave propagation through the oceans.

Poor quality laser rods in the early days of lasers, in fact, often exhibited highly irregular transverse beam patterns, with many small spots across the face of the rod, due to strong ducting effects at random locations across the rod cross-section. More controlled ducting effects are also now deliberately produced to create waveguiding effects in optical fibers or in laser rods manufactured under such trade names as \(\text{GRIN}\) (graded refractive index) or \(\text{SEL}\)-\(\text{FOC}\) rods.

Note that the spot size \(w_1\) of the trapped beam in a duct will almost always be much smaller than the half-width \(1/n_2^{1/2}\) of the duct itself, unless the duct is extremely small.

Hence it will only be the quadratic or \(n_2\equiv-\partial^2n/\partial r^2\) variation of the index close to the axis that will usually be important, even if the index profile no longer remains purely quadratic at larger radii.

Larger ducts can also trap an entire family of higher-order Hermite-gaussian or Laguerre-gaussian modes such as we have analyzed in the previous chapter. We will develop a more general theorem of gaussian ducts in future tutorial, including both quadratic index and loss variations, by using a complex \(\text{ABCD}\) matrix approach

 

 

7.  NUMERICAL BEAM PROPAGATION METHODS 

Hermite-gaussian modes, particularly the lowest-order gaussian beams, provide extremely useful tools for analyzing optical beam propagation in simple situations, especially in low-loss stable optical resonators, and in other situations where the physical modes of the problem are close to Hermite-gaussian in character, and where the effects of diffraction by hard-edged apertures are negligible or entirely absent.

We can even, if necessary, expand any arbitrary optical beam as a summation of Hermite-gaussian modes; and then calculate the propagation of the arbitrary beam by calculating the propagation of the individual Hermite-gaussian modes. 

There are many situations in real optical systems, however—such as unstable resonator problems, for example—where aperture diffraction effects play a major role.

We must then treat the effects of edge diffraction, and analyze the propagation of beams with rather arbitrary and irregular amplitude and phase profiles. Numerical calculation methods play a large role in such analyses.

The most efficient numerical methods then generally center around the use of Huygens' integral together with various "fast transform" numerical methods.

In this section we will review briefly some of the analytical and numerical tools that become important in handling these more messy situations that arise in real-world problems.

 

Paraxial Wave Propagation: The Finite Difference Approach 

We have already mentioned the possibility of calculating the forward propagation of an arbitrary optical beam by writing the paraxial wave equation in the form

\[\tag{78}\frac{\partial\tilde{u}(s,z)}{\partial z}=\frac{j}{2k}\nabla^2_t\tilde{u}(s,z),\]

and then integrating this equation forward numerically using so-called finite difference methods, first to calculate the transverse derivative of the known wavefunction at one plane, and then to step forward in \(z\) to the next plane. 

This finite-difference approach to paraxial wave propagation can be of some practical usefulness for calculating beam propagation through inhomogeneous regions, such as perturbed atmospheres or problems involving thermal blooming or an inhomogeneous laser medium. Even in such inhomogeneous situations, however, fast transform methods, properly applied, are probably still superior.

 

Huygens' Integral: Fourier Transform Interpretation 

Huygens' integral provides another straightforward way to propagate an arbitrary optical wavefront from an input plane at \(z_0\) to any later plane \(z\). In doing this, it is particularly useful to note that the one-dimensional Huygens' integral written in rectangular coordinates has exactly the form of a convolution integral. 

That is, Huygens' integral as given in Equation 25. has exactly the form of a convolution (in \(x\)) of the input field \(\tilde{u}_0(x_0,z_0)\) with a spherical wavefunction \(\text{exp}[-j\pi x^2/(z-z_0)\lambda]\) in the form

\[\tag{79}\tilde{u}(x,z)=\tilde{u}_0(x_0)*\text{exp}[-j\pi x^2_0/(z-z_0)\lambda],\]

where the symbol \(*\) indicates the convolution operation. (We leave out the constant in front for simplicity.) 

The convolution of two functions, as in Equation 79., can be accomplished, however, according to Fourier transform theory, by \(\text{(i)}\) calculating the Fourier transform of each function individually; \(\text{(ii)}\) multiplying together the two Fourier transforms point by point to get a product transform; and then \(\text{(iii)}\) inverse Fourier transforming this product transform to get the desired convolution.

Efficient fast Fourier transform algorithms then provide a very practical way of doing the required transforms in order to numerically convolve two arbitrary functions such as Equation 79. on a digital computer. (In doing Huygens' integral calculations in this fashion, the spherical function corresponding to the kernel in Huygens' integral must be transformed only once and then stored.)

This approach is generally by far the most efficient way to evaluate the Huygens- Fresnel integral numerically, in order to calculate the propagation and diffraction spreading of an arbitrary optical beam in a numerical calculation, if the work is to be done in rectangular coordinates.

 

Alternative Fourier Transform Approach

As a slightly different approach, we can also rewrite Huygens' integral for one transverse dimension in the form 

\[\tag{80}\tilde{u}(x,z)=\text{exp}\left(\frac{-j\pi x^2}{L\lambda}\right)\sqrt{\frac{j}{L\lambda}}\int^{\infty}_{-\infty}\tilde{u}'_0(x_0,z_0)\times\text{exp}[j(2\pi/l\lambda)xx_0]dx_0,\]

where \(L\equiv z-z_0\) and \(\tilde{u}'_0(x_0,z_0)\) is a modified input function given by

\[\tag{81}\tilde{u}'_0(x_0,z_0)\equiv\tilde{u}_0(x_0,z_0)e^{-j\pi x^2_0/L\lambda}.\] 

But Equation 80. simply has the mathematical form of a Fourier transform between the variables \(x\) and \(x_0\). In this form, therefore, the Huygens-Fresnel propagation integral appears as a single (scaled) Fourier transform between the input and output functions \(\tilde{u}_0\) and \(\tilde{u}\).

This transform is applied, however, to the modified function \(\tilde{u}_0(x_0,z_0)\;\text{exp}(-j\pi x^2_0/L\lambda)\), and is followed by multiplication by another factor of \(\text{exp}(-j\pi x^2/L\lambda)\).

Applying a fast Fourier transform algorithm directly to the evaluation of Equation 80. provides another related but different way to do the same propagation calculation, using now a single Fourier transform.

However, this transform is now applied to a more complex input function, because of the additional spherical wave factor \(\text{exp}(-j\pi x^2_0/L\lambda)\). The total amount of numerical work seems to come out about the same for either approach in most practical situations.

 

Fourier Transforms and Gaussian Beams

Huygen's integral in the Fresnel approximation thus has the mathematical form either of a convolution of the input wavefront against a spherical wavefunction, or of a Fourier transform of the input wavefront multiplied by a spherical wavefront. Suppose we put a gaussian-spherical beam as the input into either of these mathematical forms.

The reader should then know, or learn, that the convolution of a gaussian function with another (possibly complex) gaussian always gives still another gaussian.

A gaussian beam passing through the convolution process of Equation 79. will thus always come out again gaussian, as we already know. 

To express this same point in an alternative form, the Fourier transform of a gaussian function is always another gaussian transform (and in fact the Fourier transform of a Hermite-gaussian function is always another Hermite-gaussian function of the same order).

All of the gaussian beam properties we have derived earlier and will discuss later are thus deeply imbedded in the self-transforming properties of generalized complex gaussian functions, and their higher-order Hermite-gaussian extensions. 

 

Paraxial Plane Waves and Transverse Spatial Frequencies

The mathematical procedures that are employed in evaluating the Huygens- Fresnel integral by Fourier transform methods can be given a simple and graphic physical interpretation in terms of an expansion of the optical beam in a set of infinite plane waves traveling in slightly different directions.

The Fourier transforms that are calculated in the convolution procedure correspond, in fact, to the transformation of the beam profile into a "spatial frequency" domain, or into a

 

 

 

"\(k\)-vector space" of plane waves traveling at different directions with respect to the \(z\) axis, as illustrated earlier in Figure 2. 

Because this spatial-frequency or \(k\)-vector viewpoint is very graphic and because it may give useful physical insights into the nature of beam propagation, we will rederive this plane wave description in some detail in the following paragraphs, even though the final results we obtain will be entirely identical to what we have already obtained merely by applying Fourier theorems to the Huygens- Fresnel integral.

To carry out this derivation, we consider as fundamental building blocks a set of infinite plane waves of the form

\[\tag{82}\tilde{u}_{pw}(x,y,z)\equiv\text{exp}[-j\bf k\cdot\boldsymbol r]=\text{exp}[-j(k_xx+k_yy+k_zz)].\]

The propagation vector \(k=(k_x,k_y,k_z)\) for any such plane wave traveling at angles \(\theta_x\) and \(\theta_y\) with respect to the \(z\) axis in the \(x,z\) and \(y,z\) planes has transverse components which we can write (see Figure 13.) as

\[\tag{83}k_x=k\sin\theta_x\equiv2\pi s_x\quad\text{and}\quad k_y=k\sin\theta_y\equiv 2\pi s_y.\]

The quantities \(s_x\equiv(k/2\pi)\sin\theta_x\approx\theta_x/\lambda\) and \(s_y\equiv(k/2\pi)\sin\theta_y\approx\theta_y/\lambda\) are then the transverse spatial frequencies along the \(x\) or \(y\) axes for a plane wave traveling at a small angle \((\theta_x,\theta_y)\) away from the \(z\) axis toward the \(x\) or \(y\) directions.

The longitudinal \(k\)-vector component for a given plane-wave component can then be written, using the paraxial or Fresnel approximation, in the form

\[\tag{84}\begin{array}\;k_z=\sqrt{k_2-k_x^2-k^2_y}\approx k-(k^2_x+k^2_y)/2k\\=k-\pi\lambda(s^2_x+s^2_y).\end{array}\]

Each individual plane wave component, characterized by its angles \(\theta_x\) and \(\theta_y\), or by its spatial frequencies \(s_x\) and \(s_y\), will then have a \(z\) propagation given by

\[\tag{85}\tilde{u}_{pw}(x,y,z)-\tilde{u}_{pw}(x,y,0)\times\text{exp}[-jkz+j\pi\lambda(s^2_x+s^2_y)z].\]

Each plane-wave component thus travels with a slightly different propagation constant in the \(z\) direction, given (within the paraxial approximation) by the \(\pi\lambda(s^2_x+s_y^2)z\) factor.

 

Expansion as a Distribution of Plane Waves

We then assume that an arbitrary paraxial optical beam \(\tilde{u}(x,y,z)\) can be expanded into a plane-wave or spatial-frequency expansion in the form 

\[\tag{86}\begin{array}&\tilde{u}(x,y,z)=\int\int\tilde{U}_{pw}(s_x,s_y,0)\times\tilde{u}_{pw}(x,y,z)ds_xds_y\\\quad\quad\quad\quad\quad\quad\quad=\int\int\left[\tilde{U}_{pw}(s_x,s_y,0)e^{-jkz+j\pi\lambda(s^2_x+s^2_y)z}\right]\times e^{-j2\pi(s_xx+s_yy)}ds_xds_y\\=\int\int\tilde{U}_{pw}(s_x,s_y,z)\times e^{-j2\pi(s_xx+s_yy)}ds_xds_y,\end{array}\]

where \(\tilde{U}_{pw}(s_x,s_y,0)\) gives the complex amplitude at plane \(z=0\) of each planewave component in the beam having spatial frequencies \(s_x,s_y\), or traveling at angles \(\theta_x\approx\lambda s_x\) and \(\theta\approx\lambda s_y\).

In writing the second and third lines we have made use of the fact that each such plane-wave component propagates forward in \(z\) with a differential propagation constant given by

\[\tag{87}\tilde{U}_{pw}(s_x,s_y,z)=\tilde{U}_{pw}(s_x,s_y,0)\times e^{-jkz+j\pi\lambda(s^2_x+s^2_y)z},\]

so that each component of this spatial-frequency distribution rotates in phase by a slightly different amount as the beam propagates forward.

At any arbitrary input plane \(z=z_0\) the beam intensity pattern can thus be written as

\[\tag{88}\tilde{u}(x,y,z_0)=\int\int\tilde{U}_{pw}(s_x,s_y,z_0)\times e^{-j2\pi(s_xx+s_yy)}ds_xds_y.\] 

But this expression is exactly a two-dimensional Fourier transform between the transverse spatial coordinates \(x,y\) and the spatial frequencies \(s_x,s_y\).

If we know the input field distribution \(\tilde{u}(x,y,z_0)\) at plane \(z_0\), therefore, we can evaluate the spatial frequency distribution \(\tilde{U}_{pw}(s_x,s_y,z_0)\) by carrying out the inverse Fourier transformation given by 

\[\tag{89}\tilde{U}_{pw}(s_x,s_y,z_0)=\int\int\tilde{u}(x,y,z_0)\times e^{+j2\pi(s_xx+s_yy)}dxdy.\]

Having transformed \(\tilde{u}(x,y,z_0)\) into the spatial-frequency domain, we can then propagate this spatial-frequency distribution forward (or for that matter backward) to any other plane \(z\) by multiplying it by the phase shift factor

\[\tag{90}\tilde{U}_{pw}(s_x,s_y,z)=\tilde{U}_{pw}(s_x,s_y,z_0)\times e^{-jk(z-z_0)+j\pi\lambda(s^2_x+s^2_y)(z-z_0)}.\]

The field distribution \(\tilde{u){x,y,z)\) at the second plane can then be evaluated from the second Fourier transformation

\[\tag{91}\tilde{u}(x,y,z)=\int\int\tilde{U}_{pw}(s_x,s_y,z)\times e^{-j2\pi(s_xx+s_yy)}ds_xds_y.\]

Propagating any arbitrary paraxial wavefunction from plane \(z_0\) to plane \(z\) in free space is thus carried out using two Fourier transformations plus one simple multiplication step.

Note also that the \(e^{-jk(z-z_0)}\) term is just the on-axis or plane-wave phase shift, whereas the \(e^{j\pi\lambda(s^2_x+s^2_y)(z-z_0)}\) factor gives the differential phase rotation of each individual spatial frequency component. 

This whole approach, from Equation 86. to Equation 91., is nothing more than a physical reinterpretation of the convolution or double-Fourier-transform approach to the evaluation of Huygens' integral which we discussed in connection with Equation 79.

Looking at it in spatial-frequency terms, however, emphasizes again the importance of the small-angle or Fresnel approximation and the quadratic spatial-frequency dependence which this produces in all the exponents.

 

Huygens' Integral in Cylindrical Coordinates

Huygens' integral 19. for free space can also be written in cylindrical coordinates \((r,\theta,z)\), leading to the result 

\[\tag{92}\tilde{u}(r,\theta)=\frac{j}{L\lambda}\int^\infty_0 r_0dr_0\int^{2\pi}_0\tilde{u}_0(r_0,\theta_0)\times\\\text{exp}\left\{-j(\frac{\pi}{L\lambda})\left[r^2+r^2_0-2rr_0\cos(\theta-\theta_0)\right]\right\}d\theta.\]

Suppose the wavefunction has \(m\)-th order azimuthal symmetry, so that we can separate the radial and azimuthal variables in the form

\[\tag{93}\tilde{u}(r,\theta)=\tilde{u}_m(r)\times e^{\pm jm\theta}\]

for both \(u\) and \(\tilde{u}_0\). Huygens' integral then reduces to the simpler form

\[\tag{94}\tilde{u}_m(r)=\frac{2\pi j^{m+1}}{L\lambda}\int^\infty_0 r_0\tilde{u}_0(r_0)e^{-j(\pi/L\lambda)(r^2+r^2_0)}J_m(s\pi rr_0/L\lambda)dr_0,\]

where \(J_m\) is the \(m\)-th order Bessel function.

Huygens' integral in this situation takes the form of a Fourier-Bessel transform, more commonly called a Hankel transform. A quasi "fast Hankel transform" is then available for the carrying out of numerical propagation and diffraction calculations on optical beams in cylindrical coordinates.