Wave-optics investigation of branch-point density

Abstract. We use wave-optics simulations to investigate branch-point density (i.e., the number of branch points within the pupil-phase function) in terms of the grid sampling. The goal for these wave-optics simulations is to model plane-wave propagation through homogeneous turbulence, both with and without the effects of a finite inner scale modeled using a Hill spectrum. In practice, the grid sampling provides a gauge for the amount of branch-point resolution within the wave-optics simulations, whereas the Rytov number, Fried coherence diameter, and isoplanatic angle provide parameters to setup and explore the associated deep-turbulence conditions. Via Monte Carlo averaging, the results show that without the effects of a finite inner scale, the branch-point density grows without bound with adequate grid sampling. However, the results also show that as the inner-scale size increases, this unbounded growth (1) significantly decreases as the Rytov number, Fried coherence diameter, and isoplanatic angle increase in strength and (2) saturates with adequate grid sampling. These findings imply that future developments need to include the effects of a finite inner scale to accurately model the multifaceted nature of the branch-point problem in adaptive optics.

dislocations, 9,10 screw dislocations, 11,12 optical vortices, 13,14 etc.), which typically require the use of high-fidelity and wave-optics simulations. [15][16][17][18][19][20][21][22] For example, Voitsekhovich et al. 23 were the first to perform a wave-optics investigation of branch-point density (i.e., the number of branch points within the pupil-phase function). They did so as a function of propagation distance with the effect of a finite inner scale but for a fixed grid sampling. As previously mentioned, Fried was the first to describe the branch-point problem in terms of the hidden phase. 1 Since then, researchers have proposed several branch-point-tolerant phase reconstruction algorithms. [24][25][26][27][28][29] These algorithms, at large, have had limited degrees of success due to the multifaceted nature of the branch-point problem in adaptive optics. Most recently, members of the Starfire Optical Range at the Air Force Research Laboratory investigated the aggregate behavior of branch points, [30][31][32][33][34][35][36][37] and the goal being to relate the branch-point pairs received in a pupil-phase function to the upstream turbulence that created them. This work, in total, could inform the development of future branchpoint-tolerant phase reconstruction algorithms. To increase the fidelity of the results, future developments need to include the effects of additive-sensor noise, low signal-to-noise ratios, and subaperture-sampling requirements. [38][39][40][41][42][43][44][45] With this history in mind, reformulating the problem in terms of slope discrepancy, 46,47 especially when accounting for the effects of speckle due to rough-surface scattering, [48][49][50][51] could also inform these future developments.
This paper builds on this aforementioned history. In particular, it uses wave-optics simulations to investigate the branch-point density as a function of the grid sampling. The goal for these wave-optics simulations is to model plane-wave propagation through homogeneous turbulence, both with and without the effects of a finite inner scale modeled using a Hill spectrum. In practice, the Rytov number, Fried coherence diameter, and isoplanatic angle help to setup and explore the associated deep-turbulence conditions. These parameters provide a gauge for the amount of scintillation, turbulence-limited resolution, and anisoplanatism, respectively, within the wave-optics simulations. On the other hand, the grid sampling provides a gauge for the amount of branch-point resolution within the wave-optics simulations.
Via Monte Carlo averaging, the results show that without the effects of a finite inner scale, the branch-point density grows without bound with adequate grid sampling. Even so, as the innerscale size increases, the results also show that this unbounded growth (1) significantly decreases as the Rytov number, Fried coherence diameter, and isoplanatic angle increase in strength and (2) saturates with adequate grid sampling. These findings are encouraging from the standpoint that they could readily improve the performance of existing branch-point-tolerant phase reconstruction algorithms, as well as inform the development of future algorithms.
It is important to note that this paper also builds on the work contained in a recent conference proceeding. 52 The main difference is that this paper includes the effects of a finite inner scale. Both papers show that the branch-point density grows without bound when neglecting the effects of a finite inner scale. This finding speaks to a preconceived notion within the atmosphericpropagation research community that the branch-point density grows without bound with increasing branch-point resolution. Nonetheless, the results of this paper ultimately show that if one includes the effects of a finite inner scale, then this unbounded growth (1) significantly decreases as the associated deep-turbulence conditions become more pronounced and (2) saturates with adequate branch-point resolution within the wave-optics simulations. These findings are novel and worth sharing with the atmospheric-propagation research community. In general, (1) and (2) imply that future developments need to include the effects of a finite inner scale to accurately model the multifaceted nature of the branch-point problem in adaptive optics.
In what follows, Sec. 2 provides the background details needed to use the branch-point density as a metric of interest. Section 3 provides the wave-optics simulation details needed to setup and explore the associated deep-turbulence conditions, which assume plane-wave propagation through homogeneous turbulence. Results and discussion naturally follow in Sec. 4 with a conclusion in Sec. 5.

