A Hybrid Model for Analysis of Laser Beam Distortions Using Monte Carlo and Shack–Hartmann Techniques: Numerical Study and Experimental Results

: The hybrid model for analyzing distortions of a laser beam passed through a moderately scattering medium with the number of scattering events up to 10 is developed and investigated. The model implemented the Monte Carlo technique to simulate the beam propagation through a scattering layer, a ray-tracing technique to propagate the scattered beam to the measurements plane, and the Shack–Hartmann technique to calculate the scattered laser beam distortions. The results obtained from the developed model were conﬁrmed during the laboratory experiment. Both the numerical model and laboratory experiment showed that with an increase of the concentration value of scattering particles in the range from 10 5 to 10 6 mm − 3 , the amplitude of distortions of laser beam propagated through the layer of the scattering medium increases exponentially.


Introduction
Currently, optical diagnostic techniques are essential experimental tools employed in a large and diverse set of research and industrial applications.Based on measurements of scattering intensity and light extinction, a number of these techniques are able to determine the parameters of scattering particles, such as size, shape, and concentration [1].Examples of applications are numerous: in meteorology, the measurement of atmospheric constituents is performed by analyzing the backscattering light signal via the LIDAR (Light Detection and Ranging) technique [2].In oceanography, marine picoplankton is characterized from the laser light scattered by the illuminated particles [3].The optical radiation is actively used in various bio-medical applications as a non-invasive diagnostic tool, as well as to treat a number of diseases [4].During the past three decades, the investigation of laser light propagation in skin tissue and the human brain has been particularly extensive [4][5][6][7][8][9][10].In combustion engineering, a variety of laser-based diagnostics have been developed over more than three decades in order to determine the physical properties of fuel droplets from atomizing sprays [6].These optical measurements are of fundamental importance in the improvement of combustion efficiency and in the reduction of pollutant emissions from modern internal combustion engines and gas turbines.Another field where optical diagnostic techniques are of particular significance is free-space transmission of energy or information, namely, free-space optical communications, wireless energy transfer through the scattering and turbulent atmosphere, imaging objects through scattering atmospheric aerosol, space debris removal, etc.
One of the limitations of the optical techniques is related to the radiation attenuation and scattering.Under optically thick conditions, the scattering can significantly decrease the quality of the beam and thus introduce errors in measurements.However, it is necessary not only to simulate the propagation of the radiation through the medium but also to establish techniques to measure the beam distortions introduced by that medium.
Measuring laser beam distortions is one of the existing problems that has to be solved.Another problem that is also of great importance is compensating measured distortions.As is known, in a turbid medium, a part of a laser beam energy is absorbed while the rest is scattered.This scattered light makes objects look blurred and degrades beamfocusing efficiency.There are a wide variety of applications for which the solution of the problem is very important, i.e., wireless data transmission, pattern recognition, and medical noninvasive diagnosis, among others [11,12].The wavefront-shaping technique to focus the radiation that propagates through a strongly scattering layer was developed back in 2007 [13].The authors used the SLM and were able to properly shape the wavefront to focus the radiation through scattering objects.Various techniques have been developed since 2007 to resolve the focusing of light through scattering media [14,15].Authors of [16] showed the ability to obtain clear images of objects by means of holographic techniques that reverse a scattering process.Authors of [17] use coherence and spatial gating to produce images of optical properties of a tissue up to 9-mm deep with millimeter-scale resolution via a technique known as multispectral multiple-scattering low-coherence interferometry.Imaging of fluorescent objects when the space behind a scattering medium is limited can be achieved using the technique described in [18].SLM can be used to obtain images of objects located within or behind a scattering layer [19][20][21][22].In [23], authors demonstrated a steady-state focusing of coherent light through dynamic scattering media.In [24], authors demonstrated the use of the method of inversion of transmission matrix.They proved that the focus' peak-to-background intensity ratio can be higher compared to conventional methods.In [25], the authors proposed a new technique that allowed for the 204-fold increase of the peak-to-noise ratio.Another group of scientists proposed the system that uses the full-polarization optical digital-phase conjugation and doubles the peak-to-noise ratio [26].In [27], programmable acoustic optic deflectors were used to obtain the images of the targets in the living tissue.Bimorph deformable mirrors, as well as stacked-actuator mirrors, are also used to increase the radiation energy in the far field while focusing a laser beam through a scattering medium [28][29][30], turbulent medium [31,32], or a pinhole [33] for free-space optical communication systems.Spatial-light modulators are widely used in different areas of research; their use is not limited to the imaging and focusing applications only.Spatial-light modulators also show vast potential in such applications as micro-and nano-scale fabrication [34], and laser beam shaping [35,36].
In this particular paper, we will not concentrate on the problem of compensating distortions.Instead, we present the hybrid model for (1) The simulation of laser beam propagation through the moderately scattering medium using the Monte Carlo technique; and (2) Measuring distortions obtained by the laser beam by means of the numerical model of the Shack-Hartmann sensor.The proposed model was briefly mentioned in the authors' previous paper [37], which was devoted to the compensation of measured laser beam distortions using bimorph deformable mirrors.In this paper, we would like to describe the proposed model in detail and provide the results of the numerical analysis, as well as the experimental validation of the model.
The research is organized as follows.First, a radiative transfer equation and its possible solutions are presented.Second, a Monte Carlo technique to solve the radiative transfer equation is described; its software implementation and verification are provided.Third, the Shack-Hartmann principle is explained, and the developed numerical model is described.Fourth, the laboratory experimental setup with the Shack-Hartmann sensor is shown.A comparison and discussion of numerical and experimental results obtained are presented.

