\section{Introduction}
\label{sec:intro}
Significant theoretical and experimental effort is dedicated to mapping out the QCD phase diagram in the temperature $T$ and baryon chemical potential $\mu_B$ plane, and to search for the elusive critical point \cite{Aoki:2006we,Aoki:2006br,Borsanyi:2010bp,Bazavov:2011nk,Borsanyi:2013bia,Bazavov:2014pvz,Fodor:2018wul,Busza:2018rrf,Stephanov:2011pb,Stephanov:1998dy,Ratti:2006gh,Fukushima:2010bq,Critelli:2017oub,Luo:2017faz,Parotto:2018pwx,Grefa:2021qvt,Karthein:2021nxe,Adamczyk:2017iwn,Adamczewski-Musch:2020slf,Abelev:2013vea,Fu:2021oaw} (for recent reviews see e.g. \cite{Ding:2015ona,Ratti:2018ksb,Bzdak:2019pkr,Ratti:2021ubw}). Ultrarelativistic heavy-ion collisions are generating the deconfined phase of strongly interacting matter in the laboratory. By systematically decreasing the collision energy of the incoming nuclei, the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory is scanning the phase diagram in the so-called Second Beam Energy Scan (BESII), soon to be followed by even lower collision energies at NICA and at the GSI-FAIR accelerator.

First-principles lattice QCD simulations are available for several thermodynamic quantities, such as the equation of state at zero \cite{Borsanyi:2010cj,Borsanyi:2013bia,Bazavov:2014pvz} and small chemical potential \cite{Guenther:2017hnx,Gunther:2017sxn,Bazavov:2017dus,Borsanyi:2021sxv,Mondal:2021jxk}, QCD transition line \cite{Bellwied:2015rza,Bazavov:2018mes,Borsanyi:2020fev,Borsanyi:2021sxv}, as well as diagonal \cite{Borsanyi:2011sw,Borsanyi:2013hza,Bellwied:2013cta,Borsanyi:2014ewa,Bellwied:2015lba,Noronha-Hostler:2016rpd,Bazavov:2017tot,Borsanyi:2018grb,Bazavov:2020bjn} and off-diagonal \cite{Borsanyi:2018grb} fluctuations of conserved charges. However, they are currently limited to the low baryonic chemical potential regime, due to the fermionic sign problem. Effective models that can reproduce lattice QCD results in certain regimes of temperature and chemical potential are therefore very useful to extend the coverage of the phase diagram beyond the reach of lattice QCD.
The models are also necessary to make a connection to the common heavy-ion observables, such as the measurements of various distributions of identified hadrons.

In the low-temperature regime ($T \lesssim 160$), the bulk thermodynamics of QCD is generally well described by a multi-component gas of free hadrons and resonances \cite{Bazavov:2012jq,Borsanyi:2010bp,Vovchenko:2016rkn,Alba:2017mqu} the so-called Hadron Resonance Gas (HRG) model. This indicates that hadron interactions in this regime may be dominated by the formation of known resonances. The HRG model has been widely used to study the confined phase of QCD matter below the transition line \cite{Karsch:2003vd,Karsch:2003zq,Tawfik:2004sw,Huovinen:2009yb,Ratti:2010kj,Alba:2015iva,Huovinen:2017ogf,Bellwied:2019pxh,Vovchenko:2014pka}. Its remarkable agreement with the equation of state from first-principles lattice calculations has led to its popularity, especially in the study of the chemical freeze-out in HICs \cite{Becattini:2012xb,Alba:2014eba,Vovchenko:2015idt,Andronic:2017pug,Bellwied:2018tkc,Alba:2020jir,Bluhm:2018aei,Flor:2020fdw}.
However, with the availability of more differential observables like susceptibilities of conserved charges, discrepancies between the predictions of the fundamental theory and the HRG model have been observed \cite{Bazavov:2014xya,Borsanyi:2018grb,Bellwied:2021nrt}, specifically at temperatures $T \sim 150-160$~MeV that characterize the chemical freeze-out in heavy-ion collisions at the highest energies. In particular, lattice QCD results for the partial pressures have shown a need for more resonances in the strange sector than those which are already experimentally well known \cite{Bazavov:2014xya,Alba:2017mqu}. On the other hand, some susceptibility ratios including $\chi_4^B/\chi_2^B$ and $\chi_4^S/\chi_2^S$ suggest the need for repulsive interactions \cite{Bellwied:2015lba,Bazavov:2018mes,Alba:2017mqu}.

Several extensions of the HRG model have thus been proposed to improve the agreement with lattice QCD.
One possibility is extending the hadronic spectrum -- the model input -- to include more states not yet observed \cite{Bazavov:2014xya,Alba:2017mqu,Alba:2020jir}.
Other extensions incorporate additional, non-resonant interactions between hadrons such as excluded volume~\cite{Yen:1997rv,Andronic:2012ut,NoronhaHostler:2012ug,Bhattacharyya:2013oya,Vovchenko:2014pka,Albright:2015uua,Satarov:2016peb,Vovchenko:2017xad,Alba:2017bbr,Vovchenko:2018eod,Motornenko:2020yme}, van der Waals~\cite{Vovchenko:2016rkn,Vovchenko:2017zpj,Samanta:2017yhh,Sarkar:2018mbk,Vovchenko:2020lju}, mean field~\cite{Huovinen:2017ogf,Steinert:2018zni}, or are based on scattering phase shifts~\cite{Venugopalan:1992hy,Friman:2015zua,Vovchenko:2017drx,Lo:2017lym,Dash:2018can,Fernandez-Ramirez:2018vzu,Dash:2018mep}.
While one of the advantages of the standard HRG model is certainly the lack of free parameters, apart from the uncertainties in the hadronic spectrum, introducing additional interactions unavoidably leads to new free parameters that need to be constrained through comparison with lattice results. We propose a combination of these different corrections to the standard HRG model and investigate its agreement with several lattice results on the susceptibilities.

In this manuscript, we consider two HRG model extensions: the excluded volume interaction in the baryon sector, and the inclusion of additional particles in the hadronic list, beyond those experimentally observed. We emphasize the complementary nature of these two extensions and identify combinations of susceptibilities of conserved charges that allow to constrain them separately. The resulting HRG model considerably improves the description of lattice QCD results, in particular those for fluctuations of conserved charges.
The model can thus be useful for the analysis of freeze-out in heavy-ion collisions, in particular event-by-event fluctuations.

\section{The HRG Model and its modifications}\label{sec:HRG_ext}
\subsection{Ideal HRG}
The partial pressure for particle species $i$ in the ideal HRG model can be written as:
\begin{equation} \label{eq:press_partitionfxn}
\centering
\begin{split}
    &P_i (T, \mu_B, \mu_Q, \mu_S) = \frac{d_i T}{2\pi^2} \int_0^{\infty} (-1)^{B_i+1} k^2 dk \, \times
    \\ &\times \ln {\left[1 + (-1)^{B_i+1} \lambda_i(T, \mu_i) \exp\left(-\frac{\sqrt{k^2 + m_i^2}}{T}\right)\right]},
\end{split}
\end{equation}
where $d_i$ is the spin degeneracy factor, $B_i$ is the baryon number of species $i$, $k$ is the momentum, $\lambda_i (T,\mu_i) = \exp[(B_i \mu_B + Q_i \mu_Q + S_i \mu_S)/T]$ is the particle fugacity, and $m_i$ is the mass of species $i$.
The partial pressure can be presented as a series containing the modified Bessel function of the second kind by expanding the logarithm and integrating term-by-term:
\begin{equation} \label{eq:press_Bessfxn}
\centering
\begin{split}
    &P_i (T, \mu_B, \mu_Q, \mu_S) = \\
    &\frac{d_i T^2}{2\pi^2} \sum_{N=1}^{\infty} \frac{[(-1)^{B_i+1}]^{N+1} \lambda_i(T, \mu_i)}{N^2} m_i^2 K_2\left(N\frac{m_i}{T}\right).
\end{split}
\end{equation}
Taking the first term in the expansion corresponds to the Boltzmann approximation.
This approximation is sufficient for our purposes and will be used throughout.
The pressure attains the following convenient form:
\begin{equation} \label{eq:press_phi}
\centering
\begin{split}
    P_i (T, \mu_B, \mu_Q, \mu_S) &= d_i \tilde{\phi}(T,m_i) \lambda_i(T,\mu_i), \\
    \tilde{\phi}(T,m_i) &=  \frac{m_i^2 T^2}{2\pi^2} \, K_2(m_i/T).
\end{split}
\end{equation}
The full pressure in the ideal HRG model corresponds to the sum of partial pressures of all hadronic components.
It is convenient to group the contributions of the various hadrons in accordance with their quantum numbers.
For instance, introducing
\begin{equation} \label{eq:press_sectors}
\centering
    \tilde{\phi}(T) = \sum_{j \in \rm{sectors}} d_j \tilde{\phi}(T,m_j)
\end{equation}
allows us to identify the different sectors of the total HRG pressure broken down by various quantum numbers that are listed in Table \ref{tab:hadron_table}.
Note that all species apart from the $i = 0$ sector~(neutral particles)
have a corresponding antiparticle sector.
\begin{table}[h]
\centering
\begin{tabular}{c|c c c|c}
 $i$ & $B$ & $Q$ & $S$ & base hadron \\ \hline
 0 &  0 & 0  & 0 & $\pi^0$ \\
 1 &  0 & 1 & 0 & $\pi^+$ \\
 2 &  0 & 1 & 1 & K$^+$ \\
 3 &  0 & 0 & 1 & K$^0$ \\
 4 &  1 & 0 & 0 & n  \\
 5 &  1 & 1 & 0 & p  \\
 6 &  1 & 2 & 0 & $\Delta^{++}$  \\
 7 &  1 & -1 & 0 & $\Delta^{-}$ \\
 8 &  1 & 0 & -1 & $\Lambda$ \\
 9 &  1 & 1 & -1 & $\Sigma^+$ \\
 10 & 1 & -1 & -1 & $\Sigma^-$ \\
 11 & 1 & 0 & -2  & $\Xi^0$ \\
 12 & 1 &  -1 &  -2 & $\Xi^-$ \\
 13 & 1 & -1 & -3  & $\Omega^-$
\end{tabular}
\caption{
\label{tab:hadron_table}
The list of hadronic quantum number sets contributing to the pressure of the ideal HRG model.
The last column identifies the lowest mass hadron representing the set of quantum numbers.
}
\end{table}
Therefore, the pressure from  Eq.~(\ref{eq:press_phi}) takes the form of a truncated relativistic expansion in fugacities:
\begin{equation} \label{eq:Pfug}
    \centering
    P (T, \mu_B, \mu_Q, \mu_S) = \tilde{\phi_0}(T) + \sum_{i \neq 0} 2 \, \tilde{\phi_i}(T) \cosh \left(\mu_i/T \right),
\end{equation}
where $\mu_i=B_i \mu_B + Q_i \mu_Q + S_i \mu_S$ is the chemical potential of the corresponding $i$th sector.
Each term in Eq.~\eqref{eq:Pfug} corresponds to the partial pressure associated with the particular set of hadronic quantum numbers.

