# A weakly-parametric approach to stochastic background inference in LISA

Federico Pozzoli0009-0009-6265-584X  Dipartimento di Scienza e Alta
Tecnologia, Università dell’Insubria, via Valleggio 11, I-22100 Como, Italy
<EMAIL_ADDRESS>Riccardo Buscicchio0000-0002-7387-6754  Dipartimento di
Fisica “G. Occhialini”, Università degli Studi di Milano-Bicocca, Piazza della
Scienza 3, 20126 Milano, Italy  Institute for Gravitational Wave Astronomy &
School of Physics and Astronomy, University of Birmingham, Birmingham, B15
2TT, UK  INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126
Milano, Italy  INAF – Osservatorio Astronomico di Brera, via E. Bianchi 46,
23807, Merate, Italy  Christopher J. Moore0000-0002-2527-0213  Institute for
Gravitational Wave Astronomy & School of Physics and Astronomy, University of
Birmingham, Birmingham, B15 2TT, UK  Francesco Haardt0000-0003-3291-3704
Dipartimento di Scienza e Alta Tecnologia, Università dell’Insubria, via
Valleggio 11, I-22100 Como, Italy  INFN, Sezione di Milano-Bicocca, Piazza
della Scienza 3, 20126 Milano, Italy  INAF – Osservatorio Astronomico di
Brera, via E. Bianchi 46, 23807, Merate, Italy  Alberto
Sesana0000-0003-4961-1606  Dipartimento di Fisica “G. Occhialini”, Università
degli Studi di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy
INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy
INAF – Osservatorio Astronomico di Brera, via E. Bianchi 46, 23807, Merate,
Italy

###### Abstract

Detecting stochastic gravitational wave backgrounds (SGWBs) with The Laser
Interferometer Space Antenna (LISA) is among the mission science objectives.
Disentangling SGWBs of astrophysical and cosmological origin is a challenging
task, further complicated by the noise level uncertainties. In this study, we
introduce a Bayesian methodology to infer upon SGWBs, taking inspiration from
Gaussian stochastic processes. We investigate the suitability of the approach
for signal of unknown spectral shape. We do by discretely exploring the model
hyperparameters, a first step towards a more efficient transdimensional
exploration. We apply the proposed method to a representative astrophysical
scenario: the inference on the astrophysical foreground of Extreme Mass Ratio
Inspirals, recently estimated in [1]. We find the algorithm to be capable of
recovering the injected signal even with large priors, while simultaneously
providing estimate of the noise level.

##  I Introduction

The Laser Interferometer Space Antenna (LISA) [2] is a groundbreaking mission
for the detection of Gravitational Waves (GW) from space. Throughout the
nominal $4$ years of duration of operations, an abundance of individually
resolvable sources will be available. A primary focus of LISA will be the
detection of several tens of massive black-hole binary (MBHBs) mergers [3, 4],
spanning masses from approximately $10^{4}$ to $10^{7}$ times that of the Sun.
These mergers are anticipated to have high Signal-to-Noise Ratios (SNRs),
reaching up to several thousands. Additionally, LISA will identify several to
a few hundred extreme mass ratio inspirals (EMRIs) annually [5], along with
tens of thousands of Galactic white dwarf binaries [6].

Besides individual sources, the Universe is pervaded by a stochastic
gravitational wave background (SGWB), which arises from an incoherent
superposition of GWs originating from numerous unresolved or weak sources. At
nano-Hz frequencies, the first evidence for an unresolved GW signal of
astrophysical and/or cosmological origin [7, 8] has been recently been
reported by the European pulsar timing array [9], NANOGrav [10], the Parkes
Pulsar timing array [11] and the Chinese pulsar timing array [12]. The SGWB in
LISA may come from both cosmological and astrophysical sources. Cosmologically
[13], these waves could stem from primordial quantum fluctuations [14],
potentially amplified during cosmic inflation. Furthermore, processes such as
first-order phase transitions [15, 16], interacting cosmic strings [17], and
primordial black holes [18] may also contribute to stochastic GW emissions.
Astrophysically, the Galactic foreground [19] from double white dwarf binaries
(DWDs) is expected to dominate, particularly in the $0.5-3$ mHz range.
Additionally, both stellar-origin black-hole binaries (SOBHBs) [20] and EMRIs
have the potential to generate non-negligible backgrounds [1, 21].

The persistence of the SGWB in the collected data ties inextricably to the
estimation of instrumental noise. Various efforts have been made to construct
a null channel in LISA, that is a datastream effectively insensitive to GWs:
however, in realistic instrumental setups, this won’t be possible across the
whole sensitivity band [22].

Disentangling different SGWB components and instrumental noise provides a
challenging data-analysis task. Many approaches have been already proposed
using approaches based on stochastic-sampling algorithm such as Markov Chain
Monte Carlo (MCMC)[23] or nested sampling [24]. In particular, noise and
signal reconstruction leveraging templated inference for both [25, 16], or
just one [26, 27, 28] have been explored. In this study, we introduce a
weakly-parametric approach by implementing an inference based on stochastic
Gaussian processes, which can be applied in various contexts with great
flexibility in the SGWB spectral shape.

The paper is organized as follows. In Section II we provide a brief overview
of the putative astrophysical SGWBs potentially observable with LISA. In
Section III we describe the underlying assumptions of our analysis. In Section
IV a detailed exposition of the core statistical model is presented. In
Section V, we report our parameter estimation results on a set of simulated
SGWB signals. Finally, in Section VI, we summarize our findings and outline
the prospects for future developments towards a highly flexible inference on
SGWBs. Unless otherwise specified, summation over repeated indices is assumed
throughout.

##  II Astrophysical Backgrounds

Astrophysical SGWBs are expected to form from the superposition of GWs from
compact binary systems. However, other putative astrophysical SGWBs are
expected to arise, e.g. from asymmetric supernova explosions [29] or rapidly
rotating neutron stars [30].

Signals typically manifest as either isolated and faint, falling below the
detection threshold, or as overlapping ones. They pile up incoherently, hence
resulting in a “confused” timeseries, whose morphology is too complex to tell
individual signals apart. Moreover, the different astrophysical (and
cosmological) SGWBs will overlap in the LISA data making it difficult to study
them, individually [31]. However, the different spectral shapes of each
background component and the different sky distributions of Galactic
(localised near the plane/buldge of the Milky Way) versus extragalactic (more
isotropic) backgrounds can be used to separate them. The dominant contributors
to astrophysical SGWBs in the LISA sensitivity band are expected to originate
mainly from three populations. We summarize briefly their properties below.

_Galactic DWDs._ The white dwarf binary Galactic population is expected to be
formed by around $\tilde{\times}10^{7}$ systems [32]. A significant fraction
(up to $\sim 10\%$) contributes with GW signals in the LISA band, i.e. at GW
frequencies above $0.1\rm mHz$ [6]. Thousands of such highly monochromatic
systems will be resolved individually, some of which are already observed
through their electromagnetic counterparts; these are known as verification
binaries [33, 34, 35]. Ongoing and future missions such as GAIA [36] and the
Vera Rubin Observatory [37] have the potential to unveil up to a few thousand
more DWD systems. The largest fraction of the DWDs will remain unresolved
–only $\sim 0.1\%$ of DWDs are going to be resolved individually [38]– and
form a stochastic signal, commonly referred to as _confusion noise_. This is
expected to be larger than LISA instrumental noise the frequency range of
$0.5-3$ mHz. The bulk of such foreground arises from systems located toward
the Galactic center, hence the confusion noise is expected to be highly
anisotropic. Leveraging the SGWB sky distribution and the induced signal
modulation due to LISA motion, is an opportunity to characterize the Galactic
morphology [39], i.e. its disk, bulge, halo, and streams.

_SOBHBs._ SOBHBs will be detected by LISA during their early inspiral phase
[40]. A portion of these binaries will enter the final years of their inspiral
and eventually transition into the frequency range detectable by LIGO, Virgo,
and Kagra (LVK), thereby allowing for multiband gravitational-wave astronomy
[41, 42]. Based on current LVK constraints, it is anticipated that at least a
few of these binaries will be individually detected and characterized [43, 44,
45], while the majority will contribute to the formation of a SGWB. The
background is expected to be isotropic and power-law shaped with spectral
index $2/3$ for the energy density power spectrum (or equivalently, $-7/3$ for
the PSD), for quasi-circular binaries at leading order post-Newtonian
approximation [46]. Recently, accurate detectability of the SOBBHs background
has been demostrated [20], accounting for instrumental noise uncertainties and
coexistence of DWDs foreground.

_EMRIs._ A stellar mass object orbiting around a MBH constitute an EMRI. These
sources have not been observed through either electromagnetic or GW radiation,
yet. Nonetheless, GWs emitted by EMRIs are a primary focus for LISA, due to
their potential association with Fast Radio Bursts in the electromagnetic
spectrum [47]. Because of their very small mass ratios $\sim 10^{-6}-10^{-4}$,
EMRIs evolve slowly: a considerable number of these systems will persist in
the LISA band throughout the entire mission. A large number of orbital cycles
$\sim 10^{4}-10^{5}$ provides the opportunity to probe, in the test-particle
limit, highly curved spacetime close to supermassive black holes (SMBHs). More
recently, EMRIs have been proposed as potential targets for strong
gravitational lensing [48]. Most of the population will remain undetected,
hence a SGWB is expected to form. Unlike the SOBHBs population, the EMRI
background spectral shape prediction is uncertain due to (i) the spread of
each gravitational-wave signal over multiple harmonics of the orbital
frequency and (ii) the uncertainty related to the populations. A systematic
characterization of EMRI backgrounds, from astrophysically motivated model
populations developed in [5] was recently presented in [1]: the majority of
models yield a bright SGWB with SNRs ranging between tens and thousands,
contributing mostly to higher frequencies than the DWD background, namely
$1-10$ mHz.

##  III Stochastic signals in LISA

The LISA constellation is composed of three spacecrafts arranged in a nearly
equilateral configuration, respectively separated by $2.5\times 10^{6}\rm km$.
Spacecrafts will trail the Earth motion by an angle of about $20^{\circ}$,
lying on a plane inclined by $60$ degrees with respect to the orbital plane.
Distances between spacecrafts will be modulated by incident GWs, inducing
variations at the picometer level. Such variations will be detected by
monitoring the frequency (or phase, equivalently) of laser beams exchanged
between spacecrafts and comparing them to local reference lasers. The primary
source of noise will come from laser frequency fluctuations, many orders of
magnitude higher than the target sensitivity. To suppress it, six raw Doppler
measurements are suitably delayed and linearly combined in a set of time-delay
interferometric (TDI) variables [49]. Various TDI scheme exist in literature,
associated to different approximations of the LISA constellation orbits; in
this work we use TDI 1.5 combinations [50], assuming constant and unequal-
length LISA arms. A further linear combination of such TDI variables is used
to define approximately noise-orthogonal datastreams.

We perform inference on such data,
$\tilde{d}=(\tilde{A},\tilde{E},\tilde{T})$, the Fourier transforms of noise-
orthogonal TDI channels. We model the datastreams as the superposition of two
independent stochastic processes: the SGWB, the instrumental noise. In this
work we do not consider the contribution from the non-stationary Galactic
foreground, whose addition we leave for future study. At each frequency $f$
the correlation between data reads

$\left\langle\tilde{d}_{\alpha}(f)\tilde{d}^{*}_{\beta}(f)\right\rangle=S_{h,\alpha\beta}(f)+S_{n,\alpha\beta}(f),$
(1)

where $\alpha$ and $\beta$ denote each TDI channel. We assume throughout
isotropic, Gaussian, stationary SGWBs with equal-weighted and uncorrelated
polarizations.

The first component in the previous equation reads:

$S_{h,\alpha\beta}(f)=R_{\alpha\beta}(f)S_{h}(f).$ (2)

$S_{h}(f)$ represents the one-sided SGWB strain power spectral density (PSD)
[51], defined as

$\left\langle
h_{p}(f;\Omega)h^{*}_{p^{\prime}}(f^{\prime},\Omega^{\prime})\right\rangle\\!=\frac{1}{2}S_{h}(f)\delta(f-f^{\prime})\delta(\Omega,\Omega^{\prime})\delta_{pp^{\prime}}$
(3)

