# Statistical Estimates of the Binary Properties of Rotational Variables

Anya Phillips ID 1 and C. S. Kochanek1,2
1Department of Astronomy, The Ohio State University, 140 West 18th Avenue,
Columbus, OH, 43210, USA
2Center for Cosmology and Astroparticle Physics, The Ohio State University,
191 W. Woodruff Avenue, Columbus, OH, 43210, USA E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a model to estimate the average primary masses, companion mass
ranges, the inclination limit for recognizing a rotational variable, and the
primary mass spreads for populations of binary stars. The model fits a
population’s binary mass function distribution and allows for a probability
that some mass functions are incorrectly estimated. Using tests with synthetic
data, we assess the model’s sensitivity to each parameter, finding that we are
most sensitive to the average primary mass and the minimum companion mass,
with less sensitivity to the inclination limit and little to no sensitivity to
the primary mass spread. We apply the model to five populations of binary
spotted rotational variables identified in ASAS-SN, computing their binary
mass functions using RV data from APOGEE. Their average primary mass estimates
are consistent with our expectations based on their CMD locations ($\sim
0.75M_{\odot}$ for lower main sequence primaries and $\sim 0.9$–$1.2M_{\odot}$
for RS CVn and sub-subgiants). Their companion mass range estimates allow
companion masses down to $M_{2}/M_{1}\simeq 0.1$, although the main sequence
population may have a higher minimum mass fraction ($\sim 0.4$). We see weak
evidence of an inclination limit $\gtrsim 50^{\circ}$ for the main sequence
and sub-subgiant groups and no evidence of an inclination limit in the other
groups. No groups show strong evidence for a preferred primary mass spread. We
conclude by demonstrating that the approach will provide significantly better
estimates of the primary mass and the minimum mass ratio and reasonable
sensitivity to the inclination limit with 10 times as many systems.

###### keywords:

stars: binaries: close – stars: rotation – stars: starspots – stars:
statistics

††pubyear: 2024††pagerange: Statistical Estimates of the Binary Properties of
Rotational Variables–Statistical Estimates of the Binary Properties of
Rotational Variables

##  1 Introduction

Population-level properties of binary stars are important to understanding
their formation and evolution. At high masses (primaries $\gtrsim 3M_{\odot}$;
Moe & Di Stefano 2017), these systems can evolve to produce X-ray binaries
(Verbunt, 1993), pulsars Lorimer (2008), and Type Ia supernovae (Wang & Han,
2012), among other extreme astrophysical phenomena. Across stellar mass,
however, measuring the fundamental properties of binary systems provides
insight into binary star formation. An example is the distribution of binary
mass ratios. Observing companion masses distributed as the stellar initial
mass function (IMF) implies that each companion has evolved independently of
its primary (Tout, 1991; McDonald & Clarke, 1995). On the other hand,
correlated companion and primary masses imply that the two stars coevolved
through pre-main sequence fragmentation, fission, or mass transfer through
Roche Lobe overflow (Bonnell & Bastien, 1992; Kroupa, 1995a, b; Clarke, 1996;
Kratter & Matzner, 2006; Kouwenhoven et al., 2009; Marks & Kroupa, 2011; Bate,
2012; Moe & Di Stefano, 2017). Another point of interest is the relation
between multiplicity or close binary fraction on stellar mass (e.g., Duquennoy
& Mayor 1991; Duchêne & Kraus 2013; Mazzola et al. 2020).

Measuring binary properties also provides an observational foundation for
binary star evolution. Numerous studies have, for example, examined the
physics of tidal circularization (Koch & Hrivnak, 1981; Duquennoy & Mayor,
1991; Verbunt & Phinney, 1995; Meibom & Mathieu, 2005; Van Eylen et al., 2016;
Bashi et al., 2023), tidal synchronization (e.g. Lurie et al. 2017; Leiner et
al. 2022), and their dependence on stellar mass, binary mass ratio,
eccentricity, and orbital period. Finally, measured binary properties inform
binary population synthesis models, as the rates of different channels of
binary evolution depend on the assumed binary statistics (see Han et al. 2020
for a recent review).

Unresolved main sequence binaries can be identified photometrically, sitting
above the main sequence in luminosity, especially for unevolved primaries
($T_{\rm{eff}}<5250\ K$). Simonian et al. (2019) use main sequence rotational
variables in the Kepler to show that rapid rotation (rotation periods $<7$
days) occurs primarily in photometric binaries. For these systems, the
apparent elevation above the main sequence allows an estimate of the binary
mass ratio, but this method is best for detecting near equal mass binaries
because of the steepness of the mass-luminosity relation on the main sequence.

Eclipsing double lined spectroscopic binaries (SB2s) provide the most
information about a system’s fundamental binary properties, allowing
measurement of both companion masses and the orbit’s inclination (e.g., Torres
et al. 2010). Non-eclipsing SB2s allow measurement of the mass ratios of
individual systems, but their inclinations, and hence, also their masses, are
still unknown. Identifying SB2s is straightforward at periods $\lesssim 5$
years, but at any orbital period it is also possible to spectroscopically
disentangle binaries by fitting a binary spectral model rather than a single-
star model, recovering the mass ratio and atmospheric composition of both
companions (e.g., El-Badry et al. 2018). Binary samples with spectroscopic
signals from both stars are also biased toward equal mass companions because
high mass ratio systems typically have a high flux ratio and a secondary that
is difficult to detect.

Eclipsing single lined spectroscopic binaries (SB1s) require knowledge of the
primary mass in order to recover the secondary mass. For non-eclipsing SB1s,
we can only observe the orbital period $P$ and velocity semiamplitude $K$ to
compute the binary mass function,

$f=\frac{PK^{3}}{2\pi
G}=\frac{M_{2}^{3}\sin^{3}{i}}{(M_{1}+M_{2})^{2}}=\frac{q^{3}M_{1}\sin^{3}{i}}{(1+q)^{2}},$
(1)

where $q$ is the ratio of the secondary to primary mass $M_{2}/M_{1}$, and we
have assumed a circular orbit ($e=0$). However, the inclination is still
unknown, so we can only estimate the binary mass ratio as a function of
inclination, even with an estimate of the primary mass.

Where the primary masses are known, there are multiple statistical approaches
to deriving the distribution of binary mass ratios despite the problem of
unknown inclinations. One approach is using an averaged value of $\sin^{3}i$
(e.g., Aitken 1964), but this is valid only for very large samples of
binaries. Other approaches involve assuming a mass ratio distribution and
comparing the implied distribution in an observable like the binary mass
function to observations, either assuming randomly-oriented orbits (Jaschek &
Ferrer, 1972; Halbwachs, 1987; Trimble, 1990), or iteratively adjusting the
assumptions about the inclination and mass ratio distributions (Mazeh &
Goldberg, 1992). The underlying distribution of inclinations is also of
interest, both for binary orbits and single rotating stars. For example,
Jackson & Jeffries (2010) model the $\sin i$ distributions for G- and K-type
stars in the Pleiades and Alpha Persei clusters and find that rotational
modulation from starspots can be detected down to inclinations of
$30^{\circ}$.

For systems with main sequence primaries, color (or luminosity) is a good
proxy for primary mass, allowing use of these types of methods for SB1 orbits.
Particularly within star clusters, it is also possible to estimate the mass of
a main sequence primary with a stellar evolutionary track or isochrone. In
these cases, the binary mass ratio is of central interest. For example, Geller
et al. (2021) convert distributions in primary mass and binary mass function
to distributions in secondary mass and mass ratio for main sequence binaries
in M67 using the algorithm from Mazeh & Goldberg (1992). They find a uniform
distribution of mass ratios ranging from $\sim 0.1$–$1$, rather than IMF-
distributed companion masses. Raghavan et al. (2010) also find a roughly
uniform distribution of mass ratios for main sequence field binaries from
Hipparchos, but with an overabundance of “twin" systems (with mass ratios
$q\geq 0.95$). Pinsonneault & Stanek (2006) show that 50% of a sample of OB-
type detached eclipsing binaries in the Large Magellanic Cloud are twin
systems with $q>0.87$. For close binaries, uniformly distributed mass ratios,
or correlated companion masses, are expected and usually observed (e.g.,
Tokovinin 2000; Sana et al. 2011).

Estimating stellar mass based on luminosity is difficult for evolved
primaries, since there is no mass-luminosity relation for giants. It is
therefore more difficult to measure the binary properties of SB1s with red
giant primaries. Leiner et al. (2022) used the distribution of their sample of
heavily spotted red giant binaries (RS CVn; Hall 1976) in luminosity versus
rotation period to provide a rough upper limit on their primary masses of
$1$–$2M_{\odot}$. They assumed edge-on, circular, and tidally locked orbits,
and matched the maximum radius the primary could have before exceeding its
Roche lobe to the corresponding stellar mass on an isochrone. They also
examined a sample of sub-subgiants, an under-luminous and shorter-period
subclass of RS CVn, which tend to have lower masses and tidally synchronize
only at shorter periods ($\lesssim 10$ days) than the typical period range for
RS CVn ($\sim$tens of days).