\subsection{Extended Spectrum}
As a correction to the standard HRG model, our first extension is the incorporation of hadronic states beyond those which are well-known experimentally. When considering the partial pressures from Lattice QCD, it has been shown that the hadronic spectrum in QCD goes beyond what exists in the ordinary version of the hadronic list from the Particle Data Group (PDG) with only the most well-known states \cite{Zyla:2020zbs}.
The PDG ranks particles by how well established they are, with a rating based on a number of stars (*).
The **** states are those which are very well established, like e.g. nucleons or $\Delta(1232)$ resonances.
On the opposite side, the * states are the least established ones, for example $\Delta$(1750) and other high-mass resonances.
A previous investigation of the agreement between the partial pressures from the lattice and those obtained in the HRG model with different hadronic lists across the various quantum number sectors, found that the PDG2016+ particle list reproduced the largest number of quantities calculated on the lattice without exceeding them~\cite{Alba:2017mqu}.
This hadronic list incorporates all states from the 2016 version of the PDG~\cite{Patrignani:2016xqp}, starting from the well-known ones, all the way down to those listed as “seen”, with a confidence rating of *.
Another possibility is to incorporate the states predicted by the Quark Model~\cite{Capstick:1986bm,Godfrey:1985xj,Ebert:2009ub}, which includes an even larger number of hadrons than those contained in the PDG2016+.
As we can see in Eq. (\ref{eq:press_phi}), any additional states seek to increase the pressure of the system, i.e. a larger number of states will lead to a larger overall pressure.

It was recently pointed out in Ref. \cite{Bollweg:2021vqf} that some states calculated within the Quark Model from Ref. \cite{Capstick:1986bm} are overlapping with states later measured by the PDG \cite{Zyla:2020zbs}.
In light of this, we provide an update to the QM list first published in Ref. \cite{Alba:2020jir} that removes all the duplicate states that remained as an artifact of those early Quark Model calculations \cite{listrepo}.
We note that the calculations here denoted as QM are using this updated list and that we observed a minimal difference between the old and new list when calculating the fluctuations of interest in this study.

In this study, we aim at identifying the most suitable description of the lattice data, with both the excluded volume interactions and additional hadronic states incorporated.
In order to do so, we revisit the comparison with lattice data and investigate several different hadronic lists:
\begin{itemize}
    \item PDG2016 -- ordinary hadronic list with only the well-known states *** $-$ **** from the 2016 Particle Data Booklet;
    \item PDG2016+ -- the list containing both the established~(***$-$****) and unconfirmed~(*$-$**) states;
    \item Quark Model~(QM) -- the list which incorporates all states predicted by the Quark Model.
\end{itemize}
The latter two lists were introduced and described in detail in \cite{Alba:2017mqu, Alba:2020jir}.
We also checked that the most recent compilation of the established states from the Particle Data Group -- the PDG2020 list -- yields negligible differences compared to the PDG2016 list, thus we retain the latter list in the analysis for consistency with Refs.~\cite{Alba:2017mqu,Alba:2020jir}.
\subsection{Excluded Volume}

The next extension to the HRG model is the excluded volume model. This corresponds to including repulsive interactions between hadrons. Many versions of the EV-HRG model have been considered in the literature. Here we follow the approach introduced in Refs.~\cite{Vovchenko:2016rkn,Vovchenko:2017xad} where EV interactions are included only for baryon-baryon and antibaryon-antibaryon pairs.
This corresponds to a minimalistic EV extension that does not affect meson-meson and meson-baryon interactions, which are presumed to be dominated by resonance formation and thus already included in the HRG model.
The pressure is partitioned into contributions of non-interacting mesons and interacting baryons and antibaryons:
\begin{equation} \label{eq:pevhrg}
\centering
p = p_M^{\rm id} + p_B^{\rm ev} + p_{\bar{B}}^{\rm ev},
\end{equation}
where
\begin{equation} \label{eq:pMev}
p_{M}^{\rm id}  = \tilde{\phi}_0(T) ~ + \sum_{i \neq 0, \, i \in M} 2 \, \tilde{\phi}_i(T) \, \cosh(\mu_i/T),
\end{equation}
\begin{equation} \label{eq:pBev}
p_{B(\bar{B})}^{\rm ev}  =  \sum_{i \in B} \tilde{\phi}_i(T) \, \exp(\pm \mu_i/T) \,
\exp\left( \frac{- b \, p_{B(\bar{B})}^{\rm ev}}{T} \right)~.
\end{equation}
Here, $i \in M$ corresponds to mesons~($B_i = 0$), $i \in B$ corresponds to baryons~($B_i = 1$), $b$ is the baryon excluded volume parameter, and $\tilde{\phi}(T)$ is given in Eq. (\ref{eq:press_phi}).
Equation~\eqref{eq:pBev} can be solved in terms of the Lambert W function~\cite{NoronhaHostler:2012ug,Taradiy:2019taz}:
\begin{equation}\label{eq:pBevW}
p_{B(\bar{B})}^{\rm ev}  = \frac{T}{b} \, W[\varkappa_{B(\bar{B})}(T,\mu_B,\mu_Q,\mu_S)],~
\end{equation}
where
\begin{equation}
\varkappa_{B(\bar{B})}(T,\mu_B,\mu_Q,\mu_S) = b \, \sum_{i \in B} \tilde{\phi}_i(T) \, \exp(\pm \mu_i/T).
\end{equation}
The explicit form Eq.~\eqref{eq:pBevW} for the pressure in the EV-HRG model in terms of the Lambert W function allows us to forgo solving the transcendental equation for the pressure.

One should note that different formulations of the excluded volume HRG model exist in the literature.
In the formulation that we use the EV interactions are introduced only for baryon-baryon and antibaryon-antibaryon pairs with a common EV parameter $b$.
This is consistent with the model used in Refs.~\cite{Vovchenko:2016rkn,Vovchenko:2017xad}.
In other EV models~\cite{Andronic:2012ut,Alba:2017bbr,Motornenko:2020yme} the repulsive interactions are introduced for \emph{all} hadron pairs, and each hadronic species may be characterized by its own value of the EV parameter.
As such, these types of EV models contain many more free parameters, while ours is a minimalistic approach.
The EV effects influence the thermodynamics differently in the latter class of EV models, and the conclusions obtained in this work within the former class of EV models do not necessarily translate.

\section{Susceptibilities}
\label{sec:susc}

The two HRG model extensions introduced above -- the excluded volume corrections and the inclusion of extra states -- are complementary to one another.
This can be seen in the following way.
Adding extra resonances can be interpreted as adding attractive interactions among hadrons that lead to the formation of these resonances.
Given the fact that all the extra states have baryon number equal to either 0 or $\pm 1$, this may correspond to meson-meson and meson-baryon interactions, but not to baryon-baryon interactions.
On the other hand, the EV corrections considered here correspond to repulsive baryon-baryon interactions but not to any meson-meson or meson-baryon interactions.
Therefore, the two extensions describe different physics, and thus, they can and should be considered simultaneously.

Both extensions affect the equation of state.
For instance, the inclusion of extra states increases the pressure at a given temperature and chemical potential, while the EV interactions lead to its suppression.
It can be challenging, then, to constrain the two effects separately.
In order to achieve those constraints, we study differential observables that have recently been obtained on the lattice, namely the susceptibilities of conserved charges:
\begin{equation} \label{eq:chis}
\chi^{BQS}_{lmn} = \frac{\partial^{l+m+n}(p/T^4)}{\partial(\mu_{B}/T)^{l} \, \partial(\mu_{Q}/T)^{m} \, \partial(\mu_{S}/T)^{n}}
\end{equation}

The susceptibilities in Eq.~\eqref{eq:chis} can be calculated explicitly in the EV-HRG model with extra states by utilizing Eq.~\eqref{eq:pBevW} and the known properties of the Lambert W function.
One can then construct specific combinations of susceptibilities that are mainly sensitive to either the extra states or baryon excluded volume, but not to both.

\subsection{Extra states from baryon correlators}

Firstly,  we consider the following ratios of second order susceptibilities: $$ \chi_{11}^{BQ} / \chi_2^B, \quad \chi_{11}^{BS} / \chi_2^B.$$
By calculating them explicitly one obtains
\begin{align} \label{eq:chi11BQ}
\frac{\chi_{11}^{BQ}}{\chi_2^B} & = \frac{\sum_{j \in \rm{sectors}} \, B_j \, Q_j \, \tilde{\phi}_j (T) }{\sum_{j \in \rm{sectors}} \, B_j^2 \, \tilde{\phi}_j (T)}, \\
 \label{eq:chi11BS}
\frac{\chi_{11}^{BS}}{\chi_2^B} & = \frac{\sum_{j \in \rm{sectors}} \, B_j \, S_j \, \tilde{\phi}_j (T) }{\sum_{j \in \rm{sectors}} \, B_j^2 \, \tilde{\phi}_j (T)}.
\end{align}
The excluded volume parameter $b$ cancels out in this combination of susceptibilities.\footnote{A similar cancellation has been observed for the $\chi_{11}^{BQ}/\chi_2^B$ ratio in Ref.~\cite{Vovchenko:2017uko} in the framework of the van der Waals HRG model.}
On the other hand, the ratios are sensitive to the particle list encoded in the ``partial pressures'' $\tilde{\phi}_j$.
Therefore, these ratios can be used to constrain the hadronic spectrum.
In particular, it follows from Eqs.~\eqref{eq:chi11BQ} and~\eqref{eq:chi11BS} that $ \chi_{11}^{BQ} / \chi_2^B$ and $ \chi_{11}^{BS} / \chi_2^B$ probe the fractions of charged baryons and hyperons, respectively, compared to all baryons.

\begin{figure*}
    \centering
    \includegraphics[width=0.48\linewidth]{chi11BSchi2B.pdf}
    \includegraphics[width=0.48\linewidth]{chi11BQchi2B.pdf}
    \caption{Temperature dependence of second order susceptibility ratios $\chi_{11}^{BS}/\chi_2^B$ and $\chi_{11}^{BQ}/\chi_2^B$. Continuum extrapolated lattice results from Refs. \cite{Bellwied:2019pxh, Bellwied:2015lba} are shown by the black points with error bars, while the  calculations within the EV-HRG model are the curves for different hadronic lists. This combination of susceptibilities leads to the cancellation of the excluded volume parameter, $b$.}
    \label{fig:chi-BS-2B_no_EV_dep}
\end{figure*}

