# The $230$ GHz Variability of Numerical Models of Sagittarius A*
I. Parameter Surveys on Varying $R_{\rm Low}$

Ho-Sang Chan Croucher Scholar  Department of Physics and Institute of
Theoretical Physics, The Chinese University of Hong Kong, Shatin, N.T., Hong
Kong  JILA, University of Colorado and National Institute of Standards and
Technology, 440 UCB, Boulder, CO 80309-0440, USA  Department of Astrophysical
and Planetary Sciences, University of Colorado, 391 UCB, Boulder, CO 80309,
USA  Chi-kwan Chan Steward Observatory and Department of Astronomy, University
of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA  Data Science
Institute, University of Arizona, 1230 N. Cherry Avenue, Tucson, AZ 85721, USA
Program in Applied Mathematics, University of Arizona, 617 North Santa Rita,
Tucson, AZ 85721, USA  Ben S. Prather Los Alamos National Lab, Los Alamos, NM,
87545, USA  George N. Wong School of Natural Sciences, Institute for Advanced
Study, 1 Einstein Drive, Princeton, NJ 08540, USA  Princeton Gravity
Initiative, Princeton University, Princeton, NJ 08544, USA  Charles Gammie
Department of Physics, University of Illinois at Urbana-Champaign, 1110 West
Green Street, Urbana, IL 61801, USA  Department of Astronomy, University of
Illinois at Urbana-Champaign, 1002 West Green Street, Urbana, IL 61801, USA
NCSA, University of Illinois at Urbana-Champaign, 1205 W. Clark St., Urbana,
IL 61801, USA  Illinois Center for the Advanced Study of the Universe,
University of Illinois at Urbana-Champaign, 1110 West Green St., Urbana, IL
61801, USA

###### Abstract

The $230$ GHz lightcurves of Sagittarius A* (Sgr A*) predicted by general
relativistic magnetohydrodynamics and ray-tracing (GRRT) models in Event
Horizon Telescope Collaboration et al. (2022) have higher modulation index
$M_{\Delta T}$ compared to observations. In this series of papers, we explore
the origin of such large brightness variability. In this first paper, we
performed large GRRT parameter surveys that span from the optically thin to
the optically thick regimes, covering $R_{\rm Low}$ from $1$ to $60$. We find,
depending on the model parameters that _i_) increasing $R_{\rm Low}$ to a
higher value could systematically reduce $M_{\Delta T}$, with $M_{\Delta T}$
consistent with the observed variability of Sgr A* in some cases; and _ii_)
increasing $R_{\rm Low}$ would make $M_{\Delta T}$ increase to a higher value.
Our analysis of GRRT image snapshots unravels the large $M_{\Delta T}$ for the
$R_{\rm Low}=1$ models mainly comes from the photon rings. However, secondary
contributions from the accretion flow are also visible depending on the spin
parameter. Our work demonstrates the importance of the electron temperature
used for modelling radiatively inefficient accretion flows and places new
constraints on the ion-electron temperature ratio. A more in-depth analysis
for understanding the dependencies of $M_{\Delta T}$ on $R_{\rm Low}$ will be
performed in subsequent papers.

Radio astronomy(1338) — Radiative transfer(1335) — High energy
astrophysics(739) — Plasma astrophysics(1261) — Black hole physics(159) —
Black holes(162) — Galactic center (565)

††journal: ApJ

##  1 Introduction

The nature of the centre of our Galaxy, Sagittarius A* (Sgr A*), was a mystery
for decades. Although radio (Lo et al., 1993) and X-ray observations (Baganoff
et al., 2003), as well as stellar proper motion data (Eckart et al., 2002),
suggested that it could be a supermassive black hole, direct evidence was
lacking. It was not until the recent progress made by the Event Horizon
Telescope (EHT) Collaboration (Akiyama et al., 2022) that the mystery of Sgr
A* was unveiled. By directly imaging the photon ring, it is now confirmed that
Sgr A* is a low-luminosity accreting supermassive black hole.

Sgr A* is estimated to be accreting at a sub-Eddington rate of $\sim
10^{-5}\,M_{\odot}\mathrm{year}^{-1}$ (Quataert et al., 1999; Quataert &
Gruzinov, 2000). The low luminosity of Sgr A* indicates that it falls into the
category of radiatively inefficient accretion flows (RIAF) — the accretion
disc is believed to be geometrically thick and optically thin while
transporting energy to the black hole mainly through advection. Under such
circumstances, electrons are decoupled from the ions, and the plasma attains a
two-temperature state (Shapiro et al., 1976), because _i_) the electrons are
unable to cool efficiently by radiation and _ii_) Coulomb coupling between the
ions and the electrons is weak (Begelman, 2014). In addition, electrons could
attain a non-thermal (Özel et al., 2000) number density distributions
presumably through acceleration by shocks (Sironi & Spitkovsky, 2011),
magnetic reconnections (Werner et al., 2018), and/or turbulence (Zhdankin et
al., 2019), and electrons are unable to be thermalized easily due to weak
Coulomb coupling and low cooling efficiency. When modelling the emission from
RIAF, accurate electron distributions plus its precise temperature are
required, but both are still poorly understood111Note that some recent
progress has been made in better understanding the electron distributions in
RIAF. See, for instance, Scepi et al. (2022). To the first-order
approximation, it is reasonable to parameterize an energy partition between
the electrons and ions. We assume the ion temperature $T_{i}$ is a fraction
$R$ of the electron temperature $T_{e}$. The fraction $R$ could depend on the
local environment and/or global accretion flow properties.

Earlier work, such as that by Dexter et al. (2009), assumed $R$ was constant
over the entire domain, producing reasonable agreement with the millimetre
VLBI visibility of Sgr A*. However, it was soon found to overproduce Near-
Infrared and X-ray fluxes (Mościbrodzka et al., 2009). Following this, Chan et
al. (2015) proposed a prescription function where $R$ is a step-function of
$\beta$ at some threshold $\beta_{\rm Crit}$. This step-function prescription
helps clearly distinguish emissions coming from the disk and the jet, where
$\beta$ differs drastically. It was motivated by the fact that the electron
temperature can drastically differ from that of the ion when the density is
low, and the magnetic field is strong. Using this prescription function, Chan
et al. (2015) obtained parameter spaces where suites of accreting black hole
models agree with observations.

The electron and ion temperature partitioning used by Event Horizon Telescope
Collaboration et al. (2022) in their theoretical studies of Sgr A* reads:

$R=\frac{T_{i}}{T_{e}}=\frac{R_{\rm High}b^{2}+R_{\rm Low}}{b^{2}+1},$ (1)

where $\beta=2P/|B|^{2}$ is the plasma magnetization, and $b=\beta/\beta_{\rm
Crit}$. This prescription was first used by Mościbrodzka et al. (2016), where
$R_{\rm High}$ ($R_{\rm Low}$) represent the ion-electron temperature ratios
in the weakly (strongly) magnetized regime. This function smooths the sharp
transition of the step-function prescription by Chan et al. (2015). It was
motivated by the results of particle-in-cell simulations, which suggest that
collisionless plasma preferentially heat the ions for $\beta>1$ (see Akiyama
et al., 2019, and references therein). The power-law dependencies on $b$
guarantee that radiations from strong and weak ion-electron coupling regions
are easily distinguishable.

In Event Horizon Telescope Collaboration et al. (2022), both $\beta_{\rm
Crit}$ and $R_{\rm Low}$ are fixed at $1$. This simple description yielded
powerful results for accreting Sgr A* images that passed almost all
constraints, such as the image size, visibility amplitude morphology, M-rings
fits, infrared and near-infrared fluxes, etc., except for the $230$ GHz flux
variability. The variability is computed as $M_{\Delta T}=\sigma_{\Delta
T}/\mu_{\Delta T}$, which is the ratio of the standard deviation to the mean
of the $230$ GHz flux over a time interval $\Delta
T=530\,GMc^{-2}=3\,\mathrm{hours}$. The constraint of $M_{\Delta T}$ for Sgr
A* could not be met by most of the fiducial models in Event Horizon Telescope
Collaboration et al. (2022). The reason why the $230$ GHz flux of the GRMHD
and GRRT models is more variable than observations remains a mystery.

Figure 1: Visualization of Equation (1) with increasing $R_{\rm Low}$. In this
case, $R_{\rm High}$ is fixed at $160$, and we assume $\beta_{\rm Crit}=1$.
The blue-shaded region indicates where $R$ is sensitive to variations in
$\beta$.

Here, we investigate whether varying $R_{\rm Low}$ can reduce the variability
of the $230$ GHz flux of the theoretical models of Sgr A*. The motivation for
varying $R_{\rm Low}$ can be understood in two ways. First, it is the only
parameter that has not been studied in-depth before, and it is natural to
investigate, with a wider parameter space, if it has any effect on $M_{\Delta
T}$. Second, we can consider the functional form of Equation 1, which is
visualized in Figure 1. We find that $R$ undergoes a sharp transition around
$\beta=10^{-1}$ and $10^{1}$. Given that the major contributions to the
electromagnetic emissions at $230$ GHz are mostly plasma located near the
horizon (Dexter et al., 2020), where $\beta$ has a similar order of magnitude,
a slight fluctuation in $\beta$ could lead to a drastic change in $R$, and
thus $T_{e}$, potentially inducing high variability in the electromagnetic
emission. By varying $R_{\rm Low}$, one can narrow the gap between $R_{\rm
High}$ and $R_{\rm Low}$, making $R$ less sensitive to $\beta$ around
$10^{-1}$ and $10^{1}$, and potentially solving the variability problem. These
are the reasons that motivate us to consider the variation of $R_{\rm Low}$ in
the GRRT of GRMHD models for Sgr A*. We note that Event Horizon Telescope
Collaboration et al. (2022) increased $R_{\rm Low}$ up to $10$ and found no
systematic reduction in the $230$ GHz variabilities. Thus, to extend their
work, we will further increase $R_{\rm Low}$ to explore a wider range of
parameter spaces. This paper is structured as follows: In Section 2, we
describe our methods, namely the GRMHD simulation libraries of Sgr A*, the
tools for GRRT parameter surveys, and the parameter spaces that we consider.
In Section 3, we present the results of our parameter searches, showing how
$M_{\Delta T}$ changes as $R_{\rm Low}$ is increased for theoretical models of
Sgr A*. In fact, we find reduction of $M_{\Delta T}$, and we will compare the
reduced value with observations. In the same section, we will also unravel why
GRMHD models of Sgr A* are more variable than observations. Finally, we
conclude our studies in Section 4.

