# Estimating Magnetic Filling Factors From Simultaneous Spectroscopy and
Photometry: Disentangling Spots, Plage, and Network

T. W. Milbourne [  Department of Physics, Harvard University, Cambridge MA
02138, USA  Center for Astrophysics — Harvard & Smithsonian, Cambridge, MA
02138, USA  D. F. Phillips Center for Astrophysics — Harvard & Smithsonian,
Cambridge, MA 02138, USA  N. Langellier Department of Physics, Harvard
University, Cambridge MA 02138, USA  Center for Astrophysics — Harvard &
Smithsonian, Cambridge, MA 02138, USA  A. Mortier Astrophysics Group,
Cavendish Laboratory, J.J. Thomson Avenue, Cambridge CB3 0HE, UK  Kavli
Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK  R. D. Haywood NASA Sagan Fellow  Center for Astrophysics —
Harvard & Smithsonian, Cambridge, MA 02138, USA  Astrophysics Group,
University of Exeter, Exeter EX4 2QL, UK  S. H. Saar Center for Astrophysics —
Harvard & Smithsonian, Cambridge, MA 02138, USA  H. M. Cegla CHEOPS Fellow,
SNSF NCCR-PlanetS  UKRI Future Leaders Fellow  Observatoire de Genève,
Université de Genève, 51 chemin des Maillettes, 1290 Versoix, Switzerland
Department of Physics, University of Warwick, Coventry, CV4 7AL, UK  A.
Collier Cameron Centre for Exoplanet Science, SUPA, School of Physics and
Astronomy, University of St Andrews, St Andrews KY16 9SS, UK  X. Dumusque
Observatoire de Genève, Université de Genève, 51 chemin des Maillettes, 1290
Versoix, Switzerland  D. W. Latham Center for Astrophysics — Harvard &
Smithsonian, Cambridge, MA 02138, USA  L. Malavolta Dipartimento di Fisica e
Astronomia “Galileo Galilei”, Università di Padova, Vicolo dell‘Osservatorio
3, I-35122 Padova, Italy  J. Maldonado INAF-Osservatorio Astronomico di
Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy  S. Thompson
Astrophysics Group, Cavendish Laboratory, J.J. Thomson Avenue, Cambridge CB3
0HE, UK  A. Vanderburg  Department of Astronomy, University of Wisconsin,
Madison, WI 53706-1507, USA  C. A. Watson  Astrophysics Research Centre,
School of Mathematics and Physics, Queen’s University Belfast, BT7 1NN,
Belfast, UK  L. A. Buchhave DTU Space, National Space Institute, Technical
University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark  M.
Cecconi INAF-Fundacion Galileo Galilei, Rambla Jose Ana Fernandez Perez 7,
E-38712 Brena Baja, Spain  R. Cosentino INAF-Fundacion Galileo Galilei, Rambla
Jose Ana Fernandez Perez 7, E-38712 Brena Baja, Spain  A. Ghedina INAF-
Fundacion Galileo Galilei, Rambla Jose Ana Fernandez Perez 7, E-38712 Brena
Baja, Spain  M. Gonzalez  INAF-Fundacion Galileo Galilei, Rambla Jose Ana
Fernandez Perez 7, E-38712 Brena Baja, Spain  M. Lodi INAF-Fundacion Galileo
Galilei, Rambla Jose Ana Fernandez Perez 7, E-38712 Brena Baja, Spain  M.
López-Morales Center for Astrophysics — Harvard & Smithsonian, Cambridge, MA
02138, USA  A. Sozzetti INAF-Osservatorio Astrofisico di Torino, via
Osservatorio 20, 10025 Pino Torinese, Italy  R. L. Walsworth Department of
Physics, University of Maryland, College Park, MD 20742, USA  Department of
Electrical and Computer Engineering, University of Maryland, College Park, MD
20742, USA  Quantum Technology Center, University of Maryland, College Park,
MD 20742, USA

###### Abstract

State of the art radial velocity (RV) exoplanet searches are limited by the
effects of stellar magnetic activity. Magnetically active spots, plage, and
network regions each have different impacts on the observed spectral lines,
and therefore on the apparent stellar RV. Differentiating the relative
coverage, or filling factors, of these active regions is thus necessary to
differentiate between activity-driven RV signatures and Doppler shifts due to
planetary orbits. In this work, we develop a technique to estimate feature-
specific magnetic filling factors on stellar targets using only spectroscopic
and photometric observations. We demonstrate linear and neural network
implementations of our technique using observations from the solar telescope
at HARPS-N, the HK Project at the Mt. Wilson Observatory, and the Total
Irradiance Monitor onboard SORCE. We then compare the results of each
technique to direct observations by the Solar Dynamics Observatory (SDO). Both
implementations yield filling factor estimates that are highly correlated with
the observed values. Modeling the solar RVs using these filling factors
reproduces the expected contributions of the suppression of convective
blueshift and rotational imbalance due to brightness inhomogeneities. Both
implementations of this technique reduce the overall activity-driven RMS RVs
from 1.64 m s-1 to 1.02 m s-1, corresponding to a 1.28 m s-1 reduction in the
RMS variation. The technique provides an additional 0.41 m s-1 reduction in
the RMS variation compared to traditional activity indicators.

exoplanets — techniques: radial velocities — Sun: activity — planets and
satellites: detection

††facilities: TNG: HARPS-N, Solar Telescope††software: HARPS-N DRS, Python:
NumPy, SciPy, scikit-learn

Corresponding Author<EMAIL_ADDRESS>
##  1 Introduction

Figure 1: A representative HMI map of the three classes of active regions
considered in this work. Spots, plage, and network are identified using the
thresholding algorithm described by Haywood et al. (2016) and MH19, with the
threshold values given by Yeo et al. (2013). This algorithm is briefly
recapped in Sec. 2.3. Image taken January 1st, 2015 at 0:0:0.00 UT.

State of the art radial velocity (RV) searches for low-mass, long-period
exoplanets are limited by the effects of stellar magnetic activity. An Earth-
mass planet in the habitable zone of a Sun-like star has an RV amplitude on
the order of 10 cm s-1. However, stellar activity processes on host stars,
such as acoustic oscillations, magnetoconvection, suppression of convective
blueshift, and long-term activity cycles, can produce signals with amplitudes
exceeding 1 m s-1. A variety of techniques exist to mitigate the effect of
these processes on the measured RVs: Chaplin et al. (2019) discuss optimal
exposure times to average out acoustic oscillations; Cegla (2019) and Meunier
et al. (2017) present strategies for mitigaing the effects of granulation; and
Aigrain et al. (2012), Rajpaul et al. (2015), Langellier et al. (2020),
Haywood et al. (2020), and numerous others discuss statistically and
physically-driven techniques for removing the effects of large-scale magnetic
regions from RV measurements.

On timescales of the stellar rotation period, the apparent radial velocity is
modulated by three main types of active regions: dark sunspots; large, bright
plage; and small, bright network regions. These different regions may be
identified using full-disk solar images, as shown in Fig. 1. Milbourne et al.
(2019) (hereafter referred to as MH19) found that the large-scale photospheric
plage contribute differently to the solar suppression of convective blueshift
than the smaller network. Failure to account for this different contribution
leads to a significant RV shift over the 800 day span their of observations.
Some of the long-term variation reported in MH19 may also be attributed to
instrumental systematics: re-reducing the HARPS-N solar data with the ESPRESSO
DRS (Pepe, F. et al., 2021; Dumusque et al., 2020) reduces this shift from 2.6
$\pm$ 0.3 m s-1 to 1.6 $\pm$ 0.5 m s-1. However, the remaining RV shift can
only be fully removed by properly accounting for network regions in the
calculated activity-driven RVs. While this analysis is possible on the Sun
using high-resolution full disk images, traditional spectroscopic activity
indicators, such as the Mt. Wilson S-Index (Wilson, 1968; Linksy & Avrett,
1970) and the derivative index $\log(R^{\prime}_{HK})$ (Vaughan et al., 1978;
Noyes et al., 1984) do not differentiate between large and small active
regions. A new activity index or combination of activity indices is therefore
necessary to successfully model the suppression of convective blueshift on
stellar targets.

In this work, we demonstrate a new technique using simultaneous spectroscopy
and photometry to estimate spot, plage, and network filling factors, and
demonstrate that these filling factors may be used to model RV variations. In
Section 2, we discuss the solar data used by our technique. An analytical
implementation of the technique is described in Section 3, and a neural
network implementation is presented in Section 4. The resulting solar filling
factors, a model of the solar RVs, and possible applications to stellar
targets are analyzed in Section 5.

##  2 Measurements

Figure 2: Time series of solar observations used in this work. From top to
bottom: Mt. Wilson and HARPS-N solar telescope observations of the calcium
S-index (red); Total Solar Irradiance (TSI) from SORCE/TIM (blue); and SDO/HMI
plage, spot, and network filling factors (black). Note that the consistent
overall shapes of the S-index, TSI, and bright filling factors, and that dips
in the TSI are coincident with peaks in the spot filling factor. Observations
are taken between April 2010 through October 2017. Note the two different reds
used in the S-index plot: the darker red points correspond to measurements by
Mt. Wilson, and bright red points are from the HARPS-N solar telescope. Note
that due to an instrumental anomaly, no TSI data is availble from SORCE/TIM
from mid July 2013 until March 2014 - we therefore do not use any times in
this period in our analysis (Kopp & Lawrence, 2005; Kopp et al., 2005)

.

###  2.1 HARPS-N/Mt. Wilson Survey

We use the HARPS-N solar telescope (Dumusque et al., 2015; Phillips et al.,
2016) measurements of the S-index, as described in MH19 and Collier Cameron et
al. (2019). The S-index quantitatively represents activity-driven
chromospheric re-emission in the Calcium II H and K lines. The presence of
spots, plage, and network all increase the S-index. The solar telescope takes
exposures every five minutes while the Sun is visible. Each measurement of the
S-index has an average precision of $2.5\times 10^{-4}$, or a fractional
uncertainty of 0.0016.

Note that HARPS-N solar telescope observations began in mid 2015. To cover the
rest of the solar cycle, we use data from the Mt. Wilson S-index survey, as
presented by Egeland et al. (2017). Observations from the two instruments
overlap between July 2015 and February 2016 (JD 2457222 and JD 2457444),
allowing us to combine these time series. The solar telescope dataset is
rescaled so that the points in the overlapping time interval have the same
mean and variance as the Mt. Wilson data from the same time interval, as
described in Haywood et al. (2020). The resulting combined dataset is shown in
the top panel of Fig. 2.

####  2.1.1 HARPS-N Solar RVs

We use the HARPS-N solar telescope’s measurements of the solar RVs to assess
our ability to model realistic RV variations using our estimated filling
factors. 111We use publicly available HARPS-N solar telescope observations
reduced using the most recent ESPRESSO pipeline, available at
https://dace.unige.ch/sun/ (Phillips et al., 2016; Collier Cameron et al.,
2019; Dumusque et al., 2020), as well as our estimated values derived from the
linear and MLP techniques. The HARPS-N RVs used in this work span the period
from July 2015 to October 2017, with exposures taken every five minutes while
the Sun is visible. Each RV measurement has an average precision of 23 cm s-1.

###  2.2 SORCE