The ratios $\chi_{11}^{BQ}/\chi_2^B$ and $\chi_{11}^{BS}/\chi_2^B$ are shown in Fig.~\ref{fig:chi-BS-2B_no_EV_dep}.
The inclusion of extra states from PDG2016+ and QM leads to the enhancement of $-\chi_{11}^{BS}/\chi_2^B$ and the suppression of $\chi_{11}^{BQ}/\chi_2^B$.
This is driven by the fact that the extra states are mainly hyperons, thus their addition increases the fraction of strange baryons~(probed by $-\chi_{11}^{BS}/\chi_2^B$) and decreases the fraction of non-strange baryons~(probed by $\chi_{11}^{BQ}/\chi_2^B$).
The comparison with continuum extrapolated lattice data \cite{Bellwied:2019pxh, Bellwied:2015lba} suggests the need for additional hyperon states from PDG2016+/QM, as previously discussed in Ref. \cite{Alba:2017mqu}.
The best agreement with the lattice data is obtained for the QM list.

None of the considered particle lists allow us to describe the lattice data within errors at $T \gtrsim 155$~MeV.
There are several possibilities which might explain these deviations.
A possible explanation, which would not be captured by our model or any modifications on it, is that this temperature corresponds to the onset of deconfinement, at which new degrees of freedom (quarks) start to be liberated. Otherwise, an improvement in the agreement between lattice results and the model could be obtained through one of the following considerations.
If there are even more strange baryons than predicted by the QM, this could improve the agreement with the lattice data.
For instance, the presence of broad, high-mass Hagedorn states~\cite{Hagedorn:1965st} may have a considerable effect on the susceptibilities as one approaches the Hagedorn temperature $T_H \sim 160-180$~MeV~\cite{Lo:2015cca}.
If the Hagedorn states contain more strange baryons than non-strange baryons, this may improve the agreement with the lattice data in Fig.~\ref{fig:chi-BS-2B_no_EV_dep}.
However, it might be challenging to preserve at the same time the agreement with the individual susceptibilities rather than in the ratios alone.

Other explanations would go beyond the physics of the model employed in the present paper.
For instance, we have modeled all resonances as free particles with zero width.
On the other hand, many non-strange baryon resonances like $\Delta$'s and $N^*$'s are broad, thus a proper treatment of their spectral functions should be important.
Modeling of broad resonances is challenging, analyses in the literature based on either pion-nucleon scattering phase shifts within the S-matrix approach~\cite{Lo:2017lym} or energy-dependent Breit-Wigner widths~\cite{Vovchenko:2018fmh} indicate that partial pressures of such resonances might be overestimated in the standard HRG model.
This implies that a more involved treatment of broad resonances may lead to a suppressed $\chi_{11}^{BQ}/\chi_2^B$ ratio~(and hence an enhanced $-\chi_{11}^{BS}/\chi_2^B$) and recover the agreement with the lattice data.

Finally, the comparison with the lattice data may be affected if there is a flavor hierarchy in baryon excluded volumes.
While the excluded volume effects cancel out in $\chi_{11}^{BQ}/\chi_2^B$ and $\chi_{11}^{BS}/\chi_2^B$ ratios when a common EV parameter $b$ is used for all baryons, this would no longer be the case if excluded volumes differ between strange and non-strange baryons.
A smaller EV for strange baryons would lead to a smaller suppression of $\chi_{11}^{BS}$ relative to $\chi_{11}^{BQ}$, thus leading to an improved agreement with the lattice data.
It is possible that a combination of the three effects discussed here is at play, and it would be interesting to study these in more detail in the future.
We indeed find indications for the flavor-dependent excluded volumes in the behavior of fourth-order susceptibilities discussed in the following subsection.

\subsection{Fourth-order cumulants and excluded volume}

In addition to the extra states, we also want to place limits on the excluded volume parameter, $b$.
This can be done by considering ratios of fourth-to-second order susceptibilities.
The following three ratios are all equal in the EV-HRG model under consideration and sensitive to the EV parameter $b$:
\begin{align}
\frac{\chi^B_4}{\chi^B_2} = \frac{\chi^{BS}_{31}}{\chi^{BS}_{11}} =  \frac{\chi^{BQ}_{31}}{\chi^{BQ}_{11}}
& = \frac{1 - 8 \, W(\varkappa_B) + 6 [W(\varkappa_B)]^2}{[1 + W(\varkappa_B)]^4} \nonumber \\
& = 1 - 12 \varkappa_B + O(\varkappa_B^2).
\end{align}
In the ideal HRG model, i.e. without the EV interactions, $\varkappa_B = 0$ so that $\frac{\chi^B_4}{\chi^B_2} = \frac{\chi^{BS}_{31}}{\chi^{BS}_{11}} = \frac{\chi^{BQ}_{31}}{\chi^{BQ}_{11}} = 1$ regardless of the inclusion of any additional hadronic states.
The suppression of these ratios relative to unity, on the other hand, is directly sensitive to the EV interactions and can be used to constrain the EV parameter $b$.
Furthermore, the fact that all three ratios are predicted to be equal within the model allows us to probe the limits of validity of the model, which would be signaled by the point where the equality among these three ratios no longer holds in the lattice data.
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{4th-to-2nd-all.pdf}
    \caption{Temperature dependence of fourth-to-second order susceptibility ratios $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$, predicted to be equal in the EV-HRG model. The lattice data at finite lattice spacing $N_{\tau}=12$ from Ref. \cite{Borsanyi:2018grb} are shown as grayscale symbols with error bars, while the calculations within the EV-HRG model are shown as bands for a range of excluded volume parameter, $b$, and for different hadronic lists. The ideal HRG result is given by the horizontal line at unity.
    }
    \label{fig:chi-4-2_equiv_EV}
\end{figure}

Figure~\ref{fig:chi-4-2_equiv_EV} depicts the results of the calculation of the ratios $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$ and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ within the EV-HRG model for range $b = 0.4-1$~fm$^3$ of the EV parameter values.
The three ratios all coincide with one another, as expected, and exhibit
minimal dependence on the hadronic list utilized.
Therefore, these particular quantities are indeed sensitive mainly to the excluded volume repulsive interactions rather than to the hadronic spectrum used in the HRG model.
In Fig. \ref{fig:chi-4-2_equiv_EV} we compare the calculations in the EV-HRG model with various particle lists to lattice data at $N_{\tau} = 12$ from the Wuppertal-Budapest collaboration~\cite{Borsanyi:2018grb}, in the temperature range $T = 135-170$~MeV.
Since not all of the available lattice data are continuum extrapolated, we choose this larger lattice spacing for the comparisons and avoid constructing ratios that would be a mixture of results at finite $N_{\tau}$ and in the continuum limit.
The qualitative behavior of the three ratios is very similar in the whole temperature range considered.
Quantitatively, we see that at temperatures below 150 MeV the three ratios sit on top of each other.
This is expected, although the lattice error bars are relatively sizable at those temperatures.

Statistically significant differences between the three susceptibility ratios in the lattice data emerge at $T \gtrsim$~160~MeV.
On the one hand, this may be a reflection of the transition to the deconfined phase where a hadronic model would be expected to break down.
In the Stefan-Boltzmann limit of massless quarks one has $\chi_4^B/\chi_2^B = \chi_{31}^{BS}/\chi_{11}^{BS} \to 2 / (3\pi^2)$ while $\chi_{31}^{BQ}/\chi_{11}^{BQ} \to 0$.
Thus, the notably smaller values of $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ at $T \gtrsim 160$~MeV compared to the other two ratios might be related to the smaller Stefan-Boltzmann limit for this quantity.
On the other hand, these differences may also reflect a flavor hierarchy in baryon excluded volumes.
The ratios are predicted to be equal in the EV-HRG model if all baryons are assigned a common EV parameter $b$.
However, if for example strange baryons have a different~(smaller) excluded volume, one would expect $\chi_{31}^{BS}/\chi_{11}^{BS}$ to exhibit smaller deviations from the baseline of unity than the other ratios.
We see in Fig.~\ref{fig:chi-4-2_equiv_EV} that the separation between the strangeness and electric charge susceptibilities is such that they lie on the higher and lower ends of the $b = 0.4-1$~fm$^3$ band, respectively.
This is an indication that the strange baryons may indeed have a smaller volume than the non-strange ones.
From a phenomenological point of view, the smaller hyperon volume could reflect the fact that hyperon-hyperon interactions are mediated by the exchange of heavier mesons like $\phi$ compared to nucleon-nucleon interactions which correspond to the exchange of lighter, non-strange mesons like $\sigma$ and $\omega$.
The implementation of smaller excluded volumes for strange particles is possibly not unique. Such extensions of the EV-HRG model have been considered in Refs.~\cite{Alba:2016hwx,Vovchenko:2017zpj,Alba:2017bbr,Motornenko:2020yme} but are not considered in the present study. Of course it can be done in the future in order to model the subtle differences between $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$, for example at the chemical freeze-out in heavy-ion collisions.

\subsection{Effect of the extra states and excluded volume on other susceptibilities}

Next, we investigate the following two combinations of susceptibilities that are sensitive to the extra strange states: the kurtosis of net-strangeness fluctuations $\chi_4^S / \chi_2^S$ and the correlator $\chi_{11}^{us}$ between net numbers of up and strange quarks.
The effect of extra states on these quantities was investigated in Ref.~\cite{Alba:2017mqu} without excluded volume effects.
It was shown there that extra states improve the description of $\chi_{11}^{us}$ but in the case of the QM list spoil the agreement with the lattice data for $\chi_4^S/\chi_2^S$.
Here we investigate how these quantities are affected by the presence of baryon excluded volume in addition to extra states.

Figure~\ref{fig:chi-4S-2S} depicts the temperature dependence of $\chi_4^S/\chi_2^S$.
This quantity does not involve any $\mu_B$ derivatives, thus it does not probe the repulsive baryonic interactions as directly as e.g. $\chi_4^B/\chi_2^B$ that we considered before.
Due to that fact, we calculate over a broader range of $b = 0-1$~fm$^3$ for the EV parameter value, given by the bands in Fig.~\ref{fig:chi-4S-2S}.
Since there are no $\mu_B$ derivatives, both hyperons as well as strange mesons like kaons contribute to this quantity.
Moreover, the values of $\chi_4^S/\chi_2^S$ above unity are due to multi-strange hyperons.
Figure~\ref{fig:chi-4S-2S} shows that the excluded volume suppresses $\chi_4^S/\chi_2^S$. This effect is less pronounced at smaller temperatures, but becomes sizable at $T \gtrsim 150$~MeV.
We see that the lattice data are underestimated when using the standard PDG2016 list, and the excluded volume does not improve the agreement.
For the PDG2016+ list the agreement with the lattice data is obtained at $T \lesssim 160$~MeV for the smallest considered values of the excluded volume~(the upper part of the blue band).
This is consistent with the observation made earlier that strange baryons prefer a smaller excluded volume.
On the other hand, when the QM list is considered, which contains even more extra strange baryons, the best agreement with the lattice data is obtained for the higher value of the EV parameter, $b \simeq 1$~fm$^3$~(the lower part of the red band).
We see similar effects for  $\chi^{us}_{11}$, shown in Fig.~\ref{fig:chi-us}, with a slight over-prediction in the case of the QM list.
The PDG2016 list describes the lattice data at $T \lesssim 150$~MeV, but breaks down at higher temperatures, with no benefit from introducing the excluded volume.
The PDG2016+ and QM lists allow us to extend the agreement with the lattice data for $\chi^{us}_{11}$ to $T = 160-165$~MeV when baryon excluded volumes of up to $b \simeq 0.4$~fm$^3$ and $b \simeq 1$~fm$^3$ are used, respectively.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{chi4S2S-all.pdf}
    \caption{Temperature dependence of fourth-to-second order strangeness susceptibility ratio $\chi_4^S/\chi_2^S$. Continuum extrapolated lattice results for this quantity are given by black points with error bars \cite{Alba:2017mqu} and EV-HRG calculations are shown for the full range of the parameter $b$ for different hadronic lists.}
    \label{fig:chi-4S-2S}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{chi11us_cont_ext_b-all.pdf}
    \caption{Up-strange quark susceptibility, $\chi^{us}_{11}$, as a function of the temperature. The continuum extrapolated lattice results are given as black points with error bars \cite{Alba:2017mqu}, while the EV-HRG calculations are shown for a range of the parameter $b$ between 0.4 and 1 $\text{fm}^3$ for different hadronic lists.}
    \label{fig:chi-us}