and $R_{\alpha\beta}$ is a $3\times 3$ matrix obtained by combining $M_{\rm
TDI}(f)$ ($3\times 6$) and the single-link response matrix $G(f,t_{0})$
($6\times 6$) as follows:

$R_{\alpha\beta}(f)=M_{\alpha i,\rm TDI}(f)G_{ij}(f,t_{0})M^{*}_{j\beta,\rm
TDI}(f).$ (4)

For completeness, we provide full derivation of $G(f,t_{0})$ in appendix A.

The spectrum $S_{n,\alpha\beta}$ represents the single-sided noise power
spectra in a single laser link. Following [26], we further assume no
correlations between individual links and the same PSD in each of them. Hence,
PSDs in each channel depend uniquely on the transfer matrix $M_{\rm TDI}$ and
on a common spectral model $S_{n}(f)$.

$S_{n,\alpha\beta}(f)=\frac{1}{2}S_{n}(f)M_{\alpha i,\rm TDI}(f)M_{i\beta,\rm
TDI}^{*}(f)$ (5)

Specifically, we include in $S_{n}$ secondary noise sources, namely the test
mass (TM) and the optical metrology system (OMS) noise. The overall PSD is
$S_{n}=2S_{n}^{\rm TM}+S_{n}^{\rm OMS}$ [22], where

$\displaystyle S_{n}^{\rm TM}\left(f\right)$
$\displaystyle=A^{2}\left[1+\left(\frac{0.4\rm
mHz}{f}\right)^{2}\right]\left[1+\left(\frac{f}{8\rm mHz}\right)^{4}\right]$
$\displaystyle\times\left(\frac{1}{2\pi fc}\right)^{2}\left(10^{-30}\rm
m^{2}/\rm s^{3}\right),$ (6)

and

$S_{n}^{\rm OMS}\\!\left(f\right)\\!=\\!P^{2}\\!\left[1\\!+\\!\left(\frac{2\rm
mHz}{f}\right)^{4}\right]\\!\\!\left(\frac{2\pi
f}{c}\right)^{2}\\!\\!\\!\\!\left(10^{-24}\rm m^{2}\\!\cdot\rm s\right)$ (7)

with $A=3$ and $P=15$ [52].

Through TDI channels signal and noise PSDs, the SNR for an SGWB is readily
defined as

$\mathrm{SNR}^{2}=T\sum_{\alpha}\int^{\infty}_{0}df\frac{S^{2}_{h,\alpha\alpha}}{S^{2}_{n,\alpha\alpha}}$
(8)

with $T$ the LISA observation time [53].

##  IV Inference model

In this study we only focus on the parameter estimation of SGWB and noise,
thus we assume that all non-stochastic GW sources have been perfectly
subtracted from the data. We choose to model both spectra – introduced in Eqs.
(2) and (5) – with a weakly-parametric model, which allows us to explore a
large family of spectral shapes. We do so by using expectation values of
Gaussian processes. A Gaussian process (GP) is a stochastic model, formally
describing distributions over functions $g(x)$. It is parameterized by a mean
function $m(x)$ and positive definite covariance function $k(x,x^{\prime})$.
Realizations from such process are denoted as follows

$g(x)\sim\mathcal{GP}(m(x),k(x,x^{\prime})),$ (9)

The defining properties of GPs is that for any finite subset
$X=\left\\{x_{1}\dots x_{n}\right\\}$ of the domain for $x$, the marginal
distribution is a multivariate Gaussian

$\bm{g}(X)\sim\mathcal{N}(\bm{\mu},\bm{\Sigma}),$ (10)

with mean vector and covariance matrix defined by $\mu=m(x_{i})$ and
$\Sigma_{ij}=k(x_{i},x_{j})$, respectively. Typically, GPs are used to
flexibly incorporate some observed data in the distribution over $\bm{g}(X)$
and make predictions on new domain points $X_{*}$ as $\bm{g}(X_{*})$. The
joint probability distribution $p(\bm{g}({X}_{*}),\bm{g}({X}))$ is a
multivariate normal distribution. In order to employ this formalism as a
regression model, one needs the conditional probability
$p(\bm{g}(X_{*})|\bm{g}({X}))$. Upon suitable marginalization, the conditional
distribution is also a multivariate normal distribution

$p(\bm{g}({X}_{*})|\bm{g}({X}))=\mathcal{N}(\bm{\mu}\left({{X}_{*}|{X}}\right),\bm{\Sigma}\left({{X}_{*}|{X}}\right)),$
(11)

where

$\displaystyle\bm{\mu}\left({X_{*}|X}\right)$
$\displaystyle=\bm{\mu}(X_{*})\,+$
$\displaystyle+\bm{\Sigma}\left({X_{*}|X}\right)\bm{\Sigma}\left({X,X}\right)^{-1}(\bm{g}({X})\\!-\\!\bm{\mu}({X}))$
(12)

and

$\displaystyle\bm{\Sigma}\left({X_{*}|X}\right)$
$\displaystyle=\bm{\Sigma}\left({X_{*},X_{*}}\right)+$
$\displaystyle+\bm{\Sigma}\left({X_{*},X}\right)\bm{\Sigma}^{-1}\left({X,X}\right)\bm{\Sigma}\left({X,X_{*}}\right),$
(13)

are specified by the single kernel function $k(x,x^{\prime})$. The matrix
$\bm{\Sigma}$ is also referred to as _kernel matrix_ , and it models the
covariance between each pair of its two arguments through the definition of
the bivariate function $k$. In order to be a valid covariance for the
multivariate Gaussian in Eq. (11), the kernel matrix must be symmetric and
positive definite. A variety of kernels are available in literature, to
capture different processes peculiarities. In this work, we consider the
radial basis kernel function (RBF) defined by

$\Sigma^{\rm RBF}_{ij}\left({X,Y}\right)=k_{\rm
RBF}(x_{i},y_{j})=\exp{\left(-\frac{|x_{i}-y_{j}|^{2}}{2\sigma^{2}}\right)},$
(14)

where $\sigma$ is a model positive hyperparameter.

Figure 1: Realizations from two EGP models, with fixed baseline and knots
amplitudes sampled from their respective parameter space. Dashed gray lines
represent the knot locations. EGP realizations are shown with solid lines,
while dashed red (blue) lines correspond to the fixed power-law (bump)
baseline.

In this study, we exploit GPs in an unconventional way: we use each GP
expectation value from Eq. (12) as a proposed spectrum, and set $\bm{g}({X})$
as free parameters. The GP covariance is not used in the inference model. The
dimensionality and domain location of $X$ are to be considered
hyperparameters.

Henceforth, we model the logarithm of noise and signal PSDs as two independent
EGPs over the logarithmic frequencies, i.e.
$X=\left\\{\log_{10}f_{1},\dots,\log_{10}f_{n}\right\\}$. We illustrate both
in Secs. IV.1 and IV.2, respectively. We will refer to the parameters $X$ and
to the mean function $m(X)$ as knots and baseline, respectively.

###  IV.1 Signal

We assume the response matrix $R_{\alpha\beta}(f)$ to be exactly known. Thus,
we need to specify $S_{h}(f)$, only. We use an EGP, with baseline and knots
left free to vary simultaneously. We do so to be able to capture both global
spectral shapes and fine structure in specific frequency regions: the choice
of baseline influences globally the spectral-shape proposals, while knots
control the fine-local structure. In Figure 1, we show the proposed EGP models
for two baseline families, i.e. power-laws and _bump_ spectra. The knots
variability can compensate easily for misidentification of the truth baseline
family. In our inference, we use a power-law baseline family for
$\log_{10}S_{h}(f)$, parametrized by logarithmic amplitude $\log_{10}A$ at a
reference frequency $f_{\star}=10^{-2.5}\rm Hz$ and slope $\gamma$. Thus,
$\bm{\mu}({X}_{*})$ in Eq. (12) reads
$\log_{10}A+\gamma\log_{10}(f/f_{\star})$. Fine-structure deviations with
respect to the baseline are parameterized through knots $\bm{\delta}^{h}$,
defined by

$\displaystyle g(X)_{k}-\mu({X})_{k}$
$\displaystyle=\log_{10}(S_{h}(f_{k})\cdot\delta^{h}_{k})-\log_{10}{S_{h}(f_{k})}$
$\displaystyle=\log_{10}\delta^{h}_{k}$ (15)

The number of knots is fixed in each inference, and we choose the associated
frequencies $f_{k}$ equally log-spaced within the data frequency range.
Overall the model is summarized by functions from the parametric family

$\displaystyle\log_{10}S_{h}(f;{\log_{10}A,\gamma,\bm{\delta}^{h}})=\log_{10}A\,+$
$\displaystyle+\gamma\log_{10}\left(f/f_{\star}\right)+k(f,f_{j})k({f_{j},f_{k}})^{-1}\log_{10}\delta^{h}_{k},$
(16)

###  IV.2 Noise

Figure 2: Noise realizations from GP and EGP model. Top panel: realizations
from GP. The gray shaded region represents the $90\%$ credible interval.
Bottom panel: Realizations from EGP. The gray area represents the prior $90\%$
credible interval. While GP allows for straightforward inclusion of observed
data, EGP smoothly explores a large set of spectral shapes, with the _data_
(i.e. knots amplitudes at chosen locations shown as grey dashed vertical
lines) effectively proposed as inference parameters.

We follow a similar approach to parameterize and infer upon the noise
uncertainties. Again, we assume a well known TDI transfer matrix $M_{\rm
TDI}$, so the modelling freedom is left for the single-link noise PSD.
Specifically, we introduce a set of parameter $\log_{10}\bm{\delta}^{n}$.
Contrary to the signal case discussed in Sec. IV.1, the baseline is fixed to
the nominal PSD level from ESA’s science requirement document [52]. Thus, the
final noise model reads:

$\displaystyle\log_{10}S_{n}(f;{{\bm{\delta}}^{n}})$
$\displaystyle=\log_{10}S^{\rm SciRD}_{n}+$
$\displaystyle+k(f,f_{j})k({f_{j},f_{k}})^{-1}\log_{10}\delta^{n}_{k},$ (17)

With this parameterization, we augment the noise reference model with
considerable flexibility to vary across a large functional space. In all our
inference the injected noise in the data is generated according to the fixed
baseline, so we anticipate recovering zero values for
$\log_{10}\bm{\delta}^{n}$. In Fig.2, we show a comparison between GP- and
EGP-like instrumental noise modelling, illustrating the different degree of
flexibility each one exhibits.

###  IV.3 Likelihood

Inference is performed simultaneously on SGWB and instrumental noise. We
construct joint posteriors on parameters $\bm{\theta}$

$p(\bm{\theta}|\tilde{d})\propto\mathcal{L}(\tilde{d}|\bm{\theta})\pi(\bm{\theta})$
(18)

through stochastic sampling of the likelihood
$\mathcal{L}(\tilde{d}|\bm{\theta})$ under chosen priors $\pi(\bm{\theta})$.
To do so, we use Balrog, a large codebase for simulation and inference on LISA
signals. In this study, we use it alongside a nested sampling algorithm as
implemented in Nessai [54] to obtain each posterior and evidence. Data are
distributed according to the Gaussian likelihood [55]

$\displaystyle\log\mathcal{L}(\tilde{d}|\bm{\theta})\\!=$ (19)
$\displaystyle=-\frac{1}{2}\displaystyle\sum_{\alpha,\beta}\displaystyle\sum_{f=f_{\rm
min}}^{f_{\rm
max}}\frac{\tilde{d}_{\alpha}(f)\tilde{d}^{*}_{\beta}(f)}{S_{n,\alpha\beta}(f,\bm{\delta}^{n})+S_{h,\alpha\beta}(f,\log
A,\gamma,\bm{\delta}^{h})}\Delta f,$

with
$\bm{\theta}=(\log_{10}A,\gamma,\log_{10}\bm{\delta}^{h},\log_{10}\bm{\delta}^{n})$.
We choose uniform priors for each parameter in $\bm{\theta}$, and the
following prior ranges:

  * •