We use the Total Irradiance Monitor (TIM) onboard the Solar Radiation and
Climate Experiment (SORCE) (Rottman, 2005; Kopp & Lawrence, 2005; Kopp et al.,
2005) to measure photometry for the whole solar cycle. The total solar
irradiance (TSI) is the solar analogue of the light curves obtained by
_Kepler_ , K2, TESS, and CHEOPS (Borucki et al., 2010; Howell et al., 2014;
Ricker et al., 2014; Cessa et al., 2017), though the Sun’s proximity means it
can be observed continuously over much longer periods. The TIM level 3 data
products are averaged over 6 hours, with a precision of 0.005 $W/m^{2}$.

We expect the overall brightness of the Sun to vary with the stellar cycle.
Its relative brightness increases with the presence of plage and network, and
decreases with the presence of spots. This modulation makes the TSI, and
stellar light curves in general, useful tools for isolating the effects of
stellar magnetic activity (Aigrain et al., 2012). The time series of the TSI
is shown in the second panel of Fig. 2.

###  2.3 SDO

We use images from the Helioseismic and Magnetic Imager (HMI) instrument
onboard the Solar Dynamics Observatory (SDO, Schou et al. 2012; Pesnell et al.
2012; Couvidat et al. 2016) to independently calculate solar filling factors.
HMI measures the 6173.3 Å iron line at six points in wavelength space using
two polarizations. From these measurements, they reconstruct the Doppler shift
and magnetic field strength along with the continuum intensity, line width,
and line depth at each point on the solar disk.

Spots, plage, and network are identified on HMI images using a simple
threshold algorithm:

  * •

An HMI pixel is considered magnetically active if the radial component of the
magnetic field is over three times greater than the expected noise floor:
$|B_{r}|>3\sigma_{B_{r}}$.

  * •

Active pixels below the intensity threshold of Yeo et al. (2013), such that
$I<0.89I_{quiet}$, are labelled as spots. Here, $I_{quiet}$ is the average
intensity of inactive pixels on a given image.

  * •

Active regions exceeding the above intensity threshold that span an area ¿ 20
micro-hemispheres (that is, 20 parts per million of the visible hemisphere),
or $60\ \rm Mm^{2}$, are labelled as plage.

  * •

Active regions exceeding the intensity threshold that span an area ¡ 20 micro-
hemispheres are labelled as network.

These calculations are explained in further detail in MH19. The resulting
filling factors for each feature are plotted in the bottom three panels of
Fig. 2. We use one HMI image taken every four hours in our analysis. The
photon noise at disk center for the magnetograms and continuum intensity for
these HMI images are $\sigma_{B_{r}}=8$ G and $\sigma_{I_{c}}=0.01\%$
respectively (Couvidat et al., 2016). This corresponds to uncertainties
$<0.1\%$ in the resulting magnetic filling factors.

Since SDO/HMI allows us to perform precise direct, independent measurements of
the three filling factors of interest, we use these results as the ”ground
truth” in our analysis.

An SDO analogue does not exist for non-solar stars, so if we wish to determine
feature-specific filling factors of stellar targets we must make indirect
estimates of the filling factors using spectroscopic and photometric data. In
the next section, we discuss two processes to do so. To mitigate the effects
of acoustic oscillations, granulation, and other short-timescale activity
process, we take daily averages of each set of observations used in our
analysis. We also interpolate the HARPS-N/Mt. Wilson observations, SORCE/TIM
observations, and SDO/HMI filling factors onto a common time grid of one
observation each day, when all three instruments have measurements.

##  3 Linear Technique

###  3.1 Modelling Irradiance Variations Using Filling Factors

In MH19, the authors reproduce the observed TSI using a linear combination of
the spot and plage filling factors. Following Meunier et al. (2010), they use
the SDO/HMI derived plage and spot filling factors, and assume that the solar
irradiance follows the Stefan-Boltzmann law for blackbodies:

$TSI=\mathcal{A}\sigma\left[(1-f_{spot}-f_{bright})T_{quiet}^{4}\right.\\\
+f_{spot}(T_{quiet}+\Delta T_{spot})^{4}\\\ +\left.f_{bright}(T_{quiet}+\Delta
T_{bright})^{4}\right]$ (1)

where $\sigma$ is the Stefan-Boltzmann constant, $\mathcal{A}=(R_{\odot}/1\
\mathrm{AU})^{2}$ is a geometrical constant relating the energy emitted at the
solar surface to the energy received at Earth, $T_{quiet}$ is the quiet Sun
temperature, $\Delta T_{spot}$ and $\Delta T_{bright}$ are the effective
temperature contrasts of spots and plage/network regions, and $f_{spot}$ and
$f_{bright}$ are the HMI spot and plage/network filling factors. Expanding as
a power series yields the following approximation:

$TSI\approx\mathcal{A}\sigma T_{quiet}^{4}\left(1+4\frac{\Delta
T_{spot}}{T_{quiet}}f_{spot}+4\frac{\Delta
T_{bright}}{T_{quiet}}f_{bright}\right)$ (2)

In MH19, the authors show that HMI observations of filling factors may be used
to reproduce SORCE TSI given temperature contrasts for plage/network features
and spots and the effective temperature of the quiet Sun. In this work, we
invert the process and use the resulting effective temperatures of each type
of active region, along with the correlations between the TSI, S index, and
filling factors demonstrated in MH19, to reproduce the observed magnetic
filling factors for each type of active region. The potential stellar
applications of this technique are discussed in more detail in Sec. 5.2.

We begin by fitting the SORCE TSI to Eq. 2 using the SDO/HMI measurements of
$f_{spot}$ and $f_{bright}$. This fit yields $T_{quiet}=5769.85\pm 0.01$ K,
$\Delta T_{spot}=-525\pm 8$ K, and $\Delta T_{spot}=46.1\pm 0.5$ K. Also note
that the solar radius varies as a function of wavelength. To be consistent
with HMI, we use $R_{\odot}=695982\pm 13$ km, the solar radius measured at
6173.3 Å (Rozelot et al., 2015).

Figure 3: A plot of the SORCE TSI versus the HARPS-N/Mt. Wilson S index. The
color of each point corresponds to the value of $f_{spot}$. We see that the S
index is highly correlated with the TSI, as expected. We may use this
correlation to estimate the plage/network filling factors on the Sun. However,
increased spot coverage results in a lower TSI value for a given S index,
which will bias our estimate of bright region filling factor. This, in turn,
will result in a less accurate estimate of the spot filling factor. The black
line shows the result of the straightforward linear fit of TSI and S index,
which is biased as described above. To isolate the plage and network driven
TSI variations, we find the 50% most densely clustered points in the above
scatter plot, and fit a line to the upper boundary of this region. (This
choice in point density is arbitrary, but the resulting best-fit line is
robust to variations in this parameter.) The resulting fit line, shown in red,
is unbiased by the presence of spots.

###  3.2 Differentiating Bright and Dark Regions

The brightness of a Sun-like star may be modulated by the long-term stellar
activity cycle (_e.g._ , the 11-year solar cycle). Since $\Delta T_{bright}>0$
and $\Delta T_{spot}<0$, we see from Eq. 2 that TSI variations (on timescales
of the rotation period) below the current value from the activity cycle must
be the result of spots. Since the Sun is plage dominated (Shapiro et al.,
2016), plage and network are the primary source of variation of the TSI and
S-index, with spots making negligible contributions to the variability of the
irradiance on timescales of the solar cycle. This is also visible in comparing
the plots of the TSI, S index, $f_{plage}$, and $f_{ntwk}$ shown in Fig. 2,
and is also discussed in detail by MH19. We may therefore use a linear
transformation of the S-index to provide an initial estimate of $f_{bright}$,
and then use the TSI to estimate $f_{spot}$. The full calculation is as
follows:

(1) We begin by assuming that the S-index is directly proportional to the
total plage and network filling factor, as shown in Fig. 3 (that is,
$f_{bright,1}=m_{1}S_{HK}+b_{1}$), and that the plage and network are the
dominant drivers of TSI variation. Our first estimate of the spot filling
factor is therefore $f_{spot,1}=0$. We then estimate the values of $m_{1}$ and
$b_{1}$ by fitting the TSI as a linear transformation of the S-index:

$TSI_{fit,1}(m_{1},b_{1})=\mathcal{A}\sigma T_{quiet}^{4}\left(1+4\frac{\Delta
T_{bright}}{T_{quiet}}\left(m_{1}S_{HK}+b_{1}\right)\right).$ (3)

Note that we have included the physical constants for normalization though
they are degenerate with the fit parameters.

It is not sufficient to perform a simple linear fit to the TSI and S index. In
the above step, we model the activity-driven variations of the TSI due to the
presence of bright regions. However, as shown in Fig. 3, the presence of spots
produces scatter in this relationship. This scatter is in one direction: as
$f_{spot}$ increases, the observed TSI for a given value of the S-index
decreases. To isolate the activity-driven TSI variations due to $f_{bright}$,
we determine the 50% most densely clustered points in Fig. 3, and fit the
upper boundary of this region. The best fit line of this upper bound gives us
the values of $m_{1}$ and $b_{1}$ used above.

(2) Next, we assume any deviation from the fit above are driven by spots,
which are not included in this model. We then make a second estimate of the
spot filling factor, $f_{spot,2}$, from the residuals to the above fit:

$f_{spot,2}=\begin{dcases}\frac{TSI-
TSI_{fit,1}(m_{1},b_{1})}{4\mathcal{A}\sigma\Delta
T_{spot}T_{quiet}^{3}},&TSI-TSI_{fit,1}\leq 0\\\ 0,&TSI-
TSI_{fit,1}>0.\end{dcases}$ (4)

Essentially, any point below the line of best fit in Fig. 3 is assumed to be
due to spot-driven brightness variations. This increases the importance of
avoiding spot-driven biases in Step 1. If a simple linear fit is used in Step
1 instead of the fit to the upper boundary described above, the presence of
spots will reduce the slope of the best-fit line, which will result in an
artificially reduced $f_{bright}$ value, and will also exclude real spot-
driven variations from our calculation of $f_{spot,2}$.

(3) We determine our final estimate of $f_{bright}$ and $f_{spot}$ by fitting
the following expression to the TSI:

$TSI_{fit,2}=\mathcal{A}\sigma T_{quiet}^{4}\left(1+4\frac{\Delta
T_{spot}}{T_{quiet}}\left(a_{2}f_{spot,2}\right)\right.\\\
+\left.4\frac{\Delta
T_{bright}}{T_{quiet}}\left(m_{2}S_{HK}+b_{2}\right)\right).$ (5)

where our estimated values of $f_{bright}$ and $f_{spot}$ are given by

$f_{bright}=m_{2}S_{HK}+b_{2}$ (6)

$f_{spot}=a_{2}f_{spot,2}$ (7)

and the parameters $a_{2}$, $m_{2}$, and $b_{2}$ are determined by the above
fit. The resulting best-fit parameters derived from the solar case are given
in Table 1. Note that we do not expect $m_{2}$ and $b_{2}$ to be very
different from the parameters $m_{1}$ and $b_{1}$ found previously, nor do we
expect $a_{2}$ to be very different from 1. However, since we exclude any
negative residuals from our estimate of $f_{spot}$ in Step 2 above, we perform
this final fit in case this excluded information changes the best-fit
parameters in any way.

