Simulation of defocusing effect based on two-step ABCD algorithm while a modal decomposition

Depicting the multimode laser beam by modal decomposition can potentially assess light field variations in the fiber, during propagation. The practical engineering conditions in the lab however could not realize ideal levels, hence further research on factors influencing this method, such as defocus, is especially necessitated. The grid spacing in observation plane by Fast Fourier Transform is fixed and unchangeable within diffraction imaging, hence possibly yielding erroneous data during obtaining light field intensities. Our research resolves these issues via a Two-step ABCD algorithm, applied in the modal decomposition to characterize various guided modes at the output of multimode fibers. A direct benefit is that the image plane size can be altered, further refining laser facula clarity. Furthermore, the quantitative expressions that analyze defocus factors impacting modal decomposition are acquired. The conclusions thereby prove the modal decomposition algorithm can keep effectiveness in the range of −0.25% to 0.25% of relative defocus for low order eigenmodes, having no suitable limited band for high order eigenmodes, with reference value in engineering applications.


Introduction
Further research and development have facilitated several laser beams technology applications in engineering, hugely capitalizing on the high accuracy, efficiency and power. Laser beam propagation in step optical fiber is a subject persistently researched in recent years, such that the content is resourcefully characterized for an enhanced and precise propagation in several fields, such as fine mechanism [1,2], Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. detection [3,4] and transmission [5,6], to name a few. Several established methods are used to measure a range of light field information to identify laser beams, in corresponding conditions. The S 2 measurement spatially and spectrally quantifies the number and type of modes, to analysis the mode components [7,8]. Pu Zhou and Haibin LÜ and Liangjin Huang proposed an approach based on stochastic parallel gradient descent algorithm for a comprehensive optical fiber modal decomposition [9,10]. Maxim V. Bolshakov verified a numerical algorithm based on the cross-correlation operation of light intensity for modal power decomposition [11]. Thomas Kaiser presented a computer-generated holographic filter to characterize the laser beam, even if it had singularities [12,13]. Amongst the above-mentioned methods, the direct light field description can potentially be applied to analyze optical physical peculiarity, and comprehend laser mode interactions, namely mode competition, curve loss and laser beam quality.
Notwithstanding these advantages of laser beams, some issues persist that need to be resolved. A power increase in the optical fiber leads to nonlinear effects, such as Raman scattering and Brillouin scattering. Large-mode-area fibers can solve this issue, however not without unexpectedly irradiating the high-order modes. The mix of diverse modes worsens beam quality, decreases stability, and destroys energy density. In fact, when a laser beam comprises a few different modes, the M 2 factor is also close to one, and therefore cannot accurately depict the same. This principle also applies to optical fiber lasers. In the field of optical fiber laser [14][15][16], maintaining high beam quality during propagation is of capital importance and significant by various methods. Hence, identifying the content of modes is essential, for an insight on their interaction. Field fluctuations, directed by the shifting of laser beam mode content, are clearly observed and analyzed, particularly when disturbances occur in the wave-guided device.
Mode decomposing is primarily essential, before laser beam content can be identified. Thomas Kaiser et al (reference [13]) proposed an optical correlation analysis filter called MODAN, short for Multimode Analysis, based on computer generated hologram (CGH), to decompose and reconstruct modes in fiber. The quintessence in Thomas Kaiser's paper is its unique and straightforward mathematical procedure. Meng Lyu et al employs digital holography to measure the light field at the output of the multimode optical fiber [17]. They calculated the modal coefficients of each mode based on the modal orthonormal property. Nevertheless, the Fast Fourier Transform utilized in most modal decomposition methods, cannot alter the spacing in the frequency domain, which is limited by Nyquist's theorem. Consequently, the light field by FFT in simulation was very tiny and concentrated, to distinguish any distribution patterns. The two-step ABCD algorithm [18] integrated Fresnel propagation with the Collins formula to arbitrarily adjust the image plane sampling rate, via a single alternation of a Fast Fourier transform and its inverse. In this research, we aim to apply the Two-step ABCD algorithm in modal decomposition theory to improve its accuracy. On this basis, the influence of defocusing on modal decomposition theory is analyzed by simulation.

