COMPLEX PARAXIAL WAVE OPTICS

This is a continuation from the previous tutorial - Stable two - mirror resonators

 

A very useful generalized form of paraxial optics has been developed in recent years. This generalized form can handle paraxial wave propagation not only in free space, and in simple lenses and ducts, but also in more general types of paraxial optical systems, including cascaded multielement optical systems (cascaded sequences of paraxial optical elements), and also systems having "soft apertures" or quadratic amplitude as well as phase variations about the axis.

This more general type of paraxial wave theory can be expressed in several mathematically equivalent forms. The approach that seems most convenient describes paraxial wave propagation entirely in terms of complex ray matrices or complex \(\text{ABCD}\) matrices.

This approach is very useful for handling complicated multielement optical resonators, as well as resonators with the very useful variable-reflectivity mirrors that are now being developed. 

In the present tutorial, therefore, we develop the basic theory for this extended form of paraxial optics using complex \(\text{ABCD}\) matrices. Topics covered include soft or gaussian apertures; complex \(\text{ABCD}\) matrices; the expansion of Huygens' integral in complex \(\text{ABCD}\) matrices; and the complex Hermite-gaussian modes that appear as the eigensolutions for generalized paraxial optical systems.

In the following chapter we will apply this approach to develop a generalized analysis of paraxial optical resonators, which will treat in one formalism not only conventional stable and unstable optical resonators, but also two new and useful classes of complex stable and unstable resonators. 

 

 

1.  HUYGENS' INTEGRAL AND ABCD MATRICES

The first step in deriving a generalized form of paraxial wave optics is to show how Huygens' integral for propagation through a cascaded series of conventional optical elements can be accomplished in one step, using nothing more than the overall \(\text{ABCD}\) matrix elements for that system.

We will give a tutorial derivation in this section which applies initially only to real \(\text{ABCD}\) matrices. The generalized

 

 

 

Huygens formula that we will derive will, however, in fact be valid for both real and complex \(\text{ABCD}\) matrices, as we will show later. 

 

Huygens' Integral in Free Space 

Huygens' integral in one transverse dimension for propagation through a distance \(\text{L}\) in free space can be written in the form 

\[\tag{1}\begin{align}tilde{u}_2(x_2)&=e^{-jkL}\int^\infty_{-\infty}\tilde{K}(x_2,x_1)\tilde{u}_1(x_1)dx_1\\&=\sqrt{\frac{j}{L\lambda}}\int^\infty_{-\infty}\tilde{u}_1(x_1)\text{exp}[-jk\rho(x_1,x_2)]dx_1,\end{align}\]

where the path length \(\rho(x_1,x_2)\) for an optical ray in free space traveling from position \(x_1\) at plane \(z_1\) to position \(x_2\) at plane \(z_2=z_1+L\) is given in the paraxial approximation by 

\[\tag{2}\rho(x_1,x_2)=\sqrt{L^2+(x_2-x_1)^2}\approx L+\frac{(x_2-x_1)^2}{2L}.\] 

Huygens' integral always involves this kind of optical path length \(\rho(r_1,r_2)\), sometimes called the eikonal function, from an optical source point at \(r_1\) to an observation or field measurement point at \(r_2\). The Huygens-Presnel kernel for wave propagation in free space thus takes on the form 

\[\tag{3}\tilde{K}(x_2,x_1)=\sqrt{\frac{j}{L\lambda}}\text{exp}\left[-j\frac{\pi(x_2-x_1)^2}{L\lambda}\right]\quad\left(\begin{array}&\text{free}\;\text{space}\\\text{propagation}\end{array}\right)\]

in one transverse dimension, or a product of two such kernels in two transverse dimensions.

 

Huygens' Integral Through a General Paraxial System 

Suppose we consider instead an input wavefunction \(\tilde{u}_1(x_1)\) traveling not through free space, but through a cascaded optical system containing an arbitrary collection of real paraxial optical elements (lenses, ducts, etc.) between the 

 

 

two planes \(z_1\) and \(z_2\), as illustrated in Figure 2. These cascaded optical elements can be characterized by an overall \(\text{ABCD}\) matrix or ray matrix. We will now show (or at least give a pretty good argument) that the eikonal function and the total Huygens' integral through this complete system can be written in one step, using only the overall \(\text{ABCD}\) matrix elements for the cascaded paraxial system. 

To do this, we will again invoke Huygens' principle of viewing the wave function \(\tilde{u}(x_1)\) at each point on the input plane \(z_1\) as a source of Huygens' wavelets, and evaluating the total amplitude \(\tilde{u}_2(x_2)\) at each point on the output plane \(z_2\) by adding up the contributions at \(x_2\) from all the input points \(x_1\), taking account of the path length or phase delay from each input point \(x_1\) to each output point \(x_2\).

We must therefore be able to calculate the overall path length, or the eikonal function \(\rho(x_1,x_2)\), from each point \(x_1\) on plane \(z_1\) in Figure 2 to each point \(x_2\) on plane \(z_2\) going through the complex paraxial system.

To do this, we can try to find the net path length \(\rho(x_1,x_2)\) for that particular optical ray which, starting out from transverse displacement \(x_1\) at plane \(z_1\), will emerge with transverse displacement \(x_2\) at plane \(z_2\), as in Figure 3.

This ray will, in general, not travel along a straight line path inside the system between the two end points, but will instead follow some more complicated trajectory within the optical elements making up the \(\text{ABCD}\) system.

What will be the optical path length \(\rho(x_1,x_2)\), or the eikonal function, for this ray between the input and output points, analogous to the free-space value used in the Huygens' integral in Equation 2?

 

Optical Path Lengths, and Fermat's Principle

