Wave-optics investigation of turbulence thermal blooming interaction: I. Using steady-state simulations

Abstract. Part I of this two-part paper uses wave-optics simulations to look at the Monte Carlo averages associated with turbulence and steady-state thermal blooming (SSTB). The goal is to investigate turbulence thermal blooming interaction (TTBI). At wavelengths near 1  μm, TTBI increases the amount of constructive and destructive interference (i.e., scintillation) that results from high-power laser beam propagation through distributed-volume atmospheric aberrations. As a result, we use the spherical-wave Rytov number and the distortion number to gauge the strength of the simulated turbulence and SSTB. These parameters simplify greatly given propagation paths with constant atmospheric conditions. In addition, we use the log-amplitude variance and the branch-point density to quantify the effects of TTBI. These metrics result from a point-source beacon being backpropagated from the target plane to the source plane through the simulated turbulence and SSTB. Overall, the results show that the log-amplitude variance and branch-point density increase significantly due to TTBI. This outcome poses a major problem for beam-control systems that perform phase compensation.


Introduction
At wavelengths near 1 μm, the effects of turbulence are often more dominant than the effects of thermal blooming. Thus, at short-exposure time scales, one does not typically see fully formed crescent or half-moon irradiance patterns, as described in Fig. 1, due to thermal blooming at wavelengths near 1 μm. From a historical perspective, however, fully formed crescent or half-moon irradiance patterns do repeatedly result from the dominant effects of thermal blooming at wavelengths in the mid-wave and long-wave infrared. [1][2][3][4][5] There is, nevertheless, an interaction that occurs between turbulence and thermal blooming at wavelengths near 1 μm. [6][7][8][9][10] This so-called turbulence thermal blooming interaction (TTBI), in practice, results in an increased amount of scintillation, which is the constructive and destructive interference that results from propagating high-power laser beams through distributed-volume atmospheric aberrations or "deep turbulence." In general, the scintillation caused by turbulence results in localized hot spots that cause localized heating of the atmosphere. This localized heating results in localized defocus-like optical effects (aka localized thermal blooming), which causes more constructive and destructive interference upon propagation through the atmosphere. Given weak-turbulence conditions, for example, the resulting scintillation caused by TTBI is often analogous to that experienced with deep-turbulence conditions. Therefore, thermal blooming can have a major impact on system performance (even at wavelengths near 1 μm).
*Address all correspondence to Mark F. Spencer, E-mail: mark.spencer@osamember.org Beam-control (BC) systems, in theory, can mitigate the nonlinear optical effects induced by thermal blooming. 11 However, when one uses a single deformable mirror (DM) for phase-only compensation, analysis predicts the possibility of an instability due to positive feedback in the BC system. [12][13][14][15][16] Appropriately termed phase compensation instability (PCI), the positive feedback arises with the time-dependent development of scintillation within the propagating highpower laser beam. Recall that the localized hot spots produce defocus-like optical effects in the atmosphere. A BC system corrects for the hot spots by applying focus-like phase compensations. In turn, these phase compensations increase the strength of the localized thermal blooming, which leads to a runaway condition until something mitigates the positive feedback (e.g., wind variations in the atmosphere [17][18][19][20][21] or branch points in the phase function [22][23][24][25][26]. Whether from TTBI or PCI, an increase in the amount of scintillation leads to an increase in the amount of total destructive interference, which leads to amplitude nulls in both the real and imaginary parts of the complex optical field. These amplitude nulls cause branch points and branch cuts (i.e., 2π discontinuities) to arise in the collimated phase function of a backpropagated beacon. 27,28 This outcome poses a major problem for BC systems that perform phase compensation. In practice, these branch points and branch cuts lead to fitting error in a BC system that uses a single continuous-face-sheet DM with a high-power coating for phase-only compensation. [29][30][31][32][33] With this last point in mind, the atmospheric propagation research community needs to quantify the impacts of TTBI and PCI using the larger grid sizes allowed by modern-day, waveoptics simulations.
The following is Part I of a two-part paper on TTBI performed using AOTools and WavePropboth of which are MATLAB toolboxes written by the Optical Sciences Company. [34][35][36][37][38] In turn, this paper investigates TTBI in the presence of turbulence and SSTB. It does so via the Monte Carlo averages associated with the log-amplitude variance and the branch-point density. These metrics result from a point-source beacon being backpropagated from the target plane to the source plane through the simulated turbulence and SSTB. Part II of this two-part paper then investigates TTBI in the presence of turbulence and time-dependent thermal blooming (TDTB). Together, these papers will inform future wave-optics investigations.
It is important to note that this paper stands on its own as an independent article. However, it is also important to note that Part II of this two-part paper complements the analysis contained in this paper. In general, the steady-state assumptions contained in this paper allows us to explore the trade space in a computationally efficient manner, whereas the time-dependent assumptions contained in Part II allow us to examine the Monte Carlo averages with increased computational fidelity. For example, relative to the time-dependent results contained in Part II, the steady-state results contained in this paper provide an upper bound on both the increase in the log-amplitude Fig. 1 Thermal blooming is a nonlinear optical effect caused by the irradiance (i.e., the power per unit area) of a high-power laser beam being absorbed by molecules and aerosols in the atmosphere. This absorbed irradiance leads to an increase in the temperature of the air and a decrease in the refractive index, creating a defocus-like optical effect that blooms the beam outward and decreases its peak irradiance. The presence of a transverse wind then contributes a tilt-like optical effect, which shifts the peak irradiance off target. After the wind has cleared through the source plane at least once, the conditions for static, whole-beam, or steady-state thermal blooming (SSTB) are generally met. The end result is a fully formed crescent or half-moon irradiance pattern in the target plane.
variance and the branch-point density due to TTBI. Such results will hopefully prove fruitful in the development of next-generation scaling laws that account for the effects of TTBI. [39][40][41][42][43][44][45] In what follows, Sec. 2 contains the setup for this paper, whereas Sec. 3 explores the trade space. Results and discussion then follow in Sec. 4, with a conclusion in Sec. 5. Before moving on to the next section, it is worth mentioning that this paper builds upon the preliminary analysis presented by Murphy and Spencer in a recent conference proceeding. 46 In particular, this paper reduces the overall trade space and uses computationally efficient steady-state simulations to clearly show an increase in both the log-amplitude variance and branch-point density due to TTBI. These results serve as a novel contribution to the atmospheric propagation research community.

Setup for the Wave-Optics Simulations
The desired setup is as follows. We wish to propagate a focused high-power laser beam with a wavelength of 1 μm along a propagation path with constant atmospheric conditions. Afterward, we wish to propagate a point-source beacon of the same wavelength back along the same path through the simulated turbulence and SSTB. For this purpose, we quickly review the details associated with the split-step beam propagation method (BPM), spherical-wave Rytov number, SSTB, distortion number, and parameters of interest in the following sections. Given a common setup, this paper shares most of these sections with Part II of this two-part paper on TTBI. The analysis in this paper specifically differs in the sections on SSTB and parameters of interest. With that said, we include the shared sections in both papers for reading independence.

Split-Step Beam Propagation Method
In this paper, we make use of the commonly used split-step BPM, which simulates the propagation of monochromatic and polychromatic light through the atmosphere. [47][48][49][50] As described by Schmidt, 51 the split-step BPM divides the propagation path into independent volumes, with the atmospheric aberrations in each volume being represented by a single phase screen. The splitstep BPM makes use of angular spectrum or plane-wave spectrum propagation to vacuumpropagate the light from a source plane to the first phase screen, applies the phase screen, and repeats this process until the monochromatic light reaches a target plane. Before moving on in the analysis, it is worth mentioning that AOTools and WaveProp make use of the split-step BPM. It is also worth mentioning that AOTools and WaveProp generate the phase screens associated with turbulence similar to the approach presented in Chapter 9 of Ref. 51. In this paper, we specifically use a Kolmogorov power-spectral density and do not add subharmonics or additional tilt to the generated phase screens.

Spherical-Wave Rytov Number
Provided the Rytov approximation, the propagation of a spherical wave through turbulence has an associated path-integral expression that serves as a gauge for the amount of scintillation. 41,[51][52][53][54] Known as the spherical-wave Rytov number R sw (aka the spherical-wave Rytov parameter or spherical-wave, log-amplitude variance), this 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 1 ; 1 1 6 ; 1 8 7 where k ¼ 2π∕λ is the angular wavenumber, Z is the propagation distance, and C 2 n ðzÞ is the path-dependent refractive index structure coefficient. Given propagation paths with constant atmospheric conditions, this path-integral expression reduces to the following closed-form expression: 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 ; 9 3 R sw ¼ 0.124k 7∕6 C 2 n Z 11∕6 : For all intents and purposes, the strong scintillation regime occurs when R sw is greater than 0.25. This regime is where the Rytov approximation readily breaks down with respect to the logamplitude fluctuations 52,53 and branch points and branch cuts (i.e., 2π discontinuities) readily show up in the phase function. 27,28 Thus, in the analysis that follows, we will use R sw as a gauge for the strength of the simulated turbulence.

Steady-State Thermal Blooming
For SSTB, we assume that a transverse wind has cleared the source plane enough times so that the TDTB reaches a steady state with a fully formed crescent or half-moon irradiance pattern in the target plane (cf. Fig. 1). The heating and cooling of the air, as described by TDTB (cf. Part II), can reach a steady state if there are no variations in the exitance of the high-power laser source and ambient atmosphere. In turn, SSTB results in a static behavior of the change in the refractive index. This static behavior causes the time rate of change of the change in the refractive index to be zero. Consequently, we drop any time dependence, and the position-dependent change in the refractive index ΔnðrÞ due to SSTB results as 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 4 3 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 4 ; 1 1 6 ; 4 8 6 is the path-dependent absorbed irradiance coefficient, n 0 ðzÞ is the path-dependent ambient refractive index, αðzÞ is the path-dependent absorption coefficient, C P ðzÞ is the path-dependent specific heat at constant pressure, ρ 0 ðzÞ is the path-dependent density of air at constant pressure, T 0 ðzÞ is the path-dependent ambient temperature, v ⊥ ðzÞ is the path-dependent transverse wind speed, and i HP ðrÞ is the position-dependent, high-power laser beam irradiance. In writing Eqs. (3) and (4), we assume that the transverse wind velocity is solely in the x direction, hence the limits of integration. In so doing, we can account for the effects of SSTB using AOTools and WaveProp since these MATLAB toolboxes solve the integral expression found in Eq. (3) numerically. 34-38

Distortion Number
It is useful to describe the refraction caused by SSTB in terms of a phase error ϕðx; y; ZÞ measured in radians. For this purpose, 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 7 2 ϕðx; y; ZÞ ¼ k where again k ¼ 2π∕λ is the angular wavenumber and Z is the propagation distance. Substituting Eq. (3) into Eq. (5) results in the following relationship: 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 ; 2 0 5 ϕðx; y; ZÞ ¼ −k As shown in Eq. (6), we can characterize the radians of distortion induced by SSTB.
For focused high-power laser beams with initial power, P 0 , and initial diameter, D 0 , we can rewrite Eq. (6) in terms of a path-integral expression known as the distortion number, N D . In particular, 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 ; 1 1 2 and 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 ; 7 2 3 From the source plane at z ¼ 0 to the target plane at z ¼ Z, the optical leverage ðZ − zÞ∕Z causes D 0 to converge upon propagation so 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 ; 6 5 5 is the path-dependent beam diameter from geometrical or ray optics. Note that we include the path-dependent transmittance τðzÞ, 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 1 0 ; 1 1 6 ; 5 9 0 in the definition of N D because of extinction effects from the path-dependent extinction coefficient κðzÞ ¼ αðzÞ þ σðzÞ, where σðzÞ is the path-dependent scattering coefficient. In turn, N D provides a gauge for the radians of distortion induced by SSTB. Given propagation paths with constant atmospheric conditions, Eq. (8) simplifies, 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 1 1 ; 1 1 6 ; 4 9 7 For all intents and purposes, the strong distortion regime occurs when the distortion number N D is greater than the critical number 34 This regime is where the radians of distortion induced by SSTB give rise to significant fluctuations in the fully formed crescent or half-moon irradiance pattern (cf. Fig. 1). Thus, in the analysis that follows we will use N D as a gauge for the strength of the simulated SSTB. Table 1 contains all the parameters of interest in the wave-optics simulations. It is important to note that the wave-optics simulations used N × N grids. The side length S was the same in both the source and target planes creating unity scaling within the wave-optics simulations. By choice, the wave-optics simulations also satisfied Fresnel scaling, such that N ¼ S 2 ∕ðλZÞ, where λ is the wavelength and Z is the propagation distance. The resulting N × N grid also minimized the effects of aliasing without making the wave-optics simulations too computationally expensive.

Parameters of Interest
In addition to Fresnel scaling, we ensured that the wave-optics simulations had at least 10 pixels across the spherical-wave Fried parameter r 0;sw (also known as the spherical-wave coherence diameter or spherical-wave coherence length). 41,51-54 Given a propagation path with constant atmospheric conditions, the associated path-integral expression simplifies into the following closed-form expression: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 2 1 9 Table 2 makes use of this closed-form expression in defining the turbulence scenarios used in the wave-optics simulations. For completeness in defining the turbulence scenarios, 41,51-54 it also makes use of the isoplanatic angle, θ 0 , and the Greenwood frequency, f G , 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 1 3 ; 1 1 6 ; 1 3 8 Notice that the scenarios defined in Table 2 include spherical-wave Rytov numbers R sw from 0.1 to 1.0, both below and above the strong-scintillation regime (when R sw > 0. 25). Also notice that the scenarios defined in Table 1 include distortion numbers N D from 8 ffiffi ffi 2 p to 32 ffiffi ffi 2 p , both below and above the strong distortion regime (when N D > N C ¼ 16 ffiffi ffi 2 p ≈ 22.6). With these regimes in mind, we will explore the overall trade space in the next section. Table 2 Turbulence scenarios used in the wave-optics simulations. Recall that C 2 n is the refractive index structure coefficient, R sw is the spherical-wave Rytov number, D 0 is the aperture diameter, r 0;sw is the spherical-wave Fried parameter, θ 0 is the isoplanatic angle, λ∕D 0 is the diffraction-limited angle, and f G is the Greenwood frequency.
Specific heat at constant pressure (J/kg/K) C P 1004 Density of air at constant pressure (kg∕m 3 ) ρ 0 1.293 Ambient temperature (K) T 0 300 Transverse wind speed (m/s) v ⊥ 5

Exploration Using the Wave-Optics Simulations
In this section, we make use of the wave-optics simulations setup in the previous section to explore the trade space. The goal is to investigate TTBI in terms of the normalized power in the bucket (PIB), P N , and the peak Strehl ratio, S P , associated with a focused high-power laser beam being propagated from the source plane at z ¼ 0 to the target plane at z ¼ Z, and the log-amplitude variance, σ 2 χ , and the branch-point density, D BP , associated with a point-source beacon being backpropagated from the target plane at z ¼ Z to the source plane z ¼ 0. Here, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 6 4 3 In Eqs. (15)-(18), cylðρÞ is the cylinder 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 1 9 ; 1 1 6 ; 4 5 8 cylðρÞ ¼ D Z ¼ 2.44λZ∕D 0 is the diffraction-limited bucket diameter; i HP ðx; y; ZÞ is the focused highpower laser beam irradiance in the target plane; P 0 is again the initial power; maxf⊙g is an operator that computes the maximum value; i DL ð0;0; ZÞ is the focused on-axis, diffractionlimited, high-power laser beam irradiance in the target plane; varf⊙g is an operator that computes the spatial variance; χ PS ðx; y; 0Þ ¼ lnfjU PS ðx; y; 0Þjg is the log amplitude of the backpropagated point-source beacon in the source plane; N BP is the number of branch points in the collimated phase function of the backpropagated point-source beacon; and D 0 is again the initial diameter.
To calculate the number of branch points N BP in Eq. (18), we used the following relationship: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 2 9 2 I C ∇ϕðx; y; 0Þ · dr ¼ AE2πðN þ − N − Þ; (20) where N þ is the number of positively charged branch points and N − is the number of negatively charged branch points within the collimated phase function ϕðx; y; 0Þ. This relationship says that we can determine the location of a branch point when the line integral around the closed curve C of the gradient of the collimated phase function ∇ϕðx; y; 0Þ does not equal zero, specifically, where ∇ϕðx; y; 0Þ is a nonconservative vector field. To account for this relationship numerically, AOTools and WaveProp discretely samples the continuous integral in Eq. (20) by breaking the N × N grid into a bunch of 2 × 2 subgrids and summing up the phase derivative around each point. A positive 2π value results in a positively charged branch point and a negative 2π value results in a negatively charged branch point. In turn, to calculate N BP , we computed the sum of the total number of positive and negative branch points associated with the pixels found within ϕðx; y; 0Þ.
In what follows, we will use the metrics defined in Eqs. (15)- (18) to visualize the following topics: 1. the number of phase screens needed to simulate the effects of SSTB, 2. the focused high-power laser beam in the target plane, 3. the backpropagated point-source beacon in the source plane, and 4. the overall trade space in terms of the normalized PIB, P N , and the peak Strehl ratio, S P .
Each topic is the subject of the following sections. These sections will inform the results and discussion presented in Sec. 4.

Number of Phase Screens Needed
To determine the number of equally spaced phase screens, N PS , needed to accurately simulate SSTB using the split-step BPM, we calculated both the normalized PIB, P N , and the peak Strehl ratio, S P , as a function of N PS . In addition, we calculated the rate of change of P N and S P as a function of N PS via an operator, difff⊙g, that determines the running difference. Figure 2 shows the outcomes of these calculations. As shown, we can see that the rate of change of P N and S P go to zero (accurate to the third decimal place) when N PS ≥ 50. Thus, in the analysis that follows, we used N PS ¼ 50 to simulate the effects of SSTB.
For convenience in the wave-optics simulations, we also used 50 equally spaced phase screens to simulate the effects of turbulence. Therefore, in the analysis that follows, the phase screens used for simulating turbulence and SSTB (using the split-step BPM) were collocated along the propagation path. This choice led to small percentage errors (less than a 10th of a percentage) between the continuous and discrete calculations of the parameters found in Table 2. 51

Focused High-Power Laser Beam
To create the focused high-power laser beam, we used a series of steps starting with the creation of a positive thin lens transmittance function of circular diameter D 0 and focus Z. Assuming plane-wave illumination, we then set the exitance i 0 of the focused high-power laser beam, such that i 0 ¼ 4P 0 ∕ðπD 2 0 Þ. This series of steps created a top-hat or flat-top beam profile in the source plane with approximately 256 grid points across D 0 . As discussed above, AOTools and WaveProp then used the split-step BPM to propagate the focused high-power laser beam from the source plane to the target plane.
Based on previously published theoretical explorations, 6-10 we hypothesized that the simulated turbulence and SSTB would increase the amount of scintillation found in a propagated high-power laser beam due to TTBI. In turn, the normalized PIB P N and peak Strehl ratio (a) (b) Fig. 2 Visualization of the number of phase screens N PS needed to simulate the effects of SSTB: (a) the normalized PIB P N and (b) the peak Strehl ratio S P (blue curves). Note that the rate of change of P N and S P goes to zero around N PS ¼ 50 (orange curves).
S P for the TTBI case would typically decrease in comparison to the simulated diffractionlimited, turbulence-only, and SSTB-only cases. Figure 3 demonstrates this hypothesis to be true for one Monte Carlo realization, where R sw ¼ 0.25 and Eqs. (2) and (11), respectively]. Notice that we report values for both P N and S P at the top of each normalized irradiance subplot in Fig. 3. Also notice that these irradiance subplots are similar but different to those reported in Part II using time-dependent simulations (cf. Fig. 3 in Part II) since both papers use the same Monte Carlo realization of turbulence.

