# Anomaly Induced Supernovae Cooling: New Contribution from Axions

Sabyasachi Chakraborty<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology, Kanpur-208016, India  Aritra Gupta
<EMAIL_ADDRESS>Instituto de Física Corpuscular (IFIC), CSIC, Parc
Científic, C/Catedrático José Beltrán, 2, E-46980 Paterna, Spain  Miguel
Vanvlasselaer<EMAIL_ADDRESS>Theoretische Natuurkunde and
IIHE/ELEM, Vrije Universiteit Brussel, & The International Solvay Institutes,
Pleinlaan 2, B-1050 Brussels, Belgium

###### Abstract

Compact stellar objects like supernovae and neutron stars are believed to cool
by emitting axions predominantly via axion bremsstrahlung ($NN\to NNa$), pion
conversion ($\pi^{-}p^{+}\to Na$) and photo-production ($\gamma N\to Na$). In
this letter, we study a previously overlooked contribution to the photo-
production channel, originating from the unavoidable anomaly-induced Wess-
Zumino-Witten term
$\propto\epsilon^{\mu\nu\alpha\beta}\,F_{\mu\nu}\,\partial_{\alpha}a\,\omega_{\beta}$.
As a result, supernovae can cool up to three times faster than what was
previously thought, translating into stronger bounds on the axion decay
constant. Furthermore, the spectrum of axions emitted in the process is unique
as it is significantly harder than those originating from bremsstrahlung. As a
consequence, we find that these axions are more likely to show up in near
future water Cherenkov detectors as compared to traditional ones.

Introduction.– Searches for new physics much lighter than the electroweak
scale have gained a lot of attention in the recent past. Among all the
possible candidates, axions or axion-like particles are perhaps the most well-
motivated. Axion emerges as the pseudo-Nambu-Goldstone boson of a
spontaneously broken $U(1)$ symmetry [1, 2, 3, 4] and offers a compelling
solution to the long-standing strong-$CP$ problem [5] of the Standard Model
(SM). Moreover, axions can also be a natural dark matter candidate [6, 7, 8],
address the hierarchy problem [9, 10, 11] and play a crucial role in resolving
the matter-antimatter asymmetry [12, 13] of the Universe. Naturally, axion
furnishes complementary means to probe physics beyond the SM at multiple
frontiers.

While the axion can have a myriad of couplings with the SM particles, the most
minimal effective Lagrangian addressing the strong $CP$ problem, up to terms
required for renormalization is given as:

$\displaystyle\mathcal{L}=\mathcal{L}_{\text{SM}}+\frac{\alpha_{s}}{8\pi}\frac{a}{f_{a}}G\widetilde{G}+\frac{1}{2}\left(\partial_{\mu}a\right)^{2}-\frac{1}{2}m_{a_{0}}^{2}a^{2}\;,$
(1)

where $G$ is the field strength tensor of the gluon and $a$ is the axion with
decay constant $f_{a}$. For lighter axions,
$m_{a}\lesssim\mathcal{O}(\text{GeV})$, it is convenient to rotate away the
$aG\widetilde{G}$ coupling by a chiral transformation of the SM quark fields.
The ensuing Lagrangian generates axion-quark kinetic and mass mixing terms and
can be matched with an effective theory such as Chiral perturbation theory
($\chi$PT) containing mesons. This effective description provides a powerful
framework to study non-perturbative effects such as axion-pion mixing and has
been used to constrain large parts of axion parameter space from intensity
frontier experiments such as beam dump [14, 15, 16], rare decays of pions [17,
18, 19] and kaons [20, 21, 22, 23, 24, 25, 26, 27], etc. However, for heavy
QCD axions, $m_{a}\gtrsim\mathcal{O}(1)$ GeV, the power counting of $\chi$PT
breaks down and one generally relies on pertubative computations to probe
axions [28, 29, 30, 31], mostly from rare $B$-decays.

On the other hand, cosmic frontier provides complementary means to probe
axions as it is sensitive to lighter masses, e.g.,
$m_{a}\lesssim\mathcal{O}(10)$ MeV. For example, nuclear reactions or thermal
processes inside the stellar interior such as White Dwarf (WD), Neutron Star
(NS), Supernova (SN), etc., are potentially powerful sources of axions. The
emission of axions from such compact stellar objects might result in a more
efficient transport of energy compared to the SM neutrinos, leading to
observational changes. This lead the authors of Ref. [32, 33] to propose the
following bound on the emissivity of axions

$\displaystyle{}\frac{Q_{a}}{\rho}\lesssim
10^{19}\;\text{erg}\;s^{-1}g^{-1}\;,$ (2)

evaluated at a temperature $T=30-40$ MeV and around the nuclear saturation
density. Traditionally, axion bremsstrahlung, i.e., $NN\to NNa$ [34, 35],
where $N=(p,n)$ is a nucleon, was considered to be the most dominant channel
of axion production in the core of the SN. Although previous studies assumed
One Pion Exchange (OPE) approximation, large suppression was found while going
beyond OPE [36]. Nevertheless, using Eq.(2), earlier bounds were drawn on the
axion mass and effective axion-nucleon coupling [37, 38, 39, 40, 41, 42]. More
recently, it was argued that pion conversion [43, 44] could potentially
enhance the emissivity by a factor of few [45, 46], so is the existence of
quark matter in the core of the SN [47]. We summarize those contributions and
the relevant operators in Table 1.

Process | Coupling | Refs
---|---|---
$NN\to NNa$ | $\left(C_{aN}/2f_{a}\right)\;\bar{N}\gamma_{5}\gamma_{\mu}N\;\partial^{\mu}a$ | [34, 35]
$\pi^{-}p\to Na$ | $\left(C_{aN}/2f_{a}\right)\;\bar{N}\gamma_{5}\gamma_{\mu}N\;\partial^{\mu}a$ | [40, 45]
$N\gamma\to Na$ | $i\;\left(C_{aN\gamma}/2\right)a\;\bar{N}\gamma_{5}\sigma_{\mu\nu}N\;F^{\mu\nu}$ | [48, 49]
$N\gamma\to Na$ | $\left(\kappa/f_{a}\right)\epsilon^{\mu\nu\rho\sigma}\;F_{\mu\nu}\;\partial_{\rho}a\;\omega_{\sigma}$ | This work
Table 1: Interactions relevant for axion emission from a supernova. The
coefficient $\kappa$ has been defined in Eq. (6).

In this minimal scenario, it is natural to ask whether there exists any other
unavoidable interactions of axions or not. For example, in the SM, there
exists Wess-Zumino-Witten (WZW) interactions [50, 51, 52, 53, 54, 55, 56, 57,
58, 59] which can account for certain processes including anomalies that can
not be generated in $\chi$PT. This is because low-energy effective theories
such as $\chi$PT can often possess more symmetries such as spurious parity
than its UV counterpart QCD. Interestingly, in the presence of background
gauge fields (such as $\omega$, $\rho$ mesons, etc.), the physics of WZW
interactions are quite rich [57]. Recently it was shown that even in SM, such
interactions open up new channels for the cooling of young neutron stars [56,
59]. In this Letter, we study the implications of WZW terms in the presence of
axions leading to novel interactions.

WZW interactions.– To generate WZW interactions, one considers the
5-dimensional action [50, 51] which is invariant under chiral symmetry, where
the boundary is identified with our 4-dimensional spacetime, e.g.,

$\displaystyle S_{\rm WZW}(A_{\mu},U)=N_{c}\int_{D}\;d^{5}y\;\omega\;,$
$\displaystyle\omega=-\frac{i}{240\pi^{2}}\;\epsilon^{\mu\nu\rho\sigma\tau}\;\text{Tr}\left(\mathcal{U}_{\mu}\;\mathcal{U}_{\nu}\;\mathcal{U}_{\rho}\;\mathcal{U}_{\sigma}\;\mathcal{U}_{\tau}\right)\;.$
(3)

Here $\mathcal{U}=U^{\dagger}\partial_{\mu}U$ with
$U=\text{Exp}\left(2i\pi^{a}T^{a}/f_{\pi}\right)$. $\pi_{a}$ are the pion
fields with decay constant $f_{\pi}$ and $A_{\mu}$’s are the gauge fields. The
coefficient $N_{c}=3$ is fixed by matching with QCD. Any arbitrary subgroup of
the chiral symmetry can be gauged but only in the 4-dimension using the trial
and error method. The full result has been nicely tabulated in [51] (also see
[50, 52, 53, 54, 55]). Meson fields such as $\omega$ can be introduced as a
background vector field as prescribed in Ref. [56] by replacing
$\tilde{A_{\mu}}=A_{\mu}+B_{\mu}$ in the effective action. If the fundamental
gauge fields are vectorlike, i.e., $A_{L}=A_{R}$, one needs to add appropriate
counter-terms to maintain gauge invariance and conservation of vector current.
Finally, tracking interactions with the pion fields (see Sup. Mat. for
details), we find

$\displaystyle\mathcal{L}_{WZW}^{\pi_{0}}$
$\displaystyle\supset\;\frac{N_{c}\;\epsilon^{\mu\nu\rho\sigma}}{24\pi^{2}}F_{\mu\nu}\left[\frac{e^{2}}{4}\frac{\pi_{0}}{f_{\pi}}\;F_{\rho\sigma}+eg_{\omega}\frac{\partial_{\rho}\pi_{0}}{f_{\pi}}\;\omega_{\sigma}\right]\;.$
(4)

Starting from Eq. (1), axions can be incorporated in this framework by a
suitable chiral transformation on the SM quark fields. This generates a mass-
mixing between the axion and pion fields (shown by a cross in Fig. (1)). The
mixing angle can be approximated as[60, 61, 62]

$\theta_{\pi^{0}-a}\simeq\frac{f_{\pi}}{2f_{a}}\bigg{(}\frac{m_{d}-m_{u}}{m_{u}+m_{d}}\bigg{)}\equiv\frac{f_{\pi}}{f_{a}}C_{A}\,.$
(5)

We then trade the the mixing angle between axion and pion fields, which
generates axion-photon coupling and the desired interactions of the form

$\mathcal{L}_{\rm
WZW}^{a}\supset\frac{\kappa}{f_{a}}\;\epsilon^{\mu\nu\rho\sigma}\;F_{\mu\nu}\;\partial_{\rho}a\;\omega_{\sigma}\;,\;\kappa=\frac{C_{A}N_{c}}{24\pi^{2}}\;eg_{\omega}\;.$
(6)

$N$$N$$\gamma$$a$$\pi^{0}$$\omega$ Figure 1: Cooling of supernovae and neutron
stars via $N\gamma\to Na$ induced by axion-WZW interactions. The cross
indicates pion-axion mixing.

On the other hand, nucleons interact with the vector meson fields $\omega$ via

$\mathcal{L}_{0}=\bar{N}\left(i\cancel{\partial}-g_{\omega}\;\cancel{\omega}-m_{N}\right)N\;,$
(7)

where $g_{\omega}$ is the coupling constant related to the pion-nucleon
coupling as $g_{\omega}\sim\left(9\,g_{\pi}/5\right)\sqrt{m_{\omega}/m_{\pi}}$
[63]. For this work, we follow the conventions of [34, 35] and adopt
$g_{\pi}\approx 13.72$.

Figure 2: Left panel: Comparison of the emissivity rates from axion
bremsstrahlung and pion conversion from [46, 64] with WZW photo-production
(blue line) and medium uncertainties (blue shaded region). Right panel: The
blue region shows the constraints drawn from WZW photo-production of axions
and the hatched region represents the medium uncertainty. We recast the
previous bounds [65] on axion parameter space from bremsstrahlung and pion
conversion in yellow. Bounds on effective axion-photon coupling
$g_{a\gamma\gamma}$ [66, 67], from HB stars [68, 69], CAST [70], Sumico [67]
and light shining through wall [71] experiment is converted using Eqs. (6) and
(5) onto the $m_{a}-f_{a}$ plane by identifying
$f_{a}=\left(2\alpha_{em}C_{A}/\pi\right)\Lambda$ (see [66], Fig. 28).
Exclusion limits from beam dump experiments [14, 15, 16], rare decays of $K$
[23, 27] and $B$-mesons [28, 29, 30] are shown in magenta and grey
respectively.

Cooling via Axions.– As depicted in Fig. 1, Eq. (6) and Eq. (7) provide an
efficient production channel of axions in a nucleon rich environment. The
coupling $f_{a}$, dictates whether the axions can come out of the stellar
object or get trapped inside it. In the former case, axions can carry a part
of the internal energy of the core with them, resulting in the cooling of the
star. To estimate, we compute the emissivity $Q$, defined as the emitted
energy per unit volume and per unit time. We take into account both neutron
degenerate (D) as well as non-degenerate (ND) scenarios, mimicking the
situation inside an NS and SN respectively. The emissivity for axion photo-
production $\gamma N\to Na$ is given by

