1. Introduction

Upon its invention in 1960, the laser immediately distinguished itself from other light sources by its high quality spatial coherence, which allowed it to propagate long distances in a pencil-like beam. Though this beam quality is just one of the laser’s many exemplary properties, it has proven particularly useful to science and technology, for example in distance measurements e.g., to the moon; high-resolution imaging, cutting, machining, welding, and 3D printing; the characterization of surface profiles such as the cornea; the construction of interferometric sensors e.g., for gravitational waves; and long-distance power delivery and communication. Nonetheless, it took another 30 years for the field to settle on a standard performance metric for laser beam quality. To this end, in 1990, Siegman drew attention to the $M^{2}$ parameter [1–3], the ratio of the product of the beam size in the near and far fields to the same product for a diffraction limited beam [2]. In $M^2$ literature, it is common to write M2 instead of $M^2$ in titles as we do in this work. A reliable way to measure $M^{2}$ was codified in 1999 by the ISO standard 11146 [4], which called for transverse spatial intensity-profiles to be recorded at ten distances along the propagation axis. The $M^{2}$ parameter is now routinely measured and reported in scientific research and included in the advertised specifications of commercial lasers.

Since 1999, beam measurement tools and optical modelling have continued to advance beyond intensity spatial-profiles. There are now established methods and even commercial devices for experimentally measuring the electric field spatial-profile of a beam. These include newer methods such as off-axis holography [5,6] and direct measurement [7], as well as time-tested devices such as Shack-Hartmann sensors [8–11] and shear plates [12,13]. In addition, with advances in computation, researchers and engineers, can now perform accurate electric field modeling of a proposed device design using finite-difference time-domain, transfer-matrix, or finite element analysis methods. With such commercial or custom optical simulation software one can predict the electric field profile produced by a nano-antenna, or transmitted by a multilayer thin-film stack, or reflected by a metalens, to give some examples. The $M^{2}$ parameter is often the goal of such simulations, particularly of devices in which new optical modes are excited e.g., large area optical fibers [14–16], materials with a non-linear optical response [17], and laser cavities [18,19]. In light of these advancements, intensity spatial-profiles may not always be the ideal choice for determining $M^{2}$.

Unlike an intensity profile, all the information about a laser beam’s quality is contained in its electric field spatial-profile in any single transverse plane. This makes the numerous intensity profiles and subsequent fitting of the ISO 11146 inconveniently circuitous, particularly for optical simulations. Moreover, these numerous measurements in the ISO standard are included, in a large part, to make the procedure robust to measurement errors (e.g., background) and uncertainties, which numerical modelling does not have. On top of this, for both simulations and experiment, the ISO procedure is additionally problematic since intensity profiles must be acquired at prescribed distances from the beam waist. This requires a pre-characterization of the beam before the numerous profiles are obtained, i.e., to find the approximate location of the beam waist and the Rayleigh range. The complexity of the ISO standard has motivated work over the last decades to simplify the experimental procedure to obtain $M^{2}$; for an overview see [20,21].

Early papers on $M^2$ [22–27] focused on characterizing $M^2$ through intensity measurements. Some papers [15,24,25,28] have developed theory based on the electric field distribution to derive formulae for $M^2$ for specific beam shapes with analytic forms. Recently, other works have developed experimental methods to obtain $M^2$ from the electric field, albeit indirectly. For example, [10] demonstrated a method to determine $M^2$ using a Shack-Hartmann sensor along with near and far-field intensity profiles, and [29,30] measured the modal decomposition of a beam in order to calculate $M^2$. Other works numerically replicated one of the experimental methods to obtain $M^2$. For example, [30–32] numerically propagated the electric field to obtain the intensity profile in different planes in order to replicate the ISO method. Reference [33] did the same for three planes in combination with the use of a cylindrical lens. Such a cylindrical lens has to be carefully selected every time that a different beam is used to ensure a fair sampling of the beam. In this work, using a position-angle phase-space picture, we motivate the fundamental meaning of $M^2$ and use that to introduce a method to find $M^2$ that only requires the electric field at one plane. This eliminates the need for numerical propagation or numerical projection onto a modal basis. We numerically validate such a method for two nontrivial test beams. Our method is particularly useful in optical simulations where the complex electric field profile is known.

The rest of this paper is organized as follows. In Section 2, we start by defining $M^{2}$ and outline the method in ISO 11146 [4] (which we call the “ISO standard” from hereon). We then introduce an angle-position phase-space for a beam’s state and relate $M^{2}$ to the state’s area, which we show is conserved under paraxial propagation (i.e., the range of angles $<<{1}\;{\textrm{rad}}$). We give explicit formulas to calculate $M^{2}$, the beam waist $w_{0}$ and its location $z_0$, the beam’s angular width $\Delta \theta$, and the Rayleigh range $z_{R}$ from the electric field profile at one plane. We call this the covariance method. In Section 3, we validate our covariance method by numerically comparing its predictions for $M^{2}$ to the ISO method for two non-trivial beam shapes. In Section 4, we discuss the advantages and drawbacks of our covariance method and the prospects for generalizing it to more complicated beams. While many of the concepts in this paper have been introduced elsewhere, they have not been connected in a simple way to allow for easy use. For those solely interested in applying our covariance method, in Section 2.5 we summarize our results and provide a detailed straightforward recipe for determining $M^{2}$ from a single electric field profile of the beam anywhere.

2. Theoretical description of laser beams

2.1 Definitions and conventions

For simplicity we will consider beam profiles in one transverse dimension only, Appendix 4. contains a discussion towards the generalization to two dimensions. The $z$ axis is the propagation direction while the $x$-axis is the transverse direction, along the beam’s spatial profile. As a beam propagates along $z$, generally its width in position $w\left (z\right )$ will change while its angular width $\Delta \theta$ will be constant. We use the standard width conventions $w\left (z\right )=\Delta x\left (z\right )=2\,\sigma _{x} \left (z\right )$ and $\Delta \theta =2\,\sigma _{\theta }$, where the $\sigma$ are the standard deviation of corresponding intensity distribution (i.e., in $x$ or $\theta$). We take $z=0$ as the position of the beam waist $w_{0}$, the beam’s minimum $w\left (z\right )$ over all $z$. It is important to note that, except for Gaussians, the $w_{0}$ width here differs from the definition of the beam waist width common in Gaussian optics, the $1/e^{2}$ intensity half-width. The Rayleigh range $z_{R}$ is defined as the distance from the waist (taken as $z=0$) at which the beam has increased in width to $w\left (z_{R}\right )=\sqrt {2}w_{0}$.

2.2 Definition of the $M^{2}$ parameter