Background
In the pupil plane of an optical system, the complex-optical field, U, takes the following phasor form: (1) where A is the pupil-amplitude function and ϕ is the pupil-phase function, such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 7 0 1 ϕ ¼ ArgðUÞ ¼ tan −1 ImðUÞ ReðUÞ : (2) Here, ReðUÞ and ImðUÞ are the real and imaginary components of U. Substituting Eq. (1) into Eq. (2) results in a modulo-2π function that researchers often refer to as the "principle value" or "wrapped phase." When A ¼ 0, the argument of U is indeterminate and therefore a multivalued function. 12 This outcome corresponds to the case where ReðUÞ ¼ ImðUÞ ¼ 0 and is the reason branch points arise within ϕ due to total-destructive interference.
Researchers can determine the location of a branch point using a contour integral around the the gradient of the pupil-phase function, ∇ϕ. 1 Specifically, when the following relationship holds true: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 6 0 ∇ϕ is no longer homogeneous and ϕ is no longer a purely potential field. In Eq.
where , include what Fried refers to as spurious branch points. 1 In practice, spurious branch points are a pair of positive and negative branch points centered on immediately adjacent grid points. This outcome is undesirable because the associated branch cut remains unresolved and could be the result of an artifact within the waveoptics simulations. To address these concerns, one can use a simple search to remove the spurious branch points from the sums in Eq. (4). Empirical evidence says that spurious branch points make up a significant portion of the total number of positive and negative branch points-∼40% on average. Thus, in this paper, we present all results with spurious branch points removed.
With results in mind, the branch-point density, D BP is the metric of interest in this paper. As such E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 2 2 7 where D is the circular-pupil diameter. Simply put, D BP is the number of branch points within the pupil-phase function.

Setup and Exploration
The wave-optic simulations performed in this paper made use of the WavePlex Toolbox for MATLAB. 53 This toolbox uses the split-step beam propagation method (BPM) to simulate the propagation of monochromatic and polychromatic light through the atmosphere. [15][16][17][18][19][20][21][22] In practice, the split-step BPM divides the atmosphere into independent volumes, such that a phase screen represents the atmospheric aberrations present in a volume. Angular-spectrum propagation to each phase screen (from the source plane to the pupil plane) then represents the propagation of light through the atmosphere. With the split-step BPM in mind, this section provides the setup and exploration needed to appreciate the results presented in the next section. Table 1 contains a summary of all the parameters of interest in the wave-optics simulations. The goal for these wave-optics simulations was to model plane-wave propagation through homogeneous turbulence. For this purpose, the wave-optics simulations made use of N × N grids, where N is the grid resolution. The physical side length, S, was the same in both the source and pupil planes, allowing for unity scaling within the setup. With that said, the source plane was setup with a plane wave of unit irradiance and the pupil plane was setup with a circular-pupil diameter, D. For the minimum grid resolution used within the wave-optics simulations, the setup also satisfied critical sampling, 16 such that N min ¼ S 2 ∕ðλZÞ ¼ 512, where N min is the minimum grid resolution, λ is the wavelength, and Z is the propagation distance. The overall setup resulted in a substantial guard band ratio (GBR), where GBR ¼ S∕D. Such a GBR helped in combating the effects of aliasing, 15,16 which we determined to be visually negligible within the wave-optics simulations for all N. The overall setup also resulted in a grid sampling, δðNÞ, where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 5 0 7 δðNÞ ¼ S N :

