Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Graph-based model for adaptive simulation of beam propagation in turbulent media

Open Access Open Access

Abstract

A graph-based approach uses a triangular adaptive mesh for simulating the propagation of light beams through the atmosphere. In this approach, the atmospheric turbulence and the beam wavefront are signals in a graph, with vertices representing an irregular distribution of signal points and edges between vertices showing their relationships. The adaptive mesh provides a better representation of the spatial variations in the beam wavefront, resulting in increased accuracy and resolution compared to regular meshing schemes. The adaptability of this approach to the propagated beam characteristics makes it a versatile tool for simulating beam propagation in various turbulence conditions.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Simulation plays a crucial role in studying the propagation of light beams through random media like the atmosphere. It allows modeling and understanding of the complex interactions between light and atmospheric turbulence, which can be challenging to measure or study directly. Simulation can study the effects of turbulence on beam propagation, and explore the performance of distinctive light beams under different atmospheric conditions [14]. Simulation can also consider the optimization of optical systems in the atmosphere, such as laser communications, laser power beaming, or optical remote sensing, by identifying the optimal operating conditions and system design parameters [58].

1.1 Background

The split-step technique is a widely adopted approach for simulating the propagation of light beams through random media, including the atmosphere. This method is known for its high efficiency, precision, and ease of implementation in computer simulations. It works by dividing the total propagation distance into smaller, manageable steps and solving the wave equation for each step individually. The results for each step are then combined to arrive at the overall solution for the entire propagation distance [9,10].

The split-step technique uses the frequency domain to solve the parabolic-wave equation by representing the beam's phase and amplitude. A regular spatial mesh is required to divide the simulation space into equally sized partitions in all directions, which enables the efficient use of discrete Fourier transforms. However, implementing the spectral method for atmospheric propagation has limitations due to turbulence's significant impact on fields. A large spatial mesh with small partitions is necessary to ensure accurate numerical simulations.

In numerical solutions of the parabolic-wave equation on a mesh, the error estimate for the optical field solution is determined by [11]

$${\epsilon \le \frac{1}{M}\mathop \sum \limits_m h_m^2{\; }|{{\nabla^2}U} |,}$$
were m represents the partitions of the numerical mesh, ${h_m}$ is the corresponding partition size, and U is the optical field. The Laplacian ${\nabla ^2}U$ measures field curvature by comparing the field value to its average at nearby points. Eq. (1) shows that reducing the mesh partition size ${h_m}$ may not be necessary to minimize error. Instead, a strategy balancing partition-wise error contributions $h_m^2\; |{{\nabla^2}U} |$ is critical. ${h_m}$ should be proportional to ${|{{\nabla^2}U} |^{ - 1/2}}$ to achieve this balance, resulting in smaller partition sizes for rough (high-curvature) solutions U and larger partitions for smooth ones.

1.2 Overview of the proposed method

Irregular or adaptive meshes featuring varying partition sizes and shapes offer an effective solution to reduce errors by balancing curvature with partition size (see Fig. 1). The Delaunay triangulation method is a reliable technique for creating irregular meshes by forming triangles between points on a surface such that no point lies inside the circumcircle of any triangle [12]. Irregular triangular meshes, commonly used in computational geometry and computer graphics, are particularly useful in simulations dealing with significant wavefront curvatures or discontinuities. They can better capture wavefront behavior by using fewer partitions in specific regions while providing higher resolution where necessary, making them ideal for simulations with multiple length scales.

 figure: Fig. 1.

Fig. 1. Illustrates the transformation of the beam wavefront into a signal in a graph, where vertex ${{\boldsymbol r}_n}$ are used to represent an uneven distribution of signal points in space, and edges ${e_m}$ are used to represent the relationships between these signal points. The discrete representation of the Laplacian operator on an irregular triangular mesh takes into account the cotangents of angles ${\alpha _m}$ and ${\beta _m}$ opposite each edge ${e_m}$.

Download Full Size | PDF

Irregular triangular meshes present challenges for conventional signal processing and numerical modeling techniques due to their complex structures. Graph theory provides a solution for modeling and analyzing such meshes. It involves representing mesh points as graph vertices, interconnecting them with edges, and using an edge weight function to define the significance of the relationships between points [13,14]. This study aims to leverage graph theory and develop a graph-based approach to simulate beam propagation through turbulence using an irregular adaptive mesh and the split-step method.

The paper's organization is as follows: Section 2 explains the graph-based propagation approach. Section 3 then discusses the beam propagation method using triangular mesh adaptivity. Finally, Section 4 presents the conclusions of the study.

2. Graph-based split-step method

2.1 Split-step approach and its limitations

The split-step approach solves the parabolic-wave equation [15]