The beam or mode quality $M^{2}$, also known as the beam propagation factor, is defined in terms of the beam waist $w_{0}$ and angular width $\Delta \theta$. The product of these two widths divided by the same product for a Gaussian beam is the definition of $M^{2}$. Since for any Gaussian this product equals $\lambda /\pi$ [3] we have,

Whereas the product $w\left (z\right )\Delta \theta$ will change as the beam propagates, since $w_{0}$ and $\Delta \theta$ are independent of $z$ their product will not change. They are characteristics of the beam as a whole, as is the beam quality $M^{2}$.

As we will show, the value of $M^{2}$ is bounded by the Heisenberg uncertainty principle. As can be seen in Fig. 1(a), the relation between the widths in angle ($\Delta \theta$) and transverse wavevector ($\Delta k_{x}$) can be used to express the transverse momentum ($\Delta p_{x}$) of a photon as follows

 

Fig. 1. The evolution of paraxial beams under propagation. (a) Sketch showing the beam width along $z$ for a Gaussian beam (red, $M^2=1$) and combination of first four Hermite-Gauss modes (blue, $M^2=1.5$) with $w_0=$ 500 µm, $\lambda =$ 400 nm. The $z$-values are in units of the Gaussian’s Rayleigh distance $z_R =$ 1.96 m. (b) Normalized intensity profiles of the same beams as in (a) at three $z$-positions. (c) The phase-space ellipse representing the beam state at a general $z$-position. (d) The $x$-shear resulting from paraxial propagation of the beam. Blue cross indicates origin.

Download Full Size | PDF

Although Eq. (1) for $M^{2}$ seems quite simple, determining its value is deceptively difficult. The angular width $\Delta \theta$ can be found from the beam width far from the waist, i.e., in the far-field. In contrast, to measure the beam waist $w_{0}$ requires that one knows its $z$-location. This is typically not known a priori, and it is seldomly known with good accuracy. The ISO standard avoids this accuracy issue by recommending measuring the beam width $w\left (z\right )$ in five different $z$-planes on each side of the beam waist, half of which should lie within one Rayleigh range $z_{R}$ and the rest at least two Rayleigh ranges away from it. While one is free to choose the precise locations of the planes, to do so one still must approximately know $z_{R}$ and the location of the beam waist. A rough pre-characterization can take the place of this a priori information, but this adds complexity to the procedure. The basic definition (Eq. (1)) is hardly applicable in either experiment or simulation.

2.3 Connection between phase-space and the electric field

We now introduce a phase-space picture of the beam in order to justify a simpler method to determine $M^{2}$. Appropriately, the 2D phase-space is made up of position $x$ on one axis and angle $\theta$ on the perpendicular axis (though the latter could equally well be $k_{x}$ or $p_{x}$, as explained in the last section). In this phase-space, the beam profile at any $z$-position can be used to define a tilted ellipse (see Fig. 1(c)) which contains information about the extent of the beam’s state in phase-space. The length of the projection of the ellipse on the $x$-axis is the full-width intensity standard deviation, $w\left (z\right )=\Delta x$. Similarly, the projection on the $\theta$-axis is the angular beam width, $\Delta \theta$.

When the beam has reached its waist (taken as $z = 0$), the ellipse is aligned with the phase-space axes and its full-width in $x$ is minimized, equalling $w_{0}$. Consequently at the waist, the area $\mathcal {A}$ of the ellipse is simple to express, being proportional to its width times its height, $\mathcal {A}=\tfrac {1}{4}\pi w_{0}\Delta \theta$. In turn, using Eq. (1), we find that the beam quality is proportional to the ellipse area at the waist,

As the beam propagates away from the waist location, correlations between position and angle form, the ellipse becomes tilted and its principal axes change in length. These correlations follow from simple geometry; waves traveling at angle $\theta$ will increase in position as $x=\theta z$ in the small angle approximation, creating a correlation between $x$ and $\theta$. This shears the ellipse along the $x$-direction in phase-space [34] as illustrated in Fig. 1(d). (Conversely, momentum conservation in free space ensures that each constituent wave’s angle $\theta$ is constant.) Equivalently, free space propagation is given by the transfer matrix in Eq. (13), which describes a shear transformation of a beam in phase-space along the $x$-direction. However, since the ellipse is now tilted, the ellipse projections (divided by two) on the $x$ and $\theta$ axes are distinct from the width and height of the ellipse itself i.e., along its semi-minor and semi-major axes. Crucially, in geometry, a shear transformation always preserves area [35,36] and, thus, $M^{2}$. This implies that Eq. (4) is true for all $z$, not just at the waist. This finding has been obtained by other means in [37–39]. Thus, the problem of determining $M^{2}$ from information in a single $z$-plane reduces to finding the area of the corresponding tilted ellipse.

2.4 Covariance matrix

The covariance matrix is a powerful mathematical tool, used across physics and math, for example for the description of multimode quantum states of light (see Ref. [40] and the references therein). The covariance matrix is used to describe the beam at given $z$ and, thus, is closely related to “beam matrix” that is common in $M^2$ literature (e.g., in [26,41] and Part 3 of the ISO standard [42]), which is defined in terms of the Wigner function. Formally, the covariance matrix is simply the Weyl transform of the beam matrix. The transform ensures that it is directly computed directly from the electric field distribution $E\left (x;z\right )$ rather than indirectly, i.e., from expectation values on the Wigner function. Leaving behind that formal and somewhat opaque definition, in the following we physically motivate the covariance matrix.

We gather the parameters of a tilted ellipse centered at the origin (i.e., the ellipse equation, $ax^{2}-2bx\theta +c\theta ^{2}=ac-b^{2}$) into a symmetric matrix [43], a method from geometry known as the matrix of quadratic form. This matrix sets the ellipse aspect ratio, orientation of its axis, and its size. In the context of beam widths, the matrix of quadratic form is the following covariance matrix,

The advantage of this matrix formulation is that the determinant of the matrix $Q$ is proportional to the ellipse area, regardless of any tilt. The ellipse area is now simple to find by $\mathcal {A}=\pi \sqrt {\det Q}$. Combining this with Eq. (4) we arrive at our central idea, the determinant of the covariance matrix is directly related to $M^{2}$ via

This relation is valid at arbitrary $z$-planes as long as the (paraxial) beam propagates only through first-order optical systems (e.g., spherical mirrors and lenses) or freely in space [41].