$\displaystyle Q^{\rm WZW}_{N\gamma\to
Na}\equiv\int\frac{d^{3}\vec{p}_{\gamma}f_{\gamma}\left(p_{\gamma}\right)}{(2\pi)^{3}2E_{\gamma}}\int\frac{d^{3}\vec{p}_{N_{1}}}{(2\pi)^{3}2E_{N_{1}}}\frac{d^{3}\vec{p}_{N_{2}}}{(2\pi)^{3}2E_{N_{2}}}$
$\displaystyle
g_{\gamma}\,g_{N}\int\frac{d^{3}\vec{p}_{a}}{(2\pi)^{3}2E_{a}}E_{a}\langle|\mathcal{M}|^{2}\,\rangle
f_{N}(E_{N_{1}})\,(1-f_{N}(E_{N_{2}}))$
$\displaystyle\quad\times\,(2\pi)^{4}\delta(E_{N_{1}}-E_{N_{2}}-Q_{0})\,\delta^{3}(\vec{p}_{N_{1}}-\vec{p}_{N_{2}}-\vec{q})\,.$
(8)

where $g_{\gamma(N)}=2$ are the dof. of the photon (neutron), $f_{N}$ is the
neutron distribution function, $E_{i}$ and $\vec{p}_{i}$ are the energies and
three-momenta of the $i^{\rm th}$ species and $Q_{0}\equiv
E_{a}-E_{\gamma},\;\vec{q}\equiv\vec{p}_{a}-\vec{p}_{\gamma}$. The kinematics
of the process is similar to a fixed target experiment where the initial
neutron is assumed to be at rest. Due to the momentum transfer, the final
neutron receives a kick and we consider them to be nearly non-relativistic,
expanding in the energies $E_{N_{i}}\approx
m_{N}+\mathcal{O}(\vec{p}^{2}_{N_{i}}/2m_{N})$. The spin averaged amplitude
squared $\langle|\mathcal{M}|^{2}\rangle$ is

$\displaystyle\langle|\mathcal{M}^{\rm WZW}_{\gamma N\to Na}|^{2}\rangle$
$\displaystyle\approx\frac{8\kappa^{2}m_{N}^{2}}{f_{a}^{2}}\;\frac{g_{\omega}^{2}}{m_{\omega}^{4}}\;E_{\gamma}^{2}|\vec{p}_{a}|^{2}\;\sin^{2}\theta\;.$
(9)

Here, $\theta$ is the angle between the axion and photon three momenta. The
kinematics of the process dictates $E_{a}\sim E_{\gamma}$ when $m_{\gamma}\ll
T$. This is a reasonable approximation as in the core of SN, the photons
acquire a mass of the order of
$m_{\gamma}\sim\mathcal{O}(10)\;\text{MeV}\times(Y_{e}\rho_{14})^{1/3}$, where
$Y_{e}$ is the electron fraction. Since $T/m_{\gamma}\sim(3-4)$, we therefore
do not expect the mass of the photon to play a significant role.

Non-degenerate scenario.– For ND neutrons, the absence of Pauli blocking helps
us to further simplify Eq. (Anomaly Induced Supernovae Cooling: New
Contribution from Axions) as $1-f_{N}(p_{N})\approx 1$. Using the thermal
distribution for photons [72] and integrating the axion energy from
$(m_{a},\infty)$, we find

$\displaystyle Q_{N\gamma\to Na}^{\rm WZW,ND}\approx$
$\displaystyle\;g_{\gamma}g_{\omega}^{2}\;\frac{\kappa^{2}n_{b}}{2\pi^{3}}\;\frac{T^{2}}{f_{a}^{2}}\;\frac{m_{a}^{4}}{m_{\omega}^{4}}\;\left[4m_{a}T\;K_{3}\left(\frac{m_{a}}{T}\right)\right.$
$\displaystyle\left.+(m_{a}^{2}+35T^{2})\;K_{4}\left(\frac{m_{a}}{T}\right)\right]\;,$
(10)

where $K_{3,4}$ are the Bessel function of type three and four respectively.
Further, the number density reads as

$\displaystyle n_{b}=4.5\times
10^{-3}\rho_{15}\;\text{GeV}^{3},\quad\rho_{15}\equiv\rho/10^{15}\text{g
cm}^{-3}\,.$ (11)

In the limit of $m_{a}\to 0$, the expression in Eq.(10) further simplifies to

$\displaystyle\frac{Q_{\gamma N\to Na}^{\rm
WZW,ND}}{10^{34}\;\text{erg/s/cm}^{3}}\approx$
$\displaystyle\;5.6\;C_{A}^{2}\;g_{40}^{4}\;T_{40}^{8}\,\rho_{15}\left(\frac{10^{9}\;\text{GeV}}{f_{a}}\right)^{2}\;,$
(12)

here we define $T_{40}$ as $(T/40~{}\text{MeV})$ and
$g_{40}\equiv(g_{\omega}/40)$. The temperature dependence in Eq. (12) can be
understood as the phase space scales as $T^{3}$, the cross-section and emitted
energy as $T^{4}$ and $T$ respectively. However, the cooling argument in Eq.
(2) applies only when the axions can free-stream through the SN core and do
not get trapped[38]. The mean free path can be deduced from
$\ell=\left(n_{b}\sigma\right)^{-1}$ and can be approximated by

$\frac{L_{a}}{10\;\text{Km}}\approx\frac{\rho_{15}^{-1}}{C_{A}^{2}}\;\frac{1}{g_{40}^{2}}\left(\frac{f_{a}}{10^{9}~{}\text{GeV}}\right)^{2}\left(\frac{1~{}\text{GeV}}{E_{a}}\right)^{4}\;.$
(13)

This translates to an upper bound on the energy of the escaping axion.
Therefore, Eq. (10) was an estimate based only on free-streaming axions. We
impose this upper cut-off on the axion energy while deriving the final limits
on the axion parameter space, shown in the right panel Fig. 2.

Degenerate scenario.–The derivation of the cooling rate is more involved when
neutrons are degenerate inside a medium. We closely follow the derivation of
Ref. [59] (see also Sup. mat.) to write the emissivity in terms of the
momentum transfer and obtain,

$\displaystyle\frac{Q_{\gamma N\to Na}^{\rm
WZW,D}}{10^{34}\text{erg/s/cm}^{3}}\approx
1.7\;C_{A}^{2}\;g_{40}^{4}T_{40}^{9}\bigg{(}\frac{10^{9}~{}\text{GeV}}{f_{a}}\bigg{)}^{2}\,.$
(14)

The final emissivity is given by the minimum of the degenerate and the non-
degenerate expressions: $Q=\text{Min}[Q_{\rm D},Q_{\rm ND}]$. At the
saturation density $\rho\sim\rho_{0}\approx 0.3\times 10^{15}\text{g
cm}^{-3}$, we obtain that the rate is

$\displaystyle{}Q^{N\gamma\to Na}=\begin{cases}Q^{N\gamma\to Na}_{\rm D}\qquad
T\lesssim 40\text{MeV}\\\ Q^{N\gamma\to Na}_{\rm ND}\qquad T\gtrsim
40\text{MeV}\end{cases}$ (15)

and we conclude that for the SN, it suffices to use the ND limit. On the left
panel of Fig. 2, we show the emissivity of the photo-production and compare it
with pion conversion and bremsstrahlung.

At this point, it is worth mentioning that the high densities in the SN cores
lead to the modification of the masses and the couplings. The Brown-Rho
scaling dictates that the variation of the mass of the vector mesons follows
the neutron mass[73]:

$\displaystyle{}\frac{f_{\pi}^{\star}}{f_{\pi}}\approx\frac{m^{\star}_{\omega}}{m_{\omega}}\approx\frac{m^{\star}_{N}}{m_{N}}\Rightarrow\frac{g_{\omega}^{\star}}{g_{\omega}}\approx\frac{g_{\pi}^{\star}}{g_{\pi}}\approx\frac{g_{A}^{\star}}{g_{A}}\,.$
(16)

where the $\star$ indicates that the quantity is evaluated in the dense
environment of the NS [41, 72, 74]. The $m_{N}$ dependence is obtained from
[75, 74]. For the $g_{A}$ dependence, we use two different scaling laws: thick
line from[76, 41] and dashed line from[77, 78], $g_{A}^{\star}/g_{A}\approx
1/(1+(1/3)(m^{\star}_{N}/m_{N})(\rho/\rho_{0})^{1/3})$. In the case of the
axion bremsstrahlung, such effects were lead to only moderate uncertainties in
the coupling[79, 41]. For $\rho/\rho_{0}\in[0.5,1]$, the uncertainties in the
WZW photo production amounts to

$\displaystyle{}\frac{Q^{\star}_{N\gamma\to Na}}{Q_{N\gamma\to
Na}}\approx\frac{m_{\omega}^{4}}{g_{\omega}^{4}}\frac{g_{\omega\star}^{4}}{m_{\omega\star}^{4}}\in[0.8,1.5]\;.$
(17)

Implications.– We summarize our results in fig. 2 where we show WZW photo-
production with $T=40$ MeV can exclude regions with $f_{a}\lesssim 10^{9}$ GeV
for $m_{a}\lesssim 100$ MeV. The previous constraints, deduced from
simulations in[65], reaches $f_{a}\lesssim 0.77\times 10^{9}$ GeV for
$m_{a}\lesssim 50$ MeV. We see that the combination of the two rates would
bring a factor $\mathcal{O}(2-3)$ with respect to the previous bound from
bremsstrahlung and pion conversion only. We note in passing that axion photo-
production can also be induced via the electric dipole portal [49]. However,
the constraints are far subleading, $f_{a}/\text{GeV}\gtrsim 6\times 10^{5}$
compared to WZW process. As a consequence, we also neglect the possible cross
terms between the WZW and the electric dipole that could arise.

It is also important to mention that the landscape of QCD axion models is
quite broad. Without going into the details of such different scenarios, we
consider a minimal situation given by Eq. (1), similar to the KSVZ [80, 81]
framework. However, we find that similar conclusions hold in other types of
axion models such as DFSZ [82, 83] (see Sup. Mat.). This demonstrates the
robustness of the WZW term as far as cooling of SN is concerned.

The emissivity is related to the number density of particles emitted via

$\displaystyle{}Q_{\rm a}=\int dE_{a}E_{a}\frac{d^{2}n_{a}}{dE_{a}dt}\,,$ (18)

where the spectrum is defined as the number of axions emitted per interval of
energy and time. Comparing with Eq. (Anomaly Induced Supernovae Cooling: New
Contribution from Axions), it is straightforward to obtain the spectrum in the
ND limit and for $m_{a}\ll T$. Furthermore, the total emitted flux is well
approximated by integrating the spectrum over the radius of the star in the
optically thin limit. Taking $R_{\rm PNS}\approx 10$ Km, we obtain

$\displaystyle\frac{d^{2}N_{a}}{dE_{a}dt}\approx C_{f}$
$\displaystyle\rho_{15}\bigg{(}\frac{C_{A}10^{9}}{f_{a}/\text{GeV}}\bigg{)}^{2}g^{4}_{40}\bigg{(}\frac{E_{a}}{\text{MeV}}\bigg{)}^{6}e^{-E_{a}/T}\;,$
$\displaystyle C_{f}$ $\displaystyle=4.6\times
10^{42}\;\text{MeV}^{-1}\text{s}^{-1}\,.$ (19)

As shown in Fig. 3 by the blue curve, the total flux from the WZW photo-
production is distinctively different from axion bremsstrahlung, with a
comparatively harder spectrum.

Figure 3: Spectrum of axions due to WZW photo-production (blue) and
bremsstrahlung [74] (black).

Axion detection prospects.– Finally, we discuss the prospect of detecting such
axions in Water Cerenkov experiments [39] such as Hyper-Kamiokande with
fiducial mass 374 ktons [84]. Usually, axions emitted during future galactic
SN can interact with the detector via the channels $ap\to p\gamma$ and
$ap^{+}\to n\pi^{+}$ [49]. Around the peak of the WZW axion spectrum, i.e.,
$E_{a}\sim 200-300$ MeV, the cross section for $ap^{+}\to n\pi^{+}$ is
dominated by the $\Delta(1.23~{}\text{GeV})$ resonance and roughly translates
to $\sigma_{ap^{+}\to N\pi^{+}}\approx 10^{-25}C_{A}^{2}(f_{\pi}/f_{a})^{2}$
cm-2 [85]. Assuming $N_{t}$ target particles, the spectrum of produced
particles ($\gamma$ or $\pi$) from axions is

$\displaystyle\frac{d^{2}N_{p}}{dtdE_{p}}=\frac{N^{p}_{t}}{4\pi(d/\text{Kpc})^{2}}\times\frac{d^{2}N_{a}}{dtdE_{a}}\left(\frac{\sigma_{2\to
2}}{\text{cm}^{2}}\right),$ (20) $\displaystyle\;N^{p}_{t}\approx 2\times
10^{9}N_{a}\bigg{(}\frac{M_{\rm
detector}}{\text{kton}}\bigg{)}\bigg{(}\frac{\text{g}/\text{mol}}{m_{H_{2}O}}\bigg{)}\,.$

where $N_{a}$ is the Avogadro number, $m_{H_{2}O}=18\text{g}/\text{mol}$ the
molecular mass of water, and $d$ is the distance to the SN.

