Investigation of turbulence-tolerant free-space optical communications via multiplane light conversion

Abstract. Multiplane light conversion (MPLC) offers an alternative to adaptive optics for coupling a turbulence-corrupted free-space optical beam into single-mode fiber or waveguide. Recent published test results suggest that such a conversion device has comparable or better performance than an adaptive optics system. To better understand the device characteristics, simulations were conducted to quantify the power loss for different turbulence strengths and the number of Hermite–Gaussian modes used in the conversion process. The specific case studied is a prototype free-space laser communication system developed by the U.S. Army Research Laboratory. Simulation and statistical results of the conversion performance are reported. Analysis of a beam power combiner after MPLC is also discussed.


Introduction
Multiplane light conversion (MPLC) devices perform unitary spatial transformation 1-3 of optical beams composed of a set of orthogonal modes, such as Hermite-Gaussian (HG), to a different set of orthogonal modes via specially designed phase plates. MPLC covers a range of applications, such as spatial mode division (de)multiplexing for high-capacity optical fiber communications, 4 light-wave coherent detection, 5 and optical phase-front precompensation for communications and directed energy. 6 For free-space optical communications (FSOCs), MPLC decomposes turbulence-corrupted optical beam into a relatively small number (e.g., N ¼ 10) of low-order HG modes that are subsequently transformed into N spatially separated fundamental Gaussian modes for coupling into N single-mode fibers. Guided Gaussian beams in these fibers are then coherently combined to recover the transmitted optical signal. 7 An example of a simplified FSOC receiver schematic with a four-output MPLC followed by a 4-to-1 power combiner is shown in Fig. 1. Beam pointing and tracking devices, such as fast-steering mirrors, for maintaining optical alignment to the MPLC are not shown, but their presence before the MPLC as a part of the FSOC receiver architecture is assumed. Combining the N outputs into one singlemode fiber or waveguide eliminates the complexity of implementing multiple or N photoreceivers, especially for large N. Also, the likelihood of some of these photoreceivers becoming thermal-noise-limited is expected to be higher than the detection of the combined optical signal using a single photoreceiver. Note that the four-output MPLC and the 4-to-1 power combiner shown in Fig. 1 are for clarity and illustrative purpose only. In general, N-output MPLC and N-to-1 power combiner configurations with N ¼ 10, 15, 21, or higher are typical. In this report, the cases of N ¼ 10, 15, and 21 were investigated.
Current commercially available MPLC devices are based on microoptics. 8 Recently published FSOC test results suggest that the MPLC device has a comparable or better performance than adaptive optics. 7 Adaptive optics systems mitigate turbulence-induced beam distortion and allow free-space coupling into single-mode fiber. These systems have high performance but *Address all correspondence to Pak Cho, pak.s.cho.ctr@army.mil typically require complex subsystems, such as a wave-front sensor and two-dimensional (2D) deformable reflector array. As a result of their complexity, relatively large size, and cost, they are not commonly employed in military FSOC applications that require mobility and low size, weight, and power (SWaP). MPLC, which is passive with no moving components and has a comparatively small footprint, is an alternative to adaptive optics that potentially fulfill military requirements for FSOC.
In MPLC, the turbulence-corrupted optical beam can be expressed as a superposition or summation of many orthogonal HG spatial modes with corresponding complex coefficients c mn , as described in Sec. 7.1 in the Appendix. The transformation of HG into Gaussian modes in MPLC was achieved via free-space propagation to a series of phase plates. These specially designed phase plates impart phase-front perturbations to the input beam. As the beam propagates in the free-space region between phase plates these phase-front perturbations are converted to spatial amplitude perturbations via an equivalent 2D Fourier transform, thereby shaping the beam intensity profile. Alternating between successive phase-front perturbation induced by the phase plate and free-space propagation eventually transforms the HG mode into the desired mode, such as Gaussian beam. The phase plates profiles, number of propagations, and distance between propagations require careful designs to achieve the desired output beam profiles and distribution, 1,5 e.g., N spatially separated Gaussian beams. The MPLC described in Ref. 1, for example, consists of a single extended reflective phase plate and a spherical mirror in a multipass cavity configuration that performs simultaneous spatial unitary transformation of three HG modes into three spatially separated Gaussian beams that couple into single-mode fibers array, or vice versa. Conversion between a larger number of spatially separated Gaussian and HG beams or arbitrary input beam shapes has been demonstrated. 4,5,7 Typically, small intensity features or high-spatial frequencies in the turbulence-corrupted optical beam require more higher-order HG modes to reproduce or recover the energy of fine features of the original turbulence-corrupted optical beam. In practice, the MPLC device is limited to a relatively small N due to implementation complexity, optics and fiber coupling losses, and cost. With only a limited number of HG modes for decomposition, smaller intensity features in the original beam are not well reproduced, yielding a net power loss, i.e., P R < P, where P and P R are the optical beam powers of the original and the N-mode HG reconstructed beams, respectively. As a result, conversion loss from free-space to optical fiber due to an insufficient number of modes can be significant. Due to the random nature of turbulence, such a conversion or power loss is expected to fluctuate in a similar way as the effect of scintillation. In addition, there are static losses from optics (e.g., Fresnel loss) in the MPLC device, but they are not affected by turbulence.
To better understand the impact of C 2 n (the refractive-index structure parameter) and N (the number of HG modes used for transformation of the original beam) on power loss or the fidelity of the reconstructed beam, wave-optics simulations were performed to generate many different turbulence-corrupted beam profiles. The simulated turbulence-corrupted beam is then decomposed into N HG modes, and the corresponding N complex HG mode coefficients c mn are computed, as described in Sec. 7.1 in the Appendix. Statistics of the power loss extracted from the c mn values versus different C 2 n and N were estimated with many turbulence propagation simulation trials. The case studied here followed parameters of a prototype FSOC system designed and constructed at the U.S. Army Research Laboratory (ARL). Simulation results suggest that, for C 2 n ≤ 10 −13 m −2∕3 , the power loss or penalty below 1 dB is achievable with only 10 HG modes decomposition. The overall variation of the power loss or normalized variance is below 0.01 over 1200 simulation trials. A larger number of HG modes or lower C 2 n yields a better power loss and stability as expected. The relatively low penalty and tolerance to high C 2 n is in part due to the small receiver aperture diameter (D ¼ 46 mm) relative to the diffraction-only beam diameter (∼1.2 m) used in ARL's FSOC prototype system, i.e., a reduced D∕r 0 ratio, where r 0 is the Fried parameter. 9 At the MPLC outputs, the N transformed Gaussian beams are coupled into N single-mode fibers, which are then coherently combined into one single-mode fiber or waveguide, as depicted in Fig. 1, for N ¼ 4. A silicon-on-insulator photonic integrated circuit (PIC) beam combiner, with its small footprint, is a candidate for achieving this function. For example, a compact low-loss Y-branch or Y-junction-based 10,11 2 m -to-1 tree combiner structure (m: integer) with integrated thermo-optic phase shifters used in some optical phased array designs 12 with an edge couplers array can be readily employed to accept single-mode fibers coupling. The efficiency of such a symmetric tree combiner is largest when all input powers are the same and all input optical phases are aligned, i.e., identical or modulo of 2π. However, as shown later, the mode power or HG mode coefficients distribution can be highly nonuniform due to turbulence-induced beam distortion, yielding uneven distribution of the N Gaussian beam powers. As a result, degraded beam power combining efficiency is expected. Analysis and simulations were conducted to quantify such a degradation or penalty versus C 2 n and N for 8-, 16-, and 32-to-1 Y-branch tree combiners.

