Ray transfer matrix analysis

Ray transfer matrix analysis (also known as ABCD matrix analysis) is a mathematical form for performing ray tracing calculations in sufficiently simple problems which can be solved considering only paraxial rays. Each optical element (surface, interface, mirror, or beam travel) is described by a Template:Nowrap ray transfer matrix which operates on a vector describing an incoming light ray to calculate the outgoing ray. Multiplication of the successive matrices thus yields a concise ray transfer matrix describing the entire optical system. The same mathematics is also used in accelerator physics to track particles through the magnet installations of a particle accelerator, see electron optics.

This technique, as described below, is derived using the paraxial approximation, which requires that all ray directions (directions normal to the wavefronts) are at small angles Template:Mvar relative to the optical axis of the system, such that the approximation Template:Math remains valid. A small Template:Mvar further implies that the transverse extent of the ray bundles (Template:Mvar and Template:Mvar) is small compared to the length of the optical system (thus "paraxial"). Since a decent imaging system where this is Template:Em the case for all rays must still focus the paraxial rays correctly, this matrix method will properly describe the positions of focal planes and magnifications, however aberrations still need to be evaluated using full ray-tracing techniques.<ref>Extension of matrix methods to tracing (non-paraxial) meridional rays is described by Template:Harvp.</ref>

Matrix definition

File:RayTransferMatrixDefinitions.svg

In ray transfer (ABCD) matrix analysis, an optical element (here, a thick lens) gives a transformation between Template:Math at the input plane and Template:Math when the ray arrives at the output plane.

The ray tracing technique is based on two reference planes, called the input and output planes, each perpendicular to the optical axis of the system. At any point along the optical train an optical axis is defined corresponding to a central ray; that central ray is propagated to define the optical axis further in the optical train which need not be in the same physical direction (such as when bent by a prism or mirror). The transverse directions Template:Mvar and Template:Mvar (below we only consider the Template:Mvar direction) are then defined to be orthogonal to the optical axes applying. A light ray enters a component crossing its input plane at a distance Template:Math from the optical axis, traveling in a direction that makes an angle Template:Math with the optical axis. After propagation to the output plane that ray is found at a distance Template:Math from the optical axis and at an angle Template:Math with respect to it. Template:Math and Template:Math are the indices of refraction of the media in the input and output plane, respectively.

The ABCD matrix representing a component or system relates the output ray to the input according to <math display="block"> \begin{bmatrix}x_2 \\ \theta_2\end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix}x_1 \\ \theta_1\end{bmatrix}, </math> where the values of the 4 matrix elements are thus given by <math display="block">A = \left.\frac{x_2}{x_1} \right|_{\theta_1 = 0} \qquad B = \left.\frac{x_2}{\theta_1} \right|_{x_1 = 0},</math> and <math display="block">C = \left.\frac{\theta_2}{ x_1 } \right|_{\theta_1 = 0} \qquad D = \left.\frac{\theta_2}{\theta_1 } \right|_{x_1 = 0}.</math>

This relates the ray vectors at the input and output planes by the ray transfer matrix (Template:Dfn) Template:Math, which represents the optical component or system present between the two reference planes. A thermodynamics argument based on the blackbody radiation Template:Citation needed can be used to show that the determinant of a RTM is the ratio of the indices of refraction: <math display="block">\det(\mathbf{M}) = AD - BC = \frac{n_1}{n_2}. </math>

As a result, if the input and output planes are located within the same medium, or within two different media which happen to have identical indices of refraction, then the determinant of Template:Math is simply equal to 1.

A different convention for the ray vectors can be employed. Instead of using Template:Math, the second element of the ray vector is Template:Math,Template:Sfnp which is proportional not to the ray angle per se but to the transverse component of the wave vector. This alters the ABCD matrices given in the table below where refraction at an interface is involved.

The use of transfer matrices in this manner parallels the Template:Val matrices describing electronic two-port networks, particularly various so-called ABCD matrices which can similarly be multiplied to solve for cascaded systems.