We note first that if a ray is to enter at a specified point \(x_1\) and exit at a specified point \(x_2\), then from the ray relationship \(x_2=Ax_1+Bx'_1\) the input slope of this particular ray must be given by 

\[\tag{4}x'_1=\frac{x_2-Ax_1}{B},\] 

and the exit slope must be

\[\tag{5}x'_2=\frac{Dx_2-x_1}{B}.\]

 

 

The input ray may then be viewed as coming from an on-axis source point \(P_1\) located a distance \(R_1\) behind the input plane, as shown in Figure 3. Hence, \(R_1\) is given by 

\[\tag{6}\frac{R_1}{n_1}\equiv\frac{x_1}{x'_1}=\frac{Bx_1}{x_2-Ax_1}.\]

Upon passing through the cascaded paraxial optical system, this input ray is converted into an output ray which can again be associated with a spherical wave having an output radius of curvature \(R_2\) given by

\[\tag{7}\frac{R_2}{n_2}\equiv\frac{x_2}{x'_2}=\frac{Bx_2}{Dx_2-x_1}.\]

This output ray thus intersects the axis at point \(P_2\) located a distance \(R_2\) behind (or \(-R_2\) in front of) the plane \(z_2\) and we now have both \(R_1\) and \(R_2\) given in terms only of \(x_1\) and \(x_2\) and the \(\text{ABCD}\) matrix elements.

Now in optical terms, the points \(P_1\) and \(P_2\) in Figure 3 are conjugate points, or object-image points, in the sense that all rays leaving the axis from \(P_1\) will be focused back to the axis at \(P_2\), or vice versa.

But Fermat's principle then says that "all rays connecting two conjugate points must have the same optical path length between these two points." Applied to our system, Fermat's principle requires that

\[\tag{8}\left[\begin{align}\text{ray}\;\text{path}\;\text{from}\;P_1\;\text{to}\;P_2\\\text{through}\;X_1\;\text{and}\;X_2\end{align}\right]\equiv\left[\begin{array}&\text{ray}\;\text{path}\;\text{from}\;P_1\;\text{to}\;P_2\\\text{along}\;\text{the}\;\text{optical}\;\text{axis}\end{array}\right],\]

where we use \(X_1\) and \(X_2\) to denote the points on planes \(z_1\) and \(z_2\) where the off-axis rays intersect and leave the \(\text{ABCD}\) system.

Suppose the total optical path length through the paraxial system, from plane \(z_1\) to plane \(z_2\), for a ray traveling exactly on the axis is denoted by \(L_0\). In general this length will be given by a summation like

\[\tag{9}L_0=\sum_in_iL_i,\]

where each individual element inside the system has physical thickness \(L_i\) and index of refraction \(n_i\). If we also assume for complete generality that there are different indices of refraction \(n_1\) in the region before \(z_1\) and \(n_2\) in the region after \(z_2\), then the total optical path length along the optical axis from point \(P_1\) to \(P_2\) is given by

\[\tag{10}P_1P_2\equiv n_1R_1+L_0-n_2R_2.\] 

(Again, the minus sign is associated with the sign convention for \(R_2\).) 

On the other hand the total distance going along the off-axis ray trajectory through points \(X_1\) and \(X_2\) is given, in the paraxial approximation, by 

\[\tag{11}\begin{align}P_1X_1X_2P_2&=P_1X_1+X_1X_2+X_2P_2\\&=n_1(R^2_1+x^2_1)^{1/2}+\rho(X_1,x_2)-n_2(R^2_2+x^2_2)^{1/2}\\&\approx n_1\left(R_1+\frac{x^2_1}{2R_1}\right)+\rho(x_1,x_2)-n_2\left(R_2+\frac{x^2_2}{2R_2}\right),\end{align}\]

where \(\rho(x_1,x_2)\) is the path length inside the \(\text{ABCD}\) system. Fermat's principle then requires that the total path length from point \(P_1\) to conjugate point \(P_2\) going along either path be the same, i.e., \(P_1P_2=P_1X_1X_2P_2\).

Equating 10 and 11 and using the paraxial approximations, along with 6 and 7 for \(R_1\) and \(R_2\), then gives for the desired eikonal function

\[\tag{12}\rho(x_1,x_2)=L_0-\frac{n_1x^2_1}{2R_1}+\frac{n_2x^2_2}{2R_2}=L_0+\frac{1}{2B}(Ax^2_1-2x_1x_2+Dx^2_2).\]

The optical path length \(\rho(x_1,x_2)\) through the \(\text{ABCD}\) system from point \(x_1\) to point \(x_2\) is thus equal to the on-axis distance \(L_0\), phis an added optical length which is, within the paraxial approximation, quadratic in the displacements \(x_1\) and \(x_2\), and otherwise involves only the overall \(\text{ABCD}\) matrix elements for the system.

This added length is the generalized eikonal function, or a generalized analog to the added optical distance \((x_2-x_1)^2/2L\) that appears in the free-space Huygens' kernel between a source point at \((x_1,z_1)\) and an observation point at \((x_2,z_2)\) when \(z_2-z_1=L\).

 

Huygens' Integral Through the General ABCD System

We can make therefore a (correct) guess that the Huygens' integral for wave propagation all the way through the entire paraxial system, from plane \(z_1\) to plane \(z_2\), can be written in one step in the same general form as for the free-space situation of Equation 1, namely 

\[\tag{13}\tilde{u}_2(x_2)=e^{-jkL_0}\int\int\tilde{K}(x_2,x_1)\tilde{u}_1(x_1)dx_1,\]

but now with the Huygens kernel in one transverse dimension given by 

\[\tag{14}\tilde{K}(x_2,x_1)=\sqrt{\frac{j}{B\lambda_0}}\text{exp}\left[-j\frac{\pi}{B\lambda_0}(Ax^2_1-2x_1x_2+Dx^2_2)\right]\;\left(\begin{array}&\text{arbitrary}\\\text{ABCD}\;\text{system}\end{array}\right)\]

with \(\lambda_0\) being the optical wavelength in free space. The scale factor \(\sqrt{j/B\lambda_0}\) in front of the kernel is inserted out of necessity to conserve power, and to make the general Huygens kernel agree with the free space result when the \(\text{ABCD}\) system consists only of free space.

The matrix element \(B\) plays the same role in this kernel as the length \(L=z_2-z_1\) in the simple free-space Huygens' integral.

This form for Huygens' integral says that, within the paraxial approximation, an arbitrary optical wave can be propagated through a complete paraxial optical system, including all diffraction effects, using knowledge only of the overall \(\text{ABCD}\) coefficients of the system. (This is assuming no significant apertures or stops within the optical system, between the input and output planes.)

A more rigorous derivation to be given later will verify that this is in fact true even for more general paraxial systems having complex-valued \(\text{ABCD}\) matrices. A still more general version, applying to nonorthogonal systems as well, has been stated without proof in Equation 83 of the earlier tutorial on ray optics and ray matrices.

If there are apertures within the system, it is necessary to apply Huygens' integral in separate steps from the input up to the first aperture, from that aperture on to the next aperture, and so on.

Note also that the on-axis optical length \(L_0\) of the system is not contained in the generalized Huygens' integral, and in fact is not given by the \(\text{ABCD}\) coefficients.

This length represents a separate and independent parameter, outside of the \(\text{ABCD}\) coefficients. Also, if the paraxial system is astigmatic, we will use different \(\text{ABCD}\) matrices in the \(x\) and \(y\) directions.

 

 

2.  GAUSSIAN BEAMS AND ABCD MATRICES

Gaussian beams are the "eigenfunctions of free space", as we have pointed out in this Gaussian Beam Propagation in Ducts tutorial, and they may be the eigenfunctions of general complex paraxial wave systems also—especially since the Huygens-Fresnel integrals for free space and for a general paraxial \(\text{ABCD}\) system are so similar in form. 

To explore this point, let us therefore consider next what happens if we transmit a gaussian (or Hermite-gaussian) optical beam through a multielement paraxial optical system that is described only by its real (or possibly complex) \(\text{ABCD}\) matrix, and see if a Hermite-gaussian beam in produces a Hermitegaussian beam out. 

For this purpose consider an input beam of the form

\[\tag{15}\tilde{u}_1(x_1)=\text{exp}\left(-j\frac{\pi x^2_1}{\tilde{q}_1\lambda_1}\right)\quad\text{where}\quad\frac{1}{\tilde{q}_1}\equiv\frac{1}{R_1}-j\frac{\lambda_1}{\pi w^2_1}.\]

(Note that the \(\lambda\) used in writing \(\tilde{q}\) is always, by our definition, the wavelength in the medium where the beam is currently located.) If we put such a gaussian beam into the generalized one-dimensional Huygens' integral given in Equation 14, we obtain the integral 

\[\tag{16}\tilde{u}_2(x_2)=\sqrt{\frac{j}{B\lambda_0}}\int^{\infty}_{-\infty}\text{exp}\left[-j\frac{\pi x^2_1}{\tilde{q}_1\lambda_1}-j\frac{\pi}{B\lambda_0}(Ax^2_1-2x_1x_2+Dx^2_2)\right]dx_1,\]

where \(\lambda_0\) in this formula is the vacuum or free-space wavelength.

This integral contains only linear and quadratic powers of \(x_1\), and thus is easily evaluated using the lemma ("Siegman's lemma")

\[\tag{17}\int^{\infty}_{-\infty}e^{-ax^2-2bx}dx=\sqrt{\frac{\pi}{a}}e^{b^2/a},\]

where \(a\) and \(b\) may in general be complex, the only requirement being that \(a\) have a slightly positive real part. The result of this integration is that the beam at the output plane will still be a gaussian beam, but now in the form

\[\tag{18}\tilde{u}_2(x_2)=\sqrt{\frac{1}{A+n_1B/\tilde{q}_1}}\text{exp}\left(-j\frac{\pi}{\tilde{q}_2\lambda_2}x^2_2\right),\]

where its complex \(\tilde{q}\) parameter will have been transformed according to the relationship

\[\tag{19}\frac{\tilde{q}_2}{n_2}=\frac{A(\tilde{q}_1/n_1)+B}{C(\tilde{q}_1/n_1)+D}\]

with \(n_2\lambda_2=n_1\lambda_1\). But this is exactly the same as the ray-matrix equation 28. In other words, the complex \(\tilde{q}\) parameter for a gaussian beam can be transformed through an arbitrary real or complex paraxial system by exactly the same rule as applies to the radius of curvature \(R\) for a spherical wave using purely geometric optics. 

Equation 19 is extremely useful. It permits a gaussian beam to be propagated through multiple paraxial elements in sequence, using only the cascaded \(\text{ABCD}\) matrices for those elements.

One can either step the q value through individual elements in sequence, one by one, using the individual \(\text{ABCD}\) matrices; or alternatively we can cascade-multiply all the \(\text{ABCD}\) matrices first, and then propagate the gaussian beam through the entire system in one step using the overall \(\text{ABCD}\) matrix.

We noted earlier that if we defined the reduced radius of curvature \(\hat{R}\) for a purely spherical wave at any plane by

\[\tag{20}\hat{R}(z)\equiv\frac{R(z)}{n(z)},\]

where \(R(z)\) is the real radius of curvature at that plane, then the fundamental law of geometric optics can be simplified to

\[\tag{21}\hat{R}_2=\frac{A\hat{R}_1+B}{C\hat{R}_1+D}\quad(\text{geometric}\;\text{ray}\;\text{optics}\;\text{only}).\]

Note again that this formula is strictly valid only for purely spherical waves, that is, only in the geometric optics limit. Equation 21 does not give correct results for the radius \(R\) of a gaussian beam, except in the limit as the spot size \(w\) approaches infinity. 

Obviously we can also define a reduced \(\tilde{q}\) parameter, \(\hat{q}\), at any plane \(z\) by the corresponding definition

\[\tag{22}\frac{1}{\hat{q}}\equiv\frac{n}{\tilde{q}}\equiv\frac{n}{R}-j\frac{n\lambda}{\pi w^2}=\frac{1}{\hat{R}}-j\frac{\lambda_0}{\pi w^2},\]

where \(R(z)\) is the real radius of curvature; \(w(z)\) is the real spot size; \(\lambda (z)=\lambda_0/n(z)\) is the wavelength in the medium at that plane; and \(\lambda_0\) is the optical wavelength in vacuum.

The paraxial wave transformation rule using the \(\text{ABCD}\) matrix elements and the q gaussian beam parameters then becomes

\[\tag{23}\hat{q}_2=\frac{A\hat{q}_1+B}{C\hat{q}_1+D}\quad(\text{full}\;\text{paraxial}\;\text{wave}\;\text{optics}).\]

By using the reduced \(\hat{q}\) values we can carry out all calculations using only the \(\text{ABCD}\) matrix elements and the vacuum wavelength, with the local index of refraction \(n(z)\) coming into the calculations only when we go from reduced to real variables.

Obviously the radius \(R\) and the "complex radius" \(\tilde{q}\) scale by the index \(n\) in going from real to reduced variables, whereas the spot size \(w\) remains unchanged.

We will make extensive use of this gaussian beam transformation rule in future sections. It remains valid even for complex-valued \(\text{ABCD}\) matrices, to be introduced shortly.

 

Gaussian Beam Amplitude Transformation

Note also that the complex amplitude coefficient in front of the gaussian beam is transformed in propagating through the \(\text{ABCD}\) system in the form (for one transverse dimension) 

\[\tag{24}\frac{\tilde{u}_2(x_2=0)}{\tilde{u}_1(x_1=0)}=\sqrt{\frac{1}{A+B/\hat{q}_1}}.\]

With some algebraic manipulation this can be converted for purely real-valued \(\text{ABCD}\) elements into the form (again, for one transverse dimension) 

With some algebraic manipulation this can be converted for purely real-valued \(\text{ABCD}\) elements into the form (again, for one transverse dimension)

\[\tag{25}\frac{\tilde{u}_2(x_2=0)}{\tilde{u}_1(x_1=0)}=\sqrt{\frac{w_1}{w_2}}\text{exp}(-j\psi/2),\]

where \(\psi\) is now the phase angle defined by the complex quantity

\[\tag{26}\frac{A+B/\hat{q}_1}{|A+B/\hat{q}_1|}\equiv\text{exp}(j\psi).\]

This phase shift is a generalization of the Guoy phase shift introduced earlier.

The ratio \(\sqrt{w_1/w_2}\) is just the amplitude scaling for a one-dimensional or cylindrical wave. 

For a beam that has gaussian variations in both transverse dimensions, the ratio \(\tilde{u}_2(0)/\tilde{u}_1(0)\) given in Equation 25 must be squared, taking account also of any astigmatism in the beam parameters (i.e., \(\tilde{q}_{1x},\;\tilde{q}_{1y})\) or in the optical system \((A_x,\;A_y,\)etc.).

The phase shift \(\psi/2\), or \((\psi_x+\psi_y)/2\) for the two-transverse-dimension situation, must then be added to the \(e^{-jkL_0}\) term which is normally left out of these calculations.

 

 

3.  GAUSSIAN APERTURES AND COMPLEX ABCD MATRICES

In this section we introduce the very important new concepts of gaussian or "soft" apertures and gaussian ducts, and their description as generalized or complex paraxial elements, described by complex-valued \(\text{ABCD}\) matrices. 

 

Gaussian Apertures 

We might first recall that a thin convergent lens of focal length \(f\), gives an optical wave an added quadratic phase shift of the form 

\[\tag{27}\tilde{t}(x)\equiv\frac{\tilde{u}_2(x)}{\tilde{u}_1(x)}=\text{exp}\left(+j\frac{\pi x^2}{f\lambda}\right)\]

and is represented by an ABCD matrix of the form

\[\tag{28}\left[\begin{array}&A&B\\C&D\end{array}\right]=\left[\begin{array}&1&0\\-1/f&1\end{array}\right]\quad\left(\begin{array}\text{ABCD}\;\text{matrix}\\\text{for}\;\text{thin}\;\text{lens}\end{array}\right).\]

Let us consider instead a thin optical element consisting of a "soft aperture" (Figure 4) which has a quadratic transversely varying wave-amplitude transmission of the form

\[\tag{29}\tilde{t}(x)\equiv\frac{\tilde{u}_2(x)}{\tilde{u}_1(x)}=\text{exp}\left(-\frac{a_2x^2}{2}\right).\]

This might be the transmission function through a thin apodized filter element, as in Figure 4, or the amplitude reflection function from a planar mirror with a radially varying reflectivity, often referred to as a variable-reflectivity mirror or \(\text{VRM}\).

We will commonly refer to an optical element with this sort of radially varying transmission or reflection function as a "gaussian aperture" or a "gaussian variable-reflectivity mirror." Note that the intensity transmission through this element will be \(T{x}=|\tilde{t}(x)|^2\), and that a two-transverse-dimensional version will have \(\tilde{t}(x,y)=\tilde{t}_x(x)\tilde{t}_y(y)\) along its transverse principal axes.

 

 

The transmisson of a centered gaussian beam through this element is then given (in one transverse dimension) by 

\[\tag{30}\tilde{u}_2(x)\equiv\text{exp}\left(-j\frac{\pi x^2}{\hat{q}_2\lambda_0}\right)=\tilde{t}(x)\times\tilde{u}_1(x)=\text{exp}\left(-\frac{a_2x^2}{2}-j\frac{\pi x^2}{\hat{q}_1\lambda_0}\right).\]

Assuming that \(a_2\) is a positive quantity, the gaussian spot size is reduced, and the complex gaussian beam parameter \(\hat{q}\) is transformed in passing through a gaussian aperture according to

\[\tag{31}\frac{1}{\hat{q}_2}=\frac{1}{\hat{q}_1}-j\frac{\lambda_0a_2}{2\pi}.\]

We can put this into a more systematic form by rewriting it as

\[\tag{32}\hat{q}_2=\frac{\hat{q}_1+0}{(-j\lambda_0a_2/2\pi)\hat{q}_1+1}=\frac{A\hat{q}_1+B}{C\hat{q}_1+D}.\]

A gaussian aperture thus seems to act, at least so far as the q parameter of a gaussian beam is concerned, like a complex paraxial element, with a complex-valued \(\text{ABCD}\) matrix given by

\[\tag{33}\left[\begin{align}&A&B\\&C&D\end{align}\right]=\left[\begin{array}&1&0\\-j\lambda_0a_2/2\pi&1\end{array}\right]\quad\left(\begin{array}\text{ABCD}\;\text{matrix}\;\text{for}\\\text{"gaussian}\;\text{aperture"}\end{array}\right).\]

The gaussian aperture (with quadratic variation in the real part of the exponent) is thus very much like a thin lens or curved mirror (with quadratic variation in the imaginary part of the exponent), except that its focal length appears to be a purely imaginary quantity.

 

Gaussian Aperture Plus Thin Lens

As a slightly more general case, we might suppose that the element under consideration combines a thin lens having focal length \(f\) with a quadratic transmission variation \(a_2x^2/2\) as in Equation 27. (This could equally well be a curved variable-reflectivity mirror, with a radius of curvature \(R=2f\) and the 

 

 

same amplitude variation.) The complex \(\text{ABCD}\) matrix then becomes 

\[\tag{34}\left[\begin{align}&A&B\\&C&D\end{align}\right]=\left[\begin{array}&1&0\\-1/f-j\lambda_0a_2/2\pi&1\end{array}\right]\quad\left(\begin{array}\text{thin}\;\text{lens}\;\text{plus}\\\text{gaussian}\;\text{aperture}\end{array}\right).\]

The complex-valued matrix element shows up only in the lower left corner of this particular matrix. One can easily see, however, that if this \(\text{ABCD}\) matrix is cascaded with other purely real matrices, the complex values can easily spread to all four elements of the overall \(\text{ABCD}\) matrix.

 

Complex ABCD Matrices

So far we have only shown that a complex \(\text{ABCD}\) matrix like that in Equation 34 provides one way to interpret the propagation law for q through a gaussian aperture. In the following sections we will prove more generally, however, one of the fundamental principles of generalized or complex paraxial optics — namely, that any paraxial optical system which includes a gaussian (that is, a quadratic-exponent) variation in amplitude transmission across the axis may be considered as a complex paraxial system, and described by a complex-valued \(\text{ABCD}\) matrix or ray matrix. 

It may not be clear what it means physically to have a ray passing through a complex-valued ray matrix—at least, not in a purely geometric ray-optic analysis. We will see shortly, however, that essentially every formula developed for real paraxial systems and \(\text{ABCD}\) matrices will hold equally well for complex paraxial systems described by complex \(\text{ABCD}\) matrices.

A transverse variation in transmission (or reflection) amplitude, as described in Equation 29, is the basic condition for a complex paraxial system. Such a transverse variation may arise physically from a transverse variation in loss or absorption, or it may equally well represent a transverse variation in laser gain, as illustrated in Figure 5. The gaussian aperture described in Equation 29, for example, may equally well have the more general form

\[\tag{35}\tilde{t}(x)=\text{exp}\left[-a_0-\frac{a_2x^2}{2}\right]=\tilde{t}_0\;\text{exp}\left[-\frac{a_2x^2}{2}\right].\]

The on-axis transmission may then either be to \(\tilde{t}_0 <1\) and hence \(a_0>0\), representing an on-axis loss which increases radially; or it may be \(\tilde{t}_0>1\) and \(a_0<0\), i.e., an on-axis gain which decreases radially.

The complex \(\text{ABCD}\) matrix remains

 

 

the same, since it involves only the quadratic part of the transverse variation in either situation.

 

Other Forms of Aperture Transmission 

The complex paraxial analysis developed in this chapter will be mathematically exact (within the limits of the overall paraxial approximation) for systems having only linear or quadratic transverse variations in the exponent, as given in Equation 35. (Linear variations in phase or amplitude across the beam represent tilt, misalignment, or a displacement of the center of the beam.)

This analysis can also be extended, however, at least as a first approximation, to any transversely varying system whose transmission has a quadratic variation to first order near the optical axis. 

That is, consider any "soft" aperture whose transmission varies with transverse coordinate, at least near the axis, in the approximate form 

\[\tag{36}\tilde{t}(x)\approx t_0\times(1-a_2x^2/2),\]

where the coefficient \(a_2\) is now defined by the transverse second derivative, i.e.,

\[\tag{37}a_2\equiv-\frac{1}{t_0}\left|\frac{d^2\tilde{t}(x)}{dx^2}\right|_{x=0}\]

evaluated at the axis. This aperture may be approximated to first order by a gaussian aperture, and by a complex \(\text{ABCD}\) matrix with the same value of \(a_2\).

The complex paraxial analysis will then remain a good approximation so long as the resulting mode solutions remain confined sufficiently close to the axis so that \(\left|\frac{1}{2}a_2x^2\right|\ll1\) across the main portion of the wave.

The transverse coefficient \(a_2\) may even be a negative number, representing a radially increasing transmission, i.e.

\[\tag{38}\tilde{t}(x)=\text{exp}\left[+\frac{|a_2|x^2}{2}\right],\]

as in Figure 6, and the basic analysis will still apply equally well, with appropriate changes of sign. Questions arise, of course, as to how one realizes a radially increasing transmission function in practice, at least over any very large radius; and we will also see shortly that serious mode instability problems arise

 

 

in systems with radially increasing transmisson functions. The basic analysis, however, still remains valid regardless of the sign of \(a_2\). 

 

Gaussian Ducts 

As a generalization of the gaussian aperture, or gaussian aperture plus lens, we next introduce the concept of a complex gaussian duct.

We mean by this an arbitrary length of a dielectric medium which may have in general gaussian transverse variations of either the index of refraction \(n(x)\) and/or the loss (or gain) coefficient \(a(x)\) about the optical axis.

Such a gaussian duct will also lead to a complex \(\text{ABCD}\) matrix which is a direct extension of the ray matrix for purely real gaussian ducts. 

A gaussian duct is thus a transversely inhomogeneous medium in which the refractive index and the absorption coefficient may both vary transversely, in the forms

\[\tag{39}n(x)=n_0-\frac{1}{2}n_2x^2\quad\text{and}\quad a(x)=a_0+\frac{1}{2}a_2x^2.\]

Note that the coefficients \(a\), \(a_0\) and \(a_2\) in this situation are loss coefficients or loss factors per unit length, rather than total loss factors \(a_0\) and \(a_2\) through a discrete aperture as in the previous section.

Such a medium is obviously a complex extension of the graded-index ducts that we have considered earlier. The complex propagation constant \(k\) in such a duct varies near the axis in the form

\[\tag{40}k^2(x)=k^2_0-k_0k_2x^2\quad\text{or}\quad k(x)\approx k_0-\frac{1}{2}k_2x^2,\]

where the coefficient \(k_2\) corresponds to the transverse second derivative

\[\tag{41}k_2\equiv-\frac{d^2k(x)}{dx^2}\bigg|_{x=0}=\frac{2\pi}{\lambda}\left(\frac{n_2}{n_0}+j\frac{\lambda a_2}{2\pi}\right)\]

measured on the axis of the duct. The presence of a small uniform background loss or gain ao will also give the on-axis propagation constant a (normally very small) imaginary part, so that &o is expanded to \(k_0-ja_0\).

Note again that \(k_0\) in this section means the k value on axis, at \(x=0\), and not necessarily the free-space or vacuum value.

 

 

 

 

Differential Analysis 

Consider a small axial segment \(dz\) of such a duct, as in Figure 8. The results of the preceding section say that this small segment can be considered as a gaussian aperture plus a thin lens, with a complex \(\text{ABCD}\) matrix given by 

\[\tag{42}\boldsymbol{M}(dz)=\left[\begin{align}&1&dz/n_0\\&0&1\end{align}\right]\times\left[\begin{array}&1&0\\-n_0\bar{\gamma}^2dz&1\end{array}\right]=\left[\begin{array}&1&dz/n_0\\-n_0\bar{\gamma}^2dz&1\end{array}\right].\]

That is, the segment is essentially the combination of a small increment of free space having length \(B=dz/n_o\), plus a thin lens-gaussian aperture combination with a complex focal power given by

\[\tag{43}-n_0\bar{\gamma}^2dz=-\left(n_2+j\frac{\lambda_0a_2}{2\pi}\right)dz\]

for the incremental length \(dz\).

Now let \(\boldsymbol{M}(z)\) be the \(\text{ABCD}\) matrix for a complex (but axially uniform) duct from some arbitrary initial plane \(z_0\) up the the plane \(z\), and let \(\boldsymbol{M}(z + dz)\) be the same matrix from \(z_0\) to \(z+dz\). From the rules for cascading matrices we can then write

\[\tag{44}\boldsymbol{M}(z+dz)=\boldsymbol{M}(dz)\times\boldsymbol{M}(z),\]

where \(\boldsymbol{M}(z)\) has matrix elements \(\boldsymbol{A}(z)\), B(z), etc., and \(\boldsymbol{M}(z+dz)\) has matrix elements \(A(z+dz)\), \(B(z+dz)\), and so on. If we multiply out this matrix product and take the limit as \(dz\rightarrow 0\), we will obtain the same four differential equations as we obtained in this tutorial PARAXIAL OPTICAL RAYS AND RAY MATRICES, except now with a complex argument \(\bar{\gamma}\). Hence the complex overall \(\text{ABCD}\) matrix is given by the same result we derived in Equations 15 and 18, namely, 

\[\tag{45}\left[\begin{array}&A&B\\C&D\end{array}\right]=\left[\begin{array}&\cos\bar{\gamma}(z-z_0)&(n_0\bar{\gamma})^{-1}\sin\bar{\gamma}(z-z_0)\\-n_0\bar{\gamma\sin\bar{\gamma}(z-z_0)}&\cos\bar{\gamma}(z-z_0)\end{array}\right],\]

except that the \(\bar{\gamma}\) parameters are now all complex. This is the general complex \(\text{ABCD}\) matrix for a general gaussian dust. 

The matrix solution in Equation 45 is obviously the complex generalization of the purely real ray matrix results derived earlier for a quadratic-index duct. Since \(\bar{\gamma}\) will in general now be complex (except for the limiting situation of \(a_2\equiv\;0\;\text{and}\;n_2>0)\), the cosines and sines will now have complex arguments, and must be interpreted as combinations of trigonometric and hyperbolic functions according to the usual rules.

As one limiting situation, suppose the length \(\Delta z\equiv(z-z_0)\) of such a duct goes to zero, but the strength \(\bar{\gamma}\) of the transverse variation increases in such a way that \(\bar{\gamma}^2\Delta z\) remains finite. The \(\text{ABCD}\) matrix then simplifies to

\[\tag{46}\left[\begin{array}&A&B\\C&D\end{array}\right]\approx\left[\begin{array}&1&0\\-n_0\bar{\gamma}^2\Delta z&1\end{array}\right]=\left[\begin{array}&1&0\\-(n_2+j\lambda_0a_2/2\pi)\Delta z&1\end{array}\right].\]

This becomes just the combination of a thin converging lens (for \(n_2>0)\) and a thin gaussian aperture with quadratic coefficient \(z_2\) as discussed earlier.

We have chosen the definitions for \(n_2\) and \(a_2\) in this section such that \(n_2>0\) represents an index maximum on axis, whereas \(a_2>0\) represents a transmission maximum (i.e., a loss minimum, or a gain maximum) on axis.

These are the conditions which usually lead to confined and stable resonator or waveguide modes, and hence they are usually the conditions of most interest. Either \(n_2\) or \(a_2\) could change sign, and all of Equations 27 to 45 would remain valid, provided that we keep proper track of signs and complex phase angles in the various complex-valued formulas. 

In an orthogonal but astigmatic system a separate matrix like this can obviously be written for each of the transverse principal axes, with the appropriate and possible different values of \(\bar{\gamma}\) for the \(x\) or \(y\) transverse axes.

 

 

4.  COMPLEX PARAXIAL OPTICS

We have argued, though not rigorously proved, in the previous two sections that any general paraxial optical element can be described by a complex-valued \(\text{ABCD}\) matrix.

The most general cascaded, orthogonal, paraxial optical system can then be described by its overall \(\text{ABCD}\) matrix, which will be the matrix product of the complex matrices of the individual elements.

Imaginary parts of the \(\text{ABCD}\) matrix elements will appear in such systems in general whenever there are components that have a quadratic transverse variation of gain or loss, such as quadratic apertures or ducts. 

The primary objectives in this section are to extend this complex ray matrix picture to Huygens' integral, and then to prove in a general way that this picture is indeed a valid representation for any complex-valued quadratic or paraxial system.

 

Complex Huygens' Integral 

We developed in the first section of this chapter a generalized form of Huygens' integral for arbitrary real \(\text{ABCD}\) systems, namely, (in one transverse 

\[\tag{47}\tilde{u}_2(x_2)=\sqrt{\frac{j}{B\lambda_0}}\int^\infty_{-\infty}\tilde{u}_1(x_1)\text{exp}\left[-j\frac{\pi}{B\lambda_0}(Ax^2_1-2x_1x_2+Dx_2^2)\right]dx_1\]

in which \(\tilde{u}_1(x_1)\) is the complex input wave at plane \(z_1\); \(\tilde{u}_2(x_2)\) is the output wave at plane \(z_2\); and the \(\text{ABCD}\) elements are the overall complex ray matrix elements from \(z_1\) to \(z_2\).

This equation does not include the on-axis phase shift factor \(e^{-jk(z_2-z_1)}\). This total phase shift factor will be given in general by a sum of individual terms like

\[\tag{48}k(z_2-z_1)=\sum_ik_i(z_i-z_{i-1}),\]

where \(k_i(z_i-z_{i-1})\) is the on-axis optical path length through the \(i\)-th element. This axial phase shift factor is not uniquely determined by the \(\text{ABCD}\) matrix elements, and must be separately evaluated if its exact value is required. 

For simplicity we are writing Huygens' integral and all of the corresponding equations throughout this section in one transverse coordinate only, assuming an orthogonal optical system.

It is also assumed throughout that there are no hard apertures that will cause significant diffraction effects anywhere inside the optical system beyond the input plane, other than the soft gaussian apertures represented by the complex parts of the \(\text{ABCD}\) matrices.

 

General Proof For Complex Paraxial Systems

We now want to prove in general that the complex form 47 for Huygens' integral, and all of the other results given so far in this chapter, will apply to any arbitrary complex paraxial system. 

In order to do this, we first note that any complex paraxial optical system can be made up of only two basic elements, namely, complex gaussian propagation segments, or "ducts," such as we have previously described; and curved dielectric interfaces, having a spherical radius of curvature \(R\), between such gaussian propagation segments.

All other paraxial elements can be formed from combinations and/or limiting situations of these two basic elements. (Note that if we were to use the reduced ray matrix notation discussed in an earlier chapter, even the planar dielectric interfaces would drop out, and we could formulate everything from ducts alone.)

The general validity of complex paraxial optics can then be established by

(a) Showing that the generalized Huygens' integral in the form given in Equation 47 is valid for both of these two basic elements individually; and 

(b) Showing that Huygens' integral in this form can be cascaded through multiple elements in sequence simply by using \(\text{ABCD}\) matrix multiplication.

If Huygens' integral is valid for complex elements, then all the other results that we have established using the general Huygens' integral must also be valid in general for complex \(\text{ABCD}\) elements. 

 

Huygens' Integral for Complex Ducts

Let us now proceed with these steps. First of all, the general definition of a quadratic propagation medium or duct is any region in which the complex propagation constant \(k\) has the quadratic transverse variation 

\[\tag{49}k^2(x)=k^2_0-k_0k_2x^2.\] 

Although the constants \(k_0\) and \(k_2\) are in general complex when gain or loss are present, the imaginary part of ko is usually neglected in the \(\text{ABCD}\) formulation; and the restriction that \(|k_2x^2|\ll|k_0|\) over the volume of interest is generally assumed.

If we combine this with the usual paraxial approximation that \(|\partial^2u/\partial z^2|\ll|k\partial u/\partial z|\), the paraxial wave equation for propagation in this quadratic duct then becomes

\[\tag{50}\left[\frac{\partial^2}{\partial x^2}-2jk_0\frac{\partial}{\partial z}-k_0k_2x^2\right]\tilde{u}(z,z)=0.\]

Any valid solution for wave propagation in a gaussian duct must satisfy this wave equation. 

Suppose we now take Huygens' integral in the form given in Equation 47, together with the complex \(\text{ABCD}\) matrix in the form 

\[\tag{51}\left[\begin{array}&A(z)&B(z)\\C(z)&D(z)\end{array}\right]=\left[\begin{array}&\cos\bar{\gamma}(z-z_0)&(n_0\bar{\gamma})^{-1}\sin\bar{\gamma}(z-z_0)\\-n_0\bar{\gamma}\sin\bar{\gamma}(z-z_0)&\cos\bar{\gamma}(z-z_0)\end{array}\right],\]

where

\[\tag{52}\bar{\gamma}^2=\frac{k_2}{k_0}=\frac{n_2}{n_0}+j\frac{\lambda a_2}{2\pi},\]

and substitute all of this into the extended paraxial wave equation 50 for the duct. We can then show—after a considerable amount of algebraic labor—that this combination indeed satisfies Equation 50 exactly in any gaussian duct.

The kernel of the complex Huygens' integral also reduces to the usual Huygens- Fresnel kernel—that is, to simple spherical wavelets emanating from each source point \(x_1\)—in the limit as \(|\tilde{\gamma}|\ll1\), so that the whole solution is also valid in the free-space limit. 

Huygen's integral as given at the beginning of this section, plus the general \(\text{ABCD}\) matrix as given in Equation 51, thus provide general solutions within any gaussian duct.

It remains only to prove that these results can be extended to dielectric interfaces, and that they cascade properly.

 

Huygens' Integral for Dielectric Interfaces

A spherically curved boundary or interface between any two optical elements, and in particular between any two sections of complex gaussian duct, is the other basic paraxial optical element.

Such an interface will have a basically real \(\text{ABCD}\) matrix that is given in general by 

\[\tag{53}\left[\begin{array}&A&B\\C&D\end{array}\right]=\left[\begin{array}&1&0\\(n_2-n_1)/R&1\end{array}\right]\]

where \(n_1\) and \(n_2\) are the on-axis propagation constants on the input and output sides of the interface, and \(R\) is the radius of curvature of the interface. The sign

 

 

convention for this curvature says that \(R\) is positive for a surface concave toward the input medium or medium #1. 

The \(\text{ABCD}\) matrix for an interface is basically purely real, because \(n_1\) and \(n_2\) are essentially real quantities.

Any small imaginary part to these on-axis propagation constants, due to an average gain or loss coefficient \(a_0\), can normally be neglected. 

We must now verify that the combination of Huygens' integral plus this interface \(\text{ABCD}\) matrix reproduces what actually happens at a dielectric interface. If we insert the matrix elements given in Equation 51 directly into Huygens' integral, this will lead to a singular integral, because \(B\equiv0\) for a discontinuous interface.

We can avoid this difficulty, however, by cascading a short length \(\Delta z\) of the appropriate dielectric medium on each side of the interface, so that the integral then becomes finite. 

When we then evaluate Huygens' integral across an interface in the limit as the length \(\Delta z\) of each attached segment goes to zero, we find that the kernel of Huygens' integral acquires a Dirac delta function character; and the effect on the incident wavefront \(\tilde{u}_1(x)\) predicted by Huygens' integral becomes in the limit simply

\[\tag{54}\tilde{u}_2(x)=\tilde{u}_1(x)\text{exp}\left[-j\frac{\pi(n_2-n_1)x^2}{R\lambda}\right].\]

But this is exactly the expected physical result for a curved interface.

In carrying out beam calculations with cascaded optical elements, it is thus important to remember that even though the optical fields themselves will be continuous across an interface between two different media, the interface still has an \(\text{ABCD}\) matrix as given in Equation 53, which must be included in the calculations.

The \(\tilde{q}\) value of a gaussian beam changes discontinuously in going through such an interface, for example, because the wavelength in the medium changes, even though the spot size of the beam is continuous. The radius of curvature changes because of Snell's law. 

 

Cascading Huygens' Integral

To complete our proof, we must verify that we can cascade Huygens' integral in the \(\text{ABCD}\) matrix form off Equation 47 using matrix multiplication of the complex \(\text{ABCD}\) elements.

Propagating through any sequence of complex ducts and interfaces can then be done by first multiplying out their cascaded \(\text{ABCD}\) matrices, and then applying Huygens' integral using the product matrix for the system. 

Suppose we want to transform an arbitrary input wavefunction \(\tilde{u}_0(x_0)\) from an input plane \(z_0\) to an intermediate wavefunction \(\tilde{u}_1(x_1)\) through one paraxial

 

 

element, using Huygens' integral with some set of complex matrix elements \(\boldsymbol{M}_1=(A_1,B_1,C_1,D_1)\); and then (without any intervening aperture) transform from the intermediate wavefunction \(\tilde{u}_1(x_1)\) to an output wavefunction \(\tilde{u}_2(x_2)\) as in Figure 10, again using Huygens' integral through a second set of matrix elements \(\boldsymbol{M}_2=(A_2,B_2,C_2,D_2)\). 

Suppose we do this in straightforward fashion, by writing the cascaded Huygens' integrals

\[\tag{55}\begin{array}\tilde{u}_2(x_2)&=\int^{\infty}_{-\infty}\tilde{K}_2(x_2,x_1)\tilde{u}_1(x_1)dx_1\\&=\int_{(-\infty)}^\infty\tilde{K}_2(x_2,x_1)\left\{\int^\infty_{-\infty}\tilde{K_1}(x_1,x_0)\tilde{u}_0(x_0)dx_0\right\}dx_1\\&=\int^\infty_{-\infty}\left\{\int^\infty_{-\infty}\tilde{K}_2(x_2,x_1)\tilde{K}_1(x_1,x_0)dx_1\right\}\tilde{u}_0(x_0)dx_0,\end{array}\]

where the transition from the second to the third lines requires interchanging the order of integration between \(dx_0\) and \(dx_1\).

Doing this calculation then involves quadratic terms and products in \(x_0\), \(x_1\) and \(x_2\) that arise from the exponents of the Huygens' kernels.

It therefore requires some algebraic exercise in completing squares and using "Siegman's lemma." 

But we will find that the final result is exactly the same as we would have gotten by first multiplying together the two matrices to get a product matrix \(\boldsymbol{M}_{20}\equiv\boldsymbol{M}_2\times\boldsymbol{M}_1\) using standard matrix multiplication; writing a Huygens' kernel \(\tilde{K}_{20}(x_2,x_0)\) in the usual fashion using these matrix elements; and then evaluating Huygens' integral in one step in the form

\[\tag{56}\tilde{u}_2(x_2)=\int^\infty_{-\infty}\tilde{K}_{20}(x_2,x_0)\tilde{u}_0(x_0)dx_0.\] 

(Note as usual that the matrices must be multiplied together in reverse order.) Huygens' integral in the generalized \(\text{ABCD}\) matrix form can thus always be cascaded through multiple sections simply by multiplying together the individual complex \(\text{ABCD}\) matrices in the appropriate (i.e., reversed) order.

There is one weak mathematical constraint on this proof. Convergence of the intermediate integration over \(x_2\) in Equation 55 does require a weak condition on the matrix elements involved, to ensure that the fields do not blow up radially at the intermediate plane if any negative gaussian apertures are present.

However, this does not appear to present a significant limitation for actual physical optical systems.

 

Conclusions

All the expressions for the \(\text{ABCD}\) matrices, Huygens' integral, the transformation rule for the gaussian q parameter, and in fact all of the other results we have derived in earlier sections and chapters for real \(\text{ABCD}\) matrices, thus apply equally well to complex paraxial systems and complex ABCD matrices.

The reader will note that the fundamental formula \(AD-BC=1\) also applies to all the basic building blocks of complex paraxial analysis, and hence to any complex cascaded system as well.

This general approach, employing complex \(\text{ABCD}\) matrices in the \(x\) and \(y\) directions separately, seems to provide the most general form for analyzing propagation of any generalized paraxial wave through any generalized or complex paraxial optical system. 

There are several other more sophisticated mathematical approaches, besides the complex \(\text{ABCD}\) matrix approach, through which these same generalized paraxial optical results can derived and expressed (see References).

Some of these alternative approaches express generalized paraxial optics in terms of generalized operators, whereas others employ complex ray functions, eikonal functions, or Green's functions.

The generalized \(\text{ABCD}\) approach of this section seems to the author be the simplest approach, both to learn and to apply, as we will attempt to demonstrate in the following sections and tutorials.

 

 

5.  COMPLEX HERMITE-GAUSSIAN MODES

The ordinary Hermite-gaussian modes derived in Chapter 16 form a set of "normal modes of free space." Any such Hermite-gaussian (or Laguerre-gaussian) mode will propagate through free space retaining its Hermite-gaussian form and changing only in its radius of curvature \(R(z)\) and its scale factor or spot size \(w(z)\).

The propagation of such Hermite-gaussian modes is entirely governed by the free-space propagation law \(\tilde{q}(z)=\tilde{q}(z_0)+z-z_0\). 

As a very general extension of this free-space gaussian beam theory, a family of generalized or complex Hermite-gaussian modes will be identified in this section as normal modes for any arbitrary complex paraxial optical system.

In this section we will develop the formalism and the propagation laws for these complex Hermite-gaussian modes of any order propagating through any complex paraxial \(\text{ABCD}\) system.

The propagation of these complex Hermite-gaussian modes can be completely described, in fact, by simple transformation rules very much like the transformation rules for free space, using only the complex \(\text{ABCD}\) elements of the system.

 

Complex Hermite-Gaussian Waves 

An \(n\)-th order complex Hermite-gaussian mode of the type we wish to consider may be written in its most general form, in one transverse coordinate, as 

\[\tag{57}\tilde{u}_n(x)=\tilde{a}_n\tilde{v}^nH_n\left(\frac{\sqrt{2}x}{\tilde{v}}\right)\text{exp}\left(-j\frac{kx^2}{2\tilde{q}}\right),\]

where the three parameters \(\tilde{a}_n\), \(\tilde{v}\) and \(\tilde{q}\) all become complex quantities in the most general situation.

The complex radius of curvature \(\tilde{q}\) is still defined, just as in the free-space situation, by

\[\tag{58}\frac{1}{\tilde{q}}\equiv\frac{1}{R}-j\frac{\lambda}{\pi w^2},\]

where \(R\) and \(w\) are the purely real radius of curvature and gaussian spot size. The Hermite polynomial part of Equation 57 contains, however, a new complex Hermite scale factor or "complex spot size" \(\tilde{v}\) which appears in the argument of the Hermite polynomials.

This parameter is basically a new independent complex parameter, representing a generalization on the real gaussian spot size \(w\).

We will henceforth refer to waves of the form in which \(\tilde{v}\) is purely real and equal to w as "real" or "ordinary" Hermite-gaussian waves.

The same functions with \(\tilde{v}\) complex, or even with \(\tilde{v}\) real but different from \(w\), we will refer to as "complex" or "generalized" Hermite-gaussian waves. The distinction between these "ordinary" and "complex" Hermite-gaussians is the primary topic of this section.

We have pointed out in this tutorial - Numerical Beam Propagation Methods that there is no meaningful difference between the ordinary and the complex Hermite-gaussian waves for the two lowest-order modes \(n=0\) and \(n=1\), since the factor \(\tilde{v}\) drops out of both of these.

For modes with \(n>1\), however, a complex argument in the Hermite polynomial produces a transverse phase as well as amplitude variation in the Hermite polynomial part of the function.

As a result, the phase front of the mode is no longer purely spherical, but has additional phase distortion or wrinkling due to the complex Hermite factor as well. The transverse mode pattern also changes shape with distance for complex Hermite-gaussians with \(n\geqslant 2\).

 

Propagation of Complex Hermites Through Complex Systems

Suppose a complex Hermite-gaussian mode function \(\tilde{u}_{1n}(x_1)\) in the general form given in Equation 57, with subscripts 1 on the quantities \(\tilde{a}_{1n},k_1,\hat{q}_1\) and \(\tilde{v}_1\), is sent into a general paraxial optical system characterized by a complex \(\text{ABCD}\) matrix. (Remember that \(\hat{q}_1\) is the reduced value of the q parameter, as defined earlier.)

Propagation of this wave through the complete paraxial system can then be calculated in one step by using the complex Huygens' integral of Equation 47 in terms of the complex \(\text{ABCD}\) elements. Use of the generating function for the Hermite polynomials, namely,

\[\tag{59}\text{exp}(2xt-t^2)=\sum^\infty_{n=0}\frac{H_n(x)t^n}{n!},\]

allows this Huygens' integral to be evaluated analytically. The output function for a complex Hermite-gaussian input mode function \(\tilde{u}_{1n}(x_1)\) propagating through a complex \(\text{ABCD}\) system is found to be still a complex Hermite-gaussian mode function of exactly the same order as the input, but with in general new and different values for the parameters \(\tilde{a}_{2n},\hat{q}^2\) and \(\tilde{v}_2\). 

The output value for the \(\hat{q}\) parameter after passing through the complex paraxial system is related to the input value by exactly the same transformation rule as for real \(\text{ABCD}\) systems, namely,

\[\tag{60}\hat{q}_2=\frac{A\hat{q}_1+B}{C\hat{q}_1+D}.\]

In addition, the input and output mode amplitudes are related by the simple rule

\[\tag{61}\frac{\tilde{a}_{2n}}{\tilde{a}_{1n}}=\left(\frac{1}{A+B/\hat{q}_1}\right)^{n+1/2}.\]

Finally, the input and output values of the complex \(\tilde{v}\) parameter are related by a new, third rule, namely,

\[\tag{62}\tilde{v}^2_2=(A+B/\hat{q}_1)^2\times\tilde{v}^2_1+j\frac{4B}{k_1}(A+B/\hat{q}_1).\]

This is a new basic rule governing the propagation of the complex spot size parameter \(\tilde{v}\) through a general complex \(\text{ABCD}\) system.

The three transformation laws 60 to 62 are the most general form of the propagation rules for a complex Hermite-gaussian function through a complex paraxial optical system.

They involve the complex (reduced) radius of curvature \(\hat{q}\), the complex spot size parameter \(\tilde{v}\), and the complex amplitude (and phase) parameter an.

In terms of these parameters any input Hermite-gaussian of order n transforms into another Hermite-gaussian of the same order after passing through any complex paraxial system. The on-axis phase shift factor \(\text{exp}[-jk(z_2-z_1)]\) is not included in this transformation, but will be the same for all modes of any order \(n\).

We can then show algebraically that each of the transformation rules in Equations 60 to 62 cascades properly through multiple \(\text{ABCD}\) systems if we use as the matrix elements for the cascaded system the matrix product of the matrices of the individual cascaded elements.

Thus, the cascading property of the complex \(\text{ABCD}\) matrices holds here as it does in all other situations.

In carrying out the evaluation of Huygens' integral with a general Hermitegaussian input, there is one mathematical condition that arises, namely, it is necessary for convergence to assume that the imaginary part of the input quantity \(k(A/B+1/\hat{q}_1)\) is negative.

For real \(k\) and real \(\text{ABCD}\) elements this means that the input and output spot sizes must both be positive and finite.

For the complex situation, especially in systems that might have both radially decreasing and radially increasing loss at different points along the system, the physical meaning of this restriction is somewhat obscure.

A reasonable hypothesis is that this mathematical restriction may mean that the input wave must be limited to values that keep the output spot size bounded and the beam energy finite everywhere within the overall system.

 

Applications of the Complex Hermite-Gaussian Functions

The complex-argument Hermite-gaussian functions of Equation 57 have not yet found wide application, in optics or elsewhere, and the properties of these functions have not yet become as familiar as the conventional real-argument Hermite-gaussians.

It may be helpful, therefore, to present a few examples to illustrate the application of these general propagation rules and their connection with various earlier analyses.

 

(1) The Complex Hermite-Gaussian Scale Factor 

The complex scale factor v appearing in the generalized Hermite-gaussian functions in Equation 57 is a complex generalization of the usual gaussian spot size w.

The physical relationship between \(w\) and \(v\) in the general situation is not entirely clear.

In fact, it may be that some other way of parameterizing these functions would lead to a clearer separation between \(w\) and \(\tilde{v}\).

In any case, the transformation rule given in Equation 62 for the parameter \(\tilde{v}\) may perhaps be made to appear somewhat more plausible as follows. 

The usual propagation law for the real gaussian spot size \(w\) through a purely real \(\text{ABCD}\) system may be written, using earlier gaussian beam formulas, as

\[\tag{63}w^2_2=(A+B/R_1)^2\times w^2_1+(2B/k_1)^2(1/w^2_1),\]

where \(w_1\) and \(R_1\) are the real gaussian parameters at the input to the system and the \(\text{ABCD}\) elements are assumed purely real.

With some algebraic manipulation, this result may also be recast in the form (again, for real matrices only)

\[\tag{64}w^2_2=(A+B/\hat{q}_1)^2\times w^2_1+(4jB/k_1)(A+B/\hat{q}_1).\]

This involves the complex factor \(\hat{q}_1\) but gives a real answer. The transformation rule 62 for \(\tilde{v}\) then appears to be a fully complex generalization of these expressions, applicable to the fully complex situation.

 

(2) Propagation Eigenmodes of a Complex Quadratic Medium

As one example of the use of the general propagation rules, we can search for the propagation eigenmodes of a complex quadratic medium, i.e., a gaussian duct.

That is, we can look for a set of complex Hermite-gaussian eigenmodes \(\tilde{u}_n(x,z)\) for which the transverse field distribution will be independent of the axial distance \(z\).

Formally this means we must look for input values of \(\hat{q}_1\) and \(\tilde{v}_1\) such that the transformed values of these quantities remain unchanged, or \(\tilde{q}(z)=\hat{q}_1\) and \(\tilde{v}(z)=\tilde{v}_1\), through an arbitrary length of complex gaussian duct. 

From the transformation rules for \(\hat{q}\) and \(\tilde{v}\) plus the matrix elements of a gaussian duct, we can find that the complex Hermite-gaussian parameter values for these eigenmodes are

\[\tag{65}\hat{q}^2_2(z)=\hat{q}^2_1=-\tilde{\gamma}^{-2},\]

and

\[\tag{66}\tilde{v}^2(z)=\tilde{v}^2_1=2/k_0\tilde{\gamma},\]

where \(k_0\) is the on-axis propagation constant, \(\tilde{\gamma}\) is the complex parameter characterizing the duct as defined in Equation 43, and \(\tilde{a}_{1n}\) is an arbitrary initial constant. The phase shift and attenuation factor for one of these eigenmodes traveling a distance \(z\) down the duct is then given by 

\[\tag{67}\frac{\tilde{a}_{n}(z)}{\tilde{a}_{1n}}=\text{exp}[+j(m+1/2)\tilde{\gamma}z.\]

These results are identical with results obtained previously by Marcuse, who solved the inhomogeneous wave equation directly.

However, we will also show in a later section that modes with \(m\geq2\) may be unstable in this system, so that only the two lowest eigenmodes may be physically useful.

 

(3) Complex Gaussian Beams in Purely Real \(\text{ABCD}\) Systems

A general complex beam with a complex scale factor \(\tilde{v}\) can be launched into a purely real ray matrix system. In working with the complex \(\tilde{v}\) parameter, it is sometimes useful to define a parameter \(\tilde{V}\) as the ratio of the complex to real spot sizes, i.e., 

\[\tag{68}\tilde{V}\equiv\frac{\tilde{v}}{w},\]

so that the generalized Hermite-gaussian beam has the form

\[\tag{69}\tilde{u}_n(x)=\tilde{a}_n[\tilde{V}w]^nH_n\left(\frac{\sqrt{2}x}{\tilde{V}w}\right)\text{exp}\left(-j\frac{kx^2}{2\hat{q}}\right).\]

When \(\tilde{V}=\pm1\) the complex Hermite-gaussian reduces to a conventional real Hermite-gaussian. 

We can then show that for a purely real \(\text{ABCD}\) system the propagation rule for \(\tilde{v}\) reduces to

\[\tag{70}\frac{\tilde{V}^2_2-1}{\tilde{V}^2_1-1}=\frac{A+B/\hat{q}_1}{A+B/\hat{q}^*_{1r}},\]

which can be converted to 

\[\tag{71}\tilde{V}^2_2-1=(\tilde{V}^2_1-1)\times\text{exp}[2j\angle(A+B/\hat{q}_1)],\]

where \(\angle(A+B/\hat{q}_1)\) means the phase angle of that quantity. This shows that any optical system with real matrix elements will have an output wave with \(\tilde{V}_2=\pm1\) if and only if the input wave has \(\tilde{V}_1=\pm1\).

Therefore if we start with a purely real gaussian beam there is no optical system with purely real matrix elements that will transform it into a complex gaussian beam. 

Conversely, once we have a generalized complex gaussian beam with \(\tilde{V}^2_1\neq1\), there is no optical system with real matrix elements that can transform it back into a purely real gaussian beam.

This can be done only by a complex \(\text{ABCD}\) system. These conclusions are meaningful only for mode indices \(n\geq2\), since, as we have noted before, there is no meaningful distinction between real or complex modes for \(n=0\) or \(n=1\).

As one specific example, consider the propagation of a complex Hermite-gaussian beam in free space. Suppose that the input plane of the optical system, \(z_1\), coincides with the waist of the beam.

The matrix elements of the system out to any plane z are then \(A=D=1\), \(B=z\) and \(C=0\), and we then find

\[\tag{72}\frac{\tilde{V}^2(z)-1}{\tilde{V}^2(0)-1}=\frac{1-j(z-z_1)/zR}{1+j(z-z_1)/zR}.\]

This result is essentially identical to that of Pratesi and Ronchi, who obtained it by solving the wave equation directly. 

 

Beam Expansions In Complex Hermite-Gaussians

Any field distribution \(\tilde{u}(x)\) that satisfies the wave equation in a paraxial optical system can be expanded using as a basis set the complex Hermite-gaussian functions

\[\tag{73}\tilde{u}_n(x)=\tilde{a}_n\tilde{v}^nH_n\left(\frac{\sqrt{2}x}{\tilde{v}}x\right)\text{exp}\left(-j\frac{kx^2}{2\hat{q}}\right),\]

where \(\tilde{a}_n\) is an appropriate normalization constant, and where \(\tilde{v}\) and \(\tilde{q}\) may be arbitrarily chosen. The series expansion has the usual form

\[\tag{74}\tilde{u}(x)=\sum^\infty_{n=0}c_n\tilde{u}_n(x).\]

The general orthogonality relation for Hermite polynomials is

\[\tag{75}\int^\infty_{-\infty}H_n(\sqrt{c}x)H_m(\sqrt{c}x)\text{exp}(-cx^2)dx=\sqrt{\frac{\pi}{c}}2^nn^!\delta_{nm}\quad R_e[c]>0.\]

Using this we can show that the complex Hermite-gaussian functions \(\tilde{u}_n(x)\) are biorthogonal to a set of functions \(\phi_n(x)\) given by

\[\tag{76}\phi_n(x)=\beta_nH_n\left(\frac{\sqrt{2}x}{\tilde{v}}\right)\text{exp}\left[+j\frac{kx^2}{2\hat{q}}-\frac{2x^2}{\tilde{v}^2}\right]\]

with the normalization constant given by

\[\tag{77}\beta_n=\left(\frac{2}{\pi}\right)^{1/2}\frac{1}{2^nn!\tilde{a}_n\tilde{v}^{n+1}}\]

in the sense that \(\int^\infty_{-\infty}\phi_n(x)\tilde{u}_(x)dx=\delta_{nm}\). The coefficients in the expansion are then given by

\[\tag{78}c_n=\int^\infty_{-\infty}\phi_n(x)\tilde{u}(x)dx\]

as in the usual expansion fashion.

If the input wave \(\tilde{u}_1(x)\) to a complex \(\text{ABCD}\) system is expanded in this fashion, then the output wave \(\tilde{u}_2(x)\) from the system can be found by propagating each Hermite-gaussian wave \(\tilde{u}_n(x)\) through the system using the transformation rules of the previous section, and then reexpanding \(\tilde{u}_2(x)\) at the output using the same \(c_n\) coefficients.

In other words the coefficients \(c_n\) in the expansion do not change as the wave passes through a general paraxial system. The Hermitegaussian basis functions themselves change, as given by the transformation rules.

This invariance is one of the useful properties of the complex Hermite-gaussian expansion. 

 

 

 

 

6.  COORDINATE SCALING WITH HUYGENS' INTEGRALS 

Once we have the complex \(\text{ABCD}\) elements, it is possible, for purposes of either analysis or numerical computation, to transform the propagation of any arbitrary beam through any arbitrary complex paraxial system into an equivalent propagation of a near-collimated beam through an equivalent length of free space.

We will give both a mathematical derivation and a physical explanation of this transformation in the present section. 

 

Transformation of the General Huygens' Integral

To demonstrate this mathematically, suppose an arbitrary optical beam of approximate width \(2a_1\), which may be either diverging or converging, is sent into an arbitrary complex \(\text{ABCD}\) system, as illustrated in Figure 11. For a complicated \(\text{ABCD}\) system, with multiple internal lenses, gaussian apertures and the like, this beam will then emerge with some different approximate width \(2a_2\), and with a wavefront which may be either diverging or converging, as illustrated in Figure 11. 

The transmission of this beam through the \(\text{ABCD}\) system is given by the same Huygens' integral as in the preceding sections, namely,

\[\tag{79}\tilde{u}_2(x_2)=\sqrt{\frac{j}{B\lambda_0}}\int^{a_1}_{-a_1}\tilde{u}_1(x_1)\text{exp}\left[-j\frac{\pi}{B\lambda_0}[Ax^2_1-2x_1x_2+Dx^2_2]\right]dx_1,\]

except that we now suppose that the limits of the integral are set to the values \(\pm a_1\), where we choose the input aperture width \(2a_1\) either to equal an actual aperture in the input plane, or else just wide enough to comfortably contain all the significant portion of the input beam (including adequate allowance for "spillover"). 

Suppose we then also chose a scale width \(2a_2\) at the output, where a2 is arbitrary, but is best made wide enough to contain all the significant portion of the output beam.

We then define the ratio of these widths as a magnification scaling factor \(M\equiv a_2/a_1\), where \(\text{M}\) may in fact be either greater or less than unity, depending upon what actually happens to the particular beam in going through the paraxial system.

Let us next define normalized transverse coordinates \(y_1\) and \(y_2\) at the input and output planes by

\[\tag{80}y_1\equiv\frac{x_1}{a_1}\quad\text{and}\quad y_2\equiv\frac{x_2}{a_2}=\frac{x_2}{Ma_1},\]

and also define transformed input and output wavefunctions by

\[\tag{81}\tilde{v}_1(y_1)\equiv a_1^{1/2}\tilde{u}_1(x_1)\text{exp}\left[-j\frac{\pi(A-M)x^2_1}{B\lambda}\right]\]

and

\[\tag{82}\tilde{v}_2(y_2)\equiv a_2^{1/2}\tilde{u}_2(x_2)\text{exp}\left[+j\frac{\pi(D-M)x^2_2}{B\lambda}\right].\]

Obviously, these transformations—which depend in part upon the choice of the magnification \(\text{M}\)—amount to multiplying the input and output wavefronts by a quadratic phase transformation, as in passing through a thin lens, plus also possibly multiplying them by a gaussian aperture transmission if the \(A,B\) or \(D\) elements are complex.

The physical effect of these transformations, as we will see in the following, is essentially to convert the input beam, whether diverging or converging, into a quasi collimated input beam, and then to convert the output from the analysis back into the appropriate diverging or converging beam.

We will also sometimes refer to this transformation as "extracting out the spherical portion of the input or output wavefront."

If we make these analytical transformations and plug them into the general Huygens' integral 79, this integral is then transformed into the very simple form

\[\tag{83}\tilde{v}_2(y_2)=\sqrt{jN_c}\int^1_1\tilde{v}_1(y_1)\text{exp}\left[-j\pi N_c(y_1-y_2)^2\right]dy_1.\]

But this is exactly the normalized form of Huygens' integral for transmission through a free-space region with a Fresnel number given by

\[\tag{84}N_c\equiv\frac{a_1a_2}{B\lambda}=\frac{Ma^2_1}{B\lambda},\]

where \(N_c\) is a new kind of "collimated Fresnel number" for this propagation calculation.

 

The Collimated Fresnel Number

The formal transformation from the original Huygens' integral to the normalized free-space form in Equation 83 is mathematically valid independent of the choice we make for the output scale factor \(a_2\) or the magnification \(\text{M}\), so that the values assigned to these quantities are quite arbitrary so far as the mathematical analysis is concerned.

The most meaningful choice, however, will be for us to select the value of \(a_1\) so that \(a_1\) either matches the actual input aperture, if there is one, or else is just large enough to match the full width of the input beam.

We should then choose \(\text{M}\) so that \(a_2\) is similarly large enough to just include all the significant portions of the output beam. 

 

 

 

If we do this, the normalized functions \(\tilde{v}_1(y_1)\) and \(\tilde{v}_2(y_2)\) will then both more or less fill up the range \(\pm1\) in the input or \(y_1\) and output or \(y_2\) coordinates, i.e., the optical beam will look more or less like a collimated beam confined to the range \(\approx\pm1\) in the dimensionless y coordinate.

After the quadratic phase (and possibly quadratic amplitude) variations are extracted out by the transformations given in Equations 80 through 82, therefore, the diffraction effects of propagating through a complex \(\text{ABCD}\) system are exactly the same as propagating a quasi collimated beam through an equivalent length \(B/M\) of free space, as expressed by the Fresnel number \(N_c\).

We therefore refer to this Fresnel number Nc as the "collimated Fresnel number" for the system. 

This collimated Fresnel number is a primary measure of the amount of numerical work needed to evaluate the wave propagation integral 79 or 83 using any of the numerical techniques discussed in Section GAUSSIAN SPHERICAL WAVES (see also Section APERTURE DIFFRACTION: RECTANGULAR APERTURES).

It depends on the paraxial optical system through the parameter \(\text{B}\) on the width \(2a_1\) of the input beam; and also on the divergence or convergence of the input beam, through the magnification \(\text{M}\) that is needed to match the input and output beam widths.

 

Physical Interpretation

The preceding mathematical analysis can be given a more physical interpretation as follows. Consider again the wave propagation through an arbitrary \(\text{ABCD}\) system as in Figure 11.

This propagation can always be modeled by the collection of elements shown in Figure 12, consisting of a (possibly complex) lens with (possibly complex) focal length \(f_1\) at the input end, which converts the input beam into a quasi collimated beam; a "magnifier" (i.e., a relay imaging system of zero effective length) which magnifies this collimated beam pattern transversely by a magnification \(M_1\); a free space section of (as yet unknown) length \(L\); a second "magnifier" of power \(M_2\); and an output lens of (possibly complex) focal length \(f_2\), which gives the output beam the correct output curvature. 

The total transformation through the equivalent model of Figure 12 can then be equated to the \(\text{ABCD}\) matrix of the actual paraxial system in Figure 11 by evaluating the matrix products

\[\tag{85}\left[\begin{array}&A&B\\C&D\end{array}\right]=\left[\begin{array}&1&0\\-1/f_2&1\end{array}\right]\left[\begin{array}&M_2&0\\1/M_2&0\end{array}\right]\left[\begin{array}&1&L\\0&1\end{array}\right]\left[\begin{array}&M_1&0\\1/M_1&0\end{array}\right]\left[\begin{array}&1&0\\-1/f_1&1\end{array}\right],\]

where the \(\text{ABCD}\) coefficients on the left side are those of the actual paraxial system. If we multiply this out, we get the relationships

\[\tag{86}\begin{align}&A=M\times\left[1-\frac{L_{eq}}{f_1}\right]\\&B=M\times L_{eq}\\&D=\frac{1}{M}-\frac{ML_{eq}}{f_2},\end{align}\]

where \(L_{eq}\equiv L/M^2_1\), and where \(\text{M}\) is now set equal to the total magnification \(M_1\times M_2\) through the model system, which we will presumably adjust to match the actual transverse magnification from input to output beams.

We can then see that for any given choice of \(M\), the three independent elements \(A,B\) and \(D\) of the actual system can be matched by three adjustable parameters \(f_1\), \(f_2\) and an "equivalent length" \(L_eq\equiv L/M^2_1\) which always appears in combination in Equations 86.

If we invert these relations, we find that the equivalent free space length and the input and output focal lengths in the model of Figure 12 are given in terms of the \(M\) value and the parameters of Figure 11 by

\[\tag{87}L_{eq}=\frac{L}{M^2_1}=\frac{B}{M},\quad\frac{1}{f_1}=\frac{M-A}{B}\quad\text{and}\quad\frac{1}{f_2}=\frac{1/M-D}{B}.\]

Moreover, the collimated free-space region in the center of Figure 12—which is the only region where Huygens' integral need be evaluated—has a width \(M_1a_1\) and a length \(L=M^2_1B/M\).

Hence the Fresnel number in this free-space portion of the model is given by

\[\tag{88}N_c=\frac{(M_1a_1)^2}{L\lambda}=\frac{a^2_1}{(L/M^2_1)\lambda}=\frac{Ma^2_1}{B\lambda}=N_c.\] 

In other words, the Fresnel number in the collimated region of the model is exactly the same as the result in the purely mathematical analysis in Equation 84, independent of how we divide the total magnification \(M\) between the input magnification \(M_1\) and the output magnification \(M_2\).

The collimated Fresnel number \(N_c\)—which determines the computational difficulty of the propagation task—is indifferent to this choice.

 

More On the Optimum Magnification

Suppose further that in the system of Figure 11 the \(\text{ABCD}\) matrix elements are purely real, or nearly so, so that we may with reasonable accuracy think about real rays going through the system.

In addition, suppose further that the beam coming into the system is more or less similar to a bounded spherical wave with an input width \(2a_1\) and an approximate radius of curvature \(R_1\).

The outer edge of this beam can then be delineated by an input ray with displacement \(x_1=a_1\) and slope \(x'_1=a_1/R_1\).

The corresponding output ray, which will delineate the approximate outer edge of the output beam, will than have a displacement \(x_2\approx a_2=Ax_1+Bx'_1=(A+B/R_1)a_1\).

Hence, the approximate or geometric magnification of the transverse width for this beam going through the system will be given by

\[\tag{89}M=\frac{a_2}{a_1}\approx A+B/R_1.\]

But if we plug this magnification into Equations 87, the resulting focal lengths for the input and output lenses will be

\[\tag{90}\begin{align}&\frac{1}{f_1}=\frac{M-A}{B}\approx\frac{1}{R_1}\\&\frac{1}{f_2}=\frac{1/M-D}{B}\approx\frac{1}{R_2}.\end{align}\]

In other words, with the commonsense geometric choice for \(M\), the input and output lenses do just convert from the external diverging or converging beam profiles to an essentially collimated profile internal to the model.

In other words, with the commonsense geometric choice for \(M\), the input and output lenses do just convert from the external diverging or converging beam profiles to an essentially collimated profile internal to the model.

As the simplest possible example of this beam transformation process, let us consider a uniform-amplitude spherical wave of initial radius \(R_0\) diverging away from a hard-edged input aperture of width \(2a_1\) and traveling forward through a real distance \(L\), as in Figure 13.

The \(\text{ABCD}\) matrix for the paraxial "system" in this situation will just be the free-space matrix

\[\tag{91}\left[\begin{array}&A&B\\C&D\end{array}\right]\left[\begin{array}&1&L\\0&1\end{array}\right],\]

and the magnification from input to output will be

\[\tag{92}M=1+L/R_0=\frac{R_0+L}{R_0}.\]

This diverging beam problem will then be formally equivalent to a collimated beam problem with a collimated Fresnel number given by

\[\tag{93}N_c=\frac{Ma^2_1}{B\lambda}=\frac{(L+R_0)a^2_1}{R_0L\lambda}=\frac{a^2_1}{L_{eq}\lambda}.\]

In other words, propagation of the uniform diverging beam through distance \(L\) will be formally equivalent to propagation of a uniform collimated beam of the same width \(a_1\), through an equivalent (shorter) length \(L_{eq}\) given by 

\[\tag{94}L_{eq}-\frac{B}{M}=\frac{R_0L}{R_0+L},\]

as in Figure 13. 

In physical terms, the Fresnel diffraction profile of the diverging beam at distance \(L\) will be exactly the same as the Fresnel profile of the collimated beam at the shorter distance \(L_{eq}\), except for transverse scaling and an underlying wavefront curvature in the former.

We might crudely say that the diverging beam in the upper part of the figure suffers less Fresnel diffraction because it 

 

 

keeps getting wider transversely, and thus diffraction effects act more slowly. If the beam is converging, and has a negative value for \(R_0\), then exactly the opposite is true.

In fact, as we know, following a converging beam down to its focus or caustic point is equivalent to following the equivalent collimated beam all the way out to \(L_{eq}\rightarrow\infty\)

Use of the coordinate transformation developed in this section makes it evident that propagation and aperture diffraction calculations for diverging and converging beams traveling different physical distances are entirely equivalent, if they are described by the same collimated Fresnel number \(N_c\).

In fact, \(N_c\) is the only parameter that is relevant for any situation in which a uniform intensity plane or spherical wave, of any curvature, comes through a slit or circular aperture of width or diameter \(2a\).

The use of this transformation also makes it possible to handle the diffraction calculations for any paraxial propagation problem, through any paraxial system, with a single computer program which solves the propagation problem only through free space.

Such a program is generally best written, as we have said several times, using some sort of fast Fourier or Hankel transform algorithm.

 

 

 

 

7. SYNTHESIS AND FACTORIZATION OF ABCD MATRICES 

The previous section briefly introduced the concepts of extracting out a quadratic phase and amplitude profile from an arbitrary wavefront, and of synthesizing a given complex \(\text{ABCD}\) system out of elementary paraxial elements.

These are very useful analytical techniques, and in this section we add a few additional notes on how they might be generalized and employed. 

 

Curved Reference Planes and Self-Canceling Lens Pairs

One technique for simplifying \(\text{ABCD}\) matrix problems and paraxial propagation calculations is to extract out a quadratic phase front from the actual optical beam by multiplying the actual beam by a given spherical wave function, as we did in the previous section.

In general we do this so as to convert a quasi spherical beam that is diverging or converging into a quasi collimated beam. From another viewpoint, this is exactly the same as measuring or observing the optical beam not on a transverse reference plane but on a spherical reference surface intersecting the optical axis at the given location \(z\). 

Another way of describing this same technique physically is to think of introducing into the \(\text{ABCD}\) system at any point of interest a self-canceling lens

 

 

 

pair—that is, a pair of thin lenses of equal but opposite focal power, as shown in Figure 14. The overall \(\text{ABCD}\) matrix through these two lenses is just the identity matrix, and so introducing this lens pair will do nothing to the overall performance of the \(\text{ABCD}\) system.

Choosing a reference plane located between the two lenses, however, is equivalent to observing the optical beam at that location on a spherically curved reference surface. 

The generalization of these ideas is to extract out both a quadratic or spherical phase and a quadratic or gaussian amplitude variation from the optical wavefront, by multiplying the actual wave by an arbitrary complex-quadratic exponential, which is to say, by an arbitrary gaussian aperture.

This is equivalent to observing the optical wave on some kind of "complex-spherical" reference surface. Physically it corresponds, of course, to inserting into the \(\text{ABCD}\) system a self-canceling thin lens plus gaussian aperture pair, with one aperture having a positive and the other a negative gaussian transmission; and then again observing the optical beam at the reference plane between these two elements.

 

Transforming the ABCD Matrix

As one example of this technique, let us show how we can transform the elements of an arbitrary \(\text{ABCD}\) matrix in various ways by putting in such selfcanceling pairs before and after the \(\text{ABCD}\) system.

To do this, we visualize an optical element consisting of a thin lens plus possibly a gaussian aperture, so that this element has a total complex focal power, call it \(\tilde{p}\), given by \(\tilde{p}\equiv(1/f+j\lambda a_2/2\pi)\) in the notation of the previous sections.

We then introduce a self-canceling pair of such elements both before and after the arbitrary \(\text{ABCD}\) system, as in Figure 15.

Introducing these elements will obviously not change the physical performance of this \(\text{ABCD}\) system as a periodic focusing or roundtrip resonator matrix; and we can go through one period or one round trip from the reference planes between the self-canceling elements. 

Multiplying out the ray matrices will then give a transformed or primed ray matrix in the form

\[\tag{95}\left[\begin{array}&\hat{A}&\hat{B}\\\hat{C}&\hat{D}\end{array}\right]=\left[\begin{array}&1&0\\\tilde{p}&1\end{array}\right]\left[\begin{array}&A&B\\C&D\end{array}\right]\left[\begin{array}&1&0\\-\tilde{p}&1\end{array}\right]\\=\left[\begin{array}&=A-B\tilde{p}&B\\A\tilde{p}-B\tilde{p}^2+C-D\tilde{p}&D+B\tilde{p}\end{array}\right].\]

The ray matrix \(\text{ABCD}\) is then the ray matrix going from one reference plane between the self-canceling elements to the next such reference plane. 

 

Symmetric Matrix Form

Provided that the \(\text{B}\) element of the original \(\text{ABCD}\) matrix is not zero, we can convert the "hat" form of the \(\text{ABCD}\) matrix into a symmetrized form with \(\hat{A}=\hat{D}\) by choosing the complex curvature \(\tilde{p}\) to be 

\[\tag{96}\tilde{p}=\frac{A-D}{2B}.\]

With this choice of reference plane, the \(\text{ABCD}\) matrix will take on the symmetrized form 

\[\tag{97}\left[\begin{array}&\hat{A}&\hat{B}\\\hat{C}&\hat{D}\end{array}\right]=\left[\begin{array}&m&B\\(m^2-1)/B&m\end{array}\right]\quad\left(\begin{array}\text{symmetric}\\\text{matrix}\;\text{form}\end{array}\right).\]

Making the \(\hat{A}\) and \(\hat{D}\) elements equal will complete the square in the exponent of the Huygens' kernel of the previous section, and it will also simplify some of the resonator results of the following chapter.

Note also that for any choice of \(\tilde{p}\), the value of the \(B\) element of the matrix—which is essentially the "length" of the paraxial system—remains unchanged.

 

Canonical Matrix Form 

As another example, if we choose \(1/\tilde{p}\) to be equal to either of the eigen-\(\tilde{p}\)-values of the system, as discussed in the following chapter, namely, 

\[\tag{98}\tilde{p}=\frac{1}{B}\left[\frac{A-D}{2}\pm\sqrt{\left(\frac{A+D}{2}\right)^2-1}\right]=\frac{1}{\tilde{q}_a}\quad\text{or}\quad\frac{1}{\tilde{q}_b},\]

then the \(\text{ABCD}\) matrix will be cast into a canonical form in which \(\hat{C}\equiv0\) and \(\hat{A}\times\hat{D}=1\). The \(\hat{A}\) and \(\hat{D}\) elements will in fact become equal to the eigenvalues \(\lambda_a\) and \(\lambda_b\) of the system (next tutorial), and the \(\text{ABCD}\) matrix can be converted into the simplified form

\[\tag{99}\left[\begin{array}&\hat{A}&\hat{B}\\\hat{C}&\hat{D}\end{array}\right]=\left[\begin{array}&\lambda_a&B\\0&\lambda_b\end{array}\right]\quad\left(\begin{array}&\text{canonical}\\\text{matrix}\;\text{form}\end{array}\right),\]

or else the reverse of this with a and \(b\) interchanged.

Once it is in canonical form, the \(\text{ABCD}\) matrix can then easily be raised to the \(n\)-th power by writing 

\[\tag{100}\left[\begin{array}&\lambda_a&B\\0&\lambda_b\end{array}\right]=\left[\begin{array}&\lambda^n_a&B_n\\0&\lambda^n_b\end{array}\right],\]

where \(B_n\) is given by

\[\tag{101}\frac{B_n}{B}=\sum^{k=n}_{k=1}\lambda^{n-k}_a\lambda^{k-1}_b=\sum^{k=n}_{k=1}\lambda^{n+1-2k}_a=\sum^{k=n}_{k=1}\lambda^{n+1-2k}_b\]

These various kinds of transformation can prove useful in next tutorials.

 

 

 

Factorization of Arbitrary ABCD Matrices 

Other related questions in this same area involve how we might break down or factorize an arbitrary complex \(\text{ABCD}\) matrix into various subelements, for example, in order to synthesize an arbitrary system from a minimum number of subelements, or to find the \(n\)-th root of an arbitrary complex matrix.

These questions are closely associated with the basic mathematical properties of matrix theory; and paraxial optical matrices can provide interesting practical examples to illustrate these more general mathematical properties. Some preliminary work on this topic has been done. 

Casperson has shown, for example, as an illustration of the possibilities of these techniques, that the triple-cascaded combination of lenses and free spaces in Figure 16 has the interesting property that

\[\tag{102}\left\{\left[\begin{array}&1&0\\-3/L&1\end{array}\right]\times\left[\begin{array}&1&L\\0&1\end{array}\right]\right\}^3=I,\]

i.e., each individual section consisting of free space of length \(L\) and a thin lens of focal length \(L/3\) has an \(\text{ABCD}\) matrix corresponding to the cube root of the identity matrix. More work on this topic remains to be done.