Note that this technique only requires knowledge of the star’s distance and
radius, along with estimates of the spot and plage/network temperature
contrasts and the quiet star effective temperature - this means that it can be
used to estimate filling factors without prior knowledge of the filling
factors from full-disk images. If the plage and network features are only
being used to decorrelate activity-driven RV variations, only time series
correlated with the spot, network, and plage filling factors are needed, and
the above terms may be absorbed into the fit coefficients in Eqs. 3, 5—this is
discussed further in Sec. 5.2.

$m_{1}$ | $2.02\pm 0.07$
---|---
$b_{1}$ | $-0.31\pm 0.01$
$a_{2}$ | $0.9912\pm 0.0008$
$m_{2}$ | $2.072\pm 0.002$
$b_{2}$ | $-0.3169\pm 0.0003$
Table 1: Best-fit parameters for the linear filling factor estimation
technique. As expected, $m_{1}$ and $m_{2}$ are consistent within error bars,
as are $b_{1}$ and $b_{2}$. Similarly, $a_{2}$ is very close to 1, as
expected..

###  3.3 Differentiating the Network and Plage Filling Factor

In the discussion above, we extract $f_{bright}$, the combined plage and
network filling factor. However, we may consider these two separately by
adding a network term to Eq. 1:

$TSI=\mathcal{A}\sigma\left[(1-f_{spot}-f_{plage}-f_{ntwk})T_{quiet}^{4}\right.\\\
+f_{spot}(T_{quiet}+\Delta T_{spot})^{4}\\\ +f_{plage}(T_{quiet}+\Delta
T_{plage})^{4}\\\ +\left.f_{ntwk}(T_{quiet}+\Delta T_{ntwk})^{4}\right].$ (8)

Fitting this equation to the TSI using the SDO observed filling factors
reveals that the plage and network have distinct effective temperatures,
$\Delta T_{plage}=32\pm 1K$ and $\Delta T_{ntwk}=79\pm 2K$. This is consistent
with the intensity maps produced by HMI, which show that network regions are,
on average, indeed brighter than plage These temperature contrasts are
necessary for seperating the plage and network contributions to the filling
factor.

Setting Eq. 1 equal to Eq. 8 and expanding as a power series, we find

$f_{bright}\Delta T_{bright}\approx f_{plage}\Delta T_{plage}+f_{ntwk}\Delta
T_{ntwk}.$

Since areas are additive, we also expect

$f_{bright}=f_{plage}+f_{ntwk}.$

Combining these equations and solving for $f_{ntwk}$ and $f_{plage}$ in terms
of $f_{bright}$ yields the following expressions:

$f_{ntwk}=\frac{\Delta T_{plage}-\Delta T_{bright}}{\Delta T_{plage}-\Delta
T_{ntwk}}f_{bright}+\mathcal{B},\\\ $ (9)

and

$f_{plage}=\frac{\Delta T_{bright}-\Delta T_{ntwk}}{\Delta T_{plage}-\Delta
T_{ntwk}}f_{bright}-\mathcal{B}.$ (10)

Note the prefactors for each estimate, which simply rescale the brightness
contributions of each class of active region to account for the different
effective temperatures. Also note the offset $\mathcal{B}$ in our estimate of
$f_{plage}$: this accounts for the fact that $f_{plage}$ goes to 0 at solar
minimum, while $f_{ntwk}$ has a basal value at solar minimum. In this
analysis, the value of $\mathcal{B}$ may be found from the expected value of
$f_{bright}$ at solar minimum:

$\mathcal{B}=\frac{\Delta T_{bright}-\Delta T_{ntwk}}{\Delta T_{plage}-\Delta
T_{ntwk}}\textrm{min}\left(f_{bright}\right).$ (11)

Determining this offset therefore requires TSI and S index observations taken
at solar minimum, which may increase the observational load associated with
this technique. However, modeling the effects of network on the activity-
driven solar RV variations only requires a quantity correlated with
$f_{plage}$. The value of this offset is therefore unimportant for our
purposes.

##  4 Machine Learning Technique

While machine learning techniques are predicatively powerful, their black-box
nature makes them not physically explanatory, and therefore not necessarily
useful for some scientific applications. However, the existence of a clear
causal connection between the S-index, TSI, and filling factors makes machine
learning a strong candidate for the problem of estimating feature-specific
magnetic filling factors from spectroscopic and photometric information. We
already know the physics connecting these variables, and can therefore have
machine learning ”discover” and refine the relationships found above. A neural
network used as a universal function approximator (Cybenko, 1989) may be able
to determine subtle details of these relationships that are not incorporated
into our linear model, such as the different effects of network vs. plage, how
underlying spatial distributions of active regions affect the resulting
filling factors and activity indicators, and the correlations between spots
and plage.

Hidden Layer Sizes | $\alpha$ | $\beta$
---|---|---
(64,64) | 0.0001 | 0.001
Table 2: Hyperparameter values for MLP filling factor calculation as optimized
from cross-validation. Here $\alpha$ gives the $L_{2}$ regularization
parameter and $\beta$ is the learning rate.

We therefore compare the linear technique discussed in the previous section
with a type of neural network known as a Multilayer Perceptron (MLP, Hinton
1989). The MLP consists of an input layer, several fully-connected hidden
layers, and an output layer. It is one of the simplest neural networks that
may be used as a universal function approximator, making it ideal for this
application. We implement the MLP using the MLPRegressor class in the scikit-
learn package in Python (van Rossum, 1995; Pedregosa et al., 2011).

We train the MLP using the TSI and S-index inputs, and using the SDO plage,
spot, and network filling factors as outputs. 75% of the total available data
(taken over the whole solar cycle) is used for training, with 25% set aside to
test the performance of the trained network. The MLP uses two hidden network
layers, each with 64 neurons. We optimized the size and number of these layers
as well as the $L_{2}$ regularization parameter, $\alpha$, which combats
overfitting by constraining the size of the fit parameters as measured with an
$L_{2}$ norm; and the learning rate, $\beta$, which controls the step-size in
the parameter space search using five-fold cross validation. That is, we
randomly shuffled the training data and divided it into five groups. We
trained the network on four of these groups, and then tested the network on
the remaining group. We repeated this process using each of the five groups as
a test set to mitigate the effects of overfitting on our network, and then
repeated the entire five-fold process using each combination of network
parameters to determine which combination of hyperparameters resulted in the
best performance. The resulting values are summarized in Table 2. The network
was optimized to minimize square error using a stochastic gradient descent
algorithm (SGD), and was trained for a maximum of $10^{4}$ steps (though the
algorithm may stop training earlier once the network converges.)

Note that we may also use this MLP approach to fit the solar RVs directly
using the TSI and S-index, without first computing magnetic filling factors.
In Sec. 5.1, we compare a direct MLP fit of this form to RV models dervived
from our estimated filling factors to determine if there is any additional RV
information in the TSI and S index which is not incorporated into our filling
factor estimates.

##  5 Results and Discussion

Figure 4: Comparison of the SDO/HMI-measured magnetic filling factors (black)
to the machine learning (blue) and linear (orange) estimates derived from the
S-index and TSI. The time series for the three filling factors are plotted in
the left column. The estimated filling factors are plotted as a function of
the HMI filling factors in the right column—the grey dashed lines indicate a
slope of 1, and are meant to guide the eye. Both the linear and machine
learning techniques reproduce the directly-observed values of $f_{spot}$,
$f_{plage}$, and $f_{ntwk}$. Note that there is a slight offset between the
linear estimate of $f_{ntwk}$ and the SDO measurements. However, this offset
is well within the expected 20% - 50% definitional variations reported by
Meunier et al. (2010) Figure 5: 230 day subsets of the time series of the
SDO/HMI-observed magnetic filling factors (black), along with the MLP (blue)
and linear (orange) estimates derived from the S-index and TSI. Three subsets
are shown, taken during the middle of the stellar cycle (left), during solar
maximum (middle), and approaching solar minimum (right). Both techniques
successfully reproduce $f_{spot}$, $f_{plage}$, and $f_{ntwk}$, with
especially good performance at solar minimum.

Fig. 4 shows that both the linear and MLP-based techniques successfully
reproduce the directly-observed values of $f_{spot}$, $f_{plage}$, and
$f_{ntwk}$. Fig. 5 shows the same information as Fig. 4, but for three 230 day
regions taken in the middle of the solar cycle, at solar maximum, and at solar
minimum. We see that, again, both the linear and MLP techniques are able to
reproduce the SDO-measured values of $f_{spot}$, $f_{plage}$, and $f_{ntwk}$
at all points in the stellar activity cycle on these timescales.

Note that there is a systematic $\sim 0.004$ offset between the linear
estimates of $f_{ntwk}$ and the SDO measured values. This is likely the result
of the significant covariance between $b_{1}$ and $\mathcal{A}\sigma
T_{quiet}^{4}$ in Eq. 3. Any systematic errors in the measured values of
$R_{\odot}$ and $T_{quiet}^{4}$ will change the resulting value of $b_{1}$,
resulting in an offset in the estimated values of $f_{ntwk}$. Small changes to
these parameters can dramatically change the observed offset in $f_{ntwk}$:
artificially increasing $T_{quiet}^{4}$ by 0.15 K eliminates the offset
entirely. This is well below the precision achieved for measurement of stellar
temperatures. While we attain good precision in the solar case, in general
linear estimates of $f_{ntwk}$ should assumed to be true up to a constant
offset. As stated previously, using these filling factors to remove activity-
driven signals from RV measurements only requires values correlated with the
filling factor value, making this offset unimportant.

We also note that, while $R_{\odot}$ and $T_{quiet}^{4}$ are assumed to be
constants in our model, they do change in time as the result of physical
processes not included in our model. These quantities also vary with
wavelength: Since here we are using the Ca II H&K lines and integrated visible
intensity to reproduce filling factors measured at 6173.3 Å, uncertainties in
these parameters associated with their wavelength dependence are inevitable.
Indeed, Meunier et al. (2010) note that measured filling factors will vary by
20% to 50% as a result of these dependencies and other definitional
differences: our estimated $f_{ntwk}$ values are certainly consistent with the
SDO measured values within these margins.

In Table 3, we list the Pearson correlation coefficients between the HMI
derived filling factors and our estimates from the linear and MLP techniques.
(For the sake of consistency, note that for both the linear and MLP estimates,
we compute correlation coefficients only for results generated using the
fraction of data reserved for testing the MLP.) We see that both techniques
reproduce the information contained in the HMI filling factors, with the MLP
performing slightly better than the linear model on all three filling factors.
This may indicate that there is some additional information about the filling
factors present in the TSI and S index observations that is not being used by
the linear technique. However, given the high degree of correlation produced
by both techniques, we can use both estimates of the magnetic filling factors
to reduce the effects of activity on observed RVs.

| Spots | Plage | Network
---|---|---|---
Linear | 0.81 | 0.87 | 0.81
MLP | 0.85 | 0.87 | 0.89
Table 3: Pearson correlation coefficients between HMI ground-truth filling
factors and the linear and MLP estimates for each class of filling factor.

###  5.1 Application to Solar RVs