Backpropagated Point-Source Beacon
To create the backpropagated point-source beacon, AOTools and WaveProp used a series of steps starting with the creation of a positive thin lens transmittance function of square width 2D 0 and focus Z. Assuming Fresnel scaling, AOTools and WaveProp vacuum-propagated this positive thin lens transmittance function from the source plane to the target plane using angular spectrum or plane-wave spectrum propagation. This series of steps created a sinc-like function in the target plane with three pixels across its central lobe. As discussed above, AOTools and WaveProp then used the split-step BPM to backpropagate the point-source beacon from the target plane to the source plane. (c) (d) Fig. 3 Visualization of the focused high-power laser beam in the target plane. Here, we plot the normalized irradiance associated with (a) the diffraction-limited case, (b) the turbulence-only case, (c) the SSTB-only case, and (d) the TTBI case. We also report the normalized PIB P N and the peak Strehl ratio S P at the top of each subplot. Note that the white circles represent the diffractionlimited bucket diameter D Z .
Again, based on previously published theoretical explorations, 6-10 we hypothesized that the simulated turbulence and SSTB would increase the amount of scintillation found in a backpropagated point-source beacon due to TTBI. As a result, the log-amplitude variance σ 2 χ and branchpoint density D BP for the TTBI case would typically increase in comparison to the simulated diffraction-limited, turbulence-only, and SSTB-only cases. Figures 4 and 5 demonstrate this hypothesis to be true for one Monte Carlo realization, where R sw ¼ 0.25 and N D ¼ Eqs. (2) and (11), respectively]. Notice that we report values for both σ 2 χ and D BP at the top of each normalized irradiance subplot in Fig. 4 and each wrapped phase subplot in Fig. 5. Also notice that these irradiance subplots are similar but different to those reported in Part II using time-dependent simulations (cf. Figs. 4 and 5 in Part II) since both papers use the same Monte Carlo realization of turbulence.
To further explore the effects of TTBI, we also calculated the irrotational phase estimate, ϕ irr ðx; y; 0Þ, and rotational phase estimate,φ rot ðx; y; 0Þ. Figure 6 shows these estimates. In practice, these estimates originated from the collimated phase functions ϕðx; y; 0Þ found in Figs. 5(b) and 5(d), for the turbulence-only and TTBI cases, respectively. To calculateφ irr ðx; y; 0Þ and ϕ rot ðx; y; 0Þ from ϕðx; y; 0Þ, we used the following relationships: 27 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 6 ; 2 8 3φ rot ðx; y; 0Þ ¼ argfexpðj½ϕðx; y; 0Þ − ϕ LS ðx; y; 0ÞÞg; (22) where LSf⊙g is an operator that unwraps the phase using a least-squares algorithm and argf⊙g is an operator that extracts the wrapped phase ϕ from a phasor of the form exp½jϕ. Figure 6 shows that both the irrotational phase estimate,φ irr ðx; y; 0Þ, and the rotational phase estimate,φ rot ðx; y; 0Þ, contain more phase features due to the TTBI case [cf. Figs. 6(b) and 6(d)] compared to the turbulence-only case [cf. Figs. 6(a) and 6(c)]. This outcome says that a BC system, in general, will need more resolution and stroke to correct forφ irr ðx; y; 0Þ. It also says thatφ rot ðx; y; 0Þ will lead to an increased amount of fitting error because of the 2π discontinuities associated with branch points and branch cuts, especially when using a single continuousface-sheet DM with a high-power coating for phase-only compensation. [29][30][31][32][33] 3.4 Overall Trade Space Figure 7 shows the overall trade space in terms of the normalized PIB, P N , and the peak Strehl ratio, S P . Here, we calculated both metrics as a function of the spherical-wave Rytov number, R sw , and the distortion number, N D . As both R sw and N D increase, both P N and S P typically (b) (a) Fig. 7 Visualization of the overall trade space in terms of the normalized PIB P N and the peak Strehl ratio S P : (a) P N and (b) S P . Note that we plot both metrics as a function of the sphericalwave Rytov number, R sw , and the distortion number, N D . Also note the curves represent the averages and the error bars represent the standard deviations associated with 100 Monte Carlo realizations.
decrease, which makes good sense. Note that the curves report the averages, and the error bars report the standard deviations associated with 100 Monte Carlo realizations. Also note that the widths of these error bars are reasonably small, and the curves themselves are reasonably smooth. Thus, we will use the same number of Monte Carlo realizations in the analysis that follows.