\end{figure}

Based on the analysis of $\chi_4^S/\chi_2^S$ and $\chi^{us}_{11}$, we observe a correlation between the number of extra strange states and baryon excluded volume: the larger the number of extra states, the larger the excluded volume must be in order to recover agreement with the lattice data for these two quantities.
Thus the attractive interactions via the inclusion of extra states can be balanced by an additional repulsion in the baryon sector in these observables.
One way to break this degeneracy is to consider quantities that probe only a single one of the effects, like the $\chi_{11}^{BQ}/\chi_2^B$ and $\chi_{11}^{BS}/\chi_2^B$ ratios that we considered earlier.
However, those ratios may be sensitive to additional physics like modeling of broad resonances, as discussed above.
Instead, we study here  a couple of additional susceptibilities that probe the strangeness content of baryons.
First we look at $\chi^{BS}_{22}$, which is a doubly strange quantity and thus more sensitive to multi-strange hyperons.
This quantity, shown in the left panel of Fig.~\ref{fig:chi-22BS}, paints a picture consistent with $\chi_4^S/\chi_2^S$ and $\chi^{us}_{11}$: extra states are clearly needed to describe the lattice data, but PDG2016+ prefers a smaller excluded volume while the QM list prefers a larger one.
On the other hand, the mixed $BQS$ susceptibility containing only one strangeness derivative $\chi^{BQS}_{211}$, shown in the right panel of Fig.~\ref{fig:chi-22BS}, exhibits only mild dependence on the excluded volume
but large sensitivity to the number of extra states.
The minimal dependence on the excluded volume parameter $b$ of this fourth order susceptibility can be understood when one considers that the derivatives with respect to $Q$ and $S$ are only single derivatives.
Some baryons, e.g. $\Sigma^+$, carry these conserved charges with opposite signs~(see Table~\ref{tab:hadron_table}).
For similar reasons, the contributions of the various baryon-baryon interactions to $\chi^{BQS}_{211}$ can be either positive or negative and as such, each term making up this quantity, as evaluated from Eqs. \eqref{eq:pBevW} and \eqref{eq:chis}, can contribute with different signs.
This leads to a smaller $b$-dependence of $\chi^{BQS}_{211}$ than that of $\chi^{BS}_{22}$~(the left panel of Fig.~\ref{fig:chi-22BS}), in which case all contributions carry the same sign due to the fact that both the baryon number $B$ and the strangeness $S$ in this observable are squared.
The lattice results tend to lie in between the predictions based on the PDG2016+ and QM lists.
This indicates that the number of extra states might be underestimated in the PDG2016+ list but overestimated in the QM list.
Thus, the difference between the two lists may be taken as a systematic uncertainty in the particle list.

\begin{figure*}
    \centering
    \includegraphics[width=0.48\linewidth]{chi22BS-all.pdf}
    \includegraphics[width=0.48\linewidth]{chi211BQS-all.pdf}
    \caption{Temperature dependence of mixed fourth order susceptibilities $\chi^{BS}_{22}$ and $\chi^{BQS}_{211}$. The results at finite lattice spacing $N_{\tau}=12$ from Ref. \cite{Borsanyi:2018grb} are shown as black points with error bars. Left: Fourth-order baryon-strangeness susceptibility calculated with several different hadronic lists in the EV-HRG model for the full range of the EV parameter $b$. Right: Mixed fourth order $BQS$ susceptibility calculated with several different hadronic lists in the EV-HRG model for a range of the EV parameter $b$ between 0.4 and 1~fm$^3$.}
    \label{fig:chi-22BS}
\end{figure*}

We note that the lattice data on
the susceptibilities $\chi^{BS}_{22}$ and $\chi^{BQS}_{211}$, as well as the ratios $\chi_{31}^{BS}/\chi_{11}^{BS}$ and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$, are not continuum extrapolated but only available up to a lattice spacing of $N_{\tau} = 12$. While continuum extrapolated lattice results for $\chi_4^B/\chi_2^B$ have been shown e.g. in Ref. \cite{Borsanyi:2013hza}, in this manuscript we use $N_{\tau}=12$ also for this quantity, for a consistent comparison with the other ratios.
These data could be subject to slight alterations in the continuum limit, in particular at the lower temperatures.
Additionally, there are many more susceptibilities available from the lattice than those explored in this work.
However, in this manuscript we focused mostly on baryon and strangeness observables in order to probe the extra strange states and baryon-baryon interactions.
This allows us to avoid the most severe lattice systematics, for instance taste violation, which mainly affect calculations sensitive to pion degrees of freedom, like e.g. electric charge susceptibilities.

\section{Conclusions}\label{sec:concl}

We have investigated two common extensions of the Hadron Resonance Gas model
that implement additional attractive and repulsive interactions among hadrons.
The attractive interactions correspond to adding extra
states exceeding those measured with high confidence by the Particle Data Group, leading to additive corrections to the overall pressure in the HRG model.
This has been studied with the use of the PDG2016+ and QM particle lists.
On the other hand, we also apply excluded-volume corrections in the baryon sector,
which model the presence of repulsive core in (anti)baryon-(anti)baryon interactions.
We demonstrate that these two extensions are complementary and find support in the available first-principles lattice QCD data.

To constrain the two extensions simultaneously we constructed specific combinations of conserved charge susceptibilities that probe the two effects separately.
We show that the second order ratios, $\chi_{11}^{BS}/\chi_2^B$ and $\chi_{11}^{BQ}/\chi_2^B$, probe only the hadron spectrum but not the excluded volume.
The inclusion of additional states improves the agreement with the lattice data for these two quantities, with the best description obtained using the QM list.
Even for the QM list, however, deviations from the lattice data emerge at $T \sim 150-155$~MeV.
We argued that these deviations may necessitate a more involved modeling of broad resonances as well the possibility of smaller excluded volumes for strange baryons compared to non-strange ones.
We further studied the constraints on the hadronic spectrum by analyzing the various strangeness susceptibilities, including $\chi_4^S/\chi_2^S$, $\chi^{us}_{11}$, $\chi^{BS}_{22}$ and $\chi^{BQS}_{211}$.
The analysis of the lattice data for these quantities indicates that both the PDG2016+ and QM lists are preferable over the standard list, where the former contains most but not all the extra strange states while the QM list contains too many.
Therefore, the difference between the results using these two lists can be taken as a systematic uncertainty due to the hadron spectrum.

The fourth-to-second order ratios, $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ are shown to be mainly sensitive to the excluded volume corrections and, thus, suitable to constrain these corrections.
In the absence of the excluded volume corrections, these three ratios are equal to unity irrespective of the hadronic spectrum.
The excluded volume effects, on the other hand, suppress the ratios and make them behave similarly to the lattice data at $T \sim 155-165$~MeV.
The EV-HRG model that we use assigns a constant EV parameter $b$ for all baryons and predicts $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ to be all equal to one another at a given temperature.
The lattice data, on the other hand, reveal statistically significant differences between the three ratios at $T \gtrsim 155$~MeV, which follow a hierarchy $\chi_{31}^{BS}/\chi_{11}^{BS} > \chi_4^B/\chi_2^B > \chi_{31}^{BQ}/\chi_{11}^{BQ}$.
We argued that this hierarchy indicates a flavor dependence in the baryon excluded volumes, namely that strange baryons have generally smaller excluded volumes than non-strange baryons.
Therefore, while using a constant EV parameter in a range $b = 0.4-1$~fm$^3$ may be good enough to capture the general suppression of $\chi_4^B/\chi_2^B$, $\chi_{31}^{BS}/\chi_{11}^{BS}$, and $\chi_{31}^{BQ}/\chi_{11}^{BQ}$ ratios at $T \sim 155-165$~MeV, a more involved model is necessary to describe the subtle differences between the three.

In summary, an extended HRG model incorporating extra states via either the PDG2016+ or QM list as well as baryon excluded volume with parameter $b = 0.4-1$~fm$^3$ significantly improves the description of many lattice QCD susceptibilities at temperatures up to $T \simeq 160-165$~MeV over the standard HRG model.
We note that the results from these two lists provide an estimate of the theoretical uncertainty on the hadronic list.
This is particularly relevant for the chemical freeze-out conditions realized in heavy-ion collisions and can be used for improved modeling of event-by-event fluctuations measured in the corresponding experiments at the LHC, RHIC, and SPS.
For instance, the EV-HRG model studied here can be directly used in the generalized Cooper-Frye particlization routine developed in Ref.~\cite{Vovchenko:2020kwg}.
We summarize the performance of each list as follows. We find that the QM list performs better for the ratio $\chi_{11}^{BS}/\chi_2^B$, while both lists describe the continuum extrapolated lattice data well for $\chi_{11}^{BQ}/\chi_2^B$ as shown in Fig. \ref{fig:chi-BS-2B_no_EV_dep}.
For ratios of fourth-to-second order susceptibilities as shown in Fig. \ref{fig:chi-4-2_equiv_EV}, we find, as we expect, a small variation in the results within the EV-HRG model for different lists.
From strangeness-sensitive susceptibilities $\chi_4^S/\chi_2^S$, $\chi_{11}^{us}$ and $\chi_{22}^{BS}$, we determined that the PDG2016+ list agrees best with the continuum extrapolated lattice data for a small excluded volume $b \simeq 0.2-0.4$~fm$^3$, while the QM list finds agreement when $b \simeq 1$~fm$^3$.
For the mixed fourth order susceptibility  $\chi_{211}^{BQS}$, the optimal list for treatment of the chemical freeze-out is PDG2016+.
Given that there is still some tension between the PDG2016+ list and the QM list, our results could indicate that more work is still needed to determine the number of states in the hadronic spectrum in heavy-ion collisions.
Further improvements of the model can be achieved by considering differences in excluded volumes of strange and non-strange baryons, as well as a more involved modeling of broad resonances.