Some examples

Free space example

As one example, if there is free space between the two planes, the ray transfer matrix is given by: <math display="block"> \mathbf{S} = \begin{bmatrix} 1 & d \\ 0 & 1 \end{bmatrix} , </math> where Template:Mvar is the separation distance (measured along the optical axis) between the two reference planes. The ray transfer equation thus becomes: <math display="block"> \begin{bmatrix} x_2 \\ \theta_2 \end{bmatrix} = \mathbf{S} \begin{bmatrix} x_1 \\ \theta_1\end{bmatrix} , </math> and this relates the parameters of the two rays as: <math display="block"> \begin{aligned} x_2 &= x_1 + d\theta_1 \\ \theta_2 &= \hphantom{x_1 + d}\theta_1 \end{aligned} </math>

Thin lens example

Another simple example is that of a thin lens. Its RTM is given by: <math display="block"> \mathbf{L} = \begin{bmatrix} 1 & 0 \\ -\frac{1}{f} & 1 \end{bmatrix} , </math> where Template:Mvar is the focal length of the lens. To describe combinations of optical components, ray transfer matrices may be multiplied together to obtain an overall RTM for the compound optical system. For the example of free space of length Template:Mvar followed by a lens of focal length Template:Mvar: <math display="block">\mathbf{L}\mathbf{S} = \begin{bmatrix} 1 & 0 \\ -\frac{1}{f} & 1\end{bmatrix} \begin{bmatrix} 1 & d \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & d \\ -\frac{1}{f} & 1-\frac{d}{f} \end{bmatrix} . </math>

Note that, since the multiplication of matrices is non-commutative, this is not the same RTM as that for a lens followed by free space: <math display="block"> \mathbf{SL} = \begin{bmatrix} 1 & d \\ 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ -\frac{1}{f} & 1 \end{bmatrix} = \begin{bmatrix} 1-\frac{d}{f} & d \\ -\frac{1}{f} & 1 \end{bmatrix} . </math>

Thus the matrices must be ordered appropriately, with the last matrix premultiplying the second last, and so on until the first matrix is premultiplied by the second. Other matrices can be constructed to represent interfaces with media of different refractive indices, reflection from mirrors, etc.

Eigenvalues

A ray transfer matrix can be regarded as a linear canonical transformation. According to the eigenvalues of the optical system, the system can be classified into several classes.Template:Sfnp Assume the ABCD matrix representing a system relates the output ray to the input according to <math display="block"> \begin{bmatrix}x_2 \\ \theta_2\end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix}x_1 \\ \theta_1\end{bmatrix} =\mathbf{T}\mathbf{v} .</math>

We compute the eigenvalues of the matrix <math> \mathbf{T} </math> that satisfy eigenequation <math display="block"> [\boldsymbol{T}-\lambda I] \mathbf{v} = \begin{bmatrix} A-\lambda & B \\ C & D-\lambda \end{bmatrix} \mathbf{v} = 0 ,</math> by calculating the determinant <math display="block"> \begin{vmatrix} A-\lambda & B \\ C & D-\lambda \end{vmatrix} = \lambda^2 - (A+D) \lambda + 1 = 0 .</math>

Let <math>m = \frac{(A+D)}{2}</math>, and we have eigenvalues <math>\lambda_{1}, \lambda_{2}=m \pm \sqrt{m^{2}-1}</math>.

According to the values of <math>\lambda_{1}</math> and <math>\lambda_{2}</math>, there are several possible cases. For example:

1. A pair of real eigenvalues: <math>r</math> and <math>r^{-1}</math>, where <math>r\neq1</math>. This case represents a magnifier <math> \begin{bmatrix} r & 0 \\ 0 & r^{-1} \end{bmatrix} </math>
2. <math>\lambda_{1}=\lambda_{2}=1</math> or <math>\lambda_{1}=\lambda_{2}=-1</math>. This case represents unity matrix (or with an additional coordinate reverter) <math> \begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix} </math>.
3. <math>\lambda_{1}, \lambda_{2}=\pm1</math>. This case occurs if but not only if the system is either a unity operator, a section of free space, or a lens
4. A pair of two unimodular, complex conjugated eigenvalues <math>e^{i\theta}</math> and <math>e^{-i\theta}</math>. This case is similar to a separable Fractional Fourier Transform.

