Reconstructing the Spatial Parameters of a Laser Beam Using the Transport-of-Intensity Equation

A simple method for reconstructing the spatial parameters of a laser beam, based on the transport-of-intensity equation, is presented. Registration of cross-section intensity distributions in several planes was carried out using a single CMOS camera. The processing of the experimental measurements with the help of specialized software helped to reconstruct all of the spatial parameters, namely, the radius and position of the waist, Rayleigh length, angular divergence, quality parameter M2 The method was compared with measurements made according to the international standard ISO 11146 and showed that the difference in the spatial parameters is 10% or less, which shows good agreement.


Introduction
When characterizing conditionally static wave fields of light beams, their complex amplitude is usually represented in the form of two related components-amplitude (intensity) and phase distributions. A reliable estimation of the phase component of the field, often referred to as wavefront reconstruction [1][2][3] or a measurement of its aberrations [4,5], is usually the most difficult procedure. This is because one has to face problems of a possible ambiguity of the phase reconstruction [6], the presence of singular regions [7], a search for optimal boundary conditions [8] and limitations on the admissible dynamic range of phase gradients [9], which is typical for a particular measuring instrument.
Meanwhile, the coherent radiation formed by stable configurations of laser resonators is well described by a special class of light beams-Gaussian beams. Such beams have a field amplitude distribution in the form of Hermite-or Laguerre-Gauss functions, with a phase surface that has a paraboloid of rotation. In laser engineering and technologies, the difference of the spatial parameters of a real laser beam from an ideal Gaussian beam is characterized by the beam quality parameter M 2 .
Properties characterizing the propagation process of each laser beam are described by ten independent parameters [10,11]. The theoretical basis for finding these parameters is the method of determining second-order moments. However, due to the symmetry of most laser beams, which are widespread in many applications, their complete description requires fewer parameters, including the radius and position of the waist, Rayleigh length, angular divergence, quality parameter M 2 . When measuring the spatial and energy characteristics and parameters of laser radiation, it is necessary to form a beam with a valid waist to measure intensity distributions in at least 10 sections, of which five measurements should be taken near the waist (left and right of the waist plane). Then, for each section, the normalized moments of intensity distribution σ 2 x (z), σ 2 y (z) and the effective beam radius w e f f (z) = 2 σ 2 x (z) + σ 2 y (z) . By approximating the measurement results via a hyperbolic relationship, all the spatial characteristics of the laser beam above are calculated. However, despite the fact that the presented method is simple in its implementation, it is labor-intensive due to the need for multiple registrations of intensity distributions in the field under study.
Therefore, a number of alternative methods have been developed for more efficient M 2 measurements, including methods based on wavefront sensors [12,13], liquid crystal lens [14], computer-generated holograms [15], spatial light modulators [16], and several others.
At the same time, a method of real-time reconstruction of the complex amplitude exists for measuring the M 2 parameter by means of transport-of-intensity equation [17][18][19]. The transport-of-intensity equation (TIE) provides a non-interferometric way to obtain quantitative information on the phase of light beams by recording transverse field intensity distributions in two or more cross-sections [20,21]. The downside of the rapidity is that the measurement scheme requires two identical cameras [17], which will lead to additional errors and phase overruns, the magnitude of which will depend on the degree of identity of the selected cameras, the quality of the scheme assembly/alignment and the accuracy of the manufacture of the individual optical components. The disadvantages of the proposed approach in [18] include a large number of adjustable parameters, including the choice of the step between planes dz, the value of the shape factor for implementing the radial basis function (optimal approximation of the phase) and the Savitsky-Golay filter parameters. The determination of the optimal values of each of these parameters is a separate and very time-consuming task.
Therefore, this paper proposes a simple method for reconstructing the spatial parameters of a laser beam based on solving the TIE, using only one camera and a minimum number of optical elements. The method was compared with measurements made according to the international standard ISO 11146. Figure 1 shows a longitudinal section of an axisymmetric laser beam propagating along the z axis, with a waist diameter 2w 0 at a point with coordinate z = 0. The envelope (caustics) of the axisymmetric laser beam along the 1/e 2 level of the total flux is a singlecavity hyperboloid of rotation, and the phase surface represents a paraboloid of rotation. section, the normalized moments of intensity distribution ( ), ( ) and the effective beam radius ( ) = 2 ( ) + ( ) . By approximating the measurement results via a hyperbolic relationship, all the spatial characteristics of the laser beam above are calculated. However, despite the fact that the presented method is simple in its implementation, it is labor-intensive due to the need for multiple registrations of intensity distributions in the field under study. Therefore, a number of alternative methods have been developed for more efficient measurements, including methods based on wavefront sensors [12,13], liquid crystal lens [14], computer-generated holograms [15], spatial light modulators [16], and several others.

