4 gray 90 -1.5,-10

Complete One-Loop Corrections to
$e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$ for
Different Scenarios

S. M. Seif1 , T.A. Azim1
1Faculty of Science, Physics Department, Cairo University, Giza, Egypt.
(Dated: November 20, 2014)

###### Abstract

In the present work, the full one-loop corrections to the production of a
light neutral minimal supersymmetric standard model Higgs boson ($h^{0}$) with
a pair of lightest neutralinos (${\tilde{\chi}}_{1}^{0}$) in $e^{+}e^{-}$
collisions within the Minimal Supersymmetric Standard Model (MSSM) are
presented. The details of the renormalization scheme used are presented. Our
results also include the QED corrections as well as the weak corrections. It
is found that the contribution from the weak and QED corrections is
significant and needs to be taken into account in the future linear collider
experiments. Numerical results for two different SUSY scenarios —Higgsino and
Gaugino scenarios— for
$e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$ are
given.

PACS numbers:

##  1 Introduction

One of the main goals of the Tevatron and the LHC experimental programs was to
detect a Higgs boson. On the 4th of July 2012, the CMS and the ATLAS
experimental teams at the LHC, announced independently, that they both
discovered a previously unknown boson of mass between 125 and 127 GeV [1, 2,
3], whose behavior so far is "consistent with" a Higgs boson, and it is
confirmed likely, on March 2013, to be a Higgs boson, although yet it is
unclear which model best supports the particle or whether multiple Higgs
bosons exist. This discovery has impact on the search for particles such as
neutralino [4]. Supersymmetry (SUSY) is a novel space-time symmetry between
bosons and fermions. In realistic models, SUSY is broken at the weak scale
implying that all Standard Model (SM) particles must have superpartners with
masses in the range $\sim$ 100 – 1000 GeV that will be accessible to
colliders. In the Minimal Supersymmetric Standard Model (MSSM), there are as
many as five Higgs mass states: two scalars, h0 and H0, a pseudo scalar, A0
and a pair of charged bosons, H±, which makes the experimental search more
involved. The couplings of the Higgs bosons to the SUSY scalar fermions
$\tilde{f}$ to the charginos $\tilde{\chi}^{\pm}$ and neutralinos
$\tilde{\chi}^{0}$ depend on the soft–SUSY breaking parameters and therefore
carry information on the fundamental SUSY theory.

Since the mass of neutralinos are among the precision observables with lots of
information on the SUSY-breaking structure, the relations between the particle
masses and the SUSY parameters are important theoretical quantities for
precision calculations. In (MSSM) [5], one has four neutralinos
$\tilde{\chi}_{1}^{0}$-$\tilde{\chi}_{4}^{0}$, which are the fermion mass
eigenstates of the supersymmetric partners of the photon, the $Z^{0}$ boson,
and the neutral Higgs bosons $\textit{H}_{1,2}^{0}$. Their mass matrix depends
on the parameters $M_{1}$, $M_{2}$, $\mu$, and $\tan\beta$, where $M_{1}$ and
$M_{2}$ the SU(2) and U(1) gauge mass parameter, $\tan\beta=v_{1}/v_{2}$ with
$v_{1,2}$ the vacuum expectation values of the two neutral Higgs doublet
fields. If supersymmetry is realized in nature, neutralinos should be found in
the next generation of high energy experiments at Tevatron, LHC [6] and a
future $e^{+}e^{-}$ collider. Especially at a linear $e^{+}e^{-}$ collider, it
will be possible to perform measurements with high precision [7, 8].

In view of the experimental prospects, it is inevitable to include
higher–order terms in the calculation of the measured quantities in order to
achieve theoretical predictions matching the experimental accuracy. Former
studies on chargino–pair production [9, 10, 11] and scalar–quark decays [12]
have revealed that the Born–level predictions can be influenced significantly
by one–loop radiative corrections.

In this paper, we use on-shell renormalization scheme in the loop calculations
of the Higgs and neutralino sectors of the CP-conserving MSSM. The calculation
was performed using the FeynArts and FormCalc computer packages. All the
renormalization constants, required to determine the various counterterms for
the Higgs, neutralino and other sectors, being implemented in the MSSM version
of FeynArts [13] for completion at the one-loop level. The resulting
amplitudes were algebraically simplified using FormCalc and then converted to
a FORTRAN program. The LoopTools package was used to evaluate the one-loop
scalar and tensor integrals [14].

The paper is arranged as follows: The analytical calculations of the Born
cross section to the
$e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$
process is given in section 2, where some numerical results are shown. The
virtual, the electroweak, and the soft photonic corrections are studied in
section 3. The numerical results are presented in section 4. Finally, the
conclusions are given in section 5.

##  2 Treee Level

The total cross section at the tree level for the process of production of the
Higgs bosons in association with a neutralino-pair

$\displaystyle
e^{+}(p_{1})+e^{-}(p_{2})\rightarrow{\tilde{\chi}}_{1}^{0}(p_{3})+{\tilde{\chi}}_{1}^{0}(p_{4})+h^{0}(p_{5})$

can be written as:

$\displaystyle\sigma^{0}(e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0})=\frac{(2\pi)^{4}}{4|\overline{p}_{1}|\sqrt{s}}\int\sum\limits_{spins}{|\mathcal{M}_{0}|^{2}}d\Phi_{3}.$
(1)

The integration is performed over the three-body phase space of the final
state. The three-body phase space element $d\Phi_{3}$ is read as

$\displaystyle
d\Phi_{3}=\delta^{4}(p_{1}+p_{2}-\sum\limits_{i=3}^{5}p_{i})\prod\limits_{j=3}^{5}\frac{d^{3}P_{j}}{(2\pi)^{3}2E_{j}},$

from which we get

$\displaystyle\sigma^{0}=\frac{|\mathcal{M}_{0}|^{2}}{16(2\pi)s}\\{\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4}-p_{5})\\}$
$\displaystyle\left\\{\
\frac{d^{3}p_{3}}{E_{3}}\frac{d^{3}p_{4}}{E_{4}}\frac{d^{3}p_{5}}{E_{5}}\right\\},$
(2)

where $p_{1}$, $p_{2}$, $p_{3}$, $p_{4}$ and $p_{5}$ are the four-momenta of
the incoming and the outgoing particles, respectively.
$\lvert\mathcal{M}_{0}\rvert^{2}$ is the square of amplitudes corresponding to
the Feynman diagrams in Fig. 1 and $s=(p_{1}+p_{2})^{2}$ is the square of the
total energy in the centre of mass system. The Feynman diagrams contributing
to the production of the lightest CP–even Higgs boson in association with
neutralino pairs is shown in Fig. 1. The diagrams where the $h^{0}$ boson is
emitted from the electron and positron lines give negligible contributions. A
first class of contributions (a) is formed by diagrams where the Higgs boson
is emitted from the neutralino states, the latter being produced through
s–channel $Z$ boson exchange and left– and right–handed selectron exchange. A
second class (b) is formed by the Higgs–strahlung production process, where
the $Z$ boson is virtual and splits into two neutralinos. Finally, a third
class (c) consists of the diagrams where the Higgs boson is emitted from the
internal selectron lines. The cross section will therefore depend on the
$h^{0}$ boson couplings to both the neutralinos and sleptons [15]. In the
MSSM, the Higgs boson couplings to the neutralinos
$\tilde{\chi}_{(1,\cdots,4)}^{0}$ depend on $\tan\beta$, the higgsino
parameter $\mu$ and the gaugino —bino and wino— mass parameters $M_{1}$ and
$M_{2}$.

Figure 1: The lowest order (LO) Feynman diagrams for the
$e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$
process.

As clarified in [16, 17], all Higgs bosons in the MSSM, except the lightest
CP–even one, are too heavy to play an important role in both the current and
the near future experiments. Therefore, the present study concentrates on the
lightest Higgs boson $h^{0}$ only. The LO cross-section is studied for the
process
$e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$ as a
function of the light neutral Higgs boson mass $m_{h^{0}}$. From the
interaction Lagrangian [18, 19]

$\displaystyle\mathcal{L}_{Z^{0}\bar{e}e}$
$\displaystyle=-\frac{g}{\cos\theta_{W}}Z_{\mu}^{0}\bar{e}\gamma^{\mu}[C_{L}P_{L}+C_{R}P_{R}]e,$
$\displaystyle\mathcal{L}_{Z^{0}\tilde{\chi}_{i}^{0}\tilde{\chi}_{j}^{0}}$
$\displaystyle=\frac{g}{2\cos\theta_{W}}Z_{\mu}^{0}\bar{\tilde{\chi}}_{i}^{0}\gamma^{\mu}[O_{ij}^{\prime\prime
L}P_{L}+O_{ij}^{\prime\prime R}P_{R}]\tilde{\chi}_{j}^{0},$
$\displaystyle\mathcal{L}_{e\tilde{e}\tilde{\chi}_{i}^{0}}$
$\displaystyle=gf_{i}^{L}\bar{e}P_{R}\tilde{\chi}_{i}^{0}\tilde{e}_{L}+gf_{i}^{R}\bar{e}P_{L}\tilde{\chi}_{i}^{0}\tilde{e}_{R}+h.c.,$
$\displaystyle\mathcal{L}_{h^{0}\tilde{\chi}_{i}^{0}\tilde{\chi}_{j}^{0}}$
$\displaystyle=-\frac{g}{2}h^{0}\cos\alpha\bar{\tilde{\chi}}_{i}^{0}[Q_{ij}^{\prime\prime*}P_{L}+Q_{ij}^{\prime\prime}P_{R}]\tilde{\chi}_{j}^{0}$
$\displaystyle+\frac{g}{2}h^{0}\sin\alpha\bar{\tilde{\chi}}_{i}^{0}[S_{ij}^{\prime\prime*}P_{L}+S_{ij}^{\prime\prime}P_{R}]\tilde{\chi}_{j}^{0},$
$\displaystyle\mathcal{L}_{h^{0}ZZ}$
$\displaystyle=\frac{gm_{Z}}{2\cos\theta_{W}}Z_{\mu}Z^{\mu}h^{0}\cos(\beta-\alpha),$
$\displaystyle\mathcal{L}_{h^{0}\tilde{e}\tilde{e}}$
$\displaystyle=\frac{gm_{Z}}{\cos\theta_{W}}h^{0}\left[(I^{3L}-e_{\tilde{e}}\sin^{2}\theta_{W})\tilde{e}^{*}_{L}\tilde{e}_{L}+e_{\tilde{e}}\sin^{2}\theta_{W}\tilde{e}^{*}_{R}\tilde{e}_{R}\right]$

one obtains the following couplings:

$\displaystyle C_{L,R}$ $\displaystyle=I^{3L,R}+\sin^{2}\theta_{W},\quad
I^{3L}=-\frac{1}{2},\quad I^{3R}=0.$ $\displaystyle O_{ij}^{\prime\prime L}$
$\displaystyle=-O_{ij}^{\prime\prime
R*}=-\frac{1}{2}N_{i3}N_{j3}^{*}+N_{i4}N_{j4}^{*}.$ $\displaystyle f^{L}_{i}$
$\displaystyle=-\frac{\sqrt{2}}{2}\left(\tan\theta_{W}N_{i1}+N_{i2}\right),\quad
f_{i}^{R}=\sqrt{2}\tan\theta_{W}N_{i1}^{*}.$ $\displaystyle
gQ_{ij}^{\prime\prime}$
$\displaystyle=\frac{1}{2}[N_{i3}(gN_{j2}-g^{\prime}N_{i1})+\sqrt{2}h^{*}N_{i4}N_{j5}+(i\leftrightarrow
j)].$ $\displaystyle gS_{ij}^{\prime\prime}$
$\displaystyle=\frac{1}{2}[N_{i4}(gN_{j2}-g^{\prime}N_{i1})+\sqrt{2}h^{*}N_{i3}N_{j5}+(i\leftrightarrow
j)].$ $\displaystyle R_{ij}^{\prime\prime}$
$\displaystyle=\frac{1}{2m_{W}}[M^{*}N_{i2}N_{j2}+M^{\prime*}N_{i1}N_{j1}-M^{*}(N_{i3}N_{j4}+N_{i4}N_{j3})].$

Figure 2: Vertex Corrections. Figure 3: Box Corrections. Figure 4: Propagator
Corrections.

##  3 One-Loop Calculations

The linear electron-positron collider is considered to be the best environment
for precise studies of supersymmetric models, especially for the MSSM. From
precise measurements of masses and cross sections, the fundamental parameters
of the MSSM Lagrangian can be reconstructed [20] to shed light on the
mechanism of SUSY breaking. Adequate theoretical predictions matching the
experimental accuracy require the proper inclusion of higher-order terms in
the calculations of mass spectra, cross sections and decay rates.