\section*{Acknowledgements}
The authors would like to thank Paolo Parotto for discussions and help preparing the updated version of the QM list used in this publication.
This material is based upon work supported by the National Science Foundation under grant no. PHY1654219, by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-AC02-05CH11231 and within the framework of the Beam Energy Scan Theory (BEST) Collaboration.
V.V. acknowledges the support via the
Feodor Lynen program of the Alexander von Humboldt
foundation.


\section{Introduction}
Directed energy phased array (DEPA) systems are preferable to single aperture systems in directed energy applications that require large apertures for long range targets because of their scalability in both power and size. Scalability in power opens up a new mission space with high power requirements, and scalabilty in size allows illumination of targets at greater distances with minimal loss from diffraction. Existing monolithic single mode laser elements are limited to $\sim$10~kW in power \cite{Zervas2014_hplasers} and meter-class aperture sizes at maximum. In principle, a DEPA could be built to have any power or diameter by incorporating the appropriate number of modular subelements. Furthermore, the ability to control the phase and alignment of each laser subelement provides a means to correct for atmospheric turbulence. These features make DEPAs a potential solution to various problems involving energy supply to distant objects on the ground, in the air or in space in scenarios where sunlight is insufficient or inaccessible. For example, sending power to the moon during the lunar night could be accomplished by using an Earth-based DEPA to provide optical power to a high efficiency tuned  photovoltaic array on the moon. Power could also be sent to other extraterrestrial targets such as low Earth orbit, medium Earth orbit, geosynchronous and deep space spacecraft for power or propulsion. The latter could be done via photovoltaic conversion on the spacecraft to power high specific impulse ion engines and enable station keeping or rapid interplanetary transit. Kilometer scale arrays would be capable of diverting incoming asteroids off course by ablation, or accelerating low mass spacecraft to relativistic speeds for interstellar missions \cite{Lubin2016_roadmap,Lubin2020_ESA}. Other applications include but are not limited to deep space communications, powering of ground-based remote sensors and rovers in inhospitable locations, remote sensing LIDAR, laser machining, and experiments involving matter under extreme conditions.

\begin{figure}
    \centering
    \includegraphics[width=1\textwidth]{schematic_bigtext}
    \caption{Schematic diagram of the simulated system with coordinates. A beacon laser signal (right, blue) passes through atmospheric phase screens and illuminates the DEPA (left), which measures the reverse-propagating beacon in fiber for phase conjugation and alignment feedback. The transmit beam (red) originates from a master seed laser, which is split and amplified in fiber before propagating in free space. Here a 91 element hexagonal array is shown, but the DEPA can be scaled to any size or shape.}
\label{fig:schematic}
\end{figure}

While it is feasible to construct DEPA systems for space-based applications in space to avoid atmospheric effects and increase target visibility, cost and accessibility make it far more practical to start with ground based systems. It is therefore of great importance to quantify the systems' effectiveness in correcting for atmospheric perturbations. This study considers a ground-based continuous wave DEPA supplying power to a space-based target. This baseline DEPA is composed of a large number of circular sub-apertures of 10~cm diameter, which are hexagonally close packed so their edges touch each other. The target is equipped with a ``beacon'' laser, which illuminates the DEPA after being perturbed by the atmosphere. The atmosphere is modeled as a large number of two dimensional phase screens at different altitudes that represent the wavefront (phase) errors that a plane wave acquires as it travels in the zenith direction. The DEPA is equipped with sensors to measure the reverse-propagating beacon light and use it for phase retrieval and tip-tilt alignment. For simplicity we assume the target is stationary, making the results applicable to slow moving assets but not necessarily fast moving targets such as low Earth orbit satellites. Fig. \ref{fig:schematic} schematically represents the full system being simulated. The wavelength of both the beacon and transmit beams is assumed to be 1.064~$\mu$m to take advantage of developed ytterbium fiber laser amplification technology and low atmospheric absorption \cite{GeminiObs_transmission,Lubin2020_ESA}, though the model can easily be adapted to other wavelengths such as 1.55~$\mu$m (erbium-doped fiber amplifiers) and 2~$\mu$m (thulium-doped fiber amplifiers). This paper outlines how noise in such a system is modeled and provides simulated results for performance under various probable conditions at astronomical sites.

We utilize a recently developed DEPA beam simulation tool described in \cite{Hettel2019_spie} to calculate the transmitted beam profile. Conventional propagation methods, such as the fast Fourier transform based angular spectrum method \cite{Novotny2012_prop}, are not well suited to handle large numbers of array elements with uniquely defined perturbations on each subelement. The computation time and memory requirements for parameterizing the near field DEPA beam to sufficient resolution are significant, and any FFT based algorithm is poorly scalable with array size unless many spatial measurement points of the beam profile at the target are needed. Our method works by summing the analytic electric field solutions of each subelement at a given measurement point. This operation makes parameterizing the near field beam trivial, and for $N$ subelements with a fixed number of measurement points it requires $\mathcal{O}(N)$ calculations. While this study deals with a modest number of subelements, the incorporation of temporal dynamics makes scalability of critical importance to enable swift calculation of the beam profile at many points in time. We use this infrastructure to simulate extremely large numbers of elements ($N>10^{10}$, limited only by availability of computing resources) to calculate the beam profile at multiple points in time in parallel.

\section{Simulation pipeline}
Fig. \ref{fig:pipeline} summarizes the steps of a DEPA atmospheric simulation and how they connect to each other. First, atmospheric phase screens are generated with the Python extension SOAPY \cite{Reeves2016_soapy} using a modified Hufnagel-Valley (HV) model for atmospheric turbulence, and a modified Bufton model for wind speed (see \textcolor{urlblue}{Supplement 1}, Section 2). The modifications of both models are informed by site survey data taken for the Thirty Meter Telescope (TMT) \cite{Schock2009_tmtdata}. The phase screens are used to optically propagate the beacon to the DEPA and produce its intensity and phase as a function of space and time. The data is then mapped to the locations of the DEPA sub-apertures to produce time series for piston phase, tip-tilt error, incident beacon power, intensity reduction due to higher order aberrations, and time-of-flight error for each subelement. A model for tip-tilt servo performance is used to generate tip-tilt residual error, which is combined with incident beacon power and higher order aberrations to calculate the beacon power collected by each subelement. Beacon power and piston phase values, along with detector noise and data-based fiber amplifier phase noise, are used as inputs for simulation of the stochastic parallel gradient descent (SPGD)  \cite{cbc2013_spgd} phasing servomechanism. SPGD is an algorithm routinely used in adaptive optics systems to phase lock numerous subelements. The SPGD simulation generates residual piston phase errors for the transmit beam. The phase residuals, tip-tilt residuals, intensity reduction due to higher order aberrations, and time-of-flight effects are used as inputs for the transmit beam propagator to produce final results for the beam at the target. Each step of the process is described in further detail in \textcolor{urlblue}{Supplement 1}, with exception of the optical propagation steps in the following section.

\begin{figure}
    \centering
    \includegraphics[width=1\textwidth]{pipeline_bigtext}
    \caption{Simulation flow diagram including the various processes and noise inputs.}
\label{fig:pipeline}
\end{figure}

\section{Optical propagation of beacon and transmit beam}
\label{sec:propagation}
The beacon is modeled as a plane wave propagating from the highest altitude phase screen to the ground. The electric field is iteratively calculated at the location of each phase screen to account for scintillation effects. Scintillation arises as an effect of the atmosphere refracting or scattering light in different directions, and the light interfering with itself as it propagates. The scintillation index $\sigma_I^2$ represents the normalized intensity variance with respect to time at a given point in space, and is defined as \cite{Andrews2001_scint}:
\begin{equation}
    \sigma_I^2 = \frac{<I^2>}{<I>^2} - 1.
\end{equation}
The plane wave model of the beacon makes initial parameterization of its electric field $\mathbf{E}_\mathrm{B}$ trivial, and we require a large number of spatial measurement points to accurately map it to the DEPA sub-apertures. These features make the angular spectrum method an appropriate propagator for $\mathbf{E}_\mathrm{B}$ \cite{Novotny2012_prop}:
\begin{equation}
\label{eq:angularspec}
\mathbf{E}_\mathrm{B}(x,y,z) = \iint_\infty^\infty \hat{\mathbf{E}}_\mathrm{B}(k_x,k_y;0)e^{i[k_xx + k_yy +k_zz]}dk_xdk_y
\end{equation}
where
\begin{eqnarray}
\hat{\mathbf{E}}_\mathrm{B}(k_x,k_y;z) &=& \frac{1}{4\pi^2}\iint_\infty^\infty \mathbf{E}_\mathrm{B}(x,y,z)e^{i[k_xx + k_yy]}dxdy \\
k_z &=& \sqrt{k^2 - k_x^2 - k_y^2} \ \ \ \approx \ \ \ k - \frac{k_x^2 + k_y^2}{2k}.
\label{eq:k}
\end{eqnarray}
Here $x$, $y$, and $z$ are the spatial coordinates in Fig. \ref{fig:schematic}, $\hat{\mathbf{E}}_\mathrm{B}$ is the Fourier transform (or angular spectrum) of $\mathbf{E}_\mathrm{B}$, and $\boldsymbol{k}$ is the wave vector of the propagating light. $\hat{\mathbf{E}}_\mathrm{B}$ is propagated in the $z$ direction by the $k_zz$ term in Eq. (\ref{eq:angularspec}), and inverse Fourier transformed to solve for $\mathbf{E}_\mathrm{B}$. The paraxial approximation is employed in Eq. (\ref{eq:k}). Both $\hat{\mathbf{E}}_\mathrm{B}$ and $\mathbf{E}_\mathrm{B}$ are efficiently calculated with fast Fourier transforms. Fig. \ref{fig:beaconprop} is an example of the phase and intensity of a beacon of 1.064 $\mu$m wavelength at a given point in time after it propagates to the DEPA through an atmosphere with Fried length $r_0$ = 10~cm (for 500~nm wavelength) defined by the turbulence model discussed in \textcolor{urlblue}{Supplement 1}. All values of $r_0$ in this paper unless otherwise specified are at a wavelength of 500 nm. $\sigma_I^2 \approx0.2$, in agreement with theoretical and experimental results \cite{Thomas2005_scint,Schock2009_tmtdata}.

\begin{figure}
    \centering
        \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\textwidth]{phasescreen_bigtext}
        \caption{Simulated phase of beacon as it reaches the DEPA.}
        \label{fig:phasescreen}
    \end{subfigure}
    ~
    \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\linewidth]{beaconprofile_bigtext}
        \caption{Simulated intensity of beacon as it reaches the DEPA.}
        \label{fig:beaconprofile}
    \end{subfigure}

    \caption{Simulated spatial profiles of phase (Fig. \ref{fig:phasescreen}) and intensity (Fig. \ref{fig:beaconprofile}) of a plane wave at 1.064 $\mu$m wavelength after propagating to the DEPA through an atmosphere with $r_0$ = 10~cm (for 500~nm wavelength) defined by the turbulence model discussed in \textcolor{urlblue}{Supplement 1}.}
\label{fig:beaconprop}
\end{figure}