Methods and Approaches
At the same time, a method of real-time reconstruction of the complex amplitude exists for measuring the parameter by means of transport-of-intensity equation [17][18][19]. The transport-of-intensity equation (TIE) provides a non-interferometric way to obtain quantitative information on the phase of light beams by recording transverse field intensity distributions in two or more cross-sections [20,21]. The downside of the rapidity is that the measurement scheme requires two identical cameras [17], which will lead to additional errors and phase overruns, the magnitude of which will depend on the degree of identity of the selected cameras, the quality of the scheme assembly/alignment and the accuracy of the manufacture of the individual optical components. The disadvantages of the proposed approach in [18] include a large number of adjustable parameters, including the choice of the step between planes dz, the value of the shape factor for implementing the radial basis function (optimal approximation of the phase) and the Savitsky-Golay filter parameters. The determination of the optimal values of each of these parameters is a separate and very time-consuming task.
Therefore, this paper proposes a simple method for reconstructing the spatial parameters of a laser beam based on solving the TIE, using only one camera and a minimum number of optical elements. The method was compared with measurements made according to the international standard ISO 11146. Figure 1 shows a longitudinal section of an axisymmetric laser beam propagating along the z axis, with a waist diameter 2 at a point with coordinate z = 0. The envelope (caustics) of the axisymmetric laser beam along the 1 ⁄ level of the total flux is a singlecavity hyperboloid of rotation, and the phase surface represents a paraboloid of rotation.  Since the processes of lase-matter interaction are three-dimensional, not only is the diameter of the focused spot (waist diameter) important, but also the characteristic waist length of the beam or, as experts in laser material processing term it, "depth of focus" or "waist length". This value is determined by the Rayleigh length of the laser beam z R and is its doubled value (see Figure 1).

Methods and Approaches
In turn, the Rayleigh length is the distance from the waist of the beam at which the beam size increases √ 2 times, i.e., the radiation power density at a distance ±z R from the beam waist decreases by a factor of 2. Therefore, in the longitudinal direction, the laser beam can be divided into three parts. In the central part, for |z| ≤ z R , the transverse dimensions of the beam vary relatively little with changes in z. In the two peripheral parts (for z < −z R and z > z R ), however, the dimensions of the transverse beam increase noticeably with increasing z and are proportional to |z| for larger z. Therefore, the Rayleigh length of the laser beam also determines the so-called «near zone length» of the beam. The dependences of the laser beam radius w(z) and the radius of curvature of its wavefront R(z) on the longitudinal coordinate z are described by the following equations [22]: where z coordinate is measured from the beam waist position. In order to characterize lasers whose beam quality is far from ideal (including highpower fiber and disk lasers), laser manufacturers widely use the Beam Parameter Product (BPP), [mm·mrad]: In this relation, λ denotes the wavelength of the laser radiation; M 2 is a dimensionless parameter that determines the quality of a beam and is usually defined as the ratio of the BPP of a real laser beam to that of a diffraction-limited Gaussian beam at the same wavelength. For an ideal TEM 00 Gaussian beam, M 2 = 1, and for a real beam M 2 > 1.
Rayleigh length and angular divergence of the beam are calculated using the following equations: These equations show that for a laser beam with a given wavelength λ and parameter M 2 , one cannot simultaneously reduce the waist size and the angular divergence of the beam. Thus, decreasing the waist size results in decreasing the Rayleigh length and increasing the angular divergence. At the same time, it is possible to reduce only the waist size and angular divergence by increasing the refractive index of the medium in which the beam is formed. In addition, increasing the BPP of the beam increases the diameter of the focused spot and decreases the length of the divergence with respect to the ideal Gaussian beam, for which BPP G = λ/π.

