A method to determine the M2 beam quality from the electric field in a single plane

Laser beam quality is a key parameter for both industry and science. However, the most common measure, the M2 parameter, requires numerous intensity spatial-profiles for its determination. This is particularly inconvenient for modelling the impact of photonic devices on M2, such as metalenses and thin-film stacks, since models typically output a single electric field spatial-profile. Such a profile is also commonly determined in experiments from e.g., Shack-Hartmann sensors, shear plates, or off-axis holography. We introduce and test the validity and limitations of an explicit method to calculate M2 from a single electric field spatial-profile of the beam in any chosen transverse plane along the propagation direction.


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][2][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].A reliable way to measure M 2 was codified in 1999 by the ISO standard 11146 [4], which called for transverse spatial intensityprofiles to be recorded at ten distances along the propagation axis.The M 2 numerically validate such a method for two nontrivial test beams.
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 rad).We give explicit formulas to calculate M 2 , the beam waist w 0 and its location z 0 , the beam's angular width ∆θ, 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

Definitions and conventions
For simplicity we will consider beam profiles in one transverse dimension only, although results generalize to two dimensions.The z axis is the propagation direction while x-axis is the transverse direction, along the beam's spatial profile.As a beam propagates along z, generally its width in position w (z) will change while its angular width ∆θ will be constant.We use the standard width conventions w (z) = ∆x (z) = 2 σ x (z) and ∆θ = 2 σ θ , where the σ are the standard deviation of corresponding intensity distribution (i.e., in x or θ).We take z = 0 as the position of the beam waist w 0 , the beam's minimum w (z) 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

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 ∆θ.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 λ/π [3] we have, Whereas the product w (z) ∆θ will change as the beam propagates, since w 0 and ∆θ 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. 1a, the relation between the widths in angle (∆θ), transverse wavevector (∆k x ) and transverse momentum (∆p x ) of a photon can be expressed as where k = 2π/λ is the total wavevector magnitude, λ is the wavelength, and the small angle approximation was applied.Using Eq. ( 2), one can rewrite M 2 as where the inequality comes from the Heisenberg uncertainty principle: σ x σ px ≥ h/2 or, equivalently, ∆x∆p x ≥ 2h.Gaussian beams are minimum uncertainty states and thus saturate inequality (3), M 2 = 1, while all other beams give strict inequality.Although Eq. (1) for M 2 seems quite simple, determining its value is deceptively difficult.The angular width ∆θ 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 (z) 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 definitions (Eq.( 1)) are hardly applicable in either experiment or simulation.