Matrices for simple optical components

Element	Matrix	Remarks
Propagation in free space or in a medium of constant refractive index	<math>\begin{pmatrix} 1 & d\\ 0 & 1 \end{pmatrix} </math>	Template:Mvar = distance
Refraction at a flat interface	<math>\begin{pmatrix} 1 & 0 \\ 0 & \frac{n_1}{n_2} \end{pmatrix} </math>	Template:Math = initial refractive index Template:Math = final refractive index.
Refraction at a curved interface	<math>\begin{pmatrix} 1 & 0 \\ \frac{n_1-n_2}{R \cdot n_2} & \frac{n_1}{n_2} \end{pmatrix} </math>	Template:Mvar = radius of curvature, Template:Math for convex (center of curvature after interface) Template:Math = initial refractive index Template:Math = final refractive index.
Reflection from a flat mirror	<math> \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} </math>Template:Sfnp	Valid for flat mirrors oriented at any angle to the incoming beam. Both the ray and the optic axis are reflected equally, so there is no net change in slope or position.
Reflection from a curved mirror	<math> \begin{pmatrix} 1 & 0 \\ -\frac{2}{R_e} & 1 \end{pmatrix} </math>	<math>R_e = R\cos\theta</math> effective radius of curvature in tangential plane (horizontal direction) <math>R_e = R/\cos\theta</math> effective radius of curvature in the sagittal plane (vertical direction) Template:Mvar = radius of curvature, Template:Math for concave, valid in the paraxial approximation Template:Mvar is the mirror angle of incidence in the horizontal plane.
Thin lens	<math> \begin{pmatrix} 1 & 0 \\ -\frac{1}{f} & 1 \end{pmatrix} </math>	Template:Mvar = focal length of lens where Template:Math for convex/positive (converging) lens. Only valid if the focal length is much greater than the thickness of the lens.
Thick lens	<math>\begin{pmatrix} 1 & 0 \\ \frac{n_2-n_1}{R_2n_1} & \frac{n_2}{n_1} \end{pmatrix} \begin{pmatrix} 1 & t \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ \frac{n_1-n_2}{R_1n_2} & \frac{n_1}{n_2} \end{pmatrix}</math>	Template:Math = refractive index outside of the lens. Template:Math = refractive index of the lens itself (inside the lens). Template:Math = Radius of curvature of First surface. Template:Math = Radius of curvature of Second surface. Template:Mvar = center thickness of lens.
Single prism	<math> \begin{pmatrix} k & \frac{d}{nk} \\ 0 & \frac{1}{k} \end{pmatrix} </math>	<math>k = (\cos\psi / \cos\phi)</math> is the beam expansion factor, where Template:Mvar is the angle of incidence, Template:Mvar is the angle of refraction, Template:Mvar = prism path length, Template:Mvar = refractive index of the prism material. This matrix applies for orthogonal beam exit.<ref name=TLO>Template:Harvp</ref>
Multiple prism beam expander using Template:Mvar prisms	<math> \begin{pmatrix} M & B \\ 0 & \frac{1}{M} \end{pmatrix} </math>	Template:Mvar is the total beam magnification given by Template:Math, where Template:Mvar is defined in the previous entry and Template:Mvar is the total optical propagation distanceTemplate:Clarify of the multiple prism expander.<ref name=TLO />

Relation between geometrical ray optics and wave optics

The theory of Linear canonical transformation implies the relation between ray transfer matrix (geometrical optics) and wave optics.Template:Sfnp