##  2 Methodology

This section describes our methods of obtaining the GRMHD models, performing
the GRRT parameter surveys, and constructing the GRRT image libraries of Sgr
A* with varying $R_{\rm Low}$.

###  2.1 GRMHD Simulations

Our non-radiative GRMHD models of Sagittarius A* are obtained using the GRMHD
code kharma, which is the GPU-enabled version of iharm3d (Prather et al.,
2021). kharma solves the ideal GRMHD equations (with $c=1$) and the
divergence-less constraint of magnetic fields using the flux-constrained
transport scheme (Tóth, 2000). kharma assumes the ‘spherical-polar’ version of
the Kerr-Schild coordinates with the metric tensor given as:

$g^{\mu\nu}=\begin{pmatrix}-1+2r/\Sigma&2r/\Sigma&0&-2ar\text{sin}^{2}\theta/\Sigma\\\
2r/\Sigma&1+2r/\Sigma&0&-a\text{sin}^{2}\theta(1+2r/\Sigma)\\\ 0&0&\Sigma&0\\\
-2ar\text{sin}^{2}\theta/\Sigma&-a\text{sin}^{2}\theta(1+2r/\Sigma)&0&\text{sin}^{2}\theta(\Sigma+a^{2}\text{sin}^{2}\theta[1+2r/\Sigma])\end{pmatrix}$
(2)

where $r$ and $\theta$ are spherical distances and polar angles.
$\Sigma=r^{2}+a^{2}\text{cos}^{2}\theta$, where $a$ is the black-hole spin.
kharma uses the WENO reconstruction method (Jiang & Shu, 1996) to reconstruct
primitive variables and the Lax-Friedrichs solver to solve the Riemann problem
on cell boundaries. kharma evolves the GRMHD equations using the 2nd-order
Strong Stability Preserving Runge-Kutta Method (Gottlieb et al., 2011) with a
CFL number of $0.7$. We refer readers to Wong et al. (2022) to look for the
detailed GRMHD equations being solved and their technical details.

The kharma simulation libraries for the EHT campaign include both the Standard
and Normal Evolution (SANE) and Magnetic Arrested Disc (MAD) initial
conditions. In this work, we consider only the MAD state. The accretion disc
assumed an ideal gas equation of state with an adiabatic index of $4/3$, with
the initial conditions being a geometrically thick torus computed using the
method by Fishbone & Moncrief (1976). The initial torus has an inner radius at
$20$ and a pressure maximum at $41$, all in units of $GMc^{-2}$. To initialize
the MAD state, we set the initial magnetic field within the torus using the
vector potential (Wong et al., 2022):

$A_{\phi}=\text{max}\left[\frac{\rho}{\rho_{\rm
max}}\left(\frac{r}{r_{0}}\text{sin}\theta\right)e^{-r/400}-0.2,0\right],$ (3)

where $\rho_{\rm max}$ is the maximum density and $r_{0}$ is the inner radial
boundary of the computational domain. Magnetic fields are computed using
$\vec{B}=\nabla\times\vec{A}$. Finally, the simulation library employed the
funky modified Kerr-Schild coordinates (Prather, 2022) to overcome time step
restrictions along the pole. The radial outer boundary is located at
1,000$\,GMc^{-3}$, the simulation resolution is $288\times 128\times 128$, and
the simulation duration is 30,000 $GMc^{-2}$.

###  2.2 GRRT

We compute the images and luminosities of the GRMHD models using the open-
source GRRT code IPOLE (Mościbrodzka & Gammie, 2018), which is based on the
GRRT code GRTRANS (Dexter, 2016). In principle, IPOLE solves the polarized
radiative transfer code, solving for the Stokes parameters ($I$, $Q$, $U$,
$V$). However, the $Q$, $U$, and $V$ parameters are out of scope in this
study. Thus, to reduce computational time cost, we solve only the unpolarized
radiative transfer equation built-in IPOLE. The unpolarized, covariant
radiative transfer equation reads:

$\frac{d}{d\lambda}\left(\frac{I_{\nu}}{\nu^{3}}\right)=\frac{j_{\nu}}{\nu^{2}}-\alpha_{\nu}\frac{I_{\nu}}{\nu^{2}},$
(4)

where the subscript $\nu$ means the quantity for a given photon frequency,
$I_{\nu}$ is the intensity, $j_{\nu}$ is the emissivity, and $\alpha_{\nu}$ is
the absorptivity. $\lambda$ is the affine parameter along the photon geodesic:

$\displaystyle\frac{dx^{\nu}}{d\lambda}$ $\displaystyle=k^{\nu},$ (5)
$\displaystyle\frac{dk^{\nu}}{d\lambda}$
$\displaystyle=-\Gamma^{\nu}_{\mu\sigma}k^{\mu}k^{\sigma}.$ (6)

where $k^{\nu}$ is the photon four-vector. We use the transfer coefficients
presented in Dexter (2016) to solve Equation 4. Following the work of Event
Horizon Telescope Collaboration et al. (2022), we fix the black-hole image’s
field of view (FOV) to be $200$. The resolution of the image is taken as
$400\times 400$, which is fine enough to ensure convergence. The distance from
Earth to Sgr A* is taken as 8,178 parsec, while the mass of Sgr A* is assumed
to be $4.154\times 10^{6}\,M_{\odot}$.

Figure 2: Example contour representations of $M_{\Delta T}$ against $R_{\rm
Low}$ and $M_{\rm Unit}$ obtained through GRRT parameter surveys. This is for
the GRMHD models with $a=+0.94$. The upper right box in each sub-grid shows
the $\theta$ and $R_{\rm High}$ assumed for a particular set of GRRT. Also,
the black-dotted contour lines represent models of Sgr A* where the time-
averaged flux is $2.4\,\textrm{Jy}$, which is the observed averaged luminosity
of Sgr A*. For $R_{\rm High}=10$ and $50$, the slashed region represent the
parameter spaces where $R_{\rm High}<R_{\rm low}$, which is out of our
interest. The colour bar is scaled so that the middle (white band) represents
$M_{\Delta T}=0.1$, the upper limit of the historical distribution of
$M_{\Delta T}$ for Sgr A*. We also overlay grey-dotted contour lines to
represent regions where $M_{\Delta T}=0.1$. The contour plots for other spins
are included in Figure Set 1.

Fig. Set1. $M_{\Delta T}$ against $R_{\rm Low}$ and $M_{\rm Unit}$ for
$a=+0.5,0,-0.5,$ and $-0.94$.

Figure 3: $M_{\Delta T}$ against $R_{\rm Low}$ for the black-hole models that
satisfy the $2.4\,\textrm{Jy}$ constraint for all the black-hole spin we
considered. Each sub-figure represents results for different $a$. The sub-grid
plot in each sub-figure shows the variations of $M_{\Delta T}$ against $R_{\rm
Low}$ for a given $\theta$ across different $R_{\rm High}$. Except for $R_{\rm
High}=10$, the shape of the $M_{\Delta T}$ versus $R_{\rm Low}$ curves are
less sensitive across different $R_{\rm High}$ than $\theta$. The values of
$M_{\Delta T}$ and $R_{\rm Low}$ are obtained from the parameter search
results using Bivariate spline approximation provided by Scipy in the log-log
scale.

Figure 4: Example of time-averaged, pixel-wise luminosity plots. These images
are computed with the GRMHD snapshots within the time interval of
$\tau=(\textrm{29,465}$ – $\textrm{29,995})\,GMc^{-2}$, while spanning
different $R_{\rm High}$ and $\theta$. Here, $a=+0.94$, and we assume $R_{\rm
Low}=1$. The upper right box in each sub-grid shows the $\theta$ and $R_{\rm
High}$ assumed. The colour map is scaled down to $[0,1]$ in each sub-plot. The
luminosity plots for other spins are included in Figure Set 2.

Fig. Set2. Time-averaged, pixel-wise luminosity across different $R_{\rm
High}$ and $\theta$ (with $R_{\rm Low}=1$) for $a=+0.5,0,-0.5,$ and $-0.94$.

Figure 5: Same as Figure 4, but for the pixel-wise, time-domain standard
deviation of the luminosity. The standard deviation plots for other spins are
included in Figure Set 3.

Fig. Set3. Pixel-wise, time-domain standard deviation of the luminosity across
different $R_{\rm High}$ and $\theta$ (with $R_{\rm Low}=1$) for
$a=+0.5,0,-0.5,$ and $-0.94$.

###  2.3 Model Parameters