$${\frac{{\partial U({z,{\boldsymbol r}} )}}{{\partial \textrm{z}}} = \frac{j}{{2k}}[{{\nabla^2} + 2{k^2}\; n({z,{\boldsymbol r}} )} ]\; U({z,{\boldsymbol r}} ),}$$
by restricting the propagation to the forward direction. Equation (2) describes the scalar optical field $U({z,{\boldsymbol r}} )$, where z is a Cartesian coordinate in the direction of propagation and ${\boldsymbol r}$ represents the transverse coordinates. The transverse coordinates are not necessarily Cartesian. The free-space wavenumber is $k = 2\pi /\lambda $ and $n({z,{\boldsymbol r}} )$ describes the random fluctuations in the local refraction index resulting in optical turbulence. By expressing the fields $U({z,{\boldsymbol r}} )$ in terms of the eigenfunctions ${\phi _n}({\boldsymbol r} )$ of the transverse Laplacian operator ${\nabla ^2}$
$${U({z,{\boldsymbol r}} )= \mathop \sum \limits_n {u_n}(z ){\phi _n}({\boldsymbol r} ),}$$
with ${u_n}(z )= \smallint U({z,{\boldsymbol r}} )\phi _n^\ast ({\boldsymbol r} )\; d{\boldsymbol r}$, and where ${\nabla ^2}{\phi _n}({\boldsymbol r} )= {\lambda _n}^2{\phi _n}({\boldsymbol r} )$ and $\smallint {\phi _n}^\ast ({\boldsymbol r} ){\phi _m}({\boldsymbol r} )d{\boldsymbol r} = {\delta _{nm}}$, the split-step solution of the parabolic-wave Eq. (2) at $z + \Delta z$ is found [9]:
$${U({z + \Delta z,{\boldsymbol r}} )= exp[{\; j{\; }S({\Delta z,{\boldsymbol r}} )\; } ]\mathop \sum \limits_n {u_n}(z )\; exp\left( {\frac{j}{{2k}}{\lambda_n}^2\Delta z} \right)\; {\phi _n}({\boldsymbol r} ).}$$

This solution (4) shows that the exponential terms containing the eigenvalues ${\lambda _n}^2$ of the Laplacian ${\nabla ^2}$ propagate the optical field a step $\Delta z$ forward in a uniform media, while the exponential containing the turbulence phase fluctuations $S({\Delta z,{\boldsymbol r}} )= k\mathop \smallint \nolimits_{\Delta z} n({z,{\boldsymbol r}} )dz$ account for the inhomogeneous refractive index [16,17]. The significant result of this approach is that the extended random medium can be divided into slabs of thickness $\Delta z$, each introducing an independent random contribution $S({\Delta z,{\boldsymbol r}} )$ to the wave phase but no significant change in amplitude. Amplitude variations build up by diffraction over several $\Delta z$. For simplicity, Eq. (4) ignores a first step of uniform media propagation from $z = 0$ to $z = \Delta z/2$ necessary to initiate the splitting solution.

The turbulence phase fluctuations $S({\Delta z,{\boldsymbol r}} )$ differ from the propagation exponential terms as it does not depend on the Laplacian eigenfunction ${\phi _n}({\boldsymbol r} )$. Instead, phase fluctuations are a homogeneous, normal, and zero-averaged random field characterized by a covariance function ${\Gamma _s}({{{\boldsymbol r}_1},{{\boldsymbol r}_2}} )$. It measures the similarity between the phase fluctuations at two different spatial points and can be used to describe the turbulence affecting the wave. The turbulence phase can also admit a spectral representation

$${S({\Delta z,{\boldsymbol r}} )= \mathop \sum \limits_n {s_n}({\Delta z} ){\psi _n}({\boldsymbol r} )}$$
in terms of the eigenfunctions ${\psi _n}$ of the covariance integral equation $\smallint {\Gamma _s}({{{\boldsymbol r}_1},{{\boldsymbol r}_2}} ){\psi _n}({{{\boldsymbol r}_2}} )\; d{{\boldsymbol r}_2} = {\sigma _n}^2{\psi _n}({{{\boldsymbol r}_1}} )$, where the constants ${\sigma _n}^2$ are the corresponding eigenvalues [18]. Typically, zero-averaged circular statistics are used for the complex spectral amplitudes $\langle{s_n}\rangle$, with $\langle{s_n}\rangle = 0$, $\langle{s_n}{s_m}\rangle = 0$, and $\langle{s_n}s_m^\ast\rangle{=} {\sigma _n}^2{\delta _{nm}}$. The eigenfunctions ${\psi _n}$ give the Karhunen–Loeve basis functions of the covariance function.

The split-step approach (4) solves the parabolic-wave equation in the spectral domain. It takes advantage of the fact that optical field U and turbulence S admit representation in terms of the eigenfunctions ${\phi _n}$ and ${\psi _n}$ of the transverse Laplacian operator and the covariance function, respectively. The derivation of the splitting solution (4) is independent of the eigenfunctions, and any coordinate system is possible as long as a numerical representation of the associated eigenfunction basis is available.