Theory of modal decomposition
The optical fiber guide modes are mathematical solutions of Helmholtz and Wave equations. The transverse eigenmode expressions, varying within stipulated conditions, include linear polarization (LP) scale mode in a weak guide fiber, Hermit-Gauss (HG) mode in a resonant cavity with square mirror, and Laguerre-Gauss (LG) mode in a circular mirror cavity, collectively known as eigenmodes, and exhibiting orthogonality, totality and stability. Any incident light field, owing to its characteristics, can be expressed as superposition of eigenmodes. The research we conducted is on the premise of a weak guide fiber, wherein the refractive index of the core is approximating cladding's number, with the difference being about 0.001. This condition yields an LP scale eigenmode.
The illuminating field is designated as, where u is the Cartesian coordinate. The complex coefficient C n is the weight of modes, comprising the relative amplitude and intramode phase.
Hence, we can ascertain the modes permitted in fiber, based on the V number. Eigenmodes are vectors with usually two different polarizations that do not affect the other; degeneracy in the weakly guided fiber occurs along the polarizing direction, forming the LP scale modes. The eigenmode orthogonality is designated as where the asterisk denotes the complex conjugate. In the references, a diffractive optical element is designed as a superposition of eigenmode conjugate, which conveys the calculated frequency-Vn.
where T n is transmittance function of correlation filter. It shows the transmittance function is the sum of all eigenmode conjugate with different carried frequencies. The Fourier frequency shift theorem, shown in equation (6), isolates the different frequency-shift modes in the spatial spectrum field. In other words, the MODAN filter transmittance function is modulated via angular multiplexing. The mixture of modes passes through the regular filter element directions and arrives at the image plane, to realize mode decomposition.
We analyzed each optical path as a separate branch, so as to simultaneously and synchronously comprehend a number of channels, as also their interferometric superposition. Finally, specifying different output patterns, primarily the near-field and far-field distributions; near-field signifies the rear surface of the filter, which is associated with the intensity of all modes, whilst far-field is where the Fourier transform is achieved. We set up a lens for Fourier transform, so the back focal plane is equivalent to the far-field distribution [13].
where f is the focal length of the lens. The symbol '~' denotes the frequency domain.
It is noteworthy that the far-field description comprises several superpositions of cross correlations of each conjugated eigenmode. The orthogonality satisfies the complex amplitudes of different modes at precisely the Un = Vn * f/k spots, and as a result, patterns do not entirely conform with the initial eigenmodes. A CCD camera is placed at specific points on the rear focal plane for light intensity, which is proportional to the mode weight squared. The relative phase and M 2 factor can also be acquired via another algorithm, based on the intensity data. Generally, LP 01 is defined as basic mode, and its initial phase is zero.