Setup
As such, δðNÞ provided a gauge for the amount of branch-point resolution within the waveoptics simulations.
While Table 1 contains all the parameters of interest in the wave-optics simulations, Table 2 makes use of several path-integral expressions to define the deep-turbulence conditions. These path-integral expressions were the subject of a recent conference proceeding that discusses the limitations of deep turbulence. 54 For compactness, Table 2 contains a representative subset of these deep-turbulence conditions, which consisted of Rytov numbers ranging from 0.1 to 10.0 in increments of 0.1.
Recall that the Rytov number provides a gauge for the amount of scintillation. Given planewave propagation, 15 the path-integral expression takes the following form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 3 4 8 where C 2 n ðzÞ is the path-dependent refractive index structure coefficient and k ¼ 2π∕λ is the angular wavenumber. With homogeneous turbulence, the path-integral expression reduces to a closed-form expression, where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 2 7 3 Also recall that when the Rytov number increases above 0.1, the scintillation becomes severe and total-destructive interference gives rise to branch points in the pupil-phase function. Thus, when 0.1 < R PW ≤ 0.5, 0.5 < R PW ≤ 1, and R PW > 1, one is left with weak-to-moderate, moderate-to-strong, and strong scintillation conditions, respectively. 54 For completeness in defining the deep-turbulence conditions, Table 2 also includes values for the Fried coherence diameter and the isoplanatic angle. Given plane-wave propagation through homogeneous turbulence, 15 the path-integral expressions reduce to closed-form expressions, such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 2 3 8 for the Fried coherence diameter, and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 1 7 9 for the isoplanatic angle. For all intents and purposes, the Fried coherence diameter helps in parameterizing resolution, whereas the isoplanatic angle helps in parameterizing anisoplanatism. 55 Thus, when D∕r 0;PW > 1 and θ 0 ∕ðλ∕DÞ < 10, one is left with turbulence-limited resolution and anisoplanatic aberrations, respectively. 54 In what follows, we present an exploration of the deep-turbulence conditions provided in Table 2, both with and without the effects of a finite inner scale. In particular, we use the well-known Kolmogorov spectrum when the wave-optics simulations do not include the effects of a finite inner scale. This spectrum is an idealization with an inner-scale size, l 0 , of zero and an outer-scale size, L 0 , of infinity. Such a spectrum does not give an accurate representation of the actual energy distribution in a turbulent volume. 15 This shortcoming can be remedied by using a modified spectrum such as the Hill spectrum. 56 Specifically, this paper makes use of the Hill spectrum because it includes a small rise at high wavenumbers near 1∕l 0 , and it is readily available within the WavePlex Toolbox for MATLAB. Figures 1 and 2 show the numerical log-amplitude variance, σ 2 χ , as a function of the analytical Rytov number, R PW . In particular, Fig. 1 corresponds to the first half of the deep-turbulence conditions setup in Table 2, whereas Fig. 2 corresponds to all of the deep-turbulence conditions setup in Table 2. Both figures include this scintillation-strength exploration for the minimum grid sampling, δðN ¼ 512Þ ¼ δ min , and the maximum grid sampling, δðN ¼ 16;384Þ ¼ δ max , within the wave-optics simulations (cf. Table 1). Both figures also include this scintillation-strength exploration for values of 3.1, 6.2, 12.4, and 24.8 mm for the inner-scale size, l 0 . These values represent finite inner scales equal to one, two, four, and eight times δ min , respectively.

Exploration
Due to the limitations of the Rytov approximation, 57 the strength of the scintillation is often split into two regimes with respect to R PW . Specifically, the weak-scintillation regime occurs when R PW ≤ 0.25, and the strong-scintillation regime occurs when R PW > 0.25. One can clearly see both regimes in Figs. 1 and 2 with R PW ≈ 0.25 serving as an inflection point with respect to σ 2 χ . Beyond this inflection point, σ 2 χ heads into a saturated regime. 58 This saturation process reaches another inflection point when R PW ≈ 1. Beyond this second inflection point, σ 2 χ heads into a supersaturated regime. 58 It is important to note that Fig. 1 highlights the saturated regime when R PW > 0.25, whereas