For small momentum transfer, the energy of the final state particle is roughly
equal to the energy of the incoming axion, i.e., $E_{p}\simeq E_{a}$. For
constant cross-sections with $100-500$ MeV axions, the produced particle will
inherit the spectrum of the emitted axions. The number of pions produced
during such interaction can be obtained by integrating Eq.(20) and we find

$\displaystyle{}\frac{dN_{\pi}}{dt}\approx
6\;\rho_{15}\;C_{A}^{4}\;\bigg{(}\frac{10^{9}}{f_{a}/\text{GeV}}\bigg{)}^{4}g^{4}_{40}\,T_{40}^{7}\bigg{(}\frac{\text{Kpc}}{d}\bigg{)}^{2}$
$\displaystyle\times\bigg{(}\frac{M_{\rm
detector}}{\text{kton}}\bigg{)}\bigg{(}\frac{\text{g}/\text{mol}}{m_{H_{2}O}}\bigg{)}\,s^{-1}\,.$
(21)

A SN at a distance of 1 Kpc with axions having $f_{a}=10^{9}$ GeV, would thus
roughly produce $\sim 5$ pions in the Hyper-Kamiokande, on a timescale of 10
seconds. On the other hand, the number of photons produced via $ap\rightarrow
p\gamma$ is always negligible because of a smaller cross-section.

Summary & outlook.– In summary, we point out a new channel for axion photo-
production, i.e., $\gamma N\to Na$, mediated by the WZW term
$\epsilon^{\mu\nu\alpha\beta}F_{\mu\nu}\;\partial_{\alpha}a\;\omega_{\beta}$.
We find that for temperatures $T\gtrsim 30$ MeV, the WZW photo-production
emissivity can compete and even dominate over the usual axion bremsstrahlung
process $NN\to NNa$, as depicted in Fig. 2 This is because the three-body
phase space suppression of the bremsstrahlung competes with the loop
suppression arising from the anomaly-induced WZW term. As a result, the WZW
process can increase the SN cooling rate by a factor of 2-3 and should be
considered in the SN simulations [86]. Moreover, the spectrum of such axions
is distinctively different compared to the prototypical SN cooling processes.
As a result, WZW axions emitted from a close-by SN can leave a smoking-gun
signature in the near future Water Cherenkov experiments, producing
$\mathcal{O}(10)$ pions in the final state. Photo-production of axions can
also be probed at GlueX[87], where the photon beam of $\sim 10$ GeV hits a
fixed target and produces a plethora of mesons. As far as degenerate neutron
stars are concerned, we find that the cooling rate scales as $T^{9}$ and is
subdominant compared to $NN\to NNa$ for temperatures $T\lesssim 30$ MeV.

At this point, it is pertinent to ask whether there exist SM WZW interactions
such as $\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}\omega_{\rho}Z_{\sigma}$ that
can lead to SN cooling. Such interactions induce the process $\gamma N\to
N\nu\bar{\nu}$ resulting in the emission of a pair of neutrinos from the decay
of an off-shell $Z$ boson. It was found in [59], that the emissivity from such
interactions can play a significant role for NS with degenerate stellar cores
when $T>10-20$ MeV. For SN, WZW induced SM photo-production of neutrinos can
dominate compared to other processes around $T\sim 40$ MeV. However, further
investigation is required as the neutrinos at such temperatures might get
trapped [88]. Finally, assuming purely non-relativistic neutrons and a large
magnetic field, the axion WZW term reduces to a mixing between $a-\omega$
field or can also contribute to axion-photon oscillation in the presence of a
varying neutron density. Naturally, the phenomenological implications are
rich, interesting and need to be explored further.

SC thanks the Science and Engineering Research Board, Govt. of India
(SRG/2023/001162) for financial support. MV is supported by the “Excellence of
Science - EOS” - be.h project n.30820817, and by the Strategic Research
Program High-Energy Physics of the Vrije Universiteit Brussel. AG is supported
by the “Generalitat Valenciana” and CSIC through the GenT Excellence Program
(CIDEIG/2022/22). MV thanks IFIC (Valencia) and CSIC for their hospitality
during the completion of this project and Alberto Mariotti for comments on the
draft.

## References

  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg [1978] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * ’t Hooft [1976] G. ’t Hooft, Phys. Rev. Lett. 37, 8 (1976).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Graham  _et al._ [2015a] P. W. Graham, D. E. Kaplan, and S. Rajendran, Phys. Rev. Lett. 115, 221801 (2015a), arXiv:1504.07551 [hep-ph] .
  * Hook and Marques-Tavares [2016] A. Hook and G. Marques-Tavares, JHEP 12, 101 (2016), arXiv:1607.01786 [hep-ph] .
  * Trifinopoulos and Vanvlasselaer [2023] S. Trifinopoulos and M. Vanvlasselaer, Phys. Rev. D 107, L071701 (2023), arXiv:2210.13484 [hep-ph] .
  * Co and Harigaya [2020] R. T. Co and K. Harigaya, Phys. Rev. Lett. 124, 111602 (2020), arXiv:1910.02080 [hep-ph] .
  * Chakraborty  _et al._ [2022] S. Chakraborty, T. H. Jung, and T. Okui, Phys. Rev. D 105, 015024 (2022), arXiv:2108.04293 [hep-ph] .
  * Bjorken  _et al._ [1988] J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker, and P. Rassmann, Phys. Rev. D 38, 3375 (1988).
  * Blumlein  _et al._ [1991] J. Blumlein _et al._ , Z. Phys. C 51, 341 (1991).
  * Bergsma  _et al._ [1985] F. Bergsma _et al._ (CHARM), Phys. Lett. B 157, 458 (1985).
  * Aguilar-Arevalo  _et al._ [2019] A. Aguilar-Arevalo _et al._ (PIENU), Phys. Lett. B 798, 134980 (2019), arXiv:1904.03269 [hep-ex] .
  * Pocanic  _et al._ [2004] D. Pocanic _et al._ , Phys. Rev. Lett. 93, 181803 (2004), arXiv:hep-ex/0312030 .
  * Altmannshofer  _et al._ [2020] W. Altmannshofer, S. Gori, and D. J. Robinson, Phys. Rev. D 101, 075002 (2020), arXiv:1909.00005 [hep-ph] .
  * Georgi  _et al._ [1986] H. Georgi, D. B. Kaplan, and L. Randall, Phys. Lett. B 169, 73 (1986).
  * Bardeen  _et al._ [1987] W. A. Bardeen, R. D. Peccei, and T. Yanagida, Nucl. Phys. B 279, 401 (1987).
  * Alves and Weiner [2018] D. S. M. Alves and N. Weiner, JHEP 07, 092 (2018), arXiv:1710.03764 [hep-ph] .
  * Gori  _et al._ [2020] S. Gori, G. Perez, and K. Tobioka, JHEP 08, 110 (2020), arXiv:2005.05170 [hep-ph] .
  * Artamonov  _et al._ [2005] A. V. Artamonov _et al._ (E949), Phys. Lett. B 623, 192 (2005), arXiv:hep-ex/0505069 .
  * Lazzeroni  _et al._ [2014] C. Lazzeroni _et al._ (NA62), Phys. Lett. B 732, 65 (2014), arXiv:1402.4334 [hep-ex] .
  * Ahn  _et al._ [2019] J. K. Ahn _et al._ (KOTO), Phys. Rev. Lett. 122, 021802 (2019), arXiv:1810.09655 [hep-ex] .
  * Bauer  _et al._ [2021] M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, Phys. Rev. Lett. 127, 081803 (2021), arXiv:2102.13112 [hep-ph] .
  * Aloni  _et al._ [2019a] D. Aloni, Y. Soreq, and M. Williams, Phys. Rev. Lett. 123, 031803 (2019a), arXiv:1811.03474 [hep-ph] .
  * Chakraborty  _et al._ [2021] S. Chakraborty, M. Kraus, V. Loladze, T. Okui, and K. Tobioka, Phys. Rev. D 104, 055036 (2021), [Erratum: Phys.Rev.D 108, 039903 (2023)], arXiv:2102.04474 [hep-ph] .
  * Bertholet  _et al._ [2022] E. Bertholet, S. Chakraborty, V. Loladze, T. Okui, A. Soffer, and K. Tobioka, Phys. Rev. D 105, L071701 (2022), arXiv:2108.10331 [hep-ph] .
  * Bauer  _et al._ [2022] M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 09, 056 (2022), arXiv:2110.10698 [hep-ph] .
  * Raffelt [1990] G. G. Raffelt, physrep 198, 1 (1990).
  * Raffelt [1996] G. G. Raffelt,  _Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles_ (1996).
  * Iwamoto [1984] N. Iwamoto, Phys. Rev. Lett. 53, 1198 (1984).
  * Iwamoto [2001a] N. Iwamoto, Phys. Rev. D 64, 043002 (2001a).
  * Hanhart  _et al._ [2001] C. Hanhart, D. R. Phillips, and S. Reddy, Phys. Lett. B 499, 9 (2001), arXiv:astro-ph/0003445 .
  * Fischer  _et al._ [2016] T. Fischer, S. Chakraborty, M. Giannotti, A. Mirizzi, A. Payez, and A. Ringwald, Phys. Rev. D 94, 085012 (2016), arXiv:1605.08780 [astro-ph.HE] .
  * Burrows  _et al._ [1990] A. Burrows, M. T. Ressell, and M. S. Turner, Phys. Rev. D 42, 3297 (1990).
  * Engel  _et al._ [1990] J. Engel, D. Seckel, and A. C. Hayes, Phys. Rev. Lett. 65, 960 (1990).
  * Keil  _et al._ [1997] W. Keil, H.-T. Janka, D. N. Schramm, G. Sigl, M. S. Turner, and J. Ellis, Phys. Rev. D 56, 2419 (1997).
  * Mayle  _et al._ [1989] R. Mayle, J. R. Wilson, J. Ellis, K. A. Olive, D. N. Schramm, and G. Steigman, Physics Letters B 219, 515 (1989).
  * Turner  _et al._ [1989] M. S. Turner, H.-S. Kang, and G. Steigman, Phys. Rev. D 40, 299 (1989).
  * Turner [1992] M. S. Turner, Phys. Rev. D 45, 1066 (1992).
  * Raffelt and Seckel [1995] G. Raffelt and D. Seckel, Phys. Rev. D 52, 1780 (1995).
  * Carenza  _et al._ [2021] P. Carenza, B. Fore, M. Giannotti, A. Mirizzi, and S. Reddy, Phys. Rev. Lett. 126, 071102 (2021).
  * Choi  _et al._ [2022] K. Choi, H. J. Kim, H. Seong, and C. S. Shin, JHEP 02, 143 (2022), arXiv:2110.01972 [hep-ph] .
  * Cavan-Piton  _et al._ [2024] M. Cavan-Piton, D. Guadagnoli, M. Oertel, H. Seong, and L. Vittorio, (2024), arXiv:2401.10979 [hep-ph] .
  * Graham and Rajendran [2013] P. W. Graham and S. Rajendran, Phys. Rev. D 88, 035023 (2013), arXiv:1306.6088 [hep-ph] .
  * Lucente  _et al._ [2022] G. Lucente, L. Mastrototaro, P. Carenza, L. Di Luzio, M. Giannotti, and A. Mirizzi, Phys. Rev. D 105, 123020 (2022), arXiv:2203.15812 [hep-ph] .
  * Wess and Zumino [1971] J. Wess and B. Zumino, Phys. Lett. B 37, 95 (1971).
  * Witten [1983] E. Witten, Nucl. Phys. B 223, 422 (1983).
  * Kaymakcalan  _et al._ [1984] O. Kaymakcalan, S. Rajeev, and J. Schechter, Phys. Rev. D 30, 594 (1984).
  * Chou  _et al._ [1984] K.-c. Chou, H.-y. Guo, K. Wu, and X.-c. Song, Phys. Lett. B 134, 67 (1984).
  * Kawai and Tye [1984] H. Kawai and S. H. H. Tye, Phys. Lett. B 140, 403 (1984).
  * Pak and Rossi [1985] N. K. Pak and P. Rossi, Nucl. Phys. B 250, 279 (1985).
  * Harvey  _et al._ [2008] J. A. Harvey, C. T. Hill, and R. J. Hill, Phys. Rev. D 77, 085017 (2008), arXiv:0712.1230 [hep-th] .
  * Harvey  _et al._ [2007] J. A. Harvey, C. T. Hill, and R. J. Hill, Phys. Rev. Lett. 99, 261601 (2007), arXiv:0708.1281 [hep-ph] .
  * Hill and Hill [2007] C. T. Hill and R. J. Hill, Phys. Rev. D 76, 115014 (2007), arXiv:0705.0697 [hep-ph] .
  * Chakraborty  _et al._ [2023] S. Chakraborty, A. Gupta, and M. Vanvlasselaer, Journal of Cosmology and Astroparticle Physics 2023, 030 (2023).
  * Di Luzio  _et al._ [2020] L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Notari  _et al._ [2023] A. Notari, F. Rompineve, and G. Villadoro, Phys. Rev. Lett. 131, 011004 (2023), arXiv:2211.03799 [hep-ph] .
  * Grilli di Cortona  _et al._ [2016] G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, JHEP 01, 034 (2016), arXiv:1511.02867 [hep-ph] .
  * Downum  _et al._ [2006] C. Downum, T. Barnes, J. R. Stone, and E. S. Swanson, Phys. Lett. B 638, 455 (2006), arXiv:nucl-th/0603020 .
  * Iwamoto [2001b] N. Iwamoto, Phys. Rev. D 64, 043002 (2001b).
  * Lella  _et al._ [2023] A. Lella, P. Carenza, G. Co’, G. Lucente, M. Giannotti, A. Mirizzi, and T. Rauscher, (2023), arXiv:2306.01048 [hep-ph] .
  * Agrawal  _et al._ [2021] P. Agrawal _et al._ , Eur. Phys. J. C 81, 1015 (2021), arXiv:2102.12143 [hep-ph] .
  * Graham  _et al._ [2015b] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Ann. Rev. Nucl. Part. Sci. 65, 485 (2015b), arXiv:1602.00039 [hep-ex] .
  * Ayala  _et al._ [2014] A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi, and O. Straniero, Phys. Rev. Lett. 113, 191302 (2014), arXiv:1406.6053 [astro-ph.SR] .
  * Straniero  _et al._ [2015] O. Straniero, A. Ayala, M. Giannotti, A. Mirizzi, and I. Dominguez, in  _11th Patras Workshop on Axions, WIMPs and WISPs_ (2015) pp. 77–81.
  * Anastassopoulos  _et al._ [2017] V. Anastassopoulos _et al._ (CAST), Nature Phys. 13, 584 (2017), arXiv:1705.02290 [hep-ex] .
  * Redondo and Ringwald [2011] J. Redondo and A. Ringwald, Contemp. Phys. 52, 211 (2011), arXiv:1011.3741 [hep-ph] .
  * Payez  _et al._ [2015] A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, and A. Ringwald, JCAP 02, 006 (2015), arXiv:1410.3747 [astro-ph.HE] .
  * Brown and Rho [1991] G. E. Brown and M. Rho, Phys. Rev. Lett. 66, 2720 (1991).
  * Raffelt [2001] G. G. Raffelt, The Astrophysical Journal 561, 890 (2001).
  * 198 [1987] nphysa 473, 760 (1987).
  * Rho [1985] M. Rho, Phys. Rev. Lett. 54, 767 (1985).
  * Voskresensky [2001] D. N. Voskresensky, Lect. Notes Phys. 578, 467 (2001), arXiv:astro-ph/0101514 .
  * Fischer [2016] T. Fischer, Astron. Astrophys. 593, A103 (2016), arXiv:1608.05004 [astro-ph.HE] .
  * Mayle  _et al._ [1988] R. Mayle, J. R. Wilson, J. Ellis, K. Olive, D. N. Schramm, and G. Steigman, Physics Letters B 203, 188 (1988).
  * Kim [1979] J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ [1980] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Zhitnitsky [1980] A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980).
  * Dine  _et al._ [1981] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Abe  _et al._ [2011] K. Abe _et al._ , (2011), arXiv:1109.3262 [hep-ex] .
  * Henley and Garcia [2007] E. M. Henley and A. Garcia,  _Subatomic Physics: Third Edition_ (2007).
  * [86] S. Chakraborty, A. Gupta, and M. Vanvlasselaer, .
  * Aloni  _et al._ [2019b] D. Aloni, C. Fanelli, Y. Soreq, and M. Williams, Phys. Rev. Lett. 123, 071801 (2019b), arXiv:1903.03586 [hep-ph] .
  * Hannestad and Raffelt [1998] S. Hannestad and G. Raffelt, Astrophys. J. 507, 339 (1998), arXiv:astro-ph/9711132 .