HG Modes Decomposition
In this section, the definition of the power loss or ratio that quantifies the performance of the lowmode-number HG decomposition in MPLC is introduced. This is followed by a discussion on the spatial scale factor for HG mode required to maximize the power ratio or beam reconstruction fidelity.

Power Ratio Definition
A power ratio is defined here as the ratio of the optical power of the reconstructed (P R ) to the original (P) beam and is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 1 1 6 ; 3 2 9 Power ratio ≡ P R P ¼ RR jE R ðx; yÞj 2 dx dy RR jEðx; yÞj 2 dx dy ≤ 1; where E R ðx; yÞ and Eðx; yÞ are the reconstructed and original optical field spatial profiles after the receiver aperture, respectively. The reconstructed beam is a superposition of N HG modes with complex coefficients, c mn , computed from the original turbulence-corrupted beam at the receiver, Eðx; yÞ. It can be shown that this power ratio definition is equivalent to the conventional mode overlap integral, as described in Sec. 7.2 in the Appendix. The power ratio and its turbulence statistics investigated here provide one of several performance metrics for free-space coupling into single-mode fiber via MPLC.

Spatial Scale Factor of HG Modes
HG modes are characterized by the mode number or index ðm; nÞ and a spatial scale factor or W 0 (HG mathematical form and plots are shown in Sec. 7.3 in the Appendix). Because the spatial extend of HG mode is determined by W 0 , careful selection of its value is critical for accurately reconstructing the original beam and achieving a maximum power ratio. The selection of W 0 is not entirely obvious, especially at strong turbulences. Intuitively, W 0 should not be too different from the original beam size. However, as the product C 2 n L (L is the free-space propagation distance) increases, the Fried parameter or atmospheric coherence diameter seen by the original propagated beam reduces, leading to more likelihood of breakup of the optical beam into many smaller beam spots and patches with random sizes, i.e., the beam energy shifting to higher spatial frequencies with more finer intensity features. The beam size then becomes an increasingly ill-defined quantity. As a result, more smaller scale higher-order HG modes are required to reconstruct these smaller features, leading to a decrease of the "optimal" scale factor W opt 0 . One approach to select an optimal W 0 value is by searching for a W 0 that maximizes P R , i.e., the reconstructed beam power or the power ratio. This W opt 0 value, however, is expected to vary due to the random nature of turbulence even though C 2 n and L remain constant. Many propagation simulations are needed to extract the statistics of W opt 0 , and the ensemble average value of W opt 0 or hW opt 0 i is used to inform the design of MPLC devices and experimental testing.