Element	Matrix in geometrical optics	Operator in wave optics	Remarks
Scaling	<math>\begin{pmatrix} b^{-1} & 0\\ 0 & b\end{pmatrix} </math>	<math>\mathcal{V}[b] u(x)=u(b x)</math>
Quadratic phase factor	<math>\begin{pmatrix} 1 & 0\\ c & 1 \end{pmatrix} </math>	<math>Q[c]=\exp i \frac{k_{0}}{2} c x^{2}</math>	<math>k_0</math>: wave number
Fresnel free-space-propagation operator	<math>\begin{pmatrix} 1 & d\\ 0 & 1 \end{pmatrix} </math>	<math>\mathcal{R}[d]\left\{U\left(x_{1}\right)\right\}=\frac{1}{\sqrt{i \lambda d}} \int_{-\infty}^{\infty} U\left(x_{1}\right) e^{i \frac{k}{2 d}\left(x_{2}-x_{1}\right)^{2}} d x_1 </math>	<math>x_1 </math>: coordinate of the source <math>x_2 </math>: coordinate of the goal
Normalized Fourier-transform operator	<math>\begin{pmatrix} 0 & 1\\ -1 & 0 \end{pmatrix} </math>	<math>\mathcal{F}=\left(i \lambda_{0}\right)^{-1 / 2} \int_{-\infty}^{\infty} d x\left[\exp \left(i k_{0} p x\right)\right] \ldots </math>

Common decomposition

There exist infinite ways to decompose a ray transfer matrix <math> \mathbf{T} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} </math> into a concatenation of multiple transfer matrices. For example in the special case when <math>n_1 = n_2</math>:

1. <math> \begin{bmatrix} A & B \\ C & D \end{bmatrix}

= \left[\begin{array}{ll} 1 & 0 \\ D / B & 1 \end{array}\right]\left[\begin{array}{rr} B & 0 \\ 0 & 1 / B \end{array}\right]\left[\begin{array}{ll} 0 & 1 \\ -1 & 0 \end{array}\right]\left[\begin{array}{ll} 1 & 0 \\ A / B & 1 \end{array}\right] </math>.

1. <math> \begin{bmatrix} A & B \\ C & D \end{bmatrix}

= \left[\begin{array}{ll} 1 & 0 \\ C / A & 1 \end{array}\right]\left[\begin{array}{rr} A & 0 \\ 0 & A^{-1} \end{array}\right]\left[\begin{array}{ll} 1 & B / A \\ 0 & 1 \end{array}\right] </math>

1. <math> \begin{bmatrix} A & B \\ C & D \end{bmatrix}

= \left[\begin{array}{ll} 1 & A / C \\ 0 & 1 \end{array}\right]\left[\begin{array}{lr} -C^{-1} & 0 \\ 0 & -C \end{array}\right]\left[\begin{array}{ll} 0 & 1 \\ -1 & 0 \end{array}\right]\left[\begin{array}{ll} 1 & D / C \\ 0 & 1 \end{array}\right] </math>

1. <math> \begin{bmatrix} A & B \\ C & D \end{bmatrix}

= \left[\begin{array}{ll} 1 & B / D \\ 0 & 1 \end{array}\right]\left[\begin{array}{ll} D^{-1} & 0 \\ 0 & D \end{array}\right]\left[\begin{array}{ll} 1 & 0 \\ C / D & 1 \end{array}\right] </math>

Resonator stability

RTM analysis is particularly useful when modeling the behavior of light in optical resonators, such as those used in lasers. At its simplest, an optical resonator consists of two identical facing mirrors of 100% reflectivity and radius of curvature Template:Mvar, separated by some distance Template:Mvar. For the purposes of ray tracing, this is equivalent to a series of identical thin lenses of focal length Template:Math, each separated from the next by length Template:Mvar. This construction is known as a lens equivalent duct or lens equivalent waveguide. The Template:Abbr of each section of the waveguide is, as above, <math display="block">\mathbf{M} =\mathbf{L}\mathbf{S} = \begin{pmatrix} 1 & d \\ \frac{-1}{f} & 1-\frac{d}{f} \end{pmatrix} .</math>