The two-step ABCD algorithm
In the Two-step ABCD algorithm, a virtual medium plane was set during propagation, such that laser beams satisfy the Fresnel diffraction, from both the source to the medium, and medium to the image plane. The two-step Fresnel diffraction can realize flexible spacing grid in image plane. From modal decomposition theory, it reveals that the principle of separating eigenmodes is based on the Fourier transform with serval carried frequencies. As defocus existing in detecting distance, the function of lens is Fresnel transform, instead of the Fourier transform. Hence, the Collin's formula is more convenient and more direct way to operate both two transforms, only changing the distance between lens and detector. As a consequence, in the case of this investigation, the Two-step ABCD algorithm is more suitable way to change the sampling intervals and analyze the effect of defocus.
The medium plane herein could rest between the source and observation planes, or even be far-removed. For simplicity, we derived the field in only one instance, wherein the propagations considered are from source plane U 1 to medium plane U d , and from observation plane U 2 to medium plane U d , as indicated in figure 2. A Fresnel diffraction from U 1 (x,y) presents U d (u,v) as: where k is the wave number, and Z 1 the propagation distance.
We similarly ascertain U d (u,v) from U 2 (ξ,η), via the expression We designate ∆ 1 as sampling length in U 1 and ∆ d,1 as sampling length in U d . N * N is the grid dimension in all planes. L 1 denotes the width of the subject plane, and B x its cut-off frequency in the spectrum domain of the image plane. Laser source-Nd:YAG laser; BS-beam splitter; CCD1-detecting far-field intensity, CCD2-detecting near-field intensity. In simulation, we supposed that the laser beam is collimated and broadened input and output from optical fiber. MODAN is the filter based on self-correlation theory.
Equating the right-hand sides of equations (1) and (2) helps correlate U 2 (ξ,η) and U 1 (x,y), at the exclusion of U d (u,v). We can obtain the U 2 (ξ,η) field via the expression where F is the Fourier-transform operator, and F −1 its inverse. Equation (11) is the outcome of a two-step Fresnel diffraction for scale paraxial propagation. A similar equation for sampling distance in U d replaces ∆ 1 with ∆ 2 .
The two-step method works numerically, only with U d in the same position, starting from either U 1 or U 2 . Equations (10) and (12) must therefore be equal, so we can acquire We adjust the sampling distance ∆ 2 to ultimately obtain the variable reconstruction image size, which also connects with the distances Z 1 and Z 2 .
Here, we combine this two-step diffraction algorithm with the Collin's formula to adjust the system parameters more easily. The laser beam transferring through the optical system can be represented as the Collin's formula in polar coordinate system. In while context, the bold terms are vector and containing the angular information, such as r 1 and r 2 . We set two variables Schematic diagram of variable sampling display system. U 1 is incident field; U d is the visual middle plane; U 2 is the imaging plane. Both from U 1 and U 2 to U d are Fresnel diffraction.
where A,B,C and D are the elements of system matrix and satisfy that AD-BC = 1. It is worth to know that this is valid only for optical system with azimuthal symmetry. We refer the Collin's formula to calculate the light field in polar coordinates as follow.
where r 1 , r 2 is vector containing the amplitude and azimuth. This equation can be calculated by convolution and the output coordinate r 2 can be enlarged as Ar 2 .
we set 'm' to change the size of both source plane and image plane.
We substitute it into equation (17) and obtain where Q represent phase factor and E" represent convolution formula.
We define scale coordinate r' 2 as follow and substitute it into equation (21). We obtain It could be taken into the equation (18). Finally, we acquire the formula of output field as: where FFT is the Fourier transform and IFFT is the inverse Fourier transform. Now that we have get the expression of angular-spectrum propagation, we can test the grid spacing η 1 in the source plane and η 2 in the image plane.
From F [f 1 , r 1 ]: Above all, we can control over the observation plane grid spacing η 2 to overcome the fixed spacing grid by FFT. Moreover, adjusting the optical system matrix parameters can operate both the Fourier transform and the Fresnel transform, which analysis of defocus use the Fresnel transform.
The simulation analysis with proposed algorithm is seen as a test run for experiments to get clear results, not the main guidance. By simulation, we can try to find the suitable sample cell size and match it with the parameter of CCD camera existed in market. The proposed algorithm is helpful to know what size the sampling unit is suitable. It has reference value to choose CCD camera for experiment on aspect of resolution.