Here we take a slightly different statistical approach to previous work,
presenting a method to estimate the average primary mass, minimum mass ratio,
minimum inclination limit for recognizing rotational variables, and primary
mass spread for binaries at the population level. The method does not require
full SB2 statistics, eclipses, or a prior estimate of the primary mass. In
doing so, we aim to understand the five groups of binary spotted rotational
variable stars in the All-Sky Automated Survey for Supernovae (ASAS-SN)
variable star catalogue (Jayasinghe et al., 2018, 2019a, 2019b, 2020, 2021;
Christy et al., 2023) examined in Phillips et al. (2024). We fit circular
orbits at the ASAS-SN rotation period to sparse radial velocity (RV) curves
from the Apache Point Observatory Galactic Evolution Experiment (APOGEE;
Majewski et al. 2017) to determine their binary mass functions. We then fit
each group’s binary mass function distribution using Monte Carlo simulations
to estimate the parameters. In Phillips et al. (2024), we found that each
group consists of tidally synchronized or pseudosynchronized systems. The MS2
stars are binary main sequence stars with periods of $\lesssim 10$ days. The
G1 and G3 stars are RS CVn binaries with periods of tens of days. The G2 stars
are RS CVn and sub-subgiant binaries, with lower luminosity primaries and
shorter periods ($\sim$10 days) than G1/3. The G4 stars have intermediate
luminosities to G1-3 but with rotation periods approaching 100 days that tend
to be modestly longer than their binary orbital periods. We ignore the MS1
group since it consists of single stars.

In Section 2 we describe the data used, our method for fitting circular orbits
to determine each system’s mass function, and the Monte Carlo simulations we
use to fit the mass function distributions. We present the resulting fits and
discuss implications for each population of binaries in Section 3, concluding
and discussing future directions in Section 4.

##  2 Methods

We use RVs from the allVisit catalog of APOGEE DR17 (Abdurro’uf et al., 2022).
We require a signal to noise ratio of $>$20, that N_COMPONENTS$\neq$0, and
that there are $\geq 4$ RV measurements. With $\geq 4$ RVs we can fit a
circular orbit with at least one degree of freedom as a check of the fit. We
assume the binaries have circular orbits with the photometric period of the
rotational variable (see Lurie et al. 2017 for a discussion of tidal
synchronization in orbits with periods $\lesssim$ tens of days). This appeared
to be true for all but the G4 group, which tended to be slightly
subsynchronous. The RV curves are fit as

$v=v_{0}+a\cos(2\pi(t-t_{0})/P)+b\sin(2\pi(t-t_{0})/P)$ (2)

for $v_{0}$, $a$, and $b$. The velocity semiamplitude $K=\sqrt{a^{2}+b^{2}}$
and the orbital period $P$ provide the system’s binary mass function (Equation
1).

Because the time span used to determine periods from the ASAS-SN light curve
($\Delta T_{A}\simeq 5$ years) is shorter than the time span of some of the
APOGEE data ($\Delta T_{s}\simeq 9$ years) and because the precision of
rotational periods is limited by star spot evolution, we allow small shifts in
the period when fitting the APOGEE velocities. A period shift of $\Delta P$
leads to a phase shift of $\Delta n=\Delta T\Delta P/P^{2}$ over the period
$\Delta T$, so the period accuracy needed for APOGEE is $\Delta P\simeq
0.1P^{2}/\Delta T$ (allowing a 10% phase shift over the span of the data). We
fit circular orbits over the period range $P_{\rm{rot}}\pm 30P^{2}/\Delta
T_{A}$ in steps of $\frac{\Delta P}{10}\frac{\Delta T_{S}}{\Delta T_{A}}$.
From this best fit, we perform a refined local search over $P_{\rm{new}}\pm
P^{2}_{\rm{new}}/\Delta T_{A}$ in steps of $\Delta P_{\rm{new}}/100$. Phillips
et al. (2024) found that the reported rotation period was sometimes half the
true period for the G1 group of spotted giants, so for these stars, we also
carried out fits assuming $P_{\rm{orb}}=2P_{\rm{rot}}$.

Some systems have RV measurements that do not well constrain both extrema of
their RV curve, allowing unphysical velocity semiamplitudes and mass
functions. To filter out these cases, we find the closest RV to the minimum
and maximum velocity phases of the RV curve and the distance in phase space
between these RVs and their respective peak. The larger of these two distances
for each system is shown against the estimated mass function in Figure 1. A
vertical line at $f=0.25\ M_{\odot}$ is the mass function for a system of two
$1\ M_{\odot}$ stars in an edge-on circular orbit. The systems with high mass
functions all have large phase distances so we limit our sample for further
analysis to exclude systems with phase distances $>0.25$.

Figure 1: The maximum phase distance between an RV measurement and a velocity
extrema. The vertical line at $f=0.25\ M_{\odot}$ is the mass function of an
edge-on circular orbit of two $1M_{\odot}$ stars. We discard systems above the
horizontal line at a phase distance of $0.25$.

Figure 2 shows the mass function distribution of the five groups of rotational
variables. We empirically use $f=0.001M_{\odot}$, denoted by a vertical line
in the figure, to separate single and binary stars. Our models will correct
for the tail of binaries with $f<0.001M_{\odot}$. Of the 2133 APOGEE
rotational variables examined in Phillips et al. (2024), we are left with 225
systems after the criteria for the number and signal to noise ratios of the
RVs, the quality of the orbital solution, and the minimum mass function are
applied. There are 76, 18, 82, 26, and 23 for the MS2, G1, G2, G3, and G4
groups, respectively.

Figure 2: Distribution in binary mass function of the five groups of
rotational variables from Phillips et al. (2024). A vertical line indicates
$f=0.001M_{\odot}$, the empirical division we draw between single and binary
stars.

We fit the distributions in mass function using four parameters: the (average)
primary mass $M_{1}$, the minimum mass ratio $q_{\rm{min}}$ for a uniform
distribution of secondary masses $q_{\rm{min}}M_{1}\leq M_{2}\leq M_{1}$, a
minimum inclination limit $i_{\rm{min}}$, and a primary mass spread $r$. We
vary the inclination limit to model the selection effect that we do not expect
to find pole-on rotators. The parameter $r$ allows a uniform distribution of
primary masses from $M_{\rm{min}}$ to $M_{\rm{max}}$, where
$M_{\rm{min}}=2M_{1}/(1+r)$ and $M_{\rm{max}}=2rM_{1}/(1+r)$ so that the
$M_{1}$ parameter is the average primary mass and
$M_{\rm{max}}=rM_{\rm{min}}$.

Figure 3: The integral distributions of $f$ individually varying each model
parameter: the primary mass $M_{1}$ (leftmost), the minimum mass ratio
$q_{\rm{min}}$ (center left), the inclination limit $i_{\rm{min}}$ (center),
and the primary mass spread $r$ (center right). The rightmost panel shows the
observed distributions in mass function for each group of rotational
variables. (Leftmost panel) Distributions in $f$ assuming $r=1.0$,
$q_{\rm{min}}=1.0$, and no inclination limit for 10 values of $M_{1}$ from
0.1–1.0$M_{\odot}$. (Center left panel) Distributions in $f$ assuming
$M_{1}=1.0M_{\odot}$, $r=1.0$, and no inclination limit for 10 equally spaced
values of $q_{\rm{min}}$ from 0.1–1.0. (Center panel) Distributions in $f$
assuming $M_{1}=1.0M_{\odot}$, $r=1$, $q_{\rm{min}}=1.0$, and inclination
limits $i_{\rm{min}}=0^{\circ},10^{\circ},20^{\circ},...,80^{\circ}$. (Center
right panel) Distributions in $f$ assuming $M_{1}=1.0M_{\odot}$,
$q_{\rm{min}}=1.0$ and no inclination limit for 10 equally spaced values of
$r$ from 1.0–3.0.

Figure 3 shows integral probability distributions in mass function
systematically varying each of the parameters $M_{1}$, $q_{\rm{min}}$,
$i_{\rm{min}}$, and $r$ while holding the others constant, along with the mass
function distributions for each group of rotational variables. The average
primary mass has the greatest effect on the mass function distribution, so we
expect to be most sensitive to $M_{1}$. The effects of varying $q_{\rm{min}}$
and $i_{\rm{min}}$ are partially degenerate, so we expect the estimates of
$q_{\rm{min}}$ to be sensitive to the model’s treatment of $i_{\rm{min}}$,
with lower estimates of $q_{\rm{min}}$ when we have $i_{\rm{min}}\neq 0$.
Despite it’s noticeable effect on the distribution, we are limited in
sensitivity to low $i_{\rm{min}}$ because of our binary criteria. For example,
for $M_{1}=1M_{\odot}$ and $q=1,0.8,0.6,0.4$, and $0.2$, requiring that
$f>0.001M_{\odot}$ excludes systems with $i<9,11,13,18,$ and $34$ degrees even
before including any inclination limit $i_{\rm{min}}$. The effect of varying
$r$ is relatively small compared to the other parameters, so we expect limited
sensitivity to the primary mass spread. Because it mainly affects the high
mass function end of the distribution, we expect $r$ to also be sensitive to
mass functions that have been overestimated.

It is possible that some of the mass function estimates are simply wrong
either because of velocity sampling or a failure in one of our assumptions
(tidally locked circular orbits). This is a particular problem at the high
mass function end of the distribution where we must have $M_{1}>4f$ in order
to find a non-zero probability. We accommodate this by assuming that there is
some probability $p$ for a mass function value to be good, and a probability
$1-p$ for it to be bad. We can define a vector with $v_{i}=0$ if $f_{i}$ is
bad and $v_{i}=1$ if $f_{i}$ is good. The Bayesian probability of a model with
$M_{1}$ and $q_{\rm{min}}$ is

$P(M_{1},q|D)\propto\int
dpP(p)\sum_{v}\Pi_{v_{i}=1}pP(f_{i}|M_{1},q)\Pi_{v_{i}=0}(1-p)P_{B},$ (3)

where $P(p)$ is a prior on the probability $p$ and $P_{B}$ is the probability
for “bad" $f_{i}$ values (Press 1997). Equation 3 can be reduced to (Kochanek,
private communication to Press)

$P(M_{1},q|D)\propto P(M_{1},q)\int
dpP(p)\Pi_{i}(pP(f_{i}|M_{1},q)+(1-p)P_{B})$ (4)