Anomaly Induced Supernovae Cooling: New Contribution from Axions

Supplemental Material
Sabyasachi Chakraborty, Aritra Gupta and Miguel Vanvlasselaer

##  I Wess-Zumino-Witten Interactions

The dynamics and interactions of Goldstones bosons associated with the
spontaneous breaking of chiral symmetry is aptly described by the Chiral
Lagrangian ($\chi$PT). At the leading order, the Lagrangian consists of terms

$\mathcal{L}=\frac{f_{\pi}^{2}}{4}\text{Tr}\left(D_{\mu}U^{\dagger}D^{\mu}U+U^{\dagger}\chi+\chi^{\dagger}U\right)\;,$
(22)

where $D_{\mu}U=\partial_{\mu}U-ir_{\mu}U+iU\ell_{\mu}$. $r_{\mu}$ and
$\ell_{\mu}$ describe left and right-handed currents which can be used to
construct the vector and axial vector currents. Moreover,
$\chi=2B_{0}\mathcal{M}$, where $\chi$ is the isospin symmetry breaking
spurion with $\mathcal{M}$ denoting the quark mass matrix. Constants $B_{0}$
and $f_{\pi}$ are fixed by experiments. However, Eq. (22) possesses more
symmetries such as spurious parity and charge conjugation compared to its
underlying UV theory, i.e., QCD. As a result, $\chi$PT in Eq. 22 fails to
describe processes such as $K^{+}K^{-}\to\pi^{+}\pi^{-}\pi^{0}$ and
$\pi^{0}\to\gamma\gamma$, which are otherwise allowed in QCD. Such
interactions can be restored in the Chiral Lagrangian by the Wess-Zumino-
Witten term. In the geometric representation of the WZW form, one considers a
5-dimensional action, whose boundary is identified with 4-dimensional space in
the form

$S_{WZW}(U)=\kappa\int_{D}d^{5}y\;\omega\;,\;\text{where}\;\;\omega=-\frac{i}{240\pi^{2}}\epsilon^{\mu\nu\rho\sigma\tau}\;\text{Tr}\left[\;\mathcal{U}_{\mu}\;\mathcal{U}_{\nu}\;\mathcal{U}_{\rho}\;\mathcal{U}_{\sigma}\;\mathcal{U}_{\tau}\right]\;.$
(23)

Here, we define $\mathcal{U}_{\mu}=U^{\dagger}\partial_{\mu}U$ and
$\mathcal{U}\approx\text{Exp}\left(2i\pi^{a}T^{a}/f_{\pi}\right)$ with
$f_{\pi}=93$ MeV. Under chiral symmetry, $U$ transforms as $U\to
RUL^{\dagger}$ with $R$ and $L$ describing the $SU(2)_{R}$ and $SU(2)_{L}$
transformation matrices. Note that Eq. (23) correctly describes the process
$K^{+}K^{-}\to\pi^{+}\pi^{-}\pi^{0}$. To generate $\pi^{0}\to\gamma\gamma$ or
similar processes, one has to gauge $S_{WZW}$. However, any gauging has to be
done in the 4-dimension. In principle, any arbitrary subgroup of the WZW
action can be gauged and the resultant terms are

$S_{WZW}\left(U,A_{L},A_{R}\right)=S_{WZW}\left(U\right)+\frac{N_{c}}{48\pi^{2}}\int
d^{4}x\;\epsilon^{\mu\nu\rho\sigma}\;Z_{\mu\nu\rho\sigma}\;,$ (24)

where the full expression of $Z_{\mu\nu\rho\sigma}$ has been nicely tabulated
in Ref. [51]. However, for our work, we only track the pion fields and
therefore the relevant terms are

$\displaystyle Z_{\mu\nu\rho\sigma}(U,A_{L},A_{R})\supset i\;$
$\displaystyle\text{Tr}\left[\left\\{\left(\partial_{\mu}A_{\nu_{L}}\right)A_{\rho_{L}}+A_{\mu_{L}}(\partial_{\nu}A_{\rho_{L}})\right\\}\partial_{\sigma}U\;U^{-1}+\left\\{\left(\partial_{\mu}A_{\nu_{R}}\right)A_{\rho_{R}}+A_{\mu_{R}}(\partial_{\nu}A_{\rho_{R}})\right\\}U^{-1}\partial_{\sigma}U\right]$
$\displaystyle+$ $\displaystyle
i\;\text{Tr}\left[\left\\{(\partial_{\mu}A_{\nu_{R}})\;U^{-1}A_{\rho_{L}}\partial_{\sigma}U+A_{\mu_{L}}U^{-1}(\partial_{\nu}A_{\rho_{R}})\;\partial_{\sigma}U\right\\}\right]\;.$
(25)

Expanding in pion fields, $\partial_{\mu}U$ can be written as

$\partial_{\mu}U=i\frac{\partial_{\mu}\pi_{0}}{f_{\pi}}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}+...\,.$ (26)

This choice results in a canonically normalized kinetic term from $\pi_{0}$.
Moreover, in the SM, the gauge fields are expressed as

$\displaystyle
A_{L}=g_{2}W^{a}\frac{\tau^{a}}{2}+g_{1}W^{0}\;\begin{pmatrix}\frac{1}{6}&0\\\
0&\frac{1}{6}\end{pmatrix}\;,\;\;\;A_{R}=g_{1}W^{0}\;\begin{pmatrix}\frac{2}{3}&0\\\
0&-\frac{1}{3}\end{pmatrix}\;.$ (27)

The physical gauge fields $A$ and $Z$ are introduced after electroweak
symmetry breaking and they are related with $W^{0}$ and $W^{3}$ fields via the
weak rotation angle $\theta_{w}$ as

$W^{0}=-s_{W}Z+c_{W}A,\qquad W^{3}=c_{W}Z+s_{W}A\,.$ (28)

Using $g_{1}c_{W}=g_{2}s_{W}=e$, Eq. (27) becomes

$A_{L}=g_{2}s_{W}\frac{A}{2}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}+g_{1}c_{W}A\begin{pmatrix}\frac{1}{6}&0\\\
0&\frac{1}{6}\end{pmatrix}=eA\begin{pmatrix}\frac{2}{3}&0\\\
0&-\frac{1}{3}\end{pmatrix}=A_{R}\,,$ (29)

To introduce the vector mesons fields, we use the background field method
presented in Ref. [56] where $\widetilde{A}_{L,R}=A_{L,R}+B_{L,R}$ (the charge
as been included inside the definition of the field), including both the
fundamental and background gauge fields. As a result, we make the following
replacement $A_{L,R}\to\widetilde{A}_{L,R}$ in Eq. (25). The background field
$\omega$ is then identified as

$\displaystyle B_{V}=\frac{2g_{\omega}}{3}\;\omega\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\,,\qquad B_{A}=0\,,\;\;\;\Rightarrow
B_{L}=B_{R}=\frac{B_{V}}{2}=\frac{g_{\omega}}{3}\;\omega\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\,.$ (30)

Using the definition of $A_{L,R}$ and $B_{L,R}$ as well as $N_{c}=3$, fixed by
matching with QCD, we find from Eq. (25)

$\mathcal{L}_{\rm
WZW}\supset\frac{N_{c}}{24\pi^{2}}eg_{\omega}\frac{\partial_{\rho}\pi_{0}}{f_{\pi}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}\;\omega_{\sigma}+\frac{N_{c}}{96\pi^{2}}e^{2}\;\epsilon^{\mu\nu\rho\sigma}\;\frac{\pi_{0}}{f_{\pi}}F_{\mu\nu}F_{\rho\sigma}\;.$
(31)

The coefficient $N_{c}=3$ is fixed by matching with QCD.

##  II Computation of the emissivities for the axion

In this section, we present the computation of the emissivities of the WZW
photo-production reaction $\gamma N\to Na$ for the non-degenerate case and
then for the degenerate case.

###  II.1 Emissivity in the non-degenerate case

We study the $2\to 2$ interaction $\gamma N\to Na$ in the non-degenerate
limit. In the centre of mass frame of the collision,
$E_{CM}=E_{N_{1}}+E_{\gamma}=E_{N_{2}}+E_{a}$. The matrix element takes the
form

$\displaystyle\langle|\mathcal{M}|^{2}\rangle$
$\displaystyle=8\frac{\kappa^{2}m_{N}^{2}g_{\omega}^{2}}{f_{a}^{2}m_{\omega}^{4}}\bigg{(}E_{\gamma}^{2}|\vec{p}_{a}|^{2}\sin^{2}\theta\bigg{)}+\mathcal{O}(p_{F}^{2}/m_{N}^{2}),\qquad\kappa\equiv
C_{A}\frac{eN_{c}}{24\pi^{2}}g_{\omega}\,,$ (32)

where we neglect $m_{\gamma}$ compared to the energy scale. $\theta$ is the
angle between the incoming neutron and the outgoing axion, $p_{F}$ is the
Fermi momentum of neutrons. Also, $E_{i}$ and $\vec{p_{i}}$ are the energy and
three-momenta of the $i^{\rm th}$ species. In the CM frame

$E_{\rm CM}\equiv E_{N_{1}}+E_{\gamma}=E_{N_{2}}+E_{a}\;,\qquad
E_{\gamma}\approx E_{a}\;,\qquad E_{a}^{2}=p_{a}^{2}+M_{a}^{2}\,.$ (33)

