# On the impact of the Migdal effect in reactor CE$\nu$NS experiments

M. Atzori Corona<EMAIL_ADDRESS>Dipartimento di Fisica,
Università degli Studi di Cagliari, Complesso Universitario di Monserrato -
S.P. per Sestu Km 0.700, 09042 Monserrato (Cagliari), Italy  Istituto
Nazionale di Fisica Nucleare (INFN), Sezione di Cagliari, Complesso
Universitario di Monserrato - S.P. per Sestu Km 0.700, 09042 Monserrato
(Cagliari), Italy  M. Cadeddu<EMAIL_ADDRESS>Istituto Nazionale di
Fisica Nucleare (INFN), Sezione di Cagliari, Complesso Universitario di
Monserrato - S.P. per Sestu Km 0.700, 09042 Monserrato (Cagliari), Italy  N.
Cargioli<EMAIL_ADDRESS>Dipartimento di Fisica, Università degli
Studi di Cagliari, Complesso Universitario di Monserrato - S.P. per Sestu Km
0.700, 09042 Monserrato (Cagliari), Italy  Istituto Nazionale di Fisica
Nucleare (INFN), Sezione di Cagliari, Complesso Universitario di Monserrato -
S.P. per Sestu Km 0.700, 09042 Monserrato (Cagliari), Italy  F. Dordei
<EMAIL_ADDRESS>Istituto Nazionale di Fisica Nucleare (INFN), Sezione
di Cagliari, Complesso Universitario di Monserrato - S.P. per Sestu Km 0.700,
09042 Monserrato (Cagliari), Italy  C. Giunti<EMAIL_ADDRESS>Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Torino, Via P. Giuria
1, I–10125 Torino, Italy

###### Abstract

The search for coherent elastic neutrino nucleus scattering (CE$\nu$NS) using
reactor antineutrinos represents a formidable experimental challenge, recently
boosted by the observation of such a process at the Dresden-II reactor site
using a germanium detector. This observation relies on an unexpected
enhancement at low energies of the measured quenching factor with respect to
the theoretical Lindhard model prediction, which implies an extra observable
ionization signal produced after the nuclear recoil. A possible explanation
for this additional contribution could be provided by the so-called Migdal
effect, which however has never been observed. Here, we study in detail the
impact of the Migdal contribution to the standard CE$\nu$NS signal calculated
with the Lindhard quenching factor, finding that the former is completely
negligible for observed energies below $\sim 0.3\,\mathrm{keV}$ where the
signal is detectable, and thus unable to provide any contribution to CE$\nu$NS
searches in this energy regime. To this purpose, we compare different
formalisms used to describe the Migdal effect that intriguingly show a perfect
agreement, making our findings robust.

##  I Introduction

Coherent elastic neutrino-nucleus scattering (CE$\nu$NS) is a pure weak
neutral current low-energy process predicted by Freedman in 1973 [1] where the
neutrino interacts with the nucleus as a whole. Namely, all the nucleons
within the nucleus respond coherently to the neutrino interaction, leading to
a higher cross section than other low-energy processes involving neutrinos.
The remarkable observations of CE$\nu$NS from the COHERENT Collaboration [2,
3, 4] have opened a new era to test our knowledge of the Standard Model (SM)
of particle physics and has posed an exciting technological challenge to
develop innovative detectors capable of spotting the extremely tiny nuclear
recoils produced as a single outcome of the interaction [5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. These experiments play a crucial role
in advancing our knowledge of neutrino interactions and for their implications
for fundamental physics [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
34, 35, 36].

CE$\nu$NS requires a high neutrino flux to produce signal events above
experimental backgrounds. Among the different neutrino sources available, in
this study, we focus on CE$\nu$NS produced from reactor antineutrinos and
observed with germanium detectors. There are three main germanium detectors
currently operating, namely NCC-1701 [13, 14] (also referred to as Dresden-
II), CONUS [37] and $\nu$GEN [12], located 10.39 m, 17.1 m and 11 m away from
$2.96\;\rm{GW}_{\rm{th}}$, $3.9\;\rm{GW}_{\rm{th}}$, $3.1\;\rm{GW}_{\rm{th}}$
commercial reactors, respectively. Interestingly, the Ricochet experiment at
the ILL site, 8.8 m away from the core of the $58.3\;\rm{MW}_{\rm{th}}$
research nuclear reactor, is also aiming to measure CE$\nu$NS down to the
sub-100 eV nuclear energy recoil regime [38].

In particular, the recent first observation of CE$\nu$NS at the Dresden-II
reactor [14] has gained a lot of attention due to the broad impact of such a
result on current and future CE$\nu$NS searches, and the physics that can be
extracted within the SM and beyond [39, 40, 41, 42]. As a matter of fact, this
measurement is highly affected by the knowledge of the germanium quenching
factor (QF) at low nuclear recoil energies. The QF quantifies the reduction of
the ionization yield produced by a nuclear recoil with respect to an electron
recoil of the same energy. Indeed, the CE$\nu$NS observation by Dresden-II
depends crucially on the two new QF measurements reported in Ref. [43]. They
have been obtained from photo-neutron source measurements, so-called YBe, and
from iron-filtered monochromatic neutrons, so-called Fef [43]. However, these
two QF determinations are in contrast with and significantly higher than the
standard Lindhard prediction with the parameter $k=0.157$ [44] and other
independent experimental measurements (see e.g. Ref. [45] for a recent
measurement and the summary plot in Fig. 1). Moreover, CONUS data disfavours
quenching parameters above $k=0.27$ [9] and a recent low-energy determination
of the QF finds a good agreement with the Lindhard theory with a parameter
$k=0.162\pm 0.004$ (stat+sys) [46].

Among the possible solutions to explain the increase of the QF at very low
energies, it has been proposed that the Lindhard model might be not sufficient
to encode the full behavior of the QF in the low-energy regime. In particular,
it has been suggested [43] that the Migdal effect may play a crucial role.
This yet unobserved process111On the same day our manuscript was made public
on arXiv, another one appeared on this topic [47]. In this work, the authors
reported the first direct search for the Migdal effect in liquid xenon using
nuclear recoils produced by tagged neutron scatters. If the suppression of the
Migdal effect with respect to the theoretical prediction found in Ref. [47]
happens also in germanium, the explanation of the high QF found in Ref. [43]
through the Migdal effect is even less likely than that found in this work.
has been first proposed in 1941 [48] and recently taken into serious
consideration in the context of dark matter searches [49, 50, 51, 52, 53, 54,
55, 56] and neutrino physics [57, 58]. The Midgal effect might happen after a
nuclear recoil is induced by a neutral particle, i.e., a neutron, a neutrino,
or dark matter, due to the displacement between the recoiling nucleus and the
electronic cloud. Indeed, the atomic electrons do not immediately follow the
motion of the recoiling nucleus. Thus, in order to restore equilibrium, extra
ionization can be injected into the detector. In Refs. [59, 60] there have
been attempts to include the Migdal effect by means of some phenomenological
parameters originally introduced to go beyond the standard Lindhard theory
[61]. However, the approach exploited in these studies neglects the
microscopic physics of the Migdal process that can instead be accounted for
using existing formalisms.

The goal of this paper is to perform a detailed characterization of the Migdal
effect in a germanium detector searching for CE$\nu$NS at a reactor site. To
do so, we use different approaches to describe the microscopic physics of the
phenomenon to quantify the robustness of our findings. Finally, we discuss to
which extent the Migdal effect plays a role in the observation of CE$\nu$NS at
the Dresden-II reactor and thus, whether it can be claimed as an explanation
of the anomalous enhancement of the QF at low energies.

##  II CE$\nu$NS

The differential CE$\nu$NS cross section for a neutrino of flavor
$\mathcal{\ell}$ that scatters off a nucleus $\mathcal{N}$ with $Z$ protons
and $N$ neutrons is given by [1, 62, 63, 64, 65]

$\displaystyle\dfrac{d\sigma_{\nu_{\ell}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}}$
$\displaystyle=\dfrac{G_{\text{F}}^{2}M}{\pi}\left(1-\dfrac{MT_{\mathrm{nr}}}{2E_{\nu}^{2}}\right)\times$
$\displaystyle\left[g_{V}^{p}\left(\nu_{\ell}\right)ZF_{Z}\left(|\vec{q}|^{2}\right)+g_{V}^{n}NF_{N}\left(|\vec{q}|^{2}\right)\right]^{2},$
(1)

where $T_{\rm{nr}}$ is the nuclear recoil energy, $E_{\nu}$ is the neutrino
energy, and the term inside the square brackets is the so-called weak charge
of the nucleus $\mathcal{Q}_{W}$. In Eq. (1) $G_{\text{F}}$ is the Fermi
constant and $M$ the nuclear mass. For Ge ($Z$=32) we use $N$=(38, 40, 41, 42,
44) corresponding to the natural abundances of 0.2057 ($\mathrm{{}^{70}Ge}$),
0.2745 ($\mathrm{{}^{72}Ge}$), 0.0775 ($\mathrm{{}^{74}Ge}$), 0.3650
($\mathrm{{}^{74}Ge}$), 0.0773 ($\mathrm{{}^{76}Ge}$) [66].
$F_{Z}\left(|\vec{q}|^{2}\right)$ and $F_{N}\left(|\vec{q}|^{2}\right)$ are,
respectively, the proton and neutron form factors that parameterize the loss
of coherence for increasing values of the momentum transfer $|\vec{q}|$. In
the low-energy regime that characterizes reactor neutrino experiments, the
form factors are practically equal to one, so basically the result does not
depend on the parametrization used as well as on the values of the proton and
neutron distribution radii. The neutrino-nucleon couplings are computed using
the radiative corrections as described in Refs. [67, 64] that yield to
$g_{V}^{p}(\nu_{e})=0.0382$ and $g_{V}^{n}=-0.5117$.

For a CE$\nu$NS experiment located at a reactor site, only the $\bar{\nu}_{e}$
flavor contributes. The antineutrino spectrum $dN_{\bar{\nu}_{e}}/dE_{\nu}$ is
obtained by combining the expected spectra for $E_{\nu}>2$ MeV from Ref. [68,
69] with the low energy part determined in Ref. [70]. This parameterization is
referred to as HMVE, and it has been shown [39] that the usage of other
parameterizations [71, 72, 73] does not significantly affect the results of
the Dresden-II analysis.

When an isolated nuclear recoil occurs, the electron-equivalent (ee) nuclear
recoil energy observed in the detector, $E_{\rm{det}}$, is expressed as

$E_{\rm{det}}=f_{Q}(T_{\rm{nr}})T_{\rm{nr}},$ (2)

where $f_{Q}(T_{\rm{nr}})$ is the energy-dependent parameterization of the QF.
The theoretical CE$\nu$NS event rate is therefore

$\displaystyle\frac{dR}{dE_{\rm{det}}}=$ $\displaystyle
N_{T}(\textrm{Ge})\int_{E_{\nu}^{\rm{min}}}^{E_{\nu}^{\rm{max}}}dE_{\nu}\frac{dN_{\bar{\nu}_{e}}}{dE_{\nu}}\dfrac{d\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}}\times$
$\displaystyle\left(f_{Q}+T_{\rm{{nr}}}\frac{df_{Q}}{dT_{\rm{nr}}}\right)^{-1},$
(3)