Simulation Method and Results
Wave-optics simulations 13 were performed to generate a turbulence-corrupted 1550-nm optical beam profile at a receiver 2 km away. This was then used to compute N complex HG mode coefficients c mn , as described in Sec. 7.1 in the Appendix. These N HG modes were combined to reconstruct the original beam from which the power ratio was computed. To extract relevant statistics, 1200 independent simulation runs were executed for C 2 n values of 10 −14 , 10 −13 , and 10 −12 m −2∕3 and for N ¼ 10, 15, and 21. Using the Rytov variance 14 , the three C 2 n considered here are typically categorized as weak fluctuations (σ 2 1 ¼ 0.71 < 1), strong fluctuations (σ 2 1 > 1), and focusing/saturation regime (σ 2 1 ≫ 1), respectively. The wave-optics simulation was based on split-step propagation with the angular spectrum form of the Fresnel diffraction integral. A source was specified at the input plane and propagated across alternating instances of free space and a turbulent phase screen. The angular spectrum form of the Fresnel diffraction integral via 2D fast Fourier transform was used for the propagation calculations. Propagation was bounded by a super-Gaussian absorbing boundary to suppress energy spreading beyond the simulation grid similar to the concept of data windowing and antialiasing filter. 13 Five independent random phase screens separated by 400 m were placed along the 2-km propagation path. Accurate phase screen representation of the turbulent profile is critical for obtaining reliable scintillation data and statistics. These five random phase screens were generated using the 2D Fourier transform technique. 13 These phase screens were enhanced with subharmonic components to improve accuracy and low spatial frequency representation. 13 The modified von Kármán power spectral density (PSD) was used for the phase screen computation, taking into account both inner and outer scales and its ability to collapse to a Kolmogorov PSD. The inner and outer scales used in the simulations were 5 mm and 5 m, respectively. A Cartesian grid with 2048 × 2048 mesh points was used to represent the optical beam profiles with a 1-mm mesh size and an interpolated 29.2135-μm mesh size after beam truncation through the 46-mm receiver aperture. The quadratic phase-front envelop of the turbulence-corrupted beam due to diffraction was removed by subtracting the quadratic phasefront of a zero-turbulence propagated beam from it. Further details on wave-optics simulation can be found in Ref. 13 and the references therein.
The launched optical beam at a wavelength of 1.55 μm has a Gaussian profile that was used in all simulations. The Gaussian beam has a 1∕e 2 diameter 22.5 mm with a full beam divergence of 0.6 mrad. The diffraction-only beam diameter at 2 km is 2w 1 ¼ 1.2 m due to beam spreading without turbulence. The receiver aperture diameter is 46 mm. As a result, only a small fraction of beam area (½46∕1200 2 ¼ 0.15%) is collected by the receiver. These parameters are derived from a prototype FSOC terminal designed and constructed at ARL. Figure 2(a) shows a series of snapshots of simulated turbulence-corrupted beam (field amplitude) within a 2 m × 2 m viewing window at 2 km for C 2 n ¼ 10 −14 to 10 −12 m −2∕3 . Fried parameters computed for spherical wave (r 0 ¼ ½0.159k 2 C 2 n L −3∕5 ) along with D∕r 0 values are listed on the plots: D∕r 0 ¼ 0.49, 1.96, and 7.82 for C 2 n ¼ 10 −14 , 10 −13 , and 10 −12 m −2∕3 , respectively. D ¼ 46 mm is the receiver aperture diameter, and L ¼ 2 km is the propagation distance; these two quantities are fixed.
The beam after the 46-mm aperture is analyzed and hereafter is referred to as the original turbulence-corrupted beam to be used for HG modes decomposition and reconstruction.
Optical powers (P) after the 46-mm receiver aperture collected for 1200 simulation trials are shown in Fig. 2(b) (left to right), for C 2 n ¼ 10 −14 , 10 −13 , and 10 −12 m −2∕3 , respectively. The mean (μ) and normalized variance (σ 2 n ¼ σ 2 ∕μ 2 , σ: standard deviation) of P are shown in the legend. Note that σ 2 n increases by about a factor of 10 for C 2 n from 10 −14 to 10 −13 m −2∕3 but only by a factor of 2.4 from 10 −13 to 10 −12 m −2∕3 , indicating saturation. Figure 3 shows a series of snapshots (leftmost column) of the simulated turbulent-corrupted beam (optical field amplitude) truncated by the 46-mm diameter receiver aperture at 2 km for C 2 n values of 10 −14 , 10 −13 , and 10 −12 m −2∕3 within a viewing window slightly larger than the receiver aperture. The corresponding HG-based reconstructed beam profiles for N ¼ 10, 15, and 21 are shown in the second to fourth columns (left to right), respectively. As can be seen, the fidelity of the reconstructed beam decreases at higher C 2 n but generally improves with larger N. The W 0 values shown in the reconstructed beam plots are W opt 0 in the number of grid points, and they tend to decrease as C 2 n increases. The original beam power P defined earlier is referenced to the truncated beam power or total power over the 46-mm diameter receiver aperture.
Note that, for C 2 n ¼ 10 −14 m −2∕3 shown in the top-left plot, the beam profile is relatively flat over the receiver aperture. Decomposition of a flat profile with the first few HG modes is challenging, as can be seen in the plot next to the top-left one in which a hole or dip near the center appears to "balance" the two adjacent peaks to achieve an effective flat profile. Figure 4 shows distribution of N HG mode coefficient amplitudes, jc mn j, corresponds to reconstructed beams for C 2 n and N values depicted in Fig. 3. As expected, jc 00 j for the fundamental HG or Gaussian mode is largest at low C 2 n , but it spreads to other jc mn j with higher turbulence strengths. The general trend of jc mn j versus the mode number appears exponential-like, with more complex substructures at higher C 2 n values. Note that HG mode coefficients generally decrease with increasing C 2 n values. This is due to the reduction of optical power transmitted through the receiver aperture (46-mm diameter) as the beam distortion increases with the turbulence strength.
The power ratio or conversion loss for 1200 independent simulation trials is shown in Fig. 5 for different C 2 n and N values. The power ratio shown was obtained using ensemble averaged hW opt 0 i generated from separate sets of 1200 simulation runs, where W opt 0 was determined via power ratio maximization in each run (see Sec. 7.3 in the Appendix). It is clear that the power   Leftmost column: simulated turbulence-corrupted beam (optical field amplitude) truncated by the 46-mm diameter receiver aperture at 2 km for C 2 n ¼ 10 −14 , 10 −13 , and 10 −12 m −2∕3 from top to bottom rows, respectively. The reconstructed beam from the same turbulence-corrupted beam on each row using 10, 15, and 21 HG modes is shown in the second to fourth columns, respectively. W 0 is the spatial factor for the selected HG mode that yields the maximum power ratio.
ratio and its fluctuation improve with larger N and lower C 2 n . The mean (μ) and normalized variance (σ 2 n ) values of the power ratio over 1200 samples are also shown in the plots. Figure 6 shows σ 2 n of jc mn j versus the HG mode number. Note that σ 2 n for jc 00 j (HG mode 1) increases with C 2 n , indicating that the contribution of jc 00 j becomes less reliable at high turbulence conditions, as expected. Figure 7 shows the convergence of σ 2 n for the power ratio versus the accumulated number of simulation runs up to 1200. As can be seen, the statistics of power ratio generally stabilized around 600 or more trials. This suggests that our 1200-run statistical results should be reliable for the current investigated case.  Statistics of the power ratio, mean (μ), and normalized variance (σ 2 n ¼ ðσ∕μÞ 2 ) for different C 2 n and N over 1200 samples are shown in Fig. 8, including the case of N ¼ 36. Numbers highlighted in red are ensemble averaged power ratios in decibel units, which can be used as a guide in the link budget computation. The last column shows the ensemble average of the ratio hW opt 0 i∕R, where R ¼ 23 mm is the receiver aperture radius. Bar graphs for μ in dB units and σ 2 n versus N are shown on the right panel for a visual comparison. For the HG mode decomposition penalty to be <1 dB, MPLC with a mode number as low as 10 can tolerate a C 2 n up to about 10 −13 m −2∕3 . The power ratio fluctuation or σ 2 n is below 0.01 in this case, indicating stability and robustness of the conversion to turbulence with C 2 n ≤ 10 −13 m −2∕3 , as can be seen in the table and bar graph.