Power-law amplitude $\log_{10}A$: $[-70,-35]$

  * •

Power-law slope $\gamma$: $[-5,5]$

  * •

Signal knots amplitudes $\log_{10}\bm{\delta}^{h}$: $\left[-2,2\right]$

  * •

Noise knots amplitudes $\log_{10}\bm{\delta}^{n}$: $[-0.5,0.5]$

Additionally, we will consider a prior uniform in $\left[-5,5\right]$ for
$\log_{10}\bm{\delta}^{h}$, instrumental to the results discussion in Sec. V.

We highlight once again that the number of knots $\bm{\delta}^{h}$ and
$\bm{\delta}^{n}$ can be further optimized, as much as the remaining
hyperparameters (i.e. the kernel length scale, its functional form, and the
knots locations). This represent a potential high-degree of flexibility to our
algorithm, which can be leveraged to find unexpected spectral features upon
inference. In this study, we focus on a finite set of kernel lengths $\sigma$
and number of knots. Noise parameter estimation is instead fixed to three
knots, given the baseline choice matching the noise model used for injection
in the data.

In absence of a trans-dimensional sampling framework [56] in our codebase, our
study is constrained to perform inference for each choice of hyperparameters,
and select the most suitable via Bayesian model selection. In order to compare
the different choices, we use the log-Bayes factor

$\log_{10}\mathcal{B}^{\sigma,n}_{\sigma^{\prime},n^{\prime}}=\log_{10}\mathcal{Z}_{\sigma,n}-\log_{10}\mathcal{Z}_{\sigma^{\prime},n^{\prime}},$
(20)

where $(\sigma,n)$ and $(\sigma^{\prime},n^{\prime})$ are labels identifying
competing models with different kernel length scale and number of knots,
respectively, and

$\mathcal{Z}=\displaystyle\int
d\bm{\theta}\mathcal{L}(\tilde{d}|\bm{\theta})\pi(\bm{\theta})$ (21)

is the Bayesian evidence. In Sec. V we present our results, together with a
map of evidences: they are to be interpreted as the model marginal
likelihoods, conditioned on hyperparameters, only. This is a core building
block towards a trans-dimensional exploration of the full EGP model for SGWB
inference and its inclusion in a larger global fit scheme [57, 58].

##  V Inference Results

In this section, we present parameter estimation results performed using the
model introduced in Sec. IV. First, we use a simple toy model to control the
correct recovery of a SGWB with known spectral shape (Sec. V.1). Then, we
apply our formalism to the parameter estimation of an SGWB from a population
of unresolved EMRIs (Sec. V.2). Throughout the analysis, LISA data are
simulated for 4 years of observation time. We proceed with the assumption that
non-stochastic GWs sources and glitches have been subtracted perfectly.
Additionally, we consider the absence of data gaps, resulting in the
generation of an idealized residual dataset.

###  V.1 Toy Model

We conduct a test retrieving an EGP model with a power-law baseline and non-
zero knot amplitudes.

Figure 3:  Spectral reconstruction of two toy model SGWBs. Injected SGWBs are
denoted by black solid lines. They have SNR of 209 and 5, respectively. This
places them above and below the nominal powerlaw sensitivity curve (solid blue
line) to an SGWB with SNR of 10 after 4 years of observation. Noise injection
and prior 90$\%$ confidence intervals are shown as gray solid and purple
dashed lines, respectively. Posterior median and $90\%$ confidence intervals
on each signal –analysed independently– are denoted by red dashed lines and
light-red shared areas. The EGP flexibility captures features in the PSD
shape.

Although the injected signals lack an astrophysical interpretation, they serve
as a compelling test for our method’s capability of recovering complex
spectral shapes. We inject separately two SGWBs, with SNR of 209 and 5,
respectively. We choose the following parameters for each signal:

  * •

Power-law baseline with $\log_{10}A=-43,\gamma=-2$ (SNR=209) and
$\log_{10}A=-44.5,\gamma=-1.5$ (SNR=5);

  * •

Knots amplitudes $\log_{10}\bm{\delta}^{h}=\left[0.6,-0.2,-0.3,0.4\right]$,
identical in both cases.

We choose the parameters as representative of two SGWBs, above and below the
nominal power-law sensitivity curve [59] to an SGWB with SNR of 10 after 4
years of observation, respectively. We emphasize once more that the injected
model is an EGP with local features arising with respect to a simple power-
law. In this specific context, the PLS is not a faithful indicator of a signal
detectability, but should be interpreted as a rough reference level.

Spectral inference results and parameter posteriors are shown in Fig. 3 and
Fig. 8 in App. B, respectively. The EGP model is able to capture the spectral
shape of the high SNR injected signal, and posterior distributions are
consistent with the injected values within their $90\%$ credible intervals.
The posterior distributions for the noise parameters are consistent with zero,
as expected. If the signal SNR is too low, as for the second simulated SGWB
(SNR=5) the posterior distribution for its parameters is effectively an upper
limit, only.

###  V.2 Astrophysical Case: EMRI SGWB

We further apply our method to the astrophysical case of an SGWB from EMRIs.
We inject a background signal chosen from a set available in literature [1]:
twelve background signals from astrophysically motivated populations [5] were
generated (therein referred to as M1 to M12). The stochastic signal is
constructed through suitable processing of whole population signals, and
subtraction of resolvable sources via an iterative algorithm [60]. Each source
GW strain is computed using the state-of-the-art EMRI waveforms [61] and the
resulting LISA signal is simulated through an accurate response model.
Resulting SGWBs exhibit a wide range of amplitudes, though the majority of
them provides an SGWB potentially detectable with LISA over 4 years of
observation time.

In this study, we use M1 as our fiducial model, as it provides an intermediate
SGWB amplitude across the populations studied in [1]. Because of uncertainties
on spectral morphology and population distribution, our flexible inference
method is suitable to analyse such signal. In order to inject coherently the
signal in the Balrog simulator, we divide the A channel signal realization
obtained by [1] by the response function, following the formalism in Eq. (2).
Then, we smooth the resulting spectrum it to obtain a reference model PSD.
Subsequently, we re-generate each TDI channel data and use it for inference.

####  V.2.1 Exploring Hyperparameters Space

We perform a set of Bayesian inferences, exploring two hyperparameters for the
signal model: the number of knots, $n=3,4,\ldots,8$, and the kernel length-
scale, as integer $m_{\sigma}$ multiples of $\sigma=0.6$, with
$m_{\sigma}=\\{1,2,4,8,16,32\\}$.

Figure 4: Top Panel: Grid of logarithmic Bayes factor for different pairs of
number of signal knots ($n$) and integers of Kernel length scale $\sigma$. The
numbers are computed in relation to the model that has the highest evidence,
specifically $\left(n,m_{\sigma}\right)=(8,32)$. We observe a diagonal
pattern: as $m_{\sigma}$ decreases and $n$ increases, the model becomes
progressively less preferred. Bottom Panel: Computational Cost (i.e. Sampling
time as reported by Nessai) of parameter estimations as function of integer
multiple of $\sigma$, $m_{\sigma}$ (with $\sigma=0.6$) and the number of knots
of the signal model, $n$. The numbers are in unit of kilohours. Each inference
is run in multi-threading over 40 cores.  Figure 5: Posterior distribution for
the inference on EMRIs background, with a noise-only model. The noise
estimation employs 3 knots, with frequencies corresponding to vertical dashed
gray lines in the main panel. Bottom sub-panels show posteriors on individual
knot amplitudes. The injected spectral noise model corresponds to
$\bm{\delta}^{n}=\bm{0}$, and is denoted by solid blue lines. The noise model
can compensate with significant biases the presence of a SGWB in the data.
Nonetheless, its evidence is significantly lower when compared to the noise
and signal models, as illustrated in Fig. 4.

We evaluate each evidence through stochastic nested sampling. Results are
presented in Fig. 4 (top panel), as Bayes factors $\log_{10}\mathcal{B}$ with
respect to the highest-evidence model, i.e.
$\left(n,m_{\sigma}\right)=(8,32)$. They reveal a prevailing trend: as
$m_{\sigma}$ increases, the number of knots become less effective at
influencing each model evidence, resulting in Bayes factors increasingly close
to each as function of $m_{\sigma}$.

For completeness,in Fig. 4 (bottom panel) we also provide the computational
cost of each inference, measured in CPU core kilohours (kh). The observed
trend is the result of two competing factors: as the number of parameters
(i.e. the number of knots $n$) increases, the inference process would take
longer. However, if the increased model complexity (e.g. shorter kernel
length-scale) is not required by the simulated data, exploration of parameter
space is typically fast. As a result, we observe a pattern in computational
cost broadly similar to the log Bayes factors: as the number of knots
($\sigma$) increase (decreases), the computational cost diminishes.

As a further check we perform inference for the null-hypothesis ${\cal
H}_{0}$, where we exclusively model the data as instrumental noise. Results
reveal decisive evidence, in favour of any signal-including model
$\left(n,m_{\sigma}\right)$ with respect to the noise-only model ${\cal
H}_{0}$, with log-Bayes factors largely greater than $10^{3}$. As expected,
the spectral noise reconstruction in Fig. 5 shows significant biases to
compensate for the unmodelled SGWB signal in the data.

####  V.2.2 Spectral reconstruction

We discuss here in greater detail the signal spectral reconstruction of each
model. Our findings are best described by a comparison between the inferences
with $(n,m_{\sigma})$ equal to $(8,1)$ and $(8,32)$, whose results are shown
in top and bottom panels of both Figs. 6 and 7, respectively. We show the
reconstruction of the fiducial EMRI background with and without T channel data
in Figs. 6 and 7, respectively. The T channel is known to be less sensitive to
GWs than the A and E channel at the frequencies of our analyses, therefore
acting as a potential noise estimator. We observe however that the spectral
reconstructions with and without the T channel are broadly consistent, with
broader posteriors when the T channel is excluded, as expected. This is due to
the common noise modelling of the single link PSD spectrum, which is then
propagated coherently to each channel through $M_{\rm TDI}$ Eq. (5). Relaxing
the exact knowledge of the transfer matrix is expected to significantly
broaden the noise PSD posterior.

A comparison between top and bottom panels reveals an important result,
originating from the competition in Bayes factors between excess model
complexity and the effectiveness of our chosen parameterization. The inference
presented in this study is inherently integrated across a large frequency
band: we annotate the injected SGWB with blue values for cumulative SNRs with
respect to the true noise level on a set of representative frequencies. It is
readily apparent that $80\%$ of the $\rm SNR$ is accumulated between $1$ and
$4{\rm mHz}$: in this frequency range the reconstruction is consistent with
the injected spectrum. The dashed red lines and red shaded area denote the
posterior median and $90\%$ confidence interval respectively.

At the lowest frequencies, the preferred model (top panel,
$n,m_{\sigma}=8,32$) exhibits larger biases as compared to the disfavoured one
(bottom panel, $n,m_{\sigma}=8,1$). This is due to the narrow-band
($m_{\sigma}=1$) fluctuations being highly disfavoured by the data in the high
SNR region, where the spectrum is close to a powerlaw model, a shape already
described by the baseline parameters. On the contrary broad-band
($m_{\sigma}=32$) fluctuations can be absorbed by an adjustment of the
baseline parameters, while keeping an accurate spectral reconstruction.
However, the preference for $m_{\sigma}=32$ has an important side effect:
$m_{\sigma}$ is fixed and equal for every knot in each analysis, including the
four at the lowest frequencies. Therefore, the model with the highest evidence
is not capable of introducing local fluctuations to capture significant tilts
with respect to the baseline powerlaw. These are instead provided by the
disfavoured model.