where $N_{T}(\textrm{Ge})$ is the number of target atoms in the detector,
$E_{\nu}^{\text{min}}(T_{\text{nr}})\simeq\sqrt{MT_{\text{nr}}/2}$ and
$E_{\nu}^{\text{max}}\simeq 10$ MeV. The last term in Eq. (3) is
$dT_{\rm{nr}}/dE_{\rm{det}}$, which is needed to express the rate in terms of
the electron-equivalent nuclear recoil energy defined in Eq. (2).

##  III The Migdal effect

Despite the very intense $\bar{\nu}_{e}$ flux, the search for CE$\nu$NS from
reactor neutrinos is very challenging due to the tiny signal produced and the
still-under-debate behavior of the quenching factor at low energies. In this
scenario, it is crucial to characterize the impact of the Migdal effect to
prevent potential misinterpretations of experimental data. We first calculate
the Migdal rate using the formalism of Ibe et al. [49], which considers the
target as composed of isolated atoms. Moreover, this formalism relies on the
dipole approximation that allows one to write the Migdal transition matrix
element, $M_{fi}$, in the form

$\displaystyle M_{fi}$
$\displaystyle=\langle\psi_{f}|e^{-im_{e}\vec{v}\cdot\sum_{i=1}^{Z}\vec{r}_{i}}|\psi_{i}\rangle\simeq$
$\displaystyle\simeq-
im_{e}\vec{v}\cdot\langle\psi_{f}|\sum_{i=1}^{Z}\vec{r}_{i}|\psi_{i}\rangle$
$\displaystyle\equiv-im_{e}\vec{v}\cdot\vec{D}_{fi},$ (4)

where $\vec{r}_{i}$ is the position operator of the $Z$ electrons,
$\vec{D}_{fi}$ is the dipole matrix element, $m_{e}$ is the electron mass,
$\vec{v}$ is the nuclear recoil velocity, while $\psi_{f}$ and $\psi_{i}$ are
the wavefunctions of the final and the initial atomic states in the nucleus
rest frame. The final state wavefunctions are boosted to the rest frame of the
recoiling nucleus by a Galilean transformation and are computed using the
Dirac-Hartree-Fock method.

Under these assumptions, the differential cross section for the Migdal effect
can be written as

$\displaystyle\left(\dfrac{d\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}}\right)_{\rm{Migdal}}^{{\rm{Ibe}}\;et\;al.}=$
$\displaystyle\dfrac{G_{\text{F}}^{2}M}{\pi}\left(1-\dfrac{MT_{\mathrm{nr}}}{2E_{\nu}^{2}}\right)\mathcal{Q}_{W}^{2}\times$
$\displaystyle\left|Z_{\rm{ion}}(q_{e})\right|^{2},$ (5)

where $\left|Z_{\rm{ion}}(q_{e})\right|$ is the ionization rate of an
individual electron in the target with momentum $q_{e}$. It is defined as

$\left|Z_{\rm{ion}}(q_{e})\right|^{2}=\frac{1}{2\pi}\sum_{n,\ell}\int
dT_{e}\frac{d}{dT_{e}}p^{c}_{q_{e}}(n\ell\rightarrow T_{e}),$ (6)

where $p^{c}_{q_{e}}(n\ell\rightarrow T_{e})$ are the ionization probabilities
for an atomic electron with quantum numbers $n$ and $\ell$ that is ionized
with a final energy $T_{e}$. It should be noticed that very similar results
are expected if one relies on the probabilities calculated in Ref. [74] with
an independent approach. Indeed, the authors have shown a very good agreement
with the results obtained with the Ibe et al. formalism that is also used in
this work, demonstrating that the addition of semi-inclusive ionization
probabilities is not significant at low recoil energies for atomic germanium.

The double differential cross section for the $n\ell$ state contribution as a
function of both the electron and the nuclear recoil energy is hence

$\displaystyle\left(\dfrac{d^{2}\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}dT_{e}}\right)_{n\ell}^{{\rm{Ibe}}\;et\;al.}=$
$\displaystyle\dfrac{G_{\text{F}}^{2}M}{\pi}\left(1-\dfrac{MT_{\mathrm{nr}}}{2E_{\nu}^{2}}\right)\mathcal{Q}_{W}^{2}\times$
$\displaystyle\frac{1}{2\pi}\,\frac{d}{dT_{e}}p^{c}_{q_{e}}(n\ell\rightarrow
T_{e}).$ (7)

If the nuclear recoil is followed by a Migdal emission, the total energy
deposit of the event in the detector is

$E_{\rm{det}}=f_{Q}T_{\rm{nr}}+T_{e}+E_{n\ell},$ (8)

where the first term is the nuclear recoil energy deposit, while $T_{e}$ and
$E_{n\ell}$ account for the extra energy injected in the detector, $E_{n\ell}$
being the atomic de-excitation energy for Ge [49]. We evaluate the theoretical
event rate as a function of the detected energy, which is given by

$\displaystyle\left(\frac{dR}{dE_{\rm{det}}}\right)_{\textrm{Migdal}}^{{\rm{Ibe}}\;et\;al.}$
$\displaystyle=N_{T}(\textrm{Ge})\sum_{n,\ell}\int_{E_{\nu}^{\rm{min}}}^{E_{\nu}^{\rm{max}}}dE_{\nu}\frac{dN\bar{\nu}_{e}}{dE_{\nu}}\times$
$\displaystyle\int
dT_{e}\int_{T_{\textrm{nr}}^{\rm{min}}}^{T_{\rm{nr}}^{\rm{max}}}dT_{\rm{nr}}\left(\dfrac{d^{2}\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}dT_{e}}\right)_{n\ell}^{{\rm{Ibe\;}}et\;al.}\times$
$\displaystyle\delta(E_{\rm{det}}-f_{Q}T_{\rm{nr}}-T_{e}-E_{n\ell}),$ (9)

where we have imposed energy conservation using the Dirac $\delta$ and
$T_{\rm{nr}}$ is now constrained within the values ${T_{\rm{nr}}^{\rm{min}}}$
and ${T_{\rm{nr}}^{\rm{max}}}$ given by [57]

$\frac{\left(T_{e}+E_{n\ell}\right)^{2}}{2M}\leq
T_{\rm{nr}}\leq\frac{\left(2E_{\nu}-(T_{e}+E_{n\ell})\right)^{2}}{2(M+2E_{\nu})}.$
(10)

The rate in Eq. (9) represents the Migdal contribution summed over all the
possible $n\ell$ atomic states. The total predicted event rate is thus given
by the sum of Eq. (3) and Eq. (9).

##  IV Migdal Photo Absorption

The formalism described so far to compute the dipole matrix element for the
Migdal rate relies on the assumption that the target atom is isolated. While
this assumption is acceptable for noble elements, as for argon or xenon
detectors [58], it is expected to be less valid in semiconductors, where
solid-state effects should be considered. However, developing a first-
principle theory that goes beyond the isolated atom approximation is
challenging because of the many-body effects that need to be taken into
account. Remarkably, the formalism developed in Ref. [52] relates the
photoabsorption cross section $\sigma_{\gamma}$ to the dipole matrix element,
necessary to compute the Migdal ionization rate, without requiring any many-
body calculation. This scheme will be referred to as Migdal photoabsorption
approximation (MPA). One of the major advantages of MPA is that the
photoabsorption cross section is experimentally known, such that the Migdal
rate suffers from very small uncertainties [52], well below the precision
required in this work. MPA has been so far adopted in the context of dark
matter searches, where the power of the formalism has been proved by comparing
it to other computations for silicon and xenon [49, 50]. However, MPA has
never been exploited in the context of neutrino scattering. Here, for the
first time, we use it in this context and we compare its predictions with the
formalism of Ibe et al. for germanium detectors. Explicitly, we derive the
Migdal contribution to the CE$\nu$NS cross section under MPA as (see Appendix
A for further information)

$\displaystyle\left(\dfrac{d^{2}\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}dE_{r}}\right)_{\rm{Migdal}}^{\rm{MPA}}=$
$\displaystyle\dfrac{G_{\text{F}}^{2}M}{\pi}\left(1-\dfrac{MT_{\mathrm{nr}}}{2E_{\nu}^{2}}\right)\mathcal{Q}_{W}^{2}\times$
$\displaystyle\frac{1}{2\pi^{2}\alpha_{\rm{EM}}}\frac{m^{2}_{e}}{M}\frac{T_{\rm{nr}}}{E_{r}}\sigma_{\gamma}^{\rm{Ge}}(E_{r}),$
(11)