Filling Factor Source | $r\left(\Delta v_{conv},RV(f_{spot},f_{ntwk},f_{plage})\right)$
---|---
HMI | 0.92
Linear Estimate | 0.84
MLP Estimate | 0.83
Table 4: Pearson correlation coefficients between HMI derived estimate of the suppression of convective blueshift, $\Delta v_{conv}$ and the activity driven RVs derived from Eq. 13. The very high correlation coefficients indicate that the plage and network filling factors successfully estimate the RV contribution of the suppression of convective blueshift, as expected from MH19. Filling Factor Source | $r\left(\Delta v_{phot},RV(f_{spot})\right)$
---|---
HMI | 0.62
Linear Estimate | 0.70
MLP Estimate | 0.67
Table 5: Pearson correlation coefficients between HMI derived estimate of the photometric velocity shift, $\Delta v_{phot}$ and the activity driven RVs derived from Eq. 14. The relatively high correlation coefficients indicate that the spot filling factors successfully estimate the photometric RV shifts, as expected from MH19. | RMS (m s-1)
---|---
Full solar dataset | 1.64
Decorrelated with S index | 1.10
Decorrelated with HMI filling factors | 0.91
Decorrelated with linear filling factor estimates | 1.04
Decorrelated with MLP filling factor estimates | 1.02
Decorrelated with MLP RV estimate | 0.96
Table 6: RMS RV residuals from several models and methods. Using our estimates
of $f_{spot}$, $f_{plage}$, and $f_{ntwk}$ in Eq. 15 reduces the RMS RVs by 60
cm s-1. However, using the HMI-observed filling factors reduces the RMS
residuals by a further 13 cm s-1, indicating there is additional information
in these filling factors not captured by our estimates. A direct MLP fit to
the solar RVs, using the S index and TSI as inputs, performs better than our
estimated filling factors, but does not perform as well as the fit to HMI
filling factors. This indicates that, while our estimated filling factors are
highly correlated with the observed values, the S index and TSI alone are
insufficient to completely characterize the filling factors of each feature.

In MH19, the authors found that the HARPS-N solar radial velocities were well-
represented by a linear combination of $\Delta v_{conv}$, the suppression of
convective blueshift, and $\Delta v_{phot}$, the photometric velocity shift
due to bright and dark active regions breaking the symmetry of the solar
rotational profile:

$RV=A_{1}\Delta v_{phot}+B_{1}\Delta v_{conv}+RV_{0}.$ (12)

Here, we see if we can perform a similar reconstruction using our estimates of
the magnetic filling factors. Since the presence of active regions drives the
suppression of convective blueshift, we expect the $\Delta v_{conv}$ to be
proportional to the spot, plage, and network filling factors. Based on the
results of MH19, we also expect network and plage regions to have different
contributions to $\Delta v_{conv}$. We therefore model the suppression of
convective blueshift as:

$\Delta v_{conv}=Bf_{spot}+Cf_{plage}+Df_{network}+E.$ (13)

While plage and network occupy a greater area than spots on the Sun, and
therefore dominate the suppression of convective blueshift, the higher
brightness contrast of spots means that they drive the photometric RV shift,
$\Delta v_{phot}$. We expect $\Delta v_{phot}$ to scale with number and size
of the spots rotating across the solar surface. However, we also expect a
phase lag between $f_{spot}$ and $\Delta v_{phot}$. For a single spot moving
across the solar disk, $f_{spot}$ is at its maximum value when the spot is on
the center of the solar disk. However, the absolute value of $\Delta v_{phot}$
is maximized when the spot is at the solar limb, rotating toward or away from
the observer, and is zero when the spot is at disk center. We therefore expect
$\Delta v_{phot}$ to also depend on the derivative of the filling factor with
respect to time:

$\Delta v_{phot}\propto f_{spot}\times\left(\frac{df_{spot}}{dt}\right).$ (14)

Note that this formulation mirrors the FF′ method developed by Aigrain et al.
(2012).

By combining Eqs. 12, 13, and 14, we therefore produce a model of the solar
RVs based on our estimated feature-specific magnetic filling factors:

$RV=Af_{spot}\left(\frac{df_{spot}}{dt}\right)+Bf_{spot}+Cf_{plage}+Df_{network}+RV_{0}.$
(15)

Note that the offset $E$ in Eq. 13 has been absorbed into $RV_{0}$.

We then fit the HARPS-N solar telescope RVs to Eq. 15 using the directly-
measured SDO filling factors. The results of each fit is given in Table 7.
Note that since the linear estimates of $f_{bright}$, $f_{ntwk}$, and
$f_{plage}$ are linear transformations of the S index, as shown in Eqs. 6, 9,
and 10, we require $D=0$ to avoid degeneracies when using the filling factor
estimates derived from the linear technique. Fitting the linear estimated
filling factors to Eq. 15 without this constraint is equivalent to fitting to
$RV=Af_{spot}\left(\frac{df_{spot}}{dt}\right)+Bf_{spot}+C^{\prime}S_{HK}+D^{\prime}S_{HK}+E^{\prime}$.
We therefore set $D=0$ to avoid having degeneracies between $C$ and $D$ in our
fit. No such constraint on $D$ is necessary when considering the SDO or MLP
filling factors.

Filling Factor Source |  $A$ ($10^{5}$ m) |  $B$ (m s-1) |  $C$ (m s-1) |  $D$ (m s-1) |  $RV_{0}$ (m s-1)
---|---|---|---|---|---
SDO | $7.7\pm 1.2$ | $165\pm 76$ | $244\pm 15$ | $52\pm 27$ | $-2.0\pm 0.3$
Linear | $7.1\pm 1.3$ | $470\pm 80$ | $281\pm 11$ | 0∗ | $-1.32\pm 0.08$
MLP | $7.5\pm 1.4$ | $479\pm 218$ | $242\pm 52$ | $2\pm 90$ | $-1.3\pm 1.2$
Table 7: Best fit coefficients for fitting Eq. 15 to the HARPS-N solar RVs
using SDO filling factors, linear estimates of the filling factors, and MLP
estimates of the filling factors. Note that the estimates of $f_{plage}$ and
$f_{ntwk}$ derived from the linear technique are both linear transformations
of the S index. We therefore require $D=0$ to avoid degeneracies when using
the filling factor estimates derived from the linear technique.

To ensure our fit is indeed reproducing the suppression of convective
blueshift and the photometric RV shift, as expected, we compare the relevant
terms of Eq. 15 to the SDO/HMI estimates of these RV perturbations, as
calculated in MH19. In Tables 4 and 5, we compare the estimates of $\Delta
v_{conv}$ and $\Delta v_{phot}$ computed from Eqs. 13 and 14 using the filling
factors measured by SDO, and estimated using the linear and MLP techniques to
the values of $\Delta v_{conv}$ and $\Delta v_{phot}$ derived from HMI
observations in Haywood et al. (2016) and MH19. We see that all of the
estimated values of $\Delta v_{conv}$ are highly correlated with the HMI-
derived velocities. Our estimates of $\Delta v_{phot}$ are less correlated
with the actual photometric shift, but still show good agreement.
Interestingly, including the contributions of plage and network regions in Eq.
14 — that is, adding terms $\propto
f_{plage}\times\left(\frac{df_{plage}}{dt}\right)$ and $\propto
f_{ntwk}\times\left(\frac{df_{ntwk}}{dt}\right)$ — does not appear to increase
the correlation coefficient. However, we may still conclude that the RVs
calculated using Eq. 15 indeed do correspond to the combination of suppression
of convective blueshift and photometric RV shift described by Eq. 12.

For both the SDO measured and MLP-derived filling factors, we see $C>D$. This
is consistent with the idea that the denser magnetic interconnections
available in photospheric plages are more successful in inhibiting convection,
and thus convective blueshifts, than the sparser network magnetizations, as
suggested in MH19. Indeed, we see that, using MLP estimates, the network
contribution is consistent with zero, and using SDO observations, the network
contribution is only $\sim 2\sigma$ above zero.

The $B$ coefficients, which describe the spot contributions to the suppression
of convective blueshift, vary depending on the filling factors used. The
linear and MLP estimates of $f_{spot}$ receive a heavy weighting, while the
SDO weighting is about a factor of 3 smaller. The MLP estimates also have a
contribution only $\sim 2\sigma$ above zero. However, the Sun is a plage-
dominated star, and $f_{spot}$ is about a factor of 100 times smaller than
$f_{plage}$, as shown in Fig. 2. So, while the precise weighting of $f_{spot}$
varies based on the values used, in all cases their contribution to the
suppression of convective blueshift will be negligible compared to that of
$f_{plage}$. We may therefore conclude that, as suggested by MH19, plage
regions are the dominate contribution to the solar suppression of convective
blueshift, while spots are the dominant contribution to the photometric RV
shift: knowledge of the plage and spot filling factors are therefore
sufficient to reproduce $\Delta v_{conv}$ and $\Delta v_{phot}$ respectively.

Fitting to the HMI-observed filling factors reduces the HARPS-N solar RV
residuals from 1.64 m s-1 to 0.91 m s-1, as shown in Table 6. In comparison,
the usual technique of simply decorrelating the S-index from the RV
measurements (i.e., fitting $RV=AS_{HK}+B$) results in an RMS of only 1.10 m
s-1, indicating that spots, plage, and network regions have different
contributions to the S index, and have different effects on the suppression of
convective blueshift (Meunier et al., 2010, Fig.7).

Repeating this fit with both our linear and MLP estimates of $f_{spot}$,
$f_{ntwk}$, and $f_{plage}$ reduces the RMS RV to 1.04 m s-1 and 1.02 m s-1
respectively. This implies that, while our estimates are highly correlated
with the true values of the filling factors, there is additional information
in the true filling factors that is not captured by either technique,
resulting in less precise estimates of the convective blueshift and
photometric RV shifts. Interestingly, while the linear filling factor
estimates cannot distinguish the RV contributions of the spots and network, as
discussed above, the linear and MLP estimates result in similar RMS RVs. The
fact that the linear estimates reduce the RMS RVs below the level obtained
from the S index despite this limitation highlights the importance of spots in
our models of activity-driven RVs.

To see if it is possible for a more refined technique to extract further RV
information from our inputs, we fit an MLP directly to the solar RVs using the
S-index and TSI as inputs. This is similar to the technique proposed by de
Beurs et al. (2020), but replacing the residual cross-correlation function
with the S-index and TSI. The hyperparameters of this MLP are the same as
those given in Table 2. As before, we divide our data into training and test
sets, and the quoted residuals are derived from the test set. This fit results
in an RMS residual of 0.96 m s-1, (as shown in Table 6) indicating that there
is indeed more RV information to be gained from this set of observations.
Interestingly, however, while this RMS value is below the residuals obtained
from both sets of estimated filling factors, it is greater than the 0.91 m s-1
residuals obtained by using the direct SDO measurements of the filling
factors. This appears to indicate that, while the S index and TSI contain more
information than our linear and MLP estimates could obtain, they do not
contain all the information about the solar plage, spot, and network coverage.

This is unsurprising: we note that network regions can form from decaying
plage regions. Due to the geometry of the magnetic flux tubes associated with
these regions, a network region may rotate onto the limb, become a plage
region as it rotates onto disk center, and then become a network region again
as it rotates back onto the limb. The linear technique directly uses the
different temperature contrasts of network and plage to provide a useful first
pass at differentiating these regions, but does not capture these links
between them. That is, there are additional physical effects that further
complicate the relationship between photometry, spectroscopy, filling factors,
and RVs (Miklos et al., 2020). While the underlying behavior of the MLP is
unknown, it likely employs a similar, slightly more complex technique to
differentiate the two classes of regions. The magnetic intensification effect,
which strengthens lines in the presence of a magnetic field (e.g., Leroy 1962;
Stift & Leone 2003), has an RV signal which depends on the overall filling
factor, as well as a given line’s wavelength, effective Landé $g$ value, and
the magnetic field strength (Reiners et al., 2013). HMI monitors the
photospheric 6173.3 Å iron line: these wavelength-dependent effects mean that
the filling factors derived from HMI may not be consistent with those derived
from the chromospheric calcium H and K lines. The center-to-limb dependence of
the calcium H and K lines are different than the 6173.3 Å line as well, which
could lead to mismatches in the derived filling factors as a function of
rotational phase. More complicated linear and MLP-based filling factor
estimates could use spectroscopic measurements of additional absorption lines,
and photometric measurements integrated over different wavelength bands to
compensate for these effects, and to exploit different wavelength-dependent
contrasts of each feature to better separate these three classes of magentic
active regions.