In summary, one can evaluate $M^{2}$ using the elements of the covariance matrix Eq. (5). These expectation values use the complex field $E\left (x;z\right )$ (in analogous way to the quantum wavefunction). That is, the expectation value of a general operator $v$ acting on $E$ is defined by $\langle { v }\rangle =\tfrac {1}{n}\int {E^{*}{v}Edx}$, where $E^{*}$ is the complex conjugate of the electric field, and $n=\int {|E\left (x;z\right )|^{2}dx}$ is a normalization factor. The angle operator is $\theta =-\frac {i}{k}\frac {\partial }{\partial x}$, while $x$ is simply a multiplication by $x$. Note, the action of these two operators changes if their order changes, which may motivate why each off-diagonal in the covariance matrix incorporates both orderings. In general, the phase-space ellipse might not be centered at the origin, i.e., $\langle { x }\rangle \neq 0$ and/or $\langle {\theta }\rangle \neq 0$. We can account for this by substituting $x$ and $\theta$ in the covariance matrix (Eq. (5)) with $x- \langle {x}\rangle$ and $\theta -\langle {\theta }\rangle$, respectively. We give explicit formulae for elements $a,b,$ and $c$ for such offset beams in Appendix E. For example, the beam width can be evaluated as

In Appendix A, we analytically model the propagation of the $Q$ matrix. We bring these definitions together to provide an explicit formula for $M^{2}$ and other key beam parameters in Section 2.5.

2.5 Recipe to compute $M^{2}$ and other beam parameters from the electric field

The goal of this section is to be self-contained and provide a straightforward set of steps to calculate $M^{2}$ and other key beam parameters directly from the transverse profile of the beam’s scalar electric field $E\equiv E\left (x;z\right )$ that one has obtained in any plane $z$ along the propagation direction. To do so, we use Eq. (6) and the covariance matrix Eq. (5), which was found under the assumption that $\langle { x }\rangle = \langle { \theta }\rangle = 0$. This assumption is valid for most optical simulations since the incoming beam is ideal and usually travelling centered on the $x$-axis and because the optical device usually is symmetric about $x=0$. We give explicit formulae for elements $a,b,$ and $c$ for the more general case, offset beams, in Appendix E. Reviewing, the beam must be paraxial with angles to $z$-axis much less than 1 rad; as usual, the width convention for the beam waist $w_{0}$ and angular spread $\Delta \theta$ is twice the intensity standard-deviation (e.g., $w_{0}=2\sigma _{x})$; and the beam quality is

To evaluate this formula, one uses the covariance matrix elements, all of which are real-valued:

From the $a,b,$ and $c$ matrix elements, given above we can also find other key beam parameters, such as the beam’s angular width $\Delta \theta$, Rayleigh range $z_{R}$, beam waist $w_{0}$, and waist’s location relative to the plane of the measured field, $\Delta z$. In Eqs. (15)–(19) in Appendix A, we show that

Thus, like the beam quality, these parameters can be determined from the electric field profile at one plane.

The formulas in this section are the main results of this paper. The next section numerically verifies that they are correct, and also provides further details for the numerical evaluation of the electric field integrals in the $a,b,$ and $c$ parameters.

3. Simulation and results

In this section, we compare the $M^{2}$ found from our covariance method to the ISO standard via a numerical simulation of propagation of two beams with nontrivial electric field profiles. We calculate $E\left (x\right )$ at an initial plane ($z=0$) using the analytic form of the test beam $E\left (x\right )$ and then numerically propagate it to find $E\left (x;z\right )$ and the corresponding intensity profile $I\left (x;z\right )=\left |E\left (x;z\right )\right |^{2}$ at the ten requisite planes described below. This propagation is detailed in Appendix B and contains no paraxial approximation.

For the implementation of the ISO protocol, we calculate the standard deviation of the transverse intensity profile $I\left (x;z\right )$ to find the beam width $w\left (z\right )$ using Eq. (7). We do so at five planes within half a Rayleigh distance on either side of the beam waist and at five planes in the range of four to five Rayleigh distances. These widths are fit to Eq. (18) to find $M^{2}$.

At each of these planes, we calculate the entries, $a$, $b$, and $c$ of the covariance matrix $Q$ from $E\left (x;z\right )$ by Eq. (9) and use them to calculate $M^2$ using Eq. (6). We compare these ten values of $M^{2}$ from our covariance method, which should all be identical, to the single value from the ISO protocol. In Appendix C we give some details of the numerical implementation.

3.1 Numerical comparison of the ISO and covariance methods

For the following test beams, $\lambda =$ 400 nm, we use position range of length $L=$ 200 mm and a position grid density $\delta x=\lambda /2$. The ten planes range from $z=-7$ to 7 m.

For our first test beam, we use a linear combination of the first four Hermite-Gauss (HG) modes with complex coefficients. The electric field distribution at the initial plane for the m-th HG mode is taken to be

The test beam’s electric field is given by

For this first test beam we plot the elements of the offset $Q$ in Fig. 2(a). The element $a=w^{2}\left (z\right )/4$ (blue points) lies on the theoretical beam waist curve (Eq. (18), blue curve) and the other elements (yellow and black points) lie on the curves given by the matrix elements of the propagated $Q$ matrix (Eq. (14), yellow and black curves). This agreement confirms that our numerical propagation of $E\left (x\right )$ is accurate and our calculation of the covariance matrix $Q$ is correct.

 

Fig. 2. Analytic and numerical propagation of test beams. For a beam that is the superposition of the first four Hermite-Gauss modes given in the main text ($w_0=$ 642 µm, $\lambda =$ 400 nm $M^2=2.466$ given by Eq. (28) ): (a) points are the scaled entries of the covariance matrix for the beam numerically propagated (by Eq. (20)) to each $z$ distance and curves are the analytically propagated elements (by Eq. (14)) and (b) is the corresponding intensity profiles at three $z$-locations. For a paraxial supergaussian beam ($w_0=$ 56 µm, order = 50, $\lambda =$ 400 nm, nominal $M^2=4.0466$): (c) points are the beam widths $w_0$ calculated from the numerically propagated beam and the curve is the analytical width from Eq. (18).

Download Full Size | PDF

With the correctness of our numerical calculations established, we now compare the ISO standard method to our covariance method for an offset beam to find $M^{2}$. The $M^2$ values found from our covariance method (Eqs. (8), (22)) at the ten planes have average value $M^{2}=2.470$ with a standard deviation of $0.006$. Thus, as expected, our covariance method results in the same value of $M^{2}$ and it is independent of the plane $z$ of the $E\left (x;z\right )$ used to calculate it. Moreover, it is in good agreement with the theoretical prediction of $M^{2}=2.466$. A fit to the waists in the ten planes, the ISO method, gives $M^{2}= 2.467$ with a fit error of $2.6 \times 10^{-10}$, again in good agreement. This suggests that using solely a single plane and the covariance method will agree with the ISO method up to the second decimal place of $M^2$.