for priors $P(M_{1},q_{\rm{min}})$ on $M_{1}$ and $q_{\rm{min}}$. We take
$P_{B}=1/(M_{\rm{max}}-M_{\rm{min}})$ corresponding to the range of average
primary masses $M_{1}$ over which we perform Monte Carlo simulations (see
below). We calculate the probability distribution of the “good" fraction $p$
as

$P(p|D)\propto P(p)\int
dM_{1}dqP(M_{1},q)\Pi_{i}(pP(f_{i}|M_{1},q)+(1-p)P_{B}).$ (5)

For $r=1$, we ran models with $0.1M_{\odot}\leq M_{1}\leq 2M_{\odot}$,
$0.01\leq q\leq 1$, $0\leq p<1$ and $0^{\circ}\leq i_{\rm{min}}\leq
80^{\circ}$ with 100, 100, 100, and 9 (steps of 10 degrees) points. We also
ran models with $i_{\rm{min}}=0^{\circ}$ and $i_{\rm{min}}=30^{\circ}$ using
10 values of $r$ ranging from $r=1$ to $r=3$. We did not run a full five
dimensional grid of models including $r$ because, as expected, the data are
not sensitive to the value of $r$. We assumed uniform priors for all of the
variables.

To compute the mass function probability distribution for a model, we drew
100000 samples using uniform deviates for the allowed ranges of $M_{1}$,
$M_{2}$, and $\cos i$, and then disposed of samples where $f<0.001M_{\odot}$
to match our criteria for binarity. We sorted the remaining $N$ samples to
find the integral probability distribution $P(<f_{i})=i/N$. To estimate the
differential probability of a mass function $f$ in our data, we find the index
$i$ in the integral distribution where $f_{i}\leq f<f_{i+1}$. Using $m=500$,
we find the indices $f_{i-m}$ and $f_{i+m}$ and approximate the differential
probability as

$P(f|M_{1},q_{\rm{min}})=\frac{dP}{df}\simeq\frac{2m}{N}\frac{1}{f_{i+m}-f_{i-m}}$
(6)

(with corrections at the distribution’s edges). If $f$ lies outside the range
of the simulation data, we assign the differential probability
$P(f_{i}|M_{1},q_{\rm{min}})=0$. To get the estimated values of any parameter,
we marginalize over the other parameters and report the median and
uncertainties as the values symmetrically enclosing 90% of the probability.
The probability distributions for $M_{1}$, $q_{\rm{min}}$ and $p$ are
generally well defined so we do not show many examples of the distributions
and only report the results. The distributions for $i_{\rm{min}}$ and $r$ are
more problematic, and we show them to illustrate the issues and aid in the
discussion.

To demonstrate and test the method, we generated three synthetic data sets
with 33 synthetic mass functions, the median number of systems before imposing
the $f>0.001M_{\odot}$ criteria for binarity, and with no inclination limit
($i_{\rm{min}}=0^{\circ}$), $i_{\rm{min}}=20^{\circ}$ or
$i_{\rm{min}}=60^{\circ}$. We generated the synthetic data using
$M_{1}=1M_{\odot}$, $r=2.0$, $q_{\rm{min}}=0.1$. For the data we generated,
none of the mass functions were below our limit of $0.001M_{\odot}$. We
randomly selected two of the mass functions from each data set and reset them
to $1.0$ and $1.5M_{\odot}$, meaning that the “good" fraction for each group
is $p=0.94$. We report the resulting parameter estimates in Table 1 and show
the likelihood distributions for $M_{1}$, $q_{\rm{min}}$, and $p$ in Figure 4.
The estimates of $M_{1}$ and $p$ are consistent with the ground truth in all
cases. As expected, there is a degeneracy between $q_{\rm{min}}$ and the
inclination limit. This is particularly clear for the data set with
$i_{\rm{min}}=60^{\circ}$, where, in cases where we fix
$i_{\rm{min}}=0^{\circ}$, the estimate of $q_{\rm{min}}$ is inconsistent with
the ground truth.

Table 1: The medians and 90% confidence ranges for the average primary mass $M_{1}$ (in solar units), minimum mass ratio $q_{\rm{min}}$, and “good" fraction $p$ for three synthetic data sets generated with $M_{1}=1M_{\odot}$, $q_{\rm{min}}=0.1$, $r=2$, $p\simeq 0.94$, and either no inclination limit, $i_{\rm{min}}=20^{\circ}$, or $60^{\circ}$. We report results from each Monte Carlo fit: (1) assuming no minimum inclination and no primary mass spread, (2) marginalizing over inclination limits $i_{\rm{min}}$ but with no primary mass spread, (3) marginalizing over primary mass spread $r$ with no inclination limit ($i_{\rm{min}}=0^{\circ}$), or (4) with $i_{\rm{min}}=30^{\circ}$. Model | $M_{1}\ (M_{\odot})$ | $q_{\rm{min}}$ | $p$
---|---|---|---
Synthetic data, $i_{\rm{min}}=0^{\circ}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.98^{+0.33}_{-0.36}$ | $0.11^{+0.11}_{-0.09}$ | $0.90^{+0.07}_{-0.11}$
Varying $i_{\rm{min}}$, $r=1$ | $0.94^{+0.33}_{-0.38}$ | $0.10^{+0.11}_{-0.09}$ | $0.90^{+0.07}_{-0.12}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.89^{+0.38}_{-0.27}$ | $0.12^{+0.12}_{-0.10}$ | $0.91^{+0.06}_{-0.12}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.85^{+0.36}_{-0.25}$ | $0.09^{+0.11}_{-0.08}$ | $0.90^{+0.07}_{-0.11}$
Synthetic data, $i_{\rm{min}}=20^{\circ}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.98^{+0.31}_{-0.35}$ | $0.13^{+0.12}_{-0.11}$ | $0.90^{+0.07}_{-0.11}$
Varying $i_{\rm{min}}$, $r=1$ | $0.94^{+0.33}_{-0.38}$ | $0.11^{+0.12}_{-0.09}$ | $0.90^{+0.07}_{-0.12}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.89^{+0.38}_{-0.27}$ | $0.13^{+0.14}_{-0.11}$ | $0.91^{+0.06}_{-0.12}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.85^{+0.35}_{-0.25}$ | $0.10^{+0.11}_{-0.09}$ | $0.90^{+0.07}_{-0.11}$
Synthetic data, $i_{\rm{min}}=60^{\circ}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $1.02^{+0.31}_{-0.12}$ | $0.49^{+0.14}_{-0.24}$ | $0.91^{+0.06}_{-0.11}$
Varying $i_{\rm{min}}$, $r=1$ | $1.02^{+0.31}_{-0.12}$ | $0.41^{+0.20}_{-0.39}$ | $0.91^{+0.06}_{-0.11}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.94^{+0.36}_{-0.23}$ | $0.54^{+0.29}_{-0.24}$ | $0.91^{+0.06}_{-0.11}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $1.16^{+0.46}_{-0.25}$ | $0.09^{+0.12}_{-0.08}$ | $0.91^{+0.06}_{-0.12}$
Figure 4: The likelihood of each $M_{1}$ (left), $q_{\rm{min}}$ (center), and
$p$ (right) for the synthetic data sets allowing all inclinations (top), with
$i_{\rm{min}}=20^{\circ}$ (center), and with $i_{\rm{min}}=60^{\circ}$
(bottom). The corresponding medians and 90% confidence ranges of the integral
probability for each distribution are reported in Table 1. The ground truth
values of each parameter are marked with vertical red lines. Some of the
distributions very closely overlap.

The likelihoods for each inclination limit and primary mass spread for each
synthetic data set are shown in Figure 5.

Figure 5: The likelihood of $i_{\rm{min}}$ (left) and $r$ (right) after
marginalizing over $M_{1}$, $q_{\rm{min}}$, and $p$ for synthetic data with
$M_{1}=1\ M_{\odot},\ q_{\rm{min}}=0.1,\ r=2.0,$ and no inclination limit
(top), with $i>20^{\circ}$ (center), or with $i>60^{\circ}$ (bottom). In the
right panels, the solid (dashed) lines are for $i_{\rm{min}}=0^{\circ}$
($30^{\circ}$).

All of the data sets have maximum likelihoods at $i_{\rm{min}}=0^{\circ}$,
apparently favoring finding all systems with $f>0.001M_{\odot}$. The
distributions for the $i_{\rm{min}}=0^{\circ}$ and $i_{\rm{min}}=20^{\circ}$
data sets are nearly identical. The $i_{\rm{min}}=60^{\circ}$ data shows a
local maximum at $i_{\rm{min}}=50^{\circ}$ rather that $60^{\circ}$, but does
not recover the ground truth as the global likelihood maximum. All of the
synthetic data sets have roughly flat likelihood distributions in $r$, but the
likelihood differences are too small for parameter estimation.

##  3 Results and Discussion

Table 2 shows the results for $M_{1}$, $q_{\rm{min}}$ and $p$ from our four
standard fits for each group of rotational variables. The estimates for all
three parameters are consistent between models, and each group favors little
to no contamination by bad data. The uncertainties strongly depend on the
number of systems available for the analysis. They are very good for MS2 and
G2 with 76 and 82 systems, and significantly worse for G1, G3 and G4, with
only 18, 26, and 23 systems.