The direct MLP fit to the solar RVs and its residuals are plotted in Fig. 6.
The effects of HARPS-N cryostat cold plate warm-ups, discussed in Collier
Cameron et al. (2019); Dumusque et al. (2020), are clearly visible in the fit
residuals, indicating that the MLP is not ”learning” instrumental systematics,
and that the residual RV variations below this level are likely dominated by a
combination of instrumental systematics and activity processes not reflected
by variations in the S-index and TSI. Further work is necessary to identify
these remaining activity processes, and to disentangle them from instrumental
effects.

Figure 6:  MLP fit to the HARPS-N solar telescope data. HARPS-N RVs are shown
in black, and MLP estimates of the RVs are in blue. Fit residuals are shown in
the bottom panels: HARPS-N cryostat warm-up dates (_see text_) are indicated
with black dashed lines.

###  5.2 Application to the Stellar Case

The techniques developed in this work should be applicable to Sun-like stars
with the proper observational cadence. To reproduce properly scaled filling
factors, the linear technique requires precise knowledge of radius, effective
temperature, and distance to the target, as well as the temperature contrasts
of the plage, spots, and network. The effective temperature may be calculated
spectroscopically (Buchhave et al., 2012), while the temperature contrasts may
be assumed to be Sun-like in the case of G-class stars. The stellar radius may
then be calculated photometrically, using the spectroscopic temperature as a
prior. The stellar distance may be straightforwardly determined through
parallax measurements, while the rotation period may be obtained via
photometry or through RV measurements.

Although the techniques presented here assume a plage-dominated star, it is
straightforward to rework Eqs. 3 and 5 for a spot-dominated target: in this
case, the S index is assumed to be correlated with spot-driven variations of
the TSI, with positive deviations indicating the presence of bright, plage
regions.

While properly scaling the linear estimate of $f_{plage}$ requires
observations near stellar minimum, modelling RV variations only requires time
series which are proportional to $f_{spot}$, $f_{plage}$ and $f_{ntwk}$—in
this case, this offset is unimportant, and no additional constraint is placed
on the stellar observations. Furthermore, the values of $\Delta T_{spot}$,
$\Delta T_{plage}$ and $\Delta T_{ntwk}$ may be absorbed into the fit
coefficients in Eqs. 3, 5, and 15, further simplifying matters.

The MLP machine learning technique, in contrast, requires less knowledge of
the target star: while the mathematical and physical transformations learned
by MLP are unknown to the user, the MLP is presumably learning a more
sophisticated version of the linear technique, and implicitly ”learns” the
solar values for feature temperature contrast and quiet temperature as it
identifies higher-order correlations between the TSI, S index, and filling
factors. This makes the MLP straightforward to implement when precise contrast
values are unknown. Furthermore, since the MLP uses no timing information, it
places no constraints on the observational cadence or baseline: it only
requires simultaneous photometric and spectroscopic measurements.

However, since ground-truth filling factors can only be directly measured in
the solar case, the MLP must be trained using solar data. Its stellar
application is therefore limited to Sun-like stars (that is, stars with very
similar surface filling factors as the Sun, or possibly even only solar
twins), making it less generalizable to other stellar targets. Since stars
other than the Sun cannot be resolved spatially at high resolution, assessing
just how ”Sun-like” a target needs to be for the machine learning technique to
yield meaningful results is challenging. One possibility is to generate
synthetic stellar images for targets with a variety of spectral types,
activity levels, feature contrasts, and viewing angles using SOAP 2.0
(Dumusque et al., 2014), StarSim (Herrero, Enrique et al., 2016), or a similar
platform, computing the light curve and S-index for these images, and seeing
if an MLP trained on the solar case reproduces the filling factors expected
for each image. Such an analysis is beyond the scope of this work.

##  6 Conclusions

We assess two techniques to extract spot, plage, and network filling factors
using simultaneous spectroscopy and photometry. The first technique involves a
straightforward analytical manipulation of the S-index and TSI time series,
while the second uses a neural network machine learning technique known as a
Multilayer Perceptron (MLP) trained on ground-truth filling factors derived
from full-disk solar images. Both techniques yield filling factor estimates
which are highly correlated with values derived from full-disk solar images,
with Spearman correlation coefficients ranging from 0.81 and 0.89 from each
technique.

We show that decorrelating a nearly-three-year time series of solar RVs using
HMI-observed spot, plage, and network filling factors effectively reproduces
the expected RV variations due to the convective blueshift and rotational
imbalance due to flux inhomogeneities, reducing the residual activity-driven
RVs more than the typical technique of decorrelating using spectroscopic
activity indices alone. Fitting to HMI filling factors reduces the RV RMS from
1.64 m s-1 to 0.91 m s-1, while fitting to the S-index alone results in an RMS
variation of 1.10 m s-1. Including this additional information about spots,
plage, and network thus accounts for an additional
$\sqrt{(1.10\mathrm{m~{}s}^{-1})^{2}-(0.91\mathrm{m~{}s}^{-1})^{2}}=0.62\mathrm{m~{}s}^{-1}$
of RMS variation. The filling factor estimates from both the linear and MLP
techniques offer some improvement to the RMS residuals beyond what is obtained
from only the S-index. Decorrelating with the linear estimates reduces the RMS
variation to 1.04 m $s^{-1}$, and the MLP estimated filling factors reduces
the RMS to 1.02 m s-1.

Using a MLP trained directly on the solar RVs, we reduced the RMS to 0.96 m
s-1. While this indicates that the S-index and TSI contain more RV information
than obtained by either estimate of our filling factors, it does not lower the
RMS RVs below the 0.91 m s-1 limit obtained using direct measurements of the
magnetic filling factors. This suggests that, while our initial estimates of
$f_{spot}$, $f_{plage}$, and $f_{ntwk}$ are highly correlated with the
expected value, more information is needed to fully characterize these
feature-specific filling factors. To match the performance of the HMI filling
factors, a more sophisticated version of this technique, using additional
spectral lines and photometric bands will likely be necessary.

Both the analytical and machine learning techniques may be used to extract
filling factors on other stars: the analytical technique is more widely
generalizable, but requires detailed knowledge of the star and good temporal
sampling, ideally with observations of the target at activity minimum. The
machine learning technique, in contrast, requires no additional knowledge of
the target star, and applies no constraints on the observing schedule—however,
it is only applicable to stars with very similar filling factor properties as
the Sun.