The radiative corrections to the neutralino pair production with light neutral
Higgs boson can be classified as the following generic structure of one-loop
Feynman diagrams: The virtual vertex corrections Fig. 2, the box graph
contributions to the propagators. Fig. 3, and the self-energy contributions
Fig. 4. For the virtual particles inside loops, the complete supersymmetric
spectrum is used. Self-energy and vertex diagrams contain IR and UV
divergences. Box diagrams are ultraviolet finite, but have infrared
singularities. The IR poles in the virtual corrections are eventually canceled
by analogous singularities in the real corrections to the tree-level cross
section. The UV divergences are renormalized by introducing a suitable set of
counterterms for the renormalization of the coupling constants and the
renormalization of the external wave functions.

As a calculational frame, the on-shell renormalization scheme has been chosen
such that all particle masses are defined as pole masses, i.e., on-shell
quantities. Thus, the cross sections are directly related to the physical
masses of the external particles and the other particles entering the loops.
The theoretical basis of the calculation is outlined in the following
subsections, where the presentation is based on refs.[21, 22].

In this work, the complete set of Feynman graphs is calculated with help of
the packages FeynArts-3.6, FormCalc-7.1, and LoopTools-2.7. We have
implemented our renormalization procedure into these packages. The virtual
correction suffers from both ultraviolet (UV) and infrared (IR) singularities.
For a proper treatment of the appearing UV divergencies, counter terms are
introduced in the on-shell renormalization scheme where all particle masses
are defined as pole masses, i.e., on-shell quantities. To preserve
supersymmetry, the used regularization scheme is dimensional reduction (DR).
The loop graphs with virtual photon exchange also introduce IR singularities.
Therefore, real photon emission has to be included to obtain a finite result.

$\displaystyle\sigma^{corr}(e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0})$
$\displaystyle=\sigma^{ren}(e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0})$
$\displaystyle\quad+\sigma(e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}\gamma).$

###  3.1 Renormalization of Neutralino

The bilinear part of the Lagrangian describing the neutralino sector of the
MSSM involves the $\mu$ parameter, the soft–breaking gaugino-mass parameters
$M_{1}$ and $M_{2}$, and the Higgs vacua $v_{i}$, which are related to
$\tan\beta=v_{2}/v_{1}$ and to the $W$ mass $M_{W}=gv/2$ with
$(v_{1}^{2}+v_{2}^{2})^{1/2}$ [21, 22].

Renormalization constants are introduced for the neutralino mass matrix $Y$

$Y=\left(\begin{array}[]{cccc}M_{1}&0&-M_{Z}s_{W}\cos\beta&M_{Z}s_{W}\sin\beta\\\
0&M_{2}&M_{Z}c_{W}\cos\beta&-M_{Z}c_{W}\sin\beta\\\
-M_{Z}s_{W}\cos\beta&M_{Z}c_{W}\sin\beta&0&-\mu\\\
M_{Z}s_{W}\cos\beta&-M_{Z}c_{W}\sin\beta&-\mu&0\end{array}\right)$

and for the neutralino fields $\tilde{\chi}_{i}^{0}(i=1,\ldots,4)$ by the
following transformation

$\displaystyle Y$ $\displaystyle\rightarrow Y+\delta Y,$
$\displaystyle\omega_{L}\tilde{\chi}_{i}^{0}$
$\displaystyle\rightarrow(\delta_{ij}+\frac{1}{2}[\delta
Z_{\tilde{\chi}^{0}}]_{ij})\omega_{L}\tilde{\chi}_{j}^{0},$
$\displaystyle\omega_{R}\tilde{\chi}_{i}^{0}$
$\displaystyle\rightarrow(\delta_{ij}+\frac{1}{2}[\delta
Z_{\tilde{\chi}^{0}}]_{ij}^{*})\omega_{R}\tilde{\chi}_{j}^{0}.$ (3)

The matrix $\delta{Y}$ consists of the counterterms for the following
parameters in the mass matrix $Y$: the soft-breaking gaugino-mass parameters
$M_{1}$ and $M_{2}$, the Higgsino mass parameter $\mu$, $\tan\beta$, the $Z$
mass and the electroweak mixing angle $s_{W}=\sin\theta_{W}$,
$c_{W}=\cos\theta_{W}$. The matrix-valued renormalization constant $\delta
Z_{\tilde{\chi}^{0}}$is a general complex $4\times 4$ matrix.

Using the on-shell approach [22, 23], the pole mass of the neutralino,
$\tilde{\chi}_{1}^{0}$, is considered as an input parameters to specify the
neutralino Lagrangian in terms of physical quantities. This is equivalent to
the specification of the parameters $\mu$, $M_{1}$ and $M_{2}$, which are
related to the input masses in the same way as in LO, as a consequence of the
on-shell renormalization conditions, and to the definition of the respective
counterterms. In this way, the tree-level masses $m_{\tilde{\chi}_{1}^{0}}$ as
well as the counterterm matrix $\delta Z_{\tilde{\chi}^{0}}$, are fixed.
Furthermore, one requires that the matrix of the renormalized one-particle-
irreducible two-point vertex functions $\hat{\Gamma}_{ij}^{(2)}$ becomes
diagonal for the on-shell external momenta. This fixes the non-diagonal
entries of the field-renormalization matrix $\delta Z_{\tilde{\chi}^{0}}$;
their diagonal entries are determined by normalizing the residues of the
propagators to be unity.

###  3.2 Renormalization of Higgs Sector

As known the MSSM requires two Higgs doublets $H_{1}$ and $H_{2}$ with
opposite hypercharge $Y_{1}=-Y_{2}=-1$. The quadratic part of the Higgs
potential in the MSSM is given by:

$\displaystyle
V=m_{1}^{2}H_{1}\bar{H}_{1}+m_{2}^{2}H_{2}\bar{H}_{2}+m_{12}^{2}(\epsilon_{ab}H_{1}^{a}H_{2}^{b}+h.c.)$
$\displaystyle+\frac{1}{8}(g_{1}^{2}+g_{2}^{2})(H_{1}\bar{H}_{1}-H_{2}\bar{H}_{2})^{2}-\frac{g_{2}^{2}}{2}{\lvert
H_{1}\bar{H}_{2}\rvert}^{2},$ (4)

where $m_{12}^{2}$ is defined to be negative and
$\epsilon_{12}=-\epsilon_{21}=-1$, with soft breaking parameters
$m_{1}^{2},m_{2}^{2},m_{12}^{2}$ and $g_{1},g_{2}$ are, respectiveley, SU(2)
and U(1) gauge couplings. Decomposing each Higgs doublet field $H_{1,2}$ in
terms of its components [23, 24],

$\displaystyle H_{1}=\left(\begin{array}[]{c}H_{1}^{1}\\\
H_{1}^{2}\end{array}\right)=\left(\begin{array}[]{c}(v_{1}+\phi_{1}^{0}-i\chi_{1}^{0})/\sqrt{2}\\\
-\phi_{1}^{-}\end{array}\right)$ (9)

$\displaystyle H_{2}=\left(\begin{array}[]{c}H_{2}^{1}\\\
H_{2}^{2}\end{array}\right)=\left(\begin{array}[]{c}\phi_{2}^{+}\\\
(v_{2}+\phi_{2}^{0}-i\chi_{2}^{0})/\sqrt{2}\end{array}\right)$ (14)

with vacuum expectation values $v_{1}$, $v_{2}$.

the Higgs potential (3.2) is diagonalized by the rotations

$\displaystyle\left(\begin{array}[]{c}H^{0}\\\
h^{0}\end{array}\right)=\left(\begin{array}[]{cc}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{array}\right)\left(\begin{array}[]{c}\phi_{1}^{0}\\\
\phi_{2}^{0}\end{array}\right)$ (21)

$\displaystyle\left(\begin{array}[]{c}G^{0}\\\
A^{0}\end{array}\right)=\left(\begin{array}[]{cc}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{array}\right)\left(\begin{array}[]{c}\chi_{1}^{0}\\\
\chi_{2}^{0}\end{array}\right)$ (28)

$\displaystyle\left(\begin{array}[]{c}G^{+}\\\
H^{+}\end{array}\right)=\left(\begin{array}[]{cc}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{array}\right)\left(\begin{array}[]{c}\phi_{1}^{+}\\\
\phi_{2}^{+}\end{array}\right)$ (35)

$G^{0}$, $G^{\pm}$ describe the unphysical Goldstone modes. The spectrum of
physical states consists of: a light neutral CP-even state ($h^{0}$), a heavy
neutral CP-even state ($H^{0}$), a neutral CP-odd state ($A^{0}$), and a pair
of charged states ($H^{\pm}$).

The masses of the gauge bosons and the electromagnetic charge are determined
by

$\displaystyle M_{Z}^{2}$
$\displaystyle=\frac{1}{4}(g_{1}^{2}+g_{2}^{2})(v_{1}^{2}+v_{2}^{2}),$
$\displaystyle M_{W}^{2}$ $\displaystyle=g_{2}^{2}(v_{1}^{2}+v_{2}^{2}),$
$\displaystyle e^{2}$
$\displaystyle=\frac{g_{1}^{2}g_{2}^{2}}{g_{1}^{2}+g_{2}^{2}}.$ (36)

Thus, the potential (3.2) contains two independent free parameters, which can
conveniently be chosen as

$\displaystyle\tan\beta=\frac{v_{2}}{v_{1}},M_{A}^{2}=-m_{12}^{2}(\tan\beta+\cot\beta)$
(37)

where $M_{A}$ is the mass of the $A^{0}$ boson.

Expressed in terms of (37), the masses of the other physical states read:

$\displaystyle m_{H^{0},h^{0}}^{2}$
$\displaystyle=\frac{1}{2}\biggl{[}M_{A}^{2}+M_{Z}^{2}\pm\sqrt{(M_{A}^{2}+M_{Z}^{2})^{2}-4M_{A}^{2}M_{Z}^{2}\cos^{2}2\beta}\biggr{]}$
$\displaystyle m_{H^{+}}^{2}$ $\displaystyle=M_{A}^{2}+M_{W}^{2},$ (38)

and the mixing angle $\alpha$ in the ($H^{0},h^{0}$)-system is derived from

$\displaystyle\tan 2\alpha=\tan
2\beta\frac{M_{A}^{2}+M_{Z}^{2}}{M_{A}^{2}-M_{Z}^{2}},-\frac{\pi}{2}<\alpha\leq
0.$ (39)

Hence, masses and couplings are determined by only a single parameter more
than in the standard model.

The dependence on $M_{A}$ is symmetric under
${\tan\beta}\leftrightarrow{1/\tan\beta}$, and $m_{h^{0}}$ is constrained by:

$\displaystyle m_{h^{0}}<M_{Z}\cos 2\beta<M_{Z}.$ (40)

This simple scenario, however, is changed when radiative corrections are taken
into account.

The tree-level mass matrix $m_{0}$ of the neutral scalar system that
represents bare mass system is diagonalized by (21). Loop contributions to the
quadratic part of the potential (neglecting the $q^{2}$-dependence of the
diagrams) modify the mass matrix

$\displaystyle m_{0}\rightarrow m_{0}+\delta m=m$ (41)

Re-diagonalizing the one-loop matrix $m$ yields the corrected mass eigenvalues
$m_{H^{0},h^{0}}$, replacing (3.2), and an effective mixing angle
$\alpha_{eff}$ instead of (39). The renormalization constants [23] are defined
as follows:

$\displaystyle B_{\mu}$ $\displaystyle\rightarrow(Z_{2}^{B})^{1/2}B_{\mu},$
$\displaystyle W_{\mu}^{a}$
$\displaystyle\rightarrow(Z_{2}^{W})^{1/2}W_{\mu}^{a},$ $\displaystyle H_{i}$
$\displaystyle\rightarrow Z_{H_{i}}^{1/2}H_{i},$ $\displaystyle\psi_{j}^{L}$
$\displaystyle\rightarrow{(Z_{L}^{j})}^{1/2}\psi_{j}^{L},$
$\displaystyle\psi_{j\sigma}^{R}$
$\displaystyle\rightarrow{(Z_{R}^{j\sigma})}^{1/2}\psi_{j\sigma}^{R},$
$\displaystyle g_{2}$ $\displaystyle\rightarrow
Z_{1}^{W}(Z_{2}^{W})^{-3/2}g_{2},$ $\displaystyle g_{1}$
$\displaystyle\rightarrow Z_{1}^{B}(Z_{2}^{B})^{-3/2}g_{1},$ $\displaystyle
v_{i}$ $\displaystyle\rightarrow Z_{H_{i}}^{1/2}(v_{i}-\delta v_{i}),$
$\displaystyle m_{i}^{2}$ $\displaystyle\rightarrow
Z_{H_{i}}^{-1}(m_{i}^{2}+\delta m_{i}^{2}),$ $\displaystyle m_{12}^{2}$
$\displaystyle\rightarrow Z_{H_{1}}^{-1/2}Z_{H_{2}}^{-1/2}(m_{12}^{2}+\delta
m_{12}^{2}).$ (42)