where $\alpha_{\rm{EM}}$ is the fine structure constant and $E_{r}$ is the
energy deposit due to atomic excitation or ionization such that
$E_{\rm{det}}=f_{Q}T_{\rm{nr}}+E_{r}$. The photoabsorption cross section
$\sigma_{\gamma}^{\rm{Ge}}(E_{r})$ for Ge has been taken from Refs. [75, 76]
for $E_{r}\geq 10\;\textrm{eV}_{\textrm{ee}}$. The theoretical event rate as a
function of $E_{\rm{det}}$ is obtained by integrating over all possible
nuclear recoil and neutrino energies and imposing energy conservation.

It should be noticed that depending on the crystal scale that one is able to
probe, other effects that account for the response of multiple atoms at once
should be considered, as they have been proved to highly enhance the Migdal
rate [77, 78, 79, 80]. However, although the range of $T_{\rm{nr}}$ in Eq.
(10) includes very small values, the main contribution in current Ge CE$\nu$NS
reactor experiments comes from $T_{\rm{nr}}$ of the order of 1 $\mathrm{keV}$.
Thus, the momentum transfer is $|\vec{q}|\simeq\sqrt{2MT_{\rm{nr}}}\sim
10\;\rm{MeV}$ with a corresponding de Broglie wavelength of about 20 fm. The
latter is much smaller than the scale of the interparticle spacing in the
crystal so, in this work, we can safely neglect these effects. We point here,
that at lower ionization energies, multiple atom effects must be better taken
into consideration. Moreover, it has been suggested that in germanium a larger
amount of secondary nuclear recoils may be produced following a low-energy
primary one [43]. This, in addition to possible more complex crystal response
models [81, 82, 83], can significantly affect the microscopic description of
the Migdal effect. However, the inclusion of these possible effects, which
still need to be investigated both theoretically and experimentally, is beyond
the scope of this work.

##  V Results and Discussion

In order to compare the impact of the Migdal contribution to the standard
CE$\nu$NS rate, we consider a 1 kg Ge detector located $10\,\rm{m}$ away from
a $3\,\rm{GW}_{\rm{th}}$ reactor power plant whose $\bar{\nu}_{e}$ spectrum is
given by the HMVE parameterization. This specific configuration resembles that
of current CE$\nu$NS reactor experiments, like Dresden-II, CONUS and $\nu$GEN.

Figure 1: CE$\nu$NS theoretical differential rate for a Ge detector located
$10\,\rm{m}$ away from a $3\,\rm{GW}_{\rm{th}}$ reactor using the Lindhard QF
(black solid line). We also show the Migdal rate obtained under the Ibe et al.
formalism (solid blue line), highlighting also the contributions of the
different $n=2,3,4$ atomic shells (dashed curves), and the rate obtained with
the MPA formalism (solid red line).

Under these assumptions, in Fig. 1 we show the theoretical CE$\nu$NS rate
obtained with Eq. (3) as a function of the detected energy, adopting for the
QF the standard Lindhard parameterization. We also show the Migdal rate
determined with the Ibe et al. formalism as in Eq. (9), isolating the
contributions from the different $n$ shells, obtained by summing over all the
contributions from different orbital angular momenta in the initial state
$\ell$. We compare this result with the Migdal prediction under the MPA
scheme, finding intriguingly that the two formalisms give practically
identical results in the energy range considered. Moreover, in both cases, the
Migdal contribution is completely subdominant with respect to the CE$\nu$NS
one for energies below $\sim 0.6\;\textrm{keV}_{\rm{ee}}$, with the most
significant contribution given by the $n=2,\,3$ shells. Above $\sim
0.6\;\textrm{keV}_{\rm{ee}}$ it starts to dominate222 Note that a similar
trend is also found in Ref. [58], where a comparison between the CE$\nu$NS
rate and the Migdal contribution using the formalism of Ibe et al. has been
evaluated for xenon and argon detectors in a reactor site. We have used their
publicly available code [58] to verify our results., and it could provide the
possibility to observe CE$\nu$NS above this threshold, even if being so small
it would require extremely low levels of background.

Figure 2: Expected number of CE$\nu$NS events in the Dresden-II detector
obtained using different quenching factors, i.e. Fef (blue line), YBe (cyan
line) and Lindhard (purple line). The Migdal contribution corresponds to the
red curve, while the neutrino electron scattering ($\nu$ES) contribution is
given by the green line. In the top panel, we compare these curves with the
Dresden-II reactor ON data, while in the bottom panel, we show the Dresden-II
data residuals after background subtraction. The inset shows a comparison of
Dresden-II reactor ON (Rx-ON) [14], $\nu$GEN [12] and CONUS (C1 Run-1) [9]
data, all rescaled to the same units.

We now focus on the Dresden-II science case. Here, as already stated in the
introduction, the Migdal effect has been invoked as a possible explanation of
the enhancement measured in the Fef and YBe quenching factors at low energies
that in turn enabled the observation of CE$\nu$NS in the Dresden-II data. In
the top panel of Fig. 2, we show the Dresden-II reactor ON (Rx-ON) data along
with the standard CE$\nu$NS predictions obtained with three different QFs,
namely Lindhard, Fef and YBe. To derive these spectra we used all the
experimental information on the Dresden-II detector, including energy-smearing
effects, following Refs. [14, 40, 39]. In the bottom panel of Fig. 2 we show
the same spectra for the three QFs but we compare them to the Dresden-II data
residuals after background subtraction [14]. It is evident that only the Fef
and marginally the YBe QFs fit the excess and lead to a statistically
significant CE$\nu$NS observation for $E_{\rm{det}}\lesssim
0.3\,\mathrm{keV}_{\mathrm{ee}}$. In the same figure, we also show the Migdal
contribution using the MPA formalism. It is clear that adding the Migdal
contribution to the CE$\nu$NS Lindhard prediction is not sufficient to explain
the CE$\nu$NS Fef or YBe predictions, given that the former is completely
negligible with respect to the CE$\nu$NS signal. Moreover, we find that the
neutrino-electron scattering ($\nu$ES) [39, 84] overcomes the Migdal rate for
$E_{\rm{det}}\gtrsim 0.7\;\rm{keV}_{\rm{ee}}$, as shown in the top panel of
Fig. 2. Overall, both Migdal and $\nu$ES rates are so small that with the
current experimental precision can be overlooked in SM CE$\nu$NS searches.
However, in some scenarios of physics beyond the SM, like non-standard
properties of neutrinos, their contribution could be significantly enhanced
and thus they need to be taken into account [40, 39] to derive meaningful
limits.

In the inset of Fig. 2, for comparison purposes, we show a review of the
existing data from germanium detectors searching for CE$\nu$NS at a reactor
site, i.e. Dresden-II [14], $\nu$GEN [12] and CONUS (C1 Run-1) [9].
Interestingly, despite the fact that CONUS and $\nu$Gen have reached a much
lower background level compared to Dresden-II, they have not detected
CE$\nu$NS yet. Having a conservative threshold of $\sim
0.3\;\rm{keV}_{\rm{ee}}$, future data will be needed to confirm the Dresden-II
observation of CE$\nu$NS, which is localized below this threshold.
Nevertheless, despite the low background level reached, the Migdal
contribution is so small that it could be safely neglected also in experiments
like $\nu$GEN and CONUS, which show a good agreement with the expected
background. Similar conclusions are expected also for silicon detectors like
CONNIE [10, 15] that operate in a similar energy range.

To conclude, we have shown that the Migdal contribution is orders-of-magnitude
subdominant in the region of interest for reactor CE$\nu$NS searches with
germanium detectors, independently of the formalism used to model the Migdal
effect. Thus, the enhancement of the quenching factor at low energies found in
Ref. [14] that enabled the observation of CE$\nu$NS at the Dresden-II site
requires a different explanation than the standard Migdal effect.

###### Acknowledgements.

The authors are thankful to Simon Knapen for a fruitful discussion on the
topic. The work of C. Giunti is supported by the PRIN 2022 research grant
“Addressing systematic uncertainties in searches for dark matter”, Number
2022F2843L, funded by MIUR.

##  Appendix A Derivation of the Migdal cross section with the MPA formalism

While the standard Migdal formalism derived by Ibe et al. has been already
discussed in previous studies in the context of CE$\nu$NS interactions (see
e.g. Ref. [49] and Refs. [57, 58] with their publicly-available codes), little
literature is present for the so-called Migdal photoabsorption approximation
(MPA) as this has never been used in this context. In this Appendix, we expand
the discussion on the MPA formalism providing the reader with more details. In
this scenario, the Migdal transition matrix element under the dipole
approximation, $|M_{fi}|^{2}$, is related to the ionization probability due to
the Migdal effect. It is convenient to define the double differential cross
section for an electron neutrino that scatters off a nucleus $\mathcal{N}$ as

$\displaystyle\left(\dfrac{d^{2}\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}dE_{r}}\right)_{\rm{Migdal}}=$
$\displaystyle\left(\dfrac{d\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dT_{\mathrm{nr}}}\right)|M_{fi}|^{2},$
(12)

where $T_{\rm{nr}}$ is the nuclear recoil energy, and the matrix element is a
function of the energy deposit $E_{r}$. The square of the transition matrix
element in Eq. (12) can be recast to

$\displaystyle|M_{fi}|^{2}=\left(\frac{m_{e}}{M}\right)^{2}2MT_{\rm{nr}}\overline{D^{2}_{fi}},$
(13)

where $M$ is the nuclear mass, $m_{e}$ is the electron mass and we have used
the fact that $M|\vec{v}|=\sqrt{2MT_{\rm{nr}}}$ is the momentum transfer of
the recoiling nucleus, with $\vec{v}$ the nuclear recoil velocity. Here,
$\overline{D^{2}_{fi}}$ is the average of the squared dipole matrix element.
Under the MPA scheme, one has to relate $\overline{D^{2}_{fi}}$ with the
photoabsorption cross section $\sigma_{\gamma}$ through [52]

$\displaystyle\overline{D^{2}_{fi}}=\frac{\sigma_{\gamma}(E_{r})}{4\pi^{2}\alpha_{\rm{EM}}E_{r}},$
(14)

where $\alpha_{\rm{EM}}$ is the fine structure constant.