Power Combiner
The N Gaussian outputs of MPLC devices are coupled into N single-mode fibers, which are then coherently combined into one single-mode fiber or waveguide. A silicon-on-insulator PIC Y-branch 10,11 tree combiner with 2 m inputs and phase shifters for input phase alignment is considered and analyzed. Then, simulation results are discussed. Other PIC power splitting/ combining structures, such as directional couplers, can also be considered. The Y-branch was selected here due its compactness, low excess loss, relatively broad and flat spectral response, and fabrication tolerance based on our experience of fabricated PIC Y-branches gained from a series of silicon photonics multiproject wafer runs over the years. The performance of the power combiner or its efficiency for different C 2 n and N values was computed using the complex HG coefficients values or c mn obtained from simulations. For clarity and to elucidate the impact of turbulence, internal mode conversion losses of the MPLC and fiber coupling and other static losses are ignored. These losses are turbulence-independent and can be readily included in the combiner efficiency once they are available or measured. With these assumptions and the normalized HG modes, the input optical fields to the power combiner are equal to the HG mode coefficients or c mn , as shown in Sec. 7.3 in the Appendix. Figure 9 shows an example for a schematic of a PIC 4-to-1 Y-branch symmetric tree combiner with main components, such as edge couplers, phase shifters, and Y-branches, 10,11 to combine outputs of MPLC. Here again, N ¼ 4 was selected for clarity and illustrative purpose only. The PIC phase shifters are used to perform phase-alignment: equalizing or aligning optical phases of the input optical signals from the MPLC output fibers to maximize and stabilize the combined output optical power. These phase shifters can be based on thermo-optics using resistive heaters or plasma dispersion in the PN (P-type, N-type) junction or PIN (P-type, Intrinsic, N-type) waveguide via doping. 15 Recent reported phase shifters 16,17 based on optomechanical PIC electrostatic actuators are also candidates to meet low SWaP requirements. Electrical controls for the phase shifters are not shown for clarity. A feedback control loop (e.g., Ref. 18) that maximizes the combiner output can be used to drive the phase shifters array, which is not considered here.
As described earlier and in Sec. 7.4 in the Appendix, the Y-branch coherent combiner output is maximum when all of its inputs have the same power and same phase or are phase aligned. As shown in Fig. 4, jc mn j or mode power distribution can be highly nonuniform, yielding an uneven distribution of the N Gaussian beam powers to the combiner inputs. Also, the 2 m (m: integer) number of inputs of the Y-branch combiner does not always match N, yielding suboptimal efficiency. The impact of these impairments, along with derivation of combiner efficiency, which is defined as the ratio of the combined output power (P o ) to the original receiver aperture truncated beam power (P), is described in Sec. 7.4 in the Appendix. Figure 10 shows the simulated 8-to-1 Y-branch tree combiner efficiency without (η) and with (η ϕ ) phase alignment computed for 1200 simulation trials for different C 2 n (top to bottom rows) Fig. 9 Schematic of a 4-to-1 PIC power combiner with Y -branches that accepts outputs of MPLC and produces a combined output power. Electrical controls to the phase shifters are not shown for clarity.
and N (left to right columns). The legend shows the normalized standard deviation (σ 2 n ) and mean (μ) of η and η ϕ over 1200 samples. The numbers in red are mean efficiencies in decibels.
Efficiencies for the 16-and 32-to-1 combiners were also computed. Figure 11 shows stacked bar graphs for statistics of efficiencies, mean μ, and normalized variance σ 2 n of η (left column) and η ϕ (right column) for 8-, 16-, and 32-to-1 combiners for the investigated C 2 n and N values.  Each stacked bar graph is partitioned in three parts from left to right for the 8-, 16-, and 32-to-1 Y-branch tree combiners, respectively. As can be seen, an 8-to-1 combiner yields best efficiency for 10 HG modes, whereas a 16-to-1 combiner should be considered for 15 and 21 HG modes. Overall mean combiner efficiencies vary between −5.69 and −13.34 dB without phase alignment and between −2.36 to −5.22 dB with phase alignment for the three C 2 n and three N values. The Y-branch combiner efficiency is almost always suboptimal, as shown in Sec. 7.4 in the Appendix, due to nonequal input powers from the MPLC outputs because of turbulence-induced beam distortion. Furthermore, scalability of the combiner is likely limited due to accumulated PIC component losses at high input counts.