The cross-section for the $N\gamma\to Na$ process is given by

$\displaystyle\frac{d\sigma_{N\gamma\to Na}}{d\Omega}$
$\displaystyle=\frac{|\vec{p}_{a}|}{|\vec{p}_{\gamma}|}\frac{\langle|\mathcal{M}|^{2}\rangle}{64\pi^{2}E_{CM}^{2}}\theta\left(E_{cm}-m_{N}-m_{a}\right)\;,$
$\displaystyle=\frac{\sqrt{E_{a}^{2}-m_{a}^{2}}}{E_{\gamma}}\frac{\langle|\mathcal{M}|^{2}\rangle}{64\pi^{2}E_{CM}^{2}}\theta(E_{CM}-m_{N}-m_{a})\;,$
$\displaystyle\approx\frac{\sqrt{E_{a}^{2}-m_{a}^{2}}}{E_{\gamma}}\frac{\langle|\mathcal{M}|^{2}\rangle}{64\pi^{2}m_{N}^{2}}\theta(E_{CM}-m_{N}-m_{a})\,.$
(34)

The emissivity rate in the non-degenerate limit is expressed as

$\displaystyle Q^{\gamma N\to Na}_{\rm ND}$
$\displaystyle\approx\int\frac{g_{\gamma}d^{3}\vec{p}_{\gamma}}{(2\pi)^{3}}f_{\gamma}(p_{\gamma})\int\frac{g_{N}d^{3}\vec{p}_{N}}{(2\pi)^{3}}f_{N}\int\frac{\langle|\mathcal{M}|^{2}\rangle}{4E_{N_{1}}E_{\gamma}}(2\pi)^{4}\delta^{4}\left(\sum_{i}p_{i}\right)\frac{d^{3}\vec{p_{2}}d^{3}\vec{p_{a}}}{(2\pi)^{6}4E_{a}E_{N_{2}}}E_{a}\,.$
(35)

which can be further simplified with

$\displaystyle Q^{\gamma N\to Na}_{\rm ND}$
$\displaystyle\approx\int\frac{g_{\gamma}d^{3}\vec{p_{\gamma}}}{(2\pi)^{3}}f_{\gamma}(p_{\gamma})\int\frac{g_{N}d^{3}\vec{p}_{N}}{(2\pi)^{3}}f_{N}\sigma(E_{\gamma})E_{a}\;,\qquad\sigma(E_{\gamma})\approx\frac{1}{3\pi}\frac{\kappa^{2}g_{\omega}^{2}}{f_{a}^{2}m_{\omega}^{4}}\sqrt{1-\frac{M_{a}^{2}}{E_{\gamma}^{2}}}E_{\gamma}^{2}(E_{\gamma}^{2}-M_{a}^{2})\,,$
(36)

where $\sigma(E_{\gamma})$ is the integrated cross section. The nucleon
density is given by

$n_{B}\equiv\int\frac{d^{3}p_{N}}{(2\pi)^{3}}g_{N}f_{N},\qquad n_{B}=4.5\times
10^{6}\rho_{15}\text{ MeV}^{3}\,.\\\ $ (37)

Plugging Eq.(37) in the emissivity Eq.(36), we obtain

$\displaystyle Q^{N\gamma\to Na}_{\rm ND}$
$\displaystyle\approx\frac{g_{\gamma}n_{B}}{3\pi}\frac{\kappa^{2}g_{\omega}^{2}}{f_{a}^{2}m_{\omega}^{4}}\int^{\infty}_{m_{a}}\frac{dE_{\gamma}}{2\pi^{2}}f_{\gamma}(p_{\gamma})E_{\gamma}^{5}(E_{\gamma}^{2}-m_{a}^{2})\sqrt{1-\frac{m_{a}^{2}}{E_{\gamma}^{2}}}\,.$
(38)

The two limiting conditions are:

  1. 1.

$m_{a}/T\to 0$: In the limit $m_{a}/T\to 0$, we can approximate
$\sqrt{1-\frac{m_{a}^{2}}{E_{\gamma}^{2}}}\to 1$. This results in

$\displaystyle Q^{N\gamma\to Na}_{\rm
ND}\approx\frac{g_{\gamma}n_{B}}{3\pi^{3}}\frac{4\pi^{8}}{15}\frac{\kappa^{2}g_{\omega}^{2}}{f_{a}^{2}m_{\omega}^{4}}T^{8}\approx
27g_{\gamma}n_{B}\frac{\kappa^{2}}{f_{a}^{2}}T^{8}\;.$ (39)

Using the values and the normalization of the parameters explained in the main
text along with $g_{\gamma}=2$, we get

$Q^{N\gamma\to Na}_{\rm ND}\approx 12\times
C^{2}_{A}\bigg{(}\frac{g_{\omega}}{10}\bigg{)}^{4}\frac{n_{B}}{m_{\omega}^{4}f_{a}^{2}}T^{8}\;\Rightarrow\frac{Q^{N\gamma\to
Na}_{\rm ND}}{5.6\times 10^{34}\text{erg/s/
cm}^{3}}\approx\bigg{(}\frac{C_{A}m_{N}10^{9}}{f_{a}}\bigg{)}^{2}\times
g_{40}^{4}\rho_{15}T_{40}^{8}\,.$ (40)

  2. 2.

$m_{a}/T\gtrsim 1$: In the opposite limit $m_{a}/T\gtrsim 1$, the mass of the
axion cannot be neglected anymore and the expression becomes

$\displaystyle\frac{Q^{N\gamma\to Na}_{\rm ND}}{1.1\times 10^{31}\text{erg/s/
cm}^{3}}\approx\bigg{(}\frac{C_{A}m_{N}10^{9}}{f_{a}}\bigg{)}^{2}g_{40}^{4}\times\rho_{15}G(m_{a},T_{40})\,,$
(41)

where we defined the $G$ function in terms of the Bessel functions as

$G(m_{a},T)\equiv
3m_{a}^{4}T^{2}\bigg{(}4m_{a}TK_{3}\bigg{(}\frac{m_{a}}{T}\bigg{)}+(m_{a}^{2}+35T^{2})K_{4}\bigg{(}\frac{m_{a}}{T}\bigg{)}\bigg{)}\,.$
(42)

###  II.2 Emissivity in the degenerate case

We now evaluate the emissivity in the presence of strongly degenerate
neutrons. The full expression is provided in a convenient form as

$\displaystyle Q^{N\gamma\to Na}_{\rm
axion}=2\frac{\kappa^{2}}{f_{a}^{2}}\int\dfrac{d^{3}\vec{p}_{\gamma}}{(2\pi)^{3}}g_{\gamma}\dfrac{f_{\gamma}}{2E_{\gamma}}\int\dfrac{d^{3}\vec{p}_{a}}{(2\pi)^{3}2E_{a}}\times\bigg{[}2E_{\gamma}E_{a}(p_{\gamma}\cdot
p_{a})-(p_{\gamma}\cdot p_{a})^{2}\bigg{]}\,S(q^{\mu})E_{a}\,.$ (43)

where we define $Q_{0}\equiv
E_{2}-E_{1}=E_{\gamma}-E_{a}\;,\qquad\vec{q}=\vec{p}_{\gamma}-\vec{p}_{a}$.
The response function is defined as

$\displaystyle S(q^{\mu})\equiv$ $\displaystyle
g_{N_{1}}\int\dfrac{d^{3}\vec{p}_{N_{1}}d^{3}\vec{p}_{N_{2}}m_{N}^{2}}{(2\pi)^{6}\,E_{N_{1}}\,E_{N_{2}}}f_{N}(E_{N_{1}})(1-f_{N}(E_{N_{2}}))\times(2\pi)^{4}\times$
$\displaystyle\delta(E_{N1}+Q_{0}-E_{N_{2}})\delta^{3}(\vec{p}_{N_{1}}-\vec{p}_{N_{2}}+\vec{q})\,\Theta(E_{N_{1}}-m_{N})\,\Theta(E_{N_{2}}-m_{N})\;,$
(44)

and can be computed while identifying $z=Q_{0}/T$ and $|\vec{q}|\equiv q$. We
obtain the final form of the response function

$S(Q_{0},q)=\dfrac{M_{N}^{2}\,Q_{0}}{\pi q}\dfrac{1}{1-e^{-z}}\Theta(\mu-
E_{-}).$ (45)
$E^{-}=\sqrt{\bigg{(}m_{N}+\frac{Q_{0}}{2}\bigg{)}^{2}+\bigg{(}\frac{\sqrt{q^{2}-Q_{0}^{2}}}{2}-\frac{m_{N}Q_{0}}{\sqrt{q^{2}-Q_{0}^{2}}}\bigg{)}^{2}}-\left(m_{N}+\frac{Q_{0}}{2}\right)\,.$
(46)

which can be simplified further

$E^{-}\approx\sqrt{m_{N}^{2}+\bigg{(}\frac{m_{N}Q_{0}}{\sqrt{q^{2}-Q_{0}^{2}}}\bigg{)}^{2}}-m_{N}\,.$
(47)

To perform the integral over $p_{a}$, we make the following change of
variables

$qdq=E_{a}E_{\gamma}d\cos\theta,\qquad
dQ_{0}=dE_{a}\qquad\xrightarrow{}\quad\int\frac{d^{3}p_{a}}{(2\pi)^{3}2E_{a}}=\int\frac{dQ_{0}dqq}{(2\pi)^{2}2}\frac{1}{E_{\gamma}}\,,$
(48)

and the integral becomes

$\displaystyle Q^{N\gamma\to Na}_{\rm
D}=\frac{g_{\gamma}\kappa^{2}}{f_{a}^{2}}\int\dfrac{d^{3}p_{\gamma}}{(2\pi)^{3}}\dfrac{f_{\gamma}}{E_{\gamma}}\int\frac{dQ_{0}\,dq}{(2\pi)^{2}2E_{\gamma}}\dfrac{M_{N}^{2}\,Q_{0}}{\pi}\times\bigg{[}2E_{\gamma}E_{a}(p_{\gamma}\cdot
p_{a})-(p_{\gamma}\cdot
p_{a})^{2}\bigg{]}\dfrac{E_{\gamma}-Q_{0}}{1-e^{-Q_{0}/T}}\Theta(\mu-
E_{-})\,.$ (49)

The expression between the parenthesis in Eq.(49) can be rewritten

$\bigg{[}2E_{\gamma}E_{a}(p_{\gamma}\cdot p_{a})-(p_{\gamma}\cdot
p_{a})^{2}\bigg{]}=\bigg{[}(E_{\gamma}E_{a})^{2}-(\vec{p}_{\gamma}\cdot\vec{p}_{a})^{2}\bigg{]}\approx\frac{(q^{2}-Q_{0}^{2})}{4}\bigg{(}-(q^{2}-Q_{0}^{2})+4E_{\gamma}(E_{\gamma}-Q_{0})\bigg{)}\,.$
(50)

Figure 4: Comparison plots for the DFSZ model.

To perform the integral, we need to define the boundaries of integration.
Kinematics impose the condition that $q$ variable varies from
$q\in[Q_{0},2E_{\gamma}-Q_{0}]$. The Theta function in Eq.(49) combined with
the expression for $E_{-}$ in Eq.(47) further impose the lower bound on the
$q$ region of integration as $b|Q_{0}|$. In the end, we get

$q\in\bigg{[}b|Q_{0}|,2E_{\gamma}-Q_{0}\bigg{]}\;,\qquad
a\equiv(\mu_{N}/M_{N}+1)^{2}-1\,,\qquad b\equiv\sqrt{\frac{a+1}{a}}\,.$ (51)

Taking the typical value $\mu_{N}=300$ MeV, we find $b\approx 1.56$. We will
use those two values as a range of uncertainty for the degenerate density of
the proto-NS. This induces that the $Q_{0}$ integral should be performed
between

$Q_{0}\in\bigg{[}-\frac{2}{b-1}E_{\gamma},\frac{2}{b+1}E_{\gamma}\bigg{]}\,.$
(52)

We can further separate the parameter space into the region of

$Q_{0}\in\bigg{[}-\frac{2}{b-1}E_{\gamma},0\bigg{]}\,,\qquad
Q_{0}\in\bigg{[}0,\frac{2}{b+1}E_{\gamma}\bigg{]}\,.$ (53)

The integral over $q$ can be trivially performed and by finishing the
integration numerically, we obtain

$\displaystyle\frac{Q_{N\gamma\to Na}^{\rm
WZW,D}}{10^{34}\text{erg/s/cm}^{3}}\approx
1.7\;C_{A}^{2}\;g_{40}^{4}T_{40}^{9}\bigg{(}\frac{10^{9}~{}\text{GeV}}{f_{a}}\bigg{)}^{2}\,.$
(54)

##  III Emissivity for different axion models

In generic axion and ALPs, models, the the ALP can have a myriad more
couplings, namely directly to quarks and leptons, making the whole coupling
structure quite model-dependent. To demonstrate the model dependence, we
present the comparison between the pion conversion, the axion bremsstrahlung
and the WZW photo-production in Fig.4 using the following parameterization[60]