Table 1: List of parameters of the GRMHD models and the GRRT parameter surveys we considered.  Parameters | Range
---|---
Black-hole Spin $a$ | ($0.94$, $0.5$, $0$, $-0.5$, $-0.94$)
RHigh | ($10$, $50$, $90$, $120$, $180$)
Inclination $\theta$ (Deg) | ($0.1$, $45$, $90$, $135$, $179.9$)
$R_{\rm Low}$ | ($1$ – $\text{MIN}[60,R_{\rm High}]$)
$M_{\rm Unit}$ | ($10^{17}$ – $10^{21}$)

To investigate the effect of varying $R_{\text{Low}}$ to $M_{\Delta T}$, we
perform GRRT parameter surveys on the open science grid (Pordes et al., 2007;
Sfiligoi et al., 2009; OSG, 2006), which is a set of pools of shared computing
and data capacity for distributed high-throughput computing, supported by the
National Science Foundation and the Department of Energy. Table 1 lists the
parameters we considered in this study. In particular, we consider the time
interval $\tau=(\textrm{29,465}$ – $\textrm{29,995})\,GMc^{-2}$ of snapshots,
which covers a duration of $530\,GMc^{-2}\sim 3\,\mathrm{hours}$, to compute
$M_{\Delta T}$. This interval represents the very last sets of snapshots of
the simulations, in which inflow equilibrium near the horizon, where the $230$
GHz emissions are dominated, shall be established. In this study, we consider
$10$ sets of $R_{\rm Low}$ for a given $R_{\rm High}$, ranging from $R_{\rm
Low}=1$ to $R_{\rm Low}=\text{MIN}(60,R_{\rm High})$. We set an upper limit
for $R_{\rm Low}$ because the black-hole luminosity drops quickly at large
$R_{\rm Low}$, where the highly magnetized electron becomes cold.
Additionally, as we will show in the later section, an increase in $R_{\rm
Low}$ accompany an increase in $M_{\rm Unit}$. Thus, $R_{\rm Low}$ should not
increase further, or the assumption of Sgr A* being optically thin will be
broken. In general, GRMHD simulations of accreting systems are scale-
invariant, so there is a free parameter $M_{\rm Unit}$ that describes the
density scale of the system. We consider $10$ sets of $M_{\rm Unit}$ ranging
from $10^{17}$ to $10^{21}$ (in the log scale), which spans the optically thin
to optically thick regime. Altogether, given a certain set of $a$, $R_{\rm
High}$, $\theta$, and a particular GRMHD snapshot, we would have $10\times
10=100$ of GRRT images spanning different $R_{\rm Low}$ and $M_{\rm Unit}$.

After having the black-hole luminosities at $230$ GHz for a certain set of
$a$, $R_{\rm High}$, $\theta$, and in our considered time interval $\tau$, we
can then compute the time-domain average, variance, and hence, the modulation
index for the $230$ GHz luminosity. We could then infer how $R_{\rm Low}$
should vary against $M_{\rm Unit}$ for black-hole models that satisfy the
constraint of the Sgr A* $230$ GHz mean luminosity, which is 2.4 Jy (Wielgus
et al., 2022). In addition, we would analyze GRRT images satisfying the
luminosity constraint, and from the combined analysis, we can hopefully answer
the following:

  1. 1.

Can varying $R_{\rm Low}$ reduce the $230$ GHz variability?

  2. 2.

Why are the current GRMHD models too variable at $230$ GHz compared to
observations?

Note, however, that we assumed the vertical axis of the camera aligns with the
vertical axis of the spinning black hole.

##  3 Results and Discussion

This section presents and discusses the results of our GRRT parameter surveys,
the dependencies of $M_{\Delta T}$ on $R_{\rm Low}$ and $M_{\rm Unit}$, and
the origin of the high variabilities for the $R_{\rm Low}=1$ images.

###  3.1 Parameter Surveys

We first present the results from our parameter surveys. Here, we include
results for the GRMHD model of $a=+0.94$ in Figure 2 as a reference, while the
remaining results are shown in Figure Set 1. The black-dotted lines in the
contour plots of Figure 2 represent models consistent with the Sgr A* $230$
GHz luminosity. The grey-dotted lines indicate regions of the parameter spaces
where $M_{\Delta T}\leq 0.1$. This is also the upper limit of the historical
distribution of $M_{\Delta T}$ for Sgr A*. Depending on the black-hole spins,
the parameter space where $M_{\Delta T}\leq 0.1$ could be consistent with the
Sgr A* luminosity constraint. For instance, black holes with $a=+0.94$ ad
$a=+0.5$ show a limited sequence of Sgr A* models that could have $M_{\Delta
T}\leq 0.1$, which are obtained by increasing $R_{\rm Low}$ beyond $1$. Non-
spinning black-hole consistent with Sgr A* all have $M_{\Delta T}\leq 0.1$,
irrespective of the value of $R_{\rm Low}$. We also note that all $R_{\rm
High}=10$ Sgr A* models with $a=+0.94$ and $+0.5$ could not attain $M_{\Delta
T}\leq 0.1$ irrespective of $R_{\rm Low}$ assumed.

Black-hole with negative spins, however, differ from those of the non-negative
ones. Except for $R_{\rm High}=120,180$ viewed at $\theta=0.1\,\mathrm{Deg}$,
black-hole with $a=-0.5$ that have $M_{\Delta T}\leq 0.1$ are inconsistent
with the Sgr A* luminosity constraint. In fact, we find that no parameter
spaces are consistent with the observed $M_{\Delta T}$ for Sgr A* with
$a=-0.94$. This suggests a substantial reduction in $M_{\Delta T}$ with
varying $R_{\rm Low}$ is model-dependent and sensitive to the black-hole
parameters. This also suggests that the effects of varying $R_{\rm Low}$ to
the $230$ GHz image are different between negative and non-negative spinning
black holes, and we will compare and explain the differences in more detail in
our next papers.

In the contour sub-plots, we find an almost log-linear relationship between
$R_{\rm Low}$ and $M_{\rm Unit}$ for suites of Sgr A* models for all of the
$\theta$ and $R_{\rm High}$ we considered. This could be interpreted
qualitatively: as $R_{\rm Low}$ increases, the thermal energy partitioned to
the electrons at low $\beta$ becomes smaller, making the electrons cooler and
emitting less intense synchrotron radiation. For the total time-averaged flux
to remain at $2.4\,\textrm{Jy}$, the optical depth shall increase to make the
image larger and thus increase its overall luminosity. We note a similar
pattern between $R_{\rm Low}$ and $M_{\rm Unit}$ for other black-hole spins.

###  3.2 Parameter Dependence

We want to investigate how $M_{\Delta T}$ vary as $R_{\rm Low}$ for suites of
Sgr A* models. We extract the values of $M_{\Delta T}$ and $R_{\rm Low}$ along
the black-dotted lines in Figure 2 and present them in Figure 3. To extract
the values of $M_{\Delta T}$, we perform Bivariate spline approximations with
default parameters provided by the package Scipy, and we did it for all the
spins we considered. In sum, we find three different patterns for the
$M_{\Delta T}$ versus $R_{\rm Low}$ curves: _i_) $M_{\Delta T}$ increases as
$R_{\rm Low}$, _ii_) $M_{\Delta T}$ decreases as $R_{\rm Low}$, and _iii_)
$M_{\Delta T}$ first decrease to a local minimum and then increase. We observe
that except for the case of $R_{\rm Low}=10$, the shape of the $M_{\Delta T}$
versus $R_{\rm Low}$ curves is less sensitive to the variations of $R_{\rm
Low}$, but is more sensitive to the changes in $\theta$.

Almost every non-negative spin model, with $R_{\rm High}>10$, shows pattern
_iii_ as $R_{\rm Low}$ fist increases beyond $1$. Whether or not $M_{\Delta
T}$ reduces for negative spin models depends on $a$ and $\theta$. Also, the
$M_{\Delta T}$ versus $R_{\rm Low}$ curves for negative spins differ from
those for positive spins. For instance, none of the negative spin curves shows
strong dips in the $M_{\Delta T}$ versus $R_{\rm Low}$ relations, and they
vary more slowly and smoothly than those of the positive spins. This also
proves that the GRRT images of varying $R_{\rm Low}$ for the negative-spin
black holes fundamentally differ from those for the positive spins.

We also find that some suites of Sgr A* models with $R_{\rm Low}>10$ show a
reduction in $M_{\Delta T}$ as $R_{\rm Low}$ increases. For instance, it is
possible to achieve $M_{\Delta T}\lesssim 0.1$ for almost every model with
spin $a=+0.94$ and $a=0$. This result contrasts those presented in Event
Horizon Telescope Collaboration et al. (2022) and is discovered by considering
a more comprehensive range of parameter spaces. Thus, we successfully show
that it is possible to reduce $M_{\Delta T}$ by varying $R_{\rm Low}$. We
note, however, that $M_{\Delta T}$ remains larger than $0.1$ for almost every
model for the remaining spins. Thus, whether or not a substantial reduction in
$M_{\Delta T}$ occurs when $R_{\rm Low}$ is varied along suites of Sgr A*
models highly depends on the underlying model parameters.

###  3.3 Large $M_{\Delta T}$ Uncovered - Image Domain

Next, we would like to understand deeper about the origin of large $M_{\Delta
T}$ for the GRRT images of Sgr A* at $R_{\rm Low}=1$, the value assumed in
Event Horizon Telescope Collaboration et al. (2022). To achieve this, we
computed the time domain, pixel-wise _i_) averaged and _ii_) standard
deviation of image luminosity across all $a$, and $\theta$. We show the time-
averaged, pixel-wise luminosities in Figure 4 for the $a=+0.94$ GRMHD models.
Results for the remaining $a$ are shown in Figure Set 2. The disk is optically
thin. Hence, for $\theta=0.1,179.9\,\mathrm{Deg}$, the GRRT images show a thin
photon ring. As $\theta$ progresses towards $90$ Deg, most of the emission is
concentrated to the left of the black-hole shadow due to the Doppler beaming
effect. We find that for all the $R_{\rm Low}=1$ images, the emissions are
contributed from a thin photon ring. This is true even for other spins $a$.