Discussions
The investigations conducted here mainly focused on impact of turbulence on MPLC and power combiner performance. For clarify, turbulence-independent static losses due to MPLC optics, fiber coupling, and PIC component are not included. Other imperfections of MPLC, such as internal mode conversion loss and mode crosstalk in MPLC, are not considered to elucidate the turbulence-induced penalty.
The combiner efficiency described earlier is defined and normalized to the original beam power. Thus, the power ratio of the reconstructed beam is already included in the combiner efficiency. As a result, the combiner efficiency can be interpreted as a coupling efficiency from the free-space optical beam after the receiver aperture to a single-mode fiber or waveguide.
The power combiner analysis presented here assumed an ideal Y-branch with no excess loss and perfect 50:50 power split. In practice, both excess loss and splitting ratio deviation from ideal are nonzero and wavelength dependent. For the Y-branch combiner, the accumulated excess loss increases geometrically as m due to the cascading of Y-branches.
The feedback control of phase shifters to maximize the combiner output will require further studies, which is beyond the scope of this investigation. An efficient optimization algorithm, such as stochastic parallel gradient descent, 18 can be considered.
The fiber length difference of the connecting fiber array between MPLC and PIC generates signal distortion at the combiner output similar to intersymbol interference (ISI) in multimode fibers. For a fiber length accuracy within 5 mm or a maximum delay of about 25 ps, the maximum symbol rate would be 20 Gbaud for a maximum of a half-symbol-period delay. Digital signal processing, such as adaptive equalization, 19 can be applied to compensate for fiber differential delay-induced ISI.
The overall receiver performance, such as signal-to-noise ratio (SNR), using a Y-branch combiner were estimated and are described in Sec. 7.5 in the Appendix. To mitigate the impact of polarization variation in single-mode fibers, polarization-maintaining fibers or polarization diversity using a PIC polarization rotator-splitter 20 can be employed. The output of the combiner shown in Fig. 9 can be connected directly to a PIC photodetector integrated on the same chip, 21 e.g., leveraging mature passive and active PIC fabrication capabilities of photonics foundries, such as American Institute for Manufacturing Integrated Photonics. 22