Our second test beam is meant to approximate the top hat intensity-distribution that would be transmitted by a slit. In this way, it differs greatly from a basic Gaussian beam or even combinations of low-order HG modes. Moreover, an exact top hat has an angular intensity distribution that follows a sinc-squared function and thus has a $\Delta \theta$ that diverges and is not well-defined, which makes it a particularly challenging test case. A one-dimensional supergaussian, $E\left (x\right )=\exp [-\left (|x|/w_s\right )^N]$, of high order $N$ approximates a top hat function of full-width $2w_s$. Using Eqs. (6) and (9) we theoretically found $M^2=N\sqrt {\Gamma \left (3 / N\right )\Gamma \left (2-1/N\right )}/\Gamma \left (1 / N\right )\approx \sqrt {N/3}$, where the approximation is valid for large $N$ and $\Gamma$ is the usual gamma function. (Note, this formula differs from formulae in the literature, which are for the “cylindrical $M^2$”, e.g., [45].) The test beam is a supergaussian of order $N=50$ with equivalent slit width of $2w_s = {0.1}\;{\textrm{mm}}$ and a theoretical $M^2=4.0466$. Diffraction from such a slit has zero-to-zero angular width of $1.6\times 10^{-2}\;\rm {rad}$, putting it well within the paraxial regime.

For this second test beam, we repeat the comparison of methods that was conducted on the first beam. The beam widths in the ten planes for the supergaussian are plotted in Fig. 2(c). The ten planes of our covariance method give an average value of $M^{2}$= 4.0462 with a standard deviation of $7.0 \times 10^{-5}$. The ISO method (based on a fit to the beam waists in multiple planes) gives $M^{2}= 4.0506$ with a fit error of $1.4 \times 10^{-9}$. The ISO method’s value only agrees with theory and covariance method up to two decimal places after rounding. In contrast, our covariance method agrees with the analytical theoretical value to four decimal places, suggesting it is a more reliable method for this particularly challenging test case.

The success of our covariance method for finding $M^{2}$ with these two paraxial test beams validates its correctness and demonstrates that it stands in good agreement with the ISO protocol even for the extreme case of highly non-Gaussian beams.

4. Discussion and conclusion

This paper has focused on the simplest case, a paraxial coherent one-dimensional profile. In Appendix D, we similarly test non-paraxial beams. We show that both our covariance method and the ISO method fail for non-paraxial beams, as expected. A one-dimensional description applies to simple coherent beams whose two-dimensional profile can be written as a product, $E\left (x\right )E\left (y\right )$, e.g., two dimensional HG modes. In Appendix E, we discuss the generalization of our covariance method to a wider variety of paraxial beams, including beams with general two-dimensional profiles and incoherent beams.

In summary, we have demonstrated how one can compute the beam quality factor $M^{2}$ as well as angular width, waist width, and waist location for a beam from its electric field distribution in a single arbitrary plane. The conventional method, prescribed in the ISO standard [4], relies on finding the beam width in ten prescribed propagation planes. Since each width calculation requires three integrals of the intensity distribution, the ISO method requires thirty integrals in total. In comparison, our covariance method requires only six integrals of the electric field distribution. Moreover, it eliminates the need to fit a function, a key step in the ISO method. Further, our covariance method eliminates the need for an initial estimate of the Rayleigh range and waist location.

Since there are now many tools to measure the electric field profile of a beam, our covariance method is amenable to use in laboratories and could be incorporated in optical test and measurement products. That said, we envision our method will be particularly useful for analyzing the electric field output of optical simulations and could be directly built into common simulation software. In this way, we expect our covariance method to calculate beam quality will streamline optical research and development.

Appendix A: Paraxial propagation and beam parameters

In this section we treat paraxial propagation within the matrix formalism that we introduced in Section 2.4, which will allow us to determine other key beam parameters and connect with the ISO standard procedure. As discussed in Section 2.4, as the beam propagates in space, the covariance matrix $Q$ evolves. The covariance matrix can be propagated via the ray transfer matrix

(13)$$P=\begin{bmatrix}1 & z'\\ 0 & 1 \end{bmatrix}\text{ \,\,and the propagation law\,\, }Q\left(z+z'\right)=P\,Q\left(z\right)\,P^{T},$$

where $T$ indicates transpose. This simple formula is valid in both the Fresnel and Frauhnhofer diffraction regimes. The fact that $P$ has unit determinant ensures that the ellipse area is conserved under paraxial propagation, which justifies Eq. (6) for all $z$-planes.

We use Eq. (13) to propagate an arbitrary $Q\left (z\right )$ by a distance $z'$,

(14)$$Q\left(z+z'\right)\equiv Q'=\begin{bmatrix}a' & b'\\ b' & c' \end{bmatrix}=\begin{bmatrix}a+bz'+\left(cz'+b\right)z' & b+cz'\\ b+cz' & c \end{bmatrix}.$$

From this matrix we can retrieve the common beam parameters. First, we note that $c$ is unchanged so, trivially,

(15)$$\Delta\theta=2\sqrt{c}.$$

Next, we find the distance to the waist from $z$, the plane $Q$ was determined at, by setting $Q'$ to be diagonal, $b'=0=b+cz'$. Taking the position of the waist ($Q'$) to be $z_{0}$ so that $z'=z-z_{0}\equiv -\Delta z$, we then find

(16)$$\Delta z=\frac{b}{c}=\frac{\langle{ x\theta+\theta x}\rangle}{2\langle{\theta^{2}}\rangle},$$

where a positive $\Delta z$ means the waist is further along the beam propagation direction from the plane that $E\left (x;z\right )$ (and $Q$) was determined in. With $z'=-b/c$ substituted in $Q'$ (Eq. (14)), we find the beam waist,