Having all emissions contributed from a small area indicates that a slight
variation in the photon ring luminosity would induce large variability in the
$230$ GHz flux. Indeed, we show the time-domain, pixel-wise standard deviation
of the $230$ GHz luminosity for the $a=+0.94$ GRMHD models in Figure 5.
Results for the remaining $a$ are shown in Figure Set 3. We find that the
photon rings contribute most to the image variability. However, some minor
contributions from the accretion flow could also be seen. This is true also
for the $a=+0.5$ and $0$ models. On the other hand, the situation for negative
spins is different from that of positive ones. In addition to the photon-ring
variability, the contribution to the luminosity variability from the accretion
flow is larger for the negative spin models, and their contribution is the
largest for the $a=-0.94$ model. This could be due to the misalignment between
the black-hole and torus angular momentum, creating highly irregular accretion
flow dynamics.

We note that the image asymmetry due to the Doppler beaming effect is reduced
as the black-hole spin turns more negative. Such an effect is first mentioned
in Medeiros et al. (2022), where they consider image asymmetry up until $a=0$.
Here, we briefly extend their work to consider also negative spins. The image
becomes more symmetric for $a=-0.5$ and $a=-0.94$. This is probably due to the
retrograde motion of the negatively spinning black hole, which drags the
rotating plasma in the opposite direction to the plasma angular velocity,
reducing the angular velocity differences between the plasma on the left and
right-hand side of the event horizon and thus weakening the Doppler beaming
effect.

###  3.4 $M_{\Delta T}$ at Larger $R_{\rm High}$?

Figure 6: The variation of the minimum of $M_{\Delta T}$, obtained through
varying $R_{\rm Low}$, against a given $R_{\rm High}$. They are shown as
scatter points. This is for the particular black-hole parameter with $a=+0.94$
and $\theta=0.1\,\mathrm{Deg}$, which shows exceptionally small $M_{\Delta
T}$. In each subplot, the upper right box shows the telescope data as well as
the band filter that we use to compute the $3\,\mathrm{hours}$ variability.
The solid red horizontal lines and the red shaded represent the observed mean
and standard deviation $M_{\Delta T}$, respectively. We use Bivariate spline
approximations to obtain the blue scatter points while using direct GRRT to
compute the orange scatter points.

From Figure 3, we find that the minimum $M_{\Delta T}$ across different
$R_{\rm Low}$ for spin $a=+0.94$ and $\theta=0.1\,\mathrm{Deg}$ decreases as
$R_{\rm High}$, with the minimum of $M_{\Delta T}<0.05$ for $R_{\rm
High}=180$. This motivates us to perform extra parameter searches for the
model parameters $a=+0.94$ and $\theta=0.1\,\mathrm{Deg}$ with $R_{\rm High}$
beyond $180$ and extract the corresponding minimum $R_{\rm Low}$. We aim to
find whether $M_{\Delta T}$ could be further lowered and if they are
consistent with the observed $3\,\mathrm{hours}$ variability of Sgr A*.

In particular, we consider $R_{\rm High}=200,220$, and $240$. We performed
GRRT with varying $R_{\rm Low}$ given these $R_{\rm High}$ and to look for
models that satisfy the $2.4$ Jy constraints. We then extract the
corresponding minimum of $M_{\Delta T}$. We use two methods to achieve this.
The first is to directly look for the minimum of $M_{\Delta T}$ along the
curves generated by the Bivariate spline approximations (c.f. Figure 3). The
second is to look for the value of $R_{\rm Low}$ (and $M_{\rm Unit}$) that
produces such a $M_{\Delta T}$, and we further perform GRRT to compute the
corresponding $M_{\Delta T}$. In addition, we extract ALMA and SMA
observations of Sgr A* at $230$ GHz and compute the mean and standard
deviation of the observed $M_{\Delta T}$. To do this, we perform windowed
sampling with a window size of $3\,\mathrm{hours}$ and a sampling interval of
$0.5\,\mathrm{hours}$. We also note that the ALMA and SMA observations are
performed with the High and Low bands. We compute the $M_{\Delta T}$
separately.

These results are shown in Figure 6. We find that all $R_{\rm High}=10$ models
are inconsistent with the $M_{\Delta T}$ constraint even if we vary $R_{\rm
Low}$. $R_{\rm High}\geq 50$ models could attain $M_{\Delta T}$ consistent
with the observations. As $R_{\rm High}$ increases beyond $180$, the minimum
of $M_{\Delta T}$ could be as low as $\sim 0.03$. Also, we find that the
minimum of $M_{\Delta T}$ decreases for a larger $R_{\rm High}$, but the rate
of reduction also reduces. Indeed, the minimum of $M_{\Delta T}$ for $R_{\rm
High}=240$ is close to that of $R_{\rm High}=220$. This might suggest a
threshold value of $M_{\Delta T}$, where $M_{\Delta T}$ could not be lowered
than such value even if $R_{\rm Low}$ is considered. More comprehensive
statistical studies will help reveal whether such value exists.

##  4 Conclusion

In this paper, we performed a comprehensive set of GRRT parameter surveys
using the Open Science Grid to study the reason for the unexpectedly high
modulation index $M_{\Delta T}$ of Sgr A* predicted by theoretical models in
Event Horizon Telescope Collaboration et al. (2022). Our parameter surveys
vary the $R_{\rm Low}$ parameter in the electron temperature prescription
function used by Event Horizon Telescope Collaboration et al. (2022), which
spans the optically thin to the optically thick regimes, and covers $R_{\rm
Low}$ from $1$ to $60$. Our results contrast those in Event Horizon Telescope
Collaboration et al. (2022). We find that it is possible to reduce $M_{\Delta
T}$ to a value lower than $0.1$, the threshold values of the historical
distribution of $M_{\Delta T}$ for Sgr A*, by varying $R_{\rm Low}$, for
certain sets of black-hole spin $a$, inclination angle $\theta$, and $R_{\rm
High}$.

We find that models with $a=+0.94$, $\theta=0.1\,\mathrm{Deg}$, and $R_{\rm
High}\geq 50$ could produce $M_{\Delta T}$ consistent with the observed
$M_{\Delta T}$ of Sgr A*. We also note that the minimum of $M_{\Delta T}$
reduces for larger $R_{\rm High}$ for this particular set of $a$ and $\theta$,
and $M_{\Delta T}\sim 0.03$ for $R_{\rm High}\geq 180$. Although this might be
inconsistent with the statistical average of $M_{\Delta T}$, we note that the
minimum of $M_{\Delta T}$ across the historical distribution could be as low
as $0.03$. Thus, such a result sheds light on whether $R_{\rm Low}$ could be
an important parameter to consider if one intends to perform GRRT on suites of
Sgr A* simulation libraries.

We also analyze the GRRT images for the $R_{\rm Low}=1$ models, which are also
the values assumed in Event Horizon Telescope Collaboration et al. (2022),
across different $a$ and $\theta$. By computing the time domain, pixel-wise,
averaged, and standard deviation of the image luminosity, we find that all
$R_{\rm Low}=1$ models have their luminosity contributions coming from a thin
photon ring. We also find that the thin photon ring contributes most to the
observed variability for all $a\geq 0$ models. For negative spin models,
secondary contributions of variabilities from the accretion flow are more
visible and are due to the irregular flow dynamics induced by misalignment
between the black-hole and torus angular momentum. Our findings on the origin
of high $M_{\Delta T}$ suggest that a more serious study on the electron
temperature used to model RIAF shall be adopted. These results, together with
the findings that the behaviour of $M_{\Delta T}$ versus $R_{\rm Low}$
relations are different between the negative and non-negative spinning black
holes, also lay the foundations for our successive papers, where we would
analyze why a variation of $R_{\rm Low}$ could reduce/increase $M_{\Delta T}$
for different black-hole parameters.

Note that there are still several caveats in our preliminary study. We assumed
only $R_{\rm Low}\leq 60$ in our study. This might raise concerns about the
lack of study on the behaviour of $M_{\Delta T}$ against $R_{\rm Low}$ at a
larger $R_{\rm Low}$. However, we have already shown that $M_{\Delta T}$ would
decrease in the range of parameter spaces we are interested in. Thus, there
seems to be no immediate need to consider higher $R_{\rm Low}$. Also, we show
an almost log-linear relationship between $R_{\rm Low}$ and $M_{\rm Unit}$ for
suites of Sgr A* models. Thus, a further increase in $R_{\rm Low}$ would make
the accretion disk more optically thick. This would violate our assumption of
Sgr A* being a low-luminosity, geometrically thick, and optically thin
accretion flow. Also, a more optically thick image would have a larger size
and might be inconsistent with the image-size constraint given in Event
Horizon Telescope Collaboration et al. (2022).

Computing the effects of $R_{\rm Low}$ on $M_{\Delta T}$ for one given time
interval $\tau=(\textrm{29,465}$ – $\textrm{29,995})\,GMc^{-2}$ only might
induce concerns about whether our results are statistically significant. Also,
coarse parameter spaces that we considered in our work might pose concerns
about whether the minimum of $M_{\Delta T}$ is due to interpolation errors.
However, we show in Figure 6 that the minimum of $M_{\Delta T}$ computed by
the Bivariate spline approximations and direct GRRT are close to each other.
Nonetheless, these limitations are due to the limited computational storage
and time constraints. We plan to perform more comprehensive studies in the
future to better understand, in a more serious statistical framework, the role
of $R_{\rm Low}$ on the increase/decrease of $M_{\Delta T}$.