The complete definitions and the explicit expressions of the renormalization
constants of the other sectors: sfermion sector, MSSM parameters and fields
including those of SM as the electric charge and the masses of $W,Z$, and the
fermions and their counterterms in addition to $\tan\beta$, all these are
treated as described in [23, 26, 27], to deliver all counterterms required for
propagators and vertices appearing in the amplitudes.

###  3.3 Virtual Corrections

The set of Feynman diagrams in Fig. 1, has to be dressed by the corresponding
loop contributions containing the full particle spectrum of the MSSM. The
complete cross section at the one-loop level can be written as follows:

$\displaystyle\sigma^{1-loop}$ $\displaystyle=\sigma^{0}+\sigma^{virt}$ (43)

The virtual electroweak radiative correction to the cross section can be
written as

$\displaystyle\sigma^{virt}$
$\displaystyle=\sigma^{0}\Delta_{virt}=\frac{(2\pi)^{4}}{2|\overline{p}_{1}|\sqrt{s}}\int
d\Phi_{3}\sum\limits_{spins}Re(\mathcal{M}^{\dagger}_{0}\mathcal{M}_{virt})$
(44)

where $\sigma^{0}$ and $\mathcal{M}_{0}$ are the cross section and amplitude
at the tree-level for the
$e^{+}e^{-}\rightarrow{\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$
process, respectively. $\mathcal{M}_{virt}$ is the renormalized amplitude
involving all the one-loop electroweak Feynman diagrams and corresponding
counter-terms.

###  3.4 Real Photonic Corrections

The soft (IR) divergences in the $\mathcal{M}_{virt}$ originate from the
contributions of virtual photon exchange in loops [28]. These soft (IR)
divergencies can be cancelled by the real photon bremsstrahlung corrections in
the soft photon limit in which the cross-section for real photon emission [29]
is proportional to the Born cross-section,

$\displaystyle{\biggl{(}\frac{d\sigma}{d\Omega}\biggr{)}}_{soft}={\biggl{(}\frac{d\sigma}{d\Omega}\biggr{)}}_{Born}\delta_{soft}$
(45)

where $\delta_{soft}$ is the QED correction factor from real bremsstra- hlung
in the soft–photon approximation. The real photonic corrections are induced by
the process

$\displaystyle
e^{+}(p_{1})+e^{-}(p_{2})\rightarrow{\tilde{\chi}}_{1}^{0}(p_{3})+{\tilde{\chi}}_{1}^{0}(p_{4})+h^{0}(p_{5})+\gamma(k_{\gamma})$

where $k_{\gamma}$ denotes the photon momentum.

A real photon radiates from the electron/ positron, and can have either soft
or collinear nature. The collinear singularity is regularized by keeping
electron (positron) mass. The general phase-space-slicing method (PSS) is
adopted to separate the soft photon emission singularity from the real photon
emission processes. By using this method, the bremsstrahlung phase space is
divided into singular and non-singular regions by the soft photon cut off,
$\Delta E$, and the cross section of the real photon emission process is
decomposed into soft and hard terms [30],

$\displaystyle\sigma^{real}$ $\displaystyle=\sigma^{soft}(\Delta
E)+\sigma^{hard}(\Delta E)$
$\displaystyle=\sigma^{0}(\Delta_{soft}+\Delta_{hard}).$ (46)

where $\Delta_{soft}=\Delta_{weak}$ and $\Delta_{hard}=\Delta_{QED}$

The energy of the radiated photon in the center of mass system frame is
considered as a soft term, $\Delta_{soft}$, with radiated photon energy
$k_{\gamma}^{0}<\Delta E$, and a hard term, $\delta_{hard}$, with
$k_{\gamma}^{0}>\Delta E$ , where
$k_{\gamma}^{0}=\sqrt{{\lvert\vec{k}_{\gamma}\rvert^{2}}+m_{\gamma}^{2}}$ and
$m_{\gamma}$ is the photon mass, which is used to regulate the (IR)
divergences existing in the soft term. For practical calculations,
$\sigma^{hard}$ is divided into a collinear part, where the photon is within
an angle smaller than $\Delta\theta$ with respect to the radiating particles,
and the complementary non-collinear part [23],

$\displaystyle\sigma^{hard}=\sigma^{coll}(\Delta\theta)+\sigma^{non-
coll}(\Delta\theta).$ (47)

$\sigma^{non-coll}$ is calculated numerically with the help of
multidimensional numerical integration routines DIVONNE that based on Monte
Carlo, and CUHRE, which are both part of the CUBA-library [31], while
$\sigma^{soft}$ and $\sigma^{coll}$ are calculated analytically.

The weakness of the definition in Eq.(3.4) is the larg $\Delta E$ dependence
of the weak and QED components $\propto log(\frac{\Delta E^{2}}{s})$. However,
the sum of both is cutoff independent. Therefore, we extract the $\Delta E$
terms and the leading logarithms $log(\frac{s}{m_{e}^{2}})$, caused by
collinear soft photon emission, from the weak corrections and add them to the
QED corrections [32]. According to this definition, both corrections are now
$\Delta E$ independent. The main part of the QED corrections arises from these
leading logarithms $L_{e}$, originating from photons in the beam direction.
This leads to a large dependence on the experimental cuts and detector
specifications. Therefore, We use the structure function formalism [33] and
subtract the leading logarithmic $\mathcal{O}(\alpha)$ terms of the initial
state radiation, $\sigma^{ISR,LL}$ . After subtraction of these process
independent terms, only the non-universal QED corrections remain. The total
one-loop (renormalized) cross section$\sigma^{1-loop}$ is expressed as:

$\displaystyle\sigma^{1-loop}$
$\displaystyle=\sigma^{0}+\sigma^{virt}+\sigma^{weak}+\sigma^{QED},$

where

$\displaystyle\sigma^{weak}$ $\displaystyle=\sigma^{soft}-\tilde{\sigma}.$
$\displaystyle\sigma^{QED}$
$\displaystyle=\sigma^{hard}+\tilde{\sigma}-\sigma^{ISR,LL}.$

with

$\displaystyle\tilde{\sigma}=\frac{\alpha}{\pi}\left[\left(L_{e}-1\right)\log\frac{4\Delta
E^{2}}{s}+\frac{3}{2}L_{e}\right]\sigma^{0},$
$\displaystyle\sigma^{ISR,LL}=\frac{\alpha}{\pi}L_{e}\int\limits_{0}^{1}dx\Phi\left(x\right)\sigma^{0}\left(sx\right),$
$\displaystyle\Phi\left(x\right)=\lim\limits_{\epsilon\rightarrow
0}\left\\{\delta\left(1-x\right)\left[\frac{3}{2}+2\log(\epsilon)\right]+\theta\left(1-x-\epsilon\right)\frac{1+x^{2}}{1-x}\right\\}.$

The integrated cross section at the one–loop level, can be written in the
following way:

$\displaystyle\sigma^{1-loop}$ $\displaystyle=\sigma^{0}+\sigma^{0}\Delta,$

pointing out the relative correction

$\displaystyle\Delta$
$\displaystyle=\left(\sigma^{1-loop}-\sigma^{0}\right)/\sigma^{0},$

with respect to the Born cross section. The relative correction $\Delta$ can
be decomposed into the following parts, indicating their origin,

$\displaystyle\Delta$
$\displaystyle=\Delta_{self}+\Delta_{vertex}+\Delta_{box}+\Delta_{QED}+\Delta_{weak}.$

##  4 Numerical Results

In present work, two different scenarios are studied. In the higgsino scenario
the neutralinos are both nearly pure higgsinos and therefore the process is
dominated by the schannel $Z_{0}$ exchange. In the gaugino scenario with binos
as $\tilde{\chi}_{1}^{0}$ states, the selectron exchange diagrams play the
most important role. In the following, we distinguish between the tree-level,
the 1-loop corrections, and the conventional weak and QED corrections to the
improved tree-level as discussed in the last section. For the SM input
parameters the following values have been used $\alpha(m_{Z})=1/127.922$,
$M_{W}=80.399$ GeV, $M_{Z}=91.1876$ GeV,
$\sin^{2}\theta_{W}=1-m_{W}^{2}/m_{Z}^{2}$ $m_{t}=174.3$ GeV, $m_{b}=4.7$ GeV.
$M_{2}$ is fixed by the gaugino unification relation
$M_{2}=\frac{5}{3}\tan^{2}\theta_{W}M_{1}$, and the gluino mass is related
$M_{1}$ by
$m_{\tilde{g}}=\left(\alpha_{s}\left(m_{\tilde{g}}\right)/\alpha\right)\sin^{2}\theta_{W}M_{1}$.
In the following numerical examples, we used SoftSusy code for the mass
spectrum calculations of SUSY particles in the minimal supergravity (mSUGRA)
framework. We consider two benchmark scenarios:

  * •

for $\mu\ll M_{2}$, we have pure higgsino.

  * •

for $M_{2}\ll\mu$, we have pure gaugino.

the mass spectrum of the SUSY particles are set as shown in Tables 1 and 2 for
Higgsino and Gaugino scenarios respectively. The free parameters that have
been used in our calculations are specified as follows:

  * •

For simplicity, in the sfermions sector, all soft-SUSY breaking parameters are
assumed equal and all trilinear couplings are set to a common value $A$. The
mixing between sfermion generations is neglected,
$M_{SUSY}\equiv\tilde{M}_{L}\simeq\tilde{M}_{R}$.

  * •

The MSSM Higgs sector is parametrized by the CP-odd mass, $m_{A}$, and
$\tan\beta$, taking into account radiative corrections with the help of
FormCalc.

  * •

The chargino-neutralino sector is fixed by choosing a value for the gaugino-
mass terms $M_{1}$, $M_{2}$ and for the Higgsino-mass term $\mu$.

Numerical results for the radiative corrections to the
${\tilde{\chi}}_{1}^{0}{\tilde{\chi}}_{1}^{0}h^{0}$ production are given in
Figs. 5 and 6 for Higgsino scenario. The maximum value of the relative weak
correction is about 28%, while that for QED correction is about -57% in the
investigated parameters region. In Figs. 7 and 8, which are related to Gaugino
scenario, The maximum value of the weak correction is about 24%, while that
for QED correction is about -40% in the investigated parameters region. Thus,
for two scenarios, the weak and QED corrections have to be taken into account
in future experiments.

Table 1: The mass spectrum of the SUSY particles for Higgsino scenario Particle | Mass/[GeV] | Particle | Mass/[GeV]
---|---|---|---
$h^{0}$ | 105.341 | $\tilde{\chi}_{1}^{0}$ | 86.3240
$H^{0}$ | 700.275 | $\tilde{\chi}_{2}^{0}$ | 111.646
$A^{0}$ | 700.000 | $\tilde{\chi}_{3}^{0}$ | 200.218
$H^{\pm}$ | 704.600 | $\tilde{\chi}_{4}^{0}$ | 416.025
$\tilde{g}$ | 1063.46 | $\tilde{\nu}_{e}$ | 344.128
$\tilde{\chi}_{1}^{\pm}$ | 99.2230 | $\tilde{\nu}_{\mu}$ | 344.128
$\tilde{\chi}_{2}^{\pm}$ | 416.032 | $\tilde{\nu}_{\tau}$ | 342.105
$\tilde{e}_{L}$ | 352.582 | $\tilde{e}_{R}$ | 353.214
$\tilde{\mu}_{L}$ | 352.519 | $\tilde{\mu}_{R}$ | 353.277
$\tilde{\tau}_{L}$ | 349.346 | $\tilde{\tau}_{R}$ | 356.424
$\tilde{u}_{L}$ | 345.881 | $\tilde{u}_{R}$ | 348.268
$\tilde{d}_{L}$ | 350.860 | $\tilde{d}_{R}$ | 354.925
$\tilde{c}_{L}$ | 345.667 | $\tilde{c}_{R}$ | 348.485
$\tilde{s}_{L}$ | 350.853 | $\tilde{s}_{R}$ | 354.932
$\tilde{t}_{L}$ | 281.995 | $\tilde{t}_{R}$ | 469.650
$\tilde{b}_{L}$ | 343.316 | $\tilde{b}_{R}$ | 362.288
Table 2: The mass spectrum of the SUSY particles for Gaugino scenario Particle | Mass/[GeV] | Particle | Mass/[GeV]
---|---|---|---
$h^{0}$ | 110.985 | $\tilde{\chi}_{1}^{0}$ | 91.5340
$H^{0}$ | 393.987 | $\tilde{\chi}_{2}^{0}$ | 181.009
$A^{0}$ | 393.600 | $\tilde{\chi}_{3}^{0}$ | 359.502
$H^{\pm}$ | 401.727 | $\tilde{\chi}_{4}^{0}$ | 378.874
$\tilde{g}$ | 525.351 | $\tilde{\nu}_{e}$ | 495.905
$\tilde{\chi}_{1}^{\pm}$ | 180.516 | $\tilde{\nu}_{\mu}$ | 495.905
$\tilde{\chi}_{2}^{\pm}$ | 379.562 | $\tilde{\nu}_{\tau}$ | 493.614
$\tilde{e}_{L}$ | 501.812 | $\tilde{e}_{R}$ | 502.257
$\tilde{\mu}_{L}$ | 501.586 | $\tilde{\mu}_{R}$ | 502.483
$\tilde{\tau}_{L}$ | 495.440 | $\tilde{\tau}_{R}$ | 508.551
$\tilde{u}_{L}$ | 497.123 | $\tilde{u}_{R}$ | 498.787
$\tilde{d}_{L}$ | 500.591 | $\tilde{d}_{R}$ | 503.475
$\tilde{c}_{L}$ | 497.107 | $\tilde{c}_{R}$ | 498.807
$\tilde{s}_{L}$ | 500.557 | $\tilde{s}_{R}$ | 503.508
$\tilde{t}_{L}$ | 504.356 | $\tilde{t}_{R}$ | 548.374
$\tilde{b}_{L}$ | 484.474 | $\tilde{b}_{R}$ | 519.044

The parameters for the Higgsino Scenarioare set as: {$M_{2}$, $\mu$, $A$,
$\tan\beta$, $m_{A}$, $M_{SUSY}$}= {400 GeV, $-$100 GeV, 400 GeV, 10, 700 GeV,
350 GeV}, with SUSY mass spectrum shown in table 1.

The parameters for the Gaugino scenario are set as: {$M_{2}$, $\mu$, $A$,
$\tan\beta$, $m_{A}$, $M_{SUSY}$}= {197.6 GeV, 353.1 GeV, $-$100 GeV, 10.2,
393.6 GeV, 500 GeV}, with SUSY mass spectrum shown in table 2.

Figure 5: Total cross section as a function of $\sqrt{s}$ in the Higgsino scenario. Figure 6: Relative corrections as a function of $\sqrt{s}$ in the Higgsino scenario. Table 3: The maximum cross section in Higgsino scenario. |  $(\sigma)_{max}/$Pb |  $\sqrt{s}/$GeV
---|---|---
Born | $2.82\times 10^{-6}$ | 700
1–loop | $3.28\times 10^{-6}$ | 700
QED | $1.15\times 10^{-6}$ | 650
Weak | $3.51\times 10^{-6}$ | 700
Figure 7: Total cross section as a function of $\sqrt{s}$ in the Gaugino scenario. Figure 8: Relative corrections as a function of $\sqrt{s}$ in the Gaugino scenario. Table 4: The maximum cross section in Gaugino scenario. |  $(\sigma)_{max}/$Pb |  $\sqrt{s}/$GeV
---|---|---
Born | $8.59\times 10^{-8}$ | 1050
1–loop | $1.21\times 10^{-7}$ | 1000
QED | $5.65\times 10^{-8}$ | 850
Weak | $1.05\times 10^{-7}$ | 1000

##  5 Conclusions

The full one-loop corrections to the lightest neutralino pair production with
light neutral Higgs boson in $e^{+}e^{-}$ collisions have been calculated. The
calculation was performed in an analytical way with an independent check using
the FeynArts and FormCalc computer packages. The chosen renormalization scheme
can be used for the complete MSSM parameter space. A particular attention is
paid to an appropriate definition of the weak and QED corrections. We
extracted the non-universal QED corrections by subtracting the initial state
radiation (ISR). The full one-loop corrections are in the range of 12-20% for
Higgsino scenario, and of 40-70% for Gaugino scenario, thus they have to be
taken into account in future linear collider experiments. The maximum cross
sections are presented in Tables 3 and 4 for both scenarios.

## Acknowledgments

We would like to express our sincere gratitude to our advisor Pro. Dr. M.
Khaled Hegab, and our thanks to Pro. Dr. Samiha Abou Stiet for revising the
paper.

## References

  * [1] A. Cho, Science 337, 141 (2012)
  * [2] The ATLAS Collaboration, Phys. Lett. B 716, 1 (2012)
  * [3] The ATLAS Collaboration, Phys. Rev. D 86, 032003 (2012)
  * [4] C.Beskidt, W.de Boer, D.I.Kazakov. arXiv:1402.4650 [hep-ph]
  * [5] H. P. Nilles, Phys. Rep. 110, 1 (1984); H. E. Haber and G. L. Kane, Phys. Rep. 117, 75 (1985); R. Barbieri, Riv. Nuovo. Cim. 11, 1 (1988)
  * [6] The ATLAS Collaboration. arXiv:1403.5294 [hep-ex]
  * [7] J.A. Aguilar et al., in TESLA Technical Design Report, Part III: Physics at an $e^{+}e^{-}$ Linear Collider, 2001, ed. by R. D. Heuer, D. Miller, F. Richard, and P. M. Zerwas. arXiv:hep-ph/0106315
  * [8] C. Adolphsen et al. (International Study Group Collaboration), in International study group progress report on linear collider development, 2000, SLAC-R-559 and KEKREPORT-2000-7
  * [9] M. A. Diaz, S. F. King, D. A. Ross, Nucl. Phys. B 529, 23 (1998)
  * [10] S. Kiyoura, M. M. Nojiri, D. M. Pierce, Y. Yamada, Phys. Rev. D 58, 075002 (1998)
  * [11] T. Blank, W. Hollik, Eur. Phys .J. C 24, 619 (2002). arXiv:hep-ph/0011092
  * [12] J. Guasch, W. Hollik, J. Solà, Phys. Lett. B 510, 211 (2001)
  * [13] T. Hahn, C. Schappacher, Comp. Phys. Commun. 143, 54 (2002)
  * [14] T. Hahn, Nucl. Phys. Proc. Suppl. B 89, 231 (2000); G. J. van Olden-borgh, J. A. M. Vermaseren, Z. Phys. C 46, 425 (1990)
  * [15] A. Datta, A. Djouadi, J. L. Kneur, Phys.Lett. B 509, 299 (2001). arXiv:hep-ph/0101353
  * [16] K. Nakamura, J. Phys. G: Nucl. Part. Phys. 37, 075021 (2010)
  * [17] G. K. Egiiyan, M. Jurcisin, D. I. Kazakov, Mod. Phys. Lett. A 14, 601 (1999). arXiv:hep-ph/9807411
  * [18] J. F. Gunion, H. E. Haber, Nucl. Phys. B 272, 1 (1986)
  * [19] H. E. Haber, G.L. Kane, Phys. Rep. B 117, 75 (1985)
  * [20] M. Jurcisin, D. I. Kazakov, Mod., Phys. Lett. A 14, 671 (1999). arXiv:hep-ph/9902290
  * [21] S.Y. Choi et al., Eur. Phys. J. C 14, 535 (2000). arXiv:hep-ph/0002033; G. Moortgat-Pick et al., Eur. Phys. J. C 22, 563 (2001). arXiv:hep-ph/0108117
  * [22] T. Fritzsche, W. Hollik, Eur. Phys. J. C 24, 619 (2002). arXiv:hep-ph/0203159
  * [23] T. Fritzsche and W. Hollik, Nucl. Phys. Proc. Suppl. B 135, 102 (2004). arXiv:hep-ph/0407095; T. Fritzsche, in The International Conference on Linear Colliders, Paris, 2004. arXiv:hep-ph/0408307
  * [24] A. Dabelstein, Z. Phys. C 67, 495 (1995). arXiv:hep-ph/9409375
  * [25] A. Farzinnia, H. He, J. Ren, Phys. Lett. B 727, 141 (2013). arXiv:hep-ph/ 1308.0295
  * [26] W. Hollik, H. Rzehak, Eur. Phys. J. C 32, 127 (2003). arXiv:hep-ph/0305328
  * [27] A. Denner, Fortschr. Phys. 41, 307 (1993)
  * [28] G.’t Hooft and M. Veltman, Nucl. Phys. B 153, 365 (1979); L. Jing-Jing, M. Wen-Gan, Z. Ren-You, G. Lei, J. Yi, H. Liang, Phys. Rev. D 75, 053007 (2007)
  * [29] G W. T. Giele, E. W. N. Glover, Phys. Rev. D 46, 1980 (1992); W. T. Giele, E. W. Glover and D. A. Kosower, Nucl. Phys. B 403, 633 (1993)
  * [30] B. Y. Al-Negashi, T. A. El-Azim, I. A. M. Abdul-Magead, IREPHY 7, 259 (2013)
  * [31] T. Hahn, Comput. Phys. Commun. 168, 78 (2005). arXiv:hep-ph/0404043
  * [32] A. Denner, S. Dittmaier, Nucl. Phys B 398, 239 (1993); M. Böhm, S. Dittmaier, Nucl. Phys B 409, 3 (1993)
  * [33] W. Beenakker, A. Denner, Int. J. Mod. Phys. A 9, 4837 (1994)


# An analytical model of momentum availability for predicting large wind farm
power

Andrew Kirby1<EMAIL_ADDRESS>Thomas D. Dunstan2  Takafumi
Nishino1  1Department of Engineering Science, University of Oxford, Parks
Road, Oxford OX1 3PJ, UK 2Met Office, FitzRoy Road, Exeter EX1 3PB, UK

###### Abstract

Turbine-wake and farm-atmosphere interactions influence wind farm power
production. For large offshore farms, the farm-atmosphere interaction is
usually the more significant effect. This study proposes an analytical model
of the ‘momentum availability factor’ to predict farm-atmosphere interactions.
It models the effects of net advection, pressure gradient forcing and
turbulent entrainment. The model uses steady quasi-1D flow assumptions.
Turbulent entrainment is modelled by assuming self-similar vertical shear
stress profiles. We used the model with the ‘two-scale momentum theory’
(Nishino & Dunstan 2020, J. Fluid Mech. 894, A2) to predict the power output
of large finite-sized farms. The model compared well with existing results of
large-eddy simulations (LES) of finite wind farms in conventionally neutral
boundary layers. The model captured most of the effects of atmospheric
conditions (e.g. atmospheric boundary layer (ABL) height and free-atmosphere
stratification) on farm performance by only considering the undisturbed
vertical shear stress profile of the ABL. The model predicted the power of
staggered wind farms with a typical error of 5% or less. The developed model
provides a novel way of quickly predicting the power output of large wind
farms, including the farm blockage effects. This study also provides a new
general framework for modelling farm-scale flows. Future studies can use the
framework to better model large wind farms.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see Keyword PDF for the full list). Other
classifications will be added at the same time.