ISO Method
When processing the measurement results, it is convenient to represent the hyperbolic dependence of the laser beam radius as w(z) = √ a + bz + cz 2 [10]. The coefficients a, b, c could be found by the least squares method using the following matrix expression [23]: where w i is a beam radius in the plane with coordinate z i ; n is a number of planes in which the radiation power density distribution was measured (n ≥ 3); «−1» denotes inverse matrix notation. Using the coefficients a, b, c, we calculate the spatial parameters of the beam, the relationship of which is shown in Table 1. Table 1. Functional dependences of spatial parameters of the laser beam on a, b, c.

Parameters Equations
Radius of the laser beam waist c Position of the beam waist relative to the selected reference plane (taking into account the sign rule adopted in optics) In turn, the calculation of the curvature radius of the wavefront can be based on the solution of the TIE, in which the field phase function φ(r ⊥ , z) is related to the intensity distributions I(r ⊥ , z) by the following equation [19,24]: One of the most popular methods for solving this equation is the method based on expressing differential operators through the Fourier transform [25], that is, where r ⊥ = (x, y) is the radius-vector; F and F −1 are the direct and inverse two-dimensional Fourier transform, respectively, and k x,y = 2πν x,y , where ν x,y is the coordinate grid in the frequency domain; k = 2π/λ is the wave number. Figure 2 shows a scheme of experimental measurements and an algorithm for calculating the spatial parameters of the laser beam by solving the TIE, which consists of three main stages. At the first stage, two intensity distributions I(r ⊥ , z 1 ) and I(r ⊥ , z 3 ) are recorded using a camera. Then, at the second stage, the phase φ(r ⊥ , z) is calculated using the registered intensity distributions by solving the TIE, and then the wavefront curvature radius R z 2 is calculated using the geometric method based on the phase data [19]. At the last stage, by the least-squares method the spatial parameters of the laser beam are calculated from the values of R z 2 according to the relations given in Table 1. main stages. At the first stage, two intensity distributions ( , ) and ( , ) are recorded using a camera. Then, at the second stage, the phase ( , ) is calculated using the registered intensity distributions by solving the TIE, and then the wavefront curvature radius is calculated using the geometric method based on the phase data [19]. At the last stage, by the least-squares method the spatial parameters of the laser beam are calculated from the values of according to the relations given in Table 1. Figure 2. Schematic of the experimental measurements and algorithm for calculating the spatial parameters of the laser beam.

Experimental Part
As a rule, for the experimental measurement of all spatial parameters of the beam, an additional optical system (lens) was used and radiation intensity distributions were measured in several cross sections of the beam transformed by this lens. Therefore, in the experiment, horizontally polarized radiation from a Satsuma femtosecond laser with a second harmonic wavelength = 515.6 nm, pulse duration of 10 ps, and frequency of 10 kHz was focused, using a quartz flat-convex lens with a focal distance ′ = 550 mm. The specifications for the Satsuma femtosecond laser are shown in Table 2.
The intensity distributions were registered using a CMOS camera (1920 × 1080 pixels with dimensions 5.04 × 5.04 μm) mounted on a linear motion platform that moved the CMOS camera along the radiation propagation axis (Figure 3). A view of the two-dimensional sections of the intensity distribution in the X0Y (a), (b) and (c) planes obtained using the camera is shown in Figure 4a-c, respectively. Based on such intensity distributions, the longitudinal intensity section of the axisymmetric laser beam was reconstructed (Figure 4d).