(17)$$w_{0}=2\sqrt{a'}=2\sqrt{a-b^{2}/c}.$$

We now instead set $Q\left (z_{0}\right )$ to be at the waist and, thus, diagonal ($b=0$), and using $M^{2}$ from Eq. (1), the $a'$ element of $Q'$ gives the beam width:

(18)$$4a'=w\left(z\right)=w_{0}^{2}+\left[\frac{M^{2}\lambda}{\pi w_{0}}\right]^{2}\left(z-z_{0}\right)^{2},$$

where $z'=z-z_{0}$. This gives the Rayleigh range,

(19)$$z_{R}=\frac{\pi w_{0}^{2}}{M^{2}\lambda}=\frac{\lambda}{4\pi c}.$$

In other words, the propagation of an arbitrary paraxial beam is same as the one of a Gaussian beam magnified by $M^{2}$.

Equation (18) above is precisely the relation that the ISO standard uses to model the evolution of the beam width under spatial propagation. One first finds the intensity distributions in ten planes, then finds their widths, and then fits them to the hyperbolic function in Eq. (18). This fit yields values for $M^{2}$, the beam waist $w_{0}$ and its location $z_{0}$ [4]. A big disadvantage of this approach is that we already need to have rough values for the position of the beam waist as well as the $M^{2}$ parameter in order to compute the effective Rayleigh distance $z_{R}$ and be able to follow the steps given in the protocol.

Appendix B: Propagation simulation

For the spatial propagation of the two beams to the ten planes, we draw on the angular spectrum method [46]. The Fourier transform of $E\left (x;0\right )$ decomposes the input field into plane waves, which then can be propagated by distance $z$ via:

(20)$$E\left(x;z\right)=\mathcal{F}^{{-}1}\left[e^{ik_{z}z}\mathcal{F}\left[E\left(x;0\right)\right]\right],$$

where $k_{z}=\sqrt {k^{2}-k_{x}^{2}}$. This method is accurate for all angles for a scalar field. That is, it is valid beyond the paraxial approximation. For the implementation, we used the native Fast Fourier Transform (FFT) in Python.

Appendix C: Numerical details

We now briefly discuss some details of the numerical implementation for our covariance method. We use Python to evaluate Eq. (9). Derivatives were calculated with the Numpy gradient function, which uses the central difference method. Scipy integrate (Simpson’s rule) was used for integrals as we show in Code 1 (Ref. [47]).

The accuracy of our covariance method is mainly set by the $x$-position grid of the electric field profile $E\left (x\right )$. This directly sets the range and grid density in $x$ and indirectly sets them in $\theta$. In particular, since the two directions are related by a Fourier transform, the grid spacing $\delta \theta$ in the $\theta$-direction is set by the range $L$ in the $x$-direction, $\delta \theta \propto 1/\left (kL\right )$. Additionally, the range $\Theta$ in the $\theta$-direction is set by the grid spacing $\delta x$ in the $x$-direction, $\Theta \propto 1/\left (k\delta x\right )$.

The two goals are to have sufficient grid density to resolve the state and sufficient range to span the state in phase-space. That is, in both the $x$ and $\theta$ directions, we aim for $r$ points across the minimum width of the state and a range of $s$ times the maximum width of the state. The minimum and maximum widths in $x$ are $w_0$ and $w\left (z\right )$, which set $\delta x<w_{0}/r$ and $L>w\left (z\right )s$. The beam width is always $\Delta \theta$ in $\theta$, so $\delta \theta = \Delta \theta /r$ and $\Theta =s\Delta \theta$. By the Fourier relations from above, these set $L>\Delta \theta /\left (rk\right )$ and $\delta x<\left (s\Delta \theta \right )/k$, respectively. These four inequalities ensure sufficient resolution and range.

The inequalities determine the grid for a given beam state based on the chosen range and resolution parameters. For the latter, we suggest, $s=r=10$. Since $w_{0}$ and $\Delta \theta$ are not known a priori, one should start with a trial grid $L>10$ $w\left (z\right )$ and then calculate $w_{0}$ and $\Delta \theta$ from Eqs. (9,10). At this point, one would potentially need to revise the grid. Alternately, if $E\left (x\right )$ is from an optical simulation and computational power is not a constraint, setting $\delta x<\lambda /10$ will certainly be sufficient since it is five times better than the Nyquist limit. Satisfying the above conditions will ensure an accurate value of $M^{2}$.

Appendix D: Non-paraxial beams

Both the ISO protocol and our covariance method rely on the paraxial approximation. More problematically, the derivation of the $M^{2}$ from the uncertainty principle relies on the small-angle approximation, so it is not clear that the definition is completely general. In this section, we briefly examine the validity of both the ISO and covariance methods in the non-paraxial regime. For the following, $\lambda = {400}\;{\textrm{nm}}$ and we use position range of length $L={200}\;{\textrm{mm}}$ and a position grid density $\delta x=\lambda /2$.

To examine the non-paraxial regime, we decrease the effective slit full-width of the order $N=50$ supergaussian to ${0.5}$µm, which increases the zero-to-zero angular spread to ${1.6}{\textrm{rad}}$ putting it outside the paraxial approximation. Figure 3(b) should give an impression of the non-paraxiality of the supergaussian beam. The beam quality factor for a supergaussian is solely a function of order, so the theoretical beam quality should be the same as the $N=50$ supergaussian from Section 3.1, $M^{2} = 4.0466$. Applying the ISO method, the width measurements can still be fitted to a hyperbolic function (Fig. 3(a), top curve). However, the resulting $M^{2}$ factor more than doubles to a value of 9.28. While the covariance matrix gives the correct result for $M^2$ in the $z=0$-plane, the values quickly diverge as soon as propagation starts. At $z=$ 2.9 µm, the covariance method gives, $M^2=43.5$, which is off by a factor of five at least. The limitation to paraxial beams might be lifted when defining the beam matrix through a non-paraxial version of the Wigner function [34], but a thorough investigation would be in order here.

 

Fig. 3. (a) Width evolution of a non-paraxial supergaussian ($w_0=$ 0.56 µm, order = 50, $\lambda =$ 400 nm). Measured widths (indicated by points) do not lie on predicted curve for $M^2=4.05$. (b) Intensity plot of the same supergausssian. Red lines indicate the full angular spread of the beam $2\Delta \theta = {0.91}\;{\textrm{rad}}$.

Download Full Size | PDF

Appendix E: General two-dimensional beams

We now discuss some potential generalizations of our covariance method to a wider class of paraxial optical beams. We leave the details for future work.

Offset beams

In some optical simulations and all laboratory measurements one cannot assume that the beam is perfectly centered on the beam axis, i.e., $\langle {x}\rangle \neq 0$ and $\langle {\theta }\rangle \neq 0$. In this case, the covariance matrix takes the more general form:

(21)$$Q\left(z\right)=\begin{bmatrix} \langle{x^2}\rangle - \langle{x}^{2}\rangle & \frac{1}{2} \langle{x\theta+\theta x }\rangle -\langle{ x }\rangle \langle{ \theta}\rangle\\ \frac{1}{2}\langle{ x\theta+\theta x}\rangle -\langle{ x }\rangle \langle{ \theta}\rangle & \langle{ \theta^2 }\rangle-\langle{ \theta}^{2}\rangle \end{bmatrix}=\begin{bmatrix}a & b\\ b & c \end{bmatrix}.$$

The corresponding generalized formulas for the matrix elements are

(22)$$\begin{aligned} \\& a = \frac{1}{n}\int {{x^2}} {\left| E \right|^2}dx{\left( {\frac{1}{n}\int {x{{\left| E \right|}^2}dx} } \right)^2},\\& b = \frac{{i\lambda }}{{2\pi }}\left[ {\left( {\frac{1}{n}\int {x{{\left| E \right|}^2}} dx} \right)\left( {\frac{1}{n}\int {{E^ * }} \frac{{\partial E}}{{\partial x}}dx} \right) - \left( {\frac{1}{2} + \frac{1}{n}\int {x{E^ * }} \frac{{\partial E}}{{\partial x}}dx} \right)} \right],\\& c = {\left[ {\frac{\lambda }{{2\pi }}} \right]^2}\left[ {{{\left( {\frac{1}{n}\int {{E^ * }} \frac{{\partial E}}{{\partial x}}dx} \right)}^2} - \frac{1}{n}\int {{E^ * }} \frac{{{\partial ^2}E}}{{\partial {x^2}}}dx} \right],\\& n = {\int {\left| E \right|} ^2}dx \end{aligned}$$

As before, these can be used to find $M^2$ and other beam parameters using Eqs. (8,10).

Simple astigmatic beams

So far, we have presented our method in one transverse dimension $x$ for clarity. This one-dimensional analysis is straightforward to generalize to two-dimensional simple astigmatic beams; those that have circular and elliptical cross sections, respectively. To apply our covariance method to simple astigmatic beams, one first would find the principal axes of the elliptical transverse profile and denote these by $x$ and $y$. Along these directions, the total field is separable, $E\left (x,y;z\right )=E_{x}\left (x;z\right )E_{y}\left (y;z\right )$, and each transverse direction can be separately analyzed by our covariance method. The result would be two values, $M_{x}^{2}$ and $M_{y}^{2}$, from which one can define an effective beam quality $M_{\text {eff}}^{2}=\sqrt {M_{x}^{2}M_{y}^{2}}$. See [4] for details.

General astigmatic beams

We now propose a route to adapt our covariance method to find the beam quality of more general two-dimensional paraxial fields, $E\left (x,y;z\right )\neq E_{x}\left (x;z\right )E_{y}\left (y;z\right )$. Unlike simple astigmatic beams, the principal axes of their transverse cross sections can rotate while the beam is propagating [48]. Part 2 of the ISO standard [49] derives an effective beam quality $M_{\text {eff}}^{2}$ in terms of moments of the Wigner function in the four-dimensional phase-space. Essentially, the state is a generalized ellipse in this four-dimensional space and $M_{\text {eff}}^{2}$ is the area of that state. Analogous to the generalized $P$ matrix in Part 2 of the ISO standard, a generalized $Q_{xy}\left (z\right )$ matrix would be a 4x4 matrix:

(23)$$\ Q_{xy}\left(z\right)\equiv\left[\begin{array}{cc} Q_{x} & D_{xy}\\ D_{xy} & Q_{y} \end{array}\right].$$

Here, the block diagonals are the standard 2x2 $Q$ matrices Eq. (5): $Q_{x}\equiv Q$ and $Q_{y}$ is analogous. The $x$ and $y$ are arbitrary orthogonal transverse directions unconnected to the beam state. The two off-diagonal blocks are identical and given by

(24)$$D_{xy}\left(z\right)=\left[\begin{array}{cc} \langle{ xy }\rangle & \langle{ x\theta_{y} }\rangle\\ \langle{ y\theta_{x}}\rangle & \langle{ \theta_{x}\theta_{y}\rangle } \end{array}\right].$$

Note, the $D$ blocks involve one variable from each direction and thus do not require the two orderings (e.g., $\langle { x\theta +\theta x}\rangle$) that are essential for the 2x2 $Q$ matrix. With this generalized covariance matrix, $M_{\text {eff}}^{2}=\frac {4\pi }{\lambda }\sqrt [4]{\det Q_{xy}\left (z\right )}$. Unlike the $P$ matrix in the ISO procedure, $Q_{xy}\left (z\right )$ is found directly from the electric field $E\left (x,y;z\right )$.

Incoherent beams

We now consider generalizing our covariance method to incoherent beams. We tested our covariance method with beams that are coherent superpositions of HG modes. If the beam is an incoherent mixture of beams, one would need to use a density matrix formalism, or equivalently, a coherence matrix to represent the beam in single $z$ plane, $\rho \left (x,x'\right )\equiv \overline {E\left (x\right )E^{*}\left (x'\right )}$. Here, the bar indicates an ensemble or time average, as is standard in statistical optics. We leave the details for future work, but in short the expectation values in $Q_{xy}\left (z\right )$ would then be $\langle { v}\rangle =\mathrm {Tr}\left [v\rho \right ]$ and the beam quality $M^{2}$ would be found from $Q$ as before, from Eq. (6). This approach would complete the presented covariance method. However previous work have shown how to obtain $M^2$ for an incoherent mixture of beams [1,25], we refer the reader to such references for further details.