As expected from their locations on a color magnitude diagram, the MS2
primaries are lower main sequence stars with median masses $M_{1}\simeq
0.75M_{\odot}$. That the giant groups tend to have primary masses $\gtrsim
1M_{\odot}$ is required since less massive primaries have not yet evolved off
of the main sequence (e.g., based on a 10 Gyr PARSEC isochrone; Bressan et al.
2012; Marigo et al. 2013). The G2 subgiants favor $M_{1}\simeq 0.94M_{\odot}$.
In Phillips et al. (2024), we concluded that the G1 and G3 giants were likely
a single group of RS CVn type binaries. While our estimates of $M_{1}$ and
$q_{\rm{min}}$ for these two groups are consistent, the uncertainties are so
large that this is not a strong requirement. Analyzed separately, the G1 group
tends toward smaller primary mass estimates $M_{1}\simeq 0.8M_{\odot}$
(somewhat low for a red giant, but the uncertainties extend above
$1M_{\odot}$) compared to G3 at $M_{1}\simeq 1.1M_{\odot}$. Analyzing the G1
and G3 groups jointly gives $M_{1}\simeq 1.2M_{\odot}$, more massive than the
estimates for either of the groups on their own and the most massive of the
giant groups. The G4 group has $M_{1}\simeq 1.1M_{\odot}$.

Examining the companion mass ranges, we clearly see the expected partial
degeneracy with $i_{\rm{min}}$, as the models that vary $i_{\rm{min}}$ or fix
$i_{\rm{min}}=30^{\circ}$ all result in lower $q_{\rm{min}}$ estimates than
the models with $i_{\rm{min}}=0^{\circ}$. The MS2 stars all have a restricted
range of companion masses with $q_{\rm{min}}\simeq 0.4$ even though it would
take $q_{\rm{min}}\simeq 0.1$ to reach the Hydrogen burning limit near
$0.1M_{\odot}$ given the estimates of $M_{1}$. The giant groups all allow a
broader range of companion masses with $q_{\rm{min}}\sim 0.1$–$0.2$. In
particular, the G4 group seems to have the weakest dependence of
$q_{\rm{min}}$ on the treatment of inclination, ranging only from
$q_{\rm{min}}=0.10$–$0.12$.

Table 2: The medians and 90% confidence ranges for the average primary mass $M_{1}$ (in solar units), minimum mass ratio $q_{\rm{min}}$, and “good" fraction $p$ for the MS2-G4 groups resulting from the same four fits as in Table 1. Model | $M_{1}\ (M_{\odot})$ | $q_{\rm{min}}$ | $p$
---|---|---|---
MS2 (76 systems) |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.75^{+0.12}_{-0.13}$ | $0.43^{+0.14}_{-0.16}$ | $0.97^{+0.02}_{-0.04}$
Varying $i_{\rm{min}}$, $r=1$ | $0.75^{+0.12}_{-0.13}$ | $0.36^{+0.19}_{-0.32}$ | $0.97^{+0.02}_{-0.04}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.73^{+0.19}_{-0.19}$ | $0.47^{+0.24}_{-0.19}$ | $0.98^{+0.01}_{-0.04}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.83^{+0.21}_{-0.15}$ | $0.19^{+0.07}_{-0.16}$ | $0.98^{+0.01}_{-0.04}$
G1 (18 systems) |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.77^{+0.96}_{-0.21}$ | $0.22^{+0.28}_{-0.19}$ | $0.92^{+0.07}_{-0.16}$
Varying $i_{\rm{min}}$, $r=1$ | $0.71^{+0.94}_{-0.17}$ | $0.15^{+0.27}_{-0.13}$ | $0.92^{+0.07}_{-0.16}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.93^{+0.79}_{-0.40}$ | $0.21^{+0.32}_{-0.18}$ | $0.93^{+0.06}_{-0.16}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.91^{+0.73}_{-0.36}$ | $0.13^{+0.16}_{-0.11}$ | $0.93^{+0.06}_{-0.16}$
G2 (82 systems) |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.94^{+0.17}_{-0.10}$ | $0.21^{+0.08}_{-0.16}$ | $0.99^{+0.00}_{-0.04}$
Varying $i_{\rm{min}}$, $r=1$ | $0.89^{+0.12}_{-0.04}$ | $0.08^{+0.13}_{-0.06}$ | $0.99^{+0.00}_{-0.04}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.96^{+0.19}_{-0.12}$ | $0.21^{+0.09}_{-0.16}$ | $0.99^{+0.00}_{-0.04}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.96^{+0.19}_{-0.12}$ | $0.09^{+0.11}_{-0.07}$ | $0.99^{+0.00}_{-0.04}$
G3 (26 systems) |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $1.14^{+0.44}_{-0.21}$ | $0.26^{+0.17}_{-0.21}$ | $0.93^{+0.05}_{-0.13}$
Varying $i_{\rm{min}}$, $r=1$ | $1.12^{+0.42}_{-0.21}$ | $0.19^{+0.21}_{-0.16}$ | $0.93^{+0.05}_{-0.13}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $1.08^{+0.50}_{-0.31}$ | $0.28^{+0.19}_{-0.23}$ | $0.93^{+0.05}_{-0.13}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $1.10^{+0.46}_{-0.29}$ | $0.14^{+0.14}_{-0.12}$ | $0.93^{+0.05}_{-0.13}$
G4 (23 systems) |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $1.12^{+0.50}_{-0.19}$ | $0.12^{+0.16}_{-0.10}$ | $0.96^{+0.03}_{-0.10}$
Varying $i_{\rm{min}}$, $r=1$ | $1.10^{+0.48}_{-0.17}$ | $0.11^{+0.13}_{-0.09}$ | $0.96^{+0.03}_{-0.10}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $1.10^{+0.54}_{-0.29}$ | $0.12^{+0.18}_{-0.10}$ | $0.96^{+0.03}_{-0.11}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $1.06^{+0.48}_{-0.29}$ | $0.10^{+0.09}_{-0.09}$ | $0.96^{+0.03}_{-0.11}$
G1+G3 (44 systems) |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $1.29^{+0.31}_{-0.42}$ | $0.20^{+0.15}_{-0.17}$ | $0.94^{+0.05}_{-0.09}$
Varying $i_{\rm{min}}$, $r=1$ | $1.21^{+0.36}_{-0.36}$ | $0.16^{+0.16}_{-0.14}$ | $0.94^{+0.05}_{-0.09}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $1.08^{+0.42}_{-0.29}$ | $0.24^{+0.16}_{-0.19}$ | $0.95^{+0.04}_{-0.09}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $1.06^{+0.38}_{-0.25}$ | $0.13^{+0.10}_{-0.11}$ | $0.95^{+0.04}_{-0.09}$

Since more luminous main sequence primaries should have higher masses, we
split the MS2 group of rotators into three bins in Gaia absolute G magnitude
with equal numbers of objects per magnitude bin. The estimates for $M_{1}$,
$q_{\rm{min}}$, and $p$ for each magnitude bin are reported in Table 3. Each
bin favors little contamination by bad data. We also report mass ranges for a
10 Gyr PARSEC isochrone (Bressan et al., 2012; Marigo et al., 2013) for each
magnitude bin after doubling the luminosity at fixed color to align it with
the “binary" main sequence. As expected, the primary masses increase with
increasing luminosity, with estimates of $\sim 0.62M_{\odot}$ in the dimmest
bin and $\sim 0.90M_{\odot}$ in the brightest bin. The 90% confidence regions
for primary mass estimates are consistent with the isochrone mass ranges in
each bin.

There is a selection effect pointed out by Simonian et al. (2019) that can
affect estimates of $q_{\rm{min}}$ for the MS2 sample. If luminosity on the
main sequence scales as $L\propto M^{3}$, then the luminosity of a binary is
$L\propto(1+q^{3})$, so the distance at which a binary can be found in a flux-
limited sample scales as $(1+q^{3})^{1/2}$. For a constant density, spherical
distribution, the number of binaries would then scale as $(1+q^{3})^{3/2}$ —
the actual scaling will be weaker because most of the stars are in a disk,
which asymptotically reduces the scaling to $(1+q^{2})$, and dust. Such a
selection effect reduces the numbers of lower $q$ systems in the sample, and
thus would bias the models to higher values of $q_{\rm{min}}$. We do not
explore this in detail here because doing it correctly is a significant
project that seems unwarranted given our small sample size, which make the
statistical uncertainties on $q_{\rm{min}}$ significant, and because the
parallax error requirement on the sample (see Phillips et al. 2024) that
$\varpi/\sigma_{\varpi}>10$ also makes much of our sample volume limited
rather than flux limited. In a volume limited sample there is no such bias on
$q_{\rm{min}}$. This question will need to be addressed for larger samples.

Table 3: The medians and 90% confidence ranges for the average primary mass $M_{1}$, minimum mass ratio $q_{\rm{min}}$, and “good" fraction $p$ for the MS2 group split into three bins in Gaia absolute $G$ magnitude for the same four fits as in Tables 1 and 2. We also report the mass ranges for a 10 Gyr binary PARSEC isochrone in each magnitude bin, $M_{\rm{iso}}$. Model | $M_{1}\ (M_{\odot})$ | $q_{\rm{min}}$ | $p$
---|---|---|---
$4.78\leq M_{G}<5.91$ |  |  |
$M_{\rm{iso}}=0.71$–$0.82\ M_{\odot}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.81^{+0.63}_{-0.19}$ | $0.43^{+0.36}_{-0.38}$ | $0.95^{+0.04}_{-0.12}$
Varying $i_{\rm{min}}$, $r=1$ | $0.81^{+0.61}_{-0.19}$ | $0.13^{+0.53}_{-0.11}$ | $0.95^{+0.04}_{-0.12}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.91^{+0.50}_{-0.38}$ | $0.44^{+0.42}_{-0.37}$ | $0.96^{+0.03}_{-0.11}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $1.08^{+0.46}_{-0.33}$ | $0.09^{+0.12}_{-0.08}$ | $0.96^{+0.03}_{-0.10}$
$5.92\leq M_{G}<6.54$ |  |  |
$M_{\rm{iso}}=0.65$–$0.71\ M_{\odot}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.81^{+0.31}_{-0.21}$ | $0.40^{+0.20}_{-0.20}$ | $0.96^{+0.03}_{-0.09}$
Varying $i_{\rm{min}}$, $r=1$ | $0.81^{+0.29}_{-0.21}$ | $0.33^{+0.23}_{-0.28}$ | $0.96^{+0.03}_{-0.10}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.73^{+0.35}_{-0.21}$ | $0.46^{+0.30}_{-0.22}$ | $0.96^{+0.03}_{-0.09}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.81^{+0.36}_{-0.21}$ | $0.26^{+0.15}_{-0.21}$ | $0.96^{+0.03}_{-0.09}$
$6.55\leq M_{G}<8.57$ |  |  |
$M_{\rm{iso}}=0.47$–$0.65\ M_{\odot}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.62^{+0.25}_{-0.10}$ | $0.35^{+0.30}_{-0.23}$ | $0.97^{+0.02}_{-0.09}$
Varying $i_{\rm{min}}$, $r=1$ | $0.60^{+0.25}_{-0.08}$ | $0.23^{+0.32}_{-0.20}$ | $0.97^{+0.02}_{-0.09}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.60^{+0.31}_{-0.17}$ | $0.39^{+0.45}_{-0.25}$ | $0.97^{+0.02}_{-0.09}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.64^{+0.31}_{-0.13}$ | $0.24^{+0.14}_{-0.20}$ | $0.97^{+0.02}_{-0.09}$