The transmit beam of the DEPA is propagated using the method described in \cite{Hettel2019_spie}, which models each laser subelement as a Gaussian beam truncated by a circular aperture. The electric fields $\mathbf{E}_i$ of $N$ subelements are found with their respective analytic solutions and added together to construct the synthesized electric field $\mathbf{E}_\mathrm{T}$:
\begin{equation}
\label{eq:transmitfield}
\mathbf{E}_\mathrm{T}(x,y,z) = \sum_{i=1}^N \mathbf{E}_i(x,y,z).
\end{equation}
This method enables rapid parameterization of large DEPAs (above $10^{10}$ subelements demonstrated in \cite{Hettel2019_spie}) with perturbations in phase, amplitude, position, output angle, and focal distance. Customized parallel computing allows for their electric fields to be quickly calculated on a GPU to great accuracy. In this study, we used a GeForce GTX 1070 and Titan-Z desktop GPU with double point precision.

The atmosphere is approximated to be ``thin'' when calculating the transmit beam profile, i.e. the beam encounters all atmospheric turbulence in the near field before propagating into the far field or focal plane. This simplification is valid for the discussed simulations since the total propagation distance is much larger than the thickness of the atmosphere \cite{Brooker1985_thinatm,laserprop2005_thinatm}. Another  notable simplification is that the propagator does not model subelement aberrations of higher order than focal shifts. Losses in the main beam at the target due to these aberrations are still accounted for by reduction of the on-axis intensity of the subelements by the appropriate amount (see \textcolor{urlblue}{Supplement 1}). Since all transmit beam calculations are done in the far field of the array, any differing effects will be contained in side lobes far outside our region of interest.

\section{Numerical results for transmit beam}

\begin{figure}
    \centering
        \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\textwidth]{arrayperformance_bigtext}
        \caption{Time series of simulated measured beacon power, Strehl ratio of transmit beam, and combining efficiency of transmit beam.}
        \label{fig:timeseries}
    \end{subfigure}
    ~
    \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\linewidth]{arrayperformancezoom_bigtext}
        \caption{Same as left, but zoomed in on the phase convergence. The rise time is $\sim$2~ms, or about two cycles through the dither matrix as predicted for the SPGD algorithm in \textcolor{urlblue}{Supplement 1}.}
        \label{fig:timeserieszoom}
    \end{subfigure}

    \caption{Simulation output time series with $r_0$ = 15~cm (at 500~nm wavelength), $v_\mathrm{g}$ = 5~m/s and $\delta$ = 0~deg. Fig. \ref{fig:timeseries} shows the entire duration of the simulation, while Fig. \ref{fig:timeserieszoom} is zoomed in on the initial phase convergence. The beacon power is normalized to the amount that would be in the system with perfect alignment and no atmosphere. One can see that combined beacon power, output Strehl, and combining efficiency are well correlated with one another. With Strehl > 0.8, the beam is considered diffraction limited by the Maréchal criterion \cite{Born1975_marechal}.}
\label{fig:outputs}
\end{figure}

Simulations were performed for an array composed of $N=1024$ circular sub-apertures packed hexagonally with touching edges to form a disk. The sub-apertures were 10~cm in diameter, making the full array about 3.4~m in diameter. The transmit and beacon wavelengths were both 1.064~$\mu$m. While this is a specific configuration, the modularity of the subelements and known scaling relationships make the results applicable to a DEPA of any size (with the exception of piston phase errors) and targets at any distance. The measurement plane was set to be $4\times10^8$~m away from the array (approximately the distance to the moon) and was 500x500~m with 5~m resolution, which provides about 1\% error in estimating combining efficiency. Combining efficiency is defined as the ratio of optical power within an on-axis diffraction limited spot to the total power emitted by the array, and is a measure of the fractional power on the target. The spot diameter in angular units for a circular array of diameter $D$ and wavelength $\lambda$ is $2.44\lambda/D$, the first minimum of the Airy disk. For an array composed of circular sub-apertures which collimate Gaussian beams, the maximum achievable combining efficiency is $\sim$0.62 \cite{Hettel2019_spie,cbc2013_combining}. This limit is due to a combination of diffraction effects and power clipping at the apertures.

\begin{figure}
    \centering
    \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\linewidth]{avbeamprofile_bigtext}
        \caption{Simulated time averaged intensity of transmit beam in the target plane.}
        \label{fig:trasmitprofile}
    \end{subfigure}
    ~
    \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\linewidth]{peaklocation_bigtext}
        \caption{Peak location of simulated beam as a function of time.}
        \label{fig:peaklocation}
    \end{subfigure}

    \caption{Fig. \ref{fig:trasmitprofile} is a time average of the transmitted beam profile, normalized to its peak intensity with perfect transmission ($T$ = 1) and Strehl = 1. The color bar also functions as the y axis for the line plot. Fig. \ref{fig:peaklocation} shows the location of the peak as a function of time, normalized to the beam diameter as measured to first null. The wind is in the $x$ direction which causes larger deviations in $x$ than in $y$. Both plots are for $r_0$ = 15~cm, $v_\mathrm{g}$ = 5~m/s and $\delta$ = 0.}
\label{fig:spatialprofile}
\end{figure}

Fig. \ref{fig:outputs} shows measured beacon power, the Strehl ratio of the transmit beam, and combining efficiency of the transmit beam for Fried length $r_0$ = 15~cm (at 500~nm wavelength), ground wind speed $v_\mathrm{g}$ = 5~m/s and zenith angle $\delta$ = 0~deg as a function of time. One can see that they are in close agreement with each other, taking into account their differing normalization factors. Here the Strehl ratio is defined as the ratio of peak on-axis intensity to that with no perturbations and does not take atmospheric transmission into account. Measured beacon power, on the other hand, is normalized to the power in a system with perfect alignment and no atmosphere. The average steady state Strehl is above 0.8, which is considered to be diffraction limited by the Maréchal criterion \cite{Born1975_marechal}. Fig. \ref{fig:spatialprofile} displays the time-averaged spatial profile and peak location as a function of time for the same beam. Spatial coordinates are in angular units to make them independent of propagation distance. The beam diameter is inversely proportional to array diameter, as the main lobe of a properly phased DEPA has the same geometric properties as a single aperture of the same size and shape \cite{Hettel2019_spie}. In this case the beam diameter is approximately 0.16~arcseconds (or 308~m for our propagation distance of $4\times10^8$~m). The peak location was found by performing a Gaussian fit around the maximum intensity value. The fitting procedure increases the resolution of the location values, which allows for subtle motion to be seen in the beam. See \textcolor{urlblue}{Visualization 1} for a video of the beam's spatial mode in these conditions, and \textcolor{urlblue}{Visualization 2} for a video in the harshest conditions simulated.

Figs. \ref{fig:strehlvszen} and \ref{fig:combeffvszen} summarize Strehl and combining efficiency results for all of the simulations performed. Average Strehl and combining efficiency for each of the atmospheric models are displayed as a function of zenith angle. Note that for a given atmospheric model, $r_0$ decreases with zenith angle and the specified values of $r_0$ are for $\delta$ = 0~deg. The beam remained diffraction limited for $r_0$ down to 10~cm and $\delta$ up to 30-60~deg, depending on the model (see \textcolor{urlblue}{Visualization 3-4} for spatial mode). These conditions cover over 70\% of the site survey data taken for the TMT \cite{Schock2009_tmtdata}, a promising implication for high duty cycle applications of DEPAs at high altitude sites.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{strehlvszen_bigtext}
    \caption{Strehl as a function of zenith angle for each atmospheric condition. Error bars represent the standard deviation over the 500~ms time series. The beam stays diffraction limited (Strehl > 0.8) for $r_0\geq10$~cm and $\delta\leq$ 30-60~deg depending on atmospheric conditions, where $r_0$ labels are specified for 500~nm wavelength and $\delta$ = 0~deg.}
\label{fig:strehlvszen}
\end{figure}
\begin{figure}[t!]
    \centering
    \includegraphics[width=0.7\linewidth]{combeffvszen_bigtext}
    \caption{Combining efficiency as a function of zenith angle for each atmospheric condition. Error bars represent the standard deviation over the 500~ms time series. The dashed line is the theoretical maximum combining efficiency for an array composed of Gaussian beams transmitted through circular apertures. $r_0$ labels are specified for 500~nm wavelength and $\delta$ = 0~deg.}
    \label{fig:combeffvszen}
\end{figure}

The root-mean-square peak location as a function of zenith angle for each atmospheric condition is plotted in Fig. \ref{fig:rmspeak}. There is not a clear, intuitive relationship between $r_0$ and the magnitude of drift, as it is a function of both the strength and spatial structure of the residual phases. The beam drifts more as the linear component of the residual phase profile increases, which in some cases may be less prominent in higher turbulence than in lower turbulence. For a given atmospheric turbulence structure, however, the beam drifts more as both wind speed and zenith angle increase, as one might expect.

\begin{figure}
    \centering
    \includegraphics[width=0.7\linewidth]{rmspeak_bigtext}
    \caption{RMS peak location over 500~ms as a function of zenith angle for each atmospheric condition. The time series for $r_0$ = 15~cm, $v_\mathrm{g}$ = 5~m/s, and $\delta$ = 0~deg is plotted in Fig. \ref{fig:peaklocation}. $r_0$ labels are specified for 500~nm wavelength and $\delta$ = 0~deg.}
    \label{fig:rmspeak}
\end{figure}

It is instructive to analyze how each of the loss mechanisms in the simulation individually affect the DEPA's performance. Fig. \ref{fig:lossmechanismsgood} plots the contributing power efficiency loss factors for the $r_0$ = 15~cm, $v_\mathrm{g}$ = 10~m/s case as a function of zenith angle. Black items are fundamental losses which cannot be improved. High order aberrations, in blue, may be mitigated by decreasing the sub-aperture size. Piston phase and tip-tilt efficiency, in green, may be increased with further hardware and servomechanism development. Since the subelements are modular, all losses besides piston phase are applicable to any array size.

\begin{figure}
    \centering
    \includegraphics[width=0.7\textwidth]{lossmechanisms15cmr010mpsvg_bigtext}
    \caption{A breakdown of all of the loss mechanisms as a function of zenith angle in the simulations for $r_0$ = 15~cm and $v_\mathrm{g}$ = 10~m/s. Diffraction efficiency, atmospheric transmission, and time-of-flight (black) are all fundamental losses that cannot be avoided. High-order aberrations (blue) can only be mitigated by decreasing the sub-aperture diameter, which requires more subelements for a given array diameter. Piston phase and tip-tilt errors (green) may be reduced with ongoing hardware development. Other effects, such as those due to back reflections of the transmit beam or thermal blooming, are assumed to be sufficiently mitigated or negligible. Error bars represent the standard deviation over the 500~ms time series.}
\label{fig:lossmechanismsgood}
\end{figure}

\section{Scaling to larger arrays}