Simulation of flexible spacing grid for modal decomposition
In this part, we proved by simulation that this algorithm applied in modal decomposition can alter the spacing grid flexibly and improve resolution of output patterns. The simulation background was that the Large Mode Area fiber (LMA fiber) was excited by an Nd:YAG laser (λ 0 = 1064 nm). We assumed the waist radius w 0 was 3 millimeters, and the NA close to 0.01, which the weak guide fiber satisfies. The V = 5.167 parameter clearly establishes that the fiber supports the six lowestorder eigenmodes, namely LP 01 , LP 02 , LP 11e, LP 11o , LP 21e , and LP 21o , at the same instance, each mode exhibits two different polarizations, with o and e representing the specific angular distributions, as shown in figure 3. For simplicity, we reviewed only one polarization state, so as to decompose a known field into its spatial eigenstates.
The second stage of our simulation deployed the MODAN [13] filter for modal decomposition, as well as the Two-step ABCD algorithm to characterize the light field in variable grid coordinates. The simulation system setup, as displayed in figure 1, can be experimentally established via an easy 2 fsetup. The incident light beam is a coherent superposition of six lowest order modes, with similar polarization. The diffraction optical element MODAN, comprising corresponding mode conjugations, can be prepared via a liquid crystal spatial light modulator (SLM) as a digital hologram element. The liquid crystal SLM ensures diffractive efficiency and common applicability, relative to microlithography gratings. MODAN size is assumed as triple that of the optical fiber, and the grid count is N * N. A Fast Fourier transform was realized at the back focal plane, wherein the far-field image was a crosscorrelation of the incident beam U 1 (x 1 , y 1 ) and the eigenmode Ψ (x 1 , y 1 ), recorded via CCD 1 .On the diffraction optical elements (MODAN) rear surface, it indicated a near-field intensity distribution, which is the mix of six lowest-order mode intensities, recorded by the CCD 2 . The comparation of nearfield distribution shown as figures 4(c)-(d). We picked the desired intensity points, u = Vn * k/f, from the far-field distribution, to subsequently calculate other significant parameters. The initial parameters and calculated results were recorded in table 1 and shown as figures 4(a)-(b).
In simulation, the date of modal decomposition is presented in table 1, wherein we assumed beforehand the weights of all six modes were the square root of 1/6, such that the coefficient squares aggregate as 1. We initially simulated mode decomposition via a Fast Fourier transform and results named as initial data in next parts. The modal weights and relative phase are calculated via the Two-step ABCD algorithm subsequently, named as calculated data. The comparison diagrams of the initial data and the calculated data were shown at figures 4(a)-(b). Comparing 2 sets of data, we defined the error rate as '(initial date-calculated data)/initial date'. From table 1 and figure 4, they revealed that the calculated data conformed to the initial data within the margins of error, thereby validating the method as precise, effective and practicable.
Furthermore, the two-step ABCD algorithm accomplished FFT operation, and reconstructed the light field in variable sizes, as shown in figure 5. Both figures involve a 512 * 512 coordinate grid in the spectral domain (f x ,f y ), with variable cell sizes; however, any such adjustments should consider the corresponding image plane dimensions, as well Fourier frequency shifts. The image plane spectrum range, on the other hand, corresponds with the wavelength and detection distance. Hence, the sampling gird overall correlates with the carrier frequency, wavelength, and detection distance. With the Twostep ABCD algorithm, we can adjust the parameter 'm' in the project of Matlab software to change the grid sizes of image plane (x 2 ,y 2 ), in which the grid size was δ 2 = m * δ x . The range of image plane was x 2 = (-N/2:N/2-1) * δ 2 . So was y 2 . According to the relationship between the x 2 and f x in Fraunhofer diffraction, f x = x 2 /(λ * f), we can get the range of f x was (-N/2:N/2-1) * δ 2 /(λ * f). We eventually adjusted the grid size of image plane .
The abovementioned theory and formulas in the part of 2.1 specify that the laser beam is in accord with eigenmode orthogonality only at the point of frequency-shift position. It is worth noting that the pattern of each mode is exactly not the same as the initial six lowest-order modes, and every independently separated light spot is effectively the single-mode self-correlation pattern. Figure 5 illustrates the light far-filed patterns in simulations in the frequency domain (f x , f y ). Figure 5(a) indicates the direct FFT results are too tiny and centralized to discern any significant mode pattern or positions. Figures 5(b)-(c) present simulation of the same content based on the two-step ABCD algorithm. The far-field results in figures 5(b) and (c) exhibit size-tunable light field distributions. The grid cell length and width are identical, in the two-dimensional coordinate system. The positional shift in spatial frequency coordinates (f x , f y ) corresponds to what we initially set. The initial data and calculated date were recorded in table 1 and figure 4. From them, we can know that the modal weights and the relative phase calculated by the algorithm agreed with the initial date with a tiny error. It has proved that the two-step algorithm applied in modal decomposition can observe the light field patterns in variable size and characterize laser beam correctly.