MSC Codes  (Optional) Please enter your MSC Codes here

##  1 Introduction

Wind energy is a key technology for the renewable energy transition. To meet
future energy demands, wind energy capacity will need to increase rapidly and
individual farms will likely become larger (Veers, 2022). A key aspect for
designing wind farms is predicting their power output. However, it is
difficult to model the aerodynamics of large wind farms because of the multi-
scale nature of flows involved (Porté-Agel et al., 2020).

Traditionally, there are two main approaches to predicting wind farm
performance at a low computational cost. For the first approach, semi-
analytical ‘wake’ models predict the velocity deficit in the wake behind a
turbine (e.g. Jensen, 1983; Bastankhah & Porté-Agel, 2014). To model an entire
wind farm, individual wakes are superposed using different techniques (e.g.
Katic et al., 1986; Zong & Porté-Agel, 2020). Wake models are commonly used to
optimise the layout of turbines in a farm. However, they do not consider the
atmospheric response to wind farms and thus they perform poorly for extended
wind farms (e.g. Stevens et al., 2016). The second approach uses ‘top-down’
models (e.g. Frandsen, 1992; Frandsen et al., 2006; Calaf et al., 2010). Top-
down models consider the response of an idealised atmospheric boundary layer
(ABL) to an infinitely large wind farm. They cannot however, predict the
impact of turbine layout on farm performance. Recent studies have coupled the
wake and top-down models together (e.g. Stevens et al., 2016; Starke et al.,
2021). This approach however, is still limited by the limitations of the
constituent models, e.g., idealised ABL profiles and wake superposition
methods.

Nishino & Dunstan (2020) developed the ‘two-scale momentum theory’ to predict
the power output of large wind farms. This splits the multi-scale problem of
wind farm aerodynamics into ‘internal’ turbine/array-scale and ‘external’
farm/atmospheric-scale sub-problems. The sub-problems are coupled together
using the conservation of momentum. As farms become larger, one of the grand
challenges facing the wind energy community is to understand the farm-
atmosphere interaction (Veers, 2019). Kirby et al. (2022) performed large-eddy
simulations (LES) of 50 different periodic turbine layouts. They introduced
the concepts of turbine-scale losses due to turbine-wake interactions and
farm-scale losses caused by the atmospheric response to the whole farm. For
large offshore farms, the farm-scale losses were found to be typically more
than twice as large as the turbine-scale losses. This highlights the
importance of modelling farm-scale flows to predict the performance of future
large farms.