With use of identical, modular subelements, it is expected that piston phase losses are the only thing that will change as a function of array size. Maintaining good phase integrity may be difficult for larger systems, as the closed loop bandwidth of SPGD is inversely proportional $N$. Long fiber baselines will also add phase noise \cite{Srinivasan2018_phasenoise}. However, since the proposed DEPA system can be built with nested splitting/combining junctions between the master oscillator and the sub-apertures, it is compatible with a nested SPGD configuration \cite{Goodno2014_nspgd}. Nesting may be accomplished by running SPGD in parallel with a detector at every junction, exponentially increasing the system's closed loop bandwidth for a given acquisition frequency and number of subelements. For example, if the subelements are nested in a 5-deep tree, over one billion can be phased together with only 64 sources at each junction. A demonstration that such a configuration can overcome the expected phase perturbations would suggest that kilometer scale systems placed at high altitude sites are capable of producing diffraction limited spots through the atmosphere.

\section{Conclusion}
Detailed dynamic beam propagation simulations of a directed energy phased array transmitting through atmosphere in conditions common to astronomical sites were performed using realistic system noise derived from laboratory measurements. A recently developed beam propagation code optimized via parallelization for very large numbers of elements (>$10^{10}$) in the array made the simulation of a various scenarios possible. The results provide evidence that large aperture ($>>r_0$) ground based systems at near-infrared wavelength are capable of producing a stable diffraction limited spot in outer space if the atmosphere has $r_0\geq$ 10~cm (at 500~nm wavelength and zenith angle = 0~deg) and zenith angles $\leq$ 60~deg. These conditions are often satisfied at current astronomical quality high altitude ground based sites. With the success of nested SPGD and the option to add an arbitrary number of subelements, such systems would overcome the size and power limits of current coherent directed energy technology and enable the delivery of high flux beams over vast distances. This development would enable transformational capabilities such as powering lunar bases or propelling spacecraft for rapid interplanetary travel, ablation of incoming asteroids, and even direct photon momentum transfer to enable the first relativistic spacecraft thus enabling interstellar flyby robotic missions \cite{Lubin2016_roadmap}. Similar technology may also be used for terrestrial applications such as powering of ground-based remote sensors and rovers in inhospitable locations, remote sensing LIDAR, laser machining, and matter in extreme EOS (equation of state) conditions.

\section*{Funding}
Funding for this program comes from NASA grants NIAC Phase I DEEP-IN – 2015 NNX15AL91G and NASA NIAC Phase II DEIS – 2016 NNX16AL32G and the NASA California Space Grant NASA NNX10AT93H as well as a generous gift from the Emmett and Gladys W. fund.

\section*{Acknowledgments}
PML acknowledges support from the Breakthrough Foundation as part of the Starshot program.

\section*{Disclosures}

The authors declare no conflicts of interest.

\section*{Data availability}

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.\\
\\
See \textcolor{urlblue}{Supplement 1} for supporting content.

\section{Modeling of perturbations and adaptive optics}
\label{sec:perturbations}
Each subelement of the DEPA system is assumed to have a phase shifter to control piston phase (average phase over sub-aperture) and tip-tilt actuators to control output angle. It is also assumed that the DEPA is in a bidirectional master oscillator power amplifier (MOPA) configuration such that the beacon light can enter the system and be measured after combining at a coherent junction (see Fig. \ref{fig:schematic}). The bidirectional optical amplifier will be described in an upcoming paper. In an ideal case with no atmosphere (space environment), the inherent optical reciprocity of this bidirectional system ensures that optimizing coherently combined  beacon power at the location of the master oscillator also optimizes the power of the transmit beam at the source of the beacon. With an atmosphere (earth based array), scintillation and time-of-flight effects break this reciprocity, but the principle is the same.

\subsection{Piston phase}
\label{subsec:pistonphase}
SPGD is routinely used in adaptive optics (AO) systems where measurement of the individual sub-aperture phases is unreasonable. In this method, a beacon or laser guide star provides the wavefront reference. Elements of the AO system are dithered to measure the local wavefront gradient and adjusted actively to compensate for it. The control equation is given by \cite{cbc2013_spgd}:
\begin{equation}
    A_{i,k+1} = g_\mathrm{leakage}A_{i,k} + g\left( \frac{1}{2\sigma_\mathrm{rms}^2}\right) \left( \frac{S_+ - S_-}{S_+ + S_-} \right)D_{i,k}.
\end{equation}
Here $A_{i,k}$ represents the correction phase applied to element $i$ for iteration $k$, $D_{i,k}$ is its dither value, $S_+$ and $S_-$ are measurements of the performance metric for positive and negative dithers, $\sigma_\mathrm{rms}$ is the root mean square amplitude of $D$, and $g$ and $g_\mathrm{leakage}$ are gain terms. Our simulations use a Hadamard matrix for $D$, $\sigma_\mathrm{rms}$ = 0.01~waves ($\pi/50$~radians), $g_\mathrm{leakage}$ = 1, $g$ = 2, and an iteration frequency $f_\mathrm{i}$ = 1~MHz. $S_+$ and $S_-$ are computed by squaring the amplitude of the sum of the phase shifted electric fields collected by each subelement (see subsection \ref{subsec:power}) and represent the combined beacon power at the detector in Fig \ref{fig:schematic}. The phase convergence time $t_\mathrm{c}$ is
\begin{equation}
\label{eq:convergence}
    t_\mathrm{c} = \frac{\alpha N}{f_\mathrm{i}}
\end{equation}
for $N$ elements. $\alpha$ is a system-dependent coefficient and represents how many cycles through the dither matrix it takes for the performance metric to reach a steady state. With an orthonormal dither matrix (such as a Hadamard) it is typical to have $\alpha$$\sim$2 \cite{cbc2013_spgd}.

Atmospheric phase inputs are generated by spatially averaging the simulated beacon phase over the area occupied by each sub-aperture. These phases are added to randomly generated phases based on measured optical phase noise power spectra in fiber amplifiers \cite{Srinivasan2018_phasenoise}. Detector noise is generated with a flat power spectrum, and scaled to have a noise-equivalent power (NEP) 40~dB below the nominal beacon power per subelement in order to provide a sufficient signal-to-noise ratio. For example, if the nominal beacon power per subelement were 1~nW, this corresponds to a NEP of $10^{-13} \ \mathrm{W}/\sqrt{\mathrm{Hz}}$, which is a readily available specification from commercial sources. The residual phase errors of the SPGD simulation are then used as perturbations for the transmit beam.

\subsection{Tip-tilt}
\label{subsec:tiptilt}
The slope of the phase over an aperture will cause the beam to ``wander,'' as a phase perturbation that is linear with position is equivalent to an angular displacement of the beam in the focal plane. If the subelements of the DEPA are large fiber collimators, this angular displacement may be compensated for by moving the fiber in the focal plane of the lens. Under the paraxial approximation these displacements are calculated as:
\begin{equation}
\label{eq:tiptilt}
    \theta\cos{\phi}\boldsymbol{\hat{x}} + \theta\sin{\phi}\boldsymbol{\hat{y}} = \nabla\Phi(x,y)/k
\end{equation}
\begin{wrapfigure}{r}{0.45\linewidth}
\includegraphics[width=1\linewidth]{frequencyrejectionsmall_bigtext}
\caption{Servo rejection of tip-tilt actuators in development at University of California, Santa Barbara. The residual position error can be approximated as the setpoint passed through a first order Butterworth high pass filter with $f_\mathrm{c}$ = 13~Hz. For these simulations, we assume we reach our goal of $f_\mathrm{c}$ = 40~Hz.}
\label{fig:tiptilt}
\end{wrapfigure}
where $\theta$ is the angle from the target trajectory, $\phi$ is the azimuthal angle, $\nabla\Phi(x,y)$ is the gradient of the atmosphere-induced optical phase profile, $k$ is the wave number and $x$ and $y$ are in the horizontal plane (see Fig. \ref{fig:schematic}). The two components of $\theta$ are ``tip'' and ``tilt,'' respectively. Although Eq. (\ref{eq:tiptilt}) is for a point as opposed to a finite aperture, the finite aperture values are approximated by spatially averaging $\nabla\Phi(x,y)$ over the area occupied by each sub-aperture.

The Experimental Cosmology Group at University of California, Santa Barbara has recently developed compact and easily manufacturable actuators for tip-tilt correction \cite{Hettel2020_actuator}. These actuators use light received from the beacon for feedback by modulating the fiber position in a circular perturbation at multiple kHz frequencies and using lock-in amplification to produce a two-dimensional error signal. Experimental data shows that this system roughly acts as a first order Butterworth high pass filter for positional error with a cutoff frequency of $f_\mathrm{c}$ = 13~Hz (see Fig. \ref{fig:tiptilt}). With further design optimization we expect to bring $f_\mathrm{c}$ up to 40~Hz. To model our target actuator performance, the simulated tip and tilt time series are passed through such a filter. A 4.4~kHz tone with an amplitude of 150~nrad is added to simulate the modulation used by the servo. The amplitude of the modulation is about 0.01 of the full angle divergence of a 10~cm diameter optic at 1.064~$\mu$m wavelength, well above what was used in \cite{Hettel2020_actuator}.

\subsection{Received and transmitted power}
\label{subsec:power}
The beacon power received by each subelement $P_\mathrm{beac}$ is calculated as inputs for the SPGD simulation. First, the intensity values from the beacon propagation step ($I_\mathrm{beac}$) are scaled by the atmospheric transmission and summed over the area of the sub-aperture to find the total amount of power incident on it. The atmospheric transmission $T$ is a function of zenith angle $\delta$. For a simple slab model,

\begin{equation}
    T = 1 - \sec{\delta}(1-T_0)
\end{equation}
where $T_0$ is the transmission for $\delta = 0$. We take $T_0 = 0.946$ at 1.064~$\mu$m wavelength based on Modtran simulations for a site at 4~km altitude \cite{Lubin2020_ESA}.

The amount of power coupled into the fiber is then approximated as the power incident on the aperture multiplied by the normalized on-axis intensity in the focal plane of the lens. The normalized on-axis intensity is found by combining losses from tip-tilt error and residual high-order phase aberrations over the sub-aperture (exponential terms in Eq. (\ref{eq:beaconpower})). The residual phase $\Phi_r(x,y)$ is calculated by removing the tip and tilt components of $\Phi(x,y)$:
\begin{eqnarray}
\Phi_r(x,y) &=& \Phi(x,y) - \overline{\nabla\Phi_x} \cdot x - \overline{\nabla\Phi_y} \cdot y
\end{eqnarray}
where $\overline{\nabla\Phi_x}$ and $\overline{\nabla\Phi_y}$ are the $x$ and $y$ components of $\nabla\Phi(x,y)$ spatially averaged over the sub-aperture. The variance of the residual phase, which we will denote as $\sigma_\Phi^2$, is then plugged into the Maréchal approximation \cite{Mahajan1983_marechal,Ross2009_marechal}, and we have
\begin{equation}
    P_\mathrm{beac} = Te^{-2(\theta/\theta_{1/2})^2}e^{-\sigma_\Phi^2}\sum_\mathrm{aperture} I_\mathrm{beac}\Delta x \Delta y