Figs. 9 and 10 in Appendix further support our interpretation, where we
compare the posterior distributions for the two runs shown in Fig 6 (red
contours) with two more differing only by a broader prior for the knot
amplitudes $\bm{\delta}^{h}$, uniform in $\left[-5,5\right]$ (blue contours).
The inference results for the noise are largely unaffected and uncorrelated to
the signal parameters (and hyperparameters). However, we find that relaxing
the knot priors affects the posterior on the baseline parameters: a sign of
strong correlation. Moreover, posteriors on $\bm{\delta}^{h}$ for the
preferred runs (Fig. 9) are prior dominated, while the likelihood becomes
informative in the disfavoured models $(n,m_{\sigma})=(8,1)$ (Fig. 10). In
conclusion we argue that in a realistic scenario models more flexible than EGP
should be developed, alongside astrophysically motivated priors to avoid
introducing –even in an evidence based model selection– unwantend biases.

Figure 6: Signal reconstruction of the M1 EMRI background carried out using
the EGP parametrization. In both figures, the noise estimation employs 3
knots. Top Panel: the model features 8 knots for the signal and a kernel
length scale that is 32 times larger than $\sigma$. This configuration aligns
with the model exhibiting the highest evidence. Bottom Panel: the model
features 8 knots for the signal, and the kernel length scale is set equal to
$\sigma$. This configuration corresponds to the model with the lowest
evidence. Blue numbers denote the cumulative SNR integrated from
$10^{-4}\rm{Hz}$ up to the frequency they are located at.

####  V.2.3 T Channel Influence

We aim to investigate the impact of the T channel on parameter estimation.
Given its reduced sensitivity to the SGWB, we anticipate that excluding the T
channel from data analysis will result in less informative posterior
distributions (as depicted in Fig. 7). To be more specific, we seek to assess
whether hyperparameterization modeling remains influential, as depicted in
Fig. 4.

To achieve this, we examine two models characterized by the highest and lowest
values of evidence (or hyperlikelihood), denoted as
$\left(n,m_{\sigma}\right)=\left(8,32\right)$ and
$\left(n,m_{\sigma}\right)=\left(8,1\right)$, respectively. The calculated
logarithmic Bayes factor stands at $-10.95$, underscoring a substantial
strength of evidence in favor of one model over the other. In other words, the
outcome suggests that the trend observed in Fig. 4 is robust and not
significantly affected by variations in our ability to estimate the noise
level.

Figure 7: Analogue illustration of Fig. 6 where the inference has been
conducted without incorporating the T channel in the data. Consequently, the
resulting posterior regions exhibit greater dispersion, signifying reduced
informativeness. Blue numbers denote the cumulative SGWB SNR integrated from
$10^{-4}\rm{Hz}$ up to the frequency they are located at.

##  VI Conclusion

In this study we present a novel method to estimate parameters of an SGWB
observed by LISA. We implement into our analysis in Balrog, a broader data
analysis framework for LISA. This is a building block towards the completion
of a full global fit of deterministic and stochastic GW signals, alongside a
flexible noise model. We employ a versatile parameterization to infer the SGWB
of uncertain spectral shape using EGPs. Our method is able to capture
unexpected features in the signal shape, while assessing their significativity
in a purely Bayesian formalism. We test our algorithm using a simplified toy-
model, where the true injected model is by construction within the parameter
space explored.

Following, we apply the analysis to the expected signal from a population of
unresolved EMRIs. We leverage the model flexibility to (i)explore the
hyperparameter space and optimize them based on evidence. We also investigate
the impact of T channel on the inference results. We find that the EGP model
has great flexibility at reconstructing the SGWB spectral shape. However,
while Bayesian evidence might be suitable to select models in the high SNR
region, it may disfavour less biased models in the low SNR region, at the
price of excess complexity where it is not needed. Therefore, the model
inherent non-locality should be treated with caution when interpreting
spectral results. We also find that the T channel has a significant impact on
the inference results, and its exclusion leads to naturally less constrained
results.

Our findings suggest that parametric or weakly-parametric models (as the EGP)
should be considered a first step towards even more agnostic and realistic
approaches. Including cross-correlation terms between TDI channels is expected
to provide a more robust inference, and it is a natural extension of our work.
Finally, for specific SGWBs a parameter estimation modelling non-
stationarities and non-Gaussianities would extract valuable information from
the data, thus further reducing biases.

###### Acknowledgements.

We wish to thank Jonathan Gair, Antoine Klein, Golam Shaifullah, Germano
Nardini, Chiara Caprini, Mauro Pieroni, Nikos Karnesis, Fabio Rigamonti and
all Balrog developers for fruitful discussions. We also wish to thank the
University of Geneva for the kind hospitality in summer 2023. R.B.
acknowledges the support of the ICSC National Research Center funded by
NextGenerationEU, and the Italian Space Agency grant _Phase A activity for
LISA mission_ , Agreement n.2017-29-H.0. C.J.M. acknowledges the support of
the UK Space Agency grant, no. ST/V002813/1. A.S. acknowledges the financial
support provided under the European Union’s H2020 ERC Consolidator Grant
“Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691).
Computational resources were provided by University of Birmingham BlueBEAR
High Performance Computing facility, by CINECA through EuroHPC Benchmark
access call grant “LISAFIT” (proposal ID: EHPC-BEN-2023B08-024), and by the
Google Cloud Research Credits program with the award GCP19980904. Software: We
acknowledge usage of Mathematica [62] and of the following Python [63]
packages for modeling, analysis, post-processing, and production of results
throughout: Nessai [54], matplotlib [64], numpy [65], scipy [66].

##  Appendix A SGWB response function

For completeness and consistency of definitions, we provide here a derivation
of the one-arm LISA response –induced by a GW propagating along $\hat{k}$– as
experienced along a single link identified by the sender and receiver
spacecraft $s$ and $r$ and the unit-vector pointing from the former to the
latter, i.e. $\hat{r}_{sr}$. Henceforth, no summation over repeated $s$ and
$r$ is assumed. Each fractional armlength variation reads

$W_{sr}(t,\hat{k})=\frac{\delta L_{sr}(t)}{L_{sr}}\\\
=\frac{1}{2}\hat{r}^{a}_{sr}\hat{r}^{b}_{sr}\int^{\infty}_{-\infty}\tilde{h}_{ab}(f)\tau_{sr}(f,t,\hat{k})e^{2\pi
if(t-\Delta t)}df,$ (22)

where

$\Delta
t=\frac{1}{2}\left[\hat{k}\cdot(\vec{x}_{s}+\vec{x}_{r})+L_{sr}\right]$ (23)

and

$\tau_{sr}(f,t,\hat{k})={\rm sinc}(\pi
L_{sr}f[1-\hat{k}\cdot\hat{r}_{sr}(t)]).$ (24)

where $\vec{x}_{s}$ and $\vec{x}_{r}$ denote the sender and receiver
spacecraft, respectively. Due to the motion of the satellites, the arm
direction vector $\hat{r}$ varies with time. The GW strain can be written in
TT gauge as:

$h_{ab}(t,\hat{k})=h_{+}(t)\epsilon^{+}_{ab}(\hat{k})+h_{\times}(t)\epsilon^{\times}_{ab}(\hat{k}).$
(25)

where $\epsilon^{+,\times}$ denote a linear polarization-tensor basis. The
stochastic strain is the linear superposition of plane tensor-waves from each
direction of the sky, therefore we define an integrated response to each
polarization $p=+,\times$, as follows

$W_{sr}^{p}(t)=\int_{\hat{k}}W_{sr}^{p}(t,\hat{k})d^{2}\hat{k},$ (26)

where

$\begin{split}W_{sr}^{p}(t,\hat{k})&=\frac{1}{2}\hat{r}^{a}_{sr}\hat{r}^{b}_{sr}\int^{\infty}_{-\infty}\tilde{h}_{p}(f,\hat{k})\epsilon_{ab}^{p}(\hat{k})\tau(f,t,\hat{k})e^{2\pi
if(t-\bm{\delta}t)}df=\\\
&=\frac{1}{2}\int^{\infty}_{-\infty}\tilde{h}_{p}(f,\hat{k})\xi^{p}(\hat{k},\hat{r}_{sr})\tau(f,t,\hat{k})e^{2\pi
if(t-\bm{\delta}t)}df.\end{split}$ (27)

In Eq. (26) we introduced the antenna pattern functions
$\xi^{p}(\hat{k},\hat{r}_{sr})$. They read as:

$\begin{split}&\xi^{+}(\hat{k},\hat{r}_{sr})=(\hat{u}_{k}\cdot\hat{r}_{sr})^{2}-(\hat{v}_{k}\cdot\hat{r}_{sr})^{2}\\\
&\xi^{\times}(\hat{k},\hat{r}_{sr})=2(\hat{u}_{k}\cdot\hat{r}_{sr})(\hat{v}_{k}\cdot\hat{r}_{sr}).\end{split}$
(28)

We rewrite Eq. (27) equivalently as:

$W_{sr}^{p}(t,\hat{k})=\int^{\infty}_{-\infty}\tilde{h}_{p}(f,\hat{k})\overline{G}^{p}_{sr}(f,t,\hat{k})e^{2\pi
ift}df,$ (29)

where $\overline{G}^{p}_{sr}(f,t,\hat{k})$ reads

$\overline{G}^{p}_{sr}(f,t,\hat{k})=\frac{1}{2}\xi^{p}(\hat{k},\hat{r}_{sr})e^{-i\pi
f(L_{sr}+\hat{k}\cdot(x_{s}+x_{r}))}\tau_{sr}(f,t,\hat{k}).$ (30)

Noteworthy, $\overline{G}^{p}_{sr}(f,t,\hat{k})$ is a complex quantity unlike
$\tau(f,t,\hat{k})$. Equivalently, we recast Eq. (30) as a response in
fractional frequency domain

$G^{p}_{sr}(f,t,\hat{k})=2\pi iL_{sr}f\overline{G}^{p}_{sr}(f,t,\hat{k}).$
(31)

The Fourier transform of (29) is then

$\displaystyle W_{sr}^{p}(f,\hat{k})$
$\displaystyle=\int^{\infty}_{-\infty}dte^{-2\pi ift}W_{sr}^{p}(t,\hat{k})$
(32) $\displaystyle=\int^{\infty}_{-\infty}dte^{-2\pi
ift}\int^{\infty}_{-\infty}\tilde{h}_{p}(f^{\prime},\hat{k})\overline{G}^{p}_{sr}(f^{\prime},t,\hat{k})e^{2\pi
if^{\prime}t}df^{\prime}$
$\displaystyle=\int^{\infty}_{-\infty}dt\int^{\infty}_{-\infty}\tilde{h}_{p}(f^{\prime},\hat{k})\overline{G}^{p}_{sr}(f^{\prime},t,\hat{k})e^{-2\pi
it(f-f^{\prime})}df^{\prime}$
$\displaystyle=\int^{\infty}_{-\infty}df^{\prime}\overline{G}^{p}_{sr}(f^{\prime},f-f^{\prime},\hat{k})\tilde{h}_{p}(f^{\prime},\hat{k}),$

note that

$\overline{G}^{p}_{sr}(f^{\prime},f,\hat{k})=\int^{\infty}_{-\infty}e^{-2\pi
ift}\overline{G}^{p}_{sr}(f^{\prime},t,\hat{k})dt.$ (33)

As described in Eq.(3) an isotropic, stationary, zero-mean background is
characterized entirely by the strain covariance which allows for the
definition of a response matrix for two given links $sr$ and
$s^{\prime}r^{\prime}$ reads

$\Sigma_{sr,s^{\prime}r^{\prime}}(f)=\left<W^{p}_{sr}(f)W^{*p}_{s^{\prime}r^{\prime}}(f)\right>=\int_{\hat{k}}d^{2}k\int^{\infty}_{-\infty}df^{\prime}\overline{G}^{p}_{sr}(f^{\prime},f-f^{\prime},\hat{k})\overline{G}^{*p}_{s^{\prime}r^{\prime}}(f^{\prime},f-f^{\prime},\hat{k})S_{h}(f^{\prime})$
(34)

For an isotropic distribution of sources, integrating the product of
$\overline{G}$’s and evaluating it at a reference time $t_{0}$ –without loss
of generality given the assumption of perfect stationarity– we obtain Thus,

$\Sigma^{p}_{sr,s^{\prime}r^{\prime}}(f)=S_{h}(f)G^{p}_{sr,s^{\prime}r^{\prime}}(f),$
(35)

with