For example, in Cartesian coordinates ${\boldsymbol r} = ({x,y} )$, the Laplacian eigenfunctions become plane waves ${\phi _n} = exp({ - j\; {{\boldsymbol K}_n}\Delta {\boldsymbol r}} )$ in the spatial direction ${{\boldsymbol K}_n} = ({{K_{nx}},{K_{ny}}} )$ with eigenvalues ${\lambda _n}^2 = \; {K_n}^2$ [19]. If the covariance function only depends on the Euclidean norm $\Delta r = |{{{\boldsymbol r}_1} - {{\boldsymbol r}_2}} |$, plane waves are also the eigenfunctions ${\psi _n}({\boldsymbol r} )$ for the covariance function. The corresponding eigenvalue $\smallint {\Gamma _s}({\Delta r} )\; exp({ - j\; {K_n}\Delta r} )\; d\Delta r = {\sigma _n}^2({{K_n}} )$ is the power spectral function calculated as the Fourier transform of ${\Gamma _s}({\Delta r} )$ at frequency ${K_n}$.

The spectral method in the Cartesian coordinates solves the wave equation by representing the solutions in terms of a spatial Fourier series. The discrete Fourier transform algorithm is most efficient when applied to evenly-spaced signals, which is why it typically uses a regular spatial mesh. The mesh must be of a finite extent, with N equally spaced points ${{\boldsymbol r}_n}$ and frequencies ${{\boldsymbol K}_n}$.

2.2 Graph model for irregular meshes

As mentioned earlier, achieving precise numerical simulations of atmospheric beam propagation often requires a spatial grid with a high degree of regularity and fine meshing, which includes small and reduced elements. However, irregular meshes can potentially provide even more precise results. That is because irregular meshes, such as Delaunay triangulations, can capture more accurately the wavefront complexity and its variations in different regions.

Graphs can model irregular meshes by using N points as vertices V and M weighted edges E connecting them (refer to Fig. 1). Each vertex ${V_n}$, with n ranging from $1$ to N, is associated with a reference mesh-space coordinate ${{\boldsymbol r}_{\boldsymbol n}}$ which spans the simulation mesh, and additional quantities, such as the optical field $U = U({z,{{\boldsymbol r}_{\boldsymbol n}}} )$ and turbulence $S = S({\Delta z,{{\boldsymbol r}_{\boldsymbol n}}} )$. Each weight ${W_m}$ on the edge ${e_m}$, with m ranging from $1$ to M, can be calculated in different ways depending on the specific application and type of analysis being performed. As an example, when examining the propagation of light beams through the atmosphere, the weighting of graphs can be determined by either the edge lengths or the triangle angles within the non-uniform mesh.

In graph theory, the Laplacian matrix L is a square matrix used to represent the weighted connectivity of a graph [13]. It is the difference between the degree matrix and the adjacency matrix. The degree matrix is a diagonal matrix whose entries represent the sum of the weights of all the edges incident to a particular vertex. On the other hand, the adjacency matrix is a matrix whose entries indicate the presence or absence of edges between the vertices, and the value of each entry represents the edge weight. The Laplacian matrix incorporates the properties of a graph and appears in many graph analysis techniques. One application of the Laplacian matrix is the calculation of eigenvalues and eigenvectors, which can analyze the graph structure in the spectral domain [13]. Since the Laplacian matrix of an undirected graph is positive semidefinite, all eigenvalues are real and nonnegative and form a complete set of orthogonal eigenvectors.

The eigenvectors of the Laplacian matrix correspond to different spectral modes of the graph. The eigenvectors with smaller eigenvalues capture the global structure of the graph, while the eigenvectors with larger eigenvalues capture the local structure of the graph. The form of the eigenvectors can be complex and depends on the topology and connectivity of the graph [13]. For instances, the first eigenvector of the Laplacian matrix is the constant vector and represents the center of mass of the graph. Other eigenvectors of the Laplacian matrix describe the graph's connectivity properties. For example, the second eigenvector of the Laplacian matrix, or Fiedler vector, is often used to partition the graph. Partitioning the graph allows for the graph separation into two parts. These parts are characterized by being well-connected within themselves, but with relatively few connections between them.

In this study, a multigraph representation $G = ({V,{E^C},{E^\Gamma },{W^C},{W^\Gamma }} )$ encodes the attributes of both the propagation of optical fields U and random phase turbulence S defined over the same set of vertices. When investigating wave propagation over the graph, the cotangent Laplacian matrix ${L^C}$ [20] can define the Laplacian operator ${\nabla ^2}$ used by the wave equation. The edges ${E^C}$ in this scenario connect every vertex V to its nearest neighbors, with weights ${W^C}$ calculated adding the cotangents of the angles opposite to each edge, which allows for the discrete representation of the Laplacian operator on the irregular triangular mesh (see Fig. 1). The eigenvalues and eigenvectors of the cotangent Laplacian matrix ${L^C} = \varPhi \Lambda {\varPhi ^T}$, where $\varPhi = [{{\phi_0} \ldots \; {\phi_{N - 1}}} ]$ is the matrix of the N eigenvectors of ${L^C}$, and $\Lambda = \textrm{diag}[{{\lambda_0}^2 \ldots \; {\lambda_{N - 1}}^2} ]$ is the matrix of corresponding eigenvalues, are the base for the implementation of the splitting solution (4) of the parabolic-wave equation in the graph spectral domain and the spectral representation (3) of the fields on the graph ${u_n}$. By rearranging the fields defined on the graph vertices into a one-dimensional vector U of length N, it results