Appendix F: Derivation of $M^2$ for a general superposition of HG beams

We now derive a general formula for $M^2$ for a beam $E\left (x\right )$ that is a superposition of HG beams $E_m\left (x\right )$. This result is more general than the equations obtained in Refs. [24,25,50]. Mathematically $E_m\left (x\right )$ also describes the wave function of the energy states of a simple harmonic oscillator in quantum mechanics. Thus we use the bra-ket mathematical formalism to calculate expectation values as is standard in quantum mechanics. Consider a general state $\left | \psi \right \rangle$ described in the HG beam basis i.e., $\left | \psi \right \rangle = \sum \limits _{n = 0}^\infty {{c_n}} \left | n \right \rangle$ where $\left |{n}\right \rangle$ is the state of the HG mode of order n and $c_n$ are complex coefficients satisfying $\sum _{n=0}^{\infty } |c_n|^2 = 1.$

In the main text we found Eq. (6) as the basis of the covariance method. We use the general form for the covariance matrix elements in Appendix E, Eq. (22) to find $M^2$ from $\left | \psi \right \rangle$. It is useful to use dimensionless position $X$ and momentum $P$ (rather than angle) operators defined as follows $X = x/2\sigma _x$ and $P = p\sigma _x/\hbar = 2\pi \sigma _x\theta /\lambda.$ In terms of these dimensionless operators, Eq. (6) is rewritten as

(25)$$\begin{aligned} M^2 & = 4\sqrt{\det Q} \\ & = 4\sqrt{ \left(\Delta X\right)^2\left(\Delta P\right)^2 - \left( \langle{XP+PX }\rangle/2 - \langle{ X }\rangle \langle{P}\rangle \right)^2 }\\ & = 4\sqrt{ \left(\langle{ X^2 }\rangle-\langle{ X}^{2}\rangle\right)\left(\langle{ P^2 }\rangle-\langle{ P}^{2}\rangle\right) - \left( \langle{XP+PX }\rangle/2 - \langle{ X }\rangle \langle{P}\rangle \right)^2 }\\ & = 4\sqrt{\langle{ X^2 }\rangle\langle{ P^2 }\rangle - \langle{ X^2 }\rangle\langle{ P }^2\rangle - \langle{ X }^2\rangle\langle{ P^2}\rangle + 2 \langle{ X }\rangle\langle{ P }\rangle{\textrm{Re}} \langle{ XP }\rangle - \left({\textrm{Re}}\langle{ XP}\rangle\right)^2}, \end{aligned}$$