\label{eq:beaconpower}
\end{equation}
where $\theta_{1/2}$ is the half angle divergence of the subelement beam.

The power transmitted by each subelement is scaled by $T$ and otherwise assumed not to change, although the far field intensity reduction due to the high-order phase variance $\sigma_\Phi^2$ is accounted for by defocusing the beam by the appropriate amount. This approximation conserves the amount of transmitted power and is a reasonable assumption, since most of the residual phase variance is in the form of focusing perturbations if a Kolmogorov turbulence spectrum is assumed \cite{Noll1976_zernike}.

\subsection{Time-of-Flight}
\label{subsec:TOF}
The finite speed of light becomes a factor in DEPA performance as the atmosphere changes in the time between the beacon and transmit beams passing through it. This effect is accounted for by advancing each phase screen at height $h$ by the time-of-flight delay $\Delta t_\mathrm{TOF}$ in a slab model:
\begin{equation}
    \Delta t_\mathrm{TOF} = 2h\sec{\delta}/c
\end{equation}
where $c$ is the speed of light and $\delta$ is the angle from zenith. These phase screens are added together and subtracted from $\Phi(x,y)$. The differential is then used to produce residual errors for the transmit beam in the same manner as described in subsections \ref{subsec:pistonphase}-\ref{subsec:power}. This calculation again assumes that the atmosphere is thin compared to the target distance.

\section{Atmospheric modeling and phase screen generation}
\label{sec:atmosphere}
\subsection{Turbulence model}
The Hufnagel-Valley (HV) model is commonly used to describe the atmospheric turbulence structure function $C_n^2$ as a function of altitude \cite{Tyson2004_aoguide}. It includes three terms, which respectively correspond to a turbulence peak in the tropopause layer, an exponential decay through the troposphere, and a strong turbulence patch near the ground:
\begin{equation}
\label{eq:HV}
    C_n^2(h) = 5.94 \times 10^{-23}h^{10}\left(\frac{v_\mathrm{rms}}{27}\right)^2e^{-h} + 2.7 \times 10^{-16}e^{-2h/3} + Ae^{-10h}
\end{equation}
where $v_\mathrm{rms}$ is the root mean square wind velocity in meters per second \cite{Toyoshima2011_HVmodel}, $h$ is altitude above ground level (AGL) in kilometers, and $A$ is a tunable parameter to set the strength of turbulence at the surface. $C_n^2$ has units of $\mathrm{m}^{-2/3}$. A common condition of choice is the HV 5-7 model, where $v_\mathrm{rms}=21$~m/s and $A=1.7 \times 10^{-14}~\mathrm{m}^{-2/3}$ to yield a Fried parameter $r_0$ of 5~cm. The Fried parameter, which sets the size scale of the total integrated turbulence, is related to $C_n^2$ by
\begin{equation}
\label{eq:r0}
    r_0 = \left[0.423k^2\sec{\delta}\int_{h_1}^{h_2} C_n^2(h)dh \right]^{-3/5}
\end{equation}
for a plane wave of wave number $k$ which propagates along a path from $h_1$ to $h_2$ with a zenith angle $\delta$. Since $r_0$ is wavelength dependent, it is typical for astronomers to specify $r_0$ at 500~nm wavelength. Although $r_0$ is also dependent on $\delta$, for simplicity we specify its value at $\delta = 0$ and also specify $\delta$. The tuning parameters $A$ and $v_\mathrm{rms}$ alone unfortunately do not provide enough control over $r_0$ to simulate some low turbulence conditions of interest to us. For example, even if we set $A=0$ and $v_\mathrm{rms} = 21$~m/s we obtain $r_0$ $\sim$11~cm for 500~nm wavelength, which prevents the simulation of some common conditions at high altitude sites (see Fig. \ref{fig:wind&r0overlay}). To provide more tunability, we set $A=0$ and introduce a parameter $B$ to modify the troposphere's exponential decay term:
\begin{equation}
\label{eq:HVmod}
    C_n^2(h) = 5.94 \times 10^{-23}h^{10}\left(\frac{v_\mathrm{rms}}{27}\right)^2e^{-h} + B \times 2.7 \times 10^{-16}e^{-2h/3}.
\end{equation}
This modification is conservative in how it affects the DEPA performance for a given value of $r_0$. If we were to compare two atmospheres described by Eqs. (\ref{eq:HV}) and (\ref{eq:HVmod}) with the same $r_0$ and $A>0$, by Eq. (\ref{eq:r0}) they would have the same total integrated turbulence along the path. Therefore, in Eq. (\ref{eq:HVmod}) the turbulence that is lost from the surface patch is simply redistributed to higher altitudes. The shape of the $C_n^2$ distribution primarily affects the degree at which the beacon is scintillated and the time scale at which the perturbations change. Beacon scintillation increases as turbulence is redistributed to be further away from the DEPA, as this is a result of light propagating after it hits turbulent patches. The time scale of the perturbations is dominated by wind speed, which typically increases with altitude and reaches its maximum at the tropopause layer (see subsection \ref{sec:windmodel}). Both increasing beacon scintillation and the time scale of perturbations will worsen the DEPA performance. Thus, by ignoring ground effects and introducing $B$, we do not ease the performance requirements of the DEPA from those that would be determined by the use of a standard HV model with the same $r_0$.

\subsection{Wind model}
\label{sec:windmodel}
Wind speed as a function of altitude is commonly expressed using a Greenwood wind model, which assumes a Gaussian profile with the peak at the tropopause layer \cite{Mohr2010_windmodel,Tyson2004_aoguide}. The Bufton wind model \cite{Tyson2004_aoguide} is a specific Greenwood model where the optical system is pointed at zenith, the ground wind speed is 5~m/s, and the tropopause layer is 9.4~km AGL (Above Ground Level) with a width parameter of 4.8~km and peak velocity of 30~m/s. This study uses a modified Bufton wind model with a tunable ground wind speed $v_\mathrm{g}$:
\begin{equation}
\label{eq:bufton}
    V(h) = v_\mathrm{g} + 30\mathrm{exp}\left[ -\left(\frac{h-9.4}{4.8}\right)^2\right]
\end{equation}
which allows $v_\mathrm{rms}$ to be computed as
\begin{equation}
    v_\mathrm{rms} = \left(\frac{1}{h_2-h_1}\int_{h_1}^{h_2} V^2(h)dh \right)^{1/2}
\end{equation}
which is then plugged into Eq. (\ref{eq:HVmod}). When this model is used for $\delta > 0$, it assumes a worst case scenario of the wind being perpendicular to the propagation direction. Under these definitions, the structure and dynamics of the atmosphere are completely determined by $r_0$ and $v_\mathrm{g}$.

\subsection{Choice of atmospheric conditions}
Data collected at proposed sites for the Thirty Meter Telescope (TMT) \cite{Schock2009_tmtdata} were used to determine appropriate values of $r_0$ and $v_\mathrm{g}$. These sites include Cerro Tolar, Cerro Armazones, Cerro, Tolonchar, San Pedro Mártir, and Mauna Kea. Probability functions and cumulative distributions were produced from DIMM seeing and sonic anemometer data, shown in Fig. \ref{fig:tmtdata}. $r_0$ values are for 500~nm wavelength. Our chosen conditions to simulate are denoted by the red diamonds on Fig. \ref{fig:wind&r0overlay}. These points, which include $r_0$ = 5-15~cm and $v_\mathrm{g}$ = 5-15~m/s, were picked to cover some of the most probable conditions expected on a clear night at an astronomical site as well as a variety of poor conditions to evaluate the robustness of the system.

\begin{figure}
    \centering
        \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\textwidth]{wind_r0cumulativereduced_bigtext}
        \caption{Cumulative probabilities for ground wind speed ($v_\mathrm{g}$) and Fried parameter ($r_0$) based on TMT site survey data.}
        \label{fig:wind&r0cumulative}
    \end{subfigure}
    ~
    \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\linewidth]{wind_r0overlay_bigtext}
        \caption{Probability density function of $r_0$ and $v_\mathrm{g}$. Red diamonds show simulated conditions, and white contours show cumulative probability $p$ within the data set.}
        \label{fig:wind&r0overlay}
    \end{subfigure}

    \caption{Fig. \ref{fig:wind&r0cumulative} shows the probability of having $r_0$ above or $v_\mathrm{g}$ below the specified value based on site survey data taken for the Thirty Meter Telescope (TMT). $r_0$ values are for 500~nm wavelength. Fig. \ref{fig:wind&r0overlay} shows the probability density for a given combination of $r_0$ and $v_\mathrm{g}$ from the same data. It should be noted that some high wind speed data points were discarded for this plot because $r_0$ was not measured at these times. Red diamonds mark the conditions simulated in this paper. The white contours show the cumulative probability $p$ of achieving a given set of conditions or better, and have a spacing of 0.1. The cumulative probability in the top left corner is 1.}
\label{fig:tmtdata}
\end{figure}

\begin{figure}
    \centering
        \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\textwidth]{windinputmodels_bigtext}
        \caption{Wind speed as a function of altitude above ground level for each of the simulated conditions.}
        \label{fig:windinputs}
    \end{subfigure}
    ~
    \begin{subfigure}[b]{0.45\linewidth}
        \includegraphics[width=\linewidth]{cn2profiles_bigtext}
        \caption{$C_n^2$ as a function of altitude above ground level for each of the simulated conditions. The HV 5-7 model is also shown for reference.}
        \label{fig:cn2inputs}
    \end{subfigure}

    \caption{Wind speed (left) and $C_n^2$ (right) profiles for the simulated atmospheric conditions, based on Eqs. (\ref{eq:HVmod}) and (\ref{eq:bufton}). The HV 5-7 model is shown for reference in Fig. \ref{fig:cn2inputs}. $r_0$ values are for 500~nm wavelength and $\delta$ = 0~deg.}
\label{fig:inputmodels}
\end{figure}

\subsection{Phase screen generation}
Atmospheric phase screens were produced using SOAPY \cite{Reeves2016_soapy}, a Python extension for modeling adaptive optical systems based on AOtools \cite{Townson2019_aotools}. It produces a specified number of two dimensional phase screens that correspond to turbulence patches at different heights based on a Kolmogorov spectrum. SOAPY allows the user to input custom wind and $C_n^2$ profiles, which we computed using Eqs. (\ref{eq:HVmod}) and (\ref{eq:bufton}) and are shown in Fig. \ref{fig:inputmodels}. For our desired fidelity, we generated phase screens with 2~cm horizontal resolution, 526~m vertical resolution, and 100~$\mu$s temporal resolution. The height of each atmosphere was 20~km with 0.5~s duration. Increasing the temporal or spatial resolution did not change the results to the reported accuracy. For angles off of zenith ($\delta > 0$), the distance between screens was scaled by $\sec{\delta}$ and the phase values by $\sqrt{\sec{\delta}}$ to preserve $C_n^2$. Since the SPGD simulation requires higher time resolution, intermediate phase and power inputs from Section \ref{sec:perturbations} were generated via spline interpolation and the phase outputs were downsampled back to the original frequency.