$C_{ap}=-0.47+0.29\cos^{2}\beta-0.13\sin^{2}\beta\,,\qquad
C_{an}=-0.02+0.29\sin^{2}\beta-0.13\cos^{2}\beta\,.$ (55)

We also find that only when $\beta\ll 1$, the WZW contribution subleading.
When $\beta$ is larger, the enhancement of the emissivity is order
$\mathcal{O}(2-3)$.


Sensitivity Analysis for Random Measurement Error using Regression Calibration
and Simulation-Extrapolation

Linda Nab and Rolf H.H. Groenwold

Correspondence to Linda Nab MSc, Department of Clinical Epidemiology, Leiden
University Medical Center, Postzone C7-P, P.O. Box 9600, 2300 RC Leiden, the
Netherlands (email<EMAIL_ADDRESS>
Author affliations: Department of Clinical Epidemiology, Leiden University
Medical Center, Leiden, the Netherlands (Linda Nab, and Rolf H.H. Groenwold);
and Department of Biomedical Data Sciences, Leiden University Medical Center,
Leiden, the Netherlands (Rolf H.H. Groenwold).

This study was supported by grants from the Netherlands Organization for
Scientific Research (ZonMW-Vidi project 917.16.430) and Leiden University
Medical Center.

Conflicts of interest: none declared.

Running head: Sensitivity analysis for random measurement error

Data and code availability: The data and code used for the simulation study
have been made publicly and can be accessed via
www.github.com/LindaNab/simexvsmecor.

###### Abstract

(200/200)
Objective. Sensitivity analysis for measurement error can be applied in the
absence of validation data by means of regression calibration and simulation-
extrapolation. These have not been compared for this purpose.
Study Design and Setting. A simulation study was conducted comparing the
performance of regression calibration and simulation-extrapolation in a
multivariable model. The performance of the two methods was evaluated in terms
of bias, mean squared error (MSE) and confidence interval coverage, for
ranging reliability of the error-prone measurement (0.2-0.9), sample size
(125-1,000), number of replicates (2-10), and R-squared (0.03-0.75). It was
assumed that no validation data were available about the error-free measures,
while measurement error variance was correctly estimated.
Results. In various scenarios, regression calibration was unbiased while
simulation-extrapolation was biased: median bias was $1.4\%$ (interquartile
range (IQR): $0.8$;$2\%$), and $-12.8\%$ (IQR: $-13.2$;$-11.0\%$),
respectively. A small gain in efficiency was observed for simulation-
extrapolation (median MSE: $0.005$, IQR: $0.004$;$0.006$) versus regression
calibration (median MSE: $0.006$, IQR: $0.004$;$0.007$). Confidence interval
coverage was at the nominal level of 95% for regression calibration, and
smaller than $95\%$ for simulation-extrapolation (median coverage: $92\%$,
IQR: $85$;$94\%$).
Conclusion. In the absence of validation data, the use of regression
calibration is recommended for sensitivity analysis for measurement error.

Keywords: random measurement error, sensitivity analysis, quantitative bias
analysis, regression calibration, simulation-extrapolation

##  1 INTRODUCTION

Measurement error is common in biomedical research but often ignored [1, 2].
When ignored, measurement error can lead to considerable biases in exposure-
outcome associations [3]. Random measurement error in the exposure variable is
common in various domains of epidemiology [4, 5, 6]. Random measurement error
in an exposure variable introduces bias in the exposure-outcome association,
which is sometimes referred to as attenuation bias [7] or regression dilution
bias [8, 9].

Various methods for measurement error correction are available [10, 11, 12,
13, 14, 15, 16, 17, 18], yet application of these methods is rare in
biomedical research [3]. One possible barrier is the necessity of for instance
validation data. Validation data contain extra measurements taken to validate
error-prone measurements, which are often unavailable [5]. Validation data can
be used to estimate the measurement error model and its parameters, and
subsequently used for measurement error correction.

In the absence of validation data, regression calibration [19, 20] and
simulation-extrapolation [21], for instance, can be applied to correct for
random exposure measurement error. Both methods only require assumptions about
the variance of the random measurement error, for example based on literature
or expert knowledge. Regression calibration in the absence of validation data
is available in the R [22] package mecor for measurement error correction
[23], that implements the regression calibration described by Rosner et al.
[24]. Alternatively, simulation-extrapolation is easy to use due to its
implementation in the R package simex [25] and the simex procedure [26] in
Stata [27].

Simulation-extrapolation and regression calibration have been compared in
simulation studies for scenarios where replicate measures of the error-prone
exposure were available [5, 28]. The results of the studies by Perrier et al.
[5] and Batistatou et al. [28] were consistent and showed that, regression
calibration and simulation-extrapolation reduced bias compared to when no
measurement error correction was applied or when the replicate exposure
measures were pooled. It was also shown that application of simulation-
extrapolation generally produced more biased effect estimates than regression
calibration, especially when the reliability of the error-prone measure was
low.

Perrier et al. [5] and Batistatou et al. [28] only studied a limited number of
scenarios, e.g., large sample sizes, limited range of reliability of the
error-prone measure, and univariable models only. Further investigation is
needed of the performance of regression calibration and simulation-
extrapolation in more complex settings, as typically found in epidemiologic
research. Moreover, since the simulation studies by Perrier et al. and
Batistatou et al. focused on settings where replicate measures were available,
we aim to research how their results translate to settings were no replicate
measures, but only an estimate of the measurement error variance is available.
The quantification of the performance of the two methods in this broader range
of settings is used as the input for a framework guiding the application of
measurement error correction in the absence of validation data.

This paper is structured as follows. Section 2 reviews and applies regression
calibration and simulation-extrapolation by using two motivating examples of
error-prone systolic blood pressure measurements. In section 3, a simulation
study is described that aims to compare regression calibration and simulation-
extrapolation, and results from the simulation study are shown. Section 4
introduces a framework for conducting sensitivity analysis, also known as
quantitative bias analysis [29], for random measurement error by means of
regression calibration and simulation-extrapolation. We conclude with a
discussion of our results and recommendations in section 5.

##  2 MOTIVATING EXAMPLE

To demonstrate the use of regression calibration and simulation-extrapolation
to correct for random measurement error, we used an example about the
association between blood pressure (systolic blood pressure) and kidney
function (serum creatinine).

Example 1: Blood pressure and kidney function in adults
For the first example, we used data of the 2015-2016 cycle of the National
Health And Nutrition Examination Survey [30]. Given natural variation of blood
pressure within individuals, a single measure of blood pressure often does not
reflect the true level of blood pressure. Therefore, in epidemiological
studies, often multiple blood pressure measurements are obtained. In NHANES,
three consecutive systolic BP measurements were taken, after resting quietly
in a seated position for 5 minutes. The serum creatinine concentration was
determined once. Demographic information was collected using the family and
sample person demographics questionnaires in the home, by trained
interviewers. For this analysis, pregnant women, individuals receiving
dialysis in the past 12 months or indicating that that they had been told they
had weak/failing kidneys were excluded from analysis. Additionally, all
individuals with a systolic blood pressure greater than 180 mmHg and serum
creatinine greater than 200 $\mu$mol/L were excluded.

First, the association between systolic blood pressure and serum creatinine
was determined by only using the first systolic blood pressure measure. The
association was adjusted for age. It was found that an increase of 10 mmHg in
systolic blood pressure was associated with a 0.87 $\mu$mol/L (95% confidence
interval (CI): 0.51-1.22) increase in serum creatinine (Table 1). However, the
random measurement error in the first systolic blood pressure measure was not
considered here. In the NHANES data, the within individual variation of a
systolic blood pressure measure was on average 13.1 mmHg, which was obtained
by using the three consecutive systolic blood pressure measures. The within
individual variation of 13.1 mmHg was subsequently used to correct for the
measurement error in the first systolic blood pressure measure using
regression calibration and simulation-extrapolation, while adjusting for age.
An explanation of the two methods can be found in subsection 2.1. Using
regression calibration, we found that an increase of 10 mmHg in systolic blood
pressure was associated with an 0.93 $\mu$mol/L (CI: 0.54-1.32) increase in
serum creatinine (Table 1). Using simulation extrapolation, we found that an
increase of 10 mmHg in systolic blood pressure was associated with an 0.93
$\mu$mol/L (CI: 0.50-1.36) increase in serum creatinine (Table 1).

Example 2: Blood pressure and kidney function in pregnant women
For the second example, we used data of retrospective records of all women who
attended a tertiary maternity hospital pregnancy day assessment clinic over a
6-month period in 2014 in Australia [31]. Care always included serial, manual
blood pressure measurements every 30 min by registered midwives using aneroid
sphygmomanometers [31]. Serum creatinine and demographic data were obtained
using routinely collected data. One woman with a serum creatinine level lower
than 10 $\mu$mol/L was excluded from the analysis.

First, the association between systolic blood pressure and serum creatinine
was determined by only using the systolic blood pressure measure obtained
after 30 min. The association was adjusted for age. We found that an increase
of 10 mmHg in systolic blood pressure was associated with a 1.18 $\mu$mol/L
(CI: 0.14 - 2.23) increase in serum creatinine (Table 1). In this analysis,
the random measurement error in the single systolic blood pressure measurement
was not taken into account. Using the four consecutive blood pressure
measurements (obtained after 30, 60, 90 and 120 minutes), it was found that
the within individual variation of the systolic blood pressure measures was on
average 48.3 mmHg. The within individual variation of 48.3 mmHg was
subsequently used to correct for the measurement error in the single systolic
blood pressure measurement using regression calibration and simulation-
extrapolation, while adjusting for age. Using regression calibration, we found
that an increase of 10 mmHg was associated with a 2.04 $\mu$mol/L (CI:
0.23-3.91) increase in serum creatinine (Table 1). Using simulation-
extrapolation, we found that an increase of 10 mmHg was associated with a 1.67
$\mu$mol/L (CI: 0.01-3.25) increase in serum creatinine (Table 1).

Table 1: Effect estimates (95% confidence intervals) of the association
between blood pressure (systolic blood pressure, per 10 mmHg) and kidney
function (serum creatinine, $\mu$mol/L). The uncorrected effect estimates are
obtained using the first blood pressure measure only, the corrected estimates
are obtained by using the two consecutive measurements (example 1) and the
three consecutive measurments (example 2).

Study | Uncorrected | Regression Calibration | Simulation-extrapolation
---|---|---|---
Example 1 | 0.87 (0.51-1.22) | 0.93 (0.54-1.32) | 0.93 (0.50-1.36)
Example 2 | 1.18 (0.14-2.23) | 2.04 (0.23-3.91) | 1.67 (0.01 - 3.25)

  * •

Estimates were obtained from the national health and nutrition examination
survey (example 1, reliability of the error-prone blood pressure measurement:
0.95) [30] and the pregnancy day and assessment clinic study (example 2,
reliability of the error-prone blood pressure measurement: 0.6) [31]

###  2.1 Review of regression calibration and simulation-extrapolation

In examples 1 and 2, the measurement error in the blood pressure measures was
corrected by application of regression calibration and simulation-
extrapolation. Regression calibration started by estimating the uncorrected
effect of the error-prone first blood pressure measure on serum creatinine
given age. Subsequently, the uncorrected estimate was multiplied by the
estimated variance of the first error-prone blood pressure measurement given
age, divided by the estimated variance of the first error-prone blood pressure
measurement given age minus the measurement error variance. In examples 1 and
2, the measurement error variance was estimated using the replicate
measurements, by estimating the within individual variance and averaging over
all individuals. Alternatively, the measurement error variance could be
informed by e.g. external data or expert knowledge.

Simulation-extrapolation consists of two steps. In the simulation step, extra
measurement error was added to the error-prone first blood pressure
measurement of size 0.5, 1, 1.5, 2 times the measurement error variance. The
measurement error variance was again estimated using the replicate
measurements. Using these simulated blood pressure measurements with extra
added measurement error, the effect of blood pressure on serum creatinine
given age was estimated. This was repeated 100 times for each size of the
measurement error variance and the newly obtained estimates were averaged.
Then, in the extrapolation step, a model (e.g., linear, quadratic) was fitted
trough the effect estimates for the varying sizes of the measurement error.
The corrected effect estimate was then obtained by extrapolating the fitted
model to the situation where the measurement error was equal to 0. For a
visualisation of simulation-extrapolation, see e.g. [3].

Ninety-five percent CI’s for the corrected estimates were constructed by
bootstrap resampling.

##  3 SIMULATION STUDY

We aimed to extend the simulation studies conducted by Perrier et al. [5] and
Batistatou et al. [28], by investigating the relative performance of
regression calibration and simulation-extrapolation in scenarios different
from those studied in the two former studies. The relative performance was
studied in terms of bias, mean squared error, and confidence interval coverage
of the true effect. A general description and motivation of the scenarios
studied in our simulation study is provided in the following, the specific
parameters set in our simulation study are explained in the subsequent
subsection.

Perrier et al. and Batistatou et al. assumed relatively large sample sizes
(i.e., 3000 and 1000, respectively), only four different values for the
reliability of the error-prone exposure (i.e., 0.2 and 0.6 in the study by
Perrier et al. and 0.2, 0.5 and 0.8 in the study by Batistatou et al.) and a
small coefficient of determination for the exposure-outcome model (i.e., 0.004
and 0.0625, respectively). In addition, Perrier et al. studied the effect of
increasing the number of replicate measures on the performance of regression
calibration and simulation-extrapolation. Yet, the replicate measures were
pooled in the former study, thereby decreasing measurement error bias.
Moreover, Perrier et al. and Batistatou et al. only examined models with a
single independent variable. Our simulation study focused on relatively large
reliability of the error-prone measurement (i.e., greater or equal to 0.625),
motivated by the reliability of the blood pressure measures of our two
motivating examples (i.e., 0.95 and 0.6 in example 1 and example 2,
respectively). Additionally, the relative performance of regression
calibration vs simulation-extrapolation was studied for small sample sizes
(i.e., smaller or equal to 1000). In addition, the effect of a change in the
coefficient of determination of the outcome model was tested. Furthermore,
increasing the number of replicate measurements available was studied, without
having the advantage of pooling the replicate measurements in our analysis.
Lastly, multivariable models were studied.

###  3.1 Data generating mechanism

Inspired by the data of our example of blood pressure and kidney function in
pregnant women [31], we assume the following data generating mechanisms for
age, blood pressure (BP), error-prone blood pressure (BP∗) and creatinine:

$\displaystyle\textrm{Age}\sim\mathcal{N}(32,25),\quad\textrm{BP}|\textrm{Age}\sim\mathcal{N}(120+\gamma\textrm{Age},50),\quad\textrm{BP}^{*}\sim\mathcal{N}(\textrm{BP},\tau^{2}),$
$\displaystyle\textrm{and}\quad\textrm{Creatinine}|\textrm{BP,
Age}\sim\mathcal{N}(30+0.2\textrm{BP}+0.2\textrm{Age},\sigma^{2}).$

The above defined data generating mechanism define that the error-prone blood
pressure (BP∗) has random measurement error with measurement error variance
equal to $\tau^{2}$. In our simulation study, a ‘base scenario’ was assumed
and in the consecutive scenarios studied, we changed one of the four
parameters in the data generating mechanisms (i.e., $\gamma,\tau^{2}$ or
$\sigma^{2}$), the number of observations (i.e., $n$), or the number of
replicate measures (i.e., $k$) (see Table 2). The parameters settings of the
base scenario were inspired by our example of blood pressure and kidney
function in pregnant women [31]. In the base scenario, sample size was 500,
$\gamma=0$, $\tau^{2}=30$ and $\sigma^{2}=100$ (Table 2). Further, we assumed
that three replicate measures of the error-prone blood pressure measure were
obtained in all individuals. From the parameter settings in the base scenario,
it follows that the reliability of the error-prone measure is 0.625. The
reliability is defined as the variance in the error-free blood pressure
divided by the variance of the error-prone blood pressure. Further, in the
base scenario, the R-squared of the outcome model is 0.03, and the attenuation
due to measurement error of the effect of the error-prone blood pressure on
creatinine (given age) is equal to the reliability, i.e., 0.625.

In subsequent scenarios (see Table 2), we increased the reliability of the
error–prone measure, by decreasing the variance of the error-prone measure
(i.e., $\tau^{2}$). For completeness, 3 scenarios with low reliability were
added. Subsequently, we changed the number of observations (see under ‘Sample
size’ in Table 2). Next, we changed the number of replicate measures to study
the effect of changing the precision of the estimation of the measurement
error variance (see under ‘Number of replicates’ in Table 2). Next, we
decreased the residual error variance in the outcome model, $\sigma^{2}$ (see
under ‘R-squared’ in Table 2). Decreasing the residual error in the outcome
model, increases the R-squared of the outcome model. Lastly, we increased the
reliability of the error-prone variable, by introducing a dependency between
blood pressure and age, by changing $\gamma$ (see under ‘Covariate dependency’
in Table 2). By introducing an effect of age on blood pressure, the total
variance of the error-free blood pressure increases. Consequently, the extra
variability in the error-prone blood pressure measure due to measurement error
is relatively smaller than in the base scenario. Hence, it seems as if the
error-prone variable is more reliable, though the attenuation due to
measurement error stays constant at a rate of 0.625. For each scenario, 5000
datasets were generated.

In each generated data set, the uncorrected effect was estimated using the
first replicate measurement only. Subsequently, the corrected effect was
estimated by application of regression calibration and simulation-
extrapolation using the R package mecor [23] and simex [25], respectively. The
measurement error variance was estimated using the replicate measures. Ninety-
five percent CI’s of the corrected effects were constructed using bootstrap
resampling. Performance of the three different analyses were evaluated in
terms of bias, mean squared error (MSE), and the proportion of 95% CIs that
contained the true value of the estimand (coverage). Monte Carlo standard
errors (MCSE) were calculated for all performance measures [32], using the R
package rsimsum [33]. All code used for the simulation study is publicly
available via www.github.com/LindaNab/simexvsmecor.

Table 2: Simulation study settings

Scenarios | Parameters of Data Generating Mechanisma
---|---
| $\tau^{2}$ | $n$ | $k$ | $\sigma^{2}$ | $\gamma$
Base | 30 | 500 | 3 | 100 | 0
Reliabilityb | 200, | 500 | 3 | 100 | 0
| 100, |  |  |  |
| 50, 25, |  |  |  |
| 20, 15, |  |  |  |
| 10, 5 |  |  |  |
Sample Size | 30 | 125, | 3 | 100 | 0
|  | 250, |  |  |
|  | 1000, |  |  |
|  | 10,000 |  |  |
Number of | 30 | 500 | 2, 5, | 100 | 0
Replicates |  |  | 10 |  |
R-squaredc | 30 | 500 | 3 | 20, 5, 1 | 0
Covariate | 30 | 500 | 3 | 100 | 1, 4, 8
Dependencyc |  |  |  |  |

  * a

$\tau^{2}$: measurement error variance of the error-prone blood pressure
measurement; $n$: number of observations in the main study; $k$: number of
replicate error-prone measurements; $\sigma^{2}$: residual variance of the
outcome model; $\gamma$: association between blood pressure and age. The
attenuation in the effect of blood pressure on creatinine due to measurement
error is equal to $50/(50+\tau^{2})$

  * b

Reliability is equal to $(25\gamma^{2}+50)/(25\gamma^{2}+50+\tau^{2})$

  * c

R-squared is equal to $1-\sigma^{2}/(0.4\times 50+10+\sigma^{2})$

  * d

The effect of blood pressure on creatinine when age is not included in the
model (crude model) is equal to $0.2+5\gamma/(25\gamma^{2}+50)$

###  3.2 Results

Figure 1 shows the percentage bias, MSE and confidence interval coverage for
varying values of the reliability of the error-prone measure. The uncorrected
analysis was biased for all values of the reliability, and the percentage bias
decreased when reliability increased. Regression calibration provided unbiased
results when reliability was greater or equal to $0.33$. Simulation-
extrapolation provided biased results when reliability was smaller than $0.8$.
MSE was lower for simulation-extrapolation than for the uncorrected and
regression calibration corrected analysis when reliability was equal to 0.2,
and similar to MSE of regression calibration otherwise. Coverage of the 95%
confidence intervals was at the nominal level for the regression calibration
corrected analysis when reliability was greater than or equal to $0.33$, and
for the simulation-extrapolation corrected analysis when reliability was
greater than or equal to $0.625$.
Figure 2 shows the percentage bias, MSE and confidence interval coverage for
varying samples sizes of the main study. A sample size of 125, and 250 only
increased percentage bias minimally compared to base scenario were sample size
was 500. MSE was greater for smaller sample sizes, and MSE of the uncorrected
analysis with a sample size of 125 was smaller than the regression calibration
and simulation-extrapolation corrected analysis (0.015 vs 0.026 and 0.019,
respectively, MCSE $<0.005$). Coverage was equal to the nominal level of 95%
for regression calibration for all sample sizes, and the uncorrected analysis
was undercovered with coverage levels ranging between 91% and 45% (MCSE
$<0.01$). Coverage of the 95% confidence intervals of the simulation-
extrapolation corrected analysis was close to the nominal level of 95% except
when sample size was 1000, in which case coverage was 90% (MCSE $0.004$). A
decline in confidence interval coverage for the simulation-extrapolation
corrected analysis for larger sample sizes was confirmed by the scenario where
sample size was 10,000, in which case coverage was 53% (MCSE 0.007)(not shown
in the plots in Figure 2).
Figure 3 shows that the number of replicates had no effect on the percentage
bias, MSE and confidence interval coverage for varying number of replicates of
the error-prone measure.
Figure 4 shows that R-squared had no effect on percentage bias, and only a
minor decrease in MSE was found for increasing levels of R-Squared. In
addition, Figure 4 shows that 95% confidence interval coverage was around the
nominal level for the regression calibration corrected analysis for all values
of the R-squared. However, for the uncorrected and the simulation-
extrapolation corrected analysis, confidence interval coverage decreased for
increasing values of R-Squared. For R-Squared equal to 0.75, confidence
interval coverage decreased to 15 % and 0 % (MCSE $<0.01$) for the simulation-
extrapolation corrected and the uncorrected analysis, respectively.
In the scenarios where a dependency between the covariate age and the exposure
error-free blood pressure was introduced by changing parameter $\gamma$ in the
data generating mechanism, the reliability of the error-prone measure was
respectively 0.71, 0.94 and 0.98. However, percentage bias, MSE and confidence
interval coverage of the uncorrected and corrected analyses were equal to base
scenario (the values in base scenario are shown in e.g. Figure 1).

---
Figure 1: Performance of regression calibration (RC), simulation-extrapolation
(simex) and the analysis ignoring measurement error for varying values of the
reliability of the error-prone measure in terms of A) percentage bias; B) mean
squared error and C) coverage. For all three performance measures, Monte Carlo
standard errors were smaller than 0.01 in all scenarios. The grey points
indicate the base scenario where reliability is assumed 0.625.
---
Figure 2: Performance of regression calibration (RC), simulation-extrapolation
(simex) and the analysis ignoring measurement error for varying sample sizes
of the error-prone measure in terms of A) percentage bias; B) mean squared
error and C) coverage. For all three performance measures, Monte Carlo
standard errors were $<0.01$ in all scenarios. The grey points indicate the
base scenario where sample size is assumed 500.
---
Figure 3: Performance of regression calibration (RC), simulation-extrapolation
(simex) and the analysis ignoring measurement error for varying number of
replicates of the error-prone measure in terms of A) percentage bias; B) mean
squared error and C) coverage. For all three performance measures, Monte Carlo
standard errors were smaller than 0.01 in all scenarios. The grey points
indicate the base scenario where the number of replicates is assumed 3.
---
Figure 4: Performance of regression calibration (RC), simulation-extrapolation
(simex) and the analysis ignoring measurement error for varying R-Squared of
the outcome model A) percentage bias; B) mean squared error and C) coverage.
For all three performance measures, Monte Carlo standard errors were smaller
than 0.01 in all scenarios. The grey points indicate the base scenario where
R-Squared is assumed 0.03.