Figure 3: Comparison between the Migdal cross section (under the MPA scheme)
and the CE$\nu$NS one as a function of the observed energy $E_{\rm{det}}$,
assuming the standard Lindhard quenching factor and a neutrino energy of
$E_{\nu}=6$ MeV.

In Fig. 3 we show the comparison between the differential Migdal cross section
and the CE$\nu$NS one as a function of the observed energy $E_{\rm{det}}$, for
a neutrino energy of $E_{\nu}=6\;\rm{MeV}$ and assuming the standard Lindhard
quenching factor. By substituting Eq. (14) into Eq. (13) and using Eq. (12)
one obtains Eq. (11).

In Fig. 4 we show the germanium photoabsorption cross section extracted from
Ref. [75] as used in this work, among with the differential Migdal cross
section as a function of $E_{\rm{det}}$ for $E_{\nu}=3\;$MeV.

Figure 4: Photoabsorption cross section for Germanium taken from Ref. [75]
(purple) and differential Migdal cross section for $E_{\nu}=3\;$MeV. The dots
represent the photoabsorption cross section for each x-ray emission as listed
in Ref. [75].

The two cross sections are reported at different energy scales indicated in
the corresponding $y$ axes.

The differential Migdal cross section as a function of $E_{\rm{det}}$ is
obtained by integrating over all possible nuclear recoils and is given by

$\displaystyle\left(\dfrac{d\sigma_{\bar{\nu}_{e}\text{-}\mathcal{N}}}{dE_{\rm{det}}}\right)_{\rm{Migdal}}^{\rm{MPA}}=\dfrac{G_{\text{F}}^{2}M}{\pi}\int_{T_{\textrm{nr}}^{\rm{min}}}^{T_{\rm{nr}}^{\rm{max}}}dT_{\rm{nr}}\left(1-\dfrac{MT_{\mathrm{nr}}}{2E_{\nu}^{2}}\right)\mathcal{Q}_{W}^{2}\times$
$\displaystyle\frac{1}{2\pi^{2}\alpha_{\rm{EM}}}\frac{m^{2}_{e}}{M}\frac{T_{\rm{nr}}}{E_{\rm{det}}-f_{Q}T_{\rm{nr}}}\sigma_{\gamma}^{\rm{Ge}}(E_{\rm{det}}-f_{Q}T_{\rm{nr}}),$
(15)

where $G_{\text{F}}$ is the Fermi constant, $\mathcal{Q}_{W}$ is the weak
charge of the nucleus and $f_{Q}$ is the energy-dependent parameterization of
the quenching factor. In Fig. 4, it is possible to see that the features that
characterize the differential Migdal cross section within the MPA scheme are
tightly connected with the photoabsorption peaks.