$${{u_n} = \varPhi U.}$$

In addition, when modeling random turbulence $S({\Delta z,{{\boldsymbol r}_{\boldsymbol n}}} )$, the covariance Laplacian matrix ${L^\varGamma }$ [21] characterizes the similarity of phase fluctuations at different points within the graph. The edges ${E^\varGamma }$ connect every pair of vertices in the graph and the weights ${W^\varGamma }$ represent the covariance function ${\varGamma _s}({{{\boldsymbol r}_1},{{\boldsymbol r}_2}} )$ between those pairs. This allows the Laplacian matrix ${L^\varGamma }$ to depict the relationships between the phase fluctuations in the graph. The eigenvalues and eigenvectors of the covariance Laplacian matrix ${L^\varGamma } = \varPsi \varSigma {\varPsi ^T}$, with $\varPsi = [{{\psi_0} \ldots \; {\psi_{N - 1}}} ]$ and $\varSigma = \textrm{diag}[{{\sigma_0}^2 \ldots \; {\sigma_{N - 1}}^2} ]$, enable the generation of random phase turbulence S using the graph spectral representation (5). By defining a random vector ${s_n}$, where each element is independent and has its variance determined by the corresponding element on the diagonal of $\varSigma $, leads to

$${S = {\varPsi ^T}{s_n}.}$$

The graph-based split-step approach improves propagation simulations by utilizing eigenvectors of the cotangent and covariance Laplacian matrices. It provides a precise representation of the optical fields and defines the turbulence phase with exact statistical properties. Despite requiring frequent transformations between space and spectral domains using the matrix of eigenvectors and its inverse, as well as computing the exponential with turbulence phase fluctuations, the computational complexity of the approach is comparable to that of the spectral split-step method. The process of factorizing the Laplacian matrices is simple and accurate for small graphs, which are required for beam propagation simulations with mesh adaptivity. Furthermore, as shown below, the graph-based method employs a remeshing strategy that balances accuracy and efficiency by reducing mesh complexity while preserving its overall shape. The process of reducing the number of vertices is relatively simple and accurate for small graphs required for beam propagation simulations with mesh adaptivity.

Additionally, the method leverages the graph permutation equivalence, meaning that the order of the vertices in the graph does not impact the simulation outcomes. It abstracts away the absolute position information of the mesh nodes and encodes optical fields and turbulence into vertex features. The edges between vertices represent relationships that are invariant to their spatial positions.

Both the discrete Fourier split-step and the proposed graph-based split-step techniques aim to numerically solve the one-way parabolic-wave Eq. (2) on an unbounded domain. While the discrete Fourier transform method offers the advantage of using a computationally efficient FFT algorithm, its periodicity constraints limit its applicability. In contrast, the proposed graph-based split-step approach employs the graph Laplacian and its eigenvectors to define the simulation spectral domain, which allows for a wider range of applications without the periodicity restrictions imposed by the Fourier method.

3. Beam propagation with triangular mesh adaptivity

3.1 Triangular mesh generation and adaptivity

The split-step solution on an irregular triangular mesh is an effective tool for investigating wave propagation through random turbulence. An adaptive meshing approach adjusts the size and shape of the mesh based on the local field solution to achieve optimal results.

The remeshing strategy used in this study focuses on wavefront curvature $|{{\nabla^2}U} |$, which represents the deviation of optical wavefronts from a flat surface [22]. It quantifies the local radius of curvature of the wavefront and reflects the local focusing or defocusing effect. The optical wavefront is complex and defined by both its amplitude and phase. The strategy involves refining the mesh in areas of high wavefront curvature and simplifying it in areas of low wavefront curvature, ensuring the mesh is optimally adapted to the irregularities of the optical fields throughout the simulation, regardless of the coarseness or regularity of the initial mesh. This strategy is a development of previous work on remeshing triangle mesh surfaces [23].

At each simulation stage, the triangular mesh represents the wavefront's shape, consisting of vertices and edges that sample its complex field amplitude. The remeshing process starts with an initial mesh and dynamically adjusts its composition to capture changes in the wavefront's curvature as it propagates. The vertices and edges are added, removed, or manipulated through specific operations such as splitting, collapsing, and flipping, as shown in Fig. 2(a), to reflect the evolving shape of the wavefront and ensure accurate sampling of its complex amplitude. This results in a remeshed mesh representing the changing complex wavefront.

 figure: Fig. 2.

Fig. 2. (a) The remeshing process involves manipulating edges and vertices through various operations, such as splitting, collapsing, and flipping, to create an optimal mesh distribution. (b) The wavefront curvature is calculated by evaluating the angle $\theta $ between the normal vectors at the start and end vertices of each edge in the mesh.