##  4 SENSITIVITY ANALYSIS IN THE ABSENCE OF VALIDATION DATA

In the two examples introduced in section 2, replicate measurements of the
error-prone systolic blood pressure were available. Nevertheless, validation
data in the form of replicate measurements may not always be available. When
random measurement error in a covariate is suspected in the absence of such
validation data, sensitivity analysis could be conducted using regression
calibration or simulation-extrapolation. A general framework for conducting
sensitivity analysis for measurement error is described here, where we assume
that the input of the sensitivity analysis, i.e., the measurement error
variance and its uncertainty, are obtained from literature or expert
knowledge. Next, a distribution for the measurement error variance is assumed,
e.g., a uniform, triangular, or trapezidiol distribution [29]. Subsequently,
regression calibration or simulation-extrapolation are applied to the data for
measurement error correction, informed by the measurement error variance and
its distribution. Finally the results of the application of measurement error
correction are presented, and conclusions drawn about the sensitivity of the
results to measurement error.

###  4.1 Sensitivity analysis in the motivating example

Suppose that in the example of the relation between blood pressure and kidney
function in pregnant women discussed in section 2 (example 2), only the first
systolic blood pressure measurement was available. A 10 mmHg increase in
systolic blood pressure was associated with a 1.18 $\mu$mol/L (95% CI
0.14-2.23) increase in serum creatinine. Measurement error, however, was
suspected in the single systolic blood pressure measure and suppose the
sensitivity of the results to the measurement error was studied. Suppose it
was assumed that the variance of the measurement error in systolic blood
pressure is equal to 48 mmHg, with a minimum of 37 mmHg and a maximum of 59
mmHg. Additionally, suppose a triangular distribution was assumed for the
measurement error variance, meaning that most weight was put on 48 mmHg, and
the weight was gradually reduced until it reached the assumed minimum and
maximum level. The triangular distribution was sampled in accordance with Lash
et al. [29].
Figure 5 shows the results of the application of regression calibration and
simulation-extrapolation informed by the assumed triangular distribution. For
regression calibration, a clear pattern was obtained in the sensitivity
analysis plot shown in Figure 5. The corrected effect estimates increased for
larger values of the measurement error variance, with the effect estimates
ranging from 1.75 - 2.38, with a median of 2.03. In addition, the associated
lower limits of the confidence intervals consistently suggest an association
between blood pressure and creatinine. In comparison, simulation-extrapolation
did not show a clear pattern in the corrected effect estimates. The corrected
effect estimates ranged from 1.43 - 1.88, with a median of 1.70. Figure 5
shows that the sampling variability that is inherent to simulation-
extrapolation causes more variability in the effect estimates compared to the
variability due to measurement error. Nevertheless, the lower limits of the
associated confidence intervals again consistently suggest an association
between blood pressure and creatinine.