Radiative Transfer Equation
Estimation of radiation attenuation and multiple scattering effects requires knowledge of the radiative transfer of optical energy through the medium [1].In the field of optical diagnostics, the radiative transfer equation is commonly used to describe the photons transport, specifying a balance of optical energy distributed between the incident, outgoing, absorbing, and scattering radiation propagated through the medium.The radiative transfer equation is presented below: where I is a radiance in the particular direction (either s or s'), t is time, → r is the vector position, → s is the incident direction of propagation, µ e is an extinction coefficient [38,39], and µ s is a scattering coefficient.This is a value inverse proportional to the distance on which the collimated light flux will be reduced due to the scattering [40,41].f → s , → s is the scattering phase function derived from the appropriate scattering theory (e.g., Mie or Rayleigh theory), dΩ is the solid angle spanning → s , and c is the speed of the light in the surrounding medium.
In other words, Equation ( 1) can be explained as follows: The resultant radiance along the direction → s is a sum of the scattered radiance coming from all directions → s to the direction → s and the radiance loss due to the absorption.The total extinction represented by the first right-side term of Equation ( 1) describes two processes that contribute to the total resultant radiance, namely the radiance loss due to the scattering process and the radiance loss due to the absorption process.Though the radiative transfer equation can be applied for a wide variety of media, the analytical solution of this equation is available for rather simple circumstances with different assumptions.Thus, techniques such as the method of path integrals, diffusion approximation, light scattering by Brownian particles [42], and the method of the small-angle approximation [43] were developed.
Nevertheless, the techniques mentioned above are supposed to be used under different assumptions and simplifications.For example, the methods of light scattering by Brownian particles, as well as the diffusion approximation, only suit the small anisotropy factor-scattering media.The method of small-angle approximation does not consider the diffusively scattered light and cannot be applied to our task because the idea of the research is to analyze the influence of both quasi-ballistic and diffusive components of the scattered radiation.The method of path integrals has two variations: analytical and stochastic.Analytical method is basically the modified method of Brownian particles, whereas stochastic method is Monte Carlo method.As is known, Monte Carlo technique [44,45] is the most versatile and widely used numerical technique that provides an approximate solution of the radiative transfer equation for the vast majority of conditions.

Monte Carlo Simulation
The Monte Carlo technique is a widely used approach toward sampling probability density functions for simulating a wide range of problems.The first use of the method for photon transport in biological materials was Adams and Wilson in 1983, which considered isotropic scattering [46].Keijzer et.al. introduced anisotropic scattering into the Monte Carlo simulation of biological tissues, implementing a simulation that propagated photons using cylindrical coordinates, which introduced the Hop/Drop/Spin nomenclature for organizing the program [47].Prahl et al. reformulated the program using photon propagation based on Cartesian coordinates, which made the program much simpler to convey in written form [48]. Wang and Jacques adapted and augmented the work of Keizer and Prahl to write the Monte Carlo Multi-Layered (MCML) program that considers tissues with many planar layers with different optical properties [8].
Monte Carlo technique simulates the behavior of the elementary parts of a physical system.In regards to light propagation through a medium, this technique uses quantum nature of light and simulates the propagation of photon flux [49].
In general, there are two Monte Carlo techniques: corpuscular and wave [50].Both methods have their advantages and disadvantages.
The "Wave Monte Carlo" treats the propagation of radiation as a wave propagation.It describes the propagation of a wave through a set of screens that contains scattering particles.The resultant wave becomes stochastic due to the interference between the unscattered and scattered waves.As a result, we get the complex amplitude of the random field.An ensemble of the obtained fields for the statistically independent sets of screens is the calculated intensity distribution.This technique is particularly applied for the monodirection radiation.
In contrast, the "Corpuscular Monte Carlo" treats the propagation of radiation as a corpuscula (photon) flux.This technique performs the statistical calculation of the ensemble of photons-the average number of photons depending on the coordinate equals to the intensity distribution."Corpuscular Monte Carlo" is applied for both the weakly and strongly scattering media.Another feature of this method is that it is not supposed to model the fixed position of the scatterers.It instead models the events of the photon-particle interaction with the desired probability density.From this point of view, the Brownian movement of scatterers is considered in the simulation by default.
Choosing each method is a matter of opinion, since both show similar results [50].We have chosen the "Corpuscular Monte Carlo" technique for our simulation.
First, we simulated the initial uniform intensity distribution of a laser beam with diameter equal to 4 mm and wavelength equal to 0.65 µm.The random number generator was used to calculate the initial position of each of the 2.5 × 10 11 photons within the beam aperture.Afterward, we calculated the trajectory of propagation of each photon through the 5-mm layer of scattering medium consisted of a monodisperse spheres of 1-µm diameter and refractive index of 1.582 (refractive index of the medium was 1.33).The influence of relative refractive index is rather complicated [51].For small particles, the influence on the forward-scattering signals is mainly due to the part of the internal reflection if the relative refractive index deviates from 1.However, when the relative refractive index is close to 1, the effects on the forward-scattered light from both the surface reflection and internal reflection are great.For large particles, the contributions of the surface reflection and the internal reflection to the forward-scattered light are much weaker than the diffraction when the relative refractive index deviates from 1.When the relative refractive index is very close to 1, the effects on the forward-scattered light from the internal reflection are great.
The absorption coefficient for such conditions (i.e., for the selected wavelength and refractive index) will be negligible according to [52].Thus, no absorption will appear, and we do not have to track whether a particular photon should be excluded from the simulation.The concentration value, or a number of spheres per cubic millimeter, of the scattering suspension was varied approximately from 10 5 to 10 6 mm −3 .
In order to define the trajectory of the photon, three key parameters have to be calculated.First, free path length l-the distance between two scattering events or between two sequential interactions between the photon and a particle (Equation ( 2)).Second, scattering angle θ-the angle between the current and the new direction of photon's propagation (Equation ( 3)).Third, azimuthal angle ϕ-the angle between the projection of the new direction on the plane perpendicular to the initial direction of the photon's propagation (Equation ( 4)).
where ξ l , ξ θ , ξ ϕ -random numbers uniformly distributed within the interval [0, 1), µ s -scattering coefficient, g-anisotropy factor.The direction of photon after a scattering event is defined by scattering diagram or phase function [53,54].The phase function is a function of the probability density of scattering in the direction s' of a photon moving in the direction s, i.e., a function that characterizes the elementary act of scattering.In other words, the phase function describes the angular distribution of the intensity of light scattered by an elementary volume.The shape of the scattering diagram depends on the wavelength of the incident radiation, the size of the scatterer, and the relative refractive index [55].The phase function is also characterized by the mean cosine of the scattering angle g, or the so-called anisotropy factor [56].
Anisotropy factor g shows the degree of isotropy of scattering process.If particle scatters light equally in all directions, then g = 0; if a particle scatters more energy in the forward direction, then g > 0; if a particle scatters more energy in the backward direction, then g < 0. Comprehensive study on the dependance of anisotropy factor on the radiation wavelength and scattering particle diameter can be found in [57][58][59][60].
Scattering coefficient µ s is one of the key parameters of a scattering medium.It is measured in inverse millimeters and characterizes the scattering strength per unit length.The scattering coefficient is calculated using Mie theory [8].
After calculation of free-path length, scattering and azimuthal angles of the new position of the photon is defined.If this position was inside the scattering medium, the calculation of the next set of l, θ, ϕ is repeated.Once the photon reached the border of the medium, its parameters (e.g., total travelled optical path length, last coordinates, last direction of propagation, total number of encountered scattering events, etc.) were saved, and the photon was traced to the measurements plane (emulation of camera sensor).The procedure is repeated for each photon.When trajectories of all photons were calculated, we obtained the intensity distribution at the measurements' plane.Figure 1 shows the results of the simulation of laser beam propagation through the scattering medium with three different scatterers' concentration values-from 1.5 × 10 5 mm −3 to 7.5 × 10 5 mm −3 .
Algorithms 2023, 16, x FOR PEER REVIEW 5 of 17 where  ,  ,  -random numbers uniformly distributed within the interval [0, 1), μsscattering coefficient, g-anisotropy factor.The direction of photon after a scattering event is defined by scattering diagram or phase function [53,54].The phase function is a function of the probability density of scattering in the direction s' of a photon moving in the direction s, i.e., a function that characterizes the elementary act of scattering.In other words, the phase function describes the angular distribution of the intensity of light scattered by an elementary volume.The shape of the scattering diagram depends on the wavelength of the incident radiation, the size of the scatterer, and the relative refractive index [55].The phase function is also characterized by the mean cosine of the scattering angle g, or the so-called anisotropy factor [56].
Anisotropy factor g shows the degree of isotropy of scattering process.If particle scatters light equally in all directions, then g = 0; if a particle scatters more energy in the forward direction, then g > 0; if a particle scatters more energy in the backward direction, then g < 0. Comprehensive study on the dependance of anisotropy factor on the radiation wavelength and scattering particle diameter can be found in [57][58][59][60].
Scattering coefficient μs is one of the key parameters of a scattering medium.It is measured in inverse millimeters and characterizes the scattering strength per unit length.The scattering coefficient is calculated using Mie theory [8].
After calculation of free-path length, scattering and azimuthal angles of the new position of the photon is defined.If this position was inside the scattering medium, the calculation of the next set of l, θ, φ is repeated.Once the photon reached the border of the medium, its parameters (e.g., total travelled optical path length, last coordinates, last direction of propagation, total number of encountered scattering events, etc.) were saved, and the photon was traced to the measurements plane (emulation of camera sensor).The procedure is repeated for each photon.When trajectories of all photons were calculated, we obtained the intensity distribution at the measurements' plane.Figure 1 shows the results of the simulation of laser beam propagation through the scattering medium with three different scatterers' concentration values-from 1.5 × 10 5 mm −3 to 7.5 × 10 5 mm −3 .In order to verify our implementation of Monte Carlo technique, we compared the results of our simulation with the results obtained from "MCML" model written by In order to verify our implementation of Monte Carlo technique, we compared the results of our simulation with the results obtained from "MCML" model written by Lihong Wang and Steven Jacques [8].We considered the following parameters of simulation: halfinfinite 10-mm thick layer, 1-µm diameter spheres, 1.5 × 10 5 mm −3 concentration value, 0.9 anisotropy factor, 0.2 mm −1 scattering coefficient, 6.8-mm initial laser beam diameter, 10 × 10-mm measurement plane size.
We used three comparison criteria: 1. Number of photons transmitted through and reflected from the medium.2.

3.
Intensity distributions at the exit edge of the medium (in Table 1, we presented percentage of photons located within specified regions).

Shack-Hartmann Technique
Shack-Hartmann sensor [61][62][63][64][65][66] is well-known device that is widely used in large and diverse sets of applications, primarily to measure distortions of the wavefront of the radiation passed through different media-turbulent and/or scattering atmosphere, biological tissues, etc. Since, in our work, we consider the case of moderately scattering medium where average number of scattering events is from 0 to 10-the so-called crossover scattering regime [67,68]-we assumed (and proved it experimentally) that the Shack-Hartmann technique still could be applied to measure the scattered laser beam distortions.We will cover this point later in Section 2. 4.
In a conventional Shack-Hartmann technique, the wavefront of the incident radiation is divided into a number of sub-apertures by means of lenslet array.Lenslet array is a thin, flat base with a grid of micro lenses (lenslets) etched on the base.Each lenslet commonly has diameter from 100 to 300 µm and focal length f from 3 mm to 8 mm.Radiation passes through these lenslets and forms a set of focal spots at the measurement (sensor) plane (Figure 2a).
Since the diameter of each lenslet is rather small, it is assumed that the wavefront W is flat within a single lenslet, and the only aberration it has is tip-tilt.In case a tip-tilt equals zero (i.e., a wavefront is flat and parallel to the plane of the lenslet), radiation is focused to the center of the corresponding sub-aperture of the sensor.In case a wavefront has non-zero tip-tilt within a lenslet, then the focal spot will be displaced (S x , S y ) from the center of the sub-aperture proportional to the tip-tilt value (Figure 2b).In other words, if we measure these displacements S x and S y of the focal spot per X and Y axis, correspondingly, we will obtain the values of partial derivatives ∂W/∂x and ∂W/∂y of the wavefront W within each sub-aperture (Equation ( 5)): where  Since the diameter of each lenslet is rather small, it is assumed that the wavefront W is flat within a single lenslet, and the only aberration it has is tip-tilt.In case a tip-tilt equals zero (i.e., a wavefront is flat and parallel to the plane of the lenslet), radiation is focused to the center of the corresponding sub-aperture of the sensor.In case a wavefront has nonzero tip-tilt within a lenslet, then the focal spot will be displaced (Sx, Sy) from the center of the sub-aperture proportional to the tip-tilt value (Figure 2b).In other words, if we measure these displacements Sx and Sy of the focal spot per X and Y axis, correspondingly, we will obtain the values of partial derivatives ∂W/∂x and ∂W/∂y of the wavefront W within each sub-aperture (Equation ( 5)): where , , , -partial derivatives of the wavefront W, f-focal spot of each lenslet,  ,  -displacements of a focal spot within a sub-aperture per X and Y axis.
On the other hand, in order to analytically describe and visualize the wavefront surface, one can use the polynomial approximation; for example, B-Splines [69] or Zernike polynomials [70][71][72][73], which are commonly used in optics.These polynomials are orthogonal within the unity circle and are analytically expressed in polar coordinates as follows: Thus, the wavefront W can be presented in a form of Zernike polynomial expansion: where  -Zernike coefficient that represents the aberration value, N-number of the Zernike polynomials used, l = n − 2m.
Knowing the analytical representation of polynomials, we can calculate the partial derivatives ∂W/∂x and ∂W/∂y of a wavefront W: On the other hand, in order to analytically describe and visualize the wavefront surface, one can use the polynomial approximation; for example, B-Splines [69] or Zernike polynomials [70][71][72][73], which are commonly used in optics.These polynomials are orthogonal within the unity circle and are analytically expressed in polar coordinates as follows: Thus, the wavefront W can be presented in a form of Zernike polynomial expansion: where a nm -Zernike coefficient that represents the aberration value, N-number of the Zernike polynomials used, l = n − 2m.
Knowing the analytical representation of polynomials, we can calculate the partial derivatives ∂W/∂x and ∂W/∂y of a wavefront W: Thus, the wavefront partial derivatives ∂W/∂x and ∂W/∂y can be analytically defined using Zernike polynomials (Equation ( 8)).However, they can also be calculated using the measured displacements S x and S y of focal spots on the Shack-Hartmann sensor: Finally, we determine the overdetermined system of linear equations with the unknown coefficients a i .Solving this least square problem [74], we obtain the coefficients a i .From here on, the wavefront can be analytically described and analyzed.

Hybrid Model Implementation
Since we implemented the simulation of a laser beam propagation through a scattering layer using Monte Carlo technique, and since we know how to estimate distortions of laser radiation using Shack-Hartmann technique, we can combine these techniques and numerically estimate what distortions (in terms of focal spots displacements) the beam will obtain after passing the scattering medium.
To solve this, we upgraded our beam propagation model by adding the implementation of Shack-Hartmann technique.In order to do so, we placed the measurements plane, as well as the plane that simulates the lenslet array, at the specified distance from the output edge of the scattering layer (this distance can be varied by the user).Since Monte Carlo technique treats the light as a huge number of photons (or photon packets), we can use the ray-tracing technique and principles of geometrical optics to calculate the trajectory of propagation of photon in free space after it leaves the scattering layer.For this, we calculate the direction cosines dirCosX f , dirCosY f , and dirCosZ f of each photon that left the scattering layer: where dirCosX i , dirCosY i , and dirCosZ i -direction cosines of the photon's propagation trajectory before it left the scattering medium (before the last photon-particle interaction).After calculation of direction cosines, we traced the photon to the plane of the lenslet array.The photon fell on the particular lenslet that retraced it to the measurement plane according to the incidence angle.Photon's coordinates x f , y f , z f were calculated using the formulas: where f -focal length of the lenslet.As a result, we obtain a set of coordinates of all photons that fell on the lenslet array and then retraced to the measurement plane.Afterward, we discretize the measurements plane so it becomes more like the real camera sensor-with fixed pixel size and pixel resolution.The focal spot field also called hartmannogram image obtained on the 4.8 × 4.8-mm measurements plane after the simulation of 4-mm laser beam propagation through the scattering medium is presented on Figure 3.
It can be observed from Figure 3 that there are extra focal spots that appear out of the diameter of 4 mm (the diameter of the initial collimated laser beam that fell on the scattering layer).This effect can be explained by laser beam broadening due to light scattering.A scattered photon flux contains three components [75]: ballistic, on-axis (or "snake"); and off-axis (or diffusive) photons.Ballistic photons travel through a turbid medium without interaction with the scatterers and do not change the initial trajectory.This coherent component of the scattered light is the most important for imaging and focusing applications.On-axis photons undergo few scattering events and travel in near-forward paths along a trajectory that is close to the initial direction of the beam propagation.These photons play an important role for imaging and focusing when the thickness of the scattering medium layer increases because the number of ballistic photons decreases exponentially in this case [11].Off-axis photons scatter in all directions and form a noncoherent component of the scattered light.As noted above, ballistic photons are most valuable for focusing, but their number decreases exponentially when the layer thickness or scatterer concentration increases.When the concentration value of the scattering medium increases, the number of on-axis and diffusive photons increases as well.In addition, all of these photons together form the extra focal spots located between yellow and red circle on Figure 3.It can be observed from Figure 3 that there are extra focal spots that appear out of the diameter of 4 mm (the diameter of the initial collimated laser beam that fell on the scattering layer).This effect can be explained by laser beam broadening due to light scattering.A scattered photon flux contains three components [75]: ballistic, on-axis (or "snake"); and off-axis (or diffusive) photons.Ballistic photons travel through a turbid medium without interaction with the scatterers and do not change the initial trajectory.This coherent component of the scattered light is the most important for imaging and focusing applications.On-axis photons undergo few scattering events and travel in near-forward paths along a trajectory that is close to the initial direction of the beam propagation.These photons play an important role for imaging and focusing when the thickness of the scattering medium layer increases because the number of ballistic photons decreases exponentially in this case [11].Off-axis photons scatter in all directions and form a noncoherent component of the scattered light.As noted above, ballistic photons are most valuable for focusing, but their number decreases exponentially when the layer thickness or scatterer concentration increases.When the concentration value of the scattering medium increases, the number of on-axis and diffusive photons increases as well.In addition, all of these photons together form the extra focal spots located between yellow and red circle on Figure 3.
At this point, we could apply the Shack-Hartmann technique to estimate the distortions of scattered laser radiation.We calculated the displacements of focal spots Sx and Sy (Equation ( 5))and solved the system of linear equations (Equation ( 9)) to find the unknown coefficients  of Zernike polynomials.
During the numerical estimations, we varied the concentration value of the scattering medium from 1.3 × 10 5 to 10 6 mm −3 .We simulated the propagation of the photon flux through the scattering medium with eight particular concentration values from this range.Figure 4 shows some of hartmannogram images obtained on Shack-Hartmann during simulation of laser beam propagation through a scattering medium with different concentrations.At this point, we could apply the Shack-Hartmann technique to estimate the distortions of scattered laser radiation.We calculated the displacements of focal spots S x and S y (Equation ( 5)) and solved the system of linear equations (Equation ( 9)) to find the unknown coefficients a i of Zernike polynomials.
During the numerical estimations, we varied the concentration value of the scattering medium from 1.3 × 10 5 to 10 6 mm −3 .We simulated the propagation of the photon flux through the scattering medium with eight particular concentration values from this range.Figure 4 shows some of hartmannogram images obtained on Shack-Hartmann during simulation of laser beam propagation through a scattering medium with different concentrations.
It can be observed that the spots within the initial aperture of the beam that equaled 4 mm are presented as delta functions centered in the center of the corresponding subapertures-this is because the number of ballistic photons is much bigger than the others.The peripheral focal spots outside of the initial diameter are much wider in diameter, but they also have the clear-cut center.This is due to the impact of rather big portion of on-axis photons (Figure 4a-c).However, with the increase of the concentration of scatterers, the intensity distribution within each peripheral focal spot became uneven without the strongly marked center (Figure 4d-f).This is due to the impact of diffusive photons.
As is known, conventional Shack-Hartmann sensor allows for the measurement of the wavefront of radiation.However, when a laser beam passes through a scattering medium, it does not have the wavefront in strict physical sense since part of the radiation is scattered.It can be observed that the spots within the initial aperture of the beam that equaled 4 mm are presented as delta functions centered in the center of the corresponding subapertures-this is because the number of ballistic photons is much bigger than the others.The peripheral focal spots outside of the initial diameter are much wider in diameter, but they also have the clear-cut center.This is due to the impact of rather big portion of onaxis photons (Figure 4a-c).However, with the increase of the concentration of scatterers, the intensity distribution within each peripheral focal spot became uneven without the strongly marked center (Figure 4d-f).This is due to the impact of diffusive photons.
As is known, conventional Shack-Hartmann sensor allows for the measurement of the wavefront of radiation.However, when a laser beam passes through a scattering medium, it does not have the wavefront in strict physical sense since part of the radiation is scattered.Figure 5 explains what we actually measure with a Shack-Hartmann sensor in this work.It can be observed that the spots within the initial aperture of the beam that equaled 4 mm are presented as delta functions centered in the center of the corresponding subapertures-this is because the number of ballistic photons is much bigger than the others.The peripheral focal spots outside of the initial diameter are much wider in diameter, but they also have the clear-cut center.This is due to the impact of rather big portion of onaxis photons (Figure 4a-c).However, with the increase of the concentration of scatterers, the intensity distribution within each peripheral focal spot became uneven without the strongly marked center (Figure 4d-f).This is due to the impact of diffusive photons.
As is known, conventional Shack-Hartmann sensor allows for the measurement of the wavefront of radiation.However, when a laser beam passes through a scattering medium, it does not have the wavefront in strict physical sense since part of the radiation is scattered.Figure 5 explains what we actually measure with a Shack-Hartmann sensor in this work.As a result, a focal spot is formed at the focal plane of the lenslet.In addition, the displacement of this focal spot from the optical axis of the lenslet is proportional to the local slope of a so-called "average wavefront", which is composed of a set of, generally speaking, independent wavefronts corresponding to different point sources (i.e., scatterers in the medium).Moreover, we measure and analyze this average wavefront.
Due to the axial symmetry nature of the Mie scattering process, only symmetrical distortions of the average wavefront were obtained in the simulation.The obtained coefficients of symmetrical Zernike polynomials (#3, #8, #15, #24) are presented on Figure 6.
of, generally speaking, independent wavefronts corresponding to different point sources (i.e., scatterers in the medium).Moreover, we measure and analyze this average wavefront.
Due to the axial symmetry nature of the Mie scattering process, only symmetrical distortions of the average wavefront were obtained in the simulation.The obtained coefficients of symmetrical Zernike polynomials (#3, #8, #15, #24) are presented on Figure 6.It can be clearly seen from Figure 6 that the amplitude of distortions (values of Zernike coefficients) increased with the increase of the scatterers' concentration.The reason for that is explained by the decrease of ballistic photons and, at the same time, the increase of non-ballistic photons.Finally, it led to the displacements of the focal spot centers from its reference positions.

Hybrid Model Verification: Experimental Results and Discussion
In order to verify the developed hybrid model of the estimation of distortions of a laser beam that passed through a scattering layer, we assembled the experimental setup.The photo of the laboratory setup is presented in Figure 7.It can be clearly seen from Figure 6 that the amplitude of distortions (values of Zernike coefficients) increased with the increase of the scatterers' concentration.The reason for that is explained by the decrease of ballistic photons and, at the same time, the increase of non-ballistic photons.Finally, it led to the displacements of the focal spot centers from its reference positions.

Hybrid Model Verification: Experimental Results and Discussion
In order to verify the developed hybrid model of the estimation of distortions of a laser beam that passed through a scattering layer, we assembled the experimental setup.The photo of the laboratory setup is presented in Figure 7. Laser radiation comes from the fiber-coupled diode laser and propagates through th collimator.The collimated laser beam of 4 mm in diameter then passes through the trans parent glass cuvette with the thickness of 5 mm and dimensions of 18 × 18 mm, filled wit the scattering suspension of 1 μm-diameter polystyrene microbeads produced by Mag sphere, Inc. (Pasadena, CA, USA) [76] with the refractive index equaled to 1.582 [52] di luted in distilled water (refractive index is 1.33).The beam resulting from the cuvette i registered on the Shack-Hartmann sensor [77].According to the specification sheet on th polystyrene microbeads, the diameter variation was within 10%.The microbeads have smooth surface [78] so that no complicated shape-scattering diagram from a single mi Laser radiation comes from the fiber-coupled diode laser and propagates through the collimator.The collimated laser beam of 4 mm in diameter then passes through the transparent glass cuvette with the thickness of 5 mm and dimensions of 18 × 18 mm, filled with the scattering suspension of 1 µm-diameter polystyrene microbeads produced by Magsphere, Inc. (Pasadena, CA, USA) [76] with the refractive index equaled to 1.582 [52] diluted in distilled water (refractive index is 1.33).The beam resulting from the cuvette is registered on the Shack-Hartmann sensor [77].According to the specification sheet on the polystyrene microbeads, the diameter variation was within 10%.The microbeads have a smooth surface [78] so that no complicated shape-scattering diagram from a single microbead is produced.Shack-Hartmann sensor consisted of a digital CCD camera Basler A302fs (1/2-inch sensor with the size of the receiving area equaled to 6.4 × 4.8 mm, camera frame rate is 30 Hz) and a lenslet array (focal length-6 mm, diameter of a lenslet-150 µm, total number of lenslets-greater than 1300).The distance between the Shack-Hartmann sensor and the cuvette was set to 450 mm in order to decrease the impact of the totally diffused photons on the quality of the measurements.Since we made the experiments in a very short period after the preparation of the suspension of scattering particles, the concentration during the experimental measurements was not varied; or, at least, its variation has no impact on the measurements.Moreover, for all of our measurements, we averaged hartmannograms in order to increase the accuracy.For each experiment, we prepared the new suspension.If the cuvette with the initial high-concentration solution was not used for a long period of time, we used the ultrasonic bath to split the coagulated microbeads.
We analyze the average wavefront of scattered light in the circular aperture of 4.8 mm.The aperture's center and the cameras' center coincide.The initial laser beam diameter was 4 mm, whereas the aperture diameter was 4.8 mm.It was necessary to analyze how the non-ballistic photons impact on the distortions of scattered light.The average wavefront measured by the Shack-Hartmann sensor was approximated by Zernike polynomials, as it was done in the simulation.
Figure 8 shows the comparison of the dependence of distortion amplitude (amplitude of distortions of the averaged wavefront) on the concentration values of scatterers.Both the model and experimental curves demonstrated similar trends of increasing the amplitude of distortions while increasing the concentration values of the scatterers.Figure 8 shows a similar trend of increasing distortions of the average wavefront of the scattered laser beam with increasing concentrations of scatterers in the medium in the range from 10 5 to 10 6 mm −3 , both for the model and laboratory experiment.There are a few reasons why the model curve differs from the experimental one.The explanation is schematically represented on Figure 9. Figure 8 shows a similar trend of increasing distortions of the average wavefront of the scattered laser beam with increasing concentrations of scatterers in the medium in the range from 10 5 to 10 6 mm −3 , both for the model and laboratory experiment.There are a few reasons why the model curve differs from the experimental one.The explanation is schematically represented on Figure 9.
for the simulation (red curve) and for the experiment (purple curve).Image adapted from [37].
Figure 8 shows a similar trend of increasing distortions of the average wavefront of the scattered laser beam with increasing concentrations of scatterers in the medium in the range from 10 5 to 10 6 mm −3 , both for the model and laboratory experiment.There are a few reasons why the model curve differs from the experimental one.The explanation is schematically represented on Figure 9. First, in the model (Figure 9a), we did not reproduce the initial distribution of the intensity of the real laser source used in the experiment (Figure 9b).In fact, a generator of uniformly distributed random numbers is used to generate the initial distribution of First, in the model (Figure 9a), we did not reproduce the initial distribution of the intensity of the real laser source used in the experiment (Figure 9b).In fact, a generator of uniformly distributed random numbers is used to generate the initial distribution of photons in the model.Second, the model does not take into account a diffraction phenomenon (Figure 9c), which takes place when the real laser beam is propagating through the optical path (Figure 9d).In the experiment, as a result of diffraction, the focal spots on the Shack-Hartmann sensor in the central part of the beam are displaced from the centers of corresponding sub-apertures at high scatterer concentrations, which leads to additional defocusing.This effect was not observed in the model.For this reason, the experimental curve on Figure 8 increases faster than the model curve at high concentrations.Third, in addition to the centrally symmetric aberrations in the average wavefront, there were also other distortions (in particular, coma aberration) in the experiment (Figure 9f), which could not exist in the model (Figure 9e) in principle due to the symmetry of the nature of Mie scattering on spherical particles.This was due to the fact that the walls of the glass cuvette were not plane-parallel.We have also examined this problem in detail; the results can be found in [79].We made the measurements and determined that the cuvette walls were not absolutely parallel.The deviation was not huge, but in conjugation with the impact of the scattering medium, it led to an overall tilt and, finally, coma aberration.Thus, an additional slope was introduced into the laser beam and, as a result, the coma aberration.
It can be observed from Figure 8 that there is one more data point in the model curve (for a concentration of the order of 10 6 mm −3 ).This is due to the fact that when a laser beam propagates through a scattering medium with a particles' concentration of the order of 10 6 mm −3 , only fractions of a percent of the initial radiation energy level reach the sensor aperture.The traditional, most popular, and widespread CCD cameras have insufficient quantum efficiency and a limited exposure range to detect radiation of such a low power.Regarding the numerical simulation, we could use any concentration values, since with a statistically sufficient number of simulated photons, it is possible to obtain an image of focal spot field.We displayed this data point to show that the model curve does not move far away from the experimental curve even with an increase in the concentration of scatterers.
Thus, the data obtained in the developed model and, subsequently, in a laboratory experiment showed that with an increase in the concentration of scattering particles in the range from 10 5 to 10 6 mm −3 , the amplitude of distortions of the average wavefront of radiation propagating through the layer of the scattering medium increases.

Conclusions
In this paper, we presented the numerical model and algorithm of measurements of distortions of the laser beam that passed through a moderately scattering medium.The simulation of beam propagation through a scattering medium was implemented using the Monte Carlo simulation, the scattered beam propagation in free space from the edge of scattering medium to the lenslet array of the Shack-Hartmann sensor was implemented using ray-tracing techniques, and the measurements of beam distortions was implemented using the Shack-Hartmann technique.We measured the amplitude of distortions of the scattered beam by estimating the displacements of focal spots on the Shack-Hartmann sensor.Both the model and experiment demonstrated similar trends of exponential increase of the amplitude of distortions with the increase of concentration values of the scattering particles in the medium.

Figure 1 .
Figure 1.Results obtained from a numerical simulation.Intensity distributions and cross-sections of laser beam propagated through the scattering medium with particular concentration values.Here, FWHM (full width at half maximum) shows the broadening of the beam due to scattering.

Figure 1 .
Figure 1.Results obtained from a numerical simulation.Intensity distributions and cross-sections of laser beam propagated through the scattering medium with particular concentration values.Here, FWHM (full width at half maximum) shows the broadening of the beam due to scattering.
(x,y) ∂x , ∂W(x,y) ∂y -partial derivatives of the wavefront W, f -focal spot of each lenslet, S x , S y -displacements of a focal spot within a sub-aperture per X and Y axis.

Algorithms 2023 , 17 Figure 2 .
Figure 2. (a) Set of focal spots (focal spot field, hartmannogram) formed by a lenslet array at the sensor plane of a Shack-Hartmann sensor; (b) Principle of calculation of wavefront derivatives by means of measurements of focal-spot displacement.Image adapted from [35].

Figure 2 .
Figure 2. (a) Set of focal spots (focal spot field, hartmannogram) formed by a lenslet array at the sensor plane of a Shack-Hartmann sensor; (b) Principle of calculation of wavefront derivatives by means of measurements of focal-spot displacement.Image adapted from [35].

Algorithms 2023 , 17 Figure 3 .
Figure 3. Shack-Hartmann focal spot field (hartmannogram) obtained after the beam propagation through a scattering layer.The inner yellow circle corresponded to the initial diameter (4 mm) of the collimated beam before entering the scattering volume.The outer red circle corresponded to the area (4.8 mm diameter) of the Shack-Hartmann sensor where the distortions were analyzed.

Figure 3 .
Figure 3. Shack-Hartmann focal spot field (hartmannogram) obtained after the beam propagation through a scattering layer.The inner yellow circle corresponded to the initial diameter (4 mm) of the collimated beam before entering the scattering volume.The outer red circle corresponded to the area (4.8 mm diameter) of the Shack-Hartmann sensor where the distortions were analyzed.

Figure 5
explains what we actually measure with a Shack-Hartmann sensor in this work.

Figure 5 .
Figure 5. Explanation of the term "averaged wavefront".Initially, before the scattering medium, the beam has a flat-incident wavefront.The scatterers with which radiation interacts become point sources of secondary spherical waves with their own wave vectors

Figure 6 .
Figure 6.Symmetrical Zernike coefficient dependence on the concentration value of scattering medium.

Figure 6 .
Figure 6.Symmetrical Zernike coefficient dependence on the concentration value of scattering medium.

Algorithms 2023 , 1 Figure 7 .
Figure 7. Photo of the experimental setup with a diode laser and a Shack-Hartmann sensor.

Figure 7 .
Figure 7. Photo of the experimental setup with a diode laser and a Shack-Hartmann sensor.

Figure 8 .
Figure 8. Trend lines of the dependence of the amplitude of distortions on the concentration values for the simulation (red curve) and for the experiment (purple curve).Image adapted from [37].

Figure 8 .
Figure 8. Trend lines of the dependence of the amplitude of distortions on the concentration values for the simulation (red curve) and for the experiment (purple curve).Image adapted from [37].

Figure 9 .
Figure 9.Initial intensity distribution of the simulated (a) and real, experimental (b) laser beam; schematical representation of laser beam propagation without diffraction (c) and with diffraction (d); reconstructed interference pattern with symmetrical distortions only (e) obtained in the simulation and with both symmetrical and non-symmetrical distortions (f) obtained in the real experiment.

Figure 9 .
Figure 9.Initial intensity distribution of the simulated (a) and real, experimental (b) laser beam; schematical representation of laser beam propagation without diffraction (c) and with diffraction (d); reconstructed interference pattern with symmetrical distortions only (e) obtained in the simulation and with both symmetrical and non-symmetrical distortions (f) obtained in the real experiment.

Table 1 .
Model verification: comparison of results obtained from the developed model and MCML.

Table 1
contains results of comparison of our model with MCML.