## References

  * Freedman [1974] D. Z. Freedman, Coherent effects of a weak neutral current, Phys. Rev. D 9, 1389 (1974).
  * Akimov  _et al._ [2017] D. Akimov _et al._ (COHERENT), Observation of Coherent Elastic Neutrino-Nucleus Scattering, Science 357, 1123 (2017), arXiv:1708.01294 [nucl-ex] .
  * Akimov  _et al._ [2021] D. Akimov _et al._ (COHERENT), First Detection of Coherent Elastic Neutrino-Nucleus Scattering on Argon, Phys.Rev.Lett. 126, 012002 (2021), arXiv:2003.10630 [nucl-ex] .
  * Akimov  _et al._ [2022a] D. Akimov _et al._ (COHERENT), Measurement of the Coherent Elastic Neutrino-Nucleus Scattering Cross Section on CsI by COHERENT, Phys.Rev.Lett. 129, 081801 (2022a), arXiv:2110.07730 [hep-ex] .
  * Angloher  _et al._ [2019] G. Angloher _et al._ (NUCLEUS), Exploring $\hbox{CE}\nu\hbox{NS}$ with NUCLEUS at the Chooz nuclear power plant, Eur. Phys. J. C 79, 1018 (2019), arXiv:1905.10258 [physics.ins-det] .
  * Strauss  _et al._ [2017] R. Strauss _et al._ , The $\nu$-cleus experiment: A gram-scale fiducial-volume cryogenic detector for the first detection of coherent neutrino-nucleus scattering, Eur. Phys. J. C 77, 506 (2017), arXiv:1704.04320 [physics.ins-det] .
  * Akimov  _et al._ [2022b] D. Akimov _et al._ , The COHERENT Experimental Program, in  _Snowmass 2021_ (2022) arXiv:2204.04575 [hep-ex] .
  * Baxter  _et al._ [2020] D. Baxter _et al._ , Coherent Elastic Neutrino-Nucleus Scattering at the European Spallation Source, JHEP 02, 123, arXiv:1911.00762 [physics.ins-det] .
  * Bonet  _et al._ [2021a] H. Bonet _et al._ (CONUS), Constraints on elastic neutrino nucleus scattering in the fully coherent regime from the CONUS experiment, Phys. Rev. Lett. 126, 041804 (2021a), arXiv:2011.00210 [hep-ex] .
  * Aguilar-Arevalo  _et al._ [2016] A. Aguilar-Arevalo _et al._ (CONNIE), Results of the Engineering Run of the Coherent Neutrino Nucleus Interaction Experiment (CONNIE), JINST 11 (07), P07024, arXiv:1604.01343 [physics.ins-det] .
  * Billard  _et al._ [2017] J. Billard _et al._ , Coherent Neutrino Scattering with Low Temperature Bolometers at Chooz Reactor Complex, J. Phys. G 44, 105101 (2017), arXiv:1612.09035 [physics.ins-det] .
  * Alekseev  _et al._ [2022] I. Alekseev _et al._ ($\nu$GeN), First results of the $\nu$GeN experiment on coherent elastic neutrino-nucleus scattering, Phys. Rev. D 106, L051101 (2022), arXiv:2205.04305 [nucl-ex] .
  * Colaresi  _et al._ [2021] J. Colaresi, J. I. Collar, T. W. Hossbach, A. R. L. Kavner, C. M. Lewis, A. E. Robinson, and K. M. Yocum, First results from a search for coherent elastic neutrino-nucleus scattering at a reactor site, Phys. Rev. D 104, 072003 (2021), arXiv:2108.02880 [hep-ex] .
  * Colaresi  _et al._ [2022] J. Colaresi, J. I. Collar, T. W. Hossbach, C. M. Lewis, and K. M. Yocum, Measurement of Coherent Elastic Neutrino-Nucleus Scattering from Reactor Antineutrinos, Phys. Rev. Lett. 129, 211802 (2022), arXiv:2202.09672 [hep-ex] .
  * Aguilar-Arevalo  _et al._ [2022] A. Aguilar-Arevalo _et al._ (CONNIE), Search for coherent elastic neutrino-nucleus scattering at a nuclear reactor with CONNIE 2019 data, JHEP 05, 017, arXiv:2110.13033 [hep-ex] .
  * Agnolet  _et al._ [2017] G. Agnolet _et al._ (MINER), Background Studies for the MINER Coherent Neutrino Scattering Reactor Experiment, Nucl. Instrum. Meth. A 853, 53 (2017), arXiv:1609.02066 [physics.ins-det] .
  * Su  _et al._ [2023] C. Su, Q. Liu, and T. Liang (CE$\nu$NS@CSNS), CE$\nu$NS Experiment Proposal at CSNS (2023), arXiv:2303.13423 [physics.ins-det] .
  * Colas  _et al._ [2022] J. Colas, J. Billard, S. Ferriol, J. Gascon, and T. Salagnac (RICOCHET), Development of data processing and analysis pipeline for the RICOCHET experiment, J.Low Temp.Phys. (2022), arXiv:2111.12856 [physics.ins-det] .
  * Akimov  _et al._ [2022c] D. Y. Akimov _et al._ , The RED-100 experiment, JINST 17 (11), T11011, arXiv:2209.15516 [physics.ins-det] .
  * Wong [2015] H. T.-K. Wong, Taiwan EXperiment On NeutrinO — History and Prospects, The Universe 3, 22 (2015), arXiv:1608.00306 [hep-ex] .
  * Cadeddu  _et al._ [2021a] M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. Li, E. Picciau, and Y. Zhang, Constraints on light vector mediators through coherent elastic neutrino nucleus scattering data from COHERENT, JHEP 2101, 116, arXiv:2008.05022 [hep-ph] .
  * Cadeddu and Dordei [2019] M. Cadeddu and F. Dordei, Reinterpreting the weak mixing angle from atomic parity violation in view of the Cs neutron rms radius measurement from COHERENT, Phys.Rev. D99, 033010 (2019), arXiv:1808.10202 [hep-ph] .
  * De Romeri  _et al._ [2023] V. De Romeri, O. G. Miranda, D. K. Papoulias, G. Sanchez Garcia, M. Tortola, and J. W. F. Valle, Physics implications of a combined analysis of COHERENT CsI and LAr data, JHEP 04, 035, arXiv:2211.11905 [hep-ph] .
  * Giunti [2020] C. Giunti, General COHERENT Constraints on Neutrino Non-Standard Interactions, Phys.Rev. D101, 035039 (2020), arXiv:1909.00466 [hep-ph] .
  * Lindner  _et al._ [2017] M. Lindner, W. Rodejohann, and X.-J. Xu, Coherent Neutrino-Nucleus Scattering and new Neutrino Interactions, JHEP 1703, 097, arXiv:1612.04150 [hep-ph] .
  * Billard  _et al._ [2018] J. Billard, J. Johnston, and B. J. Kavanagh, Prospects for exploring New Physics in Coherent Elastic Neutrino-Nucleus Scattering, JCAP 1811, 016, arXiv:1805.01798 [hep-ph] .
  * Cadeddu  _et al._ [2018a] M. Cadeddu, C. Giunti, Y. F. Li, and Y. Y. Zhang, Average CsI neutron density distribution from COHERENT data, Phys.Rev.Lett. 120, 072501 (2018a), arXiv:1710.02730 [hep-ph] .
  * Cadeddu  _et al._ [2021b] M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. Li, E. Picciau, C. Ternes, and Y. Zhang, New insights into nuclear physics and weak mixing angle using electroweak probes, Phys.Rev.C 104, 065502 (2021b), arXiv:2102.06153 [hep-ph] .
  * Cadeddu  _et al._ [2018b] M. Cadeddu, C. Giunti, K. Kouzakov, Y. F. Li, A. Studenikin, and Y. Y. Zhang, Neutrino Charge Radii from COHERENT Elastic Neutrino-Nucleus Scattering, Phys.Rev. D98, 113010 (2018b), arXiv:1810.05606 [hep-ph] .
  * Cadeddu  _et al._ [2020a] M. Cadeddu, F. Dordei, C. Giunti, Y. Li, E. Picciau, and Y. Zhang, Physics results from the first COHERENT observation of CE$\nu$NS in argon and their combination with cesium-iodide data, Phys.Rev. D102, 015030 (2020a), arXiv:2005.01645 [hep-ph] .
  * Aristizabal Sierra  _et al._ [2018] D. Aristizabal Sierra, V. De Romeri, and N. Rojas, COHERENT analysis of neutrino generalized interactions, Phys.Rev. D98, 075018 (2018), arXiv:1806.07424 [hep-ph] .
  * Cadeddu  _et al._ [2020b] M. Cadeddu, F. Dordei, C. Giunti, Y. Li, and Y. Zhang, Neutrino, Electroweak and Nuclear Physics from COHERENT Elastic Neutrino-Nucleus Scattering with Refined Quenching Factor, Phys.Rev. D101, 033004 (2020b), arXiv:1908.06045 [hep-ph] .
  * Corona  _et al._ [2022] M. A. Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. Li, E. Picciau, C. A. Ternes, and Y. Zhang, Probing light mediators and $(g-2)_{\mu}$ through detection of coherent elastic neutrino nucleus scattering at COHERENT, JHEP 05, 109, arXiv:2202.11002 [hep-ph] .
  * Denton  _et al._ [2018] P. B. Denton, Y. Farzan, and I. M. Shoemaker, Testing large non-standard neutrino interactions with arbitrary mediator mass after COHERENT data, JHEP 07, 037, arXiv:1804.03660 [hep-ph] .
  * Denton and Gehrlein [2021] P. B. Denton and J. Gehrlein, A Statistical Analysis of the COHERENT Data and Applications to New Physics, JHEP 04, 266, arXiv:2008.06062 [hep-ph] .
  * Co’  _et al._ [2020] G. Co’, M. Anguiano, and A. M. Lallena, Nuclear structure uncertainties in coherent elastic neutrino-nucleus scattering, JCAP 04, 044, arXiv:2001.04684 [nucl-th] .
  * Bonet  _et al._ [2021b] H. Bonet _et al._ , Large-size sub-keV sensitive germanium detectors for the CONUS experiment, Eur. Phys. J. C 81, 267 (2021b), arXiv:2010.11241 [physics.ins-det] .
  * Augier  _et al._ [2023] C. Augier _et al._ (Ricochet), Fast neutron background characterization of the future Ricochet experiment at the ILL research nuclear reactor, Eur. Phys. J. C 83, 20 (2023), arXiv:2208.01760 [astro-ph.IM] .
  * Atzori Corona  _et al._ [2022] M. Atzori Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, C. A. Ternes, and Y. Y. Zhang, Impact of the Dresden-II and COHERENT neutrino scattering data on neutrino electromagnetic properties and electroweak physics, JHEP 09, 164, arXiv:2205.09484 [hep-ph] .
  * Coloma  _et al._ [2022] P. Coloma, I. Esteban, M. C. Gonzalez-Garcia, L. Larizgoitia, F. Monrabal, and S. Palomares-Ruiz, Bounds on new physics with data of the Dresden-II reactor experiment and COHERENT, JHEP 05, 037, arXiv:2202.10829 [hep-ph] .
  * Majumdar  _et al._ [2022] A. Majumdar, D. K. Papoulias, R. Srivastava, and J. W. F. Valle, Physics implications of recent Dresden-II reactor data, Phys. Rev. D 106, 093010 (2022), arXiv:2208.13262 [hep-ph] .
  * Denton and Gehrlein [2022] P. B. Denton and J. Gehrlein, New constraints on the dark side of non-standard interactions from reactor neutrino scattering data, Phys. Rev. D 106, 015022 (2022), arXiv:2204.09060 [hep-ph] .
  * Collar  _et al._ [2021] J. I. Collar, A. R. L. Kavner, and C. M. Lewis, Germanium response to sub-keV nuclear recoils: a multipronged experimental characterization, Phys. Rev. D 103, 122003 (2021), arXiv:2102.10089 [nucl-ex] .
  * Lindhard  _et al._ [1963] J. Lindhard, V. Nielsen, M. Scharff, and P. V. Thomsen, Integral equations governing radiation effects. (notes on atomic collisions, iii), Kgl. Danske Videnskab., Selskab. Mat. Fys. Medd. 33, 10,  (1963).
  * Albakry  _et al._ [2022] M. F. Albakry _et al._ (SuperCDMS), Ionization yield measurement in a germanium CDMSlite detector using photo-neutron sources, Phys. Rev. D 105, 122002 (2022), arXiv:2202.07043 [physics.ins-det] .
  * Bonhomme  _et al._ [2022] A. Bonhomme _et al._ , Direct measurement of the ionization quenching factor of nuclear recoils in germanium in the keV energy range, Eur. Phys. J. C 82, 815 (2022), arXiv:2202.03754 [physics.ins-det] .
  * Xu  _et al._ [2023] J. Xu _et al._ , Search for the Migdal effect in liquid xenon with keV-level nuclear recoils (2023), arXiv:2307.12952 [hep-ex] .
  * Migdal [1941] A. Migdal, Ionization of atoms accompanying $\alpha$\- and $\beta$-decay, J. Phys, USSR 4, 449 (1941).
  * Ibe  _et al._ [2018] M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, Migdal Effect in Dark Matter Direct Detection Experiments, JHEP 03, 194, arXiv:1707.07258 [hep-ph] .
  * Essig  _et al._ [2020] R. Essig, J. Pradler, M. Sholapurkar, and T.-T. Yu, Relation between the Migdal Effect and Dark Matter-Electron Scattering in Isolated Atoms and Semiconductors, Phys. Rev. Lett. 124, 021801 (2020), arXiv:1908.10881 [hep-ph] .
  * Grilli di Cortona  _et al._ [2020] G. Grilli di Cortona, A. Messina, and S. Piacentini, Migdal effect and photon Bremsstrahlung: improving the sensitivity to light dark matter of liquid argon experiments, JHEP 11, 034, arXiv:2006.02453 [hep-ph] .
  * Liu  _et al._ [2020] C. P. Liu, C.-P. Wu, H.-C. Chi, and J.-W. Chen, Model-independent determination of the Migdal effect via photoabsorption, Phys. Rev. D 102, 121303 (2020), arXiv:2007.10965 [hep-ph] .
  * Agnes  _et al._ [2023] P. Agnes _et al._ (DarkSide), Search for Dark-Matter–Nucleon Interactions via Migdal Effect with DarkSide-50, Phys. Rev. Lett. 130, 101001 (2023), arXiv:2207.11967 [hep-ex] .
  * Aprile  _et al._ [2019] E. Aprile _et al._ (XENON), Search for Light Dark Matter Interactions Enhanced by the Migdal Effect or Bremsstrahlung in XENON1T, Phys. Rev. Lett. 123, 241803 (2019), arXiv:1907.12771 [hep-ex] .
  * Albakry  _et al._ [2023] M. F. Albakry _et al._ (SuperCDMS), Search for low-mass dark matter via bremsstrahlung radiation and the Migdal effect in SuperCDMS, Phys. Rev. D 107, 112013 (2023), arXiv:2302.09115 [hep-ex] .
  * Bell  _et al._ [2021] N. F. Bell, J. B. Dent, B. Dutta, S. Ghosh, J. Kumar, and J. L. Newstead, Low-mass inelastic dark matter direct detection via the Migdal effect, Phys. Rev. D 104, 076013 (2021), arXiv:2103.05890 [hep-ph] .
  * Bell  _et al._ [2020] N. F. Bell, J. B. Dent, J. L. Newstead, S. Sabharwal, and T. J. Weiler, Migdal effect and photon bremsstrahlung in effective field theories of dark matter direct detection and coherent elastic neutrino-nucleus scattering, Phys. Rev. D 101, 015012 (2020), arXiv:1905.00046 [hep-ph] .
  * Bell  _et al._ [2022] N. F. Bell, J. B. Dent, R. F. Lang, J. L. Newstead, and A. C. Ritter, Observing the Migdal effect from nuclear recoils of neutral particles with liquid xenon and argon detectors, Phys. Rev. D 105, 096015 (2022), arXiv:2112.08514 [hep-ph] .
  * Liao  _et al._ [2021] J. Liao, H. Liu, and D. Marfatia, Coherent neutrino scattering and the Migdal effect on the quenching factor, Phys. Rev. D 104, 015005 (2021), arXiv:2104.01811 [hep-ph] .
  * Liao  _et al._ [2022] J. Liao, H. Liu, and D. Marfatia, Implications of the first evidence for coherent elastic scattering of reactor neutrinos, Phys. Rev. D 106, L031702 (2022), arXiv:2202.10622 [hep-ph] .
  * Sorensen [2015] P. Sorensen, Atomic limits in the search for galactic dark matter, Phys. Rev. D 91, 083509 (2015), arXiv:1412.3028 [astro-ph.IM] .
  * Drukier and Stodolsky [1984] A. Drukier and L. Stodolsky, Principles and applications of a neutral-current detector for neutrino physics and astronomy, Phys. Rev. D 30, 2295 (1984).
  * Barranco  _et al._ [2005] J. Barranco, O. G. Miranda, and T. I. Rashba, Probing new physics with coherent neutrino scattering off nuclei, JHEP 12, 021, arXiv:hep-ph/0508299 .
  * Atzori Corona  _et al._ [2023a] M. Atzori Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, and G. Masia, Nuclear neutron radius and weak mixing angle measurements from latest COHERENT CsI and atomic parity violation Cs data, Eur. Phys. J. C 83, 683 (2023a), arXiv:2303.09360 [nucl-ex] .
  * Cadeddu  _et al._ [2023] M. Cadeddu, F. Dordei, and C. Giunti, A view of Coherent Elastic Neutrino-Nucleus Scattering (2023), arXiv:2307.08842 [hep-ph] .
  * Berglund and Wieser [2011] M. Berglund and M. E. Wieser, Isotopic compositions of the elements 2009 (iupac technical report), Pure and Applied Chemistry 83, 397 (2011).
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2022, 083C01 (2022).
  * Huber [2011] P. Huber, On the determination of anti-neutrino spectra from nuclear reactors, Phys. Rev. C 84, 024617 (2011), [Erratum: Phys.Rev.C 85, 029901 (2012)], arXiv:1106.0687 [hep-ph] .
  * Mueller  _et al._ [2011] T. A. Mueller _et al._ , Improved Predictions of Reactor Antineutrino Spectra, Phys. Rev. C 83, 054615 (2011), arXiv:1101.2663 [hep-ex] .
  * Vogel and Engel [1989] P. Vogel and J. Engel, Neutrino electromagnetic form factors, Phys. Rev. D 39, 3378 (1989).
  * Estienne  _et al._ [2019] M. Estienne _et al._ , Updated Summation Model: An Improved Agreement with the Daya Bay Antineutrino Fluxes, Phys. Rev. Lett. 123, 022502 (2019), arXiv:1904.09358 [nucl-ex] .
  * Kopeikin  _et al._ [2000] V. I. Kopeikin, L. A. Mikaelyan, and V. V. Sinev, Search for the neutrino magnetic moment in the nonequilibrium reactor anti-neutrino energy spectrum, Phys. Atom. Nucl. 63, 1012 (2000), arXiv:hep-ph/9904384 .
  * Kopeikin [2012] V. I. Kopeikin, Flux and spectrum of reactor antineutrinos, Phys. Atom. Nucl. 75, 143 (2012).
  * Cox  _et al._ [2023] P. Cox, M. J. Dolan, C. McCabe, and H. M. Quiney, Precise predictions and new insights for atomic ionization from the Migdal effect, Phys. Rev. D 107, 035032 (2023), arXiv:2208.12222 [hep-ph] .
  * Henke  _et al._ [1993] B. Henke, E. Gullikson, and J. Davis, X-ray interactions: Photoabsorption, scattering, transmission, and reflection at e = 50-30,000 ev, z = 1-92, Atomic Data and Nuclear Data Tables 54, 181 (1993).
  * Chen  _et al._ [2015] J.-W. Chen, H.-C. Chi, K.-N. Huang, H.-B. Li, C. P. Liu, L. Singh, H. T. Wong, C.-L. Wu, and C.-P. Wu, Constraining neutrino electromagnetic properties by germanium detectors, Phys. Rev. D 91, 013005 (2015), arXiv:1411.0574 [hep-ph] .
  * Knapen  _et al._ [2021] S. Knapen, J. Kozaczuk, and T. Lin, Migdal Effect in Semiconductors, Phys. Rev. Lett. 127, 081805 (2021), arXiv:2011.09496 [hep-ph] .
  * Knapen  _et al._ [2022] S. Knapen, J. Kozaczuk, and T. Lin, python package for dark matter scattering in dielectric targets, Phys. Rev. D 105, 015014 (2022), arXiv:2104.12786 [hep-ph] .
  * Adams  _et al._ [2023] D. Adams, D. Baxter, H. Day, R. Essig, and Y. Kahn, Measuring the Migdal effect in semiconductors for dark matter detection, Phys. Rev. D 107, L041303 (2023), arXiv:2210.04917 [hep-ph] .
  * Liang  _et al._ [2022] Z.-L. Liang, C. Mo, F. Zheng, and P. Zhang, Phonon-mediated Migdal effect in semiconductor detectors, Phys. Rev. D 106, 043004 (2022), [Erratum: Phys.Rev.D 106, 109901 (2022)], arXiv:2205.03395 [hep-ph] .
  * Kahn  _et al._ [2021] Y. Kahn, G. Krnjaic, and B. Mandava, Dark Matter Detection with Bound Nuclear Targets: The Poisson Phonon Tail, Phys. Rev. Lett. 127, 081804 (2021), arXiv:2011.09477 [hep-ph] .
  * Kurinsky  _et al._ [2020] N. Kurinsky, D. Baxter, Y. Kahn, and G. Krnjaic, Dark matter interpretation of excesses in multiple direct detection experiments, Phys. Rev. D 102, 015017 (2020), arXiv:2002.06937 [hep-ph] .
  * Trickle  _et al._ [2020] T. Trickle, Z. Zhang, K. M. Zurek, K. Inzani, and S. M. Griffin, Multi-Channel Direct Detection of Light Dark Matter: Theoretical Framework, JHEP 03, 036, arXiv:1910.08092 [hep-ph] .
  * Atzori Corona  _et al._ [2023b] M. Atzori Corona, W. M. Bonivento, M. Cadeddu, N. Cargioli, and F. Dordei, New constraint on neutrino magnetic moment and neutrino millicharge from LUX-ZEPLIN dark matter search results, Phys. Rev. D 107, 053001 (2023b), arXiv:2207.05036 [hep-ph] .