where we have used $\frac {1}{2}\langle {XP+PX}\rangle = {\textrm{Re}} \langle { XP }\rangle.$

The calculation of the expectation values appearing in above equation can be done using the ladder operators $a$ and $a^\dagger$ operators commonly used in quantum mechanics. These operators are related to $X$ and $P$ according to the following equations $X = \left ( a + a^\dagger \right )/2$ and $P = i\left ( a^\dagger - a \right )/2$. The action of these operators on the $\left | n \right \rangle$ states is given by $a\left | n \right \rangle = \sqrt {n}\left | {n-1} \right \rangle$ and $a^\dagger \left | n \right \rangle = \sqrt {n+1}\left | {n+1} \right \rangle$. Expectation values can now be easily obtained, for example $\langle {P}\rangle$ is calculated as follows:

(26)$$\begin{aligned}\langle{ P }\rangle & = \langle{ i\left( a^\dagger{-} a \right)/2\rangle }\\ & = \frac{i}{2}\sum_{n=0}^{\infty} \left( c_nc_{n+1}^* - c_n^*c_{n+1} \right) \sqrt{n+1}\\ & = -\textrm{Im}\left(B\right), \end{aligned}$$

where constants $A$, $B$, $C$ will be defined explicitly at the end. While $\langle {P^2}\rangle$ is given by

(27)$$\begin{aligned}\langle{P^2}\rangle & = -\langle{a^2 + \left(a^\dagger\right)^2 - 2\langle{n}\rangle - 1}\rangle/4\\ & = -\frac{1}{2} \sum_{n=0}^{\infty} {\textrm{Re}}\left(c_nc^*_{n+2}\right)\sqrt{ \left(n+1\right)\left(n+2\right)} + 1/2\sum_{n=0}^{\infty}n|c_n|^2 + 1/4\\ & = -\frac{1}{2}{\textrm{Re}}\left(A\right) + \frac{C}{4}. \end{aligned}$$

Following a similar procedure, one obtains $\langle {X}\rangle = {\textrm{Re}} \left (B\right )$, $\langle {X^2}\rangle = \frac {1}{2}{\textrm{Re}} \left (A\right ) + \frac {C}{4}$, and ${\textrm{Re}} \left (XP\right )=-\frac {1}{2}\textrm{Im} \left (A\right )$.

Substituting these into Eq. (25) we obtain the following general expression of $M^2$ for an arbitrary superposition of HG beams:

(28)$$M^2 = \sqrt{C^2 + 8{\textrm{Re}}\left(B^2A^* \right) - 4|A|^2 - 4|B|^2C},$$

where $A^*$ is the complex conjugate of $A$ and

(29)$$\begin{aligned}A & = & \sum_{n=0}^\infty c_nc^*_{n+2}\sqrt{ \left(n+1\right)\left(n+2\right)}, \end{aligned}$$

(30)$$\begin{aligned}B & = & \sum_{n=0}^\infty c_nc^*_{n+1} \sqrt{n+1}, \end{aligned}$$

(31)$$\begin{aligned}C & = & 1 + 2\sum_{n=0}^\infty n|c_n|^2. \end{aligned}$$

Funding

Canada Research Chairs; Canada First Research Excellence Fund; Natural Sciences and Engineering Research Council of Canada; Mitacs.

Acknowledgments

We thank Orad Reshef for publicly posting a request for the method provided in this paper; the absence of solutions in the replies showed that the community was unaware of our covariance method. This research was undertaken, in part, thanks to funding from the Natural Sciences and Engineering Research Council of Canada (NSERC), RGPIN-2020-05505, the Canada Research Chairs program, the Canada First Research Excellence Fund, and the MITACS Globalink Research Internship.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. A. E. Siegman, “New developments in laser resonators,” in Optical Resonators, Vol. 1224D. A. Holmes, ed., International Society for Optics and Photonics (SPIE, 1990), p. 2–14.

2. A. E. Siegman, “Defining, measuring, and optimizing laser beam quality,” Laser Reson. Coherent Opt. Model. Technol. Appl. 1868, 2–12 (1993). [CrossRef]  

3. A. E. Siegman, “How to (maybe) measure laser beam quality,” in Diode Pumped Solid State Lasers: Applications and Issues, (Optica Publishing Group, 1998), p. MQ1.

4. “Lasers and laser-related equipment - Test methods for laser beam widths, divergence angles and beam propagation ratios - Part 1”, ISO 11146-1:2021(E).

5. J. W. Goodman and R. W. Lawrence, “Digital image formation from electronically detected holograms,” Appl. Phys. Lett. 11(3), 77–79 (1967). [CrossRef]  

6. S. Grilli, P. Ferraro, S. D. Nicola, A. Finizio, G. Pierattini, and R. Meucci, “Whole optical wavefields reconstruction by digital holography,” Opt. Express 9(6), 294–302 (2001). [CrossRef]  

7. R. Boyd, S. Lukishova, and V. Zadkov, Quantum Photonics: Pioneering Advances and Emerging Applications, Springer Series in Optical Sciences (Springer International Publishing, 2019).

8. J. Hartmann, Objektivuntersuchungen, Springer Series in Optical Sciences (Springer, 1904).

9. R. V Shack, “How to (maybe) measure laser beam quality,” in Spring Meeting of Optical Society of America 1971, vol. 656 (1904).

10. B. Schäfer and K. Mann, “Determination of beam parameters and coherence properties of laser radiation by use of an extended Hartmann-Shack wave-front sensor,” Appl. Opt. 41(15), 2809–2817 (2002). [CrossRef]  

11. R. Sharma, J. S. Ivan, and C. Narayanamurthy, “Wave propagation analysis using the variance matrix,” J. Opt. Soc. Am. A 31(10), 2185–2191 (2014). [CrossRef]  

12. T. Kreis, Digital Recording and Numerical Reconstruction of Wave Fields (John Wiley & Sons, Ltd, 2004), chap. 3, pp. 81–183.

13. U. Schnars, C. Falldorf, J. Watson, and W. Jüptner, Digital Holography and Wavefront Sensing: Principles, Techniques and Applications, EBL-Schweitzer (Springer Berlin Heidelberg, 2014).

14. S. Wielandy, “Implications of higher-order mode content in large mode area fibers with good beam quality,” Opt. Express 15(23), 15402–15409 (2007). [CrossRef]  

15. S. Liao, M. Gong, and H. Zhang, “Theoretical calculation of beam quality factor of large-mode-area fiber amplifiers,” Laser Phys. 19(3), 437–444 (2009). [CrossRef]  

16. F. Stutzki, F. Jansen, H.-J. Otto, C. Jauregui, J. Limpert, and A. Tünnermann, “Designing advanced very-large-mode-area fibers for power scaling of fiber-laser systems,” Optica 1(4), 233–242 (2014). [CrossRef]  