---
Figure 5: Sensitivity analysis for the association between blood pressure and
kidney function in pregnant women (example 2) by application of regression
calibration (panel A) and simulation-extrapolation (panel B). The uncorrected
association and 95% confidence interval is depicted with a diamond and a solid
black line, the measurement error corrected associations and 95% confidence
intervals are depicted with a square and a solid gray line. The distribution
of the measurement error variance is triangular.

##  5 DISCUSSION

This paper compared regression calibration and simulation-extrapolation for
conducting a sensitivity analysis for random measurement error in an exposure
variable. A simulation study showed that with correct assumptions about the
measurement error variance, the regression calibration corrected analysis was
generally unbiased. However, the uncorrected and simulation-extrapolation
corrected analysis were generally biased, with higher bias for lower
reliability of the error-prone exposure. The uncorrected analysis showed worse
efficiency in terms of mean squared error over the two corrected analyses,
expect for small sample sizes (i.e., 125). Compared to regression calibration,
the simulation-extrapolation corrected analysis showed similar efficiency in
terms of mean squared error, except for low reliability (i.e., 0.2 and a
sample size of 500) and small sample size (i.e., 250 and 125, and a
reliability of 0.625). In addition, regression calibration showed confidence
interval coverage of levels close to 95% in all scenarios, while the
confidence intervals of the simulation-extrapolation corrected analysis and
uncorrected analysis were often undercovered. What is more, an increment of
the R-squared value of the outcome model showed a large decrease in confidence
interval coverage of the true effect for the uncorrected and simulation-
extrapolation corrected analysis.

The results of our simulation study were in line with the results of two
previous simulation studies: the corrected analyses showed lower percentages
bias compared to the uncorrected analysis and the simulation-extrapolation
corrected analysis showed higher percentage bias compared to regression
calibration [5, 28]. However, important differences were observed. First,
simulation-extrapolation showed a small gain in efficiency over regression
calibration in some settings, which was not found in the previous simulation
studies. The sample sizes for which this gain in efficiency for simulation-
extrapolation was observed (i.e., 125, 250 and 500) were smaller than those
assumed by Perrier et al. [5] and Batistatou et al. [28] (i.e., 3000 and 1000,
respectively), which may explain the found difference. Second, our simulation
showed no effect of the number of replicates on bias. While the simulation
study by Perrier et al. showed that an increasing number of replicates reduced
bias in the corrected analyses [5]. This difference is explained by the fact
that in the study by Perrier et al., the replicate measures were pooled before
applying measurement error correction. By pooling the replicate measures with
random measurement error, the measurement error variance is reduced.
Therefore, bias decreased in the study by Perrier et al. with the availability
of more replicate measures. This effect however is solely due to pooling the
replicates measures and not due to a more precise estimate of the measurement
error variance, as was shown by our results. Of note, when replicate
measurements are obtained, pooling the replicate measurements is always
preferred above only using the replicate measurements to estimate the
measurement error variance.

Our simulation study showed that percentage bias in the uncorrected analysis
was equal to 1 minus the reliability of the error-prone measure times 100, in
line with theory [8, 9]. The reliability of an error-prone measure equals the
variance of the error-free measure divided by the variance of the error-prone
measure. For example, in Figure 1, bias in the uncorrected analysis was equal
to 80% for a reliability equal to 0.2. The uncorrected effect estimate is
equal to 0.2 times the estimand 0.2, i.e., 0.04. From that, it follows that
the bias is equal to $0.2-0.04=0.16$, which is 80% of 0.2. It is, however,
important to note that the percentage bias is not equal to 1 minus the
reliability of the error-prone measure when the total variance of the error-
free measure depends on a covariate that is also included in the outcome
model. For example, in our simulation study, the association between
creatinine and systolic blood pressure given age was estimated. When a
dependency between systolic blood pressure and age was introduced, the
reliability increased to a maximum of 0.98 while the percentage bias in the
uncorrected analysis was constant at 62.5%. A formula for the attenuation in
the effect estimate due to random measurement error in multivariable models
can be found in e.g. [13].

In our simulation study, the measurement error variance used to correct for
the measurement error was estimated using replicate measures. However, we
assumed that these replicate measures were solely available to estimate the
measurement error variance, to mimic a setting where such validation data is
not available but with unbiased estimation of the measurement error variance.
In future studies, this work could be extended to settings where the
measurement error is estimated with bias, and to settings where the
measurement error model is misspecified. Also systematic measurement error,
which was not the topic of this study, could be considered in future work.

In the example presented in section 4, the five steps of a sensitivity
analysis for exposure measurement error were described: 1) quantify the
measurement error variance and its uncertainty; 2) specify the distribution of
the measurement error variance; 3) perform measurement error correction by
means of regression calibration or simulation-extrapolation; 4) visualise the
results, and 5) draw conclusions. A sensitivity analysis using regression
calibration showed that the higher the measurement error variance, the more
the corrected effect estimate departs from the null, which is in line with the
literature [8, 9]. In the sensitivity analysis using simulation-extrapolation,
the variability in the corrected effect estimates due to the sampling
variability inherent to simulation-extrapolation exceeded the variability in
the corrected effect estimates due to the assumed measurement error variance.
Since regression calibration was generally shown unbiased and simulation-
extrapolation was generally shown biased, the application of regression
calibration for conducting sensitivity analyses is recommended. Moreover, the
advantage in efficiency of simulation-extrapolation over regression
calibration in some settings does play a less important role in sensitivity
analysis than its disadvantage in performance due to bias.

In conclusion, in the absence of validation data, regression calibration and
simulation-extrapolation are suitable for conducting sensitivity analyses for
random measurement error. Particularly, regression calibration is generally
preferred over simulation-extrapolation due to its unbiasedness in most
settings. To guide sensitivity analysis for random measurement error, a
sensitivity analysis framework was provided.

## References

  * [1] TB Brakenhoff, M Mitroiu, RH Keogh, KGM Moons, RHH Groenwold, and M van Smeden.  Measurement error is often neglected in medical literature: A systematic review.  J Clin Epidemiol, 98:89–97, 2018.
  * [2] PA Shaw, V Deffner, RH Keogh, JA Tooze, KW Dodd, H Küchenhoff, V Kipnis, and LS Freedman.  Epidemiologic analyses with error-prone exposures: review of current practice and recommendations.  Ann Epidemiol, 28(11):821–828, 2018.
  * [3] RH Keogh, PA Shaw, P Gustafson, RJ Carroll, V Deffner, KW Dodd, H Küchenhoff, JA Tooze, MP Wallace, V Kipnis, and LS Freedman.  STRATOS guidance document on measurement error and misclassification of variables in observational epidemiology: Part 1—Basic theory and simple methods of adjustment.  Stat Med, 39(16):2197–2231, 2020.
  * [4] JS Frazer, GE Barnes, V Woodcock, E Flanagan, T Littlewood, RJ Stevens, S Fleming, and HF Ashdown.  Variability in body temperature in healthy adults and in patients receiving chemotherapy: prospective observational cohort study.  J Med Eng Technol, 43(5):323–333, 2019.
  * [5] F Perrier, L Giorgis-Allemand, R Slama, and C Philippat.  Within-subject pooling of biological samples to reduce exposure misclassification in biomarker-based studies.  Epidemiology, 27(3):378–388, 2016.
  * [6] B Brunekreef, D Noy, and P Clausing.  Variability of exposure measurements in environmental epidemiology.  Am J Epidemiol, 125(5):892–898, 1987.
  * [7] C. Spearman.  The proof and measurement of association between two things.  Am J Psychol, 15(1):72–101, 1 1904.
  * [8] C Frost and S Thompson.  Correcting for regression dilution bias: comparison of methods for a single predictor variable.  J Roy Stat Soc, 163(2):173–189, 2000.
  * [9] Jennifer A. Hutcheon, Arnaud Chiolero, and James A. Hanley.  Random measurement error and regression dilution bias.  BMJ, 340:c2289, 2010.
  * [10] BG Armstrong.  Measurement error in the generalised linear model.  Commun stat-simul c, 14(3):529–544, 1 1985.
  * [11] JW Bartlett and RH Keogh.  Bayesian correction for covariate measurement error: A frequentist evaluation and comparison with regression calibration.  Stat Methods Med Res, 27(6):1695–1708, 2018.
  * [12] JP Buonaccorsi.  Measurement error: Models, methods, and applications.  Chapman & Hall/CRC, Boca Raton, FL, 2010.
  * [13] RJ Carroll, D Ruppert, LA Stefanski, and CM Crainiceanu.  Measurement error in nonlinear models: A modern perspective.  Chapman & Hall/CRC, Boca Raton, FL, 2nd edition, 2006.
  * [14] SR Cole, H Chu, and S Greenland.  Multiple-imputation for measurement-error correction.  Int J Epidemiol, 35(4):1074–1081, 2006.
  * [15] JR Cook and LA Stefanski.  Simulation-extrapolation estimation in parametric measurement error models.  J Am Stat Assoc, 89(428):1314–1328, 1994.
  * [16] P Gustafson.  Measurement error and misclassification in statistics and epidemiology: Impacts and Bayesian adjustments. Chapman & Hall/CRC, Boca Raton, FL, 2004.
  * [17] WA Fuller.  Measurement Error Models.  John Wiley & Sons, New York, NY, 1987.
  * [18] RH Keogh and IR White.  A toolkit for measurement error correction, with a focus on nutritional epidemiology.  Stat Med, 33:2137–2155, 2014.
  * [19] LJ Gleser.  Measurement error models.  Chemometr Intell Lab Syst, 10:45–57, 1991.
  * [20] RJ Carroll and LA Stefanski.  Approximate quasi-likelihood estimation in models with surrogate predictors.  J Am Stat Assoc, 85(411):652–663, 1990.
  * [21] LA Stefanski and JR Cook.  Simulation-extrapolation: The measurement error jackknife.  J Am Stat Assoc, 90(432):1247, 12 1995.
  * [22] R Core Team.  R: A Language and Environment for Statistical Computing, 2020.
  * [23] L Nab, M van Smeden, RH Keogh, and RHH Groenwold.  mecor: An R package for measurement error correction in linear regression models with a continuous outcome, 2021.
  * [24] B Rosner, D Spiegelman, and WC Willett.  Correction of logistic regression relative risk estimates and confidence intervals for random within-person measurement error.  Am J Epidemiol, 136(11):1400–1413, 1992.
  * [25] W Lederer and H Küchenhoff.  A short introduction to the SIMEX and MCSIMEX.  R news, 6(4):26–31, 1997.
  * [26] James W. Hardin, Henrik Schmiediche, and Raymond J. Carroll.  The regression-calibration method for fitting generalized linear models with additive measurement error.  The Stata Journal, 3(4):361–372, 2003.
  * [27] StataCorp.  Stata Statistical Software: Release 16, 2019.
  * [28] Evridiki Batistatou and Roseanne McNamee.  Performance of bias-correction methods for exposure measurement error using repeated measurements with and without missing data.  Stat Med, 31(28):3467–3480, 2012.
  * [29] TL Lash, MP Fox, and AK Fink.  Applying quantitative bias analysis to epidemiologic data.  Springer, New York, NY, 2009.
  * [30] Centers for Disease Control and Prevention (CDC). National Center for Health Statistics (NCHS). National Health and Nutrition Examination Survey Data. Hyattsville, MD: U.S. Department of Health and Human Services, Centers for Disease Control and Prevention, [2015-2016], [https://wwwn.cdc.gov/nchs/nhanes/continuousnhanes/default.aspx?BeginYear=2015].
  * [31] EA McCarthy, TA Carins, Y Hannigan, N Bardien, A Shub, and SP Walker.  Effectiveness and safety of 1 vs 4 h blood pressure profile with clinical and laboratory assessment for the exclusion of gestational hypertension and pre-eclampsia: a retrospective study in a university affiliated maternity hospital.  BMJ Open, 5(11):e009492–e009492, 2015.
  * [32] TP Morris, IR White, and MJ Crowther.  Using simulation studies to evaluate statistical methods.  Stat Med, 38(11):2074–2102, 2019.
  * [33] A Gasparini.  rsimsum: Summarise results from Monte Carlo simulation studies.  JOSS, 3(26):739, 6 2018.