# Modeling the Age Pattern of Fertility: An Individual-Level Approach

Daniel Ciganda  Max Planck Institute for Demographic Research Statistics
Institute, University of the Republic, Montevideo, Uruguay  Nicolas Todd
UMR7206 Eco Anthropologie, Musée de l’Homme - CNRS, Paris, France

###### Abstract

Macro-level modeling is still the dominant approach in many demographic
applications because of its simplicity. Individual-level models, on the other
hand, provide a more comprehensive understanding of observed patterns;
however, their estimation using real data has remained a challenge. The
approach we introduce in this article attempts to overcome this limitation.
Using likelihood-free inference techniques, we show that it is possible to
estimate the parameters of a simple but demographically interpretable
individual-level model of the reproductive process from a set of aggregate
fertility rates. By estimating individual-level quantities from widely
available aggregate data, this approach can contribute to a better
understanding of reproductive behavior and its driving mechanisms. It also
allows for a more direct link between individual-level and population-level
processes. We illustrate our approach using data from three natural fertility
populations.

##  1 Introduction

Modeling the age distribution of fertility rates is an essential step in a
number of demographic applications. When only a tight fit to an observed
schedule is required, as in the generation of single-age rates from grouped
data, nonparametric models, typically based on splines, tend to produce the
best results. No particular price is then placed on whether each model
coefficient may be interpreted in any meaningful way. Other applications
require models that can fit the data well while also providing a well-defined,
ideally small set of parameters that represent quantities that can be
interpreted in demographic terms. This typically involves ages summarizing
relevant characteristics of the fitted curve, such as the age at which the
curve rises above 0, the age at which it peaks, or the age at which it drops
down to zero.

The interpretability of parameters is particularly relevant in a forecasting
context, where it may prove especially useful to associate a change in the
value of the parameters with an underlying behavioral process, such as
fertility postponement or the diffusion of contraceptive methods. Forecasting
each parameter’s change in time is then possible, and potentially yields
improved forecasting of the aggregate quantities of interest.

A difficulty associated with most existing parametric models is that the
interpretation of their parameters can be elusive (Hoem et al., 1981). Even in
the best cases, the relationship between mechanisms and parameters is
ambiguous and indirect. In essence, macro-level models always operate one
level above the behaviors that drive fertility.

A natural solution to this problem is to model at the individual level. This
approach involves three steps: 1) develop a model of the reproductive process;
2) use this model to generate synthetic data and compute a set of simulated
age-specific fertility rates (ASFR); 3) estimate the parameters of the model
by minimizing the distance between simulated and observed rates.

Although the building blocks have been available for some time, this approach
has not yet been systematically explored.

The basic tools and ideas required for step 1 have long been developed by some
of the pioneers of mathematical demography (Gini, 1924; Henry, 1953; Sheps et
al., 1973, notably).

These ideas were later explored in more realistic settings using simulation
methods. The objective was then the incorporation of basic characteristics of
the reproductive process like age-dependent fecundability, or to represent
change over time in fertility rates (see: Ridley and Sheps, 1966; Barrett,
1971; Le Bras, 1993). The use of simulation methods allowed researchers to
complete step 2, by generating synthetic age-specific fertility rates from
individual-level models, drawing a direct connection between behavior and
aggregate quantities.

These contributions, however, remained largely theoretical. Parameter values
where borrowed from previous studies, or calibrated through trial and error
until the simulated data would fit a given target distribution. In other
words, estimation procedures remained too fuzzy, if not arbitrary, to decide
whether individual-level models could properly be used on observed data.

Interest in these type of models dwindled after the 1980’s, just at the same
time the first statistical methods that would enable inference on complex
simulation models were starting to emerge in other disciplines Rubin (1984).

Likelihood-free inference methods, such as Approximate Bayesian Computation
(ABC), have been extensively researched since then, providing a robust
statistical framework to estimate the parameters of computational models
Beaumont (2019). What we attempt to show in the remainder of this paper is
that these advances in statistical computing, that have proved so useful in
fields such as ecology and population genetics, also provide the missing piece
to the program outlined above. Specifically, these advances allow to infer
characteristics of individual level processes from aggregate fertility rates.
More generally, the approach we advocate here provides a clear roadmap to
fully integrate behavioral mechanisms in the modeling and forecasting of
fertility trends.

##  2 Model

The simplest form of the reproductive process is the one in which the outcomes
of the process (exact timing and number of births) are not a function of
individual preferences. Louis Henry called this “natural fertility”, in an
attempt to characterize a setting in which the reproductive process is
exclusively determined by physiological factors.

A model of the reproductive process in a natural fertility setting requires
three basic inputs: The moment when the process starts, typically defined as
the age at marriage; the risk or probability of a conception, known as
fecundability; and the length of the period in which women are not able to
conceive following childbirth, know as post-partum amenorrhea.