Numerical analysis for effect of defocus on the modal decomposition
As we established, only a precise rear focal plane distribution can realize the Fourier transform. However, engineering  applications define the focal plane based on the smallest circular spot size, whilst shifting the observation plane. General engineering applications do not provide strict operation, or high-precision elements, to calibrate the focal position. Meanwhile, the detection system defocus degenerates the performance such that the light spot would be a diffused facula. In our example, the spectrum width of each mode progressively increased, causing overlaps. These could potentially be averted via large frequency shifts; however, those are difficult to set in practical applications. The errors ensuing from defocused/overlapped modes being inevitable, we theoretically derived formulas to quantify those. The generated expression specifically helps to analyze how defocus influences characterizing, based on the modal decomposition algorithm, which in any case has a positive reference value in engineering applications.
As initially mentioned, select intensities of frequency-shift positions helped depict the image plane. We would like to know how the defocus affect the propagation of modes. It was assumed that the simulation had only a single mode in the incident field, such as LP01, and the frequency domain coordinate system was fx = x 2 /(λ * B) and fy = y 2 /(λ * B). We almost achieved maximum intensity at the frequency-shift position, even if within the defocus range. Hence, in this coordinate, we assumed defocus barely influenced frequencyshifts, and continued selecting intensities on these specific positions as well. The Collins formula helps analyze laser beam propagation in optical systems. We therefore adopted this method, and expressed the output E out as:  where E in was the lens front surface(figure 3 Fourier lens) light field, E out was the light far-field distribution tested by CCD 1 in figure 3. The input E in complex amplitude is the superposition of eigenmodes with different phases.
Taking the equation (27) into the equation (26), the integrated formula is expressed as: Each eigenmode should consider intensity normalized in polar coordinates; taking LP 01 mode as an example, where q 1 is the normalized constant. The optical system is a simple 2-f setup, and its ABCD matrix is specified as, where z is the distance between the lens and the image plane.
The derivative progression refers to the integral formulas listed herewith.
where M = µ1 + µ2 +… + µn, Re(p) > 0, Re(v + M) > 0. ψ 2 is bivariate hypergeometric functions. Generally, LP01 mode is almost referred as basic mode, defined V1 = 0. In the back focal plane, A = 0, and at the point (0,0) in spectrum domain, the output distribution according to the reference integral formula equation (29) is The result reveals that the complex distribution of LP 01 at (0,0) is satisfied with mode normalization. On the other hand, it proved the above theory of modal decomposition in the way of Collins formula. So as to other modes. When the detection plane did not fix the back focal plane(A ̸ = 0), the variable ∆f denoted defocusing, such that ∆f = z-f. According the reference integral formulas equation (36), the output Eout with defocus is , Ψ 2 is bivariate hypergeometric function. The equation (38) describe the change of intensity at frequency-shift positions, via a formula that primarily factors B/A and its power function, wherein A and B are elements in the ABCD matrix; B/A approaches infinity, when defocus approaches zero. We therefore avoided the intensity at the focus, which otherwise would render the simulation result infinite. The above formulas yield the integral result as a sum of infinite polynomials. In mathematic theory, the lengthy hybrid Bessel integration is tuff to acquire analytic solutions. For simplicity, we calculated the numerical integration of the hybrid Bessel integration in simulation. In numerical integration method, the error of results is maybe caused by the limited condition, the simplicity of sum and avoiding focus position.

Simulation with two-ABCD algorithm.
As simulating the effect of defocus on intensity with relative defocus (∆f/f) range from −0.5 to 0.5, we discover that the error of mode weights is especially large based on the intensities data. The larger relative defocus, the lager error rate. Only between −0.01 and 0.01, the data is effective for analysis. During the propagation with defocus, we picked up the intensities at the specific positions and made the LP 01 mode as basic mode. The results of simulation are as follow.
The simulation results revealed similar trends for all modes that at specific position, the defocusing was lesser, and the light intensity was stronger, establishing the figure 6(a) profile. The radios of every max light intensity were nearly close to 1 around ∆f = 0, agreed with the proportion of initial mode weight squared. On account of the carried defocus, the propagation is Fresnel diffraction not Fourier transform, and the output patterns is a superposition of the six modes Fresnel diffraction with different phase factors, such as two subimages in figure 6(a). Obviously, defocus has a significant effect on the relative light intensity, which means the modal decomposition algorithm is strict with the defocus.