Results and Discussion
This section contains results for the trade space setup in Sec. 2 and explored in Sec. 3. In particular, Fig. 8 shows results for the log-amplitude variance, σ 2 χ , and the branch-point density, D BP , both as a function of the spherical-wave Rytov number, R sw , and the distortion number, N D [cf. Eqs. (17), (18), (2), and (11), respectively]. These metrics result from a point-source beacon being backpropagated from the target plane to the source plane through the simulated turbulence and SSTB. Here, we use R sw to gauge the strength of the simulated turbulence and N D to gauge the strength of the simulated SSTB. From Fig. 8(a), we can see that σ 2 χ deviates significantly from R sw when SSTB is present, even for the lowest N D . Similarly, Fig. 8(b) shows that D BP increases significantly when SSTB increases in strength (i.e., as N D increases). Note that in Fig. 8 the curves report the averages and the error bars report the standard deviations associated with 100 Monte Carlo realizations. Also note that the widths of these error bars are reasonably small, and the curves themselves are reasonably smooth. Thus, we believe that 100 Monte Carlo realizations (i.e., the same number of realizations used in Part II of this two-part study) are adequate in quantifying the effects of TTBI.
These results clearly show that TTBI results in an increased amount of scintillation when simulating turbulence and SSTB. For example, as shown in Fig. 8(a), the log-amplitude variance σ 2 χ increases monotonically as a function of the spherical-wave Rytov number R sw . However, as the distortion number N D increases, this monotonic increase becomes increasingly saturated (i.e., it approaches an asymptotic value at a faster rate) due to TTBI. For all intents and purposes, this behavior is well documented in the atmospheric propagation literature 52 and is a known shortcoming of the Rytov approximation (cf. the solid blue curve relative to the solid black curve). What is not well documented in the atmospheric propagation literature is the fact that TTBI causes this saturation to occur at a faster rate.
The branch-point density D BP also increases monotonically, as shown in Fig. 8(b), but in general, the number of branch points increases nonlinearly up until the strong scintillation regime (cf. the solid blue curve relative to the solid black curve). After this regime is met, the number of branch points grows linearly. This behavior is also well documented in the atmospheric propagation literature 55 and poses a major problem for BC systems that perform phase compensation in the strong-scintillation regime. Again, what is not well documented is the fact that TTBI causes this nonlinear-to-linear transition to occur at a faster rate.
With the above behaviors in mind, the steady-state assumptions contained in this paper ultimately allowed us to explore the full trade space contained in Table 2 with computational efficiency. However, this computational efficiency came at the expense of computational fidelity. That is why we reserved the detailed examination of the Monte Carlo averages associated with the log-amplitude variance, σ 2 χ , and the branch-point density, D BP , to Part II of this two-part study. Relative to the time-dependent results contained in Part II, the steady-state results contained in this paper provide an upper bound on both the increase in σ 2 χ and D BP due to TTBI (cf. Fig. 8 in both papers). It is our hope that such results will prove fruitful in the development of next-generation scaling laws that account for the effects of TTBI. [39][40][41][42][43][44][45]

Conclusion
In this paper, we used wave-optics simulations to look at the Monte Carlo averages associated with turbulence and SSTB. The goal throughout was to investigate TTBI. At a wavelength near 1 μm, TTBI increases the amount of scintillation that results from high-power laser beam propagation through distributed-volume atmospheric aberrations. In turn, to help gauge the strength of the simulated turbulence and SSTB, this paper made use of the following two parameters: the spherical-wave Rytov number and the distortion number. These parameters simplified greatly, given a propagation path with constant atmospheric conditions. In addition, to help quantify the effects of TTBI, this paper made use of the following two metrics: the log-amplitude variance and branch-point density. These metrics resulted from a point-source beacon being backpropagated from the target plane to the source plane through the simulated turbulence and SSTB.
Overall, the results showed that TTBI causes the log-amplitude variance and the branch-point density to increase. These results pose a major problem for BC systems that perform phase compensation. In turn, the steady-state simulations presented in this paper will provide much needed insight into the design of future systems.