$G^{p}_{sr,s^{\prime}r^{\prime}}(f)=\int_{\hat{k}}d^{2}k\overline{G}^{p}_{sr}(f,t_{0},\hat{k})\overline{G}^{*p}_{s^{\prime}r^{\prime}}(f,t_{0},\hat{k})$
(36)

which corresponds to $G_{ij}$ in Eq. (4) upon remapping of the indices
$sr,s^{\prime}r^{\prime}$ to $ij$.

##  Appendix B Posterior distribution for the toy model

Figure 8: Posterior parameter distribution for the toy-model SGWB with SNR=209
described in Sec. V.1, obtained with an EGP model inference. The signal
injected parameters are
$\left\\{\log_{10}\bm{\delta}^{h},\log_{10}A,\gamma\right\\}=\left\\{\left[0.6,-0.2,-0.3,0.4\right],-43,-2\right\\}$.
The injected noise level matches the nominal LISA sensitivity curve [52].
Darker (lighter) shaded areas denote 90% (50%) credible regions, and solid
orange lines indicate the injected values. Posteriors are consistent with the
injected values at the 90% credible level, and the noise parameters are
consistent with zero with accuracies between $10^{-3}$ and $10^{-2}$.

##  Appendix C Posterior distributions for Model M1 SGWB

Figure 9: Posterior parameter distribution for M1 EMRI SGWB for model
$(n,m_{\sigma})=(8,32)$ using a larger and narrower prior on
$\bm{\delta}^{h}$. Posteriors shown in this plot correspond to the spectral
reconstructions in top panels of Fig. 6 and 7,respectively. They support the
result interpretation provided in Sec. V.2.2.

Figure 10: Posterior parameter distribution for M1 EMRI SGWB for model
$(n,m_{\sigma})=(8,1)$ using a larger and narrower prior on $\bm{\delta}^{h}$.
Posteriors shown in this plot correspond to the spectral reconstructions in
bottom panels of Fig. 6 and 7, respectively. They support the result
interpretation provided in Sec. V.2.2.

## References

  * Pozzoli  _et al._ [2023] F. Pozzoli, S. Babak, A. Sesana, M. Bonetti, and N. Karnesis, arXiv e-prints (2023), arXiv:2302.07043 [astro-ph.GA] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, and et al., arXiv e-prints (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Klein  _et al._ [2016] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, I. Hinder, F. Ohme, and B. Wardell, Phys. Rev. D 93, 024003 (2016), arXiv:1511.05581 [gr-qc] .
  * Sesana  _et al._ [2004] A. Sesana, F. Haardt, P. Madau, and M. Volonteri, ApJ 611, 623 (2004), arXiv:astro-ph/0401543 [astro-ph] .
  * Babak  _et al._ [2017] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Phys. Rev. D 95, 103012 (2017), arXiv:1703.09722 [gr-qc] .
  * Korol  _et al._ [2020] V. Korol, S. Toonen, A. Klein, V. Belokurov, F. Vincenzo, R. Buscicchio, D. Gerosa, C. J. Moore, E. Roebber, E. M. Rossi, and A. Vecchio, A&A 638, A153 (2020), arXiv:2002.10462 [astro-ph.GA] .
  * Antoniadis  _et al._ [2023] J. Antoniadis, P. Arumugam, S. Arumugam, P. Auclair, and et al., arXiv e-prints (2023), arXiv:2306.16227 [astro-ph.CO] .
  * Afzal  _et al._ [2023] A. Afzal, G. Agazie, A. Anumarlapudi, A. M. Archibald, and et al., ApJ 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * EPTA Collaboration  _et al._ [2023] EPTA Collaboration, InPTA Collaboration, J. Antoniadis, P. Arumugam, and et al., A&A 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE] .
  * Agazie  _et al._ [2023] G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, and et al., ApJ 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Reardon  _et al._ [2023] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, and et al., ApJ 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Xu  _et al._ [2023] H. Xu, S. Chen, Y. Guo, J. Jiang, and et al., Research in Astronomy and Astrophysics 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Caprini [2015] C. Caprini, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 610 (2015) p. 012004, arXiv:1501.01174 [gr-qc] .
  * Ricciardone [2017] A. Ricciardone, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 840 (2017) p. 012030, arXiv:1612.06799 [astro-ph.CO] .
  * Gowling and Hindmarsh [2021] C. Gowling and M. Hindmarsh, J. Cosmology Astropart. Phys. 2021, 039 (2021), arXiv:2106.05984 [astro-ph.CO] .
  * Boileau  _et al._ [2023] G. Boileau, N. Christensen, C. Gowling, M. Hindmarsh, and R. Meyer, J. Cosmology Astropart. Phys. 2023, 056 (2023), arXiv:2209.13277 [gr-qc] .
  * Auclair  _et al._ [2020] P. Auclair, J. J. Blanco-Pillado, D. G. Figueroa, A. C. Jenkins, and et al., J. Cosmology Astropart. Phys. 2020, 034 (2020), arXiv:1909.00819 [astro-ph.CO] .
  * Clesse and García-Bellido [2017] S. Clesse and J. García-Bellido, Physics of the Dark Universe 18, 105 (2017), arXiv:1610.08479 [astro-ph.CO] .
  * Nelemans [2009] G. Nelemans, Classical and Quantum Gravity 26, 094030 (2009), arXiv:0901.1778 [astro-ph.SR] .
  * Babak  _et al._ [2023] S. Babak, C. Caprini, D. G. Figueroa, N. Karnesis, P. Marcoccia, G. Nardini, M. Pieroni, A. Ricciardone, A. Sesana, and J. Torrado, J. Cosmology Astropart. Phys. 2023, 034 (2023), arXiv:2304.06368 [astro-ph.CO] .
  * Bonetti and Sesana [2020] M. Bonetti and A. Sesana, Phys. Rev. D 102, 103023 (2020), arXiv:2007.14403 [astro-ph.GA] .
  * Hartwig  _et al._ [2023] O. Hartwig, M. Lilley, M. Muratore, and M. Pieroni, Phys. Rev. D 107, 123531 (2023), arXiv:2303.15929 [gr-qc] .
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of Markov Chain Monte Carlo_ (CRC press, 2011).
  * Skilling [2004] J. Skilling, in  _Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering_, American Institute of Physics Conference Series, Vol. 735, edited by R. Fischer, R. Preuss, and U. V. Toussaint (2004) pp. 395–405.
  * Adams and Cornish [2014] M. R. Adams and N. J. Cornish, Phys. Rev. D 89, 022001 (2014), arXiv:1307.4116 [gr-qc] .
  * Baghi  _et al._ [2023] Q. Baghi, N. Karnesis, J.-B. Bayle, M. Besançon, and H. Inchauspé, J. Cosmology Astropart. Phys. 2023, 066 (2023), arXiv:2302.12573 [gr-qc] .
  * Caprini  _et al._ [2019] C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone, and G. Tasinato, J. Cosmology Astropart. Phys. 2019, 017 (2019), arXiv:1906.09244 [astro-ph.CO] .
  * Flauger  _et al._ [2021] R. Flauger, N. Karnesis, G. Nardini, M. Pieroni, A. Ricciardone, and J. Torrado, J. Cosmology Astropart. Phys. 2021, 059 (2021), arXiv:2009.11845 [astro-ph.CO] .
  * Rosca-Mead  _et al._ [2023] R. Rosca-Mead, M. Agathos, C. J. Moore, and U. Sperhake, Phys. Rev. D 107, 124040 (2023), arXiv:2302.04995 [gr-qc] .
  * Rosado [2012] P. A. Rosado, Phys. Rev. D 86, 104007 (2012), arXiv:1206.1330 [gr-qc] .
  * Boileau  _et al._ [2021] G. Boileau, A. Lamberts, N. Christensen, N. J. Cornish, and R. Meyer, arXiv e-prints (2021), arXiv:2105.06793 [gr-qc] .
  * Timpano  _et al._ [2006] S. E. Timpano, L. J. Rubbo, and N. J. Cornish, Phys. Rev. D 73, 122001 (2006), arXiv:gr-qc/0504071 [gr-qc] .
  * Stroeer and Vecchio [2006] A. Stroeer and A. Vecchio, Classical and Quantum Gravity 23, S809 (2006), arXiv:astro-ph/0605227 [astro-ph] .
  * Finch  _et al._ [2023] E. Finch, G. Bartolucci, D. Chucherko, B. G. Patterson, V. Korol, A. Klein, D. Bandopadhyay, H. Middleton, C. J. Moore, and A. Vecchio, MNRAS 522, 5358 (2023), arXiv:2210.10812 [astro-ph.SR] .
  * Kupfer  _et al._ [2023] T. Kupfer, V. Korol, T. B. Littenberg, S. Shah, and et al., arXiv e-prints (2023), arXiv:2302.12719 [astro-ph.SR] .
  * Gaia Collaboration  _et al._ [2016] Gaia Collaboration, T. Prusti, J. H. J. de Bruijne, A. G. A. Brown, and et al., A&A 595, A1 (2016), arXiv:1609.04153 [astro-ph.IM] .
  * LSST Science Collaboration  _et al._ [2009] LSST Science Collaboration, P. A. Abell, J. Allison, S. F. Anderson, and et al., arXiv e-prints (2009), arXiv:0912.0201 [astro-ph.IM] .
  * Georgousi  _et al._ [2023] M. Georgousi, N. Karnesis, V. Korol, M. Pieroni, and N. Stergioulas, MNRAS 519, 2552 (2023), arXiv:2204.07349 [astro-ph.GA] .
  * Breivik  _et al._ [2020] K. Breivik, C. M. F. Mingarelli, and S. L. Larson, ApJ 901, 4 (2020), arXiv:1912.02200 [astro-ph.GA] .
  * Sesana  _et al._ [2020] A. Sesana, A. Lamberts, and A. Petiteau, MNRAS 494, L75 (2020), arXiv:1912.07627 [astro-ph.GA] .
  * Sesana [2016] A. Sesana, Phys. Rev. Lett. 116, 231102 (2016), arXiv:1602.06951 [gr-qc] .
  * Klein  _et al._ [2022] A. Klein, G. Pratten, R. Buscicchio, P. Schmidt, C. J. Moore, E. Finch, A. Bonino, L. M. Thomas, N. Williams, D. Gerosa, S. McGee, M. Nicholl, and A. Vecchio, arXiv e-prints (2022), arXiv:2204.03423 [astro-ph.HE] .
  * Buscicchio  _et al._ [2021] R. Buscicchio, A. Klein, E. Roebber, C. J. Moore, D. Gerosa, E. Finch, and A. Vecchio, Phys. Rev. D 104, 044065 (2021), arXiv:2106.05259 [astro-ph.HE] .
  * Toubiana  _et al._ [2022] A. Toubiana, S. Babak, S. Marsat, and S. Ossokine, Phys. Rev. D 106, 104034 (2022), arXiv:2206.12439 [gr-qc] .
  * Digman and Cornish [2023] M. C. Digman and N. J. Cornish, Phys. Rev. D 108, 023022 (2023), arXiv:2212.04600 [gr-qc] .
  * Phinney [2001] E. S. Phinney, arXiv e-prints (2001), arXiv:astro-ph/0108028 [astro-ph] .
  * Li  _et al._ [2023] R.-N. Li, Z.-Y. Zhao, Z. Gao, and F.-Y. Wang, ApJ 956, L2 (2023), arXiv:2309.09448 [astro-ph.HE] .
  * Toscani  _et al._ [2023] M. Toscani, O. Burke, C. Liu, N. B. Zamel, N. Tamanini, and F. Pozzoli, arXiv e-prints , arXiv:2307.06722 (2023), arXiv:2307.06722 [astro-ph.CO] .
  * Tinto and Dhurandhar [2005] M. Tinto and S. V. Dhurandhar, Living Reviews in Relativity 8, 4 (2005).
  * Tinto  _et al._ [2004] M. Tinto, F. B. Estabrook, and J. W. Armstrong, Phys. Rev. D 69, 082001 (2004), arXiv:gr-qc/0310017 [gr-qc] .
  * Moore  _et al._ [2015] C. J. Moore, R. H. Cole, and C. P. L. Berry, Classical and Quantum Gravity 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * LISA Science Study Team [2018] LISA Science Study Team,  _LISA Science Requirements Document_, Tech. Rep. 1.0 (ESA, 2018).
  * Smith and Caldwell [2019] T. L. Smith and R. R. Caldwell, Phys. Rev. D 100, 104055 (2019), arXiv:1908.00546 [astro-ph.CO] .
  * Williams  _et al._ [2021] M. J. Williams, J. Veitch, and C. Messenger, Phys. Rev. D 103, 103006 (2021), arXiv:2102.11056 [gr-qc] .
  * Cutler and Flanagan [1994] C. Cutler and É. E. Flanagan, Phys. Rev. D 49, 2658 (1994), arXiv:gr-qc/9402014 [gr-qc] .
  * Green [1995] P. J. Green, Biometrika 82, 711 (1995).
  * Karnesis  _et al._ [2023] N. Karnesis, M. L. Katz, N. Korsakova, J. R. Gair, and N. Stergioulas, MNRAS 526, 4814 (2023), arXiv:2303.02164 [astro-ph.IM] .
  * Cornish and Littenberg [2015] N. J. Cornish and T. B. Littenberg, Classical and Quantum Gravity 32, 135012 (2015), arXiv:1410.3835 [gr-qc] .
  * Thrane and Romano [2013] E. Thrane and J. D. Romano, Phys. Rev. D 88, 124032 (2013), arXiv:1310.5300 [astro-ph.IM] .
  * Karnesis  _et al._ [2021] N. Karnesis, S. Babak, M. Pieroni, N. Cornish, and T. Littenberg, Phys. Rev. D 104, 043019 (2021), arXiv:2103.14598 [astro-ph.IM] .
  * Katz  _et al._ [2021] M. L. Katz, A. J. K. Chua, L. Speri, N. Warburton, and S. A. Hughes, Phys. Rev. D 104, 064047 (2021), arXiv:2104.04582 [gr-qc] .
  * Wolfram Research Inc. [2022] Wolfram Research Inc., Mathematica (2022).
  * Van Rossum and Drake [2009] G. Van Rossum and F. L. Drake,  _Python 3 Reference Manual_ (CreateSpace, Scotts Valley, CA, 2009).
  * Hunter [2007] J. D. Hunter, Computing in Science and Engineering 9, 90 (2007).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, and et al., Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, and et al., Nature Methods 17, 261 (2020), arXiv:1907.10121 [cs.MS] .


