# Atmospheric Newtonian noise modeling for third-generation gravitational wave
detectors

D. Brundu<EMAIL_ADDRESS>I.N.F.N, Sezione di Cagliari, Cittadella
Universitaria, 09042 Monserrato, Italy  M. Cadoni<EMAIL_ADDRESS>I.N.F.N, Sezione di Cagliari, Cittadella Universitaria, 09042 Monserrato,
Italy  Dipartimento di Fisica, Università di Cagliari, Cittadella
Universitaria, 09042 Monserrato, Italy  M. Oi<EMAIL_ADDRESS>I.N.F.N,
Sezione di Cagliari, Cittadella Universitaria, 09042 Monserrato, Italy
Dipartimento di Fisica, Università di Cagliari, Cittadella Universitaria,
09042 Monserrato, Italy  P. Olla<EMAIL_ADDRESS>I.N.F.N, Sezione di
Cagliari, Cittadella Universitaria, 09042 Monserrato, Italy  ISAC-CNR, Sezione
di Cagliari, Cittadella Universitaria, 09042 Monserrato, Italy  A. P. Sanna
<EMAIL_ADDRESS>I.N.F.N, Sezione di Cagliari, Cittadella Universitaria,
09042 Monserrato, Italy  Dipartimento di Fisica, Università di Cagliari,
Cittadella Universitaria, 09042 Monserrato, Italy

###### Abstract

The sensitivity and the frequency bandwidth of third-generation gravitational-
wave (GW) detectors are such that the Newtonian noise (NN) signals produced by
atmospheric turbulence could become relevant. We build models for atmospheric
NN that take into account finite correlation times and inhomogeneity along the
vertical direction, and are therefore accurate enough to represent a reliable
reference tool for evaluating this kind of noise. We compute the NN spectral
density from our models and compare it with the expected sensitivity curve of
the Einstein Telescope (ET) with the xylophone design. The noise signal decays
exponentially for small values of the frequency and the detector’s depth,
followed by a power-law for large values of the parameters. We find that, when
the detector is built at the earth’s surface, the NN contribution in the low-
frequency band is above the ET sensitivity curve for strong wind. Building the
detector underground is sufficient to push the noise signal under the ET
sensitivity curve, but the decrement is close to marginal for strong wind. In
light of the slow decay with depth of the NN, building the detector
underground could be only partially effective as passive noise mitigation.

††preprint: ET-0126A-22

##  I Introduction