Ho Sang (Leon) Chan acknowledges support from the Croucher Scholarship for
Doctoral Studies by the Croucher Foundation. C.-k. C. acknowledge NSF grant
OISE 17-43747 support.

## References

  * Akiyama et al. (2019) Akiyama, K., Alberdi, A., Alef, W., et al. 2019, The Astrophysical Journal Letters, 875, L5
  * Akiyama et al. (2022) —. 2022, The Astrophysical Journal Letters, 930, L14
  * Baganoff et al. (2003) Baganoff, F. K., Maeda, Y., Morris, M., et al. 2003, The Astrophysical Journal, 591, 891
  * Begelman (2014) Begelman, M. C. 2014, arXiv preprint arXiv:1410.8132
  * Chan et al. (2015) Chan, C.-K., Psaltis, D., Özel, F., et al. 2015, The Astrophysical Journal, 799, 1
  * Dexter (2016) Dexter, J. 2016, Monthly Notices of the Royal Astronomical Society, 462, 115
  * Dexter et al. (2009) Dexter, J., Agol, E., & Fragile, P. C. 2009, ApJ, 703, L142, doi: 10.1088/0004-637X/703/2/L142
  * Dexter et al. (2020) Dexter, J., Jiménez-Rosales, A., Ressler, S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 494, 4168
  * Eckart et al. (2002) Eckart, A., Genzel, R., Ott, T., & Schödel, R. 2002, Monthly Notices of the Royal Astronomical Society, 331, 917
  * Event Horizon Telescope Collaboration et al. (2022) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022, ApJ, 930, L16, doi: 10.3847/2041-8213/ac6672
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565
  * Gottlieb et al. (2011) Gottlieb, S., Ketcheson, D. I., & Shu, C.-W. 2011, Strong stability preserving Runge-Kutta and multistep time discretizations (World Scientific)
  * Jiang & Shu (1996) Jiang, G.-S., & Shu, C.-W. 1996, Journal of computational physics, 126, 202
  * Lo et al. (1993) Lo, K., Backer, D., Kellermann, K., et al. 1993, Nature, 362, 38
  * Medeiros et al. (2022) Medeiros, L., Chan, C.-K., Narayan, R., Özel, F., & Psaltis, D. 2022, The Astrophysical Journal, 924, 46
  * Mościbrodzka et al. (2016) Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, Astronomy & Astrophysics, 586, A38
  * Mościbrodzka & Gammie (2018) Mościbrodzka, M., & Gammie, C. F. 2018, Monthly Notices of the Royal Astronomical Society, 475, 43
  * Mościbrodzka et al. (2009) Mościbrodzka, M., Gammie, C. F., Dolence, J. C., Shiokawa, H., & Leung, P. K. 2009, ApJ, 706, 497, doi: 10.1088/0004-637X/706/1/497
  * OSG (2006) OSG. 2006, OSPool, OSG, doi: 10.21231/906P-4D78
  * Özel et al. (2000) Özel, F., Psaltis, D., & Narayan, R. 2000, ApJ, 541, 234, doi: 10.1086/309396
  * Pordes et al. (2007) Pordes, R., Petravick, D., Kramer, B., et al. 2007, in 78, Vol. 78, J. Phys. Conf. Ser., 012057, doi: 10.1088/1742-6596/78/1/012057
  * Prather (2022) Prather, B. 2022, PhD thesis
  * Prather et al. (2021) Prather, B. S., Wong, G. N., Dhruv, V., et al. 2021, arXiv preprint arXiv:2110.10191
  * Quataert & Gruzinov (2000) Quataert, E., & Gruzinov, A. 2000, The Astrophysical Journal, 545, 842
  * Quataert et al. (1999) Quataert, E., Narayan, R., & Reid, M. J. 1999, The Astrophysical Journal, 517, L101
  * Scepi et al. (2022) Scepi, N., Dexter, J., & Begelman, M. C. 2022, Monthly Notices of the Royal Astronomical Society, 511, 3536
  * Sfiligoi et al. (2009) Sfiligoi, I., Bradley, D. C., Holzman, B., et al. 2009, in 2, Vol. 2, 2009 WRI World Congress on Computer Science and Information Engineering, 428–432, doi: 10.1109/CSIE.2009.950
  * Shapiro et al. (1976) Shapiro, S. L., Lightman, A. P., & Eardley, D. M. 1976, ApJ, 204, 187, doi: 10.1086/154162
  * Sironi & Spitkovsky (2011) Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75, doi: 10.1088/0004-637X/726/2/75
  * Tóth (2000) Tóth, G. 2000, Journal of Computational Physics, 161, 605
  * Werner et al. (2018) Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, Monthly Notices of the Royal Astronomical Society, 473, 4840
  * Wielgus et al. (2022) Wielgus, M., Marchili, N., Martí-Vidal, I., et al. 2022, ApJ, 930, L19, doi: 10.3847/2041-8213/ac6428
  * Wong et al. (2022) Wong, G. N., Prather, B. S., Dhruv, V., et al. 2022, The Astrophysical Journal Supplement Series, 259, 64
  * Zhdankin et al. (2019) Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, Physical review letters, 122, 055101

##  Appendix A GRRT parameter surveys

In this appendix, we present the GRRT parameter survey results for $a=+0.5$
(Figure 7), $a=0$ (Figure 8), $a=-0.5$ (Figure 9), and $a=-0.94$ (Figure 10).
These figures are all included in Figure Set 1 and are presented here for
peer-reviewing convenience. They will be removed upon publication.

Figure 7: $M_{\Delta T}$ against $R_{\rm Low}$ and $M_{\rm Unit}$ for
$a=+0.5$. See the caption of Figure 2 for details.  Figure 8: $M_{\Delta T}$
against $R_{\rm Low}$ and $M_{\rm Unit}$ for $a=0$. See the caption of Figure
2 for details.  Figure 9: $M_{\Delta T}$ against $R_{\rm Low}$ and $M_{\rm
Unit}$ for $a=-0.5$. See the caption of Figure 2 for details.  Figure 10:
$M_{\Delta T}$ against $R_{\rm Low}$ and $M_{\rm Unit}$ for $a=-0.94$. See the
caption of Figure 2 for details.

##  Appendix B GRRT Images for $R_{\rm Low}=1$

This appendix presents the time-domain, pixel-wise, average and standard
deviation of the GRRT images. This include results for $a=+0.5$ (Figure 11 and
), $a=0$ (Figure 13 and ), $a=-0.5$ (Figure 15 and ), and $a=-0.94$ (Figure 17
and ). These figures are all included in Figure Sets 2 and 3 and are presented
here for peer-reviewing convenience. They will be removed upon publication.

Figure 11: Time-averaged, pixel-wise luminosity across different $R_{\rm
High}$ and $\theta$ (with $R_{\rm Low}=1$) for $a=+0.5$. See the caption of
Figure 4 for details.  Figure 12: Pixel-wise, time-domain standard deviation
of the luminosity across different $R_{\rm High}$ and $\theta$ (with $R_{\rm
Low}=1$) for $a=+0.5$. See the caption of Figure 5 for details.  Figure 13:
Time-averaged, pixel-wise luminosity across different $R_{\rm High}$ and
$\theta$ (with $R_{\rm Low}=1$) for $a=0$. See the caption of Figure 4 for
details.  Figure 14: Pixel-wise, time-domain standard deviation of the
luminosity across different $R_{\rm High}$ and $\theta$ (with $R_{\rm Low}=1$)
for $a=0$. See the caption of Figure 5 for details.  Figure 15: Time-averaged,
pixel-wise luminosity across different $R_{\rm High}$ and $\theta$ (with
$R_{\rm Low}=1$) for $a=-0.5$. See the caption of Figure 4 for details.
Figure 16: Pixel-wise, time-domain standard deviation of the luminosity across
different $R_{\rm High}$ and $\theta$ (with $R_{\rm Low}=1$) for $a=-0.5$. See
the caption of Figure 5 for details.  Figure 17: Time-averaged, pixel-wise
luminosity across different $R_{\rm High}$ and $\theta$ (with $R_{\rm Low}=1$)
for $a=-0.94$. See the caption of Figure 4 for details.  Figure 18: Pixel-
wise, time-domain standard deviation of the luminosity across different
$R_{\rm High}$ and $\theta$ (with $R_{\rm Low}=1$) for $a=-0.94$. See the
caption of Figure 5 for details.


# Sailing the CE$\nu$NS Seas of Non-Standard Neutrino Interactions with the
Coherent CAPTAIN Mills Experiment

Ian M. Shoemaker<EMAIL_ADDRESS>Center for Neutrino Physics, Department of
Physics, Virginia Tech University, Blacksburg, VA 24601, USA  Eli Welch
Center for Neutrino Physics, Department of Physics, Virginia Tech University,
Blacksburg, VA 24601, USA

###### Abstract

We study future coherent elastic neutrino-nucleus scattering (CE$\nu$NS)
modifications from a variety of possible models at the Coherent CAPTAIN Mills
(CCM) experiment at Los Alamos. We show that large regions of Non-Standard
Neutrino Interaction (NSI) parameter space will be excluded rapidly, and that
stringent new bounds on the gauge coupling in $Z^{\prime}$ models will also be
placed. As a result, CCM will be able to rule out LMA-D solutions for a large
class of models with MeV-scale mediators.

##  I Introduction