# PROFICIENCY ASSESSMENT OF L2 SPOKEN ENGLISH USING WAV2VEC 2.0

###### Abstract

The increasing demand for learning English as a second language has led to a
growing interest in methods for automatically assessing spoken language
proficiency. Most approaches use hand-crafted features, but their efficacy
relies on their particular underlying assumptions and they risk discarding
potentially salient information about proficiency. Other approaches rely on
transcriptions produced by ASR systems which may not provide a faithful
rendition of a learner’s utterance in specific scenarios (e.g., non-native
children’s spontaneous speech). Furthermore, transcriptions do not yield any
information about relevant aspects such as intonation, rhythm or prosody. In
this paper, we investigate the use of wav2vec 2.0 for assessing overall and
individual aspects of proficiency on two small datasets, one of which is
publicly available. We find that this approach significantly outperforms the
BERT-based baseline system trained on ASR and manual transcriptions used for
comparison.

Index Terms—  automatic assessment of spoken language proficiency, computer
assisted language learning

##  1 Introduction

With the increasing number of learners of English as a second language (L2)
worldwide, there has been a growing demand for automated spoken language
assessment systems both for formal settings, such as internationally
recognised language tests, and for practice situations in Computer Assisted
Language Learning (CALL). In fact, compared to human graders, automatic
graders can ensure greater consistency and speed at a lower cost, since the
recruitment and training of new human experts are expensive and can offer only
a small increase in performance [1].

In automatic assessment of L2 spoken language proficiency, input sequential
data from a learner is used to predict a grade or a level with respect to
holistic proficiency or specific aspects of proficiency [2, 3]. The input may
consist of acoustic features, recognised words, phones and/or time-aligned
information, or other information, such as fundamental frequency, extracted
directly from the audio or from automatic speech recognition (ASR)
transcriptions. Most approaches in the literature extract sets of hand-crafted
features related to specific aspects of proficiency, such as fluency [4],
pronunciation [5], prosody [6] and text complexity [7], which are then fed
into graders to predict analytic scores targeting those specific aspects. A
similar approach can be used by concatenating multiple hand-crafted features
targeting more than one aspect in order to produce overall feature sets, which
are then passed through graders to predict holistic grades, as shown in [8, 9,
10, 11]. The effectiveness of hand-crafted features for grading either single
aspects or overall proficiency heavily relies on their particular underlying
assumptions, and they risk discarding potentially salient information about
proficiency. This issue for holistic grading has been addressed by replacing
hand-crafted features with automatically derived features for holistic grading
prediction, either through an end-to-end system [12] or in multiple stages
[13, 14]. Other studies have employed graders that are trained on holistic
grades but are defined with both their inputs and topology adapted to focus on
specific aspects of proficiency, such as pronunciation [15], rhythm [16] and
text [17, 18]. In these cases, an evident limitation is the lack of
information concerning aspects of proficiency that are not included in the
input data fed to the grader, although it is possible to combine multiple
graders targeting different aspects, as shown in [19]. This is particularly
true for systems using ASR transcriptions for two reasons: first, they contain
a certain word error rate (WER) and may not faithfully render the contents of
a learner’s performance; secondly, although transcriptions might preserve some
information about pronunciation (e.g., in the ASR confidence scores), they do
not yield any information about other relevant aspects of a learner’s
performance, e.g., intonation, rhythm or prosody. Instead, they remain a
valuable resource for highly specific tasks in CALL applications, such as
spoken grammatical error correction and feedback [20].

In this work, to address these issues and limitations, we propose an approach
based on self-supervised speech representations using wav2vec 2.0 [21, 22].
Recent studies have demonstrated that self-supervised learning (SSL) is an
effective approach in various downstream tasks of speech processing
applications, such as ASR, emotion recognition, keyword spotting, speaker
diarisation and speaker identification [21, 23]. In these studies, contextual
representations were applied by means of pre-trained models. Specifically, it
has been demonstrated that such models are capable of capturing a wide range
of speech-related features and linguistic information, such as audio, fluency,
suprasegmental pronunciation, and even syntactic and semantic text-based
features for L1, L2, read and spontaneous speech [24]. In the field of CALL,
SSL has been applied to mispronunciation detection and diagnosis [25, 26, 27]
and automatic pronunciation assessment [28], but, to the best of our
knowledge, it has not been investigated for the assessment of overall spoken
proficiency nor of other specific aspects of proficiency, such as formal
correctness, communicative effectiveness, lexical richness and complexity, and
relevance.

In this work, we first test the effectiveness of wav2vec 2.0 in predicting the
holistic proficiency level of L2 English learners’ answers included in a
publicly available dataset. Subsequently, we do the same on a non-publicly
available learner corpus, which also contains annotations related to single
aspects of proficiency (i.e., pronunciation, fluency, formal correctness,
communicative effectiveness, relevance and lexical richness and complexity)
that we try to predict with specific graders. The baseline system used for
comparison is a BERT [29] model that takes transcriptions as input. We use
only manual transcriptions for our experiments on the publicly available
dataset, whereas we use both manual and ASR transcriptions for our experiments
on the second dataset. In particular, the manual transcriptions also contain
hesitations and fragments of words, which serve as proxies for pronunciation
and fluency.

Another aspect that should be highlighted is that we conduct our experiments
using a small quantity of training data and still manage to achieve promising
results on both datasets.

Section 2 describes the data used in our experiments. Section 3 illustrates
the model architectures. Section 4 shows the results of our experiments.
Finally, in Section 5, we analyse and discuss the results and reflect upon
next steps.

##  2 Data

###  2.1 ICNALE

In order to test our approach, we consider the International Corpus Network of
Asian Learners of English (ICNALE) [30], a publicly available dataset
111http://language.sakura.ne.jp/icnale/download.html comprising written and
spoken responses of English learners ranging from A2 to B2 of the Common
European Framework of Reference (CEFR) for languages [31] and partially of
native speakers. The L1s of the speakers are not indicated, but they may be
inferred from their countries of origin: China, Hong Kong, Indonesia, Japan,
South Korea, Pakistan, Philippines, Singapore, Thailand, and Taiwan. The CEFR
levels were assigned prior to collecting the data, as the ICNALE team required
all the learners to take an L2 vocabulary size test and to present their
scores in English proficiency tests such as TOEFL, TOEIC, IELTS, etc. On the
basis of these two scores, the learners were classified into proficiency
levels. The spoken section of the corpus consists of two parts: one containing
monologues and the other containing dialogues. For our experiments, we only
considered the monologues, i.e., 4332 answers lasting between 36 and 69
seconds in which learners are required to express their opinion about two
statements on the importance of having a part-time job and on smoking in
restaurants. The available metadata includes manual transcriptions of the
learners’ answers, personal information about learners’ education history, and
their assigned CEFR levels. We divide the data into a training set of 3898
answers, a development set and a test set with 217 answers each. For the
experiments on this dataset, proficiency assessment is treated as a
classification task with five classes: A2, B1_1, B1_2, B2, and native speakers
(see Table 1). To the best of our knowledge, the ICNALE spoken monologues have
only been used in [32], but in this study, the answers to the two statements
are considered and evaluated independently, so no comparison is possible. The
experiments described in [33], instead, only include a section of essays and
spoken dialogues.

| Train | Dev | Test | Total
---|---|---|---|---
A2 | 299 | 16 | 17 | 332
B1_1 | 792 | 44 | 44 | 880
B1_2 | 1681 | 94 | 93 | 1868
B2 | 586 | 33 | 33 | 652
native | 540 | 30 | 30 | 600
Total | 3898 | 217 | 217 | 4332
Table 1: Number of answers for each CEFR proficiency level in ICNALE.

###  2.2 TLT-school

In Trentino, an autonomous region in northern Italy, the linguistic competence
of Italian students has been assessed in recent years through proficiency
tests in both English and German [34], involving about 3000 students ranging
from 9 to 16 years old, belonging to four different school grades ($5^{th}$,
$8^{th}$, $10^{th}$, $11^{th}$) and three CEFR proficiency levels (A1, A2,
B1). Since our experiments are conducted only on the B1 section of the English
spoken parts of the corpus, we will not describe the section concerning the
German section, as their analysis goes beyond the scope of this paper.

After eliminating the answers containing only silence or non-speech
background, the spoken section is composed of 494 responses to 7 small talk
questions about everyday life situations. It is worth mentioning that some
answers are characterized by a number of issues (e.g., presence of words
belonging to multiple languages or presence of off-topic answers). We decided
not to eliminate these answers from the data used in the experiments. As
regards the manual transcriptions, we did not eliminate the annotations
related to spontaneous speech phenomena, such as hesitations and fragments of
words, in order not to lose any possibly existing information about
pronunciation and fluency, although we acknowledge that they cannot replace
these actual speech phenomena completely. Detailed information about the
manual transcriptions and other aspects of the corpus can be found in [34].

As for the automatic speech recognition (ASR) output text, its WER is 35.9% on
the whole spoken test data, whereas it amounts to 41.13% for the B1 subset we
used in our experiments; acoustic and language models are described in [35].

The total score ranges from 0 to 12 in the spoken section and consists of the
sum of the analytic subscores assigned by human experts for each specific
proficiency indicator assigned by the human raters (i.e., relevance, formal
correctness, lexical richness and complexity, pronunciation, fluency, and
communicative effectiveness). For each indicator, human raters could choose 0,
1 or 2 points.

For this dataset, we treated proficiency assessment as a regression task when
predicting both the holistic score from 0 to 12 and the analytic subscores
ranging from 0 to 2.