Experimental Part
As a rule, for the experimental measurement of all spatial parameters of the beam, an additional optical system (lens) was used and radiation intensity distributions were measured in several cross sections of the beam transformed by this lens. Therefore, in the experiment, horizontally polarized radiation from a Satsuma femtosecond laser with a second harmonic wavelength λ las = 515.6 nm, pulse duration of 10 ps, and frequency of 10 kHz was focused, using a quartz flat-convex lens with a focal distance f = 550 mm. The specifications for the Satsuma femtosecond laser are shown in Table 2.
The intensity distributions were registered using a CMOS camera (1920 × 1080 pixels with dimensions 5.04 × 5.04 µm) mounted on a linear motion platform that moved the CMOS camera along the radiation propagation axis (Figure 3). A view of the twodimensional sections of the intensity distribution in the X0Y (a), (b) and (c) planes obtained using the camera is shown in Figure 4a-c, respectively. Based on such intensity distributions, the longitudinal intensity section of the axisymmetric laser beam was reconstructed (Figure 4d).

ISO Method
As noted above, the envelope (caustic) of an axisymmetric laser beam in terms of the level 1 ⁄ of the total flux is a single-cavity hyperboloid. Based on the obtained intensity

ISO Method
As noted above, the envelope (caustic) of an axisymmetric laser beam in terms of the level 1/e 2 of the total flux is a single-cavity hyperboloid. Based on the obtained intensity distribution in the Y0Z plane (Figure 4d), one can approximate and represent the dependence for the laser beam radius (axis in Y) as w(z) = √ a + bz + cz 2 [10]. A total of 143 cross sections obtained in the experiment with a step of 1 mm were used to construct the envelope. The coefficients a, b, c were found using the matrix expression (4) and were equal to the following: a = 0.058 mm 2 , b = −1.307 × 10 −3 mm, c = 7.866 × 10 −6 for plane X0Z and a = 0.059 mm 2 , b = −1.458 × 10 −3 mm, c = 9.506 × 10 −6 for plane Y0Z. Then, these coefficients were used to calculate the spatial parameters of the beam, the relationship of which is given in Table 1. The numerical results are presented in Table 3.

TIE Method (Proposed Method)
To find the spatial parameters of the laser beam by the proposed method, it is first necessary to register at least three transverse intensity distributions in arbitrary planes. In this case, the extreme planes will be denoted as z 1 and z 3 , and the intermediate one between them is z 2 . For a comparative analysis of the methods, we used the same experimental data obtained for the ISO method.
At the second stage, the phase function of the field φ(r ⊥ , z 2 ) was calculated using TIE for the intensity distributions in the two extreme planes z 1 and z 3 (the distance between the planes ∆z 1−3 was chosen arbitrarily and amounted to 30 mm). The resulting phase function of the field φ(r ⊥ , z 2 ) in the intermediate plane z 2 was used for the approximation and calculation of the wavefront curvature radius R(z 2 ) [19]. The radius of curvature was calculated for the different pairs of planes z 1 and z 3 from the obtained dataset. The results are shown in Figure 5a.
Next, the values of the Rayleigh length z R were obtained for the longitudinal planes X0Z and Y0Z using Equation (2). Figure 5b shows a histogram of the obtained values with lines showing their approximation. Then, the waist radius w 0 was calculated using Equation (1). The results of the calculation are presented in Figure 5c. After that, the parameter M 2 was calculated according to Equation (4). The results are shown in Figure 5d. The comparison of the spatial parameters reconstructed by the two methods is shown in Table 3. It can be seen that the methods are in good agreement with each other.
It should be noted that in this work we investigated the performance of our method with a single-mode single-frequency laser, with parameter M 2 close to 1. We will continue our research, and in the near future we plan to study the multimode beams. The purpose of this research will be to analyze the boundaries of the method applicability. We aim to configure where, first of all, it will be necessary to measure the stability of the method to fluctuations of the output laser power. However, we can already argue that the proposed method is not inferior to the ISO 11146 one in terms of measuring M 2 close to 1 (Table 3), while being easier to implement. Additionally, the proposed method is not inferior to existing commercial solutions, for example, the Thorlabs M2MS system. distribution in the Y0Z plane (Figure 4d), one can approximate and represent the dependence for the laser beam radius (axis in Y) as ( ) = √ + + [10]. A total of 143 cross sections obtained in the experiment with a step of 1 mm were used to construct the envelope. The coefficients a, b, c were found using the matrix expression (4) and were equal to the following: = 0.058 mm 2 , = −1.307 • 10 mm, = 7.866 • 10 for plane X0Z and = 0.059 mm 2 , = −1.458 • 10 mm, = 9.506 • 10 for plane Y0Z. Then, these coefficients were used to calculate the spatial parameters of the beam, the relationship of which is given in Table 1. The numerical results are presented in Table 3.