The first direct observation of gravitational waves (GW) performed in the last
years by the LIGO–Virgo collaboration [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
represents a milestone for fundamental physics and astrophysics. The direct
detection of the GW signals generated by coalescing objects like neutron-stars
or black-hole binaries has not only provided a striking confirmation of
Einstein’s General Relativity in the strong-field regime, but has also started
the new era of multi-messenger astrophysics. The international network of
second-generation GW detectors has been further enhanced with the joining of
the KAGRA detector in 2020 [11, 12, 13]. The currently operating GW detectors
use extremely sensitive Michelson interferometers and have a sensitivity band
ranging from $10$ Hz to $10$ kHz.

Third-generation GW detectors like the Einstein Telescope (ET) [14] and Cosmic
Explorer (CE) [15] have been proposed to fully open the emerging field of GW
astrophysics and cosmology [16, 17]. Their goal is to improve the sensitivity
by a factor of $10$ and push the observation band down to $1$ Hz. These
improvements are motivated by numerous scientific reasons [17]. Black-hole
mergers could be observed at higher redshift and mass; the inspiral phase
could be detected earlier allowing for a better multi-messenger investigation
of the source; possible quantum gravity effects, e.g. quantum hair(s) for
black holes, could be detected in the ringdown phase. These are just a few
examples of the scientific relevance of third-generation GW detectors.

Improvements in sensitivity and frequency bandwidth pose formidable challenges
due to the impact of various noise sources. This is not only because the
enhancement of the sensitivity of a factor of $10$ may push the latter down
close to the noise floor, but also because we need better modeling of noise in
the lowest part of the frequency band. In fact, at few Hz, the major
limitations are expected to come from the gravitational fluctuations, also
called gravitational gradient noise or Newtonian Noise (NN) [18, 19]. NN has
two main contributions coming from seismic fields and atmospheric
perturbations. While detailed estimations of seismic NN have been
performed–also by modeling seismic sources (for a comprehensive review, see,
e.g. Ref. [18] and references therein)–atmospheric NN is instead poorly
understood. Our understanding of the atmospheric contribution to the NN in GW
detectors has, until now, been essentially based on the works of Saulson [20]
and Creighton [21] (see also Ref. [22] for an analysis of the contributions
coming from pressure fluctuations produced in turbulent flow and Ref. [22] for
an analysis of acoustic NN). Saulson considered the effect of acoustic
pressure waves. Conversely, Creighton investigated the contributions to NN of
temperature perturbations, transient atmospheric shocks and sound waves
generated by colliding objects. This lack of interest in the atmospheric
contributions to the NN is motivated by the fact that estimations predicting
NN in the $10$ Hz region are several orders of magnitudes below the
sensitivity curve of second-generation GW interferometers [18].

The situation changes drastically when one considers third-generation GW
detectors. In this case, the NN prediction derived from modeling wind-advected
temperature fluctuations, which are the dominant sources of atmospheric
density fluctuations, cannot be extended below $10\ \rm Hz$ without modifying
the models [21]. This is because the two basic assumptions adopted in [20, 21,
18] for the estimation of NN generated by temperature fluctuations, namely,
the quasi-static approximation and the homogeneity and isotropy hypothesis,
are expected to fail in that frequency band.

A second crucial issue is that some third-generation GW detectors are planned
to be built underground. Although it is generally qualitatively true that
underground construction of the interferometers will represent passive
mitigation of both seismic and atmospheric NN, quantitative results for the
dependence of the noise level on the detector depth are not presently
available. Again, this is because atmospheric NN modeling in the past was
oriented towards second-generation GW-detectors, i.e. detectors built on the
earth’s surface.

The main purpose of this paper is to improve the modeling of NN generated by
wind-advected temperature fluctuations, to make it a reliable reference tool
for the evaluation of NN for third-generation GW detectors. This is going to
be a necessary ingredient not only for a general preliminary estimation of the
noise, but also for the successive, detailed evaluation of NN by numerical
simulation of atmospheric flows in realistic conditions. This is mainly
because current atmospheric codes have a grid scale that is typically well
above that of the fluctuations expected to contribute to NN.

We will improve the quasi-static and homogeneous models of [21] by working in
two different directions. Firstly, we will go beyond the former approximation
by building models in which the NN generated by the decay of the vortices is
taken into account. Secondly, we will fully take into account the fact that
atmospheric turbulence is strongly inhomogeneous along the vertical direction,
and therefore it cannot be modeled _a priori_ within a homogeneous-isotropic
(HI) turbulence framework. We will find that, whenever vortex time-decay
dominates over the wind-advection component, the NN power spectra are
characterized by a power-law behavior. This result is only weakly dependent on
the specific form of the time correlations of the turbulence, except for its
scaling properties. This is fully expected given the multi-scale behavior of
turbulent phenomena. On the other hand, the power-law regime of NN spectra has
not been previously found in HI models of frozen turbulence [21], and its
possible impact on GW detectors was therefore completely overlooked.

We will then apply the results of our models to assess the impact of
temperature-fluctuations-induced NN on the planned ET detector. For this
purpose, we will compare the numerical NN power spectrum obtained from our
models with the expected sensitivity curve of the ET with the xylophone design
(ET-D configuration). We will also discuss the dependence of the power
spectrum on the physical parameters of our model, putting a particular
emphasis on the dependence on the detector depth $r_{0}$. In the frequency
range of interest, we will find that the noise, as a function of $r_{0}$,
decays exponentially for small values of this parameter. A $1/r_{0}^{2}$
scaling sets in at relatively large values of $r_{0}$. This scaling behavior
turns out to be generally true whenever the effect of wind advection is
negligible and $r_{0}$ is large. This implies that building the detector
underground could not provide sufficient passive mitigation of atmospheric NN.

The structure of the paper is as follows. In Section II we discuss in general
terms the modeling of NN from atmospheric temperature fluctuations beyond the
quasi-static approximation. In Section III we briefly discuss the main
features of turbulence in the planetary boundary layer (PBL). We model the NN
generated by wind-advected HI turbulence in Section IV. In Section V we build
a realistic model for NN generated by turbulence in the PBL. In Section VI we
present the power spectra derived from our models, for selected values of the
parameters, compare them with the ET-D sensitivity curve and discuss our
results. In Section VII we draw our conclusions. We leave technical details to
the appendices.

##  II Atmospheric Noise from Temperature Fluctuation

Density perturbations $\delta\rho(\mathbf{r},t)$ in the air caused by
temperature fluctuations $\delta T(\mathbf{r},t)\equiv\tilde{T}(\mathbf{r},t)$
are a major source of atmospheric NN for GW detectors. The heat in the
atmosphere generates convective turbulence, mixing pockets of cold and warm
air at all length-scales down to the millimeter. Indicating with $\bar{T}$ and
$\bar{\rho}$ the mean temperature and density of air, respectively, from the
ideal gas law at constant temperature we get
$\delta\rho(\mathbf{r},t)=-(\bar{\rho}/\bar{T})\tilde{T}(\mathbf{r},t)$.
Density perturbations generated by temperature fluctuations are typically
several orders of magnitude larger than those generated by pressure
perturbations, which disperse in the atmosphere in the form of infrasound
waves. The gravitational acceleration perturbation
$\delta\mathbf{a}(\mathbf{r}_{0},t)$ produced on the test-mass of the detector
located at $\mathbf{r}_{0}$ at time $t$ is given by

$\delta\mathbf{a}(\mathbf{r}_{0},t)=-\alpha\int
dV\frac{\tilde{T}(\mathbf{r},t)}{|\mathbf{r}-\mathbf{r}_{0}|^{3}}(\mathbf{r}-\mathbf{r}_{0}),$
(1)

where $\alpha=G\bar{\rho}/\bar{T}$ is the conversion factor from temperature
to acceleration fluctuations in the detector. The acceleration fluctuation
given by Eq. 1 must then be projected onto the detector-arm direction in order
to obtain the strain.

Physically, the time variation of the acceleration
$\delta\mathbf{a}(\mathbf{r}_{0},t)$ is the result of two different effects:
the decay of vortices and their transport by the wind with average velocity
$\mathbf{U}$ past the detector. The approach adopted in [21, 18], which
approximates turbulence as a frozen field, is based on the hypothesis that the
vortex decay time-scale (identified with the eddy turnover time) is much
larger than the typical time spent by the vortex in the vicinity of the
detector [21]. In this approximation, the main contribution to Eq. 1 is due to
the effect of frozen temperature fluctuations transported by the wind near the
detector. As explained in the introduction, it is commonly believed that this
approximation breaks down at time-scales larger than $10$ seconds [21, 18].
However, this is a rather intricate point, which deserves to be carefully
analyzed. Indeed, the frozen turbulence limit has to be defined in terms of
both the decay time of temperature correlations in the reference frame of the
wind and the effective time the turbulent structure is perturbing the
detector. We will discuss this issue in detail in Section IV.1.

Direct calculation of the acceleration fluctuation from the temperature field
using Eq. 1 is hopeless. The best we can do is to characterize the gravity
gradient noise with its spectral density (power spectrum)
$S_{g}(\omega,r_{0})$

$\begin{split}S_{g}(\omega,r_{0})=&\int{\rm d}^{3}x{\rm
d}^{3}x^{\prime}\,G_{\bf k}({\bf x},r_{0})G_{\bf k}({\bf
x}^{\prime},r_{0})C^{\text{T}}_{\omega}({\bf x},{\bf x}^{\prime})\end{split}$
(2)

where

$\displaystyle C^{\text{T}}_{\omega}({\bf x},{\bf x}^{\prime})$
$\displaystyle=\int{\rm d}t\,C^{\text{T}}({\bf x},{\bf x}^{\prime};t){\rm
e}^{{\rm i}\omega t},$ $\displaystyle\equiv\int{\rm d}t\,\langle\tilde{T}({\bf
x},t)\tilde{T}({\bf x}^{\prime},0)\rangle{\rm e}^{{\rm i}\omega t},$ (3)

and $G_{\bf k}({\bf x},r_{0})$ is the spatial Green function. In the previous
equations $\langle\tilde{T}({\bf x},t)\tilde{T}({\bf x}^{\prime},0)\rangle$ is
the autocorrelation of the temperature fluctuation field at two different
points and at two different times. The stationarity of the random process
implies this autocorrelation to be a function of the time difference. The
specific form of the Green function $G_{\bf k}({\bf x},r_{0})$, which we
report in Appendix A, depends on both the geometry of the problem and the
coordinate choice. The geometry of the system and the set of coordinates used
to perform the calculations are sketched in Fig. 1.

Figure 1:  Schematic representation of the geometry of the system adopted to
compute the Green functions. The system of cylindrical coordinates
$(x_{\perp},\varphi,x_{3})$ is highlighted, while $\bf U_{x_{3}}$ refers to
the direction of the wind, parallel to the earth’s surface. The dotted blue
arrow is the projection of the wind onto the plane of the test mass, while
$\psi$ is the angle between the detector arm, which is taken along $x_{1}$,
and the wind direction.

##  III Turbulence in the planetary boundary layer

Atmospheric turbulence is concentrated in the lowest portion of the
troposphere, a region called the planetary boundary layer (PBL) [23]. The
structure of the PBL is strongly dependent on the orography, the wind and
weather conditions, and the hour in the day.

Let us consider the structure of the PBL on a typical sunny day—which means
that the layer is unstably stratified—in the presence of strong wind. Under
these conditions, the atmospheric turbulence contribution to NN is expected to
be the maximum. For simplicity we consider the case of a horizontally uniform
PBL developing over a plain region, and assume stationarity over the time-
scales of interest.

The PBL can be subdivided into a surface layer in which mechanical stresses,
originating from the wind interaction with the earth’s surface, dominate the
dynamics, and a convection layer in which stratification is dominant. A sketch
of the PBL structure is shown in Fig. 2.

Figure 2:  Portrait (not on scale) of the PBL structure. Vertical dashed lines
indicate the zeros of the curves for turbulent velocity fluctuations $u$ and
temperature fluctuations $\tilde{T}$.

In the surface layer, the mean wind speed $U=U_{x_{3}}$ is characterized by a
logarithmic vertical profile [24]

$\displaystyle U_{x_{3}}\simeq\frac{u_{*}}{\kappa}\ln(x_{3}/z_{0}),$ (4)

where $\kappa=0.4$ is a universal constant called the von Karman constant,
$u_{*}$ is the so-called friction velocity, which gives the velocity scale of
the profile, and $z_{0}$ is an experimental parameter giving the dependence of
the profile on the roughness of the terrain. In the case of a smooth surface
$z_{0}\simeq 0.1\nu/u_{*}$, where $\nu$ is the viscosity of the fluid. For Eq.
4 to be valid, it is necessary that $x_{3}\gg z_{0}$. Values of the roughness
length $z_{0}$ for different types of terrain are listed in Table 1.

$z_{0}$ (m) | Terrain surface characteristics
---|---
1.0 | city
0.8 | forest
0.2 | bushes
0.05 | farmland (open appearance)
0.008 | mown grass
0.005 | bare soil (smooth)
0.0003 | sand surfaces (smooth)
Table 1:  Schematics of terrain types and $z_{0}$ values [25].

The friction velocity $u_{*}$ gives the amplitude of the turbulent velocity
fluctuations $u$ in the surface layer. A similar parameter $T_{*}$ can be
introduced, giving the amplitude of the temperature fluctuations. Both $u_{*}$
and $T_{*}$ are constant in the surface layer. To characterize the strength of
the wind, we shall use the mean wind speed at a reference height
$x^{\text{ref}}_{3}=10\,{\rm m}$: $U_{\text{ref}}\equiv U_{10{\rm m}}$.

At height $x_{3}$, turbulent eddies are expected to be at most $\sim x_{3}$ in
size, although they tend to be more elongated in the direction of the mean
flow than in the vertical and spanwise directions. Similar considerations are
valid for temperature fluctuations. This identifies a characteristic time-
scale of fluctuations in the surface layer

$\displaystyle\tau(x_{3})\sim x_{3}/u_{*}.$ (5)

We can compare the time-scale $\tau(x_{3})$ with the time-scale of convection
$\tau_{\text{conv}}(x_{3})\sim[\bar{T}x_{3}/(gT_{*})]^{1/2}$, where
$g=9.81\,{\rm m/s^{2}}$ is the gravitational acceleration. The transition from
the surface layer to the convection-dominated part of the PBL takes place for
$\tau({x_{3}})\sim\tau_{\text{conv}}(x_{3})$, which defines (minus) the
Obukhov length [26]

$\displaystyle L_{\text{O}}=\frac{\bar{T}u_{*}^{2}}{\kappa gT_{*}}.$ (6)

For $\bar{T}=300\,{\rm K}$, $T_{*}=1\,{\rm K}$, $U_{\text{ref}}=20\,{\rm m/s}$
and $z_{0}=0.05\,$m, corresponding to $u_{*}\simeq 1.5\,{\rm m/s}$, we get
$L_{\text{O}}\simeq 170\,{\rm m}$.

For $x_{3}>L_{\text{O}}$ the mean velocity and temperature profiles become
almost constant, and, in the presence of unstable stratification, turbulence
consists of thermal plumes whose dynamics is determined by the convection time
$\tau_{\text{conv}}(L_{\text{O}})$.

We can determine the amplitude of the velocity and temperature fluctuations in
the convective layer by requiring that, in stationary conditions, the
turbulent heat flux is independent of the height,
$u_{x_{3}}\tilde{T}_{x_{3}}\sim{\rm constant}$, where $u_{x_{3}}$ and
$\tilde{T}_{x_{3}}$ are the typical velocity and temperature fluctuations at
height $x_{3}$. From $u_{x_{3}}\sim x_{3}/\tau_{\text{conv}}(L_{\text{O}})$,
we then get, for $x_{3}>L_{\text{O}}$,

$\displaystyle
u_{x_{3}}\sim\frac{u_{*}x_{3}}{L_{\text{O}}},\qquad\tilde{T}_{x_{3}}\sim\frac{T_{*}L_{\text{O}}}{x_{3}}.$
(7)

Thermal plumes accelerate as they rise in the convective layer, and
temperature gradients concurrently diminish. The rise of the thermal plumes
stops at the top of the convective layer, at a typical height $x_{3}=1\div
2\,{\rm km}$, in the entrainment layer, where stratification becomes strongly
stable.

###  III.1 Turbulence microstructure

Turbulent structures at scale $L_{\text{O}}$ or higher are expected to
generate a contribution to the noise spectrum at frequencies much below those
of interest for GW detectors. The most significant contribution is likely to
come from fluctuations close to the earth’s surface, and from small turbulent
structures originating—through the Kolmogorov cascade—from larger eddies
higher in the boundary layer. The statistical properties of turbulence can be
quantified in terms of correlation functions for velocity fluctuations $C({\bf
x},{\bf x}^{\prime};t)$ and those for temperature fluctuations $C^{T}({\bf
x},{\bf x}^{\prime};t)$ appearing in Eq. 2. In HI turbulence, the spatial
dependence of $C({\bf x},{\bf x}^{\prime};t)$ $C^{T}({\bf x},{\bf
x}^{\prime};t)$ is on ${\bf x}-{\bf x}^{\prime}$ only, so that we have $C({\bf
x},{\bf x}^{\prime};t)=C^{\text{HI}}({\bf x}-{\bf x}^{\prime},t)$, $C^{T}({\bf
x},{\bf x}^{\prime};t)=C^{\text{HI,T}}({\bf x}-{\bf x}^{\prime},t)$. Hence,
turbulent fluctuations can be described in terms of one-time wavevector
spectra $C^{\text{HI}}_{k}=\int{\rm d}^{3}x\,C^{\text{HI}}({\bf x},0){\rm
e}^{-{\rm i}{\bf k}\cdot{\bf x}}$ and $C^{\text{HI,T}}_{k}=\int{\rm
d}^{3}x\,C^{\text{HI,T}}({\bf x},0){\rm e}^{-{\rm i}{\bf k}\cdot{\bf x}}$. At
sufficiently small scales, in the so-called inertial range, turbulent
fluctuations obey the Kolmogorov scaling [27]

$\displaystyle C^{\text{HI}}_{k}$
$\displaystyle\sim\mathcal{E}^{2/3}k^{-11/3},$ (8) $\displaystyle
C^{\text{HI,T}}_{k}$
$\displaystyle\sim\frac{\mathcal{E}_{T}}{\mathcal{E}^{1/3}}k^{-11/3}.$ (9)

The dimensional constants $\mathcal{E}$ and $\mathcal{E}_{T}$ in Eqs. 8 and 9
give the speed with which turbulent fluctuations at a given scale are
converted to smaller-scale turbulent fluctuations. Viscous forces become
dominant at sufficiently small scales, and this determines the size $\eta_{\rm
K}$ of the smallest eddies (Kolmogorov scale):

$\displaystyle\eta_{\rm K}\sim\mathcal{E}^{-1/4}\nu^{3/4},$ (10)

where $\nu\simeq 0.15\,{\rm cm^{2}/s}$ is the kinematic viscosity of air. If
$L$ is the characteristic size of the largest eddies in the flow (the so-
called integral scale of vortices), and $u_{\rm L}$ and $\tilde{T}_{\rm L}$
are their characteristic velocity and temperature scales respectively,

$\displaystyle\mathcal{E}\sim u_{\rm L}^{3}/L,\qquad\mathcal{E}_{T}\sim u_{\rm
L}\tilde{T}_{\rm L}^{2}/L,$ (11)

and Eqs. 8 and 9 will apply provided that $kL\gg 1$. For $u_{L}\sim
u_{*}=1.5\,{\rm m/s}$ and $0.1\,{\rm m}<L<L_{\text{O}}=170\,{\rm m}$, we would
have $0.02\,{\rm m^{2}/s^{3}}<\mathcal{E}<3.4\,{\rm m^{2}/s^{3}}$,
corresponding to values of the Kolmogorov length $\eta_{\rm K}$ below the
millimeter.

To any given scale $l$ we can associate a characteristic velocity scale
$\Delta_{l}u\sim|{\bf u}({\bf x}+{\bf l})-{\bf u}({\bf
x})|\sim(k^{3}C_{k}^{\text{HI}})^{-1/2}_{k\sim l^{-1}}$. From Eq. 8, then, we
can define an eddy turnover time

$\displaystyle\tau_{k}\sim\left(\frac{l}{\Delta_{l}u}\right)_{l=k^{-1}}\sim\mathcal{E}^{-1/3}k^{-2/3},$
(12)

which gives the typical lifetime of turbulent structures of size $k^{-1}$.

One may expect that at height $x_{3}$ in the PBL, turbulent microstructure at
scale $k^{-1}\ll x_{3}$ could be considered, using HI turbulence concepts, as
originating from eddies of size $L\sim x_{3}$ and velocity scale $u_{x_{3}}$.

Indeed, things are more complicated, at least in the surface layer, because
structures with $kx_{3}\gg 1$ could be part of cascades originating from
structures of size $L\sim x^{\prime}_{3}$ at height $x^{\prime}_{3}>x_{3}$.
The problem is less serious in the convective layer, where thermal plumes
continuously carry smaller eddies upwards and prevent the simultaneous
presence of multiple cascades at any given height. We shall deal with this
issue by representing turbulence as a superposition of contributions from
horizontal layers of different thickness, in which turbulence is approximated
as homogeneous and isotropic.

##  IV Newtonian noise from homogeneous isotropic turbulence

The simplest model of NN production by atmospheric temperature fluctuations we
can devise is realized by treating the PBL as an infinite layer of HI
turbulence. The correlation function $C^{\text{T}}({\bf x},{\bf
x}^{\prime};t)$ has a rather intricate structure in which space and time
dependence are intertwined. Performing a successful analysis of
$C^{\text{T}}({\bf x},{\bf x}^{\prime};t)$ heavily relies on the possibility
of introducing some kind of factorization between its space and time
dependence. The space dependence of the correlation function
$C^{\text{T}}({\bf x},{\bf x}^{\prime};t)$ in Eq. 2 is determined by the
spectrum $C^{\text{HI,T}}_{k}$ in Eq. 9. The time structure of
$C^{\text{T}}({\bf x},{\bf x}^{\prime};t)$ is more complex. Equation 12 tells
us that there is a multiplicity of decay times associated with the different
spatial scales in the problem. The situation is complicated by the so-called
Taylor sweep, which consists of turbulent structures being transported by both
the mean wind and larger turbulent structures, while being stretched and
deformed into smaller-scale turbulent structures.

The analysis of the problem is facilitated by the ordering $U>u$ which allows
us, in the first approximation, to neglect the contribution of larger eddies
to the transport of eddies at any given scale. A similar approach has been
used in the design of some turbulent closures [28]. In this approximation, it
is possible to eliminate the effect of the Taylor sweep by shifting to the
reference frame of the mean wind,

$\displaystyle\langle\tilde{T}({\bf
x},t)\tilde{T}(0,0)\rangle=\langle\tilde{T}({\bf x}-{\bf
U}t,t)\tilde{T}(0,0)\rangle_{\scriptscriptstyle U},$ (13)

where $\langle\tilde{T}({\bf x},t)\tilde{T}(0,0)\rangle_{\scriptscriptstyle
U}$ can be expressed, using Eq. 9, as a superposition of Fourier modes, each
decaying at the time-scale fixed by Eq. 12:

$\displaystyle\langle\tilde{T}({\bf x},t)$
$\displaystyle\tilde{T}(0,0)\rangle_{\scriptscriptstyle
U}=\frac{\mathcal{E}_{T}}{\mathcal{E}^{2/3}}$
$\displaystyle\times\int\frac{{\rm d}^{3}k}{(2\pi)^{3}}\frac{{\rm
d}\omega}{2\pi}k^{-11/3}h(\tau_{k}\omega){\rm e}^{{\rm i}({\bf k}\cdot{\bf
x}-\omega t)}.$ (14)

We assume the function $h$ to be symmetric, normalized to 1, and going to zero
for large values of the argument. Equation (14) tells us that, in Fourier
space, the correlations for temperature fluctuations can be expressed as the
product of $h(\tau_{k}\omega)$ times a function of $k$. We substitute now Eqs.
13 and 14 into Eq. 2, and obtain after simple algebra,

$\displaystyle S_{g}=\frac{\mathcal{E}_{T}}{\mathcal{E}^{2/3}}\int$
$\displaystyle\frac{{\rm d}^{3}k}{(2\pi)^{3}}\,k^{-13/3}|G_{\bf
k}(r_{0})|^{2}$ $\displaystyle\times h\left[\tau_{k}(\omega-{\bf k}\cdot{\bf
U})\right].$ (15)

The Green function $G_{\bf k}(r_{0})$ has been computed in Appendix A. For the
case of HI turbulence, it reads (see Eq. 38)

$\displaystyle|G_{\bf
k}(r_{0})|^{2}=\left(\frac{2\pi\alpha\cos\phi}{k}\right)^{2}{\rm
e}^{-2k_{\perp}r_{0}},$ (16)

where, we recall, ${\bf k}_{\perp}\equiv(k_{1},k_{2},0)$,
$\cos\phi=k_{1}/k_{\perp}$, $r_{0}$ is the depth of the detector, and
$\alpha=G\bar{\rho}/\bar{T}$ is the conversion factor from temperature to
acceleration fluctuations [see Eq. 1]. In Section VI we will use Eq. 15,
together with the Green function (16), as the starting point for the numerical
computation of the spectral density for NN generated by HI turbulence.

We note the small-$k$ divergence of the integrand in Eq. 15, which requires
the function $h$ to decay sufficiently fast for large values of the argument.

The fast decay of the frequency spectrum $h$ reflects the character of the
time decorrelation process as the result of the continuous stretching and
scrambling of turbulent structures by the turbulent flow. This is to be
opposed to the case of thermal fluctuations, where a microscopic mechanism
(molecular motion) is at play, causing the correlation function $C(t)$ to be
not differentiable at $t=0$, and the associated spectrum to decay like
$\omega^{-2}$ at large $\omega$. We can verify that, with such a slow decay of
frequency spectrum, the integral in Eq. 15 would be logarithmic divergent at
small $k$.

Similar difficulties in following a Langevin-equation-based approach in
turbulence have been discussed in [29], concerning subgrid modeling for large-
eddy simulations of wall flows.

###  IV.1 Frozen turbulence limit

An important limit of Eq. 15, which has been explored in [21], is that of
frozen turbulence, in which the time-decay of correlations in the reference
frame of the wind is slow and the time-decorrelation of $C^{\text{HI,T}}({\bf
x};t)$ is only a consequence of the Doppler shift induced by the mean wind.
The frozen turbulence limit is realized by approximating the function $h$ in
Eq. 15 with a Dirac delta, which requires the width of $h$, seen as a function
of $k$, to be much smaller than both $U/\omega$ and the width of $|G_{\bf
k}(r_{0})|^{2}$. Inspection of Eqs. 15 and 16 gives us the condition

$\displaystyle\tau_{\omega/U}\gg\max(r_{0}/U,\omega^{-1}),$ (17)

which means that the structures contributing to the noise must have an eddy
turnover time significantly longer than both the inverse of the frequency and
the transit time over a distance $r_{0}$. The second condition is especially
interesting: $r_{0}/U$ represents the effective time the turbulent structure
is effectively seen by the detector, and Eq. 17 tells us that, for a frozen
turbulence hypothesis to be satisfied, the depth at which the detector is
situated must not be too large. In frozen turbulence conditions, Eq. 15 takes
the form

$\displaystyle S_{g}^{ft}$
$\displaystyle=\frac{\mathcal{E}_{T}}{\mathcal{E}^{1/3}}\int\frac{{\rm
d}^{3}k}{(2\pi)^{2}}\,k^{-11/3}|G_{\bf k}(r_{0})|^{2}\,\delta(\omega-{\bf
k}\cdot{\bf U}).$ (18)

Combining Eqs. 16 and 18 produces the general expression

$\displaystyle S^{ft}_{g}=\hat{S}_{g}^{ft}(\psi,\omega
r_{0}/U)\frac{\alpha^{2}\mathcal{E}_{T}U^{8/3}}{\mathcal{E}^{1/3}\omega^{11/3}},$
(19)

where $\hat{S}_{g}^{ft}$ is dimensionless, and $\psi$ is the angle between
${\bf U}$ and the detector arm (see Fig. 1). The ratio
$\mathcal{E}_{T}/\mathcal{E}^{1/3}$ in the formula coincides with the
parameter $c_{T}^{2}$ in [21], giving $\mathcal{E}_{T}/\mathcal{E}^{1/3}=0.2\
{\rm K}^{2}\ {\rm m}^{-2/3}$.

The integral in Eq. 18 can be computed analytically. The details of the
calculations are described in Appendix B, and the final result is a
combination of hypergeometric functions (see Appendix B).

The small frequency regime of the spectrum can be directly inferred from Eq.
19. As discussed in Appendix C, $\hat{S}_{g}^{ft}$ has as a finite limit
$S_{g}^{ft}(\psi,0)=S_{g}^{ft}(\psi)$, which implies a power-law scaling in
$\omega$ for the dimensional spectrum $S_{g}^{ft}$ (this can also be verified
by taking the limit of the exact expression (B)).

The regime $\omega r_{0}/U\gg 1$ is analyzed in Appendix C and gives us the
exponential behavior

$\displaystyle
S^{ft}_{g}\sim\frac{\mathcal{E}_{T}\alpha^{2}U^{8/3}\cos^{2}\psi}{\mathcal{E}^{1/3}\omega^{11/3}}\left(\frac{U}{\omega
r_{0}}\right)^{\frac{1}{2}}\exp\left(-\frac{2\omega r_{0}}{U}\right).$ (20)

It is quite interesting to compare our Eq. 20 with the results of Ref. [21],
where also the NN generated by HI turbulence in the frozen limit has been
analyzed. We see that the exponential behavior of the noise spectrum as a
function of $\omega$ and $r_{0}$ found in Ref. [21] holds only for frequencies
$\omega$ much higher than $U/r_{0}$.

###  IV.2 Weak wind regime

The weak wind regime is realized for $\omega\tau_{\omega/U}\ll 1$, which
corresponds to approximating $h[\tau_{k}(\omega-{\bf k}\cdot{\bf U})]\simeq
h(\omega\tau_{k})$ in Eq. 15. The analysis in Appendix C produces, in this
case, the power-law behaviors

$\displaystyle
S^{ww}_{g}\sim\frac{\mathcal{E}_{T}\alpha^{2}}{\mathcal{E}^{2/3}}\times\begin{cases}k_{\omega}^{-16/3}r_{0}^{-2},&k_{\omega}r_{0}\gg
1,\\\ k_{\omega}^{-10/3},&k_{\omega}r_{0}\ll 1,\end{cases}$ (21)

where

$\displaystyle k_{\omega}=\mathcal{E}^{-1/2}\omega^{3/2}$ (22)

is the inverse size of vortices with eddy turnover time $\omega^{-1}$.

##  V Newtonian noise from turbulence in the PBL

In this section, we go beyond the isotropic and homogeneous approximation for
turbulence by building a more realistic model, in which the vertical structure
of the PBL is fully taken into account. Conversely, we continue to assume
homogeneous correlations horizontally. Following an approach described in
[30], we model turbulence in the PBL as a superposition of HI turbulence
contributions in horizontal layers $[0,\bar{x}_{3}]$, where $\bar{x}_{3}$
ranges from a minimum height $x_{3,\rm min}\sim z_{0}$—marking the transition
to the region where turbulence is strongly affected by the roughness
geometry—to a maximum height $L_{\scriptscriptstyle PBL}$ identifying the top
of the PBL. As a rule of thumb, one usually sets $x_{3,\rm min}=10z_{0}$ [24].
The latter choice, however, could cut potentially important contributions from
vortices at the lower end of the logarithmic region. In our model, therefore,
we choose a value in-between, ${x_{3,\rm min}}={\rm e}z_{0}$,111The choice is
arbitrary, but it allows a smooth transition, in the vertical wind speed
profile, from a linear behavior in the roughness layer, $U_{x_{3}}\propto
x_{3}/x_{3,\rm min}$, to the logarithmic profile above
$U_{x_{3}}\propto\ln\left({\rm e}x_{3}/x_{3,\rm min}\right)$. to identify the
smallest integral-scale temperature fluctuations contributing to the noise.
These temperature fluctuations have amplitude $T_{*}$ and eddy turnover time
$\tau_{x_{3,\rm min}^{-1}}\sim x_{3,\rm min}/u_{*}$. Thus, a smaller $x_{3,\rm
min}$ will correspond to a stronger high-frequency contribution to the noise.
Indeed, we will show in Section VI that the specific choice of $x_{3,\rm min}$
in terms of $z_{0}$ only alters the high-frequency portion of the spectral
density, which is however always below the sensitivity curve of ET for
reasonable values of the parameters (see bottom panel of Fig. 8).

The contributions to the temperature fluctuation $\tilde{T}({\bf x},t)$ from
layers with $\bar{x}_{3}\geq x_{3}$, $\tilde{T}({\bf x},t|\bar{x}_{3})$, are
assumed uncorrelated. The spectrum in each layer can be assumed to obey
Kolmogorov scaling only for $k\bar{x}_{3}\gg 1$. We must thus extend Eqs. 9
and 12 to integral scales $k\bar{x}_{3}\leq 1$. The simplest possibility is to
assume a sharp transition into the Kolmogorov scaling regime exactly at
$k\bar{x}_{3}=1$, i.e. to set

$\displaystyle C^{\text{HI,T}}_{k}(\bar{x}_{3})$
$\displaystyle\sim\frac{\tilde{T}_{\bar{x}_{3}}\bar{x}_{3}^{3}}{\max\left(1,(k\bar{x}_{3})^{11/3}\right)};$
(23) $\displaystyle\tau_{k}(\bar{x}_{3})$
$\displaystyle\sim\frac{\bar{x}_{3}}{u_{\bar{x}_{3}}\max\left(1,(k\bar{x}_{3})^{2/3}\right)},$
(24)

where Eq. 11 has been used.

For simplicity, the effect of the Taylor sweep in each layer is approximated
with that of a constant wind $U_{\bar{x}_{3}}$. The temperature correlation
resulting from the superposition of the contributions in Eqs. (23) is in the
form

$\displaystyle C^{\text{T}}_{{\bf k}_{\perp},\omega}(x_{3},x_{3}^{\prime})$
$\displaystyle=\int_{\max(x_{3},x^{\prime}_{3})}^{L_{\scriptscriptstyle
PBL}}\frac{{\rm d}\bar{x}_{3}}{\bar{x}_{3}}\int\frac{{\rm
d}^{3}k}{2\pi}C^{\text{HI,T}}_{k}(\bar{x}_{3})$ $\displaystyle\times
h[\tau_{k}(\bar{x}_{3})(\omega-{\bf k}\cdot{\bf U}_{\bar{x}_{3}})]{\rm
e}^{{\rm i}k_{3}\cdot(x_{3}-x_{3}^{\prime})}.$ (25)

As in the case of Eq. 15, the spectrum $h$ is evaluated at the Doppler-shifted
frequency $\omega-{\bf k}\cdot{\bf U}_{\bar{x}_{3}}$, generated by expressing
correlations in the reference frame of the wind in the layer $\bar{x}_{3}$.
Note the factor $\bar{x}_{3}^{-1}$ in the integral, which guarantees that all
layers have equal weight.

We substitute Eq. 25 together with Eqs. 23 and 24 and the propagator (40)
calculated in Appendix A into Eq. 2 and obtain, after straightforward algebra,

$\displaystyle S_{g}$ $\displaystyle=\frac{\alpha^{2}}{\pi}\int_{x_{3,\rm
min}}^{L_{\scriptscriptstyle PBL}}{\rm d}\bar{x}_{3}\int_{0}^{+\infty}{\rm
d}k\int_{0}^{1}{\rm d}p\,\frac{p}{\sqrt{1-p^{2}}}$ $\displaystyle\ \
\times\frac{\bar{x}_{3}^{3}\tilde{T}^{2}_{\bar{x}_{3}}A({\bf
k},\bar{x}_{3})B({\bf k},\bar{x}_{3}){\rm
e}^{-2kpr_{0}}}{U_{\bar{x}_{3}}\max[1,(k\bar{x}_{3})^{13/3}]},$ (26)

where $p=k_{\perp}/k$, and

$\displaystyle A$ $\displaystyle=\int_{0}^{2\pi}{\rm d}\phi\,\cos^{2}\phi\
h[\tau_{k}(\bar{x}_{3})(\omega-{\bf k}\cdot{\bf U}_{\bar{x}_{3}})],$ (27)
$\displaystyle B$ $\displaystyle=k^{2}\int_{0}^{\bar{x}_{3}}{\rm d}x_{3}{\rm
d}x^{\prime}_{3}\,{\rm e}^{{\rm
i}k_{3}(x_{3}-x_{3}^{\prime})-k_{\perp}(x_{3}+x_{3}^{\prime})}$
$\displaystyle=1+{\rm e}^{-2k_{\perp}x_{3}}-2{\rm
e}^{-k_{\perp}x_{3}}\cos(k_{3}\bar{x}_{3}).$ (28)

In Section VI we will use Eq. 26 as a starting point for the numerical
computation of the spectral density for NN generated by turbulence in the PBL.
We note the scaling
$\tilde{T}^{2}_{\bar{x}_{3}}/U_{\bar{x}_{3}}\propto\bar{x}_{3}^{-3}$ for
$\bar{x}_{3}>L_{\text{O}}$ in Eq. 26, which descends from Eq. 7. This suggests
that the dominant contribution to $S_{g}$ comes from the surface layer, where
$\tilde{T}^{2}_{\bar{x}_{3}}/U_{\bar{x}_{3}}=T_{*}^{2}/u_{*}$. Numerical
analysis confirms this indication. This prompts us to replace
$L_{\scriptscriptstyle PBL}\to L_{\text{O}}$ in Eq. 26, thus approximating the
noise production in the PBL with the contribution in the surface layer.

To analyze the contribution to the noise in the surface layer, it is
convenient to shift to the so-called wall units [24], which we identify with a
hat:

$\displaystyle\hat{\omega}=\frac{\omega
z_{0}}{u_{*}},\quad\hat{r}_{0}=\frac{r_{0}}{z_{0}},\quad\hat{L}_{\rm
O}=\frac{L_{\text{O}}}{z_{0}}.$ (29)

Note from Table 1, that for reasonable values of $r_{0}$ and $z_{0}$, we have
always $\hat{r}_{0}\gg 1$. In terms of wall units, Eq. 26 can be written in
the form

$\displaystyle
S_{g}=\frac{\alpha^{2}T_{*}^{2}z_{0}^{3}}{u_{*}}\hat{S}_{g}(\hat{\omega},\hat{r}_{0},\hat{L}_{\rm
O}),$ (30)

with $\hat{S}_{g}$ dimensionless.

###  V.1 Limit behaviors

The asymptotic analysis in Appendix E tells us that, in the two limits of
small and large $\hat{\omega}$, the dominant contribution to $\hat{S}_{g}$
comes from turbulent structures in the integral range ($k\bar{x}_{3}<1$) and
in the inertial range ($k\bar{x}_{3}>1$), respectively. This is associated
with the observation that the large $\omega$ contribution to $S_{g}$ is
increasingly concentrated at smaller $\bar{x}_{3}$.

The mechanism for the separation between small and large frequency behavior
can be understood by considering that the frequency $\omega$ selects, as
dominant contributors to $S_{g}$, turbulent structures with size
$k_{\omega}^{-1}$, such that either the eddy turnover time $\tau_{k}$ or the
transit time $(kU_{\bar{x}_{3}})^{-1}$ is $\sim\omega^{-1}$. A maximum
frequency for the production of integral range vortices of size
$k_{\omega}^{-1}$ is then identified by $k_{\omega}z_{0}\sim 1$, as smaller
integral-scale vortices would have to reside inside the roughness of the
terrain.

The different scaling of the contributions to $S_{g}$ from integral and
inertial range turbulent structures is associated with a transition from a
small to a large frequency behavior that is characterized by an increase in
the decay rate with frequency. Indeed, the asymptotic analysis in Appendix E
gives us, for $\omega^{-1}$ shorter than the eddy turnover time of integral
scale vortices at height $z_{0}$,

$\displaystyle\hat{S}_{g}\sim\hat{\omega}^{-8}\hat{r}_{0}^{-2},\quad\min(\hat{\omega},\hat{r}_{0})\gg
1,$ (31)

which is Eq. 58. This equation tells us that we will always have a
$1/r_{0}^{2}$ behavior of $S_{g}$ at large frequencies and for large values of
$r_{0}$.

For small $\hat{\omega}$, instead, we have two possible behaviors depending on
the magnitude of $\hat{r}_{0}$ (see Eq. 60):

$\displaystyle\hat{S}_{g}\sim\begin{cases}\hat{x}_{3,\rm{max}}^{3},&1\gg\hat{\omega}\hat{r}_{0}\\\
\hat{r}_{0}^{-2}\hat{x}_{3,\rm{max}}^{5},&\hat{\omega}\hat{r}_{0}\gg
1\gg\hat{\omega},\end{cases}$ (32)

where

$\displaystyle\hat{x}_{3,\rm{max}}=\min(\hat{\omega}^{-1},\hat{L}_{\rm O}),$
(33)

and where the ratio $r_{0}/u_{*}$ in $\hat{\omega}\hat{r}_{0}=\omega
r_{0}/u_{*}$ is the eddy turnover time of integral scale vortices at height
$r_{0}$. Thus, at small frequencies, power-law decay in $r_{0}$ of the noise
spectrum will ensue for $r_{0}\sim u_{*}/\omega$, which is the size (and the
height) of integral scale vortices with eddy turnover time $\omega^{-1}$. As
already remarked, for reasonable values of $r_{0}$ and $z_{0}$, we always have
$\hat{r}_{0}\gg 1$, which is why we have chosen to disregard the range
$\hat{\omega}\gg 1\gg\hat{r}_{0}$ in Eq. 31 in the first place. Notice,
moreover, that the $1/r_{0}^{2}$ behavior arises also in the weak-wind regime
of the HI approximation (see Eq. 21). As we show in Appendix F, this is a
general behavior arising whenever the effect of eddy decay dominates over wind
advection and $r_{0}$ is large.

##  VI Results

In this section we numerically compute and analyze the behavior of the power
spectra obtained from Eq. 2 for the three different cases discussed in this
paper: $(a)$ HI turbulence in the frozen limit, $(b)$ HI turbulence with
finite correlation time, and $(c)$ turbulence in a horizontally homogeneous
PBL. Specifically, for case $(a)$, we use the expression (18), while for cases
$(b)$ and $(c)$, as starting point for the numerical computation of the power
spectra, we use the more general Eqs. 15 and 26, respectively. The extended
set of data presented and related to this work are available online [31].
Moreover, we compare the resulting spectra with the ET-D sensitivity curve, to
investigate the possible impact of atmospheric NN on next-generation GW
detector measurements. Note that the ET-D sensitivity curve is expressed in
terms of the strain power spectrum $S_{h}$, which is related to the
acceleration power spectrum $S_{g}$ in Eq. 2 by the relation

$\displaystyle S_{h}=4\frac{S_{g}}{\omega^{4}L_{\text{arm}}^{2}},$ (34)

where $L_{\rm arm}$ is the length of the detector arm, while the factor of $4$
takes into account the number of test-masses in the detector. Indeed, one can
assume that the NN contributions on each test mass are uncorrelated, since
those arising from correlated signals between two or more test masses are
expected to be negligible [21, 18]. We rescale the modeled spectra according
to Eq. 34.

(a) HI turbulence: frozen limit—In the frozen limit, the frequency spectrum of
the temperature fluctuations (function $h$ in Eq. 15), is simply a Dirac delta
(see Eq. 18). This allows us to use the analytic expression for the integral
(18) in terms of a sum of hypergeometric functions provided in Appendix B.

Figure 3:  HI turbulence - frozen limit. Top panel: noise spectra as functions
of the frequency $f=\omega/(2\pi)$ at fixed detector depth $r_{0}$ and with
varying wind speed $U$. Bottom panel: strain spectrum for $f=2\text{ Hz}$ as a
function of the detector depth $r_{0}$ for selected values of the wind speed
$U$. In all cases we have taken $\psi=0$.

The noise curves obtained in this limit, for fixed $r_{0}$ and varying $U$,
are shown in the upper panel of Fig. 3 and compared with the sensitivity curve
of ET-D. We also show the variation of the strain spectrum at the fixed
pivotal value of the frequency of $2$ Hz, as a function of the depth $r_{0}$
(bottom panel of Fig. 3). For all the curves in Fig. 3, we have set the
coefficient $c_{T}^{2}=\mathcal{E}_{T}/\mathcal{E}^{1/3}$ in Eq. 18 equal to
its value in [21], $c_{T}^{2}=0.2\text{ K}^{2}\ \text{m}^{-2/3}$. Here, and in
the following cases (HI turbulence with finite correlation time, and
inhomogeneous turbulence), we have set $\psi=0$, which corresponds to wind
blowing parallel to the detector arm (see subsection (c) for a discussion of
the dependence of the spectra on $\psi$).

As already discussed below Eq. 20, the NN is exponentially suppressed for both
large frequencies and large detector depths. In particular, it is worth noting
that the noise curve is below the ET sensitivity curve for $r_{0}\gtrsim
30\text{ m}$ and for frequencies $f\gtrsim 5\text{ Hz}$, where
$f=\omega/(2\pi)$. Moreover, the larger the wind speed, the larger the noise
amplitude, as explicitly shown in Fig. 3. The bottom panel of Fig. 3 also
shows that, for $f=2\ \rm Hz$, if $r_{0}\lesssim 5\text{ m}$ and $U\geq
10\,{\rm m/s}$, the NN spectrum is comparable or above the sensitivity of
ET-D. The same happens if $r_{0}\lesssim 20\ \rm m$ and $U\simeq 30\ \rm m/\rm
s$. The frequency range considered in the plot is not wide enough to visualize
the power-law behavior predicted by Appendix B in the $\omega\ll U/r_{0}$
limit. We have verified that the latter arises at very small frequencies $f\ll
1\text{ Hz}$, for the values of $U$ and $r_{0}$ considered here.
Alternatively, very large values $U$ or very small values of $r_{0}$ would be
required. If the interferometer is built underground, however, $r_{0}$ is at
least $\sim 1\ {\rm m}$; on the other hand, winds stronger than $30\
\text{m/s}$ are very unlikely.

(b) HI turbulence: finite correlation time—To go beyond the frozen turbulence
approximation, an explicit form of the function $h$ in Eq. 15 must be
selected. However, such a choice is not straightforward. In fact, the
functional dependence of $h$ on the physical parameters relies on the
underlying dynamics of the turbulent structures, whose detailed knowledge is
still missing. A reasonable and likewise simple form for the temporal
correlation spectrum in Eq. 15, which satisfies the conditions given in
Section IV, is a Gaussian function

$\displaystyle h[\tau_{k}(\omega-\mathbf{k}\cdot\mathbf{U})]={\rm
e}^{-\tau_{k}^{2}(\omega-k_{\perp}U\cos(\phi-\psi))^{2}}.$ (35)

We have verified that other choices of $h$, such as that of a top-hat
function, do not produce significantly different results. Using Eq. 35, the
integral (15) has no closed-form solution, thus the noise spectra must be
computed via numerical integration.

The integrand in Eq. 15 as well as in Eq. 26 is concentrated in a tiny section
of the integration domain, whose shape is highly dependent on the parameters
involved. This makes numerical integration by quadratures cumbersome. The
integrals have thus been computed using the VEGAS algorithm [32], which
exploits a Monte Carlo technique with importance sampling and is ideally
suited for multidimensional integrals. It should be pointed out that, with
this approach, the numerical estimation is highly computational consuming
($\sim 3$ core-hour for a single value of $S_{h}$). Therefore, the integrals
have been computed by using the HYDRA framework [33], which is designed to
perform data analysis and numerical integration tasks on massively parallel
platforms, commonly used in the high-energy, particle physics community. Using
this tool, we carried out the integrals in Eqs. 15 and 26 obtaining our
results with a sub-percent precision in a reasonable amount of time.

We have repeated the analysis in Fig. 3, for different values of $U$ and
$r_{0}$, and of the new parameter $\mathcal{E}$, which sets the scale of the
eddy turnover time $\tau_{k}$ (see Eq. 12). We have considered values of
$\mathcal{E}$ in the range $0.01\
\text{m}^{2}/\text{s}^{3}\lesssim\mathcal{E}\lesssim 1\text{
m}^{2}/\text{s}^{3}$ (see discussion following Eq. 11). For consistency with
the frozen turbulence case, we have taken
$\mathcal{E}_{T}=c_{T}^{2}\mathcal{E}^{1/3}$, with $c_{T}^{2}=0.2\,{\rm
K^{2}m^{-2/3}}$. As in the previous subsection, we limited our analysis to the
case $\psi=0$.

Figure 4:  HI turbulence - finite correlation time. Top panel: noise spectra
as a function of the frequency $f$ for $r_{0}=5\text{ m}$ and $\psi=0$ and for
selected values of the parameters in the HI approximation with Gaussian time
correlations. Bottom panel: noise spectra as a function of the detector depth
$r_{0}$ for $f=2\text{ Hz}$ and $\psi=0$ and for selected values of the
parameters in the HI approximation with Gaussian time correlations.

The results are shown in Fig. 4. The curves are quite similar to those in Fig.
3. In general, the noise still decays for large detector depths and large
frequencies, while it increases with the wind speed and with $\mathcal{E}$.
The latter reflects the increase in the amplitude of the temperature
fluctuations with $\mathcal{E}_{T}=c_{T}^{2}\mathcal{E}^{1/3}$, which implies,
through Eq. 12, an increase of the characteristic frequency of the
fluctuations as well. However, a new feature, not present in the frozen
regime, is observed in both panels, namely a transition to a power-law regime
in $f$ and $r_{0}$ if $f$ is very large, $U$ is very small or $r_{0}$ very
large (see Eq. 21). At intermediate values of the parameters (take, for
instance, the curve with $U=5\ \rm m/\rm s,\,\mathcal{E}=0.1\ \rm m^{2}/\rm
s^{3}$ in the top panel of Fig. 4), the noise curves show a transition from
exponential (at low frequencies) to power-law (at large frequencies) behavior.
As expected, the noise curves show a weak dependence on $\mathcal{E}$ as long
as the dominant contribution to the noise comes from wind transport rather
than from vortex decay. This occurs for $U\gtrsim 10\ {\rm m/s}$. This feature
makes the noise curves for HI turbulence with Gaussian time correlation and
those obtained in the frozen approximation comparable in the parameter region
of interest, i.e that in which the noise curves are close to the ET-D
sensitivity curve. Indeed, it can be seen that, in this exponentially-damped
regime, the NN from turbulence with finite correlation time and frozen
turbulence have roughly the same impact on the GW detector. On the contrary,
the frequency at which the transition between the two regimes (exponential
damping and power-law scaling) occurs, strongly depends on the value of
$\mathcal{E}$. In particular, we see that the power-law behavior shifts to
larger frequencies for smaller $\mathcal{E}$.

Figure 5: Inhomogeneous turbulence. Top panel: noise spectra as a function of
the frequency $f$ for $r_{0}=5\text{ m}$ and $\psi=0$ and for selected values
of the parameters in the surface-layer model with Gaussian time correlations.
Bottom panel: noise spectra as a function of the detector depth $r_{0}$ for
$f=2\text{ Hz}$ and $\psi=0$ and for selected values of the parameters in the
surface-layer model with Gaussian time correlations.

Inspection of Fig. 4 also shows that going beyond the frozen approximation
changes only slightly the impact of both the wind speed $U$ and the depth
$r_{0}$ on power spectrum noise curves at small frequencies. The differences
between the noise curves for HI turbulence with Gaussian time correlation and
those for frozen turbulence are quite small. At a fixed wind speed of $10\ \rm
m/\rm s$, for instance, both curves cuts the ET-D sensitivity curve at
$r_{0}\approx 4-5\ \rm m$.

(c) Inhomogeneous turbulence—The noise spectrum is now described by Eq. 26. We
have proceeded as in the (b) case, and assumed Gaussian time correlations, as
described in Eq. 35. The noise spectrum has been evaluated from Eq. 26, as a
function of the experimentally accessible quantities $U_{\rm ref}$, $z_{0}$
and $T_{*}$. The quantity $T_{*}$, which gives the scale of the temperature
fluctuations, enters the expression for $S_{g}$ (26) as a scale factor and has
been fixed at the reference value $T_{*}=1\text{ K}$. For $z_{0}$ we have
taken values corresponding to situations ranging from that of bare soil to
that of a forest or a city district, as described in Table 1. As in the
previous HI turbulence cases, we have set $\psi=0$.

The noise amplitude increases with the wind speed $U_{\rm ref}$ and, for small
frequencies, with the roughness length $z_{0}$. This is not surprising,
considering that the turbulence intensity, parameterized by $u_{*}$, is
proportional to $U_{\rm ref}$ and it increases with $z_{0}$ through Eq. 4. In
the high-frequency region, instead, in correspondence with the onset of the
power-law behavior, the noise curves are characterized by $\sqrt{S_{h}}$
scaling as $z_{0}^{-3/2}$, as predicted by Eqs. 31 and 30. We also recall that
a significant contribution to high frequencies fluctuations is produced in the
region near the ground, where the parameter $x_{3,\rm min}\sim z_{0}$ plays
the role of a cutoff (see Eq. 26); we will return to this point at the end of
the section. We note that, since $T_{*}$ is the same in all curves, the
increase in the noise amplitude with $U_{\rm ref}$ and $z_{0}$ is a
consequence of the increase of the characteristic frequency of the
fluctuations only.

The curves in the top panel of Fig. 5, like those in the top panel of Fig. 4,
are characterized by a transition from an exponential to a power-law behavior
at sufficiently small values of $U_{\rm ref}$ and $z_{0}$. A transition to a
power-law, this time for $r_{0}$, is observed also in the bottom panel of Fig.
5, analogous to the one observed in Fig. 4. The crossover point in $r_{0}$
shifts to the right for large values of $U_{\rm ref}$ and $z_{0}$, which is
consistent with the predictions in the asymptotic large $\omega$ and small
$\omega$ limits provided by Eqs. 31 and 32. This is not surprising, since
$U_{\text{ref}}$ plays in the inhomogeneous model the same role played by $U$
in the HI model. For very small frequencies, moreover, the spectra are
characterized by another power-law scaling, which however can be observed only
for $f\lesssim 2\text{ Hz}$ and for large values of $U_{\text{ref}}$.

Another interesting feature is the scaling with the parameter $z_{0}$,
describing the roughness of the terrain. Indeed, the spectra scale as $\sim
z_{0}^{-3}$ for $\omega\gg u_{\ast}/z_{0}$ while, for $\omega\ll
u_{\ast}/z_{0}$, they first increase with $z_{0}$ and then they become almost
independent of this quantity (see Eqs. 31 and 32). Similarly to what happens
with the parameter $\mathcal{E}$ in the HI case, the weak dependence of the
noise curves on the parameter $z_{0}$ is a characteristic of the regime where
wind transport dominates over vortex-decay. Again, this is a nice feature,
allowing for easy comparison of the curves for inhomogeneous and homogeneous
turbulence in the regions where they are close to the ET-D sensitivity curve.

The last interesting point is the dependence of the spectra on $\psi$. Indeed,
in some regimes (see, e.g., Eq. 20), the spectra could depend strongly on this
parameter. Anyway, we expect $S_{h}$ to be the maximum for $\psi=0$. Indeed,
in Fig. 7, we show that, for reasonable values of the parameters, $S_{h}$
varies at most by a factor $\sim 2$.

On a semi-quantitative level, the results of the inhomogeneous model confirm
those obtained in the HI one (both in the general case and in the frozen
approximation limit). Indeed, for fixed values of the parameters (either $U$
or $r_{0}$), the noise spectra in the three cases have comparable orders of
magnitude, at least in the regions close to the sensitivity curve of ET-D.
When the detector is located near the earth’s surface, the power spectrum
curves cut the ET-D sensitivity and are well above it in the frequency region
$2-10\ \rm Hz$. On the other hand, in the frequency band considered here, the
noise curves go below the sensitivity curve only when the detector is located
at least $50\ \rm m$ underground (see top panel of Fig. 6). We see, however,
that for large values of $U_{\rm ref}$, i.e. $U_{\rm ref}\sim 30\ \rm m/s$,
the noise curve is only a factor $\sim 5$ below the ET-D sensitivity curve,
which is a worryingly close range, considering the fact that our models are
providing only order of magnitude estimates. The situation is confirmed if we
go to greater depths. Consistently with the general scaling $r_{0}^{-1}$ of
$\sqrt{S_{h}}$ derived in Appendix F, the noise spectra are only a factor of
$\sim 10-20$ below the sensitivity curve, for $r_{0}\sim 100-200\ \rm m$,
respectively, (see bottom panel of Fig. 6).

Figure 6: Top panel: Noise spectra for the inhomogeneous model as a function
of the frequency $f$ for $r_{0}=50\text{ m}$ and for selected values of
$U_{\text{ref}}$. Bottom panel: Noise spectra for the inhomogeneous model as
functions of the frequency for selected values of $r_{0}$ and
$U_{\text{ref}}=30\text{ m/s}$.
We fixed the values of the other parameters to: $z_{0}=0.1\text{ m}$ and
$\psi=0$.

Figure 7:  Noise spectra as a function of $\psi$ in the three cases of HI
turbulence - frozen case $(a)$; HI turbulence - finite correlation time $(b)$
and inhomogeneous turbulence $(c)$, at varying $U$ (upper two panels) and
$U_{\rm ref}$ (lowest panel).
We fixed the values of the other parameters to: $f=2\,{\rm Hz}$,
$r_{0}=5\,{\rm m}$, $\mathcal{E}=0.1\,{\rm m^{2}/s^{3}}$, $z_{0}=0.1\,{\rm
m}$.

Wavevector cutoffs and finite-size effects—For all the calculations in the
present section, we have adopted a maximum wavevector $k_{\rm max}$ equal to
the inverse of the Kolmogorov scale $\eta_{\rm K}^{-1}$ defined in Eq. 10, and
we have set the parameter $x_{3,\text{min}}={\rm e}z_{0}$ as the lower bound
of integration in Eq. 26. In the inhomogeneous turbulence case of Eq. 26, we
have carried out the integral over $\bar{x}_{3}$ up to a reference height
$L_{{\scriptscriptstyle PBL}}=2\ {\rm km}$, but the contribution at
$\bar{x}_{3}>L_{\rm O}$ turned out to be negligible.

The dependence of $S_{h}$ on the parameters $k_{\rm max}$ and
$x_{3,\text{min}}$ is rather different. In the top panel of Fig. 8, we show
the effect of lowering the cutoff $k_{\rm max}$ on the noise spectrum. We
limit our analysis to the inhomogeneous turbulence case,

Figure 8:  Inhomogeneous turbulence. Top panel: dependence of the noise
spectrum on the wavevector cutoff $k_{{\rm max}}$. Bottom panel: dependence on
the minimum height $x_{3,\rm min}$ of the vortices contributing to the noise.
Values of the parameters in the two cases: $r_{0}=5\text{ m}$,
$z_{0}=0.1\text{ m}$, $U_{\text{ref}}=10\text{ m/s}$ and $\psi=0$.

setting, as usual, $r_{0}=5\,{\rm m}$ as a reference depth at which the noise
is expected to be above the sensitivity threshold for the ET detector. We see
that, when $k_{\rm max}^{-1}$ approaches the meter range, the high-frequency
portion of the noise spectrum is increasingly damped. Inspection of the top
panel of Fig. 8 indicates that noise-damping occurs above the kink, signaling
the transition from an exponential behavior dominated by wind transport to a
power-law behavior dominated by eddy decay. However, for sufficiently small
$k_{\rm max}$, i.e. $2\pi k_{\rm max}^{-1}\sim 1-2\text{ m}$, this transition
does not occur in the frequency-band considered here, and only exponentially-
damped “frozen-like” contributions are present. This is because vortices of
this size have a very long decay time (see Eq. 12) compared to the advection
time of the wind.

A similar situation occurs, as illustrated in the bottom panel of Fig. 8,
varying $x_{3,\text{min}}$. In this case, the high-frequency portion of the
noise spectrum is increasingly damped as $x_{3,\text{min}}$ gets larger. The
result is consistent with the observation in Section V.1 that the contribution
to the high-frequency portion of the NN spectrum is produced primarily by
temperature fluctuations near the ground.

The present analysis tells us that numerical simulation of the atmospheric
flow, say, by large-eddy simulation [34], would allow direct evaluation of the
NN only for frequencies below a maximum that is an increasing function of
$k_{\rm max}$ and $x_{3,\text{min}}^{-1}$. Frequencies above this maximum
would require some kind of parameterization.

##  VII Conclusions

In this paper, we have built models for the NN generated by atmospheric
turbulence, which represent a reliable tool for the assessment of the impact
of atmospheric noise on third-generation GW detectors. This has been done by
improving previous models for temperature-fluctuation induced atmospheric NN.
Owing to their basic assumptions (frozen in time, HI turbulence) the latter
are not reliable enough in the frequency bandwidth and for the sensitivity
levels of third generation GW detectors. This improvement has been made
possible by building both models for HI turbulence with finite correlation
time for temperature fluctuations and models which also take into account the
strong inhomogeneity of turbulence along the vertical direction.

We have also computed the spectral density of the NN for the three classes of
models as a function of the characteristic physical parameters, and compared
it with the sensitivity curve of the ET detector in the xylophone
configuration. The NN signal shows two kinds of regimes, one in which it
behaves exponentially and the other in which it shows a power-law behavior.
The first appears to be the signature of the dominance of wind transport and
the exponential dependence on the depth of the detector, of the Green function
connecting temperature fluctuations and NN. Since in a wind-dominated regime,
this dependence is weighed by a characteristic scale which is the ratio of the
wind velocity and the frequency, the result is an exponential scaling of the
NN with both the depth of the detector and the frequency. Whenever we depart
from this regime (either because eddy decay becomes an important factor, or
because the depth of the detector is large), the NN dependence on the wind
velocity, frequency and detector depth is a power law. While the power spectra
scale with the velocity and the frequency with a power-law which is model
dependent, the exponent characterizing the $r_{0}$ scaling is fixed at $-2$ by
the properties of the gravitational propagator. Departures from exponential
behavior occur for high frequencies, or very small frequencies, the latter
being out of the regime of interest for ET anyway. Moreover, in the regions of
interest for ET, the NN power spectra show a very weak dependence on the
parameters of the models, apart from the wind speed $U$ and the depth $r_{0}$.
Taken together, these two facts explain why the models of HI frozen turbulence
used to date worked well in the assessment of the impact of atmospheric NN for
second-generation GW detectors [21, 18].

We have then compared the NN power spectrum calculated with our models with
the sensitivity curve of ET. We have found that the atmospheric NN
contribution is above the sensitivity curve in the low-frequency band when the
detector is placed on the surface and/or when the wind speed is relatively
large. Thus, our main result is that NN generated by atmospheric turbulence
represents an important source of noise for third-generation GW detectors,
which therefore must be taken into consideration and accurately analyzed, e.g.
using numerical simulations. This is of paramount importance, especially if
the detector has to be built on the earth’s surface.

On the other hand, we have also found that passive mitigation of atmospheric
NN is only partially effective. Although placing the detector underground
suppresses the atmospheric NN contribution, the rather weak, $1/r_{0}$, decay
of the noise signal implies that even an underground construction at great
depth may not be enough to suppress the noise completely. Even at $r_{0}\sim
200\,{\rm m}$, the noise amplitude curve, despite being always below the
sensitivity curve of ET, remains close to it in the low-frequency region. On
the other hand, our modeling of turbulence provides only order of magnitude
predictions. An inaccuracy of an order of magnitude in our estimates of the
power spectrum should therefore be taken into account. Considering this fact
and the proximity of the noise amplitude to the sensitivity of the detector,
at least in strong wind situations ($\sim 30\text{ m/s}$) (see Fig. 6),
numerical simulations of the atmospheric flow (taking into account the
orography in the detector region) and on-site measurements are advisable.

In this regard, we expect our results to be relevant for any evaluation of NN
by numerical simulation of the atmospheric flow, as atmospheric codes (such
as, e.g., large-eddy simulations [34]) have a grid scale that is typically
above that of fluctuations contributing to the NN.

## Acknowledgements

We thank the Newtonian noise research group of the SARGRAV project and in
particular Jan Harms, for helpful discussions and comments.

##  Appendix A Green function evaluation

To compute the Fourier transform of the gravitational propagator, we consider
the detector at depth $r_{0}$ below the earth’s surface, which will be
considered as a flat infinite plane. We first choose a Cartesian system of
coordinates $x_{1}$, $x_{2}$ and $x_{3}$, with origin at the test mass of the
detector, and with $x_{1}$ laying along the detector arm, while $x_{3}$ will
be along the vertical.

Computations can be more easily performed by adopting the cylindrical system
of coordinates ${\bf x}=\left(x_{\perp},\,\varphi,\,x_{3}\right)$, where
$\varphi$ is the azimuthal angle, while
$x_{\perp}=(x_{1}^{2}+x_{2}^{2})^{1/2}$ refers to the direction perpendicular
to $x_{3}$ (see Fig. 1). Moreover, we will make use of the following
expressions

$\displaystyle G({\bf x};r_{0})=-\alpha\ \partial_{x_{1}}H({\bf x};r_{0});$
(36a) $\displaystyle H({\bf
x};r_{0})=\frac{\theta(x_{3}-r_{0})}{\left(x_{\perp}^{2}+x_{3}^{2}\right)^{1/2}}.$
(36b)

The Fourier transform of Eq. 36b then reads

$\begin{split}H_{\bf k}(r_{0})=&\int_{r_{0}}^{\infty}{\rm
d}x_{3}\int_{0}^{+\infty}{\rm d}x_{\perp}x_{\perp}\int_{0}^{2\pi}{\rm
d}\varphi\\\ &\times\frac{\exp\left[-{\rm
i}\left(k_{3}x_{3}+k_{\perp}x_{\perp}\cos\varphi\right)\right]}{\left(x_{\perp}^{2}+x_{3}^{2}\right)^{1/2}}\end{split}$
(37)

where ${\bf k}=({\bf k}_{\perp},k_{3})$ and ${\bf
k}_{\perp}\equiv(k_{1},k_{2})$.

When dealing with HI turbulence, correlations have to be computed in the whole
three-dimensional Fourier space. Equations 37 and 36a together yield

$\displaystyle G_{\bf k}(r_{0})=-\frac{2\pi\alpha\cos\phi}{k_{\perp}+{\rm
i}k_{3}}{\rm e}^{-(k_{\perp}+{\rm i}k_{3})r_{0}}.$ (38)

where we have defined $\cos\phi=k_{1}/k_{\perp}$.

In more realistic cases, things have to be treated more carefully, as we might
have inhomogeneities along one or more axes. In inhomogeneous turbulence, for
instance, we have inhomogeneities along the $x_{3}$-axis due to the
inhomogeneous wind profile. Since correlations on the plane orthogonal to
$x_{3}$ are assumed to be homogeneous in our model (see Section V), the Green
function will be simply given by the $x_{\perp}$ and $\varphi$ integrals in
Eq. 37, which yield

$\displaystyle G_{\bf k_{\perp}}(x_{3},r_{0})=-2\pi{\rm
i}\alpha\theta(x_{3}-r_{0})\cos\phi\ {\rm e}^{-k_{\perp}x_{3}}$ (39)

Finally, since the effects of turbulence along the vertical will be integrated
from the earth’s surface up to infinity, we simply have to translate the
origin of the system of coordinates along the $x_{3}$-axis by $r_{0}$. In
other words, we simply shift $x_{3}\to x_{3}+r_{0}$, which yields

$\displaystyle G_{k_{\perp}}(x_{3},r_{0})=-2\pi{\rm i}\alpha\cos\phi\ {\rm
e}^{-k_{\perp}x_{3}}{\rm e}^{-k_{\perp}r_{0}}.$ (40)

##  Appendix B Analytic form of the spectrum for homogeneous isotropic NN in
the frozen approximation

In this appendix, we derive the analytic results for the spectrum (18) for HI
turbulence in the frozen approximation. We start from Eq. 18

$\displaystyle
S_{g}^{ft}=\frac{\mathcal{E}_{T}}{\mathcal{E}^{1/3}}\int\frac{{\rm
d}^{3}\mathbf{k}}{(2\pi)^{2}}k^{-11/3}|G_{\mathbf{k}}(r_{0})|^{2}\delta(\omega-\mathbf{k}\cdot\mathbf{U}).$
(41)

To compute the integral, first, we choose the geometry illustrated in Fig. 1,
and we adopt a cartesian set of coordinates in the integrated variable
$\mathbf{k}$. Moreover, we choose the wind speed to be parallel to the earth’s
surface. This simplifies the calculations since

$\displaystyle\mathbf{k}\cdot\mathbf{U}=k_{\perp}U\cos(\phi-\psi),$ (42)

where $U$ is the wind speed and $k_{\perp}=\sqrt{k_{1}^{2}+k_{2}^{2}}$. To get
rid of the angle $\psi$ in the Dirac delta, we rotate the reference frame
around the $k_{3}$-axis by an angle $\psi$ and we write
$\cos\phi=k_{1}/k_{\perp}$. Then, using Eq. 38 we get

$\displaystyle S_{g}^{ft}=$ $\displaystyle\
\frac{\mathcal{E}_{T}}{(2\pi)^{2}\mathcal{E}^{1/3}}\int_{0}^{\infty}{\rm
d}^{3}\mathbf{k}\,k^{-11/3}|G_{\mathbf{k}}|^{2}\delta(\omega-k_{1}U)$
$\displaystyle=$ $\displaystyle\
\alpha^{2}\frac{\mathcal{E}_{T}}{\mathcal{E}^{1/3}}\int_{-\infty}^{\infty}{\rm
d}k_{1}\int_{-\infty}^{\infty}{\rm d}k_{2}\int_{-\infty}^{\infty}{\rm d}k_{3}$
$\displaystyle\times
k^{-17/3}k_{\perp}^{-2}(k_{1}^{2}\cos^{2}\psi+k_{2}^{2}\sin^{2}\psi)$
$\displaystyle\times\delta(\omega-k_{1}U){\rm e}^{-2k_{\perp}r_{0}}.$ (43)

The integral over $k_{3}$ can be easily done and yields

$\displaystyle\int_{-\infty}^{\infty}\frac{{\rm d}k_{3}}{k^{17/3}}$
$\displaystyle=\sqrt{\pi}\
\frac{\Gamma\left(7/3\right)}{\Gamma\left(17/6\right)}k_{\perp}^{-14/3}.$ (44)

The integral over $k_{1}$ is instead trivial using the delta function, which
simply sets $k_{1}=\omega/U$.

The last integration in $k_{2}$ can be solved changing the integration
variable, mapping $k_{2}\to\xi=\sqrt{1+(Uk_{2}/\omega)^{2}}$. With this
substitution, the integral becomes

$S_{g}^{ft}=2\mathcal{C}\int_{1}^{\infty}{\rm
d}\xi\,\frac{\xi}{\sqrt{\xi^{2}-1}}\\\ \
\times\frac{\cos^{2}\psi+(\xi^{2}-1)\sin^{2}\psi}{\xi^{20/3}}{\rm e}^{-2\omega
r_{0}\xi/U},$ (45)

where $\mathcal{C}=\alpha^{2}\mathcal{E}_{T}\sqrt{\pi}\
\Gamma(7/3)/(\mathcal{E}^{1/3}\Gamma(17/6))\ U^{8/3}\omega^{-11/3}$. The
integral can be done analytically and gives a combination of hypergeometric
functions

$\displaystyle S_{g}^{ft}=$ $\displaystyle\
2\mathcal{C}\cos^{2}\psi\,\biggl{\\{}\frac{\sqrt{\pi}\Gamma\left(\frac{17}{6}\right)\,_{1}F_{2}\left(-\frac{7}{3};-\frac{11}{6},\frac{1}{2};x^{2}\right)}{2\Gamma\left(\frac{10}{3}\right)}+32x\Gamma\left(-\frac{17}{3}\right)$
$\displaystyle\times\biggl{[}2^{2/3}x^{14/3}\,_{1}F_{2}\left(\frac{1}{2};\frac{10}{3},\frac{23}{6};x^{2}\right)-\frac{26180\sqrt{\pi}\,_{1}F_{2}\left(-\frac{11}{6};-\frac{4}{3},\frac{3}{2};x^{2}\right)}{6561\
\Gamma\left(\frac{17}{6}\right)}\biggr{]}\biggr{\\}}$
$\displaystyle+2\mathcal{C}\sin^{2}\psi\,\biggl{\\{}\frac{\sqrt{\pi}\Gamma\left(\frac{11}{6}\right)\,_{1}F_{2}\left(-\frac{7}{3};-\frac{5}{6},\frac{1}{2};x^{2}\right)}{4\Gamma\left(\frac{10}{3}\right)}+8x\Gamma\left(-\frac{11}{3}\right)$
$\displaystyle\times\biggl{[}2^{2/3}x^{8/3}\,_{1}F_{2}\left(-\frac{1}{2};\frac{7}{3},\frac{17}{6};x^{2}\right)-\frac{55\sqrt{\pi}\,_{1}F_{2}\left(-\frac{11}{6};-\frac{1}{3},\frac{3}{2};x^{2}\right)}{243\
\Gamma\left(\frac{17}{6}\right)}\biggr{]}\biggr{\\}},$ (46)

where, for convenience, we defined $x\equiv r_{0}\omega/U$.

##  Appendix C Asymptotic results in homogeneous isotropic turbulence

###  C.1 Frozen turbulence

An explicit expression for the dimensionless spectrum $\hat{S}_{g}^{ft}$ in
Eq. 19 is obtained by substituting Eq. 38 into Eq. 18. Let us set ${\bf
U}=(U,0,0)$ and introduce dimensionless quantities $\hat{\omega}=\omega
r_{0}/U$ and $\hat{k}=kU/\omega$. We find

$\displaystyle\hat{S}^{ft}_{g}=\int\frac{{\rm d}\hat{k}_{2}{\rm
d}\hat{k}_{3}}{\hat{k}^{17/3}}\cos^{2}(\psi+\phi){\rm
e}^{-2\hat{k}_{\perp}\hat{\omega}}$ (47)

where $\hat{k}_{1}=1$ from the Dirac delta in Eq. 18. We easily verify that
the integral converges to a finite function of $\psi$ for $\hat{\omega}\to 0$.
This implies that for $\omega\ll U/r_{0}$, $S_{g}$ scales like a power in
$\omega$ and $U$ (see Eq. 19).

In the opposite limit $r_{0}\gg U/\omega$, we can approximate

$\displaystyle\hat{S}_{g}^{ft}$ $\displaystyle\simeq\int\frac{{\rm
d}\hat{k}_{2}{\rm
d}\hat{k}_{3}}{\hat{k}^{17/3}}\left[\frac{\cos^{2}\psi+\hat{k}_{2}^{2}\sin^{2}\psi}{1+\hat{k}_{2}^{2}}\right]$
$\displaystyle\times\exp\left[-\hat{\omega}(2+\hat{k}_{2}^{2})\right],$ (48)

where we have exploited $\hat{k}_{2}\sim\hat{\omega}^{-1/2}\ll 1$. Evaluating
the integral in Eq. 48 to lowest order in $\hat{\omega}^{-1}$ yields Eq. 20.

###  C.2 Weak wind regime

Let us indicate $p=k_{\perp}/k$. We can evaluate the integral in Eq. 15 in
spherical coordinates. If $h[\tau_{k}(\omega-{\bf k}\cdot{\bf U})]\simeq
h(\tau_{k}\omega)$, the integral over $\phi$ in Eq. 15 is trivial, so that we
are left with

$\displaystyle S^{ww}_{g}$
$\displaystyle\sim\frac{\mathcal{E}_{T}\alpha^{2}}{\mathcal{E}^{1/3}}\int_{0}^{+\infty}{\rm
d}k\,k^{-13/3}h(\tau_{k}\omega)\int_{0}^{1}{\rm d}p\,p\,{\rm e}^{-2pkr_{0}}$
$\displaystyle=\frac{\mathcal{E}_{T}\alpha^{2}}{\mathcal{E}^{1/3}}\int_{0}^{+\infty}{\rm
d}k\,\frac{1-(1+2kr_{0}){\rm
e}^{-2kr_{0}}}{4k^{19/3}r_{0}^{2}}h(\tau_{k}\omega)$
$\displaystyle\sim\frac{\mathcal{E}_{T}\alpha^{2}}{\mathcal{E}^{1/3}}\int_{k_{\omega}}^{+\infty}{\rm
d}k\,\frac{1-(1+2kr_{0}){\rm e}^{-2kr_{0}}}{k^{19/3}r_{0}^{2}},$ (49)

where $k_{\omega}$ is defined in Eq. 22. We can now carry out the integral in
Eq. 49 in the two limits of large and small depth and we recover Eq. 21.

##  Appendix D Time correlations for turbulence in the surface layer

Let us indicate

$\displaystyle\eta=\omega\tau_{k}(\bar{x}_{3}),\quad\zeta=k_{\perp}\tau_{k}(\bar{x}_{3})U_{\bar{x}_{3}}.$
(50)

We find the following limit behaviors for the function $A$ in Eq. 27:

  * •

If either $\zeta\ll\eta\sim 1$ or $\eta\gg 1$ and
$\zeta\ll|h(\eta-\zeta)/h^{\prime}(\eta-\zeta)|$,

$\displaystyle A(\eta,\zeta)\simeq\pi h(\eta).$ (51)

  * •

If either $\zeta\gg\max(\eta,1)$ or $\zeta>\eta\gg 1$,

$\displaystyle
A(\eta,\zeta)\sim\frac{\cos^{2}(\psi+\phi_{m})}{\zeta\sin\phi_{m}},$ (52)

where $\phi_{m}=\arccos(\eta/\zeta)$.

  * •

If $\eta\gg\max(1,\zeta)$ and
$\zeta\gg|h(\eta-\zeta)/h^{\prime}(\eta-\zeta)|$,

$\displaystyle A(\eta,\zeta)\sim\cos^{2}\psi\frac{h^{3/2}(\eta-\zeta)}{|\zeta
h^{\prime}(\eta-\zeta)|^{1/2}}.$ (53)

The result in Eq. 51 is straightforward. The result in Eq. 52 is obtained by
saddle-point approximation, expanding $h(\eta-\zeta\cos\phi)\simeq
h(\zeta\phi^{\prime}\sin\phi_{m})$, $\phi^{\prime}=\phi-\phi_{m}$, where
$\phi_{m}=\arccos(\eta/\zeta)$, and then integrating from
$-(\zeta\sin\phi_{m})^{-1}$ to $(\zeta\sin\phi_{m})^{-1}$. The result in Eq.
53 is obtained by expanding

$\displaystyle h(\eta-\zeta\cos\phi)$ $\displaystyle\simeq
h(\eta-\zeta(1-\phi^{2}/2))$ $\displaystyle\simeq
h(\eta-\zeta)+(\zeta\phi^{2}/2)h^{\prime}(\eta-\zeta)$

and then integrating by steepest descent. In order to have $\hat{A}=O(1)$, it
is sufficient that $\eta\sim\zeta\sim 1$.

##  Appendix E Asymptotic results for turbulence in the surface layer

In the integral in Eq. 26 we separate contributions from integral scale
vortices (“domain 1”, $k\bar{x}_{3}<1$) and those from inertial scale vortices
(“domain 2”, $k\bar{x}_{3}>1$). Different cutoffs act in the integral:

  * •

The term ${\rm e}^{-2kpr_{0}}$, accounting for the decay of the signal with
the depth of the detector;

  * •

The factor $k^{-13/3}$, associated with the decay of turbulent fluctuations at
small scales;

  * •

The function $A$ that filters eddies at time-scales below $\omega^{-1}$.

The analysis in Appendix D tells us that the function $A$ is surely negligible
for $\eta\gg\max(1,\zeta)$. Let us analyze the remaining regions $\eta<1$ and
$\zeta>\eta>1$.

The region $\zeta>\eta>1$ corresponds in domain 1 to

$\displaystyle\hat{\omega}\hat{x}_{3}$
$\displaystyle<\frac{\hat{k}_{\perp}\hat{x}_{3}}{\kappa}\ln({\rm
e}\hat{x}_{3})<\frac{1}{\kappa}\ln({\rm e}\hat{x}_{3})$
$\displaystyle\Rightarrow
1>\hat{k}\hat{x}_{3}>\frac{\kappa\hat{\omega}x_{3}}{\ln({\rm e}\hat{x}_{3})}$
$\displaystyle\Rightarrow\hat{\omega}<\frac{\ln({\rm
e}\hat{x}_{3})}{\kappa\hat{x}_{3}}<\frac{1}{\kappa},$ (54)

where $\hat{x}_{3}=\bar{x}_{3}/z_{0}$. The region $\eta<1$ with
$\hat{x}_{3}>1$, in turn, corresponds to $\hat{\omega}<1$. We thus reach the
conclusion that for $\hat{\omega}$ sufficiently large ($\hat{\omega}\gtrsim
2.5$), integral range vortices do not contribute to $\hat{S}_{g}$.

###  E.1 Large frequency limit

The large $\hat{\omega}$ limit confines us to domain 2. The region
$1<\eta<\zeta$ corresponds to

$\displaystyle\hat{\omega}$ $\displaystyle<\frac{p\hat{k}}{\kappa}\ln({\rm
e}\hat{x}_{3})<\frac{\hat{k}}{\kappa}\ln({\rm e}\hat{x}_{3})$
$\displaystyle\Rightarrow\hat{k}>\frac{\kappa\hat{\omega}}{\ln({\rm
e}\hat{x}_{3})},$ (55)

and the integral in $p$ in Eq. 26 is bounded in this range by

$\displaystyle p>p_{\rm min}=\frac{\kappa}{\hat{k}^{2/3}\ln({\rm
e}\hat{x}_{3})}.$ (56)

Now, unless the detector is at the surface, the term $\hat{r}_{0}$ is large,
and thus, unless $p\simeq 0$, the factor ${\rm e}^{-2p\hat{k}\hat{r}_{0}}$ in
Eq. 26 is going to be very small. The contribution to $\hat{S}_{g}$ from the
region $\zeta>\eta>1$ is in this case exponentially damped.

To get $p_{\rm min}=0$ we need to go to the region $\eta<1$, corresponding to
the condition on $k$:

$\displaystyle\hat{k}>\hat{k}_{\rm
min}\simeq\hat{\omega}^{3/2}\hat{x}_{3}^{1/2}.$ (57)

The function $B$ (28) is $O(1)$ in the whole integration domain of Eq. 26, and
we thus get

$\displaystyle\hat{S}_{g}\sim\int_{1}^{+\infty}\hat{x}_{3}^{-4/3}{\rm
d}\hat{x}_{3}\int_{\hat{k}_{\rm min}}^{+\infty}\hat{k}^{-13/3}{\rm
d}\hat{k}\int_{0}^{1}p\,{\rm d}p\ {\rm
e}^{-2p\hat{k}\hat{r}_{0}}\sim\hat{\omega}^{-8}\hat{r}_{0}^{-2}.$ (58)

###  E.2 Small frequency limit

Let us consider first domain 1. The two conditions $\eta<1$ and
$\bar{x}_{3}<L_{\text{O}}$ imply

$\displaystyle\hat{x}_{3}<\hat{x}_{3,\rm{max}}=\min(\hat{\omega}^{-1},\hat{L}_{\rm
O}).$ (59)

The contribution to $\hat{S}_{g}$ from domain 1 is therefore

$\displaystyle\hat{S}_{g1}\sim\int_{\hat{x}_{3,\rm
min}}^{\hat{x}_{3,\rm{max}}}{\rm d}\hat{x}_{3}\int_{0}^{1/\hat{x}_{3}}{\rm
d}\hat{k}\int_{0}^{1}{\rm d}p\ g_{1}(p,\hat{k},\hat{x}_{3}),$

where $g_{1}(p,\hat{k},\hat{x}_{3})=pB\hat{x}^{3}{\rm
e}^{-2p\hat{k}\hat{r}_{0}}$. The integral is concentrated at
$(\hat{x}_{3},\hat{k})\sim(\hat{x}_{3,\rm{max}},0)$. We find the limit
behaviors

$\displaystyle\hat{S}_{g1}\sim\begin{cases}\hat{x}^{3}_{3,\rm
max},&\hat{r}_{0}\lesssim\hat{x}_{3,\rm{max}}\\\ \hat{x}_{3,\rm
max}^{5}\hat{r}_{0}^{-2},&\hat{r}_{0}\gg x_{3,\rm{max}}\end{cases}$ (60)

Let us switch to domain 2 and continue to focus on the region $\eta<1$, for
which $A\sim 1$. We have now

$\displaystyle\hat{S}_{g2}\sim\int_{\hat{x}_{3,\rm
min}}^{\hat{x}_{3,\rm{max}}}{\rm
d}\hat{x}_{3}\int_{1/\hat{x}_{3}}^{+\infty}{\rm d}\hat{k}\int_{0}^{1}{\rm d}p\
g_{2}(p,\hat{k},\hat{x}_{3}),$

where $g_{2}(p,\hat{k},\hat{x}_{3})=pB\hat{k}^{-13/3}\hat{x}_{3}^{-4/3}{\rm
e}^{-2p\hat{k}\hat{r}_{0}}$. In this case the integral is concentrated at
$(\hat{x}_{3},\hat{k})\sim(\hat{x}_{3,\rm{max}},1/\hat{x}_{3,\rm{max}})$. We
can verify that $\hat{S}_{g2}\sim\hat{S}_{g1}$, and thus recover Eqs. 32 and
33.

##  Appendix F General scaling of the noise spectra with $r_{0}$

In Section IV and Section V we have shown that, in some limiting cases $S_{g}$
scales as $1/r_{0}^{2}$ (see Eqs. 21, 31 and 32). Here we prove that such a
behavior generally arises whenever vortex decay dominates over wind advection
and the detector depth is sufficiently large, independently of the chosen
turbulence model. From Eqs. 15 and 26, we can see that $S_{g}$ can be
expressed in the general form

$\displaystyle S_{g}$ $\displaystyle\sim\int_{0}^{1}{\rm
d}p\frac{p}{\sqrt{1-p^{2}}}\int_{x_{\text{min}}}^{L}{\rm
d}x\int_{0}^{2\pi}{\rm d}\phi\cos^{2}\phi$ (61)
$\displaystyle\times\int_{0}^{\infty}{\rm
d}kf(p,x,\phi,k)h[\tau_{k}(\omega-p\cos\phi kU)]F(pkr_{0}).$

Here $f$ is the spatial part of the correlation functions, $h$ contains
information about time correlations and $F$ is a function whose form depend on
the chosen geometry. In particular, $F(pkr_{0})={\rm e}^{-2pkr_{0}}$ in our
case. We assume that:

  * •

The integral (61) is convergent, which is always the case for reasonable
models;

  * •

The function $h$ satisfies the properties given in Section IV below Eq. 14,
i.e. $h(z)$ has a maximum for $z=0$ and $h(z)\to 0$ at least exponentially for
$z\to\infty$;

  * •

$f$ is regular over the whole integration domain;

  * •

$F(pkr_{0})\to 0$ for $pkr_{0}\gg 1$ sufficiently fast (at least
exponentially).

Due to the properties of $h$, we see that the contribution to the integral
over $k$ in Eq. 61 will be peaked around some value $k=\bar{k}$, whose
specific value depends on the other parameters and on the specific model. Thus

$\displaystyle S_{g}$ $\displaystyle\sim\int_{0}^{1}{\rm
d}p\frac{p}{\sqrt{1-p^{2}}}\int_{x_{\text{min}}}^{L}{\rm
d}x\int_{0}^{2\pi}{\rm d}\phi\cos^{2}\phi$ (62)
$\displaystyle\times\bar{k}f(p,x,\phi,\bar{k})h[\tau_{\bar{k}}(\omega-p\cos\phi\bar{k}U)]F(p\bar{k}r_{0}).$

If we now assume $\omega\gg\bar{k}U$, i.e. vortex decay dominates over wind
advection, we see that the function $h$ will become independent of the wind
speed. One can now performe the integrals over $x$ and $\phi$, so that

$\displaystyle S_{g}$
$\displaystyle\sim\bar{k}h(\tau_{\bar{k}}\omega)\int_{0}^{1}{\rm
d}p\frac{p}{\sqrt{1-p^{2}}}$ (63) $\displaystyle\times
g(p,x_{\text{min}},L,\bar{k})F(p\bar{k}r_{0}).$

When $r_{0}\gg 1/\bar{k}$, we see that the only non-negligible contribution to
the integral comes from the values of $p\lesssim 1/(\bar{k}r_{0})\ll 1$,
according to the last assumption above. The integral becomes then

$\displaystyle S_{g}$
$\displaystyle\sim\bar{k}h(\tau_{\bar{k}}\omega)g(0,x_{\text{min}},L,\bar{k})F(0)\int_{0}^{1/(\bar{k}r_{0})}{\rm
d}p\,p$ (64)
$\displaystyle\sim\frac{h(\tau_{\bar{k}}\omega)g(0,x_{\text{min}},L,\bar{k})F(0)}{\bar{k}r_{0}^{2}}.$

## References

  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ [2017c] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Astrophys. J. Lett. 848, L13 (2017c), arXiv:1710.05834 [astro-ph.HE] .
  * Abbott  _et al._ [2016c] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 221101 (2016c), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. D 100, 104036 (2019), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, KAGRA, VIRGO), Astrophys. J. Lett. 915, L5 (2021b), arXiv:2106.15163 [astro-ph.HE] .
  * Somiya [2012] K. Somiya (KAGRA), Class. Quant. Grav. 29, 124007 (2012), arXiv:1111.7185 [gr-qc] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (KAGRA, LIGO Scientific, Virgo, VIRGO), Living Rev. Rel. 21, 3 (2018), arXiv:1304.0670 [gr-qc] .
  * Akutsu  _et al._ [2021] T. Akutsu _et al._ (KAGRA), PTEP 2021, 05A101 (2021), arXiv:2005.05574 [physics.ins-det] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , Class. Quant. Grav. 27, 194002 (2010).
  * Abbott  _et al._ [2017d] B. P. Abbott _et al._ (LIGO Scientific), Class. Quant. Grav. 34, 044001 (2017d), arXiv:1607.08697 [astro-ph.IM] .
  * Sathyaprakash  _et al._ [2012] B. Sathyaprakash _et al._ , Class. Quant. Grav. 29, 124013 (2012), [Erratum: Class.Quant.Grav. 30, 079501 (2013)], arXiv:1206.0331 [gr-qc] .
  * Maggiore  _et al._ [2020] M. Maggiore _et al._ , JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Harms [2019] J. Harms, Living Reviews in Relativity 22, 6 (2019).
  * Harms  _et al._ [2022] J. Harms, L. Naticchioni, E. Calloni, R. De Rosa, F. Ricci, and D. D’Urso, (2022), arXiv:2202.12841 [gr-qc] .
  * Saulson [1984] P. R. Saulson, Phys. Rev. D 30, 732 (1984).
  * Creighton [2008] T. Creighton, Class. Quant. Grav. 25, 125011 (2008), arXiv:gr-qc/0007050 .
  * Cafaro and Ali [2009] C. Cafaro and S. A. Ali, (2009), arXiv:0906.4844 [gr-qc] .
  * Stull [1988] R. B. Stull,  _An introduction to boundary layer meteorology_ , Vol. 13 (Springer Science & Business Media, 1988).
  * Schlichting and Gersten [2003] H. Schlichting and K. Gersten,  _Boundary-layer theory_ (Springer Science & Business Media, 2003).
  * Troen and Lundtang Petersen [1989] I. Troen and E. Lundtang Petersen, _European Wind Atlas_ (Risø National Laboratory, 1989).
  * Obukhov [1971] A. Obukhov, Boundary-Layer Meteorology 2, 7 (1971).
  * Kolmogorov [1991] A. N. Kolmogorov, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 434, 9 (1991).
  * Belinicher and L’vov [1987] V. Belinicher and V. L’vov, Sov. Phys. JETP 66, 303 (1987).
  * Wu and He [2021] T. Wu and G. He, Physical Review Fluids 6, 054602 (2021).
  * Naguib and Wark [1992] A. Naguib and C. Wark, Journal of Fluid Mechanics 243, 541 (1992).
  * GIT [2022] “AtmosphericNN GitHub repository,” https://github.com/maurooi/AtmosphericNN (2022).
  * Peter Lepage [1978] G. Peter Lepage, Journal of Computational Physics 27, 192 (1978).
  * Alves Junior [2018] A. A. Alves Junior, “Multithreadcorner/hydra,”  (2018).
  * Stoll  _et al._ [2020] R. Stoll, J. A. Gibbs, S. T. Salesky, W. Anderson, and M. Calaf, Boundary-Layer Meteorology 177, 541 (2020).