Download Full Size | PDF

The concept of curvature is defined in terms of the directional derivative of the normal vector of a wavefront and is calculated from the differences in the normal vectors of vertices [24]. The evaluation of an edge during the remeshing process starts with the assignment of a normal vector $\hat{n}$ to each vertex, as shown in Fig. 2(b). This is achieved by summing the normal vectors of the faces associated with that vertex. Next, the angle $\theta $ between the normal vectors ${\hat{n}_1}$ and ${\hat{n}_2}$ of the start and end vertices of an edge is calculated. If the calculated angle $\theta $ exceeds a pre-determined maximum threshold ${\theta _{max}}$, it is considered as high wavefront curvature, and the edge is split. On the other hand, if the calculated angle $\theta $ falls below a minimum threshold ${\theta _{min}}$, it is deemed as low wavefront curvature and the edge is collapsed. The resulting mesh triangles are evaluated for proper shape through Delaunay checking, and edge flipping may be necessary to ensure that the triangles are well-formed. When splitting a mesh, new vertices are positioned through interpolation of their spatial coordinates and the optical wavefront at the midpoint of the edge. The weights of the new edges are also calculated in this process.

The remeshing strategy also aims to strike a balance between accuracy and efficiency by reducing the complexity of the mesh. It considers lowering the number of vertices while preserving the overall shape of the mesh. The process operates by identifying and removing vertices that have minimal impact on the mesh shape. A wavefront power-density criteria determine each vertex impact. It calculates the power density at each vertex and only removes vertices with a power density below a predefined factor K of the mean wavefront power density.

3.2 Error analysis and convergence

Figure 3 typifies the results of these studies. It reveals that the adaptive mesh approach consistently delivers better optical field quality than non-adaptive methods, even for comparable mesh sizes. It considers transmission of an initial Gaussian beam ${U_0}({\boldsymbol r} )= \textrm{exp}({ - {\; }{{|{\boldsymbol r} |}^2}/2\omega_0^2} )$, characterized by a reference beam waist ${\mathrm{\omega }_0}$, that is propagating through a turbulent media described by the turbulence phase correlation function [16]

$${{\varGamma _s}({{{\boldsymbol r}_1},{{\boldsymbol r}_2}} )= {{({|{{{\boldsymbol r}_1} - {{\boldsymbol r}_2}} |/{r_c}} )}^\alpha }.}$$

Here, the parameter $\alpha $, usually within the range of $1 < \alpha < 2$, is set to $5/3$. The wave coherence radius ${r_c}$ fully characterizes the statistics of phase fluctuations for Kolmogorov turbulence along a single propagation step $\Delta z$. The simulation considers $32$ propagation steps, with a step size of $\Delta z = 125$ m. At a wavelength $\lambda = 1550$ nm and a beam with ${\omega _0} = 5$ cm, ${r_c}$ is approximately $0.15$ m, under the moderate-to-strong turbulence conditions of the numerical example analyzed in Fig. 3. However, after the final propagation step, the corresponding wave coherence radius ${r_0}$ is only about one-tenth of ${r_c}$.

 figure: Fig. 3.

Fig. 3. (a) Statistical relationships between the measured wavefront curvature and the distribution of edge lengths h in the simulation mesh. The inset provides a detailed view of the dependencies for small edge lengths in a semi-log scale. (b) Simulation error estimate as a function of the propagation distance z, normalized to the total propagation range L. The plot compares the performance of three different threshold angles ${\theta _{max}}$. (c) Graph-based spatial spectra of the intensity I of the adaptive mesh, and the intensity error i resulting from the lack of adaptation.

Download Full Size | PDF

The initial mesh distributes $N = 1024$ points evenly across the area of a circle using a canonical Fibonacci spiral method, achieving a similar point density to a regular $32 \times 32$ square mesh. The use of a Fibonacci mesh results in regular and isotropic resolution, with an equal area for each mesh point [25]. To create the initial mesh, the points are connected with $M = 2012$ edges through Delaunay triangulation. For this example, the adaptive remeshing scheme employs curvature thresholds ${\theta _{max}} = 30^\circ $ and ${\theta _{min}} = 1^\circ $, and a power density threshold of $K ={-} 30$ dB. It adjusts the spatial resolution at every other propagation step. The study generates ${10^3}$ instances of graph-based propagated beams and collects statistics from the simulations to evaluate the performance of mesh adaptivity. The distribution of points in the simulation mesh ranges from regular to highly irregular. The number of points in the mesh dynamically adjusts from the initial count of $1024$ to an average of $2162$ during the simulation. Similarly, the number of edges in the mesh changes from the initial $2012$ to an average of $6432$.