In the aftermath of the observation of neutral currents by Gargamelle in 1973
Hasert  _et al._ (1973), it was quickly realized that the neutral currents
will induce a large coherent cross section on nuclei at low energies Freedman
(1974). In 2017, this Coherent Elastic Neutrino-Nucleus Scattering (CE$\nu$NS)
was observed by the COHERENT collaboration, yielding qualitatively new
information on neutrinos Akimov  _et al._ (2017). This is the largest cross
section at low-energies, and its utility allows for the study of precision
neutrino-nucleus interactions Scholberg (2006); Coloma  _et al._ (2017);
Shoemaker (2017); Liao and Marfatia (2017); Denton  _et al._ (2018), sterile
neutrinos Anderson  _et al._ (2012); Dutta  _et al._ (2016); Kosmas  _et al._
(2017), as well as non-neutrino BSM such as light DM deNiverville  _et al._
(2015); Ge and Shoemaker (2018); Dutta  _et al._ (2020a); Hurtado  _et al._
(2020); Dutta  _et al._ (2020b) and ALPs Dent  _et al._ (2020); Aristizabal
Sierra  _et al._ (2020).

Looking to the future, the European Spallation Source (ESS) Baxter  _et al._
(2020), the Coherent CAPTAIN-Mills (CCM) at Lujan center at Los Alamos
National Laboratory Dunton (2020), and Oak Ridge National Laboratory Akimov
_et al._ (2020) will have ton and multi-ton detectors deployed. This rapid
expansion of statistics with complementarity detector types will allow global
CE$\nu$NS data to place powerful new constraints on new physics. Recent work
has highlighted the need for a careful statistical treatment Denton and
Gehrlein (2020), the usage of copula in order to combine global data in an
effort to break NSI degeneracies Dutta  _et al._ (2020c), and combining timing
and energy information to maximize the search sensitivity to new physics
Coloma  _et al._ (2020); Dutta  _et al._ (2020a). There may be a hint of
nonzero NSI in the tension between T2K and NOvA Denton  _et al._ (2021).
Experimentally, much attention has been paid to the precise value of the
quenching at the lowest energies for CsI[Na] Dutta  _et al._ (2020a); Collar
_et al._ (2019). This can play an important role in determining the precise
sensitivity to BSM, especially for example to low-mass mediators of neutrino-
nucleus scattering.

In this work, we focus on the BSM capabilities of the CCM experiment at Los
Alamos. In particular, we will study its ability to probe NSI and light
mediator $Z^{\prime}$ models. These scenarios are naturally related to each
other given that a sufficiently heavy $Z^{\prime}$ can be integrated out
leaving behind the neutral current four-fermion effective operators of NSI.

For many decades, theorists have considered the existence of new neutral
current interactions between neutrinos and matter Wolfenstein (1978). Dubbed
Non-Standard neutrino Interactions (NSIs), this new physics can alter neutrino
production, detection, and their propagation (see the reviews Ref. Davidson
_et al._ (2003); Biggio  _et al._ (2009); Farzan and Tortola (2018); Dev
(2019)). In this work, we will focus on neutral current NSI described by the
dimension-6 effective Lagrangian,

$\mathscr{L}_{NSI}\supset
2\sqrt{2}G_{F}\epsilon_{\alpha\beta}\left(\bar{\nu}^{\alpha}\gamma^{\mu}\nu_{\beta}\right)\left(\bar{f}\gamma^{\mu}f\right),$
(1)

where $\alpha,\beta$ are neutrino flavor indices and $f$, is a charged SM
fermion, $f=\\{e,u,d\\}$. The strength of these new operators is parameterized
by the $\epsilon$ coefficients in units of $G_{F}$. In the language of EFT, we
regard the above four-fermion description as an effective description below
the energy scale $\Lambda\equiv 1/\sqrt{2\sqrt{2}\epsilon G_{F}}$.
Historically NSI studies have focused on TeV (or above) scales for the source
of new physics giving rise to Eq. (1). In this case, large $\epsilon$’s come
about from order unity couplings with EW scale masses. Given that the impact
of such non-renormalizable operators grows rapidly with energy, colliders have
placed strong bounds on NSI Berezhiani and Rossi (2002); Davidson and Sanz
(2011); Friedland  _et al._ (2012); Buarque Franzosi  _et al._ (2016); Babu
_et al._ (2020).

Some of the most model-independent bounds on NSI come from oscillation data.
In this case, NSI modifies the forward coherent scattering of neutrinos as
they travel through matter Wolfenstein (1978). As such, the presence of NSI
impacts the flavor evolution of neutrinos in matter, resulting in tight
constraints on NSI Farzan and Tortola (2018); Esteban  _et al._ (2018), and
their presence can significantly complicate the interpretation of oscillation
data (see e.g. Friedland  _et al._ (2004); Friedland and Shoemaker (2012);
Capozzi  _et al._ (2020)). Bounds of this type apply to most models so long as
the mediator is sufficiently heavy ($\gtrsim 10^{-12}~{}{\rm eV}$) Gonzalez-
Garcia  _et al._ (2007); Coloma  _et al._ (2021).

CE$\nu$NS data occupies an intermediate regime of momentum transfer where the
model dependence of realistic UV completions of NSI can be probed Shoemaker
(2017); Liao and Marfatia (2017); Denton  _et al._ (2018). Although CE$\nu$NS
data therefore probes NSI in a less model-independent way than oscillations,
this dependence on the underlying physics allows for potential sensitivity to
directly measuring the new physics scale associated with NSI Shoemaker (2017).

The remainder of this paper is organized as follows. In Sec. II we review the
capabilities of CCM and the assumptions we make regarding its performance. In
Sec. III we outline our calculational framework, and the assumed BSM
modifications to CE$\nu$NS. We find that with CCM alone, diagonal NSI
constraints will be greatly improved, although some degeneracies will remain.
We also find that simple $Z^{\prime}$ models will be ruled out as possible
LMA-D solutions Miranda  _et al._ (2006), highlighting the important
complementarity between CE$\nu$NS and oscillation data. Finally in Sec. IV we
discuss future directions and conclude in Sec. V.

##  II CCM Capabilities

A significant advance in CE$\nu$NS data will happen soon, thanks to the
Coherent CAPTAIN (Cryogenic Apparatus for Precision Tests of Argon
Interactions with Neutrinos) Mills (CCM) experiment. The experiment will take
place at Los Alamos National Lab’s Lujan Center at the Los Alamos Neutron
Science Center (LANSCE). A very high instantaneous power is crucial for
measuring the signal-to-background ratio, and makes the Lujan Center a
competitive environment for neutrino physics. While a 300 ns beam pulse width
can reduce background by a factor $\sim 10^{5}$, they expect to reach 30 ns
widths with an order of magnitude gain in background rejection.

Figure 1: Illustrative predicted event spectra at CCM for heavy mediator NSI
(left) and light mediator $Z^{\prime}$ (right).

Following Dunton (2020) we assume that the CCM will employ a fiducial 7 ton
liquid argon (LAr) detector with 3 tons of LAr veto. The active veto will be
employed to reject cosmic ray and other relevant backgrounds. Assuming a 25
keV threshold, 20 meters from the target, the neutrino flux at the detector
will be $\sim 4.5\times 10^{5}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}$ Dunton (2020).
In what follows we will assume 3 years of running.

##  III CCM Analysis Setup

We will assume that the neutrino flux is well-described by pion decay at-rest,
$\pi^{+}\rightarrow\mu^{+}+\bar{\nu}_{\mu}$, followed by
$\mu^{+}\rightarrow\nu_{\mu}+e^{-}+\bar{\nu}_{e}$. This induces a sharply
peaked muon anti-neutrino flux, resulting from two-body decay kinematics. In
contrast, the electron anti-neutrino and muon neutrino arising from muon decay
are distributed over a range of energies up to a cut-off around 50 MeV. As
such we model the fluxes as:

$\phi_{\nu_{\mu}}(E_{\nu})=\phi_{0}~{}\delta(E_{\nu}-E_{\nu_{0}})$ (2)

where $E_{\nu_{0}}\equiv\frac{m^{2}_{\pi}-m^{2}_{\mu}}{2m_{\pi}}\approx 29.8$
MeV. And

$\phi_{\bar{\nu}_{\mu}}(E_{\nu})=\phi_{0}~{}\frac{64E^{2}_{\nu}}{m_{\mu}^{3}}~{}\left(\frac{3}{4}-\frac{E_{\nu}}{m_{\mu}}\right),$
(3)
$\phi_{\nu_{e}}(E_{\nu})=\phi_{0}~{}\frac{192E^{2}_{\nu}}{m_{\mu}^{3}}~{}\left(\frac{1}{2}-\frac{E_{\nu}}{m_{\mu}}\right)$
(4)

where the $\phi_{0}$ normalization is fixed by the number of proton on target
and the distance to the detector.

Next we write the CE$\nu$NS differential cross section on a nucleus of mass
$M$ as

$\frac{d\sigma_{\alpha}}{dE_{r}}=\frac{G^{2}_{F}}{2\pi}Q^{2}_{\alpha}F^{2}(2ME_{r})M\left(2-\frac{ME_{r}}{E^{2}_{\nu}}\right)$
(5)

where $E_{r}$ is the nuclear recoil energy, and $F(Q^{2})$ is the nuclear form
factor, which we parameterize as in Engel (1991). The effective charge
$Q_{\alpha}$ in the SM is given by

$Q^{2}_{\alpha,SM}=(Zg^{V}_{p}+Ng^{V}_{n})^{2},$ (6)

where $N$ and $Z$ are the number of neutrons and protons respectively, and the
SM $Z^{0}$ couplings to proton and neutron are
$g_{p}^{V}=\frac{1}{2}-2\sin^{2}\theta_{w},$ and $g_{n}^{V}=-\frac{1}{2}$.