Summary
Analysis and a series of simulations were conducted to quantify the performance of MPLC and HG decomposition/reconstruction of turbulence-corrupted free-space optical beams for different C 2 n and N values. As a case study, simulations were conducted for a specific launched beam size, divergence, and propagation distance for a prototype FSOC system that was designed and constructed at ARL. These results provide an initial assessment on potential performance of the MPLC device with the prototype system. For example, depending on the requirements of overall loss, a design parameter space of C 2 n and N that is subject to constraints of the power loss based on simulation results can be formulated. The method and approach described here can be readily applied to other FSOC laser transceivers.
Performance analyses were also conducted for a combiner that serializes multiple parallel outputs of MPLC into a single optical signal. Efficiencies of a "traditional" Y-branch symmetric tree combiners (8-, 16-, and 32-to-1) were computed for different C 2 n and N values. The overall robust performance of the passive MPLC and PIC combiner makes it an attractive alternative to adaptive optics system for turbulence-tolerant FSOC applications.

Power Ratio Definition
The turbulence-corrupted beam optical field profile at the receiver, Eðx; yÞ, is expressed in terms of an infinite sum of orthogonal HG modes, HG mn ðx; yÞ, as follows: c mn HG mn ðx; yÞ 2 dx dy: The above inequality is due to triangular inequality for any two nonzero complex numbers a and b: ja þ bj 2 ≤ jaj 2 þ jbj 2 . Recognizing that the double-integral above is equivalent to the total optical beam power, then E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 1 ; 1 1 6 ; 1 7 0 P ≤ P R þ Δ; where P ¼ RR jEðx; yÞj 2 dx dy is the original beam power, P R ¼ RR jE R ðx; yÞj 2 dx dy is the reconstructed beam power, and Δ ¼ RR j P ∞ Nþ1 c mn HG mn ðx; yÞj 2 dx dy ≥ 0 is a residual power loss term. Because Δ is a nonnegative number, P − P R ≤ Δ ⇒ P − P R ≥ 0 ⇒ P R ≤ P; Δ ≥ 0: The above analyzed equation shows that optical power of the reconstructed beam using a finite number of HG modes is always less than the optical power of the original beam. The power ratio quantity is defined as the ratio of the two optical powers as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 1 ; 1 1 6 ; 6 7 6 Power ratio ≡ P R P ≤ 1: This quantity can be readily computed via wave-optics simulation for various turbulence scenarios and N values. It can be shown that the power ratio quantity above is identical to the conventional definition of the mode overlap integral commonly used to quantify coupling efficiency between two different optical modes, as described next.