TIE Method (Proposed Method)
To find the spatial parameters of the laser beam by the proposed method, it is first necessary to register at least three transverse intensity distributions in arbitrary planes. In this case, the extreme planes will be denoted as and , and the intermediate one between them is . For a comparative analysis of the methods, we used the same experimental data obtained for the ISO method.
At the second stage, the phase function of the field ( , ) was calculated using TIE for the intensity distributions in the two extreme planes and (the distance between the planes Δ was chosen arbitrarily and amounted to 30 mm). The resulting phase function of the field ( , ) in the intermediate plane was used for the approximation and calculation of the wavefront curvature radius ( ) [19]. The radius of curvature was calculated for the different pairs of planes and from the obtained dataset. The results are shown in Figure 5a.
Next, the values of the Rayleigh length were obtained for the longitudinal planes X0Z and Y0Z using Equation (2). Figure 5b shows a histogram of the obtained values with lines showing their approximation. Then, the waist radius was calculated using Equation (1). The results of the calculation are presented in Figure 5c. After that, the parameter was calculated according to Equation (4). The results are shown in Figure 5d. The comparison of the spatial parameters reconstructed by the two methods is shown in Table  3. It can be seen that the methods are in good agreement with each other. It should be noted that in this work we investigated the performance of our method with a single-mode single-frequency laser, with parameter close to 1. We will continue our research, and in the near future we plan to study the multimode beams. The purpose of this research will be to analyze the boundaries of the method applicability. We aim to configure where, first of all, it will be necessary to measure the stability of the method to fluctuations of the output laser power. However, we can already argue that the proposed method is not inferior to the ISO 11146 one in terms of measuring close to 1 (Table 3), while being easier to implement. Additionally, the proposed method is not inferior to existing commercial solutions, for example, the Thorlabs M2MS system.
The measurement time of parameter amounted to 30 s, which is optimal, yet with current algorithmic solutions it can be further optimized to achieve better performance. The expected improvement from the algorithm optimization can reduce the measurement time considerably. Additionally, in combination with the simple installation of the elements required for the implementation of the method, the method can become a widely used alternative solution for many photonics problems.

Conclusions
This work presents a simple method for reconstructing all of the laser beam spatial parameters based on the solution of the TIE. The phase function of the field can be quickly recovered by solving the TIE and does not require the use of a complex optical system. The method was compared with measurements made according to the international standard ISO 11146 and showed that the difference of the spatial parameters is 10% or less, which means that the methods are in good agreement with each other. The advantages of this approach are its experimental simplicity and the possibility to perform quick measurements without time-consuming numerical analyses. In addition, this method is suitable for determining the beam quality of a dynamic laser system. The measurement time of parameter M 2 amounted to 30 s, which is optimal, yet with current algorithmic solutions it can be further optimized to achieve better performance. The expected improvement from the algorithm optimization can reduce the measurement time considerably. Additionally, in combination with the simple installation of the elements required for the implementation of the method, the method can become a widely used alternative solution for many photonics problems.

Conclusions
This work presents a simple method for reconstructing all of the laser beam spatial parameters based on the solution of the TIE. The phase function of the field can be quickly recovered by solving the TIE and does not require the use of a complex optical system. The method was compared with measurements made according to the international standard ISO 11146 and showed that the difference of the spatial parameters is 10% or less, which means that the methods are in good agreement with each other. The advantages of this approach are its experimental simplicity and the possibility to perform quick measurements without time-consuming numerical analyses. In addition, this method is suitable for determining the beam quality of a dynamic laser system.