**1. Introduction**

The assessment of the effectiveness of correction for turbulent distortion of laser beams and the enhancement of correction quality are problems that have been studied for more than twenty years and are still actual at present [1–5]. Often, such problems are solved with the application of numerical methods, and in developed models, Zernike polynomials are often used to define the number of degrees of freedom of the active element of an adaptive system. Actually, a phase screen that defines the turbulent distortions of a beam is represented as a series of polynomials [4,5]. Therefore, it seems appropriate to carry out the numerical estimation of this expansion precision, to our knowledge, this problem has not yet been considered thoroughly.

Optical vortices (singular points or beam dislocations) are specific objects that appear sometimes in a wavefront of radiation [6]. Such objects are characterized by a small region of zero amplitude around the vortex and a cut in phase distribution. The cut ends in the center of a vortex. The development of optical vortices was observed in experiments with beams reflected from a rough surface [7,8]. At the same time, the authors of Refs [3,4] demonstrated in numerical [3] and laboratory experiments [4] that a beam propagating in a turbulent atmosphere acquires a complex form, and scintillations and intensity zeros appear in its amplitude distribution, so in the wavefront of such a beam, one can also expect the development of singular points. Beam control under such conditions presents considerable difficulties, for example, the negative influence of dislocations on the adaptive correction was theoretically demonstrated in Ref. [1].

The present article is a continuation of research in these areas. In the text, we discuss two problems. Firstly, in numerical experiments, we demonstrated that singular points appear in a wavefront of a Gaussian beam with an initial phase profile given by Zernike polynomials or by a screen simulating atmospheric turbulence, i.e., by a smooth function.

**Citation:** Kanev, F.; Makenova, N.; Veretekhin, I. Development of Singular Points in a Beam Passed Phase Screen Simulating Atmospheric Turbulence and Precision of Such a Screen Approximation by Zernike Polynomials. *Photonics* **2022**, *9*, 285. https://doi.org/10.3390/ photonics9050285

Received: 15 March 2022 Accepted: 19 April 2022 Published: 21 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

To increase reliability, the number and coordinates of dislocations were found using three detection algorithms built on different principles [9–12]. Thus, the number of vortices in the wavefront as a function of turbulence intensity and the distance traveled by the beam was obtained. so the reconstruction of such distribution by an adaptive mirror with a continuous surface leads to errors in the process of beam control. It is expedient to analyze this problem by numerical methods, but to do so the whole mathematical and numerical model of an adaptive system is required. This model should include many rather complex blocks, so

The second problem, the solution of which is discussed in the article, is the problem of phase screen approximation. It is known, that to realize successfully the adaptive control of a laser beam in a turbulent atmosphere requires the precise reconstruction of a reference beam phase. Often this phase has a complex form and includes singular points,

appear in a wavefront of a Gaussian beam with an initial phase profile given by Zernike polynomials or by a screen simulating atmospheric turbulence, i.e., by a smooth function. To increase reliability, the number and coordinates of dislocations were found using three detection algorithms built on different principles [9–12]. Thus, the number of vortices in the wavefront as a function of turbulence intensity and the distance traveled by

*Photonics* **2022**, *9*, x FOR PEER REVIEW 2 of 15

The second problem, the solution of which is discussed in the article, is the problem of phase screen approximation. It is known, that to realize successfully the adaptive control of a laser beam in a turbulent atmosphere requires the precise reconstruction of a reference beam phase. Often this phase has a complex form and includes singular points, so the reconstruction of such distribution by an adaptive mirror with a continuous surface leads to errors in the process of beam control. It is expedient to analyze this problem by numerical methods, but to do so the whole mathematical and numerical model of an adaptive system is required. This model should include many rather complex blocks, so at the thirst stage, we just considered the precision of a phase screen approximation by Zernike polynomials. The obtained results and characteristics of the method are presented in the current paper. at the thirst stage, we just considered the precision of a phase screen approximation by Zernike polynomials. The obtained results and characteristics of the method are presented in the current paper. **2. Materials and Methods**  The propagation of coherent laser radiation in a turbulent medium was considered on the base of the numerical simulation technique. The optical layout of the numerical experiment is shown in Figure 1. The beam had a Gaussian amplitude profile; phase modulation was realized in the plane of the emitting aperture by a surface formed as a sum of Zernike polynomials [13] or by a phase screen simulating atmospheric turbulence.

#### **2. Materials and Methods** The spectral density of the index of refraction fluctuations on this screen was given by the following equation [14]:

the beam was obtained.

The propagation of coherent laser radiation in a turbulent medium was considered on the base of the numerical simulation technique. The optical layout of the numerical experiment is shown in Figure 1. The beam had a Gaussian amplitude profile; phase modulation was realized in the plane of the emitting aperture by a surface formed as a sum of Zernike polynomials [13] or by a phase screen simulating atmospheric turbulence. The spectral density of the index of refraction fluctuations on this screen was given by the following equation [14]: 2 2 2 11/6 2 2 0 0 ( ) 0.033 ( ) exp( / ), 2 / , 5.92 / *n nL m L m C Ll* <sup>−</sup> Φ κ = κ + κ −κ κ κ = π κ = . (1) Equation (1) describes the von Karman spectrum of fluctuations, *L*0, *l*0 are the outer and inner scales of turbulence, *Сn* is the structure constant of atmospheric turbulence, related with Fried's coherence length by the formula 2 2 3/5 <sup>0</sup> 1.68( ) , *<sup>n</sup> r CkL* <sup>−</sup> = here *k* is the wave number, *L* is a path length. This formula demonstrates that 0*r* depends on the in-