Figure 3(a) displays the measured statistical relationships between wavefront curvature and the distribution of edge lengths in the mesh. The edge length h, which refers to the Euclidean distance between two linked vertices, is normalized by the wave coherence radius ${r_0}$. To improve statistical reliability, the ${10^3}$ instances of propagation were averaged. The plots show that after passing through turbulence, the measured wavefront curvature ranges from $1$ to ${10^5}$, leading to an irregular mesh structure that predicts random optical fields at various scales. The edge length distribution is accurately modeled by a beta distribution with optimal shape parameters of $\alpha = 2.4$ and $\beta = 17.5$. The mesh maintains an optimal level of detail throughout the simulation area, with high resolution in regions of high wavefront curvature and low resolution in regions of low wavefront curvature. The average edge length is slightly smaller than one eighth of the coherence radius ${r_0}$, with the longest edges reaching approximately two thirds of the coherence radius, and the shortest edges down to one hundredth of the coherence radius. This adaptive meshing approach allows for simulations at a level of resolution that is not possible with regular meshing schemes.

Figure 3(b) illustrates the simulation error estimate as a function of the propagation distance z. The simulation error, defined in Eq. (1), is computed by summing the contribution terms $h_m^2\; |{{\nabla^2}U} |$ for each edge $m = 1, \ldots ,M$, where ${h_m}$ is the edge length and ${\nabla ^2}U$ is the measured wavefront curvature. The analysis supports the intuitive expectation that in a regular mesh, with regular edge length ${h_0}$, the error $\epsilon \le 1/Mh_0^2\mathop \sum \nolimits_m |{{\nabla^2}U} |$ increases linearly with distance, in proportion to the increase in wavefront curvature.

On the other hand, the adaptive mesh adjusts the edge length optimally based on wavefront curvature, using ${h_m} = C{|{{\nabla^2}U} |^{ - 1/2}}$, which limits the error to $\epsilon \le {C^2}$. The constant ${C^2}$ observed in the simulations, which represents the upper limit of the error, mostly depends on the pre-defined maximum angle ${\theta _{max}}$ between the normal vectors at the starting and end vertices of the evaluated edge. For the standard ${\theta _{max}} = 30^\circ $ considered in the simulations, the error levels off early in the propagation. By decreasing ${\theta _{max}}$ to $20^\circ $, the measured error can be reduced, although it increases the mesh complexity. This finding demonstrates that mesh adaptivity can effectively improve the simulation accuracy.

Figure 3(c) displays the intensity error i as another metric of simulation accuracy. The intensity ${|{U({\boldsymbol r} )} |^2}$ in the adaptive mesh is written as $I = {I_0} + i$, where ${I_0}$ is the intensity in the regular mesh, and i is the error in the intensity error resulting from the lack of adaptation. Here, mesh adaptivity defines the high-resolution ground truth of this error analysis, as confirmed by the findings of Fig. 3(b). As the error is more easily examined in the transform domain, the analysis estimates the graph-based spectra of intensity ${I_n} = \varPhi I$ and intensity error ${i_n} = \varPhi i$. As shown in the figure, the error is nearly white, meaning that the errors are evenly distributed across all spatial scales. Initially, the regular mesh has errors introduced at high graph frequencies, but after propagation, they affect uniformly the entire spectrum. As expected, the intensity error is particularly prominent at high spatial frequencies. The intensity ${I_n}$ and intensity error ${i_n}$ spectra were averaged over ${10^3}$ instances to improve the accuracy of the results.

3.3 Modeling and simulation results

Figure 4 presents the results of the graph-based simulations for various types of optical beams. The parameters used in these simulations are consistent with those described in Fig. 3. The first case (a) illustrates the propagation of an off-axis Gaussian beam with an asymmetrical (ellipsoidal) shape [26]. Off-axis Gaussian beams are the starting point for creating more complex beams, such as hyperbolic and sinusoidal-Gaussian beams. This example considers a Gaussian beam that is axially displaced by ${\mathrm{\omega }_0}/4$ and has an elliptical intensity distribution rotated by an angle of $\mathrm{\pi }/6$ radians with an eccentricity of 0.6. The second case (b) focuses on the propagation of a vortex beam, having a characteristic helical phase and carrying orbital angular momentum or OAM [27]. Vortex beams have a broad range of applications, including optical free-space communication, optical remote sensing, and transmission of quantum information states. In this example, the vortex beam carries OAM of $1\hbar $ per photon along the beam axis. The third case (c) in the study investigates the propagation of a partially coherent beam composed of $45$ mutually orthogonal Laguerre-Gauss modes, all propagating simultaneously along a common central axis. Spatial-mode multiplexing, as used in near-field free-space optical communication, can significantly enhance transmission capacity by transmitting independent information streams through orthogonal modes. In such scenarios, Laguerre-Gauss modes with a chosen beam waist, based on the field coherence length in the receiver plane, are optimal for transmission with minimum degradation through turbulence [28]. In this example, the simulation parameters support the propagation of these $45$ Laguerre-Gauss modes over a $4$-km range.

 figure: Fig. 4.