Results and Discussion
The results presented in this section show that without the effects of a finite inner scale, the branch-point density grows without bound with adequate grid sampling. In addition, the results show that as the inner-scale size increases this unbounded growth (1) significantly decreases as the Rytov number, Fried coherence diameter, and isoplanatic angle increase and (2) saturates with adequate grid sampling. To make these findings manifest, we discuss the relevant trends in the figures that follow.
In support of (1), Figs. 5-7 show the branch-point density, D BP , as a function of the Rytov number, R PW , the circular-pupil diameter relative to the Fried coherence diameter, D∕r 0;PW , and the isoplanatic angle relative to the diffraction-limited half angle, θ 0 ∕ðλ∕DÞ, respectively. The legends in each subplot denote the grid resolutions, N, setup in Table 1. Each subplot also uses the same inner-scale size relative to the minimum grid sampling, l 0 ∕δ min , explored in Figs. 1 and 2. In Figs. 5-7, the plotted lines represent the Monte Carlo averages associated with 100 turbulence realizations, and the error bars represent the standard deviations. Similar to Figs. 1 and 2, the width of the error bars are small (at most 4% of the mean) and thus we believe that 100 Monte Carlo realizations are adequate in quantifying the behavior of D BP in terms of R PW , D∕r 0;PW , and θ 0 ∕ðλ∕DÞ. In Fig. 5(a), where l 0 ∕δ min ¼ 0, we see that as R PW increases, D BP increases linearly without bound when N > 2048. Figures 5(b)-5(e), on the other hand, show that when l 0 increases, this unbounded linear growth significantly decreases. In particular, when l 0 ∕δ min > 2 and N > 2048, the results for D BP tend to overlap, whereas when l 0 ∕δ min ≤ 2 and N ≤ 2048, the results for D BP tend to saturate.
In Fig. 6(a), where l 0 ∕δ min ¼ 0, we see that as D∕r 0;PW increases, D BP increases exponentially without bound when N > 2048. Figures 6(b)-6(e), on the other hand, show that when l 0 increases, this unbounded exponential growth significantly decreases. In particular, when l 0 ∕δ min > 2 and N > 2048, the results for D BP tend to overlap, whereas when l 0 ∕δ min ≤ 2 and N ≤ 2048, the results for D BP tend to saturate.
In Fig. 7(a), where l 0 ∕δ min ¼ 0, we see that as θ 0 ∕ðλ∕DÞ decreases, D BP increases exponentially without bound when N > 2048. Figures 7(b)-7(e), on the other hand, show that when l 0 increases, this unbounded exponential growth significantly decreases. In particular, when l 0 ∕δ min > 2 and N > 2048, the results for D BP tend to overlap, whereas when l 0 ∕δ min ≤ 2 and N ≤ 2048, the results for D BP tend to saturate.
The aforementioned saturation in Figs. 5-7 when N ≤ 2048 is due to inadequate sampling of the turbulent volume and is not physical. One rule of thumb says that the Fried coherence diameter relative to the grid sampling needs to be greater than ten (i.e., r 0;PW ∕δðNÞ > 10). Simply put, when N ≤ 2048, we do not always satisfy this rule of thumb within the wave-optics simulations, as shown in Figs. 5-7.  In support of (2), Fig. 8 shows the branch-point density, D BP , as a function of the grid resolution, N. The legends in each subplot denote the Rytov numbers, R PW , setup in Table 2. Each subplot also uses the same inner-scale size relative to the minimum grid sampling, l 0 ∕δ min , explored in Figs. 5-7. In Fig. 8, the plotted lines represent the Monte Carlo averages associated with 100 turbulence realizations, and the error bars represent the standard deviations. Similar to Figs. 5-7, the width of the error bars are small (at most 4% of the mean) and thus we believe that 100 Monte Carlo realizations are adequate in quantifying the behavior of D BP in terms of N.
In Fig. 8(a), where l 0 ∕δ min ¼ 0, we see that as N increases, D BP increases without bound for all R PW . Figures 8(b)-8(e), on the other hand, show that when l 0 ∕δ min increases, this unbounded growth saturates. In particular, when l 0 ∕δ min > 0 and N > 2048, the results for D BP tend to roll over, whereas when l 0 ∕δ min ≤ 2 and N ≤ 2048, the results for D BP tend to grow monotonically. This monotonic growth increases as R PW increases, which again is most likely caused by insignificant sampling of the turbulent volume. Figure 8 ultimately shows that as l 0 ∕δ min increases, D BP increasingly saturates when N > 2048. This result disagrees with a preconceived notion within the atmospheric-propagation community that the branch-point density grows without bound with increasing branch-point resolution. As shown in Fig. 8(a), this preconceived notion is the result of using the well-known Kolmogorov spectrum within the wave-optics simulations. Recall that this spectrum is an idealization with an inner-scale size, l 0 , of zero and an outer-scale size, L 0 , of infinity. Such a spectrum leads to energy distribution beyond the inertial sub range, which is not physical. Also recall that this shortcoming can be remedied using a modified spectrum such as the Hill spectrum. 56 In so doing, the wave-optics simulations more accurately model the energy distribution in a turbulent volume and lead to the results presented in Figs

Conclusions
This paper used wave-optics simulations to investigate the branch-point density in terms the grid sampling. The goal for these wave-optics simulations was to model plane-wave propagation through homogeneous turbulence, both with and without the effects of a finite inner scale modeled using a Hill spectrum. In practice, the Rytov number, Fried coherence diameter, and isoplanatic angle provided parameters to setup and explore the associated deep-turbulence conditions within the wave-optics simulations. The grid sampling, on the other hand, provided a gauge for the amount of branch-point resolution within the wave-optics simulations. Via Monte Carlo averaging, the results showed that without the effects of a finite inner scale, the branchpoint density grew without bound. Nevertheless, the results also showed that as the inner-scale size increased, this unbounded growth (1) significantly decreased as the associated deep-turbulence conditions became more pronounced and (2) saturated with adequate branch-point resolution within the wave-optics simulations.
The results of this paper imply that future developments need to include the effects of a finite inner scale to accurately model the multifaceted nature of the branch-point problem in adaptive optics. Recall that this problem tends to be the "Achilles' heel" to beam-control systems that perform deep-turbulence phase compensation. Thus, the results of this paper are encouraging from the standpoint that they could readily improve the performance of existing branch-pointtolerant phase reconstruction algorithms. The results of this paper could also inform the development of future branch-point-tolerant phase reconstruction algorithms, which take into account the effects of additive-sensor noise, low signal-to-noise ratios, and subaperture-sampling requirements. These future developments are a critical next step to improving the performance of beamcontrol systems that perform deep-turbulence phase compensation.