Template:Abbr analysis can now be used to determine the stability of the waveguide (and equivalently, the resonator). That is, it can be determined under what conditions light traveling down the waveguide will be periodically refocused and stay within the waveguide. To do so, we can find all the "eigenrays" of the system: the input ray vector at each of the mentioned sections of the waveguide times a real or complex factor Template:Mvar is equal to the output one. This gives: <math display="block"> \mathbf{M} \begin{bmatrix}x_1 \\ \theta_1\end{bmatrix} = \begin{bmatrix}x_2 \\ \theta_2\end{bmatrix} = \lambda \begin{bmatrix} x_1 \\ \theta_1 \end{bmatrix} . </math> which is an eigenvalue equation: <math display="block"> \left[ \mathbf{M} - \lambda\mathbf{I} \right] \begin{bmatrix}x_1 \\ \theta_1\end{bmatrix} = 0 , </math> where <math display="inline">\mathbf{I} = \left[\begin{smallmatrix} 1&0 \\ 0&1 \end{smallmatrix}\right]</math> is the Template:Val identity matrix.

We proceed to calculate the eigenvalues of the transfer matrix: <math display="block">\det \left[ \mathbf{M} - \lambda\mathbf{I} \right] = 0 , </math> leading to the characteristic equation <math display="block"> \lambda^2 - \operatorname{tr}(\mathbf{M}) \lambda + \det( \mathbf{M}) = 0 , </math> where <math display="block"> \operatorname{tr} ( \mathbf{M} ) = A + D = 2 - \frac{d}{f} </math> is the trace of the Template:Abbr, and <math display="block">\det(\mathbf{M}) = AD - BC = 1 </math> is the determinant of the Template:Abbr. After one common substitution we have: <math display="block"> \lambda^2 - 2g \lambda + 1 = 0 , </math> where <math display="block"> g \overset{\mathrm{def}}{{}={}} \frac{ \operatorname{tr}(\mathbf{M}) }{ 2 } = 1 - \frac{ d }{ 2 f } </math> is the stability parameter. The eigenvalues are the solutions of the characteristic equation. From the quadratic formula we find <math display="block"> \lambda_{\pm} = g \pm \sqrt{g^2 - 1} . </math>

Now, consider a ray after Template:Mvar passes through the system: <math display="block"> \begin{bmatrix}x_N \\ \theta_N \end{bmatrix} = \lambda^N \begin{bmatrix}x_1 \\ \theta_1\end{bmatrix}. </math>

If the waveguide is stable, no ray should stray arbitrarily far from the main axis, that is, Template:Mvar must not grow without limit. Suppose Template:Nowrap Then both eigenvalues are real. Since Template:Nowrap one of them has to be bigger than 1 (in absolute value), which implies that the ray which corresponds to this eigenvector would not converge. Therefore, in a stable waveguide, Template:Nowrap and the eigenvalues can be represented by complex numbers: <math display="block"> \lambda_{\pm} = g \pm i \sqrt{1 - g^2} = \cos(\phi) \pm i \sin(\phi) = e^{\pm i \phi} , </math> with the substitution Template:Math.

For <math> g^2 < 1 </math> let <math> r_+ </math> and <math> r_- </math> be the eigenvectors with respect to the eigenvalues <math> \lambda_+ </math> and <math> \lambda_- </math> respectively, which span all the vector space because they are orthogonal, the latter due to Template:Nowrap The input vector can therefore be written as <math display="block"> c_+ r_+ + c_- r_- , </math> for some constants <math> c_+ </math> and Template:Nowrap

After Template:Mvar waveguide sectors, the output reads <math display="block"> \mathbf{M}^N (c_+ r_+ + c_- r_-) = \lambda_+^N c_+ r_+ + \lambda_-^N c_- r_- = e^{i N \phi} c_+ r_+ + e^{- i N \phi} c_- r_- , </math> which represents a periodic function.

Gaussian beams