This work was primarily supported by NASA award number NNX16AD42G and the
Smithsonian Institution. The solar telescope used in these observations was
built and maintained with support from the Smithsonian Astrophysical
Observatory, the Harvard Origins of Life Initiative, and the TNG. This
material is also based upon work supported by NASA under grants No. NNX15AC90G
and NNX17AB59G issued through the Exoplanets Research Program. The research
leading to these results has received funding from the European Union Seventh
Framework Programme (FP7/2007-2013) under grant Agreement No. 313014
(ETAEARTH). The HARPS-N project has been funded by the Prodex Program of the
Swiss Space Office (SSO), the Harvard University Origins of Life Initiative
(HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of
Geneva, the Smithsonian Astrophysical Observatory (SAO), and the Italian
National Astrophysical Institute (INAF), the University of St Andrews, Queen’s
University Belfast, and the University of Edinburgh. We would like to
acknowledge the excellent discussions and scientific input by members of the
International Team, ”Towards Earth-like Alien Worlds: Know thy star, know thy
planet”, supported by the International Space Science Institute (ISSI, Bern).
This work was partially performed under contract with the California Institute
of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through
the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute
(R.D.H.). SHS is grateful for support from NASA Heliophysics LWS grant
NNX16AB79G. A.M. acknowledges support from the senior Kavli Institute
Fellowships. ACC acknowledges support from the Science and Technology
Facilities Council (STFC) consolidated grant number ST/R000824/1. X.D. is
grateful to the Branco-Weiss Fellowship for continuous support. This project
has received funding from the European Research Council (ERC) under the
European Unions’ Horizon 2020 research and innovation program (grant agreement
No. 851555). CAW acknowledges support from Science and Technology Facilities
Council grant ST/P000312/1. HMC acknowledges financial support from the
National Centre for Competence in Research (NCCR) PlanetS, supported by the
Swiss National Science Foundation (SNSF), as well as a UK Research and
Innovation Future Leaders Fellowship. We thank the entire TNG staff for their
continued support of the solar telescope project at HARPS-N.

## References

  * Aigrain et al. (2012) Aigrain, S., Pont, F., & Zucker, S. 2012, Monthly Notices of the Royal Astronomical Society, 419, 3147, doi: 10.1111/j.1365-2966.2011.19960.x
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Buchhave et al. (2012) Buchhave, L. A., Latham, D. W., Johansen, A., et al. 2012, Nature, 486, 375, doi: 10.1038/nature11121
  * Cegla (2019) Cegla, H. M. 2019, Geosciences, 9, 114, doi: 10.3390/geosciences9030114
  * Cessa et al. (2017) Cessa, V., Beck, T., Benz, W., et al. 2017, in International Conference on Space Optics — ICSO 2014, ed. Z. Sodnik, B. Cugny, & N. Karafolas, Vol. 10563, International Society for Optics and Photonics (SPIE), 468 – 476, doi: 10.1117/12.2304164
  * Chaplin et al. (2019) Chaplin, W. J., Cegla, H. M., Watson, C. A., Davies, G. R., & Ball, W. H. 2019, The Astronomical Journal, 157, 163, doi: 10.3847/1538-3881/ab0c01
  * Collier Cameron et al. (2019) Collier Cameron, A., Mortier, A., Phillips, D., et al. 2019, Monthly Notices of the Royal Astronomical Society, 487, 1082, doi: 10.1093/mnras/stz1215
  * Couvidat et al. (2016) Couvidat, S., Schou, J., Hoeksema, J. T., et al. 2016, Solar Physics, 291, 1887, doi: 10.1007/s11207-016-0957-3
  * Cybenko (1989) Cybenko, G. 1989, Mathematics of Control, Signals and Systems, 2, 303, doi: 10.1007/BF02551274
  * de Beurs et al. (2020) de Beurs, Z. L., Vanderburg, A., Shallue, C. J., et al. 2020, Identifying Exoplanets with Deep Learning. IV. Removing Stellar Activity Signals from Radial Velocity Measurements Using Neural Networks.  https://arxiv.org/abs/2011.00003
  * Dumusque et al. (2014) Dumusque, X., Boisse, I., & Santos, N. C. 2014, ApJ, 796, 132, doi: 10.1088/0004-637X/796/2/132
  * Dumusque et al. (2015) Dumusque, X., Glenday, A., Phillips, D. F., et al. 2015, ApJl, 814, L21, doi: 10.1088/2041-8205/814/2/L21
  * Dumusque et al. (2020) Dumusque, X., Cretignier, M., Sosnowska, D., et al. 2020, Three Years of HARPS-N High-Resolution Spectroscopy and Precise Radial Velocity Data for the Sun.  https://arxiv.org/abs/2009.01945
  * Egeland et al. (2017) Egeland, R., Soon, W., Baliunas, S., et al. 2017, The Astrophysical Journal, 835, 25, doi: 10.3847/1538-4357/835/1/25
  * Haywood et al. (2016) Haywood, R. D., Collier Cameron, A., Unruh, Y. C., et al. 2016, Monthly Notices of the Royal Astronomical Society, 457, 3637, doi: 10.1093/mnras/stw187
  * Haywood et al. (2020) Haywood, R. D., Milbourne, T. W., Saar, S. H., et al. 2020, Unsigned magnetic flux as a proxy for radial-velocity variations in Sun-like stars.  https://arxiv.org/abs/2005.13386
  * Herrero, Enrique et al. (2016) Herrero, Enrique, Ribas, Ignasi, Jordi, Carme, et al. 2016, A&A, 586, A131, doi: 10.1051/0004-6361/201425369
  * Hinton (1989) Hinton, G. E. 1989, Artificial Intelligence, 40, 185, doi: https://doi.org/10.1016/0004-3702(89)90049-0
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, Publications of the Astronomical Society of the Pacific, 126, 398.  http://stacks.iop.org/1538-3873/126/i=938/a=398
  * Kopp et al. (2005) Kopp, G., Heuerman, K., & Lawrence, G. 2005, Solar Physics, 230, 111, doi: 10.1007/s11207-005-7447-3
  * Kopp & Lawrence (2005) Kopp, G., & Lawrence, G. 2005, Solar Physics, 230, 91, doi: 10.1007/s11207-005-7446-4
  * Langellier et al. (2020) Langellier, N., Milbourne, T. W., Phillips, D. F., et al. 2020, Detection Limits of Low-mass, Long-period Exoplanets Using Gaussian Processes Applied to HARPS-N Solar RVs.  https://arxiv.org/abs/2008.05970
  * Leroy (1962) Leroy, J. L. 1962, Annales d’Astrophysique, 25, 127
  * Linksy & Avrett (1970) Linksy, J. L., & Avrett, E. H. 1970, Publications of the Astronomical Society of the Pacific, 82, 169.  http://stacks.iop.org/1538-3873/82/i=485/a=169
  * Meunier et al. (2010) Meunier, N., Desort, M., & Lagrange, A.-M. 2010, Astronomy and Astrophysics, 512, 3637, doi: 10.1051/0004-6361/200913551
  * Meunier et al. (2010) Meunier, N., Lagrange, A. M., & Desort, M. 2010, A&A, 519, A66, doi: 10.1051/0004-6361/201014199
  * Meunier et al. (2017) Meunier, N., Lagrange, A.-M., & Borgniet, S. 2017, A&A, 607, A6, doi: 10.1051/0004-6361/201630328
  * Miklos et al. (2020) Miklos, M., Milbourne, T. W., Haywood, R. D., et al. 2020, The Astrophysical Journal, 888, 117, doi: 10.3847/1538-4357/ab59d5
  * Milbourne et al. (2019) Milbourne, T. W., Haywood, R. D., Phillips, D. F., et al. 2019, The Astrophysical Journal, 874, 107, doi: 10.3847/1538-4357/ab064a
  * Noyes et al. (1984) Noyes, R. W., Hartmann, L. W., Baliunas, S. L., Duncan, D. K., & Vaughan, A. H. 1984, ApJ, 279, 763, doi: 10.1086/161945
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Pepe, F. et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96, doi: 10.1051/0004-6361/202038306
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Solar Physics, 275, 3, doi: 10.1007/s11207-011-9841-3
  * Phillips et al. (2016) Phillips, D. F., Glenday, A. G., Dumusque, X., et al. 2016, Proc.SPIE, 9912, 9912 , doi: 10.1117/12.2232452
  * Rajpaul et al. (2015) Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, Monthly Notices of the Royal Astronomical Society, 452, 2269, doi: 10.1093/mnras/stv1428
  * Reiners et al. (2013) Reiners, A., Shulyak, D., Anglada-Escudé, G., et al. 2013, A&A, 552, A103, doi: 10.1051/0004-6361/201220437
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 1 , doi: 10.1117/1.JATIS.1.1.014003
  * Rottman (2005) Rottman, G. 2005, Solar Physics, 230, 7, doi: 10.1007/s11207-005-8112-6
  * Rozelot et al. (2015) Rozelot, J. P., Kosovichev, A., & Kilcik, A. 2015, The Astrophysical Journal, 812, 91, doi: 10.1088/0004-637x/812/2/91
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Solar Physics, 275, 229, doi: 10.1007/s11207-011-9842-2
  * Shapiro et al. (2016) Shapiro, A. I., Solanki, S. K., Krivova, N. A., Yeo, K. L., & Schmutz, W. K. 2016, A&A, 589, A46, doi: 10.1051/0004-6361/201527527
  * Stift & Leone (2003) Stift, M. J., & Leone, F. 2003, A&A, 398, 411, doi: 10.1051/0004-6361:20021605
  * van Rossum (1995) van Rossum, G. 1995, Python tutorial, Tech. Rep. CS-R9526, Centrum voor Wiskunde en Informatica (CWI), Amsterdam.  https://ir.cwi.nl/pub/5007/05007D.pdf
  * Vaughan et al. (1978) Vaughan, A., Preston, G., & Wilson, O. 1978, PASP, 90, 267, doi: 10.1086/130324
  * Wilson (1968) Wilson, O. C. 1968, Astrophysical Journal, 153, 221, doi: 10.1086/149652
  * Yeo et al. (2013) Yeo, K. L., Solanki, S. K., & Krivova, N. A. 2013, A&A, 550, A95, doi: 10.1051/0004-6361/201220682


# Big Bang Nucleosynthesis constraints on the Energy-Momentum Squared Gravity:
The $\mathbb{T}^{2}$ model

Dukjae Jang<EMAIL_ADDRESS>School of Physics and International
Research Center for Big-Bang Cosmology and Element Genesis, Beihang
University, Beijing 100083, China  Mayukh R. Gangopadhyay
<EMAIL_ADDRESS>Centre For Cosmology and Science
Popularization (CCSP), SGT University, Gurugram, Delhi- NCR, Haryana- 122505,
India  Myung-Ki Cheoun<EMAIL_ADDRESS>Department of Physics and OMEG
Institute, Soongsil University, Seoul 156-743, Republic of Korea  Toshitaka
Kajino<EMAIL_ADDRESS>School of Physics and International Research Center
for Big-Bang Cosmology and Element Genesis, Beihang University, Beijing
100083, China,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan,
National Astronomical Observatory of Japan 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  M. Sami<EMAIL_ADDRESS> Centre For Cosmology
and Science Popularization (CCSP), SGT University, Gurugram, Delhi- NCR,
Haryana- 122505, India
Center for Theoretical Physics, Eurasian National University, Astana 010008,
Kazakhstan
Chinese Academy of Sciences,52 Sanlihe Rd, Xicheng District, Beijing.

###### Abstract

A scale-independent energy-momentum squared gravity (EMSG) allows different
gravitational couplings for different types of sources and has been proven to
have interesting implications in cosmology. In this paper, the Big Bang
Nucleosynthesis (BBN) formalism and the latest observational constraints are
being used in order to extract constraints on this class of modified gravity
models. The model has been constrained using the light element reaction rates.
Using the tight constraint from BBN on the correction term in the Friedman
equation due to EMSG, we find a significant deviation in the matter power
spectrum of the cosmic microwave background (CMB) even for a very small
allowed value of the modification parameter, thus pointing out the sensitivity
of the matter power spectrum at the small scales.

##  I Introduction

The idea of modification of Albert Einstein’s General Relativity (GR)[1, 2]
dates back to the first few months after the seminal paper published by
Einstein. The proposals were made to extend GR and incorporate it into a
larger, more unified theory. A few examples are Eddington’s theory of
connections, Weyl’s scale-independent theory, and the higher dimensional
theories of Kaluza and Klein. But almost after 108 years, the field equations
proposed by Einstein remain still the best description of how space–time
behaves on macroscopic scales. Einstein’s equations govern everything that
happens in our universe, starting from its expansion, structure formation, and
black holes to the propagation of gravitational waves. Yet, the efforts to
extend or modify GR never stopped, simply with the aim of understanding the
dynamics of dark energy (DE) and dark matter (DM). To make a comprehensive
list of such models, readers can go through the following to understand the
motives and development of such theories: [3, 4, 5, 6, 7, 8, 9]. The recent
curve ball thrown to us by the universe dubbed the $H_{0}$ tension, and the
$S_{8}$ tension also points towards the requirement of some modification in
the GR or some extension.

There exists a specific class of modified theories that permit the presence of
scalars constructed from the energy-momentum tensor $T_{\mu\nu}$ in the
action. One can see it in the $f(R,T)$ gravity, where the action involves the
scalar $T=g^{\mu\nu}T_{\mu\nu}$, which is the trace of $T_{\mu\nu}$ [10]. The
$f(R,\mathbb{T}^{2})$ model has, $\mathbb{T}^{2}\equiv T^{\mu\nu}T_{\mu\nu}$
in the action [11, 12, 13, 14]. This model inspired by phenomenological
considerations is coined as Energy-Momentum-Squared-Gravity (EMSG). A similar
term is induced on RS brane as high energy correction to Einstein equations
[15, 16, 17, 18, 19]. It should be noted that most of the modifications of
gravity, in particular scale modifications, involve extra degrees of freedom,
which need to be screened out locally by involving mechanisms such as
chameleons or Vainstein. Theories that involve chameleon screening have great
potential for late-time cosmology. However, proper screening consistent with
local gravity constraints leaves no scope for late-time acceleration caused by
large-scale modifications in this scenario. One of the interesting features of
EMSG is that, unlike most of the modified theories of gravity, it does not
involve extra degrees of freedom.

Implications of EMSG for late time acceleration have been studied in Refs.[13,
14]. The model can be constrained using observations[20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30, 31]. Chaotic inflation [32, 33] has recently been examined
in the framework of EMSG[34] and then the Primordial Black Hole and
Gravitational Wave production is being studied in [35]. It has been reported
that in this case, the model like chaotic inflation (excluded in standard
cosmology by observation) in the larger umbrella theory of EMSG falls well
within the allowed limits of Planck’18 [36].

Let us note that the stringent constraints on modified gravity models in the
early universe come from the era of Big Bang Nucleosynthesis (BBN)[37, 38, 39,
40, 41, 42, 43]. Thus, it is necessary to go for the litmus test of the BBN on
the EMSG model to make an astute comment on primordial cosmology and its
phenomenology within the framework of the given modification of gravity.

The paper is organized as follows: In section II, we have given a short
synopsis of the EMSG model and finally introduce the modified Friedmann
Equation which plays the main role in this analysis. There also the stability
criterion is also discussed. In section III, we constrained the model of EMSG
using the latest observations of BBN. Interestingly, we have shown, BBN itself
demands the negative value of the model parameter $\alpha$ which is required
for the stable solution in this framework. With the given constraints imposed
on EMSG, we calculated the matter power spectrum and reported in IV. Finally,
we conclude with our findings and future directions in this context of
modified gravity models in the last section V.

Figure 1: The cosmic expansion rate in Eq. (2) during the BBN epoch. The red-
solid, blue-dotted and black-dashed lines indicate the $H$ for
$\alpha=-10^{-12}\,{\rm cm^{6}/s^{2}/g^{2}}$, $\alpha=-10^{-13}\,{\rm
cm^{6}/s^{2}/g^{2}}$, and $\alpha=0$ (standard), respectively.

##  II Background Equations and the Model

THe EMSG model actio is given by [13, 14]:

$S=\frac{8\pi G}{2}\int d^{4}x\sqrt{-g}\left[M_{\text{\rm
p}}^{2}R-\alpha\,M_{\rm
p}^{4(1-2\beta)}(\mathbb{T}^{2})^{\beta}+2\mathcal{L}_{\rm m}\right],$ (1)

where $M_{\rm p}$ is the reduced Planck mass, $R$ is the Ricci scalar
associated with the spacetime metric $g_{\mu\nu}$, $\mathcal{L}_{\rm m}$ is
the Lagrangian density corresponding to the matter source described by the
energy-momentum tensor $T_{\mu\nu}$. In addition, $\mathbb{T}^{2}\equiv
T_{\mu\nu}T^{\mu\nu}$ is a scalar and $\alpha$ is a dimensionless constant
that determines the coupling strength of the EMSG modification.

The details of the development of the background theory of this model can be
found in [34]. The main consequence of this action is the modification of the
Friedmann equation which becomes:

$H^{2}=\frac{8\pi G}{3}\rho-\alpha\left(\frac{1}{2}p^{2}+\frac{4}{3}\rho
p+\frac{1}{6}\rho^{2}\right)$ (2)

The modification term with the constant $\alpha$ in front, remains non-
negligible deep into the radiation domination and thus can impact the BBN.
Thus, the BBN constraints become very important to look into while making any
claim in this domain. Finally, one important fact is that $\alpha<0$ is simply
from the stability criterion (taking $\beta=1$). Otherwise, the model faces
gradient instability as discussed in [34]. Not only that, from the point of
view of the ghost and gradient instability for $\beta=1$, one needs
$-2\leq\alpha\leq 0$ Interestingly, this condition also becomes very important
as we will see in the next section, BBN also demands a negative value of
$\alpha$ to match the current bounds.

##  III Constraints from BBN

Figure 2: Final abundances of D and ${\rm Y}_{p}$ as a function of $\alpha$.
We adopted the observational data of D abundance from the metal-poor
Lyman-$\alpha$ absorption, ${\rm D/H}=2.527\pm 0.030$ [49], and mass fraction
of ${\rm{}^{4}He}$ from metal-poor extra-galactic H II regions, ${\rm
Y}_{p}=0.2448\pm 0.0033$ [50]. In both panels, the red and blue boxes indicate
constrained regions by the observational data within $2\sigma$ and $4\sigma$
range, respectively.

To constrain the EMSG model on BBN, we employ the BBN calculation code [44,
45] with updated reaction rates from the JINA REACLIB Database [46]. As input
parameters, we adopt the central value of the neutron mean lifetime provided
by the Particle Data Group, $\tau_{n}=878.6\pm 0.6$ [47], and the lower limit
of the baryon-to-photon ratio, $\eta=(6.104\pm 0.058)\times 10^{-10}$, which
corresponds to the baryon density based on the $\Lambda$CDM model (TT, TE,
EE+lowE) from Planck observations of the cosmic microwave background,
$\Omega_{b}h^{2}=0.02230\pm 0.0021$ [48].

Fig. 1 illustrates the cosmic expansion rate ($H$) given by Eq. (2) for
$\alpha=-10^{-12}\,{\rm cm^{6}/s^{2}/g^{2}}$ and $\alpha=-10^{-13}\,{\rm
cm^{6}/s^{2}/g^{2}}$, comparing it with the standard formula for $\alpha=0$. A
notable deviation in the $H$ is observed in the high-temperature region.
However, as the temperature decreases, both the energy and pressure squared
terms in the correction term rapidly decrease, leading to a reduction in the
change of $H$. Consequently, in this model, $H$ becomes similar to the
standard case for temperatures $T\lesssim 10^{9}\,{\rm K}$.

Adopting the modified $H$ given by Eq. (2), we proceed with BBN calculation.
Fig. 1 depicts the final abundances as a function of $\alpha$. The negative
$\alpha$ values result in an increased $H$, leading nuclear reaction rates to
freeze out earlier. The substantial change in $H$ during the early stages
leads to a higher neutron-to-proton ratio. The increased relic neutron
abundance, in turn, increases the D abundance through the enhancement of the
${}^{1}{\rm H}(n,\gamma)^{2}{\rm H}$ reaction. The increased D abundance
further contributes to the ${}^{4}{\rm He}$ abundance via the ${}^{3}{\rm
H}(d,n)^{4}{\rm He}$ reaction. Consequently, both D and ${}^{4}{\rm He}$
abundances increase as $\alpha$ decreases. This trend aligns with findings
from other studies on the effects of a modified expansion rate on primordial
abundances [51, 52, 53, 54, 55].

From the results presented in Fig. 2, we can constrain the range of $\alpha$.
For the D/H abundance, we obtain lower limits on $\alpha$ of $-1.73\times
10^{-12}\,{\rm cm^{6}/s^{2}/g^{2}}$ ($2\sigma$) and $-2.93\times
10^{-12}\,{\rm cm^{6}/s^{2}/g^{2}}$ ($4\sigma$). For ${\rm Y}_{p}$, the
constrained region is narrower, with values of $-3.50\times 10^{-13}\,{\rm
cm^{6}/s^{2}/g^{2}}$ ($2\sigma$) and $-9.50\times 10^{-13}\,{\rm
cm^{6}/s^{2}/g^{2}}$ ($4\sigma$). For ${}^{7}{\rm Li}$, we note that there is
no consistent parameter space to explain the observational data: ${}^{7}{\rm
Li/H}=1.58\pm 0.31$ [56].

Figure 3: Theoretical prediction for the matter power spectrum in the presence
of the central values of $\alpha$ for the most stringent bound from ${\rm
Y}_{p}$ (Green), allowed by D/H abundance (Red) and the $\Lambda$CDM model
(Black).

##  IV Observational Consequences at CMB

The correction term with constant $\alpha$ creates a deviation of the matter
power spectrum $P(k,z=0)$ from the standard $\Lambda$CDM. This is depicted in
figure 3. Here, we have considered one massive neutrino with a fixed mass at
0.06eV. It is known that on large scales where $k<k_{nr}$, the matter power
spectrum $P(k,z=0)$ only depends on the matter density fraction $\Omega_{m}$
today. Here $k_{nr}$ means the free-streaming scale where the neutrinos become
non-relativistic. Whereas at the small scales $k>k_{nr}$, even the slightest
presence of the term with $\alpha$, can have a significant effect. Thus, the
matter power spectrum can be written as:

$\displaystyle
P(k,z=0)=\left<\left|\frac{\delta\rho_{CDM}+\delta\rho_{b}+\delta\rho_{SP}}{\rho_{CDM}+\rho_{b}+\rho_{SP}}\right|\right>\leavevmode\nobreak\
\leavevmode\nobreak\ $
$\displaystyle=\Omega_{m}^{-2}\left<\left|\Omega_{CDM}\delta_{CDM}+\Omega_{b}\delta_{b}+\Omega_{SP}\delta_{SP}\right|^{2}\right>\leavevmode\nobreak\
\leavevmode\nobreak\ ,$ (3)

where, $\rho_{CDM}$, $\rho_{b}$, $\rho_{SP}$ represents density fluctuation of
the cold dark matter, baryon, and the sterile particle respectively. Whereas,
$\Omega_{CDM}$, $\Omega_{b}$, $\Omega_{SP}$ represent the fractional energy
density of the cold dark matter, baryon, and the sterile particle
respectively. For the most stringent bound from ${\rm Y}_{p}$, $\alpha$ takes
the value of $-3.50\times 10^{-13}\,{\rm cm^{6}/s^{2}/g^{2}}$ ($2\sigma$). For
the D/H abundance, we obtain lower limits on $\alpha$ of $-1.73\times
10^{-12}\,{\rm cm^{6}/s^{2}/g^{2}}$ ($2\sigma$). For this value of $\alpha$,
the matter power spectrum is depicted in Fig. 3 with the green line. At large
scales, one can see there is not much deviation from $\Lambda$CDM which is
expected as the linear effects are dominant there. However, at a small scale,
one can find that we observe significant deviations even for such a small
$\alpha$. Thus, keeping the BBN constraints in mind, one can see that the
presence of $\alpha$ will affect the small-scale observations and can be
studied to check the issues in $\Lambda$CDM. The suppression we see here is a
combined effect of the massive neutrino and the negative $\alpha$ value
correction on $\delta_{CDM}$. One can compare the result with the findings of
[34] where only the CMB+BAO bounds were accounted for. Interestingly enough
even with the much smaller value of $\alpha$ compared to [34], the suppression
remains significant. On the other hand, in this analysis, we did not consider
any positive value of $\alpha$, simply from the stability conditions demand
which matches the demand from the BBN bounds as well.

##  V Conclusion

We have explored the impact of the EMSG model on BBN and its implications for
the matter power spectrum. For negative values of $\alpha$, the correction
term in the EMSG model enhances the cosmic expansion rate, depending on the
squared energy density and pressure. Constraints on the parameter $\alpha$
were obtained earlier using CMB and BAO data; here we have focused on
implications due to nucleosynthesis. Given the radiation domination due to
relativistic species during the BBN epoch, the correction term rapidly decays
over cosmic time, proportional to $T^{8}$. Consequently, the EMSG correction
specifically affects the initial stage of the BBN epoch, which leads to an
increase in the primordial abundances of D and 4He for larger $\alpha$.

Using the BBN observations, we have shown that the lower limit of $\alpha$ is
constrained to $-3.50\times 10^{-13}\,{\rm cm^{6}/s^{2}/g^{2}}$ and
$-9.50\times 10^{-13}\,{\rm cm^{6}/s^{2}/g^{2}}$ within $2\sigma$ and
$4\sigma$ ranges, respectively.

One can find that the BBN bounds are much more stringent than the combined
CMB+ BBO bounds, as reported in [22]. Of course, one can expect the bounds of
BBN to always be much more stringent, as explained. Interestingly though, even
with the tightest bounds on $\alpha\leavevmode\nobreak\ \mathcal{(}10^{-12})$,
the effect is quite prominent on the small scales of the matter power
spectrum. Of course, one path to explore in this regard is to consider if the
suppression can indicate some part of Dark Matter to be hot, following [57].
One important aspect that came of the analysis is that BBN bounds prefer the
negative value of $\alpha$ which is also the theoretical demand to avoid ghost
and gradient instability in this model with $\beta=1$.

Finally, one more time, we would argue that any modification to Einstein’s
General Relativity, which can have effects at the scales of the BBN, has to
satisfy BBN bounds, to go further, and be cultivated to be useful at a later
time. In that case, BBN will provide the tightest constraints on the model.

##  VI Acknowledgement

Work of MRG is supported by the Department of Science and Technology(DST),
Government of India under the Grant Agreement number IF18-PH-228 (INSPIRE
Faculty Award) and by the Science and Engineering Research Board(SERB), DST,
Government of India under the Grant Agreement number CRG/2020/004347(Core
Research Grant).

MS is supported by the Science and Engineering Research Board (SERB), DST,
Government of India under the Grant Agreement number CRG/2022/004120 (Core
Research Grant). MS is also partially supported by the Ministry of Education
and Science of the Republic of Kazakhstan, Grant No. AP14870191 and CAS
President’s International Fellowship Initiative(PIFI).

TK and DJ are partly supported by the National Natural Science Foundation of
China (No. 12335009), and the National Key R$\&$D Program of China
(2022YFA1602401). TK is also supported in part by Grants-in-Aid for Scientific
Research of Japan Society for the Promotion of Science (20K03958).

MKC’s work was supported by the Basic Science Research Program of the National
Research Foundation of Korea (NRF) under Grants No. 2021R1A6A1A03043957 and
No. 2020R1A2C3006177.

MRG also wants to thank the organizers of OECE,2023, Beihang University for
the hospitality received when the project was initiated.

## References

  * [1] A. Einstein, The foundation of the general theory of relativity, Annalen Phys. 49, no.7, 769-822 (1916).
  * [2] A. Einstein, Explanation of the Perihelion Motion of Mercury from the General Theory of Relativity, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys. ) 1915, 831-839 (1915).
  * [3] C. Brans and R. H. Dicke, Mach’s Principle and a Relativistic Theory of Gravitation, Physical Review 124 (1961) 925–935.
  * [4] A. H. Chamseddine and V. Mukhanov, Mimetic Dark Matter, JHEP 11 (2013) 135, [1308.5410].
  * [5] J. D. Bekenstein, Relativistic gravitation theory for the MOND paradigm, Phys. Rev. D 70 (2004) 083509, [astro-ph/0403694].
  * [6] J. W. Moffat, Scalar-tensor-vector gravity theory, JCAP 03 (2006) 004, [gr-qc/0506021].
  * [7] C. Skordis and T. Zlosnik, New Relativistic Theory for Modified Newtonian Dynamics, Phys. Rev. Lett. 127 (2021) 161302, [2007.00082].
  * [8] T. P. Sotiriou and V. Faraoni, f(R) Theories Of Gravity, Rev. Mod. Phys. 82 (2010) 451–497, [0805.1726].
  * [9] F. W. Hehl and B. Mashhoon, A Formal framework for a nonlocal generalization of Einstein’s theory of gravitation, Phys. Rev. D 79 (2009) 064028, [0902.0560].
  * [10] T. Harko, F. S. N. Lobo, S. Nojiri and S. D. Odintsov, f (R, T ) gravity, Phys. Rev. D 84 (2011) 024020, [1104.2669].
  * [11] N. Katirci and M. Kavuk, $f(R,T_{\mu\nu}T^{\mu\nu})$ gravity and Cardassian-like expansion as one of its consequences, Eur. Phys. J. Plus 129 (2014) 163, [1302.4300].
  * [12] M. Roshan and F. Shojai, Energy-Momentum Squared Gravity, Phys. Rev. D 94 (2016) 044002, [1607.06049].
  * [13] O. Akarsu, N. Katırcı and S. Kumar, Cosmic acceleration in a dust only universe via energy-momentum powered gravity, Phys. Rev. D 97 (2018) 024011, [1709.02367].
  * [14] C. V. R. Board and J. D. Barrow, Cosmological Models in Energy-Momentum-Squared Gravity, Phys. Rev. D 96 (2017) 123517, [1709.09501].
  * [15] L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 3370-3373 (1999) [arXiv:hep-ph/9905221 [hep-ph]].
  * [16] L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 4690-4693 (1999) [arXiv:hep-th/9906064 [hep-th]].
  * [17] M. Bennai, H. Chakir and Z. Sakhi, Eur. J. Phys. 9, 84-93 (2006) [arXiv:0806.1137 [gr-qc]].
  * [18] M. R. Gangopadhyay and G. J. Mathews, JCAP 03, 028 (2018) [arXiv:1611.05123 [astro-ph.CO]].
  * [19] G. Calcagni, Phys. Rev. D 69, 103508 (2004) [arXiv:hep-ph/0402126 [hep-ph]].
  * [20] E. Nazari, F. Sarvi and M. Roshan, Generalized Energy-Momentum-Squared Gravity in the Palatini Formalism, Phys. Rev. D 102 (2020) 064016, [2008.06681].
  * [21] S. Bahamonde, M. Marciu and P. Rudra, Dynamical system analysis of generalized energy-momentum-squared gravity, Phys. Rev. D 100 (2019) 083511, [1906.00027].
  * [22] O. Akarsu, N. Katirci, S. Kumar, R. C. Nunes and M. Sami, Cosmological implications of scale-independent energy-momentum squared gravity: Pseudo nonminimal interact Phys. Rev. D 98 (2018) 063522, [1807.01588].
  * [23] O. Akarsu and N. M. Uzun, Cosmological models in scale-independent energy-momentum squared gravity, Phys. Dark Univ. 40 (2023) 101194, [2301.11204].
  * [24] E. Nazari, Light bending and gravitational lensing in energy-momentum-squared gravity, Phys. Rev. D 105 (2022) 104026, [2204.11003].
  * [25] N. Nari and M. Roshan, Compact stars in Energy-Momentum Squared Gravity, Phys. Rev. D 98 (2018) 024031, [1802.02399].
  * [26] A. Kazemi, M. Roshan, I. De Martino and M. De Laurentis, Jeans analysis in energy–momentum-squared gravity, Eur. Phys. J. C 80 (2020) 150, [2001.04702].
  * [27] E. Nazari, M. Roshan and I. De Martino, Constraining energy-momentum-squared gravity by binary pulsar observations, Phys. Rev. D 105 (2022) 044014, [2201.08578].
  * [28] O. Akarsu, E. Nazari and M. Roshan, Relativistic binary systems in scale-independent energy-momentum squared gravity, 2302.04682.
  * [29] O. Akarsu, J. D. Barrow, S. C¸ ıkınto˘glu, K. Y. Ek¸si and N. Katırcı, Constraint on energy-momentum squared gravity from neutron stars and its cosmological implications, Phys. Rev. D 97 (2018) 124017, [1802.02093].
  * [30] M. Faraji, N. Rashidi and K. Nozari, Inflation in energy-momentum squared gravity in light of Planck2018, Eur. Phys. J. Plus 137 (2022) 593, [2107.13547].
  * [31] C. Ranjit, P. Rudra and S. Kundu, Constraints on Energy–Momentum Squared Gravity from cosmic chronometers and Supernovae Type Ia data, Annals Phys. 428 (2021) 168432, [2010.02753]
  * [32] A. D. Linde, Chaotic inflation, Physics Letters B 129 (1983) 177–181.
  * [33] A. D. Linde, A new inflationary universe scenario: a possible solution of the horizon, flatness, homogeneity, isotropy, and primordial monopole problems, Physics Letters B 108 (1982) 389–393.
  * [34] S. A. Hosseini Mansoori, F. Felegary, M. Roshan, O. Akarsu and M. Sami, $T^{2}$\- inflation: Sourced by energy–momentum squared gravity., Phys. Dark Univ. 42, 101360 (2023), [arXiv:2306.09181 [gr-qc]].
  * [35] S. A. Hosseini Mansoori, F. Felegray, A. Talebian and M. Sami, JCAP 08, 067 (2023) [arXiv:2307.06757 [astro-ph.CO]].
  * [36] BICEP, Keck collaboration, P. A. R. Ade et al., Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations Phys. Rev. Lett. 127 (2021) 151301, [2110.00483
  * [37] J. Bernstein, L. S. Brown, and G. Feinberg, COSMOLOGICAL HELIUM PRODUCTION SIMPLIFIED., Rev. Mod. Phys. 61 (1989) 25. I, IV [44] E. W. Kolb and M. S. Turner, The Early Universe, vol. 69. 1990.
  * [38] K. A. Olive, G. Steigman, and T. P. Walker, Primordial nucleosynthesis: Theory and observations., Phys. Rept. 333 (2000) 389–407, arXiv:astro-ph/9905320.
  * [39] R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, Big Bang Nucleosynthesis: 2015., Rev. Mod. Phys. 88 (2016) 015004, arXiv:1505.01076 [astro-ph.CO].
  * [40] J. D. Barrow, S. Basilakos, and E. N. Saridakis, Big Bang Nucleosynthesis constraints on Barrow entropy., Phys. Lett. B 815 (2021) 136134, arXiv:2010.00986 [gr-qc].
  * [41] P. Asimakis, S. Basilakos, N. E. Mavromatos, and E. N. Saridakis, Big bang nucleosynthesis constraints on higher-order modified gravities., Phys. Rev. D 105 (2022) no. 8, 084010, arXiv:2112.10863 [gr-qc].
  * [42] N. Sasankan, M. R. Gangopadhyay, G. J. Mathews and M. Kusakabe, New observational limits on dark radiation in braneworld cosmology., Phys. Rev. D 95, no.8, 083516 (2017), [arXiv:1607.06858 [astro-ph.CO]].
  * [43] K. Ichiki, M. Yahiro, T. Kajino, M. Orito and G. J. Mathews, Observational constraints on dark radiation in brane cosmology., Phys. Rev. D 66, 043521 (2002), [arXiv:astro-ph/0203272 [astro-ph]].
  * [44] L. Kawano, “Let’s go: Early universe. 2. Primordial nucleosynthesis: The Computer way,” FERMILAB-PUB-92-004-A.
  * [45] M. S. Smith, L. H. Kawano and R. A. Malaney, Astrophys. J. Suppl. 85, 219-247 (1993).
  * [46] R. H. Cyburt et al., Astrophys. J. Suppl. Ser. 189, 240 (2010).
  * [47] R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022).
  * [48] N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)].
  * [49] R. J. Cooke, M. Pettini and C. C. Steidel, Astrophys. J. 855, no.2, 102 (2018).
  * [50] E. Aver et al., Mon. Not. Roy. Astr. Soc. 510, 1, 373 (2022).
  * [51] M. Kusakabe, S. Koh, K. S. Kim and M. K. Cheoun, Phys. Rev. D 93, no.4, 043511 (2016).
  * [52] D. Jang, M. Kusakabe and M. K. Cheoun, Phys. Rev. D 97, no.4, 043005 (2018).
  * [53] N. Sasankan, M. R. Gangopadhyay, G. J. Mathews and M. Kusakabe, Int. J. Mod. Phys. E 26, no.08, 1741007 (2017) [arXiv:1706.03630 [astro-ph.CO]]
  * [54] G. J. Mathews, M. Gangopadhyay, N. Sasankan, K. Ichiki and T. Kajino, AIP Conf. Proc. 1947, no.1, 020014 (2018) [arXiv:1711.04873 [astro-ph.CO]].
  * [55] G. Mathews, M. Kusakabe, M. Gangopadhyay, T. Kajino and N. Sasankan, EPJ Web Conf. 184, 01011 (2018)
  * [56] L. Sbordone et al., Astron. Astrophys. 522, A26 (2010).
  * [57] S. Bhattacharya, S. Das, K. Dutta, M. R. Gangopadhyay, R. Mahanta and A. Maharana, Phys. Rev. D 103, no.6, 063503 (2021) [arXiv:2009.05987 [astro-ph.CO]].