More generally, we will consider two BSM scenarios which modify the effective
charge $Q_{\alpha}$. First we will examine conventional NSI in the heavy
mediator limit, as in Eq. 1. Thus in the heavy mediator limit the effective
charge entering into the cross section (Eq. 5) becomes.

$Q^{2}_{\alpha}=[Z(g_{p}^{V}+2\epsilon_{\alpha\alpha}^{uV}+\epsilon_{\alpha\alpha}^{dV})+N(g_{n}^{V}+\epsilon_{\alpha\alpha}^{uV}+2\epsilon_{\alpha\alpha}^{dV})]^{2}+\sum_{\beta\neq\alpha}[Z(2\epsilon_{\alpha\beta}^{uV}+\epsilon_{\alpha\beta}^{dV})+N(\epsilon_{\alpha\beta}^{uV}+2\epsilon_{\alpha\beta}^{dV})]^{2},$
(7)

where the first squared term includes those NSI coefficients which interfere
with the SM contribution, and the second squared terms include all the off-
diagonal $\epsilon_{\alpha\beta}$ which have no SM analog. For simplicity in
this paper we will assume flavor-diagonal NSI, leaving flavor off-diagonal
terms for future work.

Figure 2: 90$\%$ CL sensitivity projections on NSI parameters for up-quark
couplings from the CCM experiment. The constraints on down-type couplings are
qualitatively similar. The right panel shows a close-up of the allow region in
the immediate vicinity of the SM.

On the other hand, new physics contributions to CE$\nu$NS can be hidden from
us at low and accessible mass scales if it is very feebly coupled Farzan
(2015); Farzan and Shoemaker (2016); Farzan and Heeck (2016). Consider for
example the following phenomenological $Z^{\prime}$ model giving rise to NSI

$\mathscr{L}\supset
g_{q}Z^{\prime}_{\mu}\bar{q}\gamma^{\mu}q+g_{\nu}Z^{\prime}_{\mu}\bar{\nu}\gamma^{\mu}\nu$
(8)

this yields Eq.(1) when the energies of interest are below $\Lambda\equiv
m_{Z^{\prime}}/\sqrt{g_{\nu}g_{q}}$. For simplicity here we assume that the
$Z^{\prime}$ couples equally to all quark flavors. We will however examine the
coupling to electron- and muon-flavor neutrinos separately. With these
assumptions the effective charge entering into Eq. 5 is

$Q^{2}_{\alpha}=[Z(g^{V}_{p}+\frac{3g^{2}}{2\sqrt{2}G_{F}(2ME_{r}+M^{2}_{Z^{\prime}})})+N(g^{V}_{n}+\frac{3g^{2}}{2\sqrt{2}G_{F}(2ME_{r}+M^{2}_{Z^{\prime}})})]^{2}$
(9)

Finally we compute the differential rate of events per unit detector mass as

$\frac{dR_{\alpha}}{dE_{R}}=\int\phi_{\alpha}\frac{d\sigma_{\alpha}}{dE_{R}}dE_{\nu},$
(10)

with $\phi_{\alpha}$ the fluxes given in Eqs. 2, 3, and 4, and the cross
section $d\sigma_{\alpha}/dE_{R}$ given in Eq. 5.

To be conservative we assume that the number of background events will be
20$\%$ of the signal events. Flux normalization uncertainties and detector
efficiencies and calibrations will be important systematic uncertainties. As a
rough treatment, we take these to be a 10$\%$ overall uncertainty on the flux
normalization. These assumptions are conservatively based on COHERENT Akimov
_et al._ (2018), though it is expected that CCM will benefit from an enhanced
background rejection capability at the Lujan facility.

We display the nuclear recoil spectra for both the NSI and $Z^{\prime}$ models
in Fig. 1. As can be seen, heavy mediator NSI produces a spectrum broadly
similar to the SM, though differing of course in normalization. This is
expected given that no novel energy dependence arises in heavy mediator NSI
(see Eq. 7). The $Z^{\prime}$ model on the other hand, includes the explicit
energy dependence in Eq. 9 arising from the $Z^{\prime}$ propagator. As a
result the spectrum appearing in the right panel of Fig. 1 differs in shape
compared to the SM expectation. As previously pointed out, such a feature can
potentially be used to determine the $Z^{\prime}$ mass in the event of a
discovery Shoemaker (2017).

Next we compute the sensitivity CCM can expect to both heavy mediator NSI and
light $Z^{\prime}$ models. In Fig. 2 we examine the sensitivity to the
diagonal elements $\epsilon_{ee}$ and $\epsilon_{\mu\mu}$. Unlike the off-
diagonal NSI coefficients, these benefit from interference with the SM
CE$\nu$NS contributions and are as a result more tightly constrained. We see
that there are four allowed regions, one of which (bottom left) which
encompass the SM ($\epsilon=0$). We note that although the regions we find are
quite different than those found for COHERENT Liao and Marfatia (2017), they
are similar to those found for a projection study of COHERENT before the
experiment ran Coloma  _et al._ (2017). We note that the apparent NSI
sensitivity at CCM will be sufficiently strong to rule out the LMA-D region
(see e.g. Fig. 7 from Ref. Coloma  _et al._ (2017)), which is off the plot in
our Fig. 2.

Figure 3: 90$\%$ CL constraints on $Z^{\prime}$ completions of NSI for
electron neutrino couplings (left) and muon neutrino couplings (right). Note
the small allowed strip of parameter space in each case (see text for
discussion). For reference we also display the BBN Kamada and Yu (2015)
exclusions and LMA-D region.

Finally, we display the expected sensitivity to the $Z^{\prime}$ model in Fig.
3 with only the electron or muon coupling turned on in each panel. We have
also included the LMA-D region, and the bounds from BBN which are strong at
low mediator masses Kamada and Yu (2015). We have focused on $Z^{\prime}$
models with equal up and down-quark couplings, but note that it has recently
been argued that more general scenarios may remove some of the tension between
LMA-D and CE$\nu$NS data Esteban  _et al._ (2018); Coloma  _et al._ (2020);
Esteves Chaves and Schwetz (2021).

##  IV Discussion and Future Directions

In the future we plan to study the utility of multiple detectors with
different target nuclei to break degeneracies in NSI couplings which arise in
light mediator models. Although NSI coupling degeneracies are studied in many
scenarios, they seem not to have been studied in light mediator models which
effectively lead to momentum-dependent NSI couplings Shoemaker (2017). As
illustrated in Fig. 3, CE$\nu$NS plays an important in the ability to
constrain neutrino oscillations in matter (as in the LMA-D solution), and can
shed light on light $Z^{\prime}$ models that impact BBN Kamada and Yu (2015).
For example, CCM found evidence of $Z^{\prime}$ inside the BBN highlighted
region of Fig. 3, it would indicate that a non-trivial cosmological history
played a key role in ameliorating tension between the new light degree of
freedom at BBN.

Additionally, the power of CE$\nu$NS to probe models non-minimal sterile
neutrino models has not received much attention. For example, non-minimal
Sterile Neutrinos have been shown to be a viable explanation of the XENON1T
excess Bally  _et al._ (2020); Shoemaker  _et al._ (2020) and may have a
connection to the observed Hubble tension Berbig  _et al._ (2020). Similar
models can also be probed at CE$\nu$NS experiments via neutrino-to-sterile
upscattering. Moreover, in sterile neutrino models with enhanced baryon
couplings, the HNLs may be directly produced in the target Batell  _et al._
(2016), though this has not been studied in light of present CE$\nu$NS data.
We plan to return to these topics in future work.

##  V Conclusions

We have examined the near term sensitivity of the Coherent Captain Mills
experiment to a variety of BSM neutrino modifications. First we focused on the
sensitivity to neutral current contact operators, i.e. NSI. There we found
that the high-statistics which CCM will offer gives it a significant
opportunity to break strong NSI parameter degeneracies. We have focused on the
flavor diagonal NSI couplings, $\epsilon_{ee}$ and $\epsilon_{\mu\mu}$, but we
expect CCM to play an important role in constraining the off-diagonal elements
$\epsilon_{e\alpha}$ and $\epsilon_{\mu\alpha}$ as well. Then we examined a
simplified $Z^{\prime}$ model of neutrino-nucleus interactions. We found that
CCM will significantly improve over COHERENT’s bounds on the gauge. As such
this will allow it to close off a small window of parameter space which
previously allowed for LMA-D to be compatible with BBN Denton  _et al._
(2018). It is likely that a viable LMA-D solution in the scenario of Ref.
Farzan (2015) with MeV-scale mediators will be ruled out by CCM, though the
sub-eV scenario cannot be tested at CE$\nu$NS type experiments. Further,
viable LMA-D solutions based on purely leptonic gauge couplings offer an
alternative route for sufficiently suppressing CE$\nu$NS constraints Coloma
_et al._ (2021). Lastly, a difference in up- and down-quark couplings can
allow for destructive interference on particular nuclei, which can also
ameliorate the tension between CE$\nu$NS and LMA-D Esteves Chaves and Schwetz
(2021).

_Acknowledgements-_ We are very grateful for discussions with Bhaskar Dutta
and Adrian Thompson. This work is supported by the U.S. Department of Energy
under the award number DE-SC0020250.