We represent the reproductive experience of a cohort of women from birth to
age 50 using a discrete-time simulation model, in which time advances at fixed
increments of one month. The age (in months) at marriage for a woman $i$ is
simulated from a lognormal distribution with mean $\mu_{m}$ and standard
deviation $\sigma_{m}$.

Every month, married women who are neither pregnant nor amenorrheic are
exposed to the risk of a conception, which depends on the woman’s age. Whether
woman $i$ conceives in a given month when aged $x^{*}$ is the outcome of a
Bernoulli trial with probability $\phi(x^{*})$.

We model the age-pattern of fecundability using two polynomial bases:

$\phi(x)=\phi_{1}\cdot 3x_{s}(1-x_{s})^{2}+\phi_{2}\cdot 3x_{s}^{2}(1-x_{s})$
(1)

Where $\phi(x)$ is fecundability at age $x$ in the reproductive window
(assumed 10 to 50 years), and $x_{s}=\frac{x-10}{50-10}$, the age scaled to
belong to $[0,1]$. These basis functions are two of the four Bernstein
polynomials of degree 3. The other two [$x_{s}\mapsto(1-x_{s})^{3}$ and
$x_{s}\mapsto x_{s}^{3}$] have no use when representing fecundability, which
is 0 at ages 10 and 50.

Different combinations of $\phi_{1}$ and $\phi_{2}$ generate different age
profiles for the risk of conception, as seen in Figure 1. $\phi_{1}$ has a
stronger influence on the level at which fecundability peaks, while $\phi_{2}$
largely controls the pace at which the risk of conception decays, until
permanent sterility is reached. While the specific values of $\phi_{1}$ and
$\phi_{2}$ are not informative on their own, the resulting age profile of
fecundability produce by a given combination reflects both the behavioral
(such as the frequency of intercourse) and biological mechanisms that
influence the risk of conception.

Figure 1: Age Distribution of Fecundability for Different Combinations of
Parameter Values. Although it depends on only two parameters, the model used
to represent the evolution of the risk of conception across the reproductive
life course captures a wide range of patterns.

A conception brings about a state of non-susceptibility (no risk of
conception) equal to the duration of the pregnancy, 9 months, plus the
duration of post-partum amenhorrea, $\delta$, that must be estimated together
with the other parameters of the model.

The model outputs individual reproductive trajectories, from which a set of
simulated age-specific fertility rates ${}_{s}f(x)$ can be computed. We then
leverage the distance between simulated and observed rates, ${}_{o}f(x)$, to
estimate the parameters of our model, namely
$(\mu_{m},\sigma_{m},\phi_{1},\phi_{2},\delta)$. This procedure is explained
in more detail in Section 3.5, after we describe the data used in our analysis
in the next section.

##  3 Data

To illustrate our approach, we use data from three high quality natural
fertility populations: Hutterites communities from 1860 to 1914, the French-
Canadian population from 1700 - 1750, and a subsample of the Louis Henry
Enquête on the French population in the XVIIIth century. All three provide
full maternity histories as well as information on the dates of other relevant
events in a woman’s lifecourse like marriage, death, or death of spouse.

###  3.1 Hutterites

The Hutterites are a Anabaptist community, originated in the XIVth century in
the Tyrolean Alps. After a long history of migration in Europe and Asia, they
relocated to North America at the end of the XIXth century where they still
live today in self-sustained, largely isolated colonies.

Like many religious communities, the Hutterites oppose the use of birth
control methods (Lee and Brattrud, 1967; Ingoldsby and Stanton, 1988). What
made them stand out for demographic analysis was how strictly they adhered to
these beliefs to at least the second half of the 20th century. The marital
fertility of the Hutterite cohorts born until the early 1900s remained close
to the theoretical maximum, with an average of around 10 children per woman,
providing researchers with an exceptional opportunity to study the
reproductive process under natural fertility conditions (Eaton and Mayer,
1953).

Another reason why the Hutterites became the gold standard for natural
fertility research was their custom of keeping detailed family records. These
records, personally checked for consistency by colony preachers, were made
available for various scientific studies in the 1950s and 60s (Eaton and
Mayer, 1953; Mange, 1964). These data were extended in the course of one of
these studies through follow-up interviews, which resulted in complete
maternity histories for 562 families Sheps (1965).

###  3.2 French-Canadian

BALSAC is a longitudinal population database that contains information on
individuals and families who have lived in Quebec from the 17th century to the
present. The information used in this article, that concerns cohorts born
before 1750, was gathered using family reconstitution methods from parish
registers (Vézina and Bournival, 2020). Given the high quality of these
records, the information on XVII Century Quebec populations has become an
important reference in the literature on natural fertility (see: Clark et al.,
2020; Larsen and Yan, 2000; Eijkemans et al., 2014, among others).

###  3.3 Henry

Another important dataset for natural fertility research is the Enquête Louis
Henry. This dataset was created by Louis Henry and his coworkers as an effort
to draw a nationally representative sample of XVIIIth century France
population. It contains data on life events for a sample of 378 parishes for
the period 1670 to 1829 from parish registers and administrative records – the
’anonymous sample’. The families in a sub-sample of 40 rural parishes – the
’nominative sample’ – have been fully reconstituted (seguy2001, for a detailed
description of the Enquête, see).

###  3.4 Sample Sizes

To keep things simple, we model a process without death or union dissolution.
Therefore, we restrict our data to “intact” marriages, i.e., marriages that do
not dissolve by death or separation before the woman reaches age 50. Table 1
contains the remaining sample sizes after we exclude censored trajectories and
records with incomplete information. Records with partial information were
recovered in some cases by imputing the month or day of events. The inclusion
of imputed data does not significantly affects our results.

While our simulated data belongs to a single birth cohort, in order to get
reliable estimates our empirical data is obtained from multiple cohorts. To
minimize potential biases we kept a range of birth cohorts that are as
homogeneous as possible with respect to the age at first marriage and the
total fertility rate. The cohorts included for each population are also
displayed in Table 1

| Nr. Marriages | Nr. Births | Cohorts
---|---|---|---
Hutterites | 161 | 1726 | 1860 - 1914
XVIII Century Quebec | 14303 | 110772 | 1722 - 1730
XVII - XVIII Century France | 3235 | 18623 | 1680 - 1760
Table 1: Sample Information for the Three Natural Fertility Datasets

Figure 2 shows the observed age-specific fertility rates for our three natural
fertility populations. As expected, the Hutterites show consistently high
fertility from age 20, with an average number of children per woman in these
cohorts of 10.76. The Hutterite data is also unique in terms of its pattern,
with a visible hump around ages 24-26, followed by a decline that is slow
until age 35 and accelerates afterwards. The French-Canadian cohorts have the
second largest fertility in the three populations, with an average of 7.7
children per woman, followed by the French cohorts with an average of 5.75
children per woman. The shape of the ASFR is more conventional in these last
two populations, characterized by a rounded top, which in the case of the
French-Canadian cohorts is found around ages 28 and in the case of the French
cohorts around age 32. This is not surprising considering the mean ages at
marriages for these cohorts which is 20.4 for the Hutterites, 21.1 for the
French-Canadian, and 24.1 for the French cohorts

Figure 2: Age Specific Fertility Rates for Three Natural Fertility
Populations.

###  3.5 Estimation

As previously mentioned, our primary goal is to demonstrate that the lack of
individual-level data does not have to be an impediment to understanding the
behaviors and mechanisms underlying fertility processes.

When only aggregate data is available, the probability of individual
observations cannot be calculated, and model parameters cannot be estimated
using a maximum likelihood approach. Instead, we rely on Approximate Bayesian
Computation (ABC), a popular likelihood-free approach for the estimation of
computational models.

The objective of ABC is to get a sample approximately drawn from the posterior
distribution without ever explicitly calculating the likelihood function.
Instead, ABC algorithms simulate data from the model for a set of parameter
values, and keep combinations of parameter values that yield simulated data
that is close to the observed data, according to a predefined distance metric.

In the basic ABC rejection algorithm we use, parameter samples are drawn from
the prior distribution (Tavaré et al., 1997; Pritchard et al., 1999). While
this method is straightforward to implement, it requires a large number of
simulations to obtain a reasonable sample from the approximate posterior.

To improve the efficiency of the approach, the parameter values accepted in
the first step can be modeled as a function of the simulated data to obtain a
better approximation of the posterior distribution, a procedure known in the
ABC literature as “regression adjustment” (Blum, 2018). For this second step,
we use a random forest (RF), as described in the pseudocode below.

Algorithm:  ABC + Regression Adjustment

Input: Prior distribution $\pi$, model $\mathcal{M}$ for simulation, observed
data $y_{\text{obs}}$, distance threshold $\epsilon$.

Output: Adjusted samples from the posterior distribution
$\\{\tilde{\theta}_{1},\tilde{\theta}_{2},\ldots,\tilde{\theta}_{N}\\}$.

begin

Draw parameter samples, $\theta^{*}$, from $\pi$

For each $\theta^{*}$, simulate data, $y^{*}$, from $\mathcal{M}(\theta^{*})$

Compute distance between $y^{*}$ and $y_{\text{obs}}$

Accept $\theta^{*}$ if the distance is less than $\epsilon$

With accepted parameters $\\{\theta_{1},\theta_{2},\ldots,\theta_{N}\\}$ and
corresponding simulated data $\\{y_{1},y_{2},\ldots,y_{N}\\}$:

begin

Train a Random Forest $m$ that maps $y$ to $\theta$

Adjust the values of the accepted parameters:

$\tilde{\theta}_{i}=m(y_{\text{obs}})+(\theta_{i}-m(y_{i}))$

end

Return $\\{\tilde{\theta}_{1},\tilde{\theta}_{2},\ldots,\tilde{\theta}_{N}\\}$
as samples from the adjusted posterior distribution.

end

###  3.6 Sample Size

The comparison between simulated and observed data in ABC typically relies on
summary statistics, as using the raw data directly is usually impractical or
impossible.