The same matrices can also be used to calculate the evolution of Gaussian beamsTemplate:Sfnp propagating through optical components described by the same transmission matrices. If we have a Gaussian beam of wavelength Template:Nowrap radius of curvature Template:Mvar (positive for diverging, negative for converging), beam spot size Template:Mvar and refractive index Template:Mvar, it is possible to define a complex beam parameter Template:Mvar by:<ref name=Lei/> <math display="block"> \frac{1}{q} = \frac{1}{R} - \frac{i\lambda_0}{\pi n w^2} . </math>

(Template:Mvar, Template:Mvar, and Template:Mvar are functions of position.) If the beam axis is in the Template:Mvar direction, with waist at Template:Math and Rayleigh range Template:Mvar, this can be equivalently written as<ref name=Lei>{{#invoke:citation/CS1|citation |CitationClass=web }} especially Chapter 5Template:Self-published source</ref> <math display="block"> q = (z - z_0) + i z_R .</math>

This beam can be propagated through an optical system with a given ray transfer matrix by using the equationTemplate:Explain: <math display="block"> \begin{bmatrix} q_2 \\ 1 \end{bmatrix} = k \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix}q_1 \\ 1 \end{bmatrix} , </math> where Template:Mvar is a normalization constant chosen to keep the second component of the ray vector equal to Template:Math. Using matrix multiplication, this equation expands as <math display="block">\begin{aligned} q_2 &= k (A q_1 + B) \\ 1 &= k (C q_1 + D)\,.\end{aligned}</math>

Dividing the first equation by the second eliminates the normalization constant: <math display="block"> q_2 =\frac{Aq_1+B}{Cq_1+D} ,</math>

It is often convenient to express this last equation in reciprocal form: <math display="block"> \frac{ 1 }{ q_2 } = \frac{ C + D/q_1 }{ A + B/q_1 } . </math>

Example: Free space

Consider a beam traveling a distance Template:Mvar through free space, the ray transfer matrix is <math display="block">\begin{bmatrix} A & B \\ C & D \end{bmatrix} = \begin{bmatrix} 1 & d \\ 0 & 1 \end{bmatrix} .</math> and so <math display="block">q_2 = \frac{A q_1+B}{C q_1+D} = \frac{q_1+d}{1} = q_1+d</math> consistent with the expression above for ordinary Gaussian beam propagation, i.e. Template:Nowrap As the beam propagates, both the radius and waist change.

Example: Thin lens

Consider a beam traveling through a thin lens with focal length Template:Mvar. The ray transfer matrix is <math display="block">\begin{bmatrix}A&B\\C&D\end{bmatrix}=\begin{bmatrix}1&0\\-1/f&1\end{bmatrix}.</math> and so <math display="block">q_2 =\frac{Aq_1+B}{Cq_1+D} = \frac{q_1}{-\frac{q_1}{f}+1} </math> <math display="block">\frac{1}{q_2} = \frac{-\frac{q_1}{f} + 1}{q_1} = \frac{1}{q_1} - \frac{1}{f} .</math> Only the real part of Template:Math is affected: the wavefront curvature Template:Math is reduced by the power of the lens Template:Math, while the lateral beam size Template:Mvar remains unchanged upon exiting the thin lens.

Higher rank matrices

Methods using transfer matrices of higher dimensionality, that is Template:Val, Template:Val, and Template:Val, are also used in optical analysis.Template:Sfnmp In particular, Template:Val propagation matrices are used in the design and analysis of prism sequences for pulse compression in femtosecond lasers.<ref name=TLO />

See also

• Transfer-matrix method (optics)
• Linear canonical transformation

References

• Template:Cite journal
• Template:Cite book
• Template:Cite book
• Template:Cite book
• Template:Cite book
• Template:Cite journal
• Template:Cite conference
• Template:Cite journal
• Template:Cite book
• Template:Cite book

Further reading

• Template:Cite book

External links

• Thick lenses (Matrix methods)
• ABCD Matrices Tutorial Provides an example for a system matrix of an entire system.
• ABCD Calculator An interactive calculator to help solve ABCD matrices.