## References

  * Hasert  _et al._ (1973) F. J. Hasert _et al._ (Gargamelle Neutrino), Phys. Lett. B 46, 138 (1973).
  * Freedman (1974) D. Z. Freedman, Phys. Rev. D 9, 1389 (1974).
  * Akimov  _et al._ (2017) D. Akimov _et al._ (COHERENT), Science 357, 1123 (2017), arXiv:1708.01294 [nucl-ex] .
  * Scholberg (2006) K. Scholberg, Phys. Rev. D 73, 033005 (2006), arXiv:hep-ex/0511042 .
  * Coloma  _et al._ (2017) P. Coloma, P. B. Denton, M. Gonzalez-Garcia, M. Maltoni, and T. Schwetz, JHEP 04, 116 (2017), arXiv:1701.04828 [hep-ph] .
  * Shoemaker (2017) I. M. Shoemaker, Phys. Rev. D 95, 115028 (2017), arXiv:1703.05774 [hep-ph] .
  * Liao and Marfatia (2017) J. Liao and D. Marfatia, Phys. Lett. B 775, 54 (2017), arXiv:1708.04255 [hep-ph] .
  * Denton  _et al._ (2018) P. B. Denton, Y. Farzan, and I. M. Shoemaker, JHEP 07, 037 (2018), arXiv:1804.03660 [hep-ph] .
  * Anderson  _et al._ (2012) A. Anderson, J. Conrad, E. Figueroa-Feliciano, C. Ignarra, G. Karagiorgi, K. Scholberg, M. Shaevitz, and J. Spitz, Phys. Rev. D 86, 013004 (2012), arXiv:1201.3805 [hep-ph] .
  * Dutta  _et al._ (2016) B. Dutta, Y. Gao, R. Mahapatra, N. Mirabolfathi, L. E. Strigari, and J. W. Walker, Phys. Rev. D 94, 093002 (2016), arXiv:1511.02834 [hep-ph] .
  * Kosmas  _et al._ (2017) T. Kosmas, D. Papoulias, M. Tortola, and J. Valle, Phys. Rev. D 96, 063013 (2017), arXiv:1703.00054 [hep-ph] .
  * deNiverville  _et al._ (2015) P. deNiverville, M. Pospelov, and A. Ritz, Phys. Rev. D 92, 095005 (2015), arXiv:1505.07805 [hep-ph] .
  * Ge and Shoemaker (2018) S.-F. Ge and I. M. Shoemaker, JHEP 11, 066 (2018), arXiv:1710.10889 [hep-ph] .
  * Dutta  _et al._ (2020a) B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, and L. E. Strigari, Phys. Rev. Lett. 124, 121802 (2020a), arXiv:1906.10745 [hep-ph] .
  * Hurtado  _et al._ (2020) N. Hurtado, H. Mir, I. M. Shoemaker, E. Welch, and J. Wyenberg, Phys. Rev. D 102, 015006 (2020), arXiv:2005.13384 [hep-ph] .
  * Dutta  _et al._ (2020b) B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, L. E. Strigari, and A. Thompson, (2020b), arXiv:2006.09386 [hep-ph] .
  * Dent  _et al._ (2020) J. B. Dent, B. Dutta, D. Kim, S. Liao, R. Mahapatra, K. Sinha, and A. Thompson, Phys. Rev. Lett. 124, 211804 (2020), arXiv:1912.05733 [hep-ph] .
  * Aristizabal Sierra  _et al._ (2020) D. Aristizabal Sierra, V. De Romeri, L. Flores, and D. Papoulias, (2020), arXiv:2010.15712 [hep-ph] .
  * Baxter  _et al._ (2020) D. Baxter _et al._ , JHEP 02, 123 (2020), arXiv:1911.00762 [physics.ins-det] .
  * Dunton (2020) E. Dunton, “Searching for sterile neutrinos at coherent captain mills.” https://indico.cern.ch/event/844613/contributions/3618342/attachments/1942075/3220598/Dunton_MagnificentCEvNS_2019.pdf (2020).
  * Akimov  _et al._ (2020) D. Akimov _et al._ (COHERENT), Phys. Rev. D 102, 052007 (2020), arXiv:1911.06422 [hep-ex] .
  * Denton and Gehrlein (2020) P. B. Denton and J. Gehrlein, (2020), arXiv:2008.06062 [hep-ph] .
  * Dutta  _et al._ (2020c) B. Dutta, R. F. Lang, S. Liao, S. Sinha, L. Strigari, and A. Thompson, JHEP 20, 106 (2020c), arXiv:2002.03066 [hep-ph] .
  * Coloma  _et al._ (2020) P. Coloma, I. Esteban, M. Gonzalez-Garcia, and M. Maltoni, JHEP 02, 023 (2020), [Addendum: JHEP 12, 071 (2020)], arXiv:1911.09109 [hep-ph] .
  * Denton  _et al._ (2021) P. B. Denton, J. Gehrlein, and R. Pestes, Phys. Rev. Lett. 126, 051801 (2021), arXiv:2008.01110 [hep-ph] .
  * Collar  _et al._ (2019) J. Collar, A. Kavner, and C. Lewis, Phys. Rev. D 100, 033003 (2019), arXiv:1907.04828 [nucl-ex] .
  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
  * Davidson  _et al._ (2003) S. Davidson, C. Pena-Garay, N. Rius, and A. Santamaria, JHEP 03, 011 (2003), arXiv:hep-ph/0302093 .
  * Biggio  _et al._ (2009) C. Biggio, M. Blennow, and E. Fernandez-Martinez, JHEP 08, 090 (2009), arXiv:0907.0097 [hep-ph] .
  * Farzan and Tortola (2018) Y. Farzan and M. Tortola, Front. in Phys. 6, 10 (2018), arXiv:1710.09360 [hep-ph] .
  * Dev (2019) _Neutrino Non-Standard Interactions: A Status Report_, Vol. 2 (2019) arXiv:1907.00991 [hep-ph] .
  * Berezhiani and Rossi (2002) Z. Berezhiani and A. Rossi, Phys. Lett. B 535, 207 (2002), arXiv:hep-ph/0111137 .
  * Davidson and Sanz (2011) S. Davidson and V. Sanz, Phys. Rev. D 84, 113011 (2011), arXiv:1108.5320 [hep-ph] .
  * Friedland  _et al._ (2012) A. Friedland, M. L. Graesser, I. M. Shoemaker, and L. Vecchi, Phys. Lett. B 714, 267 (2012), arXiv:1111.5331 [hep-ph] .
  * Buarque Franzosi  _et al._ (2016) D. Buarque Franzosi, M. T. Frandsen, and I. M. Shoemaker, Phys. Rev. D 93, 095001 (2016), arXiv:1507.07574 [hep-ph] .
  * Babu  _et al._ (2020) K. S. Babu, P. S. B. Dev, S. Jana, and A. Thapa, JHEP 03, 006 (2020), arXiv:1907.09498 [hep-ph] .
  * Esteban  _et al._ (2018) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, I. Martinez-Soler, and J. Salvado, JHEP 08, 180 (2018), [Addendum: JHEP 12, 152 (2020)], arXiv:1805.04530 [hep-ph] .
  * Friedland  _et al._ (2004) A. Friedland, C. Lunardini, and C. Pena-Garay, Phys. Lett. B 594, 347 (2004), arXiv:hep-ph/0402266 .
  * Friedland and Shoemaker (2012) A. Friedland and I. M. Shoemaker, (2012), arXiv:1207.6642 [hep-ph] .
  * Capozzi  _et al._ (2020) F. Capozzi, S. S. Chatterjee, and A. Palazzo, Phys. Rev. Lett. 124, 111801 (2020), arXiv:1908.06992 [hep-ph] .
  * Gonzalez-Garcia  _et al._ (2007) M. C. Gonzalez-Garcia, P. C. de Holanda, E. Masso, and R. Zukanovich Funchal, JCAP 01, 005 (2007), arXiv:hep-ph/0609094 .
  * Coloma  _et al._ (2021) P. Coloma, M. C. Gonzalez-Garcia, and M. Maltoni, JHEP 01, 114 (2021), arXiv:2009.14220 [hep-ph] .
  * Miranda  _et al._ (2006) O. G. Miranda, M. A. Tortola, and J. W. F. Valle, JHEP 10, 008 (2006), arXiv:hep-ph/0406280 .
  * Engel (1991) J. Engel, Phys. Lett. B 264, 114 (1991).
  * Farzan (2015) Y. Farzan, Phys. Lett. B 748, 311 (2015), arXiv:1505.06906 [hep-ph] .
  * Farzan and Shoemaker (2016) Y. Farzan and I. M. Shoemaker, JHEP 07, 033 (2016), arXiv:1512.09147 [hep-ph] .
  * Farzan and Heeck (2016) Y. Farzan and J. Heeck, Phys. Rev. D 94, 053010 (2016), arXiv:1607.07616 [hep-ph] .
  * Akimov  _et al._ (2018) D. Akimov _et al._ (COHERENT), (2018), arXiv:1803.09183 [physics.ins-det] .
  * Kamada and Yu (2015) A. Kamada and H.-B. Yu, Phys. Rev. D 92, 113004 (2015), arXiv:1504.00711 [hep-ph] .
  * Esteves Chaves and Schwetz (2021) M. Esteves Chaves and T. Schwetz, (2021), arXiv:2102.11981 [hep-ph] .
  * Bally  _et al._ (2020) A. Bally, S. Jana, and A. Trautner, Phys. Rev. Lett. 125, 161802 (2020), arXiv:2006.11919 [hep-ph] .
  * Shoemaker  _et al._ (2020) I. M. Shoemaker, Y.-D. Tsai, and J. Wyenberg, (2020), arXiv:2007.05513 [hep-ph] .
  * Berbig  _et al._ (2020) M. Berbig, S. Jana, and A. Trautner, Phys. Rev. D 102, 115008 (2020), arXiv:2004.13039 [hep-ph] .
  * Batell  _et al._ (2016) B. Batell, M. Pospelov, and B. Shuve, JHEP 08, 052 (2016), arXiv:1604.06099 [hep-ph] .