$$\Phi\_n(\kappa) = 0.033 \text{C}\_n^2 (\kappa\_L^2 + \kappa^2)^{-11/6} \exp(-\kappa^2/\kappa\_m^2), \quad \kappa\_L = 2\pi/L\_0, \kappa\_m = 5.92/l\_0. \tag{1}$$

$$\Phi\_n^{\text{ext}} = \begin{cases} 0 & \text{if } n = 1 \\ 1 & \text{if } n = 2 \end{cases}$$

**Figure 1.** Layout of numerical experiments. **Figure 1.** Layout of numerical experiments.

From the screen to the plane of observation the beam propagated under conditions of free diffraction. Propagation was described by the wave equation [15]: <sup>1</sup> 2 , *gr E E ik E zv t* <sup>⊥</sup> ∂ ∂ + =Δ ∂ ∂ (2) Equation (1) describes the von Karman spectrum of fluctuations, *L*0, *l*<sup>0</sup> are the outer and inner scales of turbulence, *C<sup>n</sup>* is the structure constant of atmospheric turbulence, related with Fried's coherence length by the formula *r*<sup>0</sup> = 1.68(*C* 2 *n k* <sup>2</sup>*L*) −3/5 , here *k* is the wave number, *L* is a path length. This formula demonstrates that *r*<sup>0</sup> depends on the intensity of the index of refraction fluctuations, wavelength, and on the distance passed by a beam.

and the fast Fourier transform was employed to solve it [16]. In Equation (2), *x* and *y* are From the screen to the plane of observation the beam propagated under conditions of free diffraction. Propagation was described by the wave equation [15]:

$$2ik\left(\frac{\partial E}{\partial z} + \frac{1}{v\_{gr}}\frac{\partial E}{\partial t}\right) = \Delta\_\perp E\_\prime \tag{2}$$