Figure 6 shows the likelihood distributions for the inclination limit and
primary mass spread (both assuming $i_{\rm{min}}=0^{\circ}$ and
$i_{\rm{min}}=30^{\circ}$). We discuss the structure of likelihood differences
across $i_{\rm{min}}$ and $r$, but note that the differences are not
particularly significant. With the exception of the G2 group, there is no
strong indication of an inclination bias (up to the insensitivity we discussed
in Section 2). G1, G3 (both separated and combined), and G4 favor finding most
if not all systems with $f>0.001M_{\odot}$ regardless of inclination, with
maximum likelihoods $i_{\rm{min}}=10^{\circ}$. MS2 has a maximum likelihood at
$i_{\rm{min}}=0^{\circ}$, but its likelihood does not monotonically decrease
with increasing $i_{\rm{min}}$. Instead, it mirrors the behavior of the
synthetic data in the bottom left panel of Figure 5, with a local maximum at
$i_{\rm{min}}=50^{\circ}$, suggesting a potential inclination limit $\gtrsim
30^{\circ}$. The G2 group has a maximum likelihood at
$i_{\rm{min}}=60^{\circ}$, preferring to be viewed quasi-equatorially.

For the values of $r$, all but the G2 group have roughly flat distributions,
with little distinction between assuming $i_{\rm{min}}=0^{\circ}$ and
$i_{\rm{min}}=30^{\circ}$. The MS2 group has maximum likelihoods at $r=3.0$
and $r=2.11$ when we assume $i_{\rm{min}}=0^{\circ}$ and
$i_{\rm{min}}=30^{\circ}$, respectively. The G1, G3 (both separated and
combined), and G4 groups have maximum likelihoods at $r=3$ regardless of
whether we assume $i_{\rm{min}}=0^{\circ}$ or $30^{\circ}$. The G2 group
favors no spread in primary mass, with a maximum likelihood at $r=1$.

Figure 6: The likelihood of $i_{\rm{min}}$ (left) and $r$ (right) after
marginalizing over $p$, $M_{1}$, and $q_{\rm{min}}$ for each group of
rotational variables. In the right panels, the solid (dashed) lines are for
$i_{\rm{min}}=0^{\circ}$ ($30^{\circ}$).

##  4 Conclusions

We presented a method to estimate the average primary mass, secondary mass
range, inclination limit, primary mass spread, and fraction of correctly
estimated binary mass functions for populations binary stars, and applied it
to five previously identified populations of binary rotational variables.
Given that we recover consistent values of the good data probability $p$ with
the ground truth in all cases of synthetic data, we are reasonably confident
that we correctly estimate it. Each group of rotational variables favors
little to no contamination by incorrectly estimated mass functions.

We are also reasonably confident in our estimates of average primary mass,
which we also recover correctly in the synthetic data tests. In the rotational
variables, we find MS2 stars to have average primary masses of $\sim
0.75M_{\odot}$, as expected from their location on the CMD. We estimate
different average primary masses for G1 ($\sim 0.8M_{\odot}$) and G3 ($\sim
1.10M_{\odot}$) when analyzing them separately, and higher masses ($\sim
1.2M_{\odot}$) when analyzing them jointly. The uncertainties on the separated
G1 and G3 mass estimates are large enough that we are still unable to
definitively say whether they are physically different groups, or a single
group as we posited in Phillips et al. (2024). Regardless, all of these mass
estimates are consistent with the $1$–$2M_{\odot}$ upper limit on RS CVn
primaries from Leiner et al. (2022). The G2 stars are somewhat less massive
than the G1 and G3 groups analyzed jointly, at $\sim 0.94M_{\odot}$,
consistent with their characterization as sub-subgiants which should be less
massive than typical RS CVn. For the G4 group we estimate average primary
masses intermediate to those in G1-3, around $1.10M_{\odot}$. We also confirm
that our estimates of $M_{1}$ increase with increasing luminosity for the MS2
main sequence binaries.

For populations without a strong inclination bias, we are also relatively
confident in our estimates of minimum mass ratio. We correctly recover
$q_{\rm{min}}$ in most cases of synthetic data, failing to do so only when we
assume $i_{\rm{min}}=0^{\circ}$ for synthetic data constructed with a high
inclination limit. The estimates of $q_{\rm{min}}$ are partially degenerate
with the model’s treatment of inclination, and this effect becomes stronger at
higher true inclination limits in synthetic data. The MS2 stars seem to have a
restricted range of companion masses, with estimates of $q_{\rm{min}}\sim
0.4$, though we have reason to believe they have an inclination limit above
our insensitivity threshold of $\sim 20^{\circ}$ (see below) and therefore may
actually have companions down to lower mass ratios. The giant groups have
$q_{\rm{min}}\simeq 0.1$–$0.2$, with the weakest dependence of $q_{\rm{min}}$
on $i_{\rm{min}}$ for the G4 group.

We are not sensitive to $i_{\rm{min}}\lesssim 20$–$30^{\circ}$ because of our
minimum mass function criteria (excluding $f<0.001$). The G1 and G3 (analyzed
both separately and jointly) and G4 stars all show maximum likelihood
inclination limits near $0^{\circ}$, favoring finding all systems with
$f>0.001M_{\odot}$ and indicating no strong inclination bias. For higher
inclination limits, our tests with synthetic data show that we are only
somewhat sensitive to $i_{\rm{min}}$, observing local maxima in likelihood at
high inclinations, but with a maximum likelihood at $i_{\rm{min}}=0^{\circ}$.
Such a shape of the $i_{\rm{min}}$ likelihood distribution may be useful as a
diagnostic for whether it is necessary to consider parameter estimates from
models varying $i_{\rm{min}}$. This is the case for the MS2 group which has a
local maximum at $i_{\rm{min}}=40^{\circ}$, much like the synthetic data with
$i_{\rm{min}}=60^{\circ}$, indicating that it may have a true inclination
limit $\gtrsim 30$–$50^{\circ}$.

We are not sensitive to the primary mass spread and are unable to recover the
correct value for $r$ as the maximum likelihood in the synthetic data. The
likelihoods in $r$ are relatively flat both for the synthetic data and for
each group of rotational variables, providing no evidence for a preferred
primary mass spread.