Mode Overlap Integral
The mode overlap integral is defined as follows: 23 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 2 ; 1 1 6 ; 5 3 7 where P and P R are optical powers of the original and reconstructed beams, respectively. The integral in the numerator using HG mode decomposition is given by The above is the same as P R because E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 2 ; 1 1 6 ; 2 6 5 where HG mn ðx; yÞ is normalized as described below, yielding RR jHG mn ðx; yÞj 2 dx dy ¼ 1. This results in the following E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 2 ; 1 1 6 ; 1 9 0 yielding a power ratio the same as η.

HG Modes
The HG mode 24  where δ ij is the Kronecker delta. For fundamental or Gaussian mode, w 0 is simply the 1∕e 2 radius of the intensity profile because jHG 00 ðx; yÞj 2 ¼ A 2 exp½−2ð x 2 þy 2 w 2 0 Þ. Figure 12 shows jHG mn ðx; yÞj for the first 36 modes. Note that the numbering of the mode in relation to the ðm; nÞ index shown in the figures is adopted here. The w 0 is same as the spatial scale factor W 0 described in the main text. Optimal value of W 0 or W opt 0 was obtained by maximizing the power ratio. Figure 13 shows the ratio W opt 0 ∕R, where R ¼ 23 mm is the receiver aperture radius obtained for each 1200 simulation run for different C 2 n and N. As can be seen, W opt 0 tends to reduce with higher N and C 2 n values. The legend shows the ratio of the ensemble averaged W opt 0 to R or hW opt 0 i∕R. With the normalized HG modes defined above, the MPLC outputs to the combiner are derived as follows. The N HG modes are transformed to a series of spatially separated Gaussian modes at the MPLC output: where MfHG mn ðx; yÞg ¼ q mn G mn ðx; yÞ and G mn ðx; yÞ is an mn-dependent spatially shifted Gaussian mode and Mf·g is a transform (free-space propagation through a series of phase plates) that converts HG mn ðx; yÞ to G mn ðx; yÞ. The complex coefficient q mn is associated with the transform from HG mn ðx; yÞ to G mn ðx; yÞ via Mf·g; it may also include any losses from the optics. Because q mn is a static loss mostly from the MPLC implementation loss and is independent of turbulence, it was assumed here that all q mn ¼ 1 and G mn ðx; yÞ is normalized, i.e., RR G mn ðx; yÞdx dy ¼ 1. Furthermore, assuming no coupling loss from the MPLC outputs to the fiber arrays or fiber arrays to the PIC chip, the j'th input signal optical field (E Sj ) to the power combiner is simply given by the HG modes coefficients: E Sj ¼ c mn .