and the fast Fourier transform was employed to solve it [16]. In Equation (2), *x* and *y* are coordinates in the plane normal to the direction of propagation; *z* is a coordinate along this direction; ∆<sup>⊥</sup> = *∂* <sup>2</sup>/*∂x* <sup>2</sup>+ *∂* <sup>2</sup>/*∂y* 2 is the Laplace operator; *v*gr is the group velocity of the beam. The coordinates *x* and *y,* and parameter *r*<sup>0</sup> were normalized to the beam initial radius *a*0, and coordinate *z* to the diffraction length *z<sup>d</sup>* = *ka*<sup>2</sup> 0 . All computer applications used in investigations of singular beam propagation were developed by the authors of the paper with Visual C++ language [17].

The following set of parameters was used to characterize the optical field in the plane of observation: Power-in-the-bucket (PIB).

$$J(t) = \frac{1}{P\_0} \iint \rho(\mathbf{x}, y) I(\mathbf{x}, y, t) d\mathbf{x} dy. \tag{3}$$

This parameter is proportional to the beam power incident in an aperture of radius *S<sup>t</sup>* . In Equation (3) *P*<sup>0</sup> is the total power of the beam and

$$\rho(\mathfrak{x}, \mathfrak{y}) = \exp[- (\mathfrak{x}^2 + \mathfrak{y}^2) / S\_t^2]$$

is an aperture function.

The shift of the beam gravity center along axis *x*:

$$X\_{\mathfrak{C}} = \frac{1}{P\_0 a\_0} \iint \mathbf{x} I(\mathbf{x}, y, t) d\mathbf{x} dy,\tag{4}$$

and *y*:

$$Y\_{\varepsilon} = \frac{1}{P\_{0}a\_{0}} \iint y I(\mathbf{x}, y, t) d\mathbf{x} dy. \tag{5}$$

The effective radius of the beam

$$R\_{Eff} = \left\{ \frac{1}{P\_0 a\_0} \iint \left( \mathbf{r}\_\perp - \mathbf{r}\_c \right)^2 I(x, y, t) dx dy \right\}^{1/2}. \tag{6}$$

To obtain complete information concerning beam distortions, in addition to the abovedescribed functions, we also registered the quantity and coordinates of optical vortices developed in the wavefront. The detailed description of the algorithms constructed to localized singular points of the wavefront was given in Ref. [9]. Algorithms were built with the use of the following special features of vortex radiation:


$$\Gamma(\mathfrak{a}) = \oint\_{P} \mathfrak{a}(\mathfrak{r}, z) d\mathfrak{r} \tag{7}$$

is equal to ±2*πn* if an optical vortex falls in an integration contour. The second algorithm was based on this property. Here *n* is a vortex topological charge and *P* is the perimeter of the integration contour.

3. A vortex is a point of intersections of isolines. Isolines should be drawn in distributions of real and imaginary parts of beam complex amplitude, magnitudes of corresponding functions were equal to zero along the line. On this property, the third algorithm was developed.

According to the assessments presented in Ref. [9], the highest precision of registration was achieved with algorithms 2 and 3. Computer processing of interference patterns (algorithm No. 1) does not provide high resolution, so only a small number of vortices can be localized with the application of this method.

#### **3. Obtained Results**

In the present paragraph, the results characterizing the development of singular points in a wavefront of radiation are presented. The phase profile of radiation has been formed as a sum of Zernike polynomials [14], or specified by a phase screen, simulating atmospheric turbulence. Furthermore we illustrated the quality of approximation of such a screen.

Changes in Gaussian beam amplitude are shown in Figures 2–5, while its initial phase was formed by a single polynomial. In the pictures, we can see the phase distribution a screen.

a screen.

a screen.

of radiation (Figure 2a), cross-sections of this distribution (Figure 2b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figure 2c). tion of radiation (Figure 2a), cross-sections of this distribution (Figure 2b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figure 2c). tion of radiation (Figure 2a), cross-sections of this distribution (Figure 2b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figure 2c). tion of radiation (Figure 2a), cross-sections of this distribution (Figure 2b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figure 2c).

Changes in Gaussian beam amplitude are shown in Figures 2–5, while its initial phase was formed by a single polynomial. In the pictures, we can see the phase distribu-

Changes in Gaussian beam amplitude are shown in Figures 2–5, while its initial phase was formed by a single polynomial. In the pictures, we can see the phase distribu-

Changes in Gaussian beam amplitude are shown in Figures 2–5, while its initial phase was formed by a single polynomial. In the pictures, we can see the phase distribu-

In the present paragraph, the results characterizing the development of singular points in a wavefront of radiation are presented. The phase profile of radiation has been formed as a sum of Zernike polynomials [14], or specified by a phase screen, simulating atmospheric turbulence. Furthermore we illustrated the quality of approximation of such

In the present paragraph, the results characterizing the development of singular points in a wavefront of radiation are presented. The phase profile of radiation has been formed as a sum of Zernike polynomials [14], or specified by a phase screen, simulating atmospheric turbulence. Furthermore we illustrated the quality of approximation of such

In the present paragraph, the results characterizing the development of singular points in a wavefront of radiation are presented. The phase profile of radiation has been formed as a sum of Zernike polynomials [14], or specified by a phase screen, simulating atmospheric turbulence. Furthermore we illustrated the quality of approximation of such

(**b**)

*Photonics* **2022**, *9*, x FOR PEER REVIEW 4 of 15

*Photonics* **2022**, *9*, x FOR PEER REVIEW 4 of 15

*Photonics* **2022**, *9*, x FOR PEER REVIEW 4 of 15

**Figure 2.** Phase specification by polynomial No. 5 (third-order astigmatism). Phase surface (**a**) and corresponding amplitude distribution of the beam passed path of 0.1 diffraction length (**c**) are shown by grayscale images. Cross-sections of the phase surface are presented in picture (**b**). Cross-section 1 was cut through the center of the beam, cross-section 2 corresponds to the cut shifted on the beam radius from the center in a positive direction of the *y*-axis, and number 3—to the cut shifted in a negative direction. **Figure 2.** Phase specification by polynomial No. 5 (third-order astigmatism). Phase surface (**a**) and corresponding amplitude distribution of the beam passed path of 0.1 diffraction length (**c**) are shown by grayscale images. Cross-sections of the phase surface are presented in picture (**b**). Cross-section 1 was cut through the center of the beam, cross-section 2 corresponds to the cut shifted on the beam radius from the center in a positive direction of the *y*-axis, and number 3—to the cut shifted in a negative direction. **Figure 2.** Phase specification by polynomial No. 5 (third-order astigmatism). Phase surface (**a**) and corresponding amplitude distribution of the beam passed path of 0.1 diffraction length (**c**) are shown by grayscale images. Cross-sections of the phase surface are presented in picture (**b**). Cross-section 1 was cut through the center of the beam, cross-section 2 corresponds to the cut shifted on the beam radius from the center in a positive direction of the *y*-axis, and number 3—to the cut shifted in a negative direction. **Figure 2.** Phase specification by polynomial No. 5 (third-order astigmatism). Phase surface (**a**) and corresponding amplitude distribution of the beam passed path of 0.1 diffraction length (**c**) are shown by grayscale images. Cross-sections of the phase surface are presented in picture (**b**). Cross-section 1 was cut through the center of the beam, cross-section 2 corresponds to the cut shifted on the beam radius from the center in a positive direction of the *y*-axis, and number 3—to the cut shifted in a negative direction.

**Figure 3.** The same as in Figure 2, but the phase was specified by polynomial No. 9 (trefoil). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). **Figure 3.** The same as in Figure 2, but the phase was specified by polynomial No. 9 (trefoil). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). **Figure 3.** The same as in Figure 2, but the phase was specified by polynomial No. 9 (trefoil). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). **Figure 3.** The same as in Figure 2, but the phase was specified by polynomial No. 9 (trefoil). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**).

**Figure 4.** The same as in Figure 2, the phase was specified by polynomial No. 7 (coma). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). **Figure 4.** The same as in Figure 2, the phase was specified by polynomial No. 7 (coma). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). **Figure 4.** The same as in Figure 2, the phase was specified by polynomial No. 7 (coma). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). **Figure 4.** The same as in Figure 2, the phase was specified by polynomial No. 7 (coma). Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. Three cross-sections coincide (**b**). *Photonics* **2022**, *9*, x FOR PEER REVIEW 5 of 15

**Figure 5.** The same as in Figure 2, the phase was specified by polynomial No. 20. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. **Figure 5.** The same as in Figure 2, the phase was specified by polynomial No. 20. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown.

In the pictures, we can see the phase distribution of radiation (Figures 2a–5a),

Usually, an approximation of a screen simulating atmospheric turbulence is performed by the sum of polynomials [14], so it is feasible to consider influence induced not only by discrete components but also by the total sum of components. Corresponding examples are shown in Figures 6 and 7; in Figure 6, summation was realized from the first to sixth (trefoil) polynomials and in Figure 7 to the seventh (coma). All coefficients

(**а**) (**b**) (**c**)

**Figure 6.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 6. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**)

**Figure 7.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 7. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**)

Registered distributions of singular points are presented in Figures 8 and 9, corresponding phase profiles were set by such polynomials as coma (Figure 8) and trefoil (Figure 9). To improve the reliability of results dislocations were localized by three algo-

(**b**) (**c**)

were the same and equal to one.

are shown.

are shown.

(**а**)

rithms built on different principles.

In the pictures, we can see the phase distribution of radiation (Figures 2a, 3a, 4a and 5a), cross-sections of this distribution (Figures 2b, 3b, 4b and 5b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figures 2c, 3c, 4c and 5c). In the pictures, we can see the phase distribution of radiation (Figures 2a–5a), cross-sections of this distribution (Figures 2b–5b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figures 2c–5c). cross-sections of this distribution (Figures 2b–5b), and amplitude distribution of a beam passed a distance of 0.1 diffraction length (Figures 2c–5c). Usually, an approximation of a screen simulating atmospheric turbulence is per-

**Figure 5.** The same as in Figure 2, the phase was specified by polynomial No. 20. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown.

In the pictures, we can see the phase distribution of radiation (Figures 2a–5a),

(**а**) (**b**) (**c**)

**Figure 5.** The same as in Figure 2, the phase was specified by polynomial No. 20. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown.

(**а**) (**b**) (**c**)

*Photonics* **2022**, *9*, x FOR PEER REVIEW 5 of 15

*Photonics* **2022**, *9*, x FOR PEER REVIEW 5 of 15

Usually, an approximation of a screen simulating atmospheric turbulence is performed by the sum of polynomials [14], so it is feasible to consider influence induced not only by discrete components but also by the total sum of components. Corresponding examples are shown in Figures 6 and 7; in Figure 6, summation was realized from the first to sixth (trefoil) polynomials and in Figure 7 to the seventh (coma). All coefficients were the same and equal to one. Usually, an approximation of a screen simulating atmospheric turbulence is performed by the sum of polynomials [14], so it is feasible to consider influence induced not only by discrete components but also by the total sum of components. Corresponding examples are shown in Figures 6 and 7; in Figure 6, summation was realized from the first to sixth (trefoil) polynomials and in Figure 7 to the seventh (coma). All coefficients were the same and equal to one. formed by the sum of polynomials [14], so it is feasible to consider influence induced not only by discrete components but also by the total sum of components. Corresponding examples are shown in Figures 6 and 7; in Figure 6, summation was realized from the first to sixth (trefoil) polynomials and in Figure 7 to the seventh (coma). All coefficients were the same and equal to one.

**Figure 6.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 6. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. **Figure 6.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 6. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. **Figure 6.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 6. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown.

(**b**) (**c**) **Figure 7.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 7. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) **Figure 7.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 7. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown. **Figure 7.** The same as in Figure 2, the phase was specified by the sum of polynomials from 1 to 7. Phase surface (**a**), cross-sections of phase surface (**b**) and corresponding amplitude distribution (**c**) are shown.

Registered distributions of singular points are presented in Figures 8 and 9, corresponding phase profiles were set by such polynomials as coma (Figure 8) and trefoil (Figure 9). To improve the reliability of results dislocations were localized by three algo-Registered distributions of singular points are presented in Figures 8 and 9, corresponding phase profiles were set by such polynomials as coma (Figure 8) and trefoil (Figure 9). To improve the reliability of results dislocations were localized by three algorithms built on different principles. Registered distributions of singular points are presented in Figures 8 and 9, corresponding phase profiles were set by such polynomials as coma (Figure 8) and trefoil (Figure 9). To improve the reliability of results dislocations were localized by three algorithms built on different principles. *Photonics* **2022**, *9*, x FOR PEER REVIEW 6 of 15

**Figure 8.** Distributions of singular points in a beam with a phase profile set by polynomial No. 8 (coma). The dislocations are shown by circles with «+» or «−» signs. To register optical vortices three algorithms were used. Firstly, vortices were localized by processing distribution of wavefront gradients (**a**); secondly, as a point of intersection of isolines (**b**); and with application of interferometric algorithm (**c**). The number of the found singular points *N*dsl = 98 (**а**), *N*dsl = 92 (**b**) и *N*dsl = 10 (**c**). The normalized path length was equal to 0.1. **Figure 8.** Distributions of singular points in a beam with a phase profile set by polynomial No. 8 (coma). The dislocations are shown by circles with «+» or «−» signs. To register optical vortices three algorithms were used. Firstly, vortices were localized by processing distribution of wavefront gradients (**a**); secondly, as a point of intersection of isolines (**b**); and with application of interferometric algorithm (**c**). The number of the found singular points *N*dsl = 98 (**a**), *N*dsl = 92 (**b**) *N*dsl = 10 (**c**). The normalized path length was equal to 0.1.

(**а**) (**b**) (**c**)

**Figure 9.** The same as in Figure 8, but the phase was specified by polynomial No. 9 (trefoil). *N*dsl =

**Figure 10.** Number of singular points (*N*dsl) in a beam wavefront with a phase specified by different polynomials (*N*z). Registration was carried out in a region with a radius equal to the initial radius of

The number of the found vortices as a function of the distance traveled by a beam is presented in Figures 11 and 12. The initial phase was formed by such polynomials as coma and trefoil (Figure 11) or set by a phase screen representing atmospheric turbulence

a beam (curve 1), to two initial radii (2), to three initial radii (3). *Z* = 0.1.

The dependence of dislocation quantity on polynomial number is shown in Figure 10. The radius of a region where detection was performed changed in the range from 1 to

are shown.

rithms built on different principles.

58 (**а**), *N*dsl = 43 (**b**) и *N*dsl = 7 (**c**).

1.4 initial radii of a beam.

(Equation (1)).

(**a**) (**b**) (**c**)

(**a**) (**b**) (**c**)

*Photonics* **2022**, *9*, x FOR PEER REVIEW 6 of 15

**Figure 8.** Distributions of singular points in a beam with a phase profile set by polynomial No. 8 (coma). The dislocations are shown by circles with «+» or «−» signs. To register optical vortices three algorithms were used. Firstly, vortices were localized by processing distribution of wavefront gradients (**a**); secondly, as a point of intersection of isolines (**b**); and with application of interferometric algorithm (**c**). The number of the found singular points *N*dsl = 98 (**а**), *N*dsl = 92 (**b**) и *N*dsl = 10

**Figure 8.** Distributions of singular points in a beam with a phase profile set by polynomial No. 8 (coma). The dislocations are shown by circles with «+» or «−» signs. To register optical vortices three algorithms were used. Firstly, vortices were localized by processing distribution of wavefront gradients (**a**); secondly, as a point of intersection of isolines (**b**); and with application of interferometric algorithm (**c**). The number of the found singular points *N*dsl = 98 (**а**), *N*dsl = 92 (**b**) и *N*dsl = 10

**Figure 9.** The same as in Figure 8, but the phase was specified by polynomial No. 9 (trefoil). *N*dsl = 58 (**а**), *N*dsl = 43 (**b**) и *N*dsl = 7 (**c**). **Figure 9.** The same as in Figure 8, but the phase was specified by polynomial No. 9 (trefoil). *N*dsl = 58 (**a**), *N*dsl = 43 (**b**) *N*dsl = 7 (**c**). **Figure 9.** The same as in Figure 8, but the phase was specified by polynomial No. 9 (trefoil). *N*dsl = 58 (**а**), *N*dsl = 43 (**b**) и *N*dsl = 7 (**c**).

The dependence of dislocation quantity on polynomial number is shown in Figure 10. The radius of a region where detection was performed changed in the range from 1 to 1.4 initial radii of a beam. The dependence of dislocation quantity on polynomial number is shown in Figure 10. The radius of a region where detection was performed changed in the range from 1 to 1.4 initial radii of a beam. The dependence of dislocation quantity on polynomial number is shown in Figure 10. The radius of a region where detection was performed changed in the range from 1 to 1.4 initial radii of a beam.

(**c**). The normalized path length was equal to 0.1.

(**c**). The normalized path length was equal to 0.1.

The number of the found vortices as a function of the distance traveled by a beam is presented in Figures 11 and 12. The initial phase was formed by such polynomials as coma and trefoil (Figure 11) or set by a phase screen representing atmospheric turbulence The number of the found vortices as a function of the distance traveled by a beam is presented in Figures 11 and 12. The initial phase was formed by such polynomials as coma and trefoil (Figure 11) or set by a phase screen representing atmospheric turbulence (Equation (1)). The number of the found vortices as a function of the distance traveled by a beam is presented in Figures 11 and 12. The initial phase was formed by such polynomials as coma and trefoil (Figure 11) or set by a phase screen representing atmospheric turbulence (Equation (1)). *Photonics* **2022**, *9*, x FOR PEER REVIEW 7 of 15

**Figure 11.** Dependence of dislocation number on distance *Z* passed by the beam. Calculation was performed for a beam with phase specified by two comas (polynomial numbers and number of curves are 7 and 8) and by trefoils (polynomials and curves are 6 and 9). The radius of registration region is 1.4 of the initial radius of the beam. **Figure 11.** Dependence of dislocation number on distance *Z* passed by the beam. Calculation was performed for a beam with phase specified by two comas (polynomial numbers and number of curves are 7 and 8) and by trefoils (polynomials and curves are 6 and 9). The radius of registration region is 1.4 of the initial radius of the beam.

**Figure 12.** Dependence of dislocation number on distance *Z* passed by the beam. Calculation was performed for several realizations of random phase screens, numbers of which are printed near

The results presented above demonstrate that the development of optical vortices is possible in a beam with an initially smooth phase profile. Let us now consider the possibility of such a profile approximation by Zernike polynomials. In the first of the two problems, the phase was represented by the sum of polynomials and approximated by the sum of polynomials, in the second, the phase was set by a turbulent screen and approximated by polynomials. In the process of approximation calculation of polynomial

The numbers of polynomials entering the basis of approximation are larger than that

The screen and approximation basis include the same number of functions (Tables

The number of polynomials in the basis is lesser than in the phase screen (Table 5). The initial phase profile and the surface obtained as a result of approximation were formed by the same functions, so to assess the precision of calculations we can compare coefficients of the members entering the sum and deviation of one surface from another. Corresponding parameters (coefficients *C*ZN, here, *N* is a polynomial number) were put in Tables 1–5. The difference between the two surfaces was calculated according to the

coefficients was realized with the least-mean-square method [15].

in the phase screen (Table 1; 12 and 9 polynomials correspondingly).

The first problem we divided into three variants:

curves.

2–4).

formula

region is 1.4 of the initial radius of the beam.

curves.

**Figure 12.** Dependence of dislocation number on distance *Z* passed by the beam. Calculation was performed for several realizations of random phase screens, numbers of which are printed near **Figure 12.** Dependence of dislocation number on distance *Z* passed by the beam. Calculation was performed for several realizations of random phase screens, numbers of which are printed near curves.

**Figure 11.** Dependence of dislocation number on distance *Z* passed by the beam. Calculation was performed for a beam with phase specified by two comas (polynomial numbers and number of curves are 7 and 8) and by trefoils (polynomials and curves are 6 and 9). The radius of registration

The results presented above demonstrate that the development of optical vortices is possible in a beam with an initially smooth phase profile. Let us now consider the possibility of such a profile approximation by Zernike polynomials. In the first of the two problems, the phase was represented by the sum of polynomials and approximated by the sum of polynomials, in the second, the phase was set by a turbulent screen and approximated by polynomials. In the process of approximation calculation of polynomial The results presented above demonstrate that the development of optical vortices is possible in a beam with an initially smooth phase profile. Let us now consider the possibility of such a profile approximation by Zernike polynomials. In the first of the two problems, the phase was represented by the sum of polynomials and approximated by the sum of polynomials, in the second, the phase was set by a turbulent screen and approximated by polynomials. In the process of approximation calculation of polynomial coefficients was realized with the least-mean-square method [15].

coefficients was realized with the least-mean-square method [15]. The first problem we divided into three variants:

The first problem we divided into three variants: The numbers of polynomials entering the basis of approximation are larger than that The numbers of polynomials entering the basis of approximation are larger than that in the phase screen (Table 1; 12 and 9 polynomials correspondingly).

in the phase screen (Table 1; 12 and 9 polynomials correspondingly). The screen and approximation basis include the same number of functions (Tables 2–4). **Table 1.** The screen was set by 9 polynomials and approximated by 12 polynomials. Dimensions of computational grid are 256 × 256 nodes.


The screen and approximation basis include the same number of functions (Tables 2–4).

**Table 2.** The screen was set by 9 polynomials and approximated by 9 polynomials. Dimensions of computational grid are 256 × 256 nodes.


**Table 3.** The screen was set by 12 polynomials and approximated by 12 polynomials. Dimensions of computational grid are 256 × 256 nodes.


**Table 4.** The screen was set by 12 polynomials and approximated by 12 polynomials. Dimensions of computational grid are 2048 × 2048 nodes.


The number of polynomials in the basis is lesser than in the phase screen (Table 5).

**Table 5.** The screen was set by 9 polynomials and approximated by 8 polynomials. Polynomial coefficient No. 9 was reduced to 0.2. Dimensions of computational grid are 1024 × 1024 nodes.


The initial phase profile and the surface obtained as a result of approximation were formed by the same functions, so to assess the precision of calculations we can compare coefficients of the members entering the sum and deviation of one surface from another. Corresponding parameters (coefficients *C*ZN, here, *N* is a polynomial number) were put in Tables 1–5. The difference between the two surfaces was calculated according to the formula

$$\varepsilon\_{\rm Ph} = \sum\_{i,j}^{M} \sqrt{\left(A\_{ij} - B\_{ij}\right)^2} \Big/ \sum\_{i,j}^{M} \sqrt{\left(A\_{ij}\right)^2} \tag{8}$$

here, *Aij* is the value of the given function in node *i, j*; *Bij* is the value of the function obtained as a result of approximation; *M* × *M* are dimensions of the computational grid.

To increase the consistency of the solution, we compared the parameters of two beams with Gaussian amplitude distributions. The phase of the first was set by the screen formed by the sum of polynomials (Figure 1), while for the phase of the second, the surface taken obtained was in the approximation procedure. Both beams passed the same distance, after that their parameters were calculated and compared. Calculation of the following

Number of the column

Values of parameters corresponding to the given phase

Values obtained as a result of approximation

Number of the

Values of parameters corresponding to the given phase

Values obtained as a result of approximation

Number of the column

Values of parameters corresponding to the given phase

Values obtained as a result of approximation

> parameters was performed: power-in-the-bucket *J* (Equation (3)), shifts *X*<sup>C</sup> and *Y*<sup>C</sup> of beam gravity center along OX and OY axes (Equations (4) and (5)), the effective radius of the beam *R*Eff (Equation (6)). Furthermore, we calculated the difference between amplitude distributions. To do so, Equation (8) was used, but instead of phase profiles, we substituted into the formula corresponding amplitude distributions. **Parameters** *C***Z1** *C***Z2** *C***Z3** *C***Z4** *C***Z5 …** *C***Z8** *J* **εPh ε<sup>A</sup>** *R***Eff** *X***c** *Y***<sup>c</sup>** 1 2 3 4 5 6 7 8 9 10 11 12 13 1.00 1.00 1.00 1.00 1.00 … 1.00 0.52 0.00 0.00 2.35 0.56 −0.41

*Photonics* **2022**, *9*, x FOR PEER REVIEW 9 of 15

**Parameters** *C***Z1** *C***Z2** *C***Z3** *C***Z4** *C***Z5 …** *C***Z12** *J* **εPh ε<sup>A</sup>** *R***Eff** *X***c** *Y***<sup>c</sup>**

column 1 2 3 4 5 6 7 8 9 10 11 12 13

**Parameters** *C***Z1** *C***Z2** *C***Z3** *C***Z4** *C***Z5 …** *C***Z12** *J* **εPh ε<sup>A</sup>** *R***Eff** *X***c** *Y***<sup>c</sup>**

of computational grid are 256 × 256 nodes.

of computational grid are 2048 × 2048 nodes.

**Table 3.** The screen was set by 12 polynomials and approximated by 12 polynomials. Dimensions

1.00 1.00 1.50 1.00 0.50 … 1.00 0.52 0.00 0.00 2.35 0.07 −0.07

19.27 −3.9 −1452.1 1.00 −1413.1 … 1.00 0.18 4.10 0.69 3.21 0.44 0.39

1 2 3 4 5 6 7 8 9 10 11 12 13

1.00 1.00 1.50 1.00 0.50 … 1.00 0.53 0.00 0.00 2.25 0.06 −0.07

0.95 0.77 952.22 1.00 960.12 … 1.00 0.54 0.22 0.14 2.24 0.06 −0.07

The third variant is illustrated by the results given in Table 5 and in Figure 13.

**Table 5.** The screen was set by 9 polynomials and approximated by 8 polynomials. Polynomial coefficient No. 9 was reduced to 0.2. Dimensions of computational grid are 1024 × 1024 nodes.

**Table 4.** The screen was set by 12 polynomials and approximated by 12 polynomials. Dimensions

Results obtained in the first variant are presented in Table 1. The basis of approximation was formed by 12 polynomials, and 9 polynomials were included in the phase screen.

In the second variant, an equal number of polynomials were included in the basis and screen. In Table 2, the date is presented for the screen formed by 9, and in Tables 3 and 4 by 12 polynomials. 46.61 1.00 17.83 1.00 17.83 … 1.00 0.13 8.19 0.69 3.47 1.95 −0.36

The third variant is illustrated by the results given in Table 5 and in Figure 13.

**Figure 13.** A beam amplitude registered at the distance of 0.1 diffraction length (**a**) and its initial phase (**b**) represented as the sum of the first nine polynomials. All polynomial coefficients except No. 9 were equal to one; the ninth coefficient was equal to 0.2. Amplitude (**c**) and phase (**d**) profiles obtained as a result of approximation with the use of eight polynomials. **Figure 13.** A beam amplitude registered at the distance of 0.1 diffraction length (**a**) and its initial phase (**b**) represented as the sum of the first nine polynomials. All polynomial coefficients except No. 9 were equal to one; the ninth coefficient was equal to 0.2. Amplitude (**c**) and phase (**d**) profiles obtained as a result of approximation with the use of eight polynomials.

The next problem is the approximation of the phase screen simulating turbulent fluctuations of the index of refraction. As in the previous example, propagation was analyzed by two beams with Gaussian amplitude distributions. The phase of the first was set by the phase screen, while the phase of the second was obtained as a result of approximation. Parameters of beams calculated in the plane of observations were presented in Tables 6–8; in Figures 14–16, we can see amplitude distributions and phase profiles of the beams.

**Table 6.** Characteristics of the beam in the plane of observations. Parameters of numeric experiment were the same as in the previous example (Figure 14).


**Table 7.** Characteristics of the beam in the plane of observations calculated with increased turbulence intensity. Parameters of numeric experiment were the same as in the previous example (Figure 15).


Values of parameters corre-

Values obtained as a result of approximation

Values obtained as a result of approximation

Values of parameters corre-

Values obtained as a result of approximation

Values of parameters corre-

Values obtained as a result of approximation


the beams.

*Photonics* **2022**, *9*, x FOR PEER REVIEW 10 of 15

**Table 8.** Characteristics of the beam in the plane of observations. Parameters of numeric experiment were the same as in the previous example (Figure 16). *Photonics* **2022**, *9*, x FOR PEER REVIEW 10 of 15

The next problem is the approximation of the phase screen simulating turbulent fluctuations of the index of refraction. As in the previous example, propagation was analyzed by two beams with Gaussian amplitude distributions. The phase of the first was set by the phase screen, while the phase of the second was obtained as a result of approximation. Parameters of beams calculated in the plane of observations were presented in Tables 6–8; in Figures 14–16, we can see amplitude distributions and phase profiles of

(**d**) with initial phase presented in Figure (**c**). Dimensions of computation grid are 1024 × 1024, *r*0 =

**Figure 14.** Precision of phase screen approximation. Phase distribution specified by the screen simulating atmospheric turbulence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**). Dimensions of computation grid are 1024 × 1024, *r*0 = 0.2. **Figure 14.** Precision of phase screen approximation. Phase distribution specified by the screen simulating atmospheric turbulence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**). Dimensions of computation grid are 1024 × 1024, *r*<sup>0</sup> = 0.2. 0.208 0.291 0.986 0.043 0.019 0.077 0.075 0.107

0.208 0.291 0.986 0.043 0.019 0.077 0.075 0.107 **Figure 15.** Increase in turbulence intensity (*r*0 = 0.1). Other parameters of numeric experiment are the same as in Figure 14. Phase distribution specified by the screen simulating atmospheric turbulence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**). **Figure 15.** Increase in turbulence intensity (*r*<sup>0</sup> = 0.1). Other parameters of numeric experiment are the same as in Figure 14. Phase distribution specified by the screen simulating atmospheric turbulence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**). sponding to the given phase 0 0 0.983 0.076 0.035 0.038 0.092 0.100 0.208 0.544 0.980 0.076 0.035 0.085 0.078 0.114

lence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**). **Table 7.** Characteristics of the beam in the plane of observations calculated with increased turbulence intensity. Parameters of numeric experiment were the same as in the previous example (Fig-**Figure 16.** Results obtained with computational grid of larger (2048 × 2048 nodes) dimensions. Other parameters are the same as in the previous example (Figure 15, Table 6). Phase distribution specified by the screen simulating atmospheric turbulence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**). **Figure 16.** Results obtained with computational grid of larger (2048 × 2048 nodes) dimensions. Other parameters are the same as in the previous example (Figure 15, Table 6). Phase distribution specified by the screen simulating atmospheric turbulence (**a**), amplitude distribution of a beam with initial phase profile specified by this screen (**b**), surface obtained as a result of approximation (**c**), amplitude of a beam (**d**) with initial phase presented in Figure (**c**).

**Table 8.** Characteristics of the beam in the plane of observations. Parameters of numeric experi-

0.244 0.675 0.955 0.076 −0.135 0.093 0.122 0.153

*4.1. Scintillations of Amplitude and Emergence of Singular Points in the Wavefront of a Beam* 

The first block of results included in paragraph 3 of the current paper illustrates the fact that a singular wavefront can be obtained if the initial phase profile of a beam is set to be Zernike polynomials. To prove this fact, the beam with a given phase passed some distance in a non-aberrating medium, its amplitude and phase profiles were calculated in the plane of observations (the optical layout of the corresponding numerical experiment is shown in Figure 1), and singular points were registered by three different algorithms. Figures 2–7 show that points of zero intensity appear in such a beam, the localization of optical vortices is illustrated in Figures 8 and 9. If a phase profile is formed by a sum of polynomials with equal coefficients, amplitude distribution is influenced mainly by the last polynomial entering the sum. This can be seen from the comparison of Figures 3 and 6. In the first case, the phase was set by trefoil, and in the second by the sum including polynomials from the first to trefoil. All polynomial coefficients were the same and equal

Calculations were performed with different intensities of atmospheric distortions and with the use of different calculation grids. In all numeric experiments, the inner scale of turbulence was larger than the diameter of the beam. Even under such conditions, the quality of the screen approximation was not high, but with a smaller inner scale, the ob-

ment were the same as in the previous example (Figure 16).

tained results were completely unsatisfying.

*with a Phase Set by Zernike Polynomials* 

**Parameters εPh εAmp** *J Sh***<sup>X</sup>** *Sh***<sup>Y</sup>** *REff***<sup>X</sup>** *REff***<sup>Y</sup>** *REff*

sponding to 0 0 0.965 0.076 −0.149 0.039 0.081 0.090

ure 15).

**4. Discussion** 

Calculations were performed with different intensities of atmospheric distortions and with the use of different calculation grids. In all numeric experiments, the inner scale of turbulence was larger than the diameter of the beam. Even under such conditions, the quality of the screen approximation was not high, but with a smaller inner scale, the obtained results were completely unsatisfying.

#### **4. Discussion**