The volume of data available for this type of analysis should increase rapidly
increase with the SDSS Milky Way Mapper (Kollmeier et al., 2017) extension of
the APOGEE survey and Gaia DR4. We explored the expected improvement by
repeating each of the synthetic mass function tests with 10 times as much data
(including 10 times as many “bad" mass functions). We report the results for
$M_{1}$, $q_{\rm{min}}$, and $p$ from each model in Table 4. The primary mass
estimates are consistent with the ground truth value of $1M_{\rm{\odot}}$, are
more uniform across models, and have significantly smaller uncertainties
compared to those in Table 1, especially for the data with
$i_{\rm{min}}=0^{\circ}$ and $20^{\circ}$. In the minimum mass fractions, we
still see the dependence of $q_{\rm{min}}$ on inclination. For the
$i_{\rm{min}}=0^{\circ}$ and $i_{\rm{min}}=20^{\circ}$ data sets, even the
values of $q_{\rm{min}}$ where the model assumes no inclination limit are
consistent with the ground truth value of $q_{\rm{min}}=0.1$. For the
$i_{\rm{min}}=60^{\circ}$ data set, the values of $q_{\rm{min}}$ when the
model fixes $i_{\rm{min}}=0^{\circ}$ ($q_{\rm{min}}\sim 0.46$) are
inconsistent with the ground truth, but given that we will be better able to
constrain the inclination limits with more data (see below), this is less of a
cause for concern.

Table 4: The estimates and 90% confidence ranges for the average primary mass $M_{1}$ (in solar units), minimum mass ratio $q_{\rm{min}}$, and “good" fraction $p$ for three 10 times larger synthetic data sets generated with $M_{1}=1M_{\odot},\ q_{\rm{min}}=0.1,\ r=2$, and either $i_{\rm{min}}=0^{\circ}$, $i_{\rm{min}}=20^{\circ}$, or $i_{\rm{min}}=60^{\circ}$ for the same 4 fits as in Tables 1–3. Model | $M_{1}$ | $q_{\rm{min}}$ | $p$
---|---|---|---
Synthetic data, $i_{\rm{min}}=0$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $0.96^{+0.08}_{-0.06}$ | $0.11^{+0.06}_{-0.09}$ | $0.92^{+0.03}_{-0.03}$
Varying $i_{\rm{min}}$, $r=1$ | $0.96^{+0.08}_{-0.06}$ | $0.10^{+0.06}_{-0.08}$ | $0.92^{+0.03}_{-0.03}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.96^{+0.08}_{-0.08}$ | $0.11^{+0.07}_{-0.09}$ | $0.92^{+0.03}_{-0.03}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.94^{+0.08}_{-0.10}$ | $0.08^{+0.05}_{-0.06}$ | $0.92^{+0.03}_{-0.03}$
Synthetic data, $i_{\rm{min}}=20^{\circ}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $1.02^{+0.06}_{-0.10}$ | $0.13^{+0.06}_{-0.11}$ | $0.92^{+0.03}_{-0.03}$
Varying $i_{\rm{min}}$, $r=1$ | $1.00^{+0.06}_{-0.10}$ | $0.11^{+0.07}_{-0.09}$ | $0.92^{+0.03}_{-0.03}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $0.98^{+0.08}_{-0.08}$ | $0.14^{+0.05}_{-0.11}$ | $0.92^{+0.03}_{-0.03}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $0.96^{+0.08}_{-0.10}$ | $0.08^{+0.06}_{-0.06}$ | $0.92^{+0.03}_{-0.03}$
Synthetic data, $i_{\rm{min}}=60^{\circ}$ |  |  |
$i_{\rm{min}}=0^{\circ},\ r=1$ | $1.10^{+0.06}_{-0.06}$ | $0.32^{+0.08}_{-0.06}$ | $0.92^{+0.03}_{-0.02}$
Varying $i_{\rm{min}}$, $r=1$ | $1.10^{+0.06}_{-0.08}$ | $0.07^{+0.05}_{-0.06}$ | $0.92^{+0.03}_{-0.02}$
Varying $r$, $i_{\rm{min}}=0^{\circ}$ | $1.08^{+0.06}_{-0.10}$ | $0.33^{+0.11}_{-0.06}$ | $0.92^{+0.03}_{-0.02}$
Varying $r$, $i_{\rm{min}}=30^{\circ}$ | $1.12^{+0.08}_{-0.06}$ | $0.08^{+0.04}_{-0.07}$ | $0.92^{+0.03}_{-0.02}$

Figure 7 shows the likelihoods of $i_{\rm{min}}$ and $r$ for the fits to the
longer synthetic data tests. For the inclination limits, the data with no
inclination limit and with $i_{\rm{min}}=20^{\circ}$ both have maximum
likelihoods at at $i_{\rm{min}}=20^{\circ}$, and decrease steeply in
likelihood at higher inclination limits. The $i_{\rm{min}}=60^{\circ}$ data
set has a maximum likelihood at $i_{\rm{min}}=50^{\circ}$, meaning that with
$\sim 10$ times as much data, we will be better able to constrain
$i_{\rm{min}}$.

The likelihoods for the primary mass spreads are still flat compared to the
likelihood distributions for $i_{\rm{min}}$. All fits favor no primary mass
spread, with maximum likelihoods at or near $r=1$, though the $i_{\rm{min}}=0$
and $i_{\rm{min}}=20$ data sets show hints at local maxima near $r=2$. We
checked whether the inclusion of “bad" mass functions was affecting our
estimates of $r$ by generating another synthetic data set with no bad points
and looking at the results with $p\equiv 1$. This yielded nearly identical
results. Thus, even with ten times the data it will be difficult to constrain
$r$.

Figure 7: The likelihood of $i_{\rm{min}}$ (left) and $r$ (right) after
marginalizing over $p$, $M_{1}$, and $q_{\rm{min}}$ for a larger synthetic
data set compared to the one described in Section 2 with $M_{1}=1\ M_{\odot},\
q_{\rm{min}}=0.1,\ r=2.0,$ and no inclination limit (top), with $i>20^{\circ}$
(center), or with $i>60^{\circ}$ (bottom). In the right panels, the solid
(dashed) lines are for $i_{\rm{min}}=0^{\circ}$ ($30^{\circ}$).

These tests show that with increasingly large samples, the constraining power
of this approach increases significantly, both in terms of better constraints
on the variables that were reasonably well-constrained for the present sample
($M_{1}$ and $q_{\rm{min}}$) but also will begin to constrain the variables
that were not well-constrained ($i_{\rm{min}}$ and $r$). SDSS Milky Way Mapper
(Kollmeier et al., 2017) is likely to provide an order of magnitude increase,
and Gaia DR4 is likely to provide an even larger expansion. The ASAS-SN sample
of rotational variables has a limited magnitude overlap with good Gaia
spectroscopic observations, so expanding the rotational variable sample being
used to brighter magnitudes (e.g., the samples from Oelkers et al. 2018, Chen
et al. 2020, Reinhold et al. 2023, or Claytor et al. 2024) can provide a still
larger increase. The method can also be applied to any population of SB1
binaries providing only mass functions.

## Acknowledgements

The authors thank Las Cumbres Observatory and its staff for their continued
support of ASAS-SN. AP thanks Dominick Rowan and Dr. Christine Mazzola Daher
for their guidance in using and interpreting APOGEE RV data, and Profs.
Krzysztof Stanek and Marc Pinsonneault for their helpful comments in preparing
the manuscript. CSK is supported by National Science Foundation grants
AST-2307385 and 2407206.

ASAS-SN is funded in part by the Gordon and Betty Moore Foundation through
grants GBMF5490 and GBMF10501, and by the Alfred P. Sloan Foundation through
grant G-2021-14192 to the Ohio State University.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions.

SDSS-IV acknowledges support and resources from the Center for High
Performance Computing at the University of Utah. The SDSS website is
www.sdss4.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

## Data Availability

All data used in this study are public. The catalogue of ASAS-SN variable
stars is available on the ASAS-SN variable stars database (https://asas-
sn.osu.edu/variables). APOGEE DR17 radial velocities are available in the
allvisit file at https://www.sdss4.org/dr17/irspec/spectro_data/.

## References

  * Abdurro’uf et al. (2022) Abdurro’uf et al., 2022, ApJS, 259, 35
  * Aitken (1964) Aitken R. G., 1964, The Binary Stars. Dover Publications
  * Bashi et al. (2023) Bashi D., Mazeh T., Faigler S., 2023, MNRAS, 522, 1184
  * Bate (2012) Bate M. R., 2012, MNRAS, 419, 3115
  * Bonnell & Bastien (1992) Bonnell I., Bastien P., 1992, ApJ, 401, 654
  * Bressan et al. (2012) Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127
  * Chen et al. (2020) Chen X., Wang S., Deng L., de Grijs R., Yang M., Tian H., 2020, ApJS, 249, 18
  * Christy et al. (2023) Christy C. T., et al., 2023, MNRAS, 519, 5271
  * Clarke (1996) Clarke C. J., 1996, MNRAS, 283, 353
  * Claytor et al. (2024) Claytor Z. R., van Saders J. L., Cao L., Pinsonneault M. H., Teske J., Beaton R. L., 2024, ApJ, 962, 47
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Duquennoy & Mayor (1991) Duquennoy A., Mayor M., 1991, A&A, 248, 485
  * El-Badry et al. (2018) El-Badry K., et al., 2018, MNRAS, 476, 528
  * Geller et al. (2021) Geller A. M., Mathieu R. D., Latham D. W., Pollack M., Torres G., Leiner E. M., 2021, AJ, 161, 190
  * Halbwachs (1987) Halbwachs J. L., 1987, A&A, 183, 234
  * Hall (1976) Hall D. S., 1976, in Fitch W. S., ed., Astrophysics and Space Science Library Vol. 60, IAU Colloq. 29: Multiple Periodic Variable Stars. p. 287, doi:10.1007/978-94-010-1175-4_15
  * Han et al. (2020) Han Z.-W., Ge H.-W., Chen X.-F., Chen H.-L., 2020, Research in Astronomy and Astrophysics, 20, 161
  * Jackson & Jeffries (2010) Jackson R. J., Jeffries R. D., 2010, MNRAS, 402, 1380
  * Jaschek & Ferrer (1972) Jaschek C., Ferrer O., 1972, PASP, 84, 292
  * Jayasinghe et al. (2018) Jayasinghe T., et al., 2018, MNRAS, 477, 3145
  * Jayasinghe et al. (2019a) Jayasinghe T., et al., 2019a, MNRAS, 485, 961
  * Jayasinghe et al. (2019b) Jayasinghe T., et al., 2019b, MNRAS, 486, 1907
  * Jayasinghe et al. (2020) Jayasinghe T., et al., 2020, MNRAS, 491, 13
  * Jayasinghe et al. (2021) Jayasinghe T., et al., 2021, MNRAS, 503, 200
  * Koch & Hrivnak (1981) Koch R. H., Hrivnak B. J., 1981, AJ, 86, 438
  * Kollmeier et al. (2017) Kollmeier J. A., et al., 2017, arXiv e-prints, p. arXiv:1711.03234
  * Kouwenhoven et al. (2009) Kouwenhoven M. B. N., Brown A. G. A., Goodwin S. P., Portegies Zwart S. F., Kaper L., 2009, A&A, 493, 979
  * Kratter & Matzner (2006) Kratter K. M., Matzner C. D., 2006, MNRAS, 373, 1563
  * Kroupa (1995a) Kroupa P., 1995a, MNRAS, 277, 1491
  * Kroupa (1995b) Kroupa P., 1995b, MNRAS, 277, 1507
  * Leiner et al. (2022) Leiner E. M., Geller A. M., Gully-Santiago M. A., Gosnell N. M., Tofflemire B. M., 2022, ApJ, 927, 222
  * Lorimer (2008) Lorimer D. R., 2008, Living Reviews in Relativity, 11, 8
  * Lurie et al. (2017) Lurie J. C., et al., 2017, AJ, 154, 250
  * Majewski et al. (2017) Majewski S. R., et al., 2017, AJ, 154, 94
  * Marigo et al. (2013) Marigo P., Bressan A., Nanni A., Girardi L., Pumo M. L., 2013, MNRAS, 434, 488
  * Marks & Kroupa (2011) Marks M., Kroupa P., 2011, MNRAS, 417, 1702
  * Mazeh & Goldberg (1992) Mazeh T., Goldberg D., 1992, ApJ, 394, 592
  * Mazzola et al. (2020) Mazzola C. N., et al., 2020, MNRAS, 499, 1607
  * McDonald & Clarke (1995) McDonald J. M., Clarke C. J., 1995, MNRAS, 275, 671
  * Meibom & Mathieu (2005) Meibom S., Mathieu R. D., 2005, ApJ, 620, 970
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Oelkers et al. (2018) Oelkers R. J., et al., 2018, AJ, 155, 39
  * Phillips et al. (2024) Phillips A., Kochanek C. S., Jayasinghe T., Cao L., Christy C. T., Rowan D. M., Pinsonneault M., 2024, MNRAS, 527, 5588
  * Pinsonneault & Stanek (2006) Pinsonneault M. H., Stanek K. Z., 2006, ApJ, 639, L67
  * Press (1997) Press W. H., 1997, in Unsolved Problems in Astrophysics. pp 49–60 (arXiv:astro-ph/9604126), doi:10.48550/arXiv.astro-ph/9604126
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Reinhold et al. (2023) Reinhold T., Shapiro A. I., Solanki S. K., Basri G., 2023, A&A, 678, A24
  * Sana et al. (2011) Sana H., Le Bouquin J. B., De Becker M., Berger J. P., de Koter A., Mérand A., 2011, ApJ, 740, L43
  * Simonian et al. (2019) Simonian G. V. A., Pinsonneault M. H., Terndrup D. M., 2019, ApJ, 871, 174
  * Tokovinin (2000) Tokovinin A. A., 2000, A&A, 360, 997
  * Torres et al. (2010) Torres G., Andersen J., Giménez A., 2010, A&ARv, 18, 67
  * Tout (1991) Tout C. A., 1991, MNRAS, 250, 701
  * Trimble (1990) Trimble V., 1990, MNRAS, 242, 79
  * Van Eylen et al. (2016) Van Eylen V., Winn J. N., Albrecht S., 2016, ApJ, 824, 15
  * Verbunt (1993) Verbunt F., 1993, ARA&A, 31, 93
  * Verbunt & Phinney (1995) Verbunt F., Phinney E. S., 1995, A&A, 296, 709
  * Wang & Han (2012) Wang B., Han Z., 2012, New Astron. Rev., 56, 122


# Population-based wind farm monitoring based on a spatial autoregressive
approach

W. Lin  K. Worden and E.J. Cross  Dynamics Research Group, University of
Sheffield, Sheffield S1 3JD, UK<EMAIL_ADDRESS>
###### Abstract

An important challenge faced by wind farm operators is to reduce operation and
maintenance cost. Structural health monitoring provides a means of cost
reduction through minimising unnecessary maintenance trips as well as
prolonging turbine service life. Population-based structural health monitoring
can further reduce the cost of health monitoring systems by implementing one
system for multiple structures (i.e. turbines). At the same time, shared data
within a population of structures may improve the predictions of structural
behaviour. To monitor turbine performance at a population/farm level, an
important initial step is to construct a model that describes the behaviour of
all turbines under normal conditions. This paper proposes a population-level
model that explicitly captures the spatial and temporal correlations (between
turbines) induced by the wake effect. The proposed model is a Gaussian
process-based spatial autoregressive model, named here a GP-SPARX model. This
approach is developed since (a) it reflects our physical understanding of the
wake effect, and (b) it benefits from a stochastic data-based learner. A case
study is provided to demonstrate the capability of the GP-SPARX model in
capturing spatial and temporal variations as well as its potential
applicability in a health monitoring system.

##  1 Introduction

In light of the current climate situation and power crisis, urgent actions are
needed to accelerate the energy transition from fossil fuels to renewables.
Wind power is a key type of renewable source and currently constitutes 9% of
the world’s electricity. Although the offshore wind sector has been growing
rapidly over the past year, the rate of offshore wind installations needs to
increase by a factor of four in order to meet the objectives of net zero by
2050 [1]. The rapid expansion in offshore wind power leads to significant
challenges as well as opportunities. Given that the cost of operation and
maintenance in offshore wind is significantly higher than its onshore
counterpart, the development of predictive maintenance strategies for wind
farms becomes more crucial than ever.

Population-based structural health monitoring (PBSHM) provides a means of
developing predictive maintenance strategies while accounting for the
interactions between the turbines across a wind farm. In contrast to the
single-structure approach in conventional data-based structural health
monitoring (SHM), PBSHM focusses on data sharing across a population of
structures to improve the predictability of models [2, 3, 4, 5, 6, 7, 8, 9].
When sharing data between structures, one has to acknowledge that no matter
how similar the structures are in terms of geometry, material and topology,
the data collected from different structures are likely to be different
because of environmental and operational variations (EOV) and manufacturing
tolerances. In PBSHM, there are currently two approaches to address this
inter-structure difference. One of them is based on domain adaptation, which
removes inter-structure difference by homogenising data from various
structures [10, 11]. The other approach is to model the inter-structure
correlations with the help of the physical understanding of EOV [12, 9], which
is the focus of this study.

The monitoring of turbine behaviours across an offshore wind farm is a unique
problem in many aspects. Firstly, the turbines in a wind farm are of the same
model, meaning that they are expected to behave in the same way under the same
environmental and operational conditions. Secondly, the turbines are
positioned close to each other, sharing the same environment. Currently, the
most commonly used wind turbines are the horizontal-axis ones, which create a
wake region downstream of each turbine rotor. The wind environment around each
turbine is affected by the existence of wakes, therefore, the turbines in a
wind farm can affect each other’s behaviour through the environment. That is,
a spatial pattern exists in the wind environment and, thus, turbine behaviours
across a wind farm. This spatial pattern changes with time-varying wind
directions, resulting in spatio-temporal correlations across the farm.
Capturing these correlations may enable a monitoring method to predict which
turbines are more likely to be damaged in a wind farm as well as to maximise
the overall power production across the farm.

In this paper, a two-level method is proposed to model the spatio-temporal
correlations across a wind farm. The first level consists of several Gaussian
process-based spatial autoregressive models (GP-SPARX), each capturing the
spatial pattern corresponding to a specific wind direction. In the second
level, (simple) criteria are designed that determine when to switch on/off the
first-level models. A case study of a simulated wind farm is used to
demonstrate the capability of the proposed switching GP-SPARX model, which
also sheds light on possible ways to improve this method.

The layout of this paper is as follows. Section 2 explains the motivation of
the proposed two-level method, by defining the wind farm modelling problem in
detail. Section 3 provides descriptions for the proposed switching GP-SPARX
model. Section 4 presents the case study and the prediction results obtained
using the proposed method, with the concluding remarks given in Section 5.

##  2 Wind farm modelling as a two-level problem

For offshore wind farms with horizontal-axis turbines, the wake effect has a
large impact on turbine power production. Wake is the region behind a rotor
where the wind field is disturbed by the rotating turbine, often characterised
by decreased mean wind speed and increased turbulence [13]. A downwind turbine
placed in the wake of an upwind turbine tends to produce less power than in
free-stream conditions. The wake-induced power loss can amount to 80% of the
free-stream power in compactly spaced wind farms such as Lillgrund [14].

The existence of wakes gives rise to a spatial pattern across a wind farm. As
the wind approaches a first-row turbine, the turbine extracts optimal power
from the wind, leaving a wake downstream. For a wake-shadowed turbine, both
its power output and subsequent wake are affected by its upstream wake
environment. As a result, the turbine behaviours along a path of wake
progression are spatially autoregressive. The correlation between neighbouring
turbines along a wake path is highly nonlinear because of the complex wind-
turbine interactions – the wind drives the rotation in a turbine, which then
disrupts downstream wind flow. To this end, the farm-wise spatial pattern,
determined by wake progression paths, is nonlinear and spatial autoregressive
by nature. A Gaussian process-based spatial autoregressive model is,
therefore, suitable to capture this pattern.

The farm-wise spatial pattern varies with time. In a wind farm, although the
turbines are programmed to face the incoming wind, they might face slightly
different directions due to local turbulence. Since the paths of wake
progression are highly sensitive to turbine angles, a change in any turbine
angle will result in a different spatial pattern. Therefore, the time-varying
turbine angles give rise to time-varying spatial patterns. A second-level
model is needed to determine how the spatial autoregressive wake progression
changes in response to turbine angles.

In this respect, the modelling of a wind turbine array can be considered a
two-level problem: the first level concerns the spatial autoregressive wake
progression across the farm; the second level captures how the spatial
autoregressive pattern changes with time-varying turbine angles.

##  3 The switching GP-SPARX model

###  3.1 GP-SPARX

Given a fixed set of wake paths, the first level of the model captures the
nonlinear, spatial autoregressive variations across a wind farm through a GP-
SPARX model [15, 16, 17],

$\displaystyle u(s,t)=f\Bigl{(}u_{\infty}(t),\ w(i,s)u(i,t)\Bigr{)},\quad
i=1,\dots,S$ (1) $\displaystyle
f\left(\cdot\right)\sim\mathcal{GP}\left(\mathbf{0},k_{SE}\left(\cdot,\cdot\right)\right)$

where $u(s,t)$ denotes the wind speed at location $s$ and time $t$.
$u_{\infty}(t)$ is the exogenous input that represents the free-stream wind
speed. The spatial correlation between turbines is indicated in the weighted
autoregressive term $w(i,s)u(i,t)$; for every turbine $i$ in the wind farm of
$S$ turbines, the weight $w(i,s)=1$ when the output location $s$ is in the
wake of turbine $i$, otherwise $w(i,s)=0$. It ensures that any non-zero
autoregressive term corresponds to an upstream neighbour that provides wake
shadowing effects to turbine $s$. To capture the nonlinearity in wake
progression, the function $f\left(\cdot\right)$ is assumed to be a Gaussian
process with a mean zero and a squared exponential covariance function. For
more details on GP regression, the reader is referred to [18].

###  3.2 The switching model

As the turbine angles change, different wake patterns result. The idea here is
to train a limited number of GP-SPARX models, each corresponding to a spatial
pattern, and use a higher-level model to switch between these GP-SPARX models
when applied to a testing data set. In this preliminary study, both the
selection of training patterns and the model switching are done manually, as a
cost-effective way of evaluating the method before developing more
complicated/automated alternatives.

$u(s,t)=\begin{cases}f_{1}\Bigl{(}u_{\infty}(t),\
w(i,s,\phi)u(i,t)\Bigr{)}&\underline{\phi}_{1}\leq\phi<\bar{\phi}_{1}\\\
\qquad\qquad\vdots\\\ f_{n}\Bigl{(}u_{\infty}(t),\
w(i,s,\phi)u(i,t)\Bigr{)}&\underline{\phi}_{n}\leq\phi<\bar{\phi}_{n}\end{cases},\
i=1,\dots,S$ (2)

In Equation 2, the output $y$ is obtained by switching between $n$ trained GP-
SPARX models $f_{1}(\cdot),\dots,f_{n}(\cdot)$, according to criteria based on
the wind direction $\phi$. Note that both the function $f(\cdot)$ and the
weight $w(i,s,\phi)$ change with wind direction. $\underline{\phi}$ and
$\bar{\phi}$ denote the lower and upper bounds of $\phi$, respectively. It
indicates that a specific spatial pattern only applies to a certain range of
wind directions.

##  4 A case study of a simulated wind farm

To evaluate the efficacy of the suggested model form, data from a simulated
wind farm are used, such that we can assume no inadequate prediction is caused
by actual anomalies that may be present in operational data.

The simulated wind farm is shown in Figure 1, where the wake patterns
associated with two incoming wind directions are illustrated. Each blue dash-
dotted line indicates a wake path, along which the wake propagates from an
upstream turbine to a downstream turbine. The wake path is determined by both
the free-stream wind direction and the directions of local air vortices.

The corresponding spatial patterns of mean wind speed are shown in Figure 2.
The wind speeds are simulated using a third-order polynomial NARX model [17],
which provides results that are in line with the most commonly known
analytical wake models [13].

(a)

(b)

Figure 1: Wake patterns in various wind directions.

(a)

(b)

Figure 2: Mean wind speed map in various wind directions.

###  4.1 Training and testing data

Figure 3 demonstrates the training and testing data coverage in terms of wind
direction $\phi$. In this case study, four wake patterns, associated with four
different wind directions (blue arrows in Figure 3), are used to train four
first-level GP-SPARX models. The testing data, on the contrary, covers the
entire range of wind directions, shown as the orange circle arrow in Figure 3.
Testing predictions are made using the four trained first-level models, by
switching between them at predefined boundaries (black dashed lines in Figure
3). For simplicity, the boundaries are set as the midpoints between the
training angles.

Figure 3: Training and testing data coverage in wind direction.

###  4.2 Results

Previously in [17], the authors have demonstrated that the GP-SPARX model is
able to capture the nonlinear correlations across a wind farm with
satisfactory accuracy and outperform deterministic spatial NARX and ARX models
in accuracy. To carry the study forward, the focus of this paper is to
evaluate the efficacy of the second-level switching model in predicting the
time-varying spatial patterns across a wind farm.

In Figure 4, model prediction errors are presented in polar coordinates to
demonstrate the angles at which the two-level model provides satisfactory
accuracy or requires improvement. Two main observations can be obtained from
this error map. Firstly, high errors occur in regions away from the training
data. With the training angles indicated as black dashed lines in Figure 4, it
is seen that the predictions at training regions are relatively accurate. In
situations when a considerable amount of error is obtained close to a training
region, the error reduces rapidly as one approaches the training region. Away
from the training regions, the model starts to extrapolate (in terms of wind
direction), thus, there is a tendency for reduced accuracy. Secondly,
prediction accuracy varies in the regions away from the training data, giving
rise to spikes of high errors. These spikes indicate the wind directions (i.e.
wake patterns) to be included as part of the training data, to improve the
overall accuracy.

Figure 4: Map of squared errors in polar coordinates.

##  5 Conclusions

This paper demonstrates a population-based method that focusses on modelling
the correlations between structures, for the application of wind farm
monitoring. One of the main themes of the paper is to encourage the use of
physical understanding to simplify modelling approaches. For example, the
understanding of wake paths provides fixed spatial lags for model training,
which circumvents the optimisation problem of spatial lag selection. With the
rapid growth of computing power, it is crucial not to misspend the
computational resources when building data-based models. One way to achieve
this is to create parsimonious models by exploring physics. In the context of
wind farm monitoring, the proposed two-level model is potentially able to
predict the spatio-temporal correlations across a wind farm, given only one
spatial reference (i.e. weather station measurements).

This paper provides a preliminary study using an unrefined model form,
intending to evaluate the efficacy of the proposed method as well as
discovering possible ways to improve the model. It is concluded that a
switching GP-SPARX model may be able to accurately capture the time-varying
spatial patterns in a wind farm if the training regions and switching
boundaries are chosen correctly. Hence, a possible next step is to use an
automatic switching model, such as a treed GP [19], which learns from data how
to select the training regions.

## Acknowledgements

The authors would like to acknowledge the support of the EPSRC, particularly
through grant reference numbers EP/R004900/1, EP/S001565/1 and EP/R003645/1.
For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

## References

  * [1] Lee J and Zhao F 2022 Global wind report 2022 Tech. rep. Global Wind Energy Council
  * [2] Worden K, Cross E, Dervilis N, Papatheou E and Antoniadou I 2015 IFAC 48(21) 1–17
  * [3] Bull L, Gardner P, Gosliga J, Maguire A, Campos C, Rogers T, Haywood-Alexander M, Dervilis N, Cross E and Worden K 2020 Towards population-based structural health monitoring, Part I: Homogeneous populations and forms Proceedings of IMAC XXXVIII – the 38th International Modal Analysis Conference, Houston, TX
  * [4] Gosliga J, Gardner P, Bull L, Dervilis N and Worden K 2020 Towards population-based structural health monitoring, Part II: Heterogeneous populations and structures as graphs Proceedings of IMAC XXXVIII – the 38th International Modal Analysis Conference, Houston, TX
  * [5] Gosliga J, Gardner P, Bull L, Dervilis N and Worden K 2020 Towards population-based structural health monitoring, Part III: Graphs, networks and communities Proceedings of IMAC XXXVIII – the 38th International Modal Analysis Conference, Houston, TX
  * [6] Gardner P and Worden K 2020 Towards population-based structural health monitoring, Part IV: Heterogeneous populations, matching and transfer Proceedings of IMAC XXXVIII – the 38th International Modal Analysis Conference, Houston, TX
  * [7] Wickramarachchi C, Brennan D, Lin W, Harvey D, Cross E and Worden K 2021 Towards population-based structural health monitoring, Part V: Network and Databases Proceedings of IMAC-XXXIX – the 39th International Modal Analysis Conference
  * [8] Worden K 2020 Towards population-based structural health monitoring, Part VI: Structures as geometry Proceedings of IMAC XXXVIII – the 38th International Modal Analysis Conference, Houston, TX
  * [9] Lin W, Worden K, Maguire A and Cross E 2020 Towards population-based structural health monitoring, Part VII: EoV fields: environmental mapping Proceedings of IMAC XXXVIII – the 38th International Modal Analysis Conference, Houston, TX
  * [10] Gardner P, Bull L, Dervilis N and Worden K 2022 Journal of Civil Structural Health Monitoring
  * [11] Poole J, Gardner P, Dervilis N, Bull L and Worden K 2022 On the application of partial domain adaptation for PBSHM Proceedings of EWSHM 2022 – European Workshop on Structural Health Monitoring, Palermo, Italy vol 270
  * [12] Lin W, Worden K, Maguire A and Cross E in press Data-Centric Engineering
  * [13] Archer C, Vasel-Be-Hagh A, Yan C, Wu S, Pan Y, Brodie J and Maguire A 2018 Applied Energy 226 1187–1207
  * [14] Dahlberg J A 2009 Assessment of the Lillgrund windfarm: Power performance and wake effects Tech. rep. Vattenfall Vindkraft AB
  * [15] Worden K, Becker W, Rogers T and Cross E 2018 Mechanical Systems and Signal Processing 104 188–223
  * [16] Rogers T 2019 Towards Bayesian System Identification: With Application to SHM of Offshore Structures Ph.D. thesis Department of Mechanical Engineering, University of Sheffield
  * [17] Lin W, Worden K and Cross E 2022 A spatial autoregressive approach for wake field prediction across a wind farm Proceedings of EWSHM 2022 – European Workshop on Structural Health Monitoring, Palermo, Italy vol 270
  * [18] Rasmussen C and Williams C 2006 Gaussian Processes for Machine Learning (Cambridge MA: The MIT Press)
  * [19] Worden K, Cross E and Brownjohn J 2013 Surrogate-Based Modeling and Optimization 337–358