In order to verify that these subscores are, in fact, targeting different
aspects of proficiency, we investigated the relationships between them with a
repeated measures design. Since our data violate both the sphericity and
normality assumptions required for rANOVA, we performed the Friedman test
[36], which is the non-parametric equivalent of rANOVA and determines whether
there are any statistically significant differences in ranks between the
distributions of multiple paired groups. As we obtained a significant _p_
-value, we found that there are significant differences among the subscores.
In order to determine exactly which pairs of subscores are significantly
different, we performed post-hoc multiple comparisons using the Nemenyi test
[37] (see Figure 1). All paired comparisons show significant differences (_p_
-value<0.05), except the pairs formal correctness-lexical richness and
complexity, formal correctness-pronunciation, pronunciation-communicative
effectiveness, and fluency-communicative effectiveness. The absence of
significant differences between the subscores related to correctness and those
related to lexical richness and complexity seem to be consistent with the fact
that: a) a poorer and simpler lexis tends to have fewer errors and b) the
human expert may have actually overlapped the two indicators due to the
presence of lexical errors, which can be connected either to poor formal
correctness or bad use of vocabulary. Similarly, pronunciation and correctness
do not show significant differences, because pronunciation errors might have
been incorporated in one or the other indicator. On the other hand, the
indicator related to communicative effectiveness intersects almost by
definition with those related to fluency and pronunciation. This has been
particularly true in recent years, as pronunciation tends to be assessed in
terms of general goals such as intelligibility [38] and communicative
effectiveness [39] rather than closeness to native English, especially in
English as an International Language (EIL) contexts [40]. Regardless of all
these considerations, one also has to take into account that halo effect [41]
might bias the scores to a certain extent.

Fig. 1: Heatmap of the results of the post-hoc Nemenyi test on the TLT-school
dataset.

Since every utterance was scored by only one expert, it was not possible to
evaluate any kind of agreement among experts. Note that the CEFR levels were
assigned before the tests and should be considered as expected proficiency
levels, whereas the test scores are effectively representing each learner’s
performance in the exam. We divided the data into three sets: a training set
of 322, a development set of 85 and a test set of 87 samples.

| Train | Dev | Test | Total
---|---|---|---|---
0-3 | 74 | 14 | 27 | 115
3-6 | 73 | 20 | 17 | 110
6-9 | 77 | 20 | 11 | 108
9-12 | 98 | 31 | 32 | 161
Total | 322 | 85 | 87 | 494
Table 2: Number of answers for each score range in TLT-school.

##  3 Model architectures

###  3.1 wav2vec2-based graders

Wav2vec 2.0 encodes speech audio through a multilayer convolutional neural
network (CNN). After encoding, masking is applied to spans of the resulting
latent representations, which are fed into a transformer to build
contextualised representations. Gumbel softmax is employed to compute the
contrastive loss on which the model is trained, and speech representations are
learned from this training. For our models, we initialised the configuration
and processor from a version provided by the HuggingFace Transformer Library
[42]222huggingface.co/patrickvonplaten/wav2vec2-base. After the learners’
answers are fed into the model, wav2vec 2.0 provides contextualised
representations. In order to handle representations of different audio
lengths, we used a mean pooling method to concatenate 3D representations into
2D representations. These are subsequently passed through a Dense layer of 768
units, a Dropout layer and, finally, through an output layer. We tried
different architectures and hyperparameters and, finally, we opted for those
described in the following paragraphs.

ICNALE: the task is multi-class classification, therefore the output layer has
5 units and softmax as activation function. The grader uses cross entropy as
loss function. Training uses the AdamW optimiser [43] with batch size 4,
gradient accumulation step 2, dropout 0.2, and learning rate 1e-5. The grader
is trained for 8 epochs.

TLT-school - holistic score: in this case, assessment is treated as a
regression task, therefore the output layer has 1 unit and a linear activation
function. The loss function is mean squared error (MSE). The grader is trained
for 12 epochs using AdamW optimiser with batch size 4, gradient accumulation
step 2, dropout 0.2, and learning rate 5e-5.

TLT-school - analytic subscores: we trained 6 different graders for each
aspect of proficiency. The batch size and gradient accumulation steps are the
same as the holistic grader, whereas the other hyperparameters are reported in
Table 3.

| Epochs | Learning rate | Dropout
---|---|---|---
Relevance | 13 | 5e-6 | 0.1
Correctness | 19 | 2e-6 | 0.1
Lexical | 12 | 4e-6 | 0.1
Pronunciation | 8 | 1e-5 | -
Fluency | 6 | 8e-6 | 0.1
Communicative | 10 | 1e-5 | 0.1
Table 3: Hyperparameters of the individual wav2vec2-based graders.

As has been said, the first component of wav2vec 2.0 consists of a stack of
CNN layers that are used to extract acoustically meaningful - but contextually
independent - features from the raw speech signal. This part of the model has
already been sufficiently trained during pretraining and does not need to be
fine-tuned. Therefore, we froze the parameters of the feature extractor for
all our experiments.

###  3.2 BERT-based graders

The baseline systems used for comparison are BERT-based models in the version
provided by the HuggingFace Transformer Library [42] 333huggingface.co/bert-
base-uncased. The models take a sequence of token embeddings, i.e., of the
answers provided by the learners as inputs. Each token is transformed into a
vector representation and then passed to BERT’s encoder layer. We use the
[CLS] token state and feed it to a classification or regression head,
depending on the nature of the task. Similarly to what we did with our
wav2vec2-based models, we kept the BERT layer frozen. We tried various
hyperparameters and architectures and we chose the ones described in the
following paragraphs.

ICNALE: the classification head consists of a stack of three Dense layers of
768 units, a stack of three Dense layers of 128 units, and an output layer of
5 units with softmax as activation function. The model is trained for 600
epochs with batch size set to 256 using Adam optimiser [44] with learning rate
set to 5e-5 and cross entropy as loss function. The maximum sequence length is
256.

TLT-school - holistic score: the regression head has the same intermediate
layers as the model used on ICNALE and an output layer of 1 unit with linear
activation function. The model is trained for 800 epochs on the manual
transcriptions and for 150 epochs on the ASR transcriptions. The batch size is
set to 256 and the maximum sequence length to 64. Training uses Adam optimiser
with learning rate 2e-5 and dropout 0.2.

TLT-school - analytic subscores: analogously to what we did with the
wav2vec2-based graders, we trained a BERT-based grader for each proficiency
indicator. The architecture and hyperparameters of the individual graders are
the same as the holistic grader except for the dropout rate of the fluency
grader, which is set to 0.4.

The performance of the graders trained on ICNALE is evaluated using accuracy
and weighted F1 score, whereas for the evaluation of both the holistic and
individual graders trained on TLT-school we considered Pearson’s correlation
coefficient (PCC), Spearman’s rank coefficient (SRC) and mean squared error
(MSE).

##  4 Experiments and results

###  4.1 Results on ICNALE

We started our series of experiments from the multi-class classifiers trained
on ICNALE. First, we trained and evaluated the BERT-based baseline grading
system. Secondly, we tried the wav2vec2-based grader. Table 4 shows the
results of the two graders in terms of accuracy and weighted F1 score on the
ICNALE test set. Figure 2 illustrates the confusion matrices of each CEFR
proficiency level for the two graders. As can be seen, the wav2vec-based
grader significantly outperforms the BERT baseline across all proficiency
levels. Specifically, it performs best on B1_2 and on the class related to
native speakers. While the reason for its results on the latter can be
attributed to a clear gap between L1 and low/mid levels of L2 English, such as
the ones featured in the dataset, the performance on the first is probably due
to the greater amount of training data available compared to the other CEFR
proficiency levels. Similarly but conversely, this can be inferred from the
relatively worse performance on A2, which is the least represented proficiency
level (see Table 1).

| Accuracy(%) | Weighted F1
---|---|---
BERT | 53.45 | 0.50
wav2vec2 | 77.88 | 0.77
Table 4: Results on the ICNALE test set of the BERT-based and wav2vec2-based
graders in terms of accuracy and weighted F1 score.

(a)

(b)

Fig. 2: Confusion matrices of CEFR proficiency levels for the two graders
(predicted labels on X-axis, true labels on the Y-axis) on the ICNALE test
set.

###  4.2 Results on TLT-school

We continued our experiments using the TLT-school data, starting from the
holistic graders. In this case, we compared the performance of the
wav2vec2-based model and the BERT-based model using the manual transcriptions
and the ASR output text. Table 5 reports the results on the TLT-school test
set. It can be observed that the wav2vec2-based approach significantly
outperforms BERT both on the manual and ASR transcriptions across all metrics.

| PCC | SRC | MSE
---|---|---|---
BERT-ASR | 0.749 | 0.743 | 9.877
BERT-manual | 0.857 | 0.863 | 6.110
wav2vec2 | 0.927 | 0.933 | 2.297
Table 5: Results on TLT-school test set (holistic score) of the BERT-based
grader (manual and ASR transcriptions) and the wav2vec2-based grader in terms
of PCC, SRC and MSE.

Subsequently, in order to verify the impact of our two approaches on
individual aspects of proficiency, we trained our models on the analytic
subscores of each proficiency indicator in the TLT-school test set: relevance,
formal correctness, lexical richness and complexity, pronunciation, fluency,
and communicative effectiveness. For this part of our experiments, we did not
consider the ASR transcriptions, but we only focused on the upper bound
provided by the manual transcriptions. We report the results in Table 3. As
can be seen, the use of wav2vec 2.0 significantly improves on the performance
of BERT. In particular, as expected, the more exquisitely speech-related
analytic subscores (i.e., pronunciation, fluency, and communicative
effectiveness) are best predicted using wav2vec 2.0, despite the presence of
proxies of fluency and pronunciation in the manual transcriptions.
Furthermore, the wav2vec2-based graders remarkably outperform the BERT
baselines when predicting the formal correctness subscore and the relevance
subscore. The only subscore on which the predictions of both approaches appear
fairly aligned is the one related to lexical richness and complexity, despite
the BERT baseline performance being lower in terms of PCC and MSE. This is
quite an expectable result, given that BERT is trained on a large quantity of
textual data, and lexical richness and complexity are competences that should
be typically constructed and evaluated starting from text, be it written or
spoken.

(a)

(b)

Fig. 3: Comparison of the BERT-based and wav2vec2-based individual graders in
terms of PCC and SRC.

Given such results, we continued our analysis focusing on this indicator.
Specifically, we wanted to analyse the performance of the wav2vec2- and BERT-
based graders across scores, in order to understand whether these two
approaches perform differently across proficiency and therefore could be
complementary to each other. Figure 4 illustrates the MSE variation across
scores applying Gaussian kernel smoothing with sigma set to 0.5.

We found that the BERT-based model performs moderately better in the middle
scores, whilst the approach including wav2vec 2.0 has a significantly lower
MSE for the lowest and highest scores. To verify their complementarity, we
merged them using: a) a shallow combination: we calculated the average of the
scores predicted by each grader; b) a deep combination: we concatenated the
two hidden representations in BERT and wav2vec 2.0, and we fed them through a
small network consisting of a Dense layer of 16 units, Dropout layer with
dropout rate set to 0.5 and a final output layer. The resulting network was
trained for 3000 epochs with learning rate 5e-5 and batch size 512. In both
cases, an interesting improvement across all considered metrics is observed,
as shown in Table 6.

Fig. 4: MSE variation of the wav2vec2-based and BERT-based (manual transcription) graders across scores for lexical richness and complexity. | PCC | SRC | MSE
---|---|---|---
BERT | 0.846 | 0.855 | 0.217
wav2vec2 | 0.863 | 0.851 | 0.178
BERT+wav2vec2 (shallow) | 0.885 | 0.876 | 0.164
BERT+wav2vec2 (deep) | 0.883 | 0.864 | 0.166
Table 6: Results on TLT-school test set (lexical richness and complexity) of
the BERT-based (manual transcriptions), the wav2vec2-based graders, and their
combinations in terms of PCC, SRC and MSE.

##  5 Conclusions and future works