Fig. 4. Examples of graph-based atmospheric beam propagation using triangular mesh adaptivity. The figure displays both the amplitude (left) and phase (right) wavefronts for different types of beams, including (a) an off-axis Gaussian beam with an asymmetrical (ellipsoidal) shape, (b) a vortex beam with characteristic helical phase, and (c) a partially coherent beam composed of 45 mutually orthogonal Laguerre-Gauss modes.

Download Full Size | PDF

Figure 4 presents both the amplitude and phase wavefronts, showcasing the effectiveness of the graph-based approach in adapting to the specific features of the propagated beams. The asymmetrical off-axis Gaussian beam shows that with the increase of propagation distance, beam spreading, power scintillation, and phase aberrations exhibit a transverse coordinate dependency, which the adaptive mesh effectively captures. On the other hand, the vortex beam has an unstable phase singularity and changes in the OAM state, which the adaptive mesh can explain. Lastly, the partially coherent beam in the atmosphere is subject to phase distortions and amplitude fluctuations that lead to a degradation of coherence and limit the number of independent modes. The adaptive mesh approach proves valuable in this scenario by providing a more precise depiction of modal coupling.

4. Conclusion

The results of this study provide several important insights into the effectiveness of the graph-based approach for simulating light propagation through atmospheric turbulence. First, the approach models accurately turbulence and the beam wavefront as signals in a graph. By using a split-step graph spectral method, the technique accurately simulates beam propagation through the turbulent media. Second, using an adaptive mesh is a critical aspect of this approach. The adaptive mesh distributes nodes unevenly and adjusts the mesh based on wavefront curvature during propagation, leading to a more accurate representation of the optical fields. The results show that even a relatively small number of nodes in the mesh could provide high-resolution detail in areas of high wavefront curvature and low resolution in places of low wavefront curvature. Third, the error analysis shows that the graph-based beam propagation with triangular mesh adaptivity provides improved beam propagation quality compared to non- adaptive methods. It highlights the effectiveness of this approach in simulating beam propagation through atmospheric turbulence.

Finally, the proposed graph-based approach can adapt to the specific features of the propagated beams. The technique works with wavefronts of arbitrary amplitude or phase and under diverse turbulence conditions, providing a more precise depiction of the random optical fields.

Funding

Agencia Estatal de Investigación (PID2020-118410RB-C21).

Disclosures

The author declares no conflicts of interest.

Data availability

No data were generated or analyzed in the presented research.

References

1. M. A. Cox, N. Mphuthi, I. Nape, N. Mashaba, L. Cheng, and A. Forbes, “Structured Light in Turbulence,” IEEE J. Sel. Top. Quantum Electron. 27(2), 1–21 (2021). [CrossRef]  

2. Y. Zhu, Y. Zhang, and Z. Hu, “Spiral spectrum of Airy beams propagation through moderate-to-strong turbulence of maritime atmosphere,” Opt. Express 24(10), 10847–10857 (2016). [CrossRef]  

3. M. Charnotskii, “Wave propagation modeling with non-Markov phase screens,” J. Opt. Soc. Am. A 33(4), 561–569 (2016). [CrossRef]  

4. R. Frehlich, “Simulation of laser propagation in a turbulent atmosphere,” Appl. Opt. 39(3), 393–397 (2000). [CrossRef]  

5. E. Villaseñor, M. He, Z. Wang, R. Malaney, and M. Z. Win, “Enhanced uplink quantum communication with satellites via downlink channels,” IEEE Trans. Quantum Eng. 2, 1–18 (2021). [CrossRef]  

6. C. Robert, J. Conan, and P. Wolf, “Impact of turbulence on high-precision ground-satellite frequency transfer with two-way coherent optical links,” Phys. Rev. A 93(3), 033860 (2016). [CrossRef]  

7. A. Belmonte, “Feasibility study for the simulation of beam propagation: consideration of coherent lidar performance,” Appl. Opt. 39(30), 5426–5445 (2000). [CrossRef]  

8. J. A. Fleck Jr., J. R. Morris, and M. D. Feit, “Time-dependent propagation of high energy laser beams through the atmosphere,” Appl. Phys. 10(2), 129–160 (1976). [CrossRef]  

9. M. Z. M. Jenu and D. H. O. Bebbington, “Simulation of optical propagation in a turbulent atmosphere using the split step method,” in 8th International Conference on Antennas and Propagation (1993), pp. 817–822.

10. M. Spivack and B. J. Uscinski, “The split-step solution in random wave propagation,” J. Compt. Appl. Math. 27(3), 349–361 (1989). [CrossRef]  

11. K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Computational Differential Equations (Cambridge University Press, 2019).

12. M. Berg, O. Cheong, M. Kreveld, and M. Overmars, Computational Geometry: Algorithms and Applications (Springer, 2008).

13. J. L. Gross, J. Yellen, and M. Anderson, Graph Theory and Its Applications (Chapman and Hall, 2018).

14. A. Sandryhaila and J. M. F. Moura, “Discrete Signal Processing on Graphs,” IEEE Trans. Signal Process. 61(7), 1644–1656 (2013). [CrossRef]  

