## Appendix A Simulation of Atmospheric Turbulence  
We want to confirm that the polarization of light propagating through atmospheric turbulence is well-preserved. Recall from Section 1 that atmospheric turbulence causes fluctuations in the index of refraction \(n(\vec{r})\) at position \(\vec{r}\) that can be modeled by Equation (1). Throughout this appendix, all vectors are three-dimensional. These fluctuations introduce a depolarizing term into the wave equation describing light in a turbulent region as shown in Equation (2). We choose to use the COMSOL Ray Optics software [37] to simulate light propagating through turbulence because it numerically solves [PERSON]'s equations and therefore captures the effects of the depolarization term in Equation (2). Here we consider three values for the initial DOP and perform 1000 simulations for each DOP. In each simulation, the initial SOP was set to the arbitrary state used in Figure 4. The fluctuations of the refractive index are generated assuming the Von Karman spectrum [18]. \[\phi_{n}(\vec{x})\propto\frac{0.033C_{n}^{2}}{(\vec{x}^{2}+\kappa_{0}^{2})^{11/6}}e ^{-\frac{\vec{x}^{2}_{n}}{\kappa_{0}^{2}}}, \tag{10}\]  
where \(\vec{x}\) is the spatial frequency, \(C_{n}^{2}\) is the structure constant, and \(\kappa_{0}\) and \(\kappa_{m}\) are the spatial frequencies corresponding to the inner and outer scales of the turbulence respectively. The kernel \(G(\vec{r})\) that defines the statistics of the fluctuations in position space is then given by [45; 46]. \[G(\vec{r})=\mathcal{F}_{\vec{x}\rightarrow\vec{r}}^{-1}[\sqrt{\phi_{n}(\vec{x })}], \tag{11}\]  
where \(\mathcal{F}_{\vec{x}\rightarrow\vec{r}}^{-1}\) denotes the inverse Fourier transform. By convolving \(G(\vec{r})\) with random uniform noise, fluctuations in \(n(\vec{r})\) with the desired statistics can be generated [45; 46]. If the vectors in Equations (10) and (11) are restricted to two dimensions, the approach for generating the thin phase screens for the aptly named phase screen method is recovered [47; 48; 49; 50; 51]. Here, we use three-dimensional vectors to generate a volume with the appropriate fluctuations in \(n(\vec{r})\). As before, the code for generating these fluctuations and performing the simulations can be found on Github [38]. The results are shown in Figure 10 and demonstrate that the DOP and SOP are preserved as expected from previous work [19]. In this appendix, we examine the experimentally observed and simulated intensity distributions. Figure A2 displays histograms of the intensities obtained from the simulation and experiment for light that is initially vertically polarized with DOP = 0.8. The experimental intensity distribution here is from the measurement of \(I(0^{\circ},90^{\circ})\). The simulation results are cutoff above \(I=2\) and below \(I=0.03\) to better reflect the dynamic range of the detector while maintaining the structure of the distribution. These histograms demonstrate that our simulation yields intensity fluctuations similar to those seen in the experiment. We also present the probability density against the normalized intensity for the simulated and experimental results in Figure A3. For both sets of results, we fit the data to the gamma and log-normal distributions' probability density functions and plot them over the histogram. These distributions, commonly used for modeling turbulence due to bubbles in water [29], are characterized by shape \(s\), location \(l\), and scale \(c\) parameters, while the quality of the fit is given by the coefficient of determination \(R^{2}\). For the simulated data, we obtain \(s=0.805\), \(l=0.005\), \(c=0.145\), and \(R^{2}=0.350\) for the gamma distribution and \(s=1.74\), \(l=0.005\), \(c=0.064\), and \(R^{2}=0.331\) for the log-normal distribution. For the experimental data, we obtain \(s=1.35\), \(l=0.028\), \(c=0.171\), and \(R^{2}=0.584\) for the gamma distribution and \(s=0.794\), \(l=0.001\), \(c=0.194\) and \(R^{2}=0.455\) for the log-normal distribution.