The 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 θ on the perpendicular axis (though the latter could equally well be k x or p x , as explained in the last section).In this phasespace, the beam profile at any z-position can be visualized as a tilted ellipse (see Fig. 1c).The length of the projection of the ellipse on the x-axis is the full-width intensity standard deviation, w (z) = ∆x.Similarly, the projection on the θ-axis is the angular beam width, ∆θ .
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 .The area A of the ellipse is proportional to its width times its height, A = 1 4 πw 0 ∆θ.In turn, using Eq. ( 1), we find that the beam quality is proportional to the ellipse area, As the beam propagates away from the waist location, correlations between position and angle form and the ellipse becomes tilted.This follows from simple geometry; waves traveling at at angle θ will increase in position as x = θz in the small angle approximation, creating a correlation between x and θ.This shears the ellipse along the x-direction in phase-space [32] as seen in Fig. 1d.(Conversely, momentum conservation in free space ensures that each constituent wave's angle θ is constant.)However, since the ellipse is now tilted, these widths are distinct from the width and height of the ellipse itself i.e., along its semiminor and semi-major axes.Crucially, in geometry, a shear transformation always preserves area [33,34] and, thus, M 2 .This implies that Eq. ( 4) is true for all z, not just at the waist [35][36][37].Thus, the problem of determining M 2 from information in a single z-plane reduces to finding the area of the corresponding tilted ellipse.

The covariance matrix
We gather the parameters of a tilted ellipse centered at the origin (i.e., the ellipse equation, ax 2 − 2bxθ + cθ 2 = ac − b 2 ) into a symmetric matrix [38], 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, where is analogous to an average or expectation value but evaluated using the complex distribution E (x; z), the transverse profile of the beam's electric field at a plane z.We give further detail on calculating at the end of this subsection and explicit expressions for Q will be given in Section 2.5.For now, we point out that the diagonals of Q are the variances, 2 /4 for a beam with x = θ = 0.When the beam has reached its waist z = 0, the ellipse is aligned with the phase-space axes and Q is diagonal.As explained in the previous subsection, at other z, correlations exist between angle and position.These are the off-diagonal covariance terms in Q.The Q matrix is closely related to beam matrices defined in terms of the Wigner function; namely, it is the Weyl transform of the matrices in [39] and Part 3 of the ISO standard [40].Unlike those matrices, Q is directly computed from the electric field distribution E (x; z) in a single z-plane.
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 A = π √ 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 [39].In summary, one can evaluate M 2 using the elements of the covariance matrix Eq. ( 5).These expectation values use the complex field E (x; z) (in analogous way to the quantum wavefunction).That is, the expectation value of a general operator v acting on E is defined by v = 1 n E * vEdx, where E * is the complex conjugate of the electric field, and n = |E (x; z) | 2 dx is a normalization factor.The angle operator is θ = − i k ∂ ∂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 centred at the origin, i.e., x = 0 and/or θ = 0. We can account for this by substituting x and θ in the covariance matrix (Eq.( 5)) with x − x and θ − θ , respectively.We give explicit formulae for elements a, b, and c for such offset beams in Appendix D. 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.

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 ≡ E (x; z) 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 x = θ = 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 D.
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 ∆θ is twice the intensity standard-deviation (e.g., w 0 = 2σ x ); and the beam quality is To evaluate this formula, one uses the covariance matrix elements, all of which are real-valued: where we have used ∂/∂x (xE) = x∂E/∂x + E to simplify b.One can save some computational effort my noticing that the second and first integrals in a and c, respectively, also appear in b.
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 ∆θ, Rayleigh range z R , beam waist w 0 , and waist's location relative to the plane of the measured field, ∆z.In Eqs.(14)(15)(16)(17)(18) 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.

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 (x) at an initial plane (z = 0) using the analytic form of the test beam E (x) and then numerically propagate it to find E (x; z) and the corresponding intensity profile I (x; z) = |E (x; z)| 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 (x; z) to find the beam width w (z) 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. ( 17) to find M 2 .
At each of these planes, we calculate the entries, a, b, and c of the covariance matrix Q from E (x; z) 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.

Numerical comparison of the ISO and covariance methods
For the following test beams, λ = 400 nm, we use position range of length L = 200 mm and a position grid density δx = λ/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 Figure 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, λ = 400 nm M 2 = 2.466 given by Eq. ( 27) ): (a) points are the scaled entries of the covariance matrix for the beam numerically propagated (by Eq. ( 19)) to each z distance and curves are the analytically propagated elements (by Eq. ( 13)) and (b) is the corresponding intensity profiles at three z-locations.For a paraxial supergaussian beam (w 0 = 56 µm, order = 50, λ = 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. (17).
the initial plane for the m-th HG mode is taken to be where w 0,m is the beam waist of the mode and H m the m-th order Hermite polynomial.Although each HG mode is centered the total beam will generally be offset from the x = 0 axis, so Eqs.(20,21) in Appendix E must be used.In addition, the z location of the waist of a superposition will not generally be at the waist location of the individual HG modes.Accommodating these possibilities, we derive a completely general theoretical value of M 2 for a superposition of HG modes with complex coefficients c n in Appendix F, with Eq. ( 27) as the final result.
For this first test beam we plot the elements of the offset Q in Fig. 2a.The element a = w 2 (z) /4 (blue points) lies on the theoretical beam waist curve (Eq.( 17), 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.( 13), yellow and black curves).This agreement confirms that our numerical propagation of E (x) is accurate and our calculation of the covariance matrix Q is correct.
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,21)) 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 (x; z) 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 × 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 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 sincsquared function and thus has a ∆θ that diverges and is not well-defined, which makes it a particularly challenging test case.A one-dimensional supergaussian, E (x) = exp − (|x|/w s ) N , of high order N approximates a top hat function of full-width 2w s .Using Eqs. ( 6) and ( 9) we theoretically found , where the approximation is valid for large N and Γ is the usual gamma function.(Note, this formula differs from formulae in the literature, which are for the "cylindrical M 2 ", e.g., [42].) The test beam is a supergaussian of order N = 50 with equivalent slit width of 2w s = 0.1 mm and a theoretical M 2 = 4.0466.Diffraction from such a slit has zero-to-zero angular width of 1.6 × 10 −2 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. 2c.The ten planes of our covariance method give an average value of M 2 = 4.0462 with a standard deviation of 7.0 × 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 × 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.

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 (x) E (y), e.g., two dimensional HG modes.In Appendix E, we generalize 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.
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.

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 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. ( 12) to propagate an arbitrary Q (z) by a distance z , From this matrix we can retrieve the common beam parameters.First, we note that c is unchanged so, trivially, 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 ≡ −∆z, we then find where a positive ∆z means the waist is further along the beam propagation direction from the plane that E (x; z) (and Q) was determined in.With z = −b/c substituted in Q Eq. ( 13), we find the beam waist, We now instead set Q (z 0 ) to be the waist and, thus, diagonal (b = 0), and using M 2 from Eq. ( 1), the a element of Q gives the beam width: where z = z − z 0 .This gives the Rayleigh range, In other words, the propagation of an arbitrary paraxial beam is same as the one of a Gaussian beam magnified by M 2 .Equation (17) 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. ( 17).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.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 nonparaxial version of the Wigner function [32], but a thorough investigation would be in order here.

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 beam is perfectly centered on the beam axis, i.e., x = 0 and θ = 0.In this case, the covariance matrix takes the more general form: The corresponding generalized formulas for the matrix elements are 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 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 (x, y; z) = E x (x; z) E y (y; z), and each transverse direction can be separately analyzed by our covariance method.The result would be two values, M 2 x and M 2 y , from which one can define an effective beam quality M 2 eff = M 2 x M 2 y .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 (x, y; z) = E x (x; z) E y (y; z).Unlike simple astigmatic beams, the principal axes of their transverse cross sections can rotate while the beam is propagating [44].Part 2 of the ISO standard [45] derives an effective beam quality M 2 eff 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 2 eff is the area of that state.Analogous to the generalized P matrix in Part 2 of the ISO standard, a generalized Q xy (z) matrix would be a 4x4 matrix: Here, the block diagonals are the standard 2x2 Q matrices Eq. ( 5): Note, the D blocks involve one variable from each direction and thus do not require the two orderings (e.g., xθ + θx ) that are essential for the 2x2 Q matrix.With this generalized covariance matrix, M 2 eff = 4π λ 4 det Q xy (z).Unlike the P matrix in the ISO procedure, Q xy (z) is found directly from the electric field E (x, y; z).

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, ρ (x, x ) ≡ E (x) E * (x ).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 (z) would then be v = Tr [vρ] 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,24], we refer the reader to such references for further details.

Figure 1 :
Figure 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, λ = 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.
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 D xy (z) = xy xθ y yθ x θ x θ y .