15. S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskii, “Principles of Statistical Radiophysics,” in Wave Propagation through Random Media (Springer, 1989).

16. J. Goodman, Statistical Optics (John Wiley & Sons, 2015).

17. R. Lane, A. Glindemann, and J. Dainty, “Simulation of a Kolmogorov phase screen,” Waves Random Media 2(3), 209–224 (1992). [CrossRef]  

18. H. Stark and J. W. Woods, Probability, Random Processes, and Estimation Theory for Engineers (Prentice Hall, 2011).

19. R. N. Bracewell, The Fourier Transform and Its Applications (McGraw-Hill, 1999).

20. A. I. Bobenko and B. A. Springborn, “A discrete Laplace-Beltrami operator for simplicial surfaces,” Discrete Comput. Geom. 38(4), 740–756 (2007). [CrossRef]  

21. M. Belkin and P. Niyogi, “Laplacian Eigenmaps for Dimensionality Reduction and Data Representation,” Neural Computation 15(6), 1373–1396 (2003). [CrossRef]  

22. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 2019).

23. P. Alliez, M. Attene, C. Gotsman, and G. Ucelli, “Recent advances in remeshing of surfaces,” in Shape Analysis and Structuring (Springer, 2007).

24. S. Rusinkiewicz, “Estimating curvatures and their derivatives on triangle meshes,” in 2nd International Symposium on 3D Data Processing, Visualization and Transmission (2004), pp. 486–493.

25. R. Swinbank and R. J. Purser, “Fibonacci grids: A novel approach to global modelling,” Q. J. R. Meteorol. Soc. 132(619), 1769–1793 (2006). [CrossRef]  

26. Y. Baykal and H. T. Eyyuboğlu, “Off-axis-Gaussian beam intensity in random medium,” Opt. Commun. 269(1), 30–38 (2007). [CrossRef]  

27. B. Rodenburg, M. P. J. Lavery, M. Malik, M. N. O’Sullivan, M. Mirhosseini, D. J. Robertson, M. Padgett, and R. W. Boyd, “Influence of atmospheric turbulence on states of light carrying orbital angular momentum,” Opt. Lett. 37(17), 3735–3737 (2012). [CrossRef]  

28. A. Belmonte and J. M. Kahn, “Optimal modes for spatially multiplexed free-space communication in atmospheric turbulence,” Opt. Express 29(26), 43556–43566 (2021). [CrossRef]  

Data availability

No data were generated or analyzed in the presented research.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

Fig. 1.
Fig. 1. Illustrates the transformation of the beam wavefront into a signal in a graph, where vertex ${{\boldsymbol r}_n}$ are used to represent an uneven distribution of signal points in space, and edges ${e_m}$ are used to represent the relationships between these signal points. The discrete representation of the Laplacian operator on an irregular triangular mesh takes into account the cotangents of angles ${\alpha _m}$ and ${\beta _m}$ opposite each edge ${e_m}$.
Fig. 2.
Fig. 2. (a) The remeshing process involves manipulating edges and vertices through various operations, such as splitting, collapsing, and flipping, to create an optimal mesh distribution. (b) The wavefront curvature is calculated by evaluating the angle $\theta $ between the normal vectors at the start and end vertices of each edge in the mesh.
Fig. 3.
Fig. 3. (a) Statistical relationships between the measured wavefront curvature and the distribution of edge lengths h in the simulation mesh. The inset provides a detailed view of the dependencies for small edge lengths in a semi-log scale. (b) Simulation error estimate as a function of the propagation distance z, normalized to the total propagation range L. The plot compares the performance of three different threshold angles ${\theta _{max}}$. (c) Graph-based spatial spectra of the intensity I of the adaptive mesh, and the intensity error i resulting from the lack of adaptation.
Fig. 4.
Fig. 4. Examples of graph-based atmospheric beam propagation using triangular mesh adaptivity. The figure displays both the amplitude (left) and phase (right) wavefronts for different types of beams, including (a) an off-axis Gaussian beam with an asymmetrical (ellipsoidal) shape, (b) a vortex beam with characteristic helical phase, and (c) a partially coherent beam composed of 45 mutually orthogonal Laguerre-Gauss modes.

Equations (8)

Equations on this page are rendered with MathJax. Learn more.

ϵ 1 M m h m 2 | 2 U | ,
U ( z , r ) z = j 2 k [ 2 + 2 k 2 n ( z , r ) ] U ( z , r ) ,
U ( z , r ) = n u n ( z ) ϕ n ( r ) ,
U ( z + Δ z , r ) = e x p [ j S ( Δ z , r ) ] n u n ( z ) e x p ( j 2 k λ n 2 Δ z ) ϕ n ( r ) .
S ( Δ z , r ) = n s n ( Δ z ) ψ n ( r )
u n = Φ U .
S = Ψ T s n .
Γ s ( r 1 , r 2 ) = ( | r 1 r 2 | / r c ) α .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.