Unlike most large wind farm models, the two-scale momentum theory does not
assume any specific profiles of the ABL (such as a logarithmic law), allowing
for the external sub-problem to be modelled in various manners. For the
external modelling, it is crucial to accurately predict how the momentum
available to the farm site increases due to the presence of the turbines,
i.e., the momentum availability. Recent studies (Patel et al., 2021; Legris et
al., 2023) used ‘twin’ numerical weather prediction (NWP) simulations to
calculate the momentum availability and thus predict how farm-scale losses
(including the so-called farm blockage losses) change with atmospheric
conditions. However, this approach would be too computationally expensive for
wind farm optimisation unless a sufficiently long set of twin NWP (such as
those reported by van Stratum et al. (2022)) has already been conducted for
candidate wind farm sites.

In the present study we propose a simple analytical model to predict the
momentum availability for large wind farms. This model, together with the two-
scale momentum theory, allows us to predict the power of a large wind farm
analytically. The model is derived using quasi-1D control volume analysis and
assuming self-similar vertical shear stress profiles. In section 2, we
summarise the two-scale momentum theory and the key wind farm parameters.
Section 3 presents the derivation of the new momentum availability factor
model. In section 4 we compare the predictions of farm power production with
existing finite wind farm LES. The model is discussed in section 5 and
concluding remarks are given in section 6.

##  2 Two-scale momentum theory

By considering the conservation of momentum for a control volume with and
without a wind farm present, Nishino & Dunstan (2020) derived the non-
dimensional farm momentum (NDFM) equation:

$C_{T}^{*}\frac{\lambda}{C_{f0}}\beta^{2}+\beta^{\gamma}=M$ (1)

where $\beta$ is the farm wind-speed reduction factor defined as $\beta\equiv
U_{F}/U_{F0}$ (with $U_{F}$ defined as the average wind speed in the nominal
farm-layer of height $H_{F}$, and $U_{F0}$ is the farm-layer-averaged speed
without the presence of the turbines); $C_{T}^{*}$ is the (farm-averaged)
‘internal’ turbine thrust coefficient defined as
$C_{T}^{*}\equiv\sum_{i=1}^{n}T_{i}/\frac{1}{2}\rho U_{F}^{2}nA$ (where
$T_{i}$ is the thrust of turbine $i$, $A$ is the rotor swept area and $n$ is
the number of turbines in the farm); $\lambda$ is the array density defined as
$\lambda\equiv nA/S_{F}$ (where $S_{F}$ is the wind farm area); $C_{f0}$ is
the natural friction coefficient of the surface defined as
$C_{f0}\equiv\langle\tau_{w0}\rangle/\frac{1}{2}\rho U_{F0}^{2}$ (where
$\tau_{w0}$ is the undisturbed bottom shear stress); $\gamma$ is the bottom
friction exponent defined as
$\gamma\equiv\log_{\beta}(\langle\tau_{w}\rangle/\langle\tau_{w0}\rangle)$
(assumed to be 2 in this study); and $M$ is the momentum availability factor
given by:

$M=\frac{X-C-\left[\frac{\partial p}{\partial
x_{F}}\right]-\frac{\partial[\rho U]}{\partial
t}}{X_{0}-C_{0}-\left[\frac{\partial p_{0}}{\partial
x_{F0}}\right]-\frac{\partial[\rho U_{0}]}{\partial t}}$ (2)

where $U$ is the velocity in the hub-height wind direction (i.e. streamwise
direction) $x_{F}$, $X$ represents the net streamwise momentum injection
through the top and side boundaries of the control volume (due to advection
and Reynolds stress), $C$ is the streamwise component of the Coriolis force
averaged over the control volume, $\partial p/\partial{x_{F}}$ is the pressure
gradient in the direction $x_{F}$ and the subscript $0$ refers to values
without the turbines present. In this study we use a fixed definition for the
farm-layer height of $H_{F}=2.5H_{hub}$ where $H_{hub}$ is the turbine hub-
height, following Kirby et al. (2022).

The theoretical framework Nishino & Dunstan (2020) used to derive equation 1
is known as the ‘two-scale momentum theory’. The left-hand side of equation 1
is expected to depend primarily on turbine-scale or ‘internal’ conditions.
This includes turbine layout and operating conditions, and the hub-height wind
speed and direction. The right-hand side of equation 1 is assumed to depend on
‘external’ farm-scale conditions.

The farm wind-speed reduction factor $\beta$ can be calculated using 1 for a
set of $C_{T}^{*}$, $\gamma$ and $M$. Using $\beta$, the average turbine power
coefficients can be calculated using

$C_{p}=\beta^{3}C_{p}^{*}$ (3)

where $C_{p}$ is the (farm-averaged) turbine power coefficient defined as
$C_{p}\equiv\sum_{i=1}^{n}P_{i}/\frac{1}{2}\rho U_{F0}^{3}nA$ ($P_{i}$ is
power of turbine $i$ in the farm) and $C_{p}^{*}$ is the (farm-averaged)
‘internal’ turbine power coefficient defined as
$C_{p}^{*}\equiv\sum_{i=1}^{n}P_{i}/\frac{1}{2}\rho U_{F}^{3}nA$.

##  3 Momentum availability factor model

The momentum availability factor $M$ describes the increase in momentum
supplied to the farm site due to the presence of the turbines (equation 4). We
introduce new variables $M_{F0}$ defined as the total net momentum flux into
the farm control volume without the turbines present, and $M_{F}$ as the total
net momentum flux with the turbines. $\Delta M_{F}$ is the change in momentum
flux due to the presence of the turbines, defined as $\Delta M_{F}\equiv
M_{F}-M_{F0}$.

$\Delta M_{F}$ can be decomposed into the contributions from different
physical mechanisms. For this study, we will consider only the contributions
from net momentum advection, pressure gradient forcing (PGF) and turbulent
entrainment (equation 4). Equation 4 provides a new framework for modelling
the impact of farm-scale flows.

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{M_{F}}{M_{F0}}=1+\frac{\Delta M_{F}}{M_{F0}}$ (4)
$\displaystyle\mbox{}=1+\frac{\Delta M_{F,Advection}}{M_{F0}}+\frac{\Delta
M_{F,PGF}}{M_{F0}}+\frac{\Delta M_{F,Entrainment}}{M_{F0}}$

In this study we propose simple analytical models for the components of
equation 4. For each component we consider a control volume of height $H_{F}$
around the farm (note that this control volume can be different from the
control volume used to derive equation 1). The height of the control volume
can be chosen arbitrarily but choosing $H_{F}$ allows terms to be linked to
farm wind-speed reduction factor $\beta$. We use a steady quasi-1D analysis
for the advection and PGF terms. The entrainment term is modelled by assuming
self-similar vertical shear stress profiles (above the top turbine-tip).
However, the modelling framework is not specific to the analytical models
proposed in the next sections. Equation 4 can be a starting point for future
studies using more sophisticated approaches for each component.

###  3.1 Net momentum advection

Figure 1 shows a rectangular wind farm control volume. $L$ is the length of
the farm in the hub-height wind direction and $W$ is the farm width. Using a
quasi-1D approach, we define the spanwise- and vertically-averaged wind speed
throughout the control volume as $U_{F0}\beta_{local}(x)$. In our notation
$x=0$ is the front of the farm and $x=L$ is the rear. Note that in figure 1a,
$\beta_{local}(0)$ is not exactly equal to 1 since the wind speed decelerates
upstream of the farm due to the farm blockage effect. Our proposed model of
net momentum advection can therefore capture the impact of farm blockage.

Figure 1: Control volume analysis for net momentum advection calculation a)
side view and b) front view.

By considering the conservation of mass of an elemental control volume (shown
by the bold box in figure 1), the mass flux out of the farm control volume at
position $x$ can be expressed as

$\displaystyle\dot{m}_{out}(x)=-\rho
H_{F}WU_{F0}\frac{\mathrm{d}\beta_{local}(x)}{\mathrm{d}x}.$ (5)

Note that this could be mass flux out of the top or sides of the farm control
volume. The momentum flux into the control volume at the front surface (i.e.
$x=0$) is given by

$\dot{m}_{in,front}U_{in,front}=\rho U_{F0}^{2}\beta_{local}(0)^{2}H_{F}W$ (6)

noting that a positive value represents a net inflow of momentum. The net
momentum flux through rear surface (i.e. $x=L$) is given by

$\dot{m}_{in,rear}U_{in,rear}=-\rho U_{F0}^{2}\beta_{local}(L)^{2}H_{F}W.$ (7)

The momentum flux through the top surface is the integral of the mass flux and
streamwise velocity at each position $x$, i.e.

$\displaystyle\dot{m}_{in,top}U_{in,top}=-\int_{0}^{L}\dot{m}_{out}(x)U_{F0}\beta_{local}(x)\,\mathrm{d}x.$
(8)

We substitute 5 into 8 and integrate to obtain

$\displaystyle\dot{m}_{in,top}U_{in,top}=\rho
U_{F0}^{2}H_{F}W\left[\frac{1}{2}\beta_{local}(L)^{2}-\frac{1}{2}\beta_{local}(0)^{2}\right].$
(9)

The total net momentum inflow into the control volume is given by the sum of
equations 6, 7 and 9, i.e.,

$\Delta M_{F,Advection}=\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left[\beta_{local}(0)^{2}-\beta_{local}(L)^{2}\right].$ (10)

This expression relates the change in momentum advection to the velocity at
the front and rear of the farm.

###  3.2 Pressure gradient forcing

LES results in the literature have shown that large wind farms can induce
additional pressure gradients across them (e.g., Allaerts & Meyers, 2017; Wu &
Porté-Agel, 2017; Lanzilao & Meyers, 2022). When farms induce additional
pressure gradients, a velocity reduction at the front of the farm is commonly
observed. This velocity reduction at the front of the farm is commonly known
as ‘wind-farm blockage’. Typically, a velocity increase and pressure reduction
at the rear of the farm also occurs. Figure 2 shows example streamwise
variations of $\beta_{local}(x)$ and pressure.

Figure 2: Example variations of a) local farm wind-speed reduction factor
$\beta_{local}(x)$ and b) pressure $p(x)$ with streamwise location.

In this section we propose a simple model for $\Delta M_{F,PGF}$. We apply
Bernoulli’s equation as an approximation to link the pressure increase to the
velocity reduction at the front of the farm. Firstly, we assume that, in the
upstream of the farm, the bottom friction force is balanced by the natural
PGF, so that the increase in pressure can be linked to the velocity reduction.
We also neglect changes in gravitational potential energy. It is assumed that
changes in pressure are uniform up to the control volume height $H_{F}$.

Bernoulli’s equation is applied on a streamline from far upstream to the front
of the farm. The increase in pressure force on the front surface of the
control volume is therefore given by:

$\Delta P_{front}H_{F}W=\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left[1-\beta_{local}(0)^{2}\right].$ (11)

We then assume that the pressure decrease at the rear of the farm is equal to
the pressure increase at the front, even though this is a strong assumption
(in reality, this depends on the ‘near-wake’ characteristics of the farm).
Therefore the change in pressure force at the rear surface is given by:

$\Delta P_{rear}H_{F}W=-\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left[1-\beta_{local}(0)^{2}\right].$ (12)

$\Delta M_{F,PGF}$ is then given by:

$\Delta M_{F,PGF}=(\Delta P_{front}-\Delta P_{rear})H_{F}W=\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left[2-2\beta_{local}(0)^{2}\right].$ (13)

Combing advection and PGF terms gives:

$\Delta M_{F,Advection}+\Delta M_{F,PGF}=\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left[2-\beta_{local}(0)^{2}-\beta_{local}(L)^{2}\right].$
(14)

Generally, the velocity at the front of the farm will be close to the
undisturbed value. The velocity at the rear of the farm will typically be
close to the farm-averaged value. As a simple first-order approach, we assume
that the velocity reduction at the front of the farm, $1-\beta_{local}(x)$, is
equal to the increase at the rear (relative to the farm-averaged wind speed),
$\beta_{local}(L)-\beta$. As such we can say that
$\beta_{local}(0)^{2}+\beta_{local}(L)^{2}\approx 1+\beta^{2}$. Applying this
approximation to equation 14:

$\Delta M_{F,Advection}+\Delta M_{F,PGF}=\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left(1-\beta^{2}\right).$ (15)

Dividing by the initial momentum supply $M_{F0}$ gives:

$\frac{\Delta M_{F,Advection}}{M_{F0}}+\frac{\Delta
M_{F,PGF}}{M_{F0}}=\frac{\frac{1}{2}\rho
U_{F0}^{2}H_{F}W\left(1-\beta^{2}\right)}{\tau_{w0}LW}=\frac{1}{C_{f0}}\frac{H_{F}}{L}\left(1-\beta^{2}\right).$
(16)

Although this is an approximation, equation 16 shows that when advection and
PGF terms are combined, the dependence on farm inlet and outlet velocities
disappears. Interestingly, this suggests that the impact of farm-scale flows
may not directly depend on the wind speed at the front of the farm.

###  3.3 Turbulent entrainment

Wind farms increase the turbulent mixing within the ABL. This increases the
momentum entrainment into the farm (Stevens & Meneveau, 2017). This mechanism
supplies momentum to the control volume through the shear stress at the top
surface. $\Delta M_{F,Entrainment}$ can be expressed as

$\frac{\Delta
M_{F,Entrainment}}{M_{F0}}=\frac{\tau_{t}LW-\tau_{t0}LW}{\tau_{w0}LW}=\frac{\tau_{t}-\tau_{t0}}{\tau_{w0}}$
(17)

where $\tau_{t}$ is the shear stress at the top of the control volume.

The ABL can be stratified with different profiles of potential temperature
which can change farm performance (Porté-Agel et al., 2020). In this study we
consider conventionally neutral boundary layers (CNBLs) because many wind farm
LES studies adopted ABLs with these profiles. CNBLs have self-similar vertical
shear stress profiles (Liu et al. (2021)). Therefore, if the surface stress
and the boundary layer height are known, then the shear stress at any height
can be determined. Note that the height is normalised by
$h=h_{0.05}/(1-0.05^{2/3})$ where $h_{0.05}$ is the height where the shear
stress is 5% of the surface value.

The self-similarity of shear stress profiles can be applied to large wind
farms. Abkar & Porté-Agel (2013) performed horizontally periodic LES of CNBLs
with and without turbines present. Figure 3a shows the vertical profiles of
the shear stress (note this is the stress in the hub-height wind direction,
$x_{F}$). Above the turbine top tip ($126.5m$ in Abkar & Porté-Agel (2013))
the vertical profiles have a similar shape. This suggests that the stress
profiles above the turbines is equivalent to an ‘empty’ CNBL with a higher
surface stress.

The black dotted lines in figure 3a show the stress profiles above the
turbines extrapolated to the surface (using a second order polynomial
regression). This corresponds to the total bottom stress, i.e., surface shear
stress and turbine thrust. When the stress profiles of Abkar & Porté-Agel
(2013) are normalised by the total bottom stress and new ABL height, they fall
onto the same curve (figure 3b). This shows that the CNBLs containing wind
farms also follow the same self-similar shear stress profile. Figure 3b shows
that wind farms increase the total bottom stress and boundary layer height of
CNBLs.

Figure 3: a) Vertical profiles of shear stress from horizontally periodic LES
with and without the farm present (Abkar & Porté-Agel, 2013) b) normalised
vertical profiles.

The self-similarity of stress profiles can be used to predict momentum
entrainment into large finite-sized wind farms. As a first order approach, we
consider the vertical shear stress profile horizontally-averaged across the
farm. Figure 4 shows a schematic of the shear stress profile with and without
the farm present. Here, $h_{0}$ is the ABL height without the farm present
averaged across the farm site. $h$ is the averaged ABL height with the farm
present. With the farm present, the shear stress is scaled by $M$ and the
heights scaled by $h/h_{0}$.

Figure 4: Schematic of vertical shear stress profiles with and without the
wind farm.

To determine $\tau_{t}$, we need to calculate the undisturbed shear stress at
a height of $H_{F}h_{0}/h$ (denoted as $\tau_{a}$). This is found by linearly
interpolating the stress between the top of the control volume and the
surface, i.e.,

$\tau_{a}=\tau_{w0}-(\tau_{w0}-\tau_{t0})\frac{h_{0}}{h}.$ (18)

Whilst the shear stress profile is not strictly linear, it is approximately
linear away from the top of the ABL. $\tau_{t}$ can then be expressed in terms
of $M$ and the undisturbed shear stress profile

$\displaystyle\tau_{t}$ $\displaystyle=$ $\displaystyle
M\tau_{a}=M\left(1-\frac{h_{0}}{h}\right)\tau_{w0}+M\frac{h_{0}}{h}\tau_{t0}.$
(19)

Equation 19 can be substituted into equation 17 to calculate $\Delta
M_{F,Advection}$,

$\displaystyle\frac{\Delta
M_{F,Entrainment}}{M_{F0}}=\frac{(\tau_{t}-\tau_{t0})LW}{\tau_{w0}LW}$
$\displaystyle=$ $\displaystyle
M+M\frac{h_{0}}{h}\left(\frac{\tau_{t0}}{\tau_{w0}}-1\right)-\frac{\tau_{t0}}{\tau_{w0}}$
(20)

The different sources of momentum increase can be summed to calculate $M$.
Subsituting equations 16 and 20 into 4 we obtain

$\displaystyle
M=1+\frac{1}{C_{f0}}\frac{H_{F}}{L}\left(1-\beta^{2}\right)+M+M\frac{h_{0}}{h}\left(\frac{\tau_{t0}}{\tau_{w0}}-1\right)-\frac{\tau_{t0}}{\tau_{w0}}$
(21)

which can be rearranged to give

$\displaystyle
M=\frac{1+\frac{1}{C_{f0}}\frac{H_{F}}{L}\left(1-\beta^{2}\right)-\frac{\tau_{t0}}{\tau_{w0}}}{\frac{h_{0}}{h}\left(1-\frac{\tau_{t0}}{\tau_{w0}}\right)}.$
(22)

Using equation 4 we can simply sum the different sources of momentum increase.
The result is a single equation to model farm-scale flows and their impact on
farm performance. Equation 22 is an approximate expression for $M$ as a
function of $\beta$ and $h/h_{0}$. Different models could be used to predict
$h/h_{0}$ in this formula for $M$. In the next section we present a simple
first order model for $h/h_{0}$ as a function of $\beta$.

###  3.4 Boundary layer height increase

In this section we present a simple model for increase in ABL height $h/h_{0}$
in response to a wind farm. We use a quasi-1D analysis of a boundary layer
flow over a farm (see figure 5). Here $\delta(x)$ is the ABL height at
steamwise position $x$. We assume there is no mass flux between the boundary
layer and free atmosphere above. The quasi-1D approach neglects the effect of
flow bypass around the sides of the farm.

Figure 5: Quasi-1D control volume analysis for ABL flow over a wind farm.

We consider the conservation of mass across a control volume encompassing a
section of the ABL (in bold in figure 5). The velocity averaged over the
height of the ABL at position $x$ is denoted $U_{A}(x)$. We assume that the
streamwise variation of $U_{A}(x)$ is the same as the streamwise variation of
$\beta_{local}(x)$. This is strictly not true, but could give a reasonable
first-order approximation of $U_{A}(x)$. Therefore
$U_{A}(x)=U_{A0}\beta_{local}(x)$ where $U_{A0}$ is the velocity averaged
throughout the undisturbed ABL. The mass flux in through the left hand side of
the control volume is given by:

$\displaystyle\dot{m}_{in}=\rho U_{A0}\beta_{local}(x)\delta(x)W.$ (23)

Taking the limit as $\Delta x$ approaches $0$, the mass flux out of the right
hand side is given by:

$\displaystyle\dot{m}_{out}$ $\displaystyle=$
$\displaystyle\rho\left[U_{A0}\beta_{local}(x)+U_{A0}\frac{\mathrm{d}\beta_{local}(x)}{\mathrm{d}x}\mathrm{d}x\right]\left[\delta(x)+\frac{\mathrm{d}\delta(x)}{\mathrm{d}x}\mathrm{d}x\right]W$
$\displaystyle\dot{m}_{out}$ $\displaystyle=$ $\displaystyle\rho
U_{A0}\beta_{local}(x)\delta(x)W+\rho
U_{A0}\beta_{local}(x)\frac{\mathrm{d}\delta(x)}{\mathrm{d}x}\mathrm{d}xW+$
(24) $\displaystyle\rho
U_{A0}\frac{\mathrm{d}\beta_{local}(x)}{\mathrm{d}x}\mathrm{d}x\delta(x)W$

whilst neglecting second order terms. From the mass conservation,
$\dot{m}_{in}=\dot{m}_{out}$,

$\displaystyle\rho
U_{A0}\beta_{local}(x)\frac{\mathrm{d}\delta(x)}{\mathrm{d}x}W$
$\displaystyle=$ $\displaystyle-\rho
U_{A0}\delta(x)\frac{\mathrm{d}\beta_{local}(x)}{\mathrm{d}x}W$
$\displaystyle\frac{\mathrm{d}\delta(x)}{\delta(x)}$ $\displaystyle=$
$\displaystyle-\frac{\mathrm{d}\beta_{local}(x)}{\beta_{local}(x)}$ (25)

Both sides of equation 25 can be integrated from far upstream to a given
position, resulting in

$\displaystyle\delta(x)=\frac{h_{0}}{\beta_{local}(x)}$ (26)

using the condition that at far upstream $\delta(x)=h_{0}$ and
$\beta_{local}(x)=1$. To find $h$, $\delta(x)$ is averaged between $x=0$ and
$x=L$ (as $h$ is the ABL height averaged across the wind farm site):

$\displaystyle\frac{h}{h_{0}}=\frac{1}{L}\int_{0}^{L}\frac{\mathrm{d}x}{\beta_{local}(x)}\approx\frac{1}{\beta}$
(27)

where the relationship between $\beta$ and $\beta_{local}(x)$ is given by
$\beta=\int_{0}^{L}\beta_{local}(x)\mathrm{d}x$. If $\beta_{local}(x)$ is
constant, the left hand side of equation 27 is equal to $1/\beta$, which we
assume to be generally true. In reality, the more $\beta_{local}(x)$ varies
the less accurate this assumption becomes. If $\beta_{local}(x)$ is close to
zero at any point the approximation in equation 27 becomes less accurate.

To check the validity of this approximation, we used data from finite wind
farm LES performed by Wu & Porté-Agel (2017). We take the hub-height wind
speed normalised by the inflow value as a proxy for $\beta_{local}(x)$. Note
that $H_{F}$ has been defined so that the hub-height wind speed is
approximately equal to the farm-layer-averaged speed (Kirby et al., 2022). For
the four farms simulated by Wu & Porté-Agel (2017), the maximum error of this
assumption was about 0.4%. Therefore the approximation is reasonable for
realistic profiles of $\beta_{local}(x)$.

| $\Gamma=1K/km$ | $\Gamma=5K/km$
---|---|---
| Staggered | Aligned | Staggered | Aligned
$1/\beta$ | 1.250 | 1.200 | 1.349 | 1.281
$\frac{1}{L}\int_{0}^{L}\frac{\mathrm{d}x}{\beta_{local}(x)}$ | 1.253 | 1.202 | 1.354 | 1.285
Percentage error | 0.244% | 0.160% | 0.414% | 0.248%
Table 1: Percentage errors of approximation in equation 27 using data from Wu
& Porté-Agel (2017)

Substituting equation 27 into equation 22, we get the following expression for
$M$

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{1+\frac{1}{C_{f0}}\frac{H_{F}}{L}\left(1-\beta^{2}\right)-\frac{\tau_{t0}}{\tau_{w0}}}{\beta\left(1-\frac{\tau_{t0}}{\tau_{w0}}\right)}.$
(28)

Equation 28 is a single analytical equation to predict the impact of farm-
scale flows on farm performance. It includes the effects of net advection, PGF
and turbulent entrainment. The only environmental parameter needed is the
undisturbed shear stress profile. Equation 28 can be used with the two-scale
momentum theory to predict farm power. In section 4 we compare the predictions
of farm power output with the results from finite wind farm LES results
reported in the literature.

###  3.5 Approximate expression for wind extractability factor $\zeta$

Patel et al. (2021) used twin NWP simulations to calculate $M$ for a realistic
wind farm site. They found, for most cases, that $M$ varied almost linearly
with the farm induction factor $(1-\beta)$. As such, $M$ can be modelled as

$\displaystyle M=1+\zeta(1-\beta)$ (29)

where $\zeta$ is the ‘wind extractability factor’ (which is identical to what
Nishino & Dunstan (2020) originally proposed as ‘momentum response factor’).
Kirby et al. (2022) reported that $\zeta$ varied with atmospheric conditions
and decreased exponentially with wind farm size. The analytical model of $M$
presented earlier can be used to derive an approximate expression for $\zeta$,
as the equation 28 can be expressed as

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{1}{\beta}+\frac{\frac{1}{C_{f0}}\frac{H_{F}}{L}}{1-\frac{\tau_{t0}}{\tau_{w0}}}\frac{(1-\beta)(1+\beta)}{\beta}.$
(30)