Power Combiner Analysis
The performance of the power combiner is analyzed, and a combiner efficiency is defined here for the Y-branch tree combiner. Consider a single ideal lossless Y-branch with two input optical fields, E 1 ¼ jE 1 je iϕ 1 and E 2 ¼ jE 2 je iϕ 2 , as shown in Fig. 14. The combined output optical field is given as 15 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 4 ; 1 1 6 ; 2 1 8 The combined output power P o is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 4 ; 1 1 6 ; 1 6 1 where P 1 ¼ jE 1 j 2 and P 2 ¼ jE 2 j 2 are the input optical powers. Consider the case in which the input phases are aligned such that ϕ 1 ¼ ϕ 2 , i.e., constructive interference in the fundamental mode at the output waveguide. The combiner efficiency is defined as follows: and N. Note that W opt 0 tends to reduce with higher N and C 2 n values. R ¼ 23 mm is the receiver aperture radius.
Note that, if P 2 ¼ P 1 , the combined output P o ¼ 2P 1 and η ¼ 1. Figure 14 shows a plot of η versus the ratio P 2 ∕P 1 . The efficiency reaches maximum at 1, where the output power is simply the sum of the two input powers only when both input powers are equal. When one of the input power is zero, the combined output power is half of the nonzero input power.
The above analysis of a single Y-branch combiner can be extended to 2 m -to-1 combiner. For N MPLC outputs to be combined using a 2 m -to-1 Y-branch combiner, the combined output optical field is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 4 ; 1 1 6 ; 6 0 8 where E j ¼ jE j je iϕ j is the complex optical field to the j'th input of the Y-branch combiner. The expression minðN; 2 m Þ means selecting whichever number is smaller in the parenthesis. For the phase-aligned case (i.e., all ϕ j are identical) using the compensating phase shifters, the combined output field is simply the sum of the input field amplitudes given as Due to the 1∕ð ffiffi ffi 2 p Þ m factor in the Y-branch combiner, unity combiner efficiency can only be achieved when all input amplitudes are identical and all phases are aligned.

SNR Analysis
Analysis of SNR for the Y-branch power combiner with direct detection is described.
The SNR of a Y-branch combiner with direct-detection using a single photodetector is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 5 ; 1 1 6 ; 3 4 5 where R is the detector responsivity in A/W, I d is the detector dark current in A, P YB is the combined optical power at the Y-branch combiner output in watts, B is the noise bandwidth in Hz, i th 2 is the detector thermal noise PSD in A 2 ∕Hz, and q is the electronic charge. The first, Fig. 14 (a) Schematic of a Y -branch with two input optical fields and a combined output optical field. (b) Combiner efficiency η versus input power ratio P 2 ∕P 1 for the case of aligned phase or ϕ 1 ¼ ϕ 2 ¼ 0. Note that the efficiency reaches maximum at 1, where the output power is simply the sum of the two input powers only when both input powers are equal.
second, and third term in the denominator represents shot, dark current, and thermal noise power of the detector, respectively. P YB is related to P rec or the optical power collected at the receiver aperture via the efficiency of the MPLC and Y-branch combiner described in the main text. For the ideal case and no turbulence, P YB ¼ P rec , the SNR is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 7 . 5 ; 1 1 6 ; 4 2 3 SNR ¼ ðRP rec Þ 2 2qRP rec B þ 2qI d B þ i th 2 B : Figure 15 shows the computed SNR versus P rec with Y-branch for thermal noise current i th ¼ 0.7 and 7 pA∕ p Hz. The other parameters are R ¼ 0.7 A∕W at 1550 nm, 21 I d ¼ 1.16 nA, 21 N ¼ 15, and B ¼ 1 GHz. Fig. 15 Computed SNR in dB versus received signal optical power in dBm at the receiver or P rec using Y -branch combiner for detector thermal noise current density 0.69465 and 6.9949 pA∕ p Hz.