We use the set of age-specific fertility rates for our three natural fertility
populations as our target summary statistics. As shown in Table 1, the number
of individual trajectories used to compute these rates varies widely between
the three populations. To adequately represent the uncertainty associated with
each sample size, we estimate each population separately, matching the number
of marriages in our observed data with the number of marriages used when
simulating from the model.

###  3.7 Prior choice

Given that one of our primary goals is to determine whether or not it is
possible to infer the parameters of our individual model from information
contained on a set of age-specific fertility rates, we use non informative
priors for our model’s parameters. More specifically, the results presented in
the following sections are obtained using uniform priors in a range containing
reasonably implausible values for the five parameters in our model.

##  4 Results

###  4.1 Model Validation

We use cross-validation to assess whether the parameters in our model are
correctly identified from a set of age-specific fertility rates. This
procedure makes use of the set of parameter values $\theta^{*}$ and their
corresponding simulated rates $y^{*}$. From this set we randomly select the
$i^{th}$ simulation as our validation set. The simulated rates in the
validation set are taken as pseudo-observed values. Taking into account all
simulations except the validation set, the model parameters are then estimated
using the procedure described in the previous section. This process is
repeated for a subset of 100 randomly chosen simulations and the prediction
error is computed as:

$\varepsilon_{\text{p}}=\sum_{i}\left(\frac{(\theta^{\ast}_{i}-\theta_{i})^{2}}{\text{Var}(\theta_{i})}\right)$

The results of the cross-validation exercise are shown in Figure 3(e). The
three sets of results correspond to the different sample sizes used in the
simulations. The results obtained with the largest sample size show that the
estimation procedure consistently returns parameter values $\theta$ that are
very close to the true parameters that generated the data $\theta^{*}$. As
expected, simulations results obtained with a smaller number of marriages are
associated with large preduction errors, which shows the estimation approach
correctly reflects sample size uncertainty.

The parameter with the smallest associated prediction errors is $\mu_{m}$ and
the parameter with the largest associated errors is $\delta$. This is not
unexpected as a change in the value of $\mu_{m}$ produces a shift of the
entire distribution of rates that is very distinct from the effect of a value
change in any other parameter. A change in the value of $\delta$, on the other
hand, produces an effect on the rates that overlaps to a certain extent to the
effect of a change in the value of $\phi_{1}$ (see appendix A for some
animations illustrating the effect of each parameters on the schedule of age-
specific fertility rates).

[figure]capposition=top

Figure 3: Results of Model Cross-Validation. The location of each of the red
dots displayed in the figure is defined by the mean of the estimated posterior
distribution in the $y$ axis and the value of the parameter that generated the
simulated results used as pseudo-observed data during the validation process.

(a) $\phi_{1}$

(b) $\phi_{2}$

(c) $\delta$

(d) $\mu_{m}$

(e) $\sigma_{m}$

###  4.2 Model Fit

Figure 4(c) depicts the model’s posterior predictive distribution against
observed ASFR for the three populations analyzed. The model accurately
captures the general characteristics of the data, such as the shape, location
of peaks, and rate of decline of the three distributions. In fact, the fit is
very satisfactory at most ages, except for the rates at the very end of the
reproductive age window (ages 47 to 50), which the model tends to slightly
overestimate. Credible intervals reflect the higher degree of uncertainty with
which observed rates are estimated in the Hutterite’s case.

[figure]capposition=top

Figure 4: Observed vs. Simulated Age-Specific Fertility Rates with 95%
Credible Intervals.

(a) Hutterites

(b) XVIII Century Quebec

(c) XVII - XVIII Century France

##  5 Discussion

The need to represent the age-distribution of fertility rates has
overwhelmingly been met with solutions at the aggregate level. These solutions
traded parameter interpretability in exchange for statistical formalization.
Thanks to the approach we introduced in this article, this compromise is no
longer necessary.

Our individual-level model is certainly a simple, even crude, description of
the reproductive process. Yet, it fits a wide range of fertility schedules
just as well as previous purely statistical models, with the major strength
that it explicitly models the behavioral and biological mechanism that give
rise to these schedules.

Our approach also allows for the transparent incorporation of additional
information via prior distributions when such information is available. This
provides a very favorable contrast to the trail-and-error calibration
procedures commonly used in micro-level models.

We believe these ideas have the potential to make contributions in at least
two fields. In biodemography, e.g. in the study of non-human primates, our
approach could be used to classify species in a meaningful way, by age at
onset of sexual reproduction, postpartum amenorrhea duration and
fecundability, rather than abstract properties of asfr curves.

In demography, we hope it inspires the exploration of more complex models with
the aim of understanding the microfoundations of present and future fertility
trends.

##  Appendix A Appendix

Figure 5: Effect of $\phi_{1}$ on Simulated Age-Specific Fertility Rates

[loop, controls=true,width=10cm]6plots/phibeta1/plotseq180

Figure 6: Effect of $\phi_{2}$ on Simulated Age-Specific Fertility Rates

[loop, controls=true,width=10cm]6plots/phibeta2/plotseq180

Figure 7: Effect of $\mu_{m}$ on Simulated Age-Specific Fertility Rates

[loop,controls=true,width=10cm]6plots/mum/plotseq180

Figure 8: Effect of $\sigma_{m}$ on Simulated Age-Specific Fertility Rates

[loop,controls=true,width=10cm]6plots/sigmam/plotseq180

Figure 9: Effect of $\delta$ on Simulated Age-Specific Fertility Rates

[loop, controls=true,width=10cm]6plots/nsp/plotseq180

##  Appendix B Appendix

[figure]capposition=bottom

| $\phi_{1}$ | $\phi_{2}$ | $\delta$
---|---|---|---
Hutterites |  |  |
French Canadian |  |  |
Historic France |  |  |
Figure 10: Adjusted Posterior Distributions

[figure]capposition=bottom

| $\mu_{m}$ | $\sigma_{m}$
---|---|---
Hutterites |  |
French Canadian |  |
Historic France |  |
Figure 11: Adjusted Posterior Distributions

## References

  * Barrett (1971) Barrett, J. C. (1971).  Use of a fertility simulation model to refine measurement techniques.  Demography 8(4), 481–490.
  * Beaumont (2019) Beaumont, M. A. (2019).  Approximate bayesian computation.  Annual review of statistics and its application 6, 379–403.
  * Blum (2018) Blum, M. G. (2018).  Regression approaches for abc.  Handbook of approximate Bayesian computation, 71–85.
  * Clark et al. (2020) Clark, G., N. Cummins, and M. Curtis (2020).  Twins support the absence of parity-dependent fertility control in pretransition populations.  Demography 57(4), 1571–1595.
  * Eaton and Mayer (1953) Eaton, J. W. and A. J. Mayer (1953).  The social biology of very high fertility among the hutterites. the demography of a unique population.  Human biology 25(3), 206.
  * Eijkemans et al. (2014) Eijkemans, M. J., F. Van Poppel, D. F. Habbema, K. R. Smith, H. Leridon, and E. R. Te Velde (2014).  Too old to have children? lessons from natural fertility populations.  Human Reproduction 29(6), 1304–1312.
  * Gini (1924) Gini, C. (1924).  Premières recherches sur la fécondabilité de la femme.  In North-Holland: (Ed.), Proceedings of the International Mathematical Congress., Volume Vol. 2., Toronto.
  * Henry (1953) Henry, L. (1953).  Fondements théoriques des mesures de la fécondité naturelle.  Revue de l’Institut International de Statistique / Review of the International Statistical Institute 21(3), 135–151.
  * Hoem et al. (1981) Hoem, J. M., D. Madsen, J. L. Nielsen, E.-M. Ohlsen, H. O. Hansen, and B. Rennermalm (1981).  Experiments in modelling recent danish fertility curves.  Demography 18(2), 231–244.
  * Ingoldsby and Stanton (1988) Ingoldsby, B. B. and M. E. Stanton (1988).  The hutterites and fertility control.  Journal of Comparative Family Studies 19(1), 137–142.
  * Larsen and Yan (2000) Larsen, U. and S. Yan (2000).  The age pattern of fecundability: an analysis of french canadian and hutterite birth histories.  Social biology 47(1-2), 34–50.
  * Le Bras (1993) Le Bras, H. (1993).  Simulation of change to validate demographic analysis. Oxford England Clarendon Press 1993.
  * Lee and Brattrud (1967) Lee, S. and A. Brattrud (1967).  Marriage under a monastic mode of life: A preliminary report on the hutterite family in south dakota.  Journal of Marriage and the Family, 512–520.
  * Mange (1964) Mange, A. P. (1964).  Growth and inbreeding of a human isolate.  Human Biology 36(2), 104–133.
  * Pritchard et al. (1999) Pritchard, J. K., M. T. Seielstad, A. Perez-Lezaun, and M. W. Feldman (1999).  Population growth of human y chromosomes: a study of y chromosome microsatellites.  Molecular biology and evolution 16(12), 1791–1798.
  * Ridley and Sheps (1966) Ridley, J. C. and M. C. Sheps (1966).  An analytic simulation model of human reproduction with demographic and biological components.  Population Studies 19(3), 297–310.
  * Rubin (1984) Rubin, D. B. (1984).  Bayesianly justifiable and relevant frequency calculations for the applied statistician.  The Annals of Statistics, 1151–1172.
  * Sheps (1965) Sheps, M. C. (1965).  An analysis of reproductive patterns in an american isolate.  Population studies 19(1), 65–80.
  * Sheps et al. (1973) Sheps, M. C., J. A. Menken, and A. P. Radick (1973).  Mathematical models of conception and birth.  University of Chicago Press Chicago.
  * Tavaré et al. (1997) Tavaré, S., D. J. Balding, R. C. Griffiths, and P. Donnelly (1997).  Inferring coalescence times from dna sequence data.  Genetics 145(2), 505–518.
  * Vézina and Bournival (2020) Vézina, H. and J.-S. Bournival (2020).  An overview of the balsac population database: past developments, current state and future prospects.  Historical Life Course Studies 9, 114–129.