17. H. Zhou, A. Jin, Z. Chen, B. Zhang, X. Zhou, S. Chen, J. Hou, and J. Chen, “Combined supercontinuum source with >200 W power using a 3 × 1 broadband fiber power combiner,” Opt. Lett. 40(16), 3810–3813 (2015). [CrossRef]  

18. N. Hodgson and T. Haase, “Beam parameters, mode structure and diffraction losses of slab lasers with unstable resonators,” Opt. Quantum Electron. 24(9), S903–S926 (1992). [CrossRef]  

19. C. Jauregui, C. Stihler, and J. Limpert, “Transverse mode instability,” Adv. Opt. Photonics 12(2), 429–484 (2020). [CrossRef]  

20. G. F. Marshall and G. E. Stutz, Handbook of optical and laser scanning, (Taylor & Francis, 2012).

21. A. Forbes, and Laser beam propagation: generation and propagation of customized light, (CRC Press, 2014).

22. S. Lavi, R. Prochaska, and E. Keren, “Generalized beam parameters and transformation laws for partially coherent light,” Appl. Opt. 27(17), 3696–3703 (1988). [CrossRef]  

23. A. E. Siegman, “Defining the effective radius of curvature for a nonideal optical beam,” IEEE J. Quantum Electron. 27(5), 1146–1148 (1991). [CrossRef]  

24. H. Weber, “Some historical and technical aspects of beam quality,” Opt. Quantum Electron. 24(9), S861–S864 (1992). [CrossRef]  

25. K. M. Du, G. Herziger, P. Loosen, and F. Rühl, “Coherence and intensity moments of laser light,” Opt. Quantum Electron. 24(9), S1081–S1093 (1992). [CrossRef]  

26. G. Nemes and A. E. Siegman, “Measurement of all ten second-order moments of an astigmatic beam by the use of rotating simple astigmatic (anamorphic) optics,” J. Opt. Soc. Am. A 11(8), 2257–2264 (1994). [CrossRef]  

27. C. Gao, H. Weber, and M. Gao, “Characterization of laser beams by using intensity moments,” in ICO20: Lasers and Laser Technologies, Vol. 6028 (SPIE, 2005), p. 428–435.

28. H. Yoda, P. Polynkin, and M. Mansuripur, “Beam quality factor of higher order modes in a step-index fiber,” J. Lightwave Technol. 24(3), 1350–1355 (2006). [CrossRef]  

29. O. A. Schmidt, C. Schulze, D. Flamm, R. Brüning, T. Kaiser, S. Schröter, and M. Duparré, “Real-time determination of laser beam quality by modal decomposition,” Opt. Express 19(7), 6741–6748 (2011). [CrossRef]  

30. D. Flamm, C. Schulze, R. Brüning, O. A. Schmidt, T. Kaiser, S. Schröter, and M. Duparré, “Fast M² measurement for fiber beams based on modal analysis,” Appl. Opt. 51(7), 987–993 (2012). [CrossRef]  

31. Y.-z. Du, G.-y. Feng, H.-r. Li, Z. Cai, H. Zhao, S.-h. Zhou, and M. Duparré, “Real-time determination of beam propagation factor by Mach-Zehnder point diffraction interferometer,” Opt. Commun. 287, 1–5 (2013). [CrossRef]  

32. Y. Du, Y. Fu, and L. Zheng, “Complex amplitude reconstruction for dynamic beam quality M² factor measurement with self-referencing interferometer wavefront sensor,” Appl. Opt. 55(36), 10180–10186 (2016). [CrossRef]  

33. M. Du, L. Loetgering, K. S. E. Eikema, and S. Witte, “Measuring laser beam quality, wavefronts, and lens aberrations using ptychography,” Opt. Express 28(4), 5022–5034 (2020). [CrossRef]  

34. M. Alonso, “Wigner functions in optics: describing beams as ray bundles and pulses as particle ensembles,” Adv. Opt. Photonics 3(4), 272–365 (2011). [CrossRef]  

35. M. Crampin and F. Pirani, Applicable Differential Geometry, Lecture note series (Cambridge University Press, 1986).

36. M. Hazewinkel, Encyclopaedia of Mathematics, no. v. 1 in Encyclopaedia of Mathematics (Springer Netherlands, 2012).

37. V. Dodonov and O. Man’ko, “Universal invariants of paraxial optical beams,” Computer Optics 1, 65–68 (1989).

38. R. Simon and N. Mukunda, “Optical phase space, Wigner representation, and invariant quality parameters,” J. Opt. Soc. Am. A 17(12), 2440–2463 (2000). [CrossRef]  

39. V. V. Dodonov and O. V. Man’ko, “Universal invariants of quantum-mechanical and optical systems,” J. Opt. Soc. Am. A 17(12), 2403–2410 (2000). [CrossRef]  

40. C. Fabre and N. Treps, “Modes and states in quantum optics,” Rev. Mod. Phys. 92(3), 035005 (2020). [CrossRef]  

41. M. J. Bastiaans, J. Flueckiger, L. Chrostowski, and D. Ratner, “Second-order moments of the Wigner distribution function in first-order optical systems,” in Optical Society of America Annual Meeting, (Optica Publishing Group, 1991).

42. Lasers and laser-related equipment - Test methods for laser beam widths, divergence angles and beam propagation ratios - Part 3 (ISO 11146-3: 2004), ISO 11146-3:2004(E).

43. A. Pettofrezzo, Matrices, Transformations (Prentice-Hall, 1966).

44. H. Weber, “Propagation of higher-order intensity moments in quadratic-index media,” Opt. Quantum Electron. 24(9), S1027–S1049 (1992). [CrossRef]  

45. S. Luo and B. Lü, “M2 factor and kurtosis parameter of super-Gaussian beams passing through an axicon,” Optik 114(5), 193–198 (2003). [CrossRef]  

46. C. S. Adams and I. G. Hughes, Optics f2f: From Fourier to Fresnel (Oxford University Press, 2018).

47. M. Griessmann, A. Martinez-Becerril, and J. Lundeen, “M2CovarianceMethod,” figshare (2023), https://doi.org/10.6084/m9.figshare.22661410.

48. A. Letsch and A. Giesen, “Characterization of a general astigmatic laser beam by measuring its ten second order moments,” in Optical Design and Engineering II, vol. International Society for Optics and Photonics5962 (SPIE, 2005).

49. Lasers and laser-related equipment - Test methods for laser beam widths, divergence angles and beam propagation ratios - Part 2: General astigmatic beams (ISO 11146-2: 2021), ISO 11146-2:2021(E).

50. K. Floettmann, “Coherent superposition of orthogonal Hermite–Gauss modes,” Opt. Commun. 505, 127537 (2022). [CrossRef]