Comparison of results from two methods.
We showed the effect of defocusing on modal decomposition detailly in defocusing range of 1% to −1%. The mode weights error percentage was selected as the element to reflect the result of modal decomposition. Figures 6(b)-(g) detailly highlighted error percentage of mode weights profiles for each eigenmode based on two methods. The comparison of the two methods was not only helpful to analyze their respective characteristics, but also obtained more accurate defocusing influence rule. According to results, the trend of intensity for all modes with the increase of relative defocus at specific positions in spectrum domain was same as figure 6(a), while the profiles for low-order modes and high-order modes were not identical.
From the figures, it revealed that both in two methods, at the position of ∆f = 0, back focal plane, the mode decomposition and characterizing can be realized exactly, which were consistent with patterns in figure 5. When there was defocus error in testing distance, the result of mode decomposition would be affected terribly, in which the range keeping precisely correct radio did not appear. the results observed both in two methods helped establish that the defocus significantly influenced the modal characterizing algorithm. Beyond this, the effect of defocusing was variant with the order of mode and worse in high-order mode. From the numerical analysis, it can be discovered that within a small range of relative defocus, from −0.25% to 0.25%, the error of mode weights can almost around 0 for LP 01 , LP 11o , LP 11e . However, as for higher-order mode, the order is higher, the range of keeping error rate around zero is smaller.
The principle of this phenomenon is that the Fourier transformation is changed as Fresnel transformation because of the introduction of defocus, which leads to larger spot width for each mode and adjacent patterns overlap. Consequently, the light intensity at special positions contain one part generated by mode orthogonality and the other part generated by overlap. The intensity data have lost the correctness for calculating the mode weights and relative phase. Furthermore, the fact is that higher order eigenmodes by Fraunhofer diffraction are bottle beam, wherein the transvers image is concentric circles with different light and shade. The power of light is mostly concentrated on the outer ring, not equally in every ring. According this, the light intensities will up and down, maybe which is difficult to find an accurate universal law for defocus effect. But expect for zero point in error percentage of mode weights at ∆f = 0, others are no physical meaning, they are just a coincidence about pure number. Overall, we infer that defocus would have worse effect on modal decomposition algorithm, especially for higher-order modes.   5,0), the results will out of calculation of Matlab software, stopping program running. So, the error can be decreased but not disappearing. Second, we only pick up the intensity of single point, while the size of single point is out of the resolution of CCD camera. Limited by the CCD camera resolution, when we want to measure the light intensity at (1.5,0) in frequency domain, the actual data displayed by the camera is the mean value of the light intensity within a small radius centered on (1.5, 0). Picking up single point intensity does not be realized in engineering application, maybe it contributes to much error of intensity data. The second method, numerical analysis, the mathematical model of eigenmode was Bessel function in my work, whose integral result is not analytic solution. In the reference integral formula (equations (34-36)), it is an infinite integral region. However, software Matlab cannot work out the integral in infinite region. For the feasibility of integral operation, we set a few limited conditions and approximate expressions. So is the sum of infinite polynomials (equation 36). It is reducible to the sum of limited polynomials. Even if the operation process has been simplified much, the running time by Matlab software is too long, which cannot be applied in instant feedback optical system.
Combined with the above reasons, the Two-step ABCD algorithm is more rapid and effective to analyze the defocus effect in the relative defocus range of −0.25% to 0.25%. Beyond this range, defocusing has a terrible effect on modal decomposition, and we also try to minimize defocusing in engineering applications. This range also indicates that finding accurate back focal plane is tuff and vital work in absolute coordinate.

Conclusion
This research has demonstrated that the Two-step ABCD algorithm in simulation can realize modal decomposition and characterization of laser modes based on correlation filer, allowing image plane size modifications. Choosing the suitable grid size needs to consider the carried frequency and the size of light filed for variable optical setup situations. It is possible to observe the light facula more clearly and ensure the correctness and accuracy of light intensities at desired positions. Combining the algorithm and numerical analysis, we analyze the effect of defocus on modal decomposition based on correlation filter and discover that error of mode weights increases with the defocus. If the relative defocus between −0.25% and −0.25%, the correctness of modal decomposition and characterization can still be kept for LP 01 , LP 11o , LP 11e . The higher-order eigenmodes does not have the similar limited band of relative defocus, which effect of defocus is worse than lower-order eigenmodes. This research verifies feasibility and effectiveness of Two-step ABCD algorithm and analyzes effect of defocus by simulation and numerical analysis, with a view to provide a more reliable data and references for the engineer applications later.
In nature, the Bessel beam has the feature of focal shift, which brings more problems to define the back focal plane. In later work, we will continue studying and researching these error problems and solutions to improve accuracy and applicability of modal decomposition based on correlation filter for more engineer situations.