Transcriptions of L2 speaking tests are not always easy to obtain, and, even
when they are available through ASR systems, they generally contain errors and
do not provide information about strictly speech-related aspects of
proficiency, such as pronunciation and fluency. This paper considers whether
it is possible to use wav2vec 2.0 representations to assess L2 spoken English
proficiency both holistically and analytically, even when a small quantity of
data is available. First, we found that this approach significantly
outperforms the BERT baseline system trained on manual transcriptions of the
ICNALE dataset in the task of CEFR level classification. Secondly, we
investigated the use of wav2vec 2.0 for a regression task on the B1 section of
TLT-school targeting holistic scores. In this case, we also obtained
significant improvements on the BERT baseline trained on ASR and manual
transcriptions. Finally, we tested this approach on subscores related to
individual aspects of proficiency (i.e., relevance, formal correctness,
lexical richness and complexity, pronunciation, fluency, and communicative
effectiveness) using manual transcriptions only, and we found that the
wav2vec2-based graders remarkably outperform the BERT-based baseline systems
across all proficiency indicators. For lexical richness and complexity, i.e.,
the only subscore on which the two strategies had a similar performance, we
found that two types of combination of the two models show an interesting
improvement, suggesting a certain complementarity. Future works will involve
other approaches for comparison, e.g., a grader with hand-crafted features,
and other types of combination, considering both shallow and deep fusion
methods. Furthermore, specific types of combinations will be investigated
according to each aspect of proficiency, e.g., a concatenation of the question
prompt and learner’s answer for the subscore related to relevance along the
lines of [45]. In this work, we have focused on answer-level assessment, but
we plan to investigate whether it is possible to evolve this approach in order
to provide feedback about specific parts of an answer, e.g., analysing the
local attention representations. Finally, we did not fine-tune wav2vec 2.0 on
L2 English data before using it for proficiency assessment, therefore we
envisage the investigation of this additional step to further improve our
results.

## References

  * [1] M. Zhang,  “Contrasting automated and human scoring of essays,”  R&D Connections, , no. 21, pp. 1–11, 2013.
  * [2] N. Iwashita, A. Brown, T. McNamara, and S. O’Hagan,  “Assessed levels of second language speaking proficiency: How distinct?,”  Applied Linguistics, vol. 29, no. 1, pp. 24–49, 2008.
  * [3] N.H. De Jong, M.P. Steinel, A.F. Florijn, R. Schoonen, and J.H. Hulstijn,  “Facets of speaking proficiency,”  Studies in Second Language Acquisition, , no. 34, pp. 5–34, 2012\.
  * [4] H. Strik and C. Cucchiarini,  “Automatic assessment of second language learners’ fluency,”  in Proc. International Congress of Phonetic Sciences (ICPhS) 1999, 1999.
  * [5] L. Chen, K. Evanini, and X. Sun,  “Assessment of non-native speech using vowel space characteristics,”  in Proc. 2010 IEEE Spoken Language Technology Workshop, 2010, pp. 139–144.
  * [6] E. Coutinho, F. Hönig, Y. Zhang, S. Hantke, A. Batliner, E. Nöth, and B. Schuller,  “Assessing the prosody of non-native speakers of English: Measures and feature sets,”  in Proc. 10th International Conference on Language Resources and Evaluation (LREC’16), 2016.
  * [7] S. Bhat and S. Yoon,  “Automatic assessment of syntactic complexity for spontaneous speech scoring,”  Speech Communication, vol. 67, pp. 42–57, 2015.
  * [8] P. Müller, F. De Wet, C. Van Der Walt, and T. Niesler,  “Automatically assessing the oral proficiency of proficient L2 speakers.,”  in Proc. Workshop on Speech and Language Technology for Education (SLaTE), 2009, pp. 29–32.
  * [9] S. Crossley and D. McNamara,  “Applications of text analysis tools for spoken response grading,”  Language Learning & Technology, vol. 17, no. 2, pp. 171–192, 2013\.
  * [10] Y. Wang, M.J.F. Gales, K.M. Knill, K. Kyriakopoulos, A. Malinin, R.C. van Dalen, and M. Rashid,  “Towards automatic assessment of spontaneous spoken English,”  Speech Communication, vol. 104, pp. 47–56, 2018.
  * [11] Z. Liu, G. Xu, T. Liu, W. Fu, Y. Qi, W. Ding, Y. Song, C. Guo, C. Kong, S. Yang, et al.,  “Dolphin: a spoken language proficiency assessment system for elementary education,”  in Proc. The Web Conference 2020, 2020, pp. 2641–2647.
  * [12] L. Chen, J. Tao, S. Ghaffarzadegan, and Y. Qian,  “End-to-end neural network based automated speech scoring,”  in Proc. 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2018, pp. 6234–6238.
  * [13] K. Takai, P. Heracleous, K. Yasuda, and A. Yoneyama,  “Deep learning-based automatic pronunciation assessment for second language learners,”  in Proc. International Conference on Human-Computer Interaction, 2020, pp. 338–342.
  * [14] S. Cheng, Z. Liu, L. Li, Z. Tang, D. Wang, and T.F. Zheng,  “ASR-Free pronunciation assessment,”  in Proc. Interspeech 2020, 2020, pp. 3047–3051.
  * [15] K. Kyriakopoulos, K.M. Knill, and M.J.F. Gales,  “A deep learning approach to assessing non-native pronunciation of English using phone distances,”  in Proc. Interspeech 2018, 2018, pp. 1626–1630.
  * [16] K. Kyriakopoulos, K.M. Knill, and M.J.F. Gales,  “A deep learning approach to automatic characterisation of rhythm in non-native English speech,”  in Proc. Interspeech 2019, 2019, pp. 1836–1840.
  * [17] X. Wang, K. Evanini, Y. Qian, and M. Mulholland,  “Automated scoring of spontaneous speech from young learners of english using transformers,”  in 2021 IEEE Spoken Language Technology Workshop (SLT), 2021, pp. 705–712.
  * [18] V. Raina, M.J.F. Gales, and K.M. Knill,  “Universal adversarial attacks on spoken language assessment systems,”  in Proc. Interspeech 2020, 2020, pp. 3855–3859.
  * [19] S. Bannò, B. Balusu, M. J. F. Gales, K. M. Knill, and K. Kyriakopoulos,  “View-specific assessment of L2 spoken English,”  in Proc. Interspeech 2022, 2022, pp. 4471–4475.
  * [20] Y. Lu, S. Bannò, and M.J.F. Gales,  “On assessing and developing spoken ’grammatical error correction’ systems,”  in Proceedings of the 17th Workshop on Innovative Use of NLP for Building Educational Applications (BEA 2022), 2022, pp. 51–60.
  * [21] A. Baevski, H. Zhou, A. Mohamed, and M. Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  in NeurIPS 2020, 2020, pp. 1–12.
  * [22] W.N. Hsu, A. Sriram, A. Baevski, T. Likhomanenko, Q. Xu, V. Pratap, J. Kahn, A. Lee, R. Collobert, G. Synnaeve, and M. Auli,  “Robust wav2vec 2.0: Analyzing domain shift in self-supervised pre-training,”  in Proc. Interspeech 2021, 2021, pp. 721–725.
  * [23] S.-W Yang, P.-H. Chi, Y.-S. Chuang, C.-I. J. Lai, K. Lakhotia, Y.Y. Lin, A.T. Liu, J. Shi, X. Chang, G.-T. Lin, T.-H. Huang, W.-C. Tseng, K.-T. Lee, D.-R. Liu, Z. Huang, S. Dong, S.-W. Li, S. Watanabe, A. Mohamed, and H.-Y. Lee,  “SUPERB: Speech Processing Universal PERformance Benchmark,”  in Proc. Interspeech 2021, 2021, pp. 1194–1198.
  * [24] J. Shah, Y.K. Singla, C. Chen, and R. Ratn Shah,  “What all do audio transformer models hear? Probing Acoustic Representations for Language Delivery and its Structure,”  arXiv e-prints, p. arXiv:2101.00387, 2021.
  * [25] L. Peng, K. Fu, B. Lin, D. Ke, and J. Zhan,  “A Study on Fine-Tuning wav2vec2.0 Model for the Task of Mispronunciation Detection and Diagnosis,”  in Proc. Interspeech 2021, 2021, pp. 4448–4452.
  * [26] M. Wu, K. Li, W.-K. Leung, and H. Meng,  “Transformer Based End-to-End Mispronunciation Detection and Diagnosis,”  in Proc. Interspeech 2021, 2021, pp. 3954–3958.
  * [27] X. Xu, Y. Kang, S. Cao, B. Lin, and L. Ma,  “Explore wav2vec 2.0 for Mispronunciation Detection,”  in Proc. Interspeech 2021, 2021, pp. 4428–4432.
  * [28] E. Kim, J.-J. Jeon, H. Seo, and H. Kim,  “Automatic pronunciation assessment using self-supervised speech representation learning,”  in Proc. Interspeech 2022, 2022, pp. 1411–1415.
  * [29] J. Devlin, M. Chang, L. Kenton, and K. Toutanova,  “BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding,”  arXiv e-prints, p. arXiv:1810.04805, 2018.
  * [30] S. Ishikawa,  “A new horizon in learner corpus studies: The aim of the ICNALE project,”  in Corpora and language technologies in teaching, learning and research, G. Weir, S. Ishikawa, and K. Poonpon, Eds., pp. 3–11. University of Strathclyde Press, 2011.
  * [31] Council of Europe,  Common European Framework of Reference for Languages: Learning, Teaching, Assessment,  Cambridge University Press, Cambridge, 2001.
  * [32] Z. Zhou, S. Vajjala, and S.V. Mirnezami,  “Experiments on Non-native Speech Assessment and its Consistency,”  in NLP4CALL 2019, 2019, pp. 86–92.
  * [33] S. Bannò and M. Matassoni,  “Cross-corpora experiments of automatic proficiency assessment and error detection for spoken English,”  in Proceedings of the 17th Workshop on Innovative Use of NLP for Building Educational Applications (BEA 2022), 2022, pp. 82–91.
  * [34] R. Gretter, M. Matassoni, S. Bannò, and D. Falavigna,  “TLT-school: a corpus of non native children speech,”  in Proceedings of the 12th Language Resources and Evaluation Conference, 2020.
  * [35] R. Gretter, M. Matassoni, K. Allgaier, S. Tchistiakova, and D. Falavigna,  “Automatic assessment of spoken language proficiency of non-native children,”  in IEEE International Conference on Acoustics, Speech and Signal Processing, 2019.
  * [36] M. Friedman,  “The use of ranks to avoid the assumption of normality implicit in the analysis of variance,”  Journal of the American Statistical Association, vol. 32, no. 200, pp. 675–701, 1937.
  * [37] P.B. Nemenyi,  Distribution-free multiple comparisons,  Ph.D. thesis, Princeton University, 1963.
  * [38] J.M. Levis,  Intelligibility, Oral Communication, and the Teaching of Pronunciation,  Cambridge University Press, Cambridge, 2018.
  * [39] M.C. Pennington and P. Rogerson-Revell,  English Pronunciation Teaching and Research: Contemporary Perspectives,  Palgrave Macmillan, London, 2019.
  * [40] E.-L. Low,  Pronunciation for English as an International Language: from Research to Practice,  Routledge, London, 2014.
  * [41] C.M. Myford and E.W. Wolfe,  “Detecting and measuring rater effects using many-facet rasch measurement: Part i,”  Journal of Applied Measurement, vol. 4, no. 4, pp. 386–422, 2003\.
  * [42] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz, J. Davison, S. Shleifer, P. von Platen, C. Ma, Y. Jernite, J. Plu, C. Xu, T. Le Scao, S. Gugger, M. Drame, Q. Lhoest, and A. Rush,  “Transformers: State-of-the-art natural language processing,”  in Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations, 2020, pp. 38–45.
  * [43] I. Loshchilov and F. Hutter,  “Decoupled weight decay regularization,”  in ICLR 2017, 2017.
  * [44] D.P. Kingma and J. Ba,  “Adam: a method for stochastic optimization,”  in International Conference on Learning Representations, 2014.
  * [45] Y. Qian, R. Ubale, M. Mulholland, K. Evanini, and X. Wang,  “A prompt-aware neural network approach to content-based scoring of non-native spontaneous speech,”  in 2018 IEEE Spoken Language Technology Workshop (SLT), 2018, pp. 979–986.