This can than be expressed in terms of farm induction factor $(1-\beta)$,
e.g.,

$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{1}{1-(1-\beta)}+\frac{\frac{1}{C_{f0}}\frac{H_{F}}{L}}{1-\frac{\tau_{t0}}{\tau_{w0}}}\frac{2(1-\beta)-(1-\beta)^{2}}{1-(1-\beta)}.$
(31)

Figure 6: Linear approximations for terms a) $\frac{1}{1-(1-\beta)}$ and b)
$\frac{2(1-\beta)-(1-\beta)^{2}}{1-(1-\beta)}$ in $M$ expression (equation
31).

The functions of $(1-\beta)$ in equation 31 are approximately linear for a
realistic range of $(1-\beta)$ (see figure 6). The analytical model of $M$
(equation 28) developed in this study therefore predicts the quasi-linear
momentum response seen in NWP simulations. We use linear interpolation for
$(1-\beta)$ between $0$ and $0.2$ (a realistic range) to find linear
approximations to the functions of $(1-\beta)$. As such $M$ can be
approximated by:

$\displaystyle M_{approx}$ $\displaystyle=$ $\displaystyle
1+1.18(1-\beta)+\frac{\frac{2.18}{C_{f0}}\frac{H_{F}}{L}}{1-\frac{\tau_{t0}}{\tau_{w0}}}(1-\beta).$
(32)

Equation 32 can be used to derive an approximate expression for $\zeta$, i.e.,

$\displaystyle\zeta_{approx}$ $\displaystyle=$ $\displaystyle
1.18+\frac{\frac{2.18}{C_{f0}}\frac{H_{F}}{L}}{1-\frac{\tau_{t0}}{\tau_{w0}}}$
(33)

using the expression $\zeta=(M-1)/(1-\beta)$. Equation 33 predicts the
exponential decrease of $\zeta$ with farm size reported by Kirby et al.
(2022). The discrepancy between $\zeta_{approx}$ and the ‘true’ $\zeta$
(obtained without the linear approximations) is a function only of
$\zeta_{approx}$ and $(1-\beta)$. Figure 7a shows the percentage error of
$\zeta_{approx}$ for a realistic range of $\zeta_{approx}$ and $(1-\beta)$.
The maximum error is approximately 10% and this occurs for low values of
$(1-\beta)$. However, for low values of $(1-\beta)$, the value of $M_{approx}$
is relatively insensitive to the value of $\zeta_{approx}$. Figure 7b shows
the percentage error of $M_{approx}$ which is typically less than 5% for
realistic farms.

Figure 7: Percentage error of a) $\zeta_{approx}$ and b) $M_{approx}$.

##  4 Comparison with finite wind farm LES

The two-scale momentum theory presented in section 2 can be used to predict
wind farm power. For arrays of actuator-discs (or aerodynamically ideal
turbines operating below the rated wind speed), the turbine power coefficient
can be calculated using

$C_{p}=\beta^{3}{C_{T}^{*}}^{\frac{3}{2}}{C_{T}^{\prime}}^{-\frac{1}{2}}$ (34)

where $C_{T}^{\prime}\equiv T_{i}/\frac{1}{2}\rho U_{T,i}^{2}A$ is a turbine
resistance coefficient describing the turbine operating conditions ($U_{T,i}$
is the streamwise velocity averaged over the rotor swept area of turbine $i$).
$C_{T}^{*}$ generally depends on the turbine layout, but its upper limit may
be predicted by using an analogy to the classical actuator-disc theory
(Nishino, 2016), as

$C_{T}^{*}=\frac{16C_{T}^{\prime}}{(4+C_{T}^{\prime})^{2}}.$ (35)

Since this model neglects turbine-scale flow interactions, it can be used to
predict an upper limit to farm performance (Kirby et al., 2022). Using this
model with $\beta=1$, the power coefficient of an isolated turbine
$C_{p,Betz}$ can be obtained as

$C_{p,Betz}=\frac{64C_{T}^{\prime}}{(4+C_{T}^{\prime})^{3}}.$ (36)

It is important to note that equations 34 and 35 are only strictly valid when
$U_{T,i}$ is the same for all turbines. However, in this study we only
consider wind farms in which $C_{T}^{\prime}$ is the same for all turbines,
and in this case, equations 34 and 35 can be applied in an approximate manner
when $U_{T,i}$ varies throughout the farm.

We used the new model of $M$ (equation 28) with the two-scale momentum theory
to predict wind farm power. See appendix A for a brief summary of how to use
the analytical model to predict farm power. We compared the predictions
against finite wind farm LES results from Wu & Porté-Agel (2017), Allaerts &
Meyers (2017) and Lanzilao & Meyers (2022). We selected these studies because
they simulated large wind farms (longer than 15km in the streamwise
direction). They published the undisturbed vertical shear stress profile which
is a required input to the model of $M$. These studies were also selected
because the results are normalised by the power of an isolated turbine rather
than the front row turbine power. The 3 studies performed LES of 4 farms with
staggered turbine layouts and 6 farms with aligned layouts.

We used the analytical model to predict the farm-averaged power normalised by
the isolated turbine power, i.e., $C_{p}/C_{p,Betz}$. Allaerts & Meyers (2017)
and Lanzilao & Meyers (2022) used an actuator-disc no rotation model for the
turbines with $C_{T}^{\prime}=1.33$. For these studies we used $C_{T}^{*}$
from the actuator-disc theory (equation 35) as an upper limit for $C_{T}^{*}$
in 1 and then used the expression $C_{p}/C_{p,Betz}=\beta^{3}$. For this
$C_{T}^{\prime}$ value, the actuator-disc theory gives $C_{T}^{*}=0.75$. Wu &
Porté-Agel (2017) used an actuator-disc with rotation implementation for the
turbines. We used the published thrust coefficient curves (Wu & Porté-Agel,
2015) as an upper limit to $C_{T}^{*}$. Kirby et al. (2022) reported that
turbine-wake interactions reduced $C_{T}^{*}$ by up to about 20%. For the
lower limit of farm performance, we used a 20% reduced $C_{T}^{*}$ value in
equation 1 and then the expression $C_{p}/C_{p,Betz}=\beta^{3}0.8^{1.5}$
(assuming that $C_{T}^{\prime}$ is unaffected by the turbine layout). Note
that this expression comes from equation 34 using a $C_{T}^{*}$ value which is
80% of the upper limit.

Lanzilao & Meyers (2022) normalised the turbine power by the power of an
imaginary turbine row 10km upstream of the farm. However, small reductions in
wind speed were observed 10km upstream for the stratified boundary layer. We
instead normalised farm powers from this study by the imaginary upstream power
for the neutrally stratified case. This is because the neutrally stratified
case had a much smaller reduction in wind speed upstream of the farm. The
imaginary upstream turbine power of the staggered and aligned farms in the
stratified boundary layer were 6.24MW and 6.17MW respectively (Lanzilao,
2023). For the staggered farm in the neutrally stratified case it was 6.61MW
(Lanzilao, 2023).

A comparison between the model predictions and LES is shown in figure 8.
Generally, the staggered LES results are close to the upper limits predicted
by the two-scale model. The aligned layouts are close to the lower limit of
the predictions, except for the results from $\Gamma 1$ in Wu & Porté-Agel
(2017) where the LES power is higher than the model predictions. The
analytical model captures the impact of atmospheric conditions on farm
performance. The model predicts well the decrease in farm performance with
decreasing inversion layer height observed by Allaerts & Meyers (2017). For
the results of Lanzilao & Meyers (2022), the staggered wind farms had very
similar performances in the neutral and conventionally neutral boundary
layers. The two-scale predictions are exactly the same for these two cases
because the shear stress profiles were identical.

Figure 8: Comparison of farm-average performance predicted by the two-scale
model and results from LES.

The model gave accurate predictions of farm power for staggered turbine
layouts. We made direct comparisons with staggered turbine layouts because
they have a $C_{T}^{*}$ value close to actuator-disc predictions. For the
aligned turbine layouts, the exact value of $C_{T}^{*}$ depends on turbine
spacing therefore a direct comparison is not currently possible. Table 2 shows
the percentage error in farm power for the four staggered farms. When
normalised by the farm power, the error was typically around 5% or less.
Normalising by the isolated turbine power gave errors of around 3% or less.

|  Wu et. al. 2017 |  Lanzilao et. al. 2022
---|---|---
| $\Gamma=1K/km$ | $\Gamma=5K/km$ | CS | NS
$\lambda$ | 0.0160 | 0.0160 | 0.0314 | 0.0314
$C_{f0}$ | 0.00557 | 0.00501 | 0.00314 | 0.00314
$H_{CV}/L$ | 0.00875 | 0.00875 | 0.0188 | 0.0188
$\tau_{t0}/\tau_{w0}$ | 0.571 | 0.391 | 0.475 | 0.475
$\beta_{model}$ | 0.847 | 0.816 | 0.821 | 0.821
$M_{model}$ | 2.41 | 2.40 | 5.74 | 5.74
$C_{p,model}/C_{p,Betz}$ | 0.607 | 0.543 | 0.554 | 0.544
$C_{p,LES}/C_{p,Betz}$ | 0.639 | 0.521 | 0.585 | 0.570
$(C_{p,model}-C_{p,LES})/C_{p,LES}$ | -5.00% | 4.33% | -5.39% | -2.87%
$(C_{p,model}-C_{p,LES})/C_{p,Betz}$ | -3.20% | 2.25% | -3.15% | -1.63%
Table 2: Model parameters and percentage errors in predicting the average
power of staggered wind farm LES.

##  5 Discussion

This study proposes a novel framework for modelling farm-scale flows. Equation
4 is a general expression for the farm momentum increase as the sum of
contributions from different mechanisms. This allows each mechanism to be
modelled separately. The contributions are summed to model net momentum
advection (including farm blockage), PGF and entrainment. Equation 4 could
provide a framework for combining different farm-scale flow models in the
future.

We have proposed simple analytical models for each mechanism increasing farm
momentum supply. The models make quasi-1D and steady flow assumptions. The
simple momentum availability factor model can be coupled with the two-scale
momentum theory to predict farm power. The predictions for staggered farms had
a typical error of around 5% or less. The model predictions are for the farm-
averaged power rather than just the power in the fully developed region.
Therefore, the model includes the effect of the development region and the
farm size. Overall, this suggests that wind farm performance in complex
atmospheric conditions can be largely explained using a simple quasi-1D
approach.

The only environmental parameter required as an input to the model is the
undisturbed vertical shear stress profile of the ABL. The model seems to
capture most of the impact of free-atmosphere stratification and boundary
layer height. This suggests that these factors impact farm performance largely
by changing the vertical profile of shear stress in the ABL. As far as the
authors are aware, this is a new way of understanding the role of atmospheric
conditions on farm efficiency. The proposed model also captures the effect of
farm blockage. This is a physics-based approach without tuning or empirical
coefficients.

Our model derivation suggests the effects of farm blockage and increased PGF
could (at least partly) counteract. A strong free-atmosphere stratification
reduces the wind speed in front of the farm (Wu & Porté-Agel, 2017; Lanzilao &
Meyers, 2022). This reduces the net advection of momentum into the farm. When
this occurs, an additional pressure gradient is induced across the farm. This
increases the momentum supply to the farm. In our model, the combined impact
of the two effects is independent of the wind speed at the front of the farm
(equation 16). This suggests the negative farm blockage effect could be
counteracted by the increased PGF. For cases NS and CS in Lanzilao & Meyers
(2022), the atmospheric conditions are kept constant except for the free-
atmosphere stratification. The CNBL showed a much greater wind farm blockage
yet both cases had a similar farm-averaged power. Hence, the interaction
between farm blockage and PGF seems important and requires further
investigation in the future.

Our model can provide physical insight into how farm efficiency changes with
atmospheric conditions. Figure 9 shows the predicted values of $M$ and $\beta$
for the LES studies considered. The model predicts an approximately linear
relationship between $M$ and $\beta$ as observed by Patel et al. (2021). Each
atmospheric condition has a different gradient which is equal to the wind
extractability factor $\zeta$. The LES results of Wu & Porté-Agel (2017) had a
much lower value of $\zeta$ because of the high bottom surface roughness
length of 0.05m (note, however, that this does not mean a higher surface
roughness leads to a lower farm power, because a higher surface roughness
gives a lower effective array density $\lambda/C_{f0}$, which usually leads to
a higher farm power (Nishino & Dunstan, 2020)). For the results of Allaerts &
Meyers (2017), S1 had the highest initial inversion layer and S4 the lowest.
Figure 9 shows that the momentum increase (or the value of $M$) is relatively
insensitive to the initial inversion layer height. However, a lower inversion
layer requires a greater wind speed reduction to increase the momentum supply
to a certain level (which depends on the ‘internal’ conditions, such as the
turbine layout and the value of $C_{T}^{\prime}$). The lower wind speed
results in S4 having a lower farm efficiency than S1 (see figure 8).

Figure 9: Predicted values of momentum availability factor ($M$) and farm
wind-speed reduction factor ($\beta$).

A taller ABL can increase the farm momentum supply with a smaller wind speed
reduction. To illustrate this, we consider the momentum response of two ABLs;
a low ABL height (with $H_{F}/h_{0}=0.75$) and a high ABL height (with
$H_{F}/h_{0}=0.25$). For both ABLs we use a typical value for large offshore
farms of $\frac{1}{C_{f0}}\frac{H_{F}}{L}=7.5$ and set $M=4$ (which physically
means that the ‘internal’ conditions of the farm are adjusted to acheive
$M=4$). The value of $\beta$ can be calculated using equation 28 for both
scenarios. The new and undisturbed shear stress profiles are shown in figure
10. The low ABL has a smaller entrainment of momentum at the top of the
control volume. The top of the control volume is closer to the top of the
boundary layer where the velocity shear is lower. Therefore the low ABL is
less efficient at entraining momentum into the farm. With less momentum
entrainment, more momentum has to be supplied through advection and PGF in
order to acheive the set level of $M$ (see figure 11). These mechanisms
increase the momentum supply by reducing the farm-averaged wind speed
(equation 16). In taller ABLs, more momentum is supplied to the farm through
entrainment, and the advection and PGF supply less momentum. Taller ABLs can
therefore sustain higher wind speeds in the farm and thus have a higher farm
efficiency.

Figure 10: Contribution of entrainment to farm momentum availability factor
$M$ for a) a low initial ABL height and b) a high initial ABL height.

Figure 11: Decomposition of farm momentum availability factor $M$ into
entrainment, and advection and PGF terms for a) a low initial ABL height and
b) a high initial ABL height. Note that the heights of green bars in this
figure correspond to the width of green bars in figure 10.

This study only presents a preliminary comparison with wind farm LES. To fully
validate the proposed model, a comparison with a larger set of LES would be
required. A limitation of the proposed model is that the quasi-1D analysis
captures only the first-order effects. However, the model still gave a close
agreement with LES results. More detailed flow physics were not considered
(e.g. gravity wave excitation) but could be included in future studies. Future
work will also aim to predict $C_{T}^{*}$ for different turbine layouts and
operating conditions. As an example, Legris et al. (2023) used a wake model to
predict $C_{T}^{*}$ for different farm designs. Future studies will also focus
on the farm-atmosphere interaction under more realistic atmospheric
conditions.

It should also be noted that, in reality, the assumption of scale separation
is not expected to be strictly valid. Therefore, the entrainment component of
the model may need some parameters that depend on the internal conditions
(e.g. turbine layout). However, given the close agreement of the current model
and LES results, such an inter-scale effect does not seem significant. This
means that the scale separation assumption introduced by Nishino & Dunstan
(2020) is approximately valid for large wind farms, at least for the
prediction of farm-average power.

##  6 Conclusions

In this study, we proposed an analytical model of the momentum availability
factor for large wind farms. This is a simple model to predict the impact of
farm-scale flows on farm performance. This study considered changes in net
advection, PGF and turbulent entrainment as the factors contributing to the
momentum availability. We derived the model using a steady quasi-1D control
volume analysis and assuming self-similar vertical shear stress profiles. The
only environmental parameters required as input to the model are the
undisturbed vertical shear stress profiles. We used the new model with the
two-scale momentum theory (Nishino & Dunstan, 2020) to predict large farm
power production. The model compared well against existing LES of finite wind
farms in CNBLs. The model captured the impact of atmospheric conditions (i.e.
boundary layer height and free-atmosphere stratification) and farm size on
farm performance. A direct comparison with LES of staggered farms showed a
typical error of 5% or less.

This study also provides a new general framework for modelling farm-scale
flows. The momentum availability factor is expressed as the sum of
contributions from different physical mechanisms. As such, the different
mechanisms can be modelled separately. We used first-order analytical models
for the different components. Despite the simplicity, the model predicts farm
power with close agreement to finite wind farm LES. Only one equation needs to
be solved analytically to calculate the farm wind-speed reduction factor
$\beta$ and thus predict farm performance. It also provides a physics-based
method to predict farm blockage effects and the impact of atmospheric
stability. Future studies could use more advanced models for the different
components of the momentum availability factor. This could provide more
accurate predictions of wind farm power production.

[Acknowledgements]The first author (AK) acknowledges the NERC-Oxford Doctoral
Training Partnership in Environmental Research (NE/S007474/1) for funding and
training. We also thank Luca Lanzilao, Prof. Johan Meyers and Prof. Richard
Stevens for providing additional data from their studies.

[Funding]This work was supported by the Natural Environmental Research Council
(NERC) award NE/S007474/1.

[Declaration of Interests]The authors report no conflict of interest.

[Author ORCID]A. Kirby, https://orcid.org/0000-0001-8389-1619; T. Nishino,
https://orcid.org/0000-0001-6306-7702.

[Author contributions]A.K. and T.N. derived the analytical model. A.K. made a
comparison of the model and existing LES results in the literature, and A.K.
wrote the paper with corrections from T.N. and T.D.D.

##  Appendix A Model formulation for predicting wind farm power

In this section we present a summary of how to predict farm power using the
proposed model. Firstly, equation 37 is solved for farm wind-speed reduction
factor $\beta$:

$C_{T}^{*}\frac{\lambda}{C_{f0}}\beta^{2}+\beta^{\gamma}=\frac{1+\frac{1}{C_{f0}}\frac{H_{F}}{L}\left(1-\beta^{2}\right)-\frac{\tau_{t0}}{\tau_{w0}}}{\beta\left(1-\frac{\tau_{t0}}{\tau_{w0}}\right)}$
(37)

where $L$ is the length of the farm in the hub-height wind
direction,$\tau_{t0}$ is the undisturbed shear stress at the top of the farm
layer (of height $H_{F}$) and we assume $\gamma=2$. Definitions of all other
terms can be found in section 2. Using the calculated value of $\beta$, the
farm-averaged power coefficient $C_{p}$ can be calculated as:

$C_{p}=\beta^{3}{C_{T}^{*}}^{\frac{3}{2}}{C_{T}^{\prime}}^{-\frac{1}{2}}.$
(38)

The definitions of $C_{p}$, $C_{T}^{*}$ and $C_{T}^{\prime}$ are summarised in
table 3.

Symbol | Name | Formula |  Note
---|---|---|---
$C_{p}$ | Average turbine power coefficient | $\frac{\sum_{i=1}^{n}P_{i}}{\frac{1}{2}\rho U_{F0}^{3}nA}$ |  $U_{F0}$ is the undisturbed farm-layer-averaged wind speed.
$C_{p}^{*}$ | Average ‘internal’ turbine power coefficient | $\frac{\sum_{i=1}^{n}P_{i}}{\frac{1}{2}\rho U_{F}^{3}nA}$ |  $U_{F}$ is the farm-layer-averaged wind speed with the turbines present.
$C_{T}^{\prime}$ | Turbine resistance coefficient | $\frac{T_{i}}{\frac{1}{2}\rho U_{T,i}^{2}A}$ |  $U_{T}$ is the wind speed averaged over the rotor swept area.
$C_{T}^{*}$ | Average ‘internal’ turbine thrust coefficient | $\frac{\sum_{i=1}^{n}T_{i}}{\frac{1}{2}\rho U_{F}^{2}nA}$ |  Determined by turbine operating conditions and layout.
Table 3: Summary of turbine power and thrust coefficients.

## References

  * Abkar & Porté-Agel (2013) Abkar, M. & Porté-Agel, F. 2013 The effect of free-atmosphere stratification on boundary-layer flow and power output from very large wind farms. Energies 6, 2338–2361.
  * Allaerts & Meyers (2017) Allaerts, D. & Meyers, J. 2017 Boundary-layer development and gravity waves in conventionally neutral wind farms. Journal of Fluid Mechanics 814, 95–130.
  * Bastankhah & Porté-Agel (2014) Bastankhah, M. & Porté-Agel, F. 2014 A new analytical model for wind-turbine wakes. Renewable Energy 70, 116–123.
  * Calaf et al. (2010) Calaf, M., Meneveau, C. & Meyers, J. 2010 Large eddy simulation study of fully developed wind-turbine array boundary layers. Physics of Fluids 22, 1–16.
  * Frandsen (1992) Frandsen, S. 1992 On the wind speed reduction in the center of large clusters of wind turbines. Journal of Wind Engineering and Industrial Aerodynamics 39, 251–265.
  * Frandsen et al. (2006) Frandsen, S., Barthelmie, R., Pryor, S., Rathmann, O., Larsen, S., Højstrup, J. & Thøgersen, M. 2006 Analytical modelling of wind speed deficit in large offshore wind farms. Wind Energy 9, 39–53.
  * Jensen (1983) Jensen, N. O. 1983 A note on wind generator interaction. Risø-M-2411 Risø National Laboratory Roskilde .
  * Katic et al. (1986) Katic, I., Hojstrup, J. & Jensen, N. O. 1986 A simple model for cluster efficiency. Proceedings of the European wind energy association conference and exhibition, Rome, Italy pp. 407–409.
  * Kirby et al. (2022) Kirby, A., Nishino, T. & Dunstan, T. D. 2022 Two-scale interaction of wake and blockage effects in large wind farms. Journal of Fluid Mechanics 953, A39.
  * Lanzilao (2023) Lanzilao, L. 2023 personal communication.
  * Lanzilao & Meyers (2022) Lanzilao, L. & Meyers, J. 2022 Effects of self-induced gravity waves on finite wind-farm operations using a large-eddy simulation framework. Journal of Physics: Conference Series 2265, 22043.
  * Legris et al. (2023) Legris, L., Pahus, M. L., Nishino, T. & Perez-Campos, E. 2023 Prediction and mitigation of wind farm blockage losses considering mesoscale atmospheric response. Energies 16, 386\.
  * Liu et al. (2021) Liu, L., Gadde, S. N. & Stevens, R. J. A. M. 2021 Universal wind profile for conventionally neutral atmospheric boundary layers. Phys. Rev. Lett. 126, 104502.
  * Nishino (2016) Nishino, T. 2016 Two-scale momentum theory for very large wind farms. Journal of Physics: Conference Series 753, 032054\.
  * Nishino & Dunstan (2020) Nishino, T. & Dunstan, T. D. 2020 Two-scale momentum theory for time-dependent modelling of large wind farms. Journal of Fluid Mechanics 894, A2.
  * Patel et al. (2021) Patel, K., Dunstan, T. D. & Nishino, T. 2021 Time-dependent upper limits to the performance of large wind farms due to mesoscale atmospheric response. Energies 14, 6437.
  * Porté-Agel et al. (2020) Porté-Agel, F., Bastankhah, M. & Shamsoddin, S. 2020 Wind-turbine and wind-farm flows: A review. Boundary-Layer Meteorology 174, 1–59.
  * Starke et al. (2021) Starke, G. M., Meneveau, C., King, J. R. & Gayme, D. F. 2021 The area localized coupled model for analytical mean flow prediction in arbitrary wind farm geometries. Journal of Renewable and Sustainable Energy 13.
  * Stevens et al. (2016) Stevens, R. J. A. M., Gayme, D. F. & Meneveau, C. 2016 Generalized coupled wake boundary layer model: applications and comparisons with field and les data for two wind farms. Wind Energy 19, 2023–2040.
  * Stevens & Meneveau (2017) Stevens, R. J. A. M. & Meneveau, C. 2017 Flow structure and turbulence in wind farms. Annual Review of Fluid Mechanics 49, 311–339.
  * van Stratum et al. (2022) van Stratum, B., Theeuwes, N., Barkmeijer, J., van Ulft, B. & Wijnant, I. 2022 A one-year-long evaluation of a wind-farm parameterization in harmonie-arome. Journal of Advances in Modeling Earth Systems 14, e2021MS002947.
  * Veers (2019) Veers, P. et al. 2019 Grand challenges in the science of wind energy. Science 366, eaau2027.
  * Veers (2022) Veers, P. et al. 2022 Grand challenges: wind energy research needs for a global energy transition. Wind Energy Science 7, 2491–2496.
  * Wu & Porté-Agel (2017) Wu, K. L. & Porté-Agel, F. 2017 Flow adjustment inside and around large finite-size wind farms. Energies 10, 2164\.
  * Wu & Porté-Agel (2015) Wu, Y. T. & Porté-Agel, F. 2015 Modeling turbine wakes and power losses within a wind farm using les: An application to the horns rev offshore wind farm. Renewable Energy 75, 945–955.
  * Zong & Porté-Agel (2020) Zong, H. & Porté-Agel, F. 2020 A momentum-conserving wake superposition method for wind farm power prediction. Journal of Fluid Mechanics 889, A8.