# Influence of the gravitational darkening effect on the spectrum of a hot,
rapidly rotating neutron star

Agnieszka Majczyna National Centre for Nuclear Research, ul. Andrzeja Sołtana
7, 05-400 Otwock, Poland;  Jerzy Madej  Astronomical Observatory, University
of Warsaw, Al. Ujazdowskie 4, 00-478 Warszawa, Poland;  Agata Różańska
Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, ul.
Bartycka 18, 00-716 Warszawa, Poland;  Mirosław Należyty  Astronomical
Observatory, University of Warsaw, Al. Ujazdowskie 4, 00-478 Warszawa, Poland;

###### Abstract

In this paper, we discuss the influence of the gravitational darkening effect
on the emergent spectrum of a fast-rotating, flattened neutron star. Model
atmosphere codes always calculate spectra of emergent intensities and fluxes
emitted from the unit surface on the star in plane-parallel geometry. Here we
took a step beyond that and calculated a small sample grid of theoretical
spectra integrated over the distorted surface of a sample rotating neutron
star seen by a distant observer at various inclination angles. We assumed
parameters like two dimensionless angular velocities $\bar{\Omega}^{2}=0.30$
and 0.60, the effective temperature of a nonrotating star $T_{\rm
eff}=2.20\times 10^{7}\,$K, the logarithm of the surface gravity of a
spherical star $\log(g)=14.40$ (cgs), and inclination angles from
$i=0^{\circ}$ to $i=90^{\circ}$ with step $\Delta i=10^{\circ}$. We assumed
that the atmosphere consists of a mixture of hydrogen and helium with $M_{\rm
H}=0.70$ and $M_{\rm He}=0.30$. At each point on the neutron star surface, we
calculated true intensities for local values of parameters ($T_{\rm eff}$ and
$\log(g)$), and these monochromatic intensities are next integrated over the
whole surface to obtain the emergent spectrum. In this paper, we compute for
the first time theoretical spectra of the fast-rotating neutron star. Our work
clearly shows that the gravitational darkening effect strongly influences the
spectrum and should be included in realistic models of the atmospheres of
rotating neutron stars.

radiative transfer – gravitational darkening effect – stars: neutron

††journal: ApJ

##  1 Introduction

A star distorted by either the centrifugal force caused by rotation or by
tidal forces from a secondary star in a close binary system radiates away
radiative energy in a nonsymmetric way and exhibits variations of the
radiation power over its surface. Consequently, the effective temperature on
the stellar surface is not uniform. The effect was recognized by Hugo von
Zeipel in 1924 (von Zeipel, 1924). He investigated the radiative stellar
interior in hydrostatic equilibrium. The main conclusion of his work was that
the locally emergent bolometric flux of the radiation of a distorted star is
proportional to the local gravity. Assuming of Stefan-Boltzmann law, it could
be written as $T_{\rm eff}\sim g^{1/4}$, known as the von Zeipel law. The
results of von Zeipel are valid only for stars with fully radiative envelopes
and hydrostatic equilibrium. This law must be written in the more general form
as $T_{\rm eff}\propto g_{\rm eff}^{\beta}$, where $\beta$ is the
gravitational darkening exponent (GDE). If the envelope of the star is not
fully radiative, the gravitational darkening exponent differs from 1/4.

The literature on the theory of stellar atmospheres involves the following
terms: the gravitational darkening and the limb darkening. We direct the
attention of the reader to the basic difference between these two terms. Limb
darkening is an optical effect seen both in spherical and distorted stars
where the central part of the star’s disk is brighter than the star’s limb.
This effect denotes that the specific intensity of the radiation emitted at
some point on the surface is the largest in the vertical direction and
decreases in directions approaching the horizontal plane for a single value of
the effective temperature at this point. Gravitational darkening denotes that
the value of the effective temperature (total luminosity) varies over the
surface of the distorted star. The distortion of the star can be caused either
by centrifugal forces due to rotation or by tidal interactions between two
stars in a close binary system.

Gravitational darkening effects were studied for different aspects and for
different types of objects (see e.g. Djurašević et al., 2003, 2006; Domiciano
de Souza et al., 2014). Theoretical models of gravitational darkening for
different types of objects were also developed.Lucy (1967) investigated a
fully convective envelope and found that $\beta\sim 0.08$. This is a very weak
dependence of the effective temperature on the effective gravity. The results
of Lucy (1967) research were confirmed by Webbink (1976), who used a simple
opacity law. The author presented an approximation, which under the assumption
of constant entropy of the envelope, gave results consistent with those
obtained by Lucy (1967). Note, that all these above results were obtained
under assumptions of rigid and small rate of rotation. In the case of neutron
stars, the latter assumption is not strictly valid.

Recently Claret (2015) investigated gravitational darkening for neutron stars
with various rotational models assumed. The author assumed a spherical shape
of the neutron star and explored the dependence of $\beta$ on the temperature
of the star. There is a weak dependence of the gravitational darkening
exponent on the temperature below 6$\times 10^{7}$ K in all assumed models of
the rotation. In the above paper, the author was concerned about the
calculation of the gravitational darkening exponent and its dependence on the
rotation law. The author did not calculate any model spectra of a rotating
star, either normal or neutron star.

The equilibrium configurations of a slowly rotating neutron star for different
central densities $\rho_{\rm c}$ and circular angular velocities $\Omega$ were
investigated by Belvedere et al. (2014). The most interesting result of this
paper was presented in their Figure 16, where the authors have shown the
dependence of the eccentricity on the frequency of the neutron star rotation.
What is important, for frequency typical for neutron stars in X-ray bursters
the eccentricity is of the order of 0.6, which means that the ratio of the
polar radius to the equatorial radius $R_{\rm p}/R_{\rm eq}\sim 0.8$. We
expect that such flattening should give rise to visible changes in the
emerging spectrum.

The effective surface acceleration of the rotating neutron star was calculated
by AlGendy & Morsink (2014) by using useful approximate formulae. On the
surface of an oblate star, the effective gravity acceleration is decreased by
the centrifugal force at the equator and increased by the centrifugal
flattening of the star. The authors have shown useful formulae on the
effective gravity for the approxi mation of fast and slow rotation. These
formulae weakly depend on the chosen equation of state of the neutron star
matter. In the slow rotation approximation, the effective gravity depends only
on the compactness ratio $M/R_{\rm eq}$, the dimensionless square of the
angular velocity $\Omega^{2}R_{\rm eq}^{3}/GM$, and the latitude on the
surface of the star. For spin frequency $\sim$600 Hz difference of the
effective gravity between the pole and equator is about 20%. This effect is
important not only for the spectrum calculation but, also for it influence on
the Eddington luminosity of a neutron star. AlGendy & Morsink (2014) have
shown that for the assumed $\bar{\Omega^{2}}=0.1$, the expected difference
between the Eddington luminosity at the pole and at the equator is about 10%.
This is crucial because X-ray bursts at the Eddington limits (PRE \-
photospheric radius expansion bursts) are used as the standard candles
(Galloway et al., 2008).

Suleimanov et al. (2020) investigated the spectra of slowly rotating neutron
stars for a few models of the emission, from the simplest blackbody to their
theoretical model atmospheres. Using their models, the authors calculated the
dilution factor $w$ and the spectral hardening factor $f_{c}$ to obtain local
specific intensities, which are next integrated over the surface of the star.
Different rotational frequencies and inclination angles were investigated. The
shape of the star and the local gravity are calculated under the approximation
of slow rotation. General relativistic effects were also included.

We note here that Jaroszynski (1986) modeled the external gravitational field
and surface shape of neutron stars rotating with the single angular velocity
$\Omega=0.27GM/R^{3}$. Jaroszynski (1986) investigated the influence of the
rotation on the observed luminosity and spectrum of the rotating neutron star.
For the assumed angular velocity, eccentricity, and other basic parameters of
the star, spectra have been calculated. Even if the blackbody model was used
for the spectrum, the light bending was included. The ray-tracing procedure
was included by the author. The main conclusion of Jaroszynski (1986) was that
a two-temperature blackbody shape is needed to fully describe the spectrum of
the rotating star. The author also showed that a small, 10%, changing
eccentricity cannot be seen in the emerging spectrum.

In this paper, we show how fast rotation and therefore the gravitational
darkening effect influences the shape of the spectrum of the oblate neutron
star. Emergent spectra of the rotating neutron star seen by a distant observer
are calculated by the ATM24 code, which solves the full radiative transfer
including Compton scattering (see e.g. Majczyna et al., 2005; Madej et al.,
2017). We calculated models with different parameters like an eccentricity,
angular velocities, or inclination angle. We show that these spectra differ
from each other and from the spectra of the spherical star. Therefore, the
gravitational darkening effect should be included in realistic models of the
neutron star atmospheres.

##  2 Model atmospheres

###  2.1 Gravitational darkening

To calculate model spectra including the gravitational darkening effect, we
need to define the geometry of the problem, calculate the effective gravity
over the surface of the oblate star, and the intensity spectra at each point
on the star. We used the spherical coordinate system where $\theta=0^{\circ}$
denotes the latitude of the equator.

In the first approximation, we assume that the shape of the rotating neutron
star could be described as an ellipsoid. The radius of this solid is

$R(\theta)=R_{\rm eq}\sqrt{1-e^{2}\cos(90^{\circ}-\theta)},$ (1)

where $\theta$ is the latitude in spherical coordinates, $R_{\rm eq}$ is the
equatorial radius, whereas $R_{\rm p}$ is the polar radius, and $e$ is
eccentricity, defined as $e=\sqrt{(R_{\rm eq}^{2}-R_{\rm p}^{2})/R_{\rm
eq}^{2}}$.

The effective gravity on the surface of the oblate star depends on the
latitude angle. To describe this dependence, we use the approximate formula
given by AlGendy & Morsink (2014), Eq. 50 in their work, in the regime of the
fast rotation. In the original equation, the colatitude angle was used,
whereas in the present paper, we use latitude angle $\theta$, therefore, we
replace $\theta^{\prime}=90^{\circ}-\theta$

$\displaystyle g(\theta)/g_{\rm 0}$ $\displaystyle=$ $\displaystyle 1+(c_{\rm
e}\bar{\Omega^{2}}+d_{\rm e}\bar{\Omega^{4}}+f_{\rm
e}\bar{\Omega^{6}})\sin^{2}(90^{\circ}-\theta)$ (2) $\displaystyle+$
$\displaystyle(c_{\rm p}\bar{\Omega^{2}}+d_{\rm p}\bar{\Omega^{4}}+f_{\rm
p}\bar{\Omega^{6}}-d_{\rm
60}\bar{\Omega^{4}})\cos^{2}(90^{\circ}-\theta)+d_{\rm
60}\bar{\Omega^{4}}\cos(90^{\circ}-\theta),$

where $g_{\rm 0}$ is the gravity of a nonrotating star, and the following
coefficients with $x=M/R_{\rm eq}$ are

$\displaystyle c_{\rm e}$ $\displaystyle=$ $\displaystyle-0.791+0.776x,$
$\displaystyle d_{\rm e}$ $\displaystyle=$ $\displaystyle-1.315x+2.431x^{2},$
$\displaystyle f_{\rm e}$ $\displaystyle=$ $\displaystyle-1.172x,$
$\displaystyle c_{\rm p}$ $\displaystyle=$ $\displaystyle 1.138-1.431x,$
$\displaystyle d_{\rm p}$ $\displaystyle=$ $\displaystyle 0.653x-2.864x^{2},$
$\displaystyle f_{\rm p}$ $\displaystyle=$ $\displaystyle 0.975x,$
$\displaystyle d_{60}$ $\displaystyle=$ $\displaystyle 13.67x-27.13x^{2}.$

The dimensionless angular velocity $\bar{\Omega}$ is defined as

$\bar{\Omega}=\Omega\left(\frac{R_{\rm eq}^{3}}{GM}\right)^{1/2}.$ (4)

Figure 1: Distribution of the effective temperature over the surface of a
neutron star for different dimensionless angular velocities:
$\bar{\Omega}^{2}=0.10$ (left panel), $\bar{\Omega}^{2}=0.40$ (middle panel),
and $\bar{\Omega}^{2}=0.60$ (right panel). We assumed the effective
temperature and surface gravity of an undistorted star: $T_{\rm
eff}=2.20\times 10^{7}\,$K and $\log(g)=14.40\,$(cgs), respectively, and
Gravitational Darkening Exponent (GDE) $\beta=0.25$, the eccentricity of the
rotating star $e=0.70$ and inclination angle $i=30^{\circ}$.

The local effective gravity and therefore the local effective temperature are
strong functions of the dimensionless angular velocity, which we have shown in
the example of the distribution of the local temperature in Fig. 1. In the
figure, we use the same scale factor for converting values of the temperature
into the gray scale. For low values of the $\bar{\Omega}$, the local effective
gravity or temperature varies with the latitude angle $\theta$ much slowly
than for higher values of $\bar{\Omega}$.

Figure 2: Dependence of the effective gravity (left panel) and the effective
temperature (right panel) on the latitude angle $\theta$ for various
dimensionless angular velocities. We assume the logarithm of the surface
gravity of a spherical star $\log(g)=14.40$(cgs), the effective temperature of
the spherical star $T_{\rm eff}=2.20\times 10^{7}\,$K, GDE $\beta=0.25$ and
various dimensionless angular velocities $\bar{\Omega}^{2}=0.10$ (blue line),
0.40 (red line), and 0.60 (black line).

Figure 2 shows the dependence of the effective surface gravity on the latitude
angle. We assume the gravity of a spherical star $\log(g)=14.40\,$(cgs). As
was shown in the Fig. 2, the surface gravity on the pole is not equal to the
assumed surface gravity of a spherical star $\log(g)=14.40$. In our
calculations, we include the quadrupole moment caused by the mass distribution
in the neutron star interior. This moment does not vanish on the pole and
causes greater gravity than the gravity of a nonrotating neutron star.

Figure 3: Dependence of the effective temperature on the latitude angle
$\theta$ for various GDEs: $\beta=0.25$ denoted by the black line,
$\beta=0.15$ by the red one, and $\beta=0.10$ by the blue one. We assume the
effective temperature of a nonrotating, spherical star, $T_{\rm
eff}=2.20\times 10^{7}\,K$, and the logarithm of the surface gravity
$\log(g)=14.40$(cgs) and dimensionless angular velocity
$\bar{\Omega}^{2}=0.60$.

The temperature at a given point on the surface of the neutron star is as
follows:

$T(\theta)=T_{\rm eff}\left(\frac{g(\theta)}{g_{0}}\right)^{1/4},$ (5)

where $T_{\rm eff}$ is the effective temperature of the spherical, nonrotating
neutron star, $g_{0}$ is the surface gravity of the spherical star, and
$g(\theta)$ is the surface gravity at a given point of the distorted star
defined by Eq.2.

In the above Eq.5, GDE is equal to 0.25 because we assumed that the atmosphere
is radiative. As was shown by (e.g., Lucy, 1967), this exponent is different
for nonradiative or not fully radiative envelopes. Very recently Claret (2021)
calculated self-consistently the GDE for DA and DB white dwarfs. The author
showed that this exponent depends not only on the temperature but also on the
opacity. In our work we did not calculate GDE. To show how different
parameters influence the emergent spectrum, we assume values of the GDE and
other parameters. In the Fig. 3, we presented the distribution of the
effective temperature over the surface of the oblate star for various GDEs. We
assume the dimensionless angular velocity $\bar{\Omega}^{2}=0.60$ and the
ratio of the mass to the equatorial radius of the neutron star $x=0.195$ and
$\beta=0.10$, $0.15$, and $0.25$.

###  2.2 ATM24

We calculate intensity spectra with the ATM24 code. Our model was described in
several papers, e.g., Madej (1991); Madej et al. (2004); Majczyna et al.
(2005). The accuracy of our code was recently improved (Madej et al., 2017;
Vincent et al., 2018).

Our code solves the radiative transfer equation in the plane-parallel
geometry. In hot atmospheres, Compton scattering plays a crucial role in
forming the emerging spectrum and therefore should be treated very carefully.
We investigated the scattering of X-ray photons with energies approaching the
electron rest mass energy on free, thermal, relativistic electrons. We allow
for any large energy exchange between photons and electrons during a single
scattering. We assume the equation of state of ideal gas in local
thermodynamical equilibrium (LTE). Note that the Compton scattering provides
non-LTE effects in the model atmospheres. When calculating opacities, we
include energy-dependent free-free and bound-free opacities form all elements
with Z$\leq$30 (atomic number) and bound-bound opacities from hydrogen,
helium, and iron.

The equation of transfer was adopted from Pomraning (1973, see also, ) and has
the following form:

$\displaystyle\mu\,\frac{dI_{\nu}}{d\tau_{\nu}}=I_{\nu}-{\frac{k_{\nu}}{k_{\nu}+\sigma_{\nu}}}B_{\nu}-\left(1-{\frac{k_{\nu}}{k_{\nu}+\sigma_{\nu}}}\right)J_{\nu}$
(6) $\displaystyle+$
$\displaystyle\left(1-{\frac{k_{\nu}}{k_{\nu}+\sigma_{\nu}}}\right)J_{\nu}\int\limits_{0}^{\infty}\Phi(\nu,\nu^{\prime})\left(1+{\frac{c^{2}}{2h{\nu^{\prime}}^{3}}}J_{\nu^{\prime}}\right)d\nu^{\prime}+$
$\displaystyle-$
$\displaystyle{\frac{k_{\nu}}{k_{\nu}+\sigma_{\nu}}}\left(1+{\frac{c^{2}}{2h\nu^{3}}}J_{\nu}\right)$
$\displaystyle\times$
$\displaystyle\int\limits_{0}^{\infty}\Phi(\nu,\nu^{\prime})J_{\nu^{\prime}}{\left({\frac{\nu}{\nu^{\prime}}}\right)}^{3}\exp\left[-{\frac{h(\nu-\nu^{\prime})}{{\rm
k}T}}\right]d\nu^{\prime},$

where $k_{\nu}$ and $\sigma_{\nu}$ denote coefficients of absorption and
electron scattering, respectively. $I_{\nu}$ is the energy-dependent specific
intensity, $J_{\nu}$ is the mean intensity of radiation, and $z$ is the
geometrical depth in the atmosphere. We formulated the exact angle-averaged
redistribution function $\Phi(\nu,\nu^{\prime})$ and the Compton scattering
cross section $\sigma(\nu\to\nu^{\prime},\vec{n}\cdot\vec{n^{\prime}})$ (for
details, see Madej et al. (2017); Nagirner & Poutanen (1993, 1994)). The
equation of transfer is solved together with the hydrostatic and radiative
equilibrium equations because we assume that the atmosphere is static and that
only photons transport the energy accumulated in the core of the neutron star.
We neglect the influence of the magnetic field, accretion, and the effect of
electron degeneracy, which is unimportant in the hot atmospheres. The limb-
darkening/-brightening effects are fully taken into account in our
computations.

Figure 4: Theoretical spectra of the unit surface calculated by the ATM24 code
for different values of the effective temperatures and surface gravities. In
the left panel, we show the theoretical spectra for the effective temperatures
$T_{\rm eff}=2.20\times 10^{7}\,$K and various surface gravities from
$\log(g)=14.30\,$(cgs) to $\log(g)=14.50\,$(cgs). In the right panel we show
spectra for the fixed surface gravity $\log(g)=14.40\,$(cgs) and various
effective temperatures from $T_{\rm eff}=2.16\times 10^{7}\,$K to $T_{\rm
eff}=2.26\times 10^{7}\,$K.

We assumed that the rotating neutron star is seen at some inclination to the
equatorial plane. Input parameters of the ATM24 code are the effective
temperature, surface gravity, and the chemical composition. We assumed that
the reference nonrotating star has the effective temperature $T_{\rm
eff}=2.20\times 10^{7}\,$K and the surface gravity $\log(g)=14.40\,$(cgs). The
atmosphere consists of hydrogen and helium with the cosmic mass abundance
$M_{\rm H}=0.70$ and $M_{\rm He}=0.30$.

In X-ray bursters, especially during the burst and soon after, we expect high
temperatures of the gas, so $T_{\rm eff}=2.20\times 10^{7}\,$K is a reasonable
assumption. Values of the effective temperature and the logarithm of the
surface gravity of the spherical star were taken arbitrarily. The chemical
composition of X-ray burst sources is very poorly constrained. It is
frequently assumed that the atmosphere consists of hydrogen and helium in
arbitrarily chosen proportions. For example, Galloway et al. (2006)
investigated Eddington-limited bursts from 4U 1636-536. For this source, a
hydrogen-rich atmosphere with $M_{\rm H}=0.70$ as well as a pure helium
atmosphere were investigated. Therefore, we assume chemical composition with
cosmic hydrogen abundance.

The rotating star is oblate, therefore the gravity and the effective
temperature are different at various latitude angles. We calculated the local
effective gravity and the local effective temperature according to the
formulae given by Equations 2 and 5. These local values of the gravity and the
effective temperature are input parameters for the calculation of
monochromatic specific intensities at various latitude angles. These
intensities were next integrated over the whole surface to obtain the spectrum
as seen by a distant observer.

For a given eccentricity of the distorted star, we were able to compute the
distance from the center of gravity for each patch on the surface and its area
(in cm2), its latitude, and the azimuth as seen by the observer.
Simultaneously, we determined the cosine of the projection angle as seen by
the observer. We determined the effective gravity as a function of the
latitude. Therefore, we know all variables that are necessary to integrate the
emergent intensity over various patches to obtain the final monochromatic
luminosity of the distorted star.

In our ATM24 code, we include the limb-darkening effect and now we add the
gravitational darkening effect. An example of flux spectra of the spherical
star calculated by the ATM24 code is presented in Figure 4. In the left panel,
we compare the spectra for the same chemical compositions given above and the
same effective temperatures but different surface gravities from
$\log(g)=14.30$ up to 14.50 (cgs). With increasing surface gravity, spectra
become softer, and their maximum shifts toward lower energies. The right panel
shows our model spectra calculated for the same surface gravity and different
effective temperatures from $T_{\rm eff}=2.16\times 10^{7}$ up to $2.24\times
10^{7}\,$K. With increasing temperature, the maximum of the spectra shifts
toward higher energies.

##  3 Spectra of the rotating neutron star

We computed a small grid of theoretical spectra of the rotating neutron stars
(56 model spectra). Because the calculations are very time consuming, we
divided the star into 18x36 points in latitude and azimuthal angles. As was
mentioned above, we assumed parameters of the rotating star like the effective
temperature of the spherical star $T_{\rm eff}=2.20\times 10^{7}\,$K, surface
gravity of the nonrotating star $\log(g)=14.40\,$(cgs), chemical composition
($M_{\rm H}=0.70$ and $M_{\rm He}=0.30$) and mass to radius ratio $x=0.195$.
We computed our theoretical spectra for different values of the other
parameters: the dimensionless angular velocity $\bar{\Omega}^{2}=0.30$ and
$0.60$, and few inclination angles from $i=0^{\circ}$ up to $90^{\circ}$ with
step $\Delta i=10^{\circ}$. We note that the both assumed values of the
dimensionless angular velocity are relevant for the fast-rotation limit.

The eccentricity of the rotating neutron star depends on the angular velocity
of the rotation. The formula that connects these two variables depends on the
equation of state of the matter building up neutron star. Such approximate
formulae were proposed by various authors (see e.g., AlGendy & Morsink, 2014;
Silva et al., 2021). In our paper, we used the relation from Silva et al.
(2021).

$e=0.251+0.935\bar{\Omega}^{2}+0.709x+0.030\bar{\Omega}^{2}x-0.472\bar{\Omega}^{4}-2.427x^{2}$
(7)

Here, we quoted Eq.19 from Silva et al. (2021) using variables defined in our
paper. For the assumed $x=0.195$ and $\bar{\Omega}^{2}=0.30$ and $0.60$, the
eccentricity of the flattened star is equal to $e=0.54$ and 0.69,
respectively. These values of the eccentricity were used during the
calculations of theoretical spectra of rotating neutron star. In some cases
where we show the influence of the shape of the star on the emergent spectrum,
we thread the eccentricity as a free parameter, independent of the angular
velocity (see Fig. 5).

Assuming $\bar{\Omega}^{2}=0.30$, and 0.60, $x=0.195$ and two representative
equations of state (AlGendy & Morsink, 2014) calculated the frequency of the
neutron star rotation (see Tables 2 and 3 in their paper). For both
dimensionless angular velocities they obtained high frequencies, e.g., for
$\bar{\Omega}^{2}=0.30$, $M=1.80\,M_{\odot}$ and $R=13.6\,$km (EOS HLPS2) the
frequency is equal to $\nu=842\,$Hz whereas for $\bar{\Omega}^{2}=0.60$,
$M=1.99\,M_{\odot}$ and $R=15.1\,$km (EOS HLPS2), $\nu=1084\,$Hz. The highest
measured frequency of the rotation is equal to 716 Hz in the case of pulsar
PSR J1748-2446ad (Hessels et al., 2006). For X-ray bursters, the fastest spin
of the neutron star equals to $\nu=620\,$Hz for 4U 1608-522 (Muno et al.,
2002). The maximal value of the frequency of the neutron star rotation depends
on the equation of state, and it is as high as 1200 Hz (see e.g., Haensel et
al., 2007).

Figure 5: The top panel shows the theoretical spectra of the flattened neutron
star for a given parameters: $T_{\rm eff}=2.20\times 10^{7}\,$K,
$\log(g)=14.40\,$(cgs), $i=30^{\circ}$, $\bar{\Omega}^{2}=0.60$ and different
values of the eccentricity of the star $e=0.10$ (red line), 0.20 (magenta),
0.30 (blue), and 0.69 (black). The lower panel shows the ratio of the
monochromatic flux ($\mathcal{F}_{0.30}$) emitted by the ellipsoidal star with
$e=0.30$ and flux ($\mathcal{F}_{0.69}$) emitted by a star with $e=0.69$.

Figure 5 (top panel) shows theoretical spectra of the rotating star with
different eccentricities for the same inclination angle $i=30^{\circ}$. We
note that these spectra have very similar shapes. It may appear that these
spectra differ only in the normalization. For this reason, in the lower panel
of the Figure 5 we show the ratio of the monochromatic flux emitted by the
star with $e=0.30$ and the flux of the star with $e=0.69$. If these spectra
differ only by the normalization, this ratio should be constant across the
whole range of photon energies, which is not the case presented in the figure.
Of course, for lower differences of the eccentricity, the ratio is closer to
constant so the spectra have more similar shapes. It is especially important
for methods of mass and radius determination that used the distance to the
source (see e.g., the review by Bhattacharyya, 2010). Our mass and radius
determination method (e.g., Majczyna & Madej, 2005) is independent of the
distance; therefore value of the obtained normalization does not affect our
results.

Figure 6: Distribution of the local effective temperature over the surface of
the neutron star seen by a distant observer. We assume parameters of our
reference model for $T_{\rm eff}$ and $\log(g)$ and the gravitational
darkening exponent $\beta=0.25$, dimensionless angular velocity
$\bar{\Omega}^{2}=0.60$ and different inclination angles $i=10^{\circ}$ (left
panel), $i=40^{\circ}$ (middle panel), and $i=80^{\circ}$ (right panel).

Figure 6 shows the temperature distribution on the star seen at different
inclination angles $i=30^{\circ}$ (left panel), $i=60^{\circ}$ (middle), and
$i=80^{\circ}$ (right panel). The last panel of this figure presents the star
seen almost at the pole on, where the temperature has the highest value.

Figure 7: Theoretical spectra of the rotating neutron star for various
inclination angles. We assumed parameters $T_{\rm eff}=2.20\times 10^{7}\,$K,
$\log(g)=14.40\,$(cgs), $\beta=0.25$, and two values of the dimensionless
angular velocity $\bar{\Omega}^{2}=0.30$ (left panel) and
$\bar{\Omega}^{2}=0.60$ (right panel).

The spectral shape significantly depends on the inclination angle as was shown
in the Figure 7. An inclination angle equal to $i=90^{\circ}$ denotes that an
observer sees the star at the pole, whereas $i=0^{\circ}$ is in the equatorial
plane. As the inclination angle increases the observer sees hotter and hotter
areas, but the spectrum does not become harder and its maximum shifts toward
lower energies. This slightly unintuitive effect is connected with the
different contributions of spots with different temperatures when the
inclination angle changes. Hotter areas located at the poles are smaller than
those colder ones located at the equatorial plane.

Figure 8: Theoretical spectra of the rotating star for different effective
temperatures of the spherical neutron star. We fixed following the parameters:
$\log(g)=14.40$, $\beta=0.25$, $i=30^{\circ}$, and $\bar{\Omega}^{2}=0.30$;
$e=0.54$ (left panel) and $\bar{\Omega}^{2}=0.60$; $e=0.69$ (right panel). The
black solid curve denotes spectrum of the oblate star calculated for assumed
temperatures of the spherical star $T_{\rm eff}=2.20\times 10^{7}\,$K whereas
blue ones assume $T_{\rm eff}=1.60\times 10^{7}\,$K and red $T_{\rm
eff}=1.00\times 10^{7}\,$K. Black, blue, and red dashed curves denote spectra
of the spherical stars with appropriate temperatures.

Figure 8 shows spectra emitted by the rotating star calculated for various
effective temperatures of the spherical star. The left panel shows spectra of
the neutron star rotating with $\bar{\Omega}^{2}=0.30$, whereas the right one
is rotating with $\bar{\Omega}^{2}=0.60$. In both cases, we assume the
logarithm of gravity of the nonrotating star $\log(g)=14.40\,$(cgs). The
dashed line denotes spectra of the spherical star. This star has various
values of temperature ($T_{\rm eff}=1.0\times 10^{7}\,$K, $1.6\times
10^{7}\,$K and $2.2\times 10^{7}\,$K respectively) constant over the surface.
The maxima of spectra with gravitational darkening effect are shifted toward
higher energies, and spectra are harder than these emitted by the spherical
star.

Table 1: Bolometric luminosities and effective temperatures of the rotating,
flattened neutron star with two dimensionless angular velocities
$\bar{\Omega}^{2}$ and a few inclination angles $i$. We assumed the effective
temperature of the spherical star $T_{\rm sph}=2.20\times 10^{7}\,$K. The
bolometric luminosity of the spherical star is equal to $L_{\rm
sph}=8.179\times 10^{37}\,{\rm erg}\,{\rm s}^{-1}$ and color temperature
$T_{c}=3.31\times 10^{7}\,$K. We assume that the equatorial radius of the
oblate star and the radius of the spherical star is equal to $R=7\,$km.
$\bar{\Omega}^{2}=0.30$ ($e=0.54$)
---
Model | $L_{\rm bol}\times 10^{37}$ | $L_{\rm bol}/L_{\rm sph}$ |  $T_{\rm eff}\times 10^{7}\,$K | $T_{\rm eff}/T_{\rm sph}$ |  $T_{c}\times 10^{7}\,$K | $f_{\rm c}$
$i=0.0$ | 6.365 | 0.78 | 2.1165 | 0.962 | 3.31199 | 1.565
$i=30$ | 7.103 | 0.87 | 2.1753 | 0.989 | 3.31199 | 1.523
$i=90$ | 9.159 | 1.12 | 2.3181 | 1.054 | 3.36002 | 1.449
$\bar{\Omega}^{2}=0.60$ ($e=0.69$)
Model | $L_{\rm bol}\times 10^{37}$ | $L_{\rm bol}/L_{\rm sph}$ |  $T_{\rm eff}\times 10^{7}\,$K | $T_{\rm eff}/T_{\rm sph}$ |  $T_{c}\times 10^{7}\,$K | $f_{\rm c}$
$i=0.0$ | 4.805 | 0.59 | 2.0007 | 0.909 | 3.55911 | 1.779
$i=30$ | 6.164 | 0.75 | 2.1292 | 0.968 | 3.50826 | 1.648
$i=90$ | 9.827 | 1.20 | 2.3926 | 1.088 | 3.40870 | 1.425

Table 1 presents examples of a few values of the effective and color
temperatures for models with two values of the dimensionless angular
velocities $\bar{\Omega}^{2}=0.30$ and $0.60$ and a few inclination angles
$i$. As to be expected, stars seen in the equatorial plane ($i=0^{\circ}$)
have the lowest effective temperatures ($T_{\rm eff}=2.117\times 10^{7}\,$K
and $T_{\rm eff}=2.001\times 10^{7}\,$K for $\bar{\Omega}^{2}=0.30$ and
$0.60$, respectively), the lowest color temperatures, and the lowest
bolometric fluxes. On the contrary, stars seen pole on are brighter and
apparently hotter. The last column of the Table 1 contains spectral hardening
factor $f_{c}=T_{c}/T_{\rm eff}$ (here, $T_{c}$ is the color temperature
whereas $T_{\rm eff}$ is the effective temperature). The color temperature is
calculated from the Wien displacement law $h\nu_{\rm max}=2.82kT$. The value
of the hardening factor depends on the inclination angle as well as on the
dimensionless angular velocity. The largest values of $f_{\rm c}$ were
obtained for the star seen pole on whereas the lowest ones were for the star
seen in the equatorial plane. For the spherical star with reference parameters
$T_{\rm eff}=2.20\times 10^{7}\,$K and $\log(g)=14.40\,$(cgs), the spectral
hardening factor is equal to $f_{\rm c}=1.50$ and differs from values for the
flattened star given in the last column of the table.

Figure 9: Spectral hardening factor $f_{c}$ for a flattened star with our
reference parameters seen at different inclination angles $i$ and two sets of
parameters $\bar{\Omega}^{2}=0.30$ and $e=0.54$ (blue line labeled 0.30) and
$\bar{\Omega}^{2}=0.60$ and $e=0.69$ (red line labeled 0.60).

Figure 9 presents the dependence of the spectral hardening factor for a
rotationally flattened neutron star with two values of the dimensionless
angular velocity $\bar{\Omega}^{2}=0.30$ and $0.60$ on inclination angles from
$i=0^{\circ}$ (star seen in the equatorial plane) up to $i=90^{\circ}$ (seen
pole on). For stars with higher values of the angular velocity of the rotation
and therefore an eccentricity, the spectral hardening factor has larger values
and changes more rapidly with the inclination angle than in the case of the
stars with a lower value of $\bar{\Omega}^{2}$.

##  4 Summary and discussion

In our paper, we presented spectra of a rotationally distorted neutron star
considering also the case of fast rotation. We have shown the substantial
influence of both the dimensionless angular velocity $\bar{\Omega}^{2}$ and
the inclination angle $i$ on the emitted spectra.

Using the ATM24 code we calculated a small sample grid of theoretical spectra
integrated over the distorted surface of a sample rotating neutron star seen
by a distant observer at various inclination angles. We assumed the following
parameters: two dimensionless angular velocities $\bar{\Omega}^{2}=0.30$ and
$0.60$, the effective temperature of the nonrotating star $T_{\rm
eff}=2.20\times 10^{7}\,$K, the logarithm of the surface gravity of the
spherical star $\log(g)=14.40\,$(cgs) and inclination angles from
$i=0^{\circ}$ to $i=90^{\circ}$ with the step $\Delta i=10^{\circ}$. We assume
that the atmosphere is a mixture of hydrogen and helium with $M_{\rm H}=0.70$
and $M_{\rm He}=0.30$.

The effective gravity and the effective temperature (determined from von
Zeipel law) are strong functions of the dimensionless angular velocity
$\bar{\Omega}^{2}$ (see Eq.2). Therefore, $\bar{\Omega}^{2}$ belongs to
parameters that play a crucial role in the formation of the emergent spectrum.
The eccentricity of the oblate star depends on $\bar{\Omega}^{2}$ and was
obtained by the approximate (independent of the equation of state) formula
from Silva et al. (2021).

Suleimanov et al. (2020) have shown the spectra of the rotating neutron star
as a function of the inclination angle and the angular velocity of the star.
They applied a diluted blackbody spectrum to describe the local emission
corrected by the hardening factor. The authors assumed slow-rotation
approximation during the calculation of the effective gravity. Their Constant
Relative Flux model is the same as our model of a flattened neutron star with
the effective temperature proportional to the local effective gravity as
$T_{\rm eff}(\theta)\sim g(\theta)^{0.25}$. Therefore, our model corresponds
to the upper panel of the Figure 6 in their paper. The lower panel of Fig.6 in
their paper simply ignores the von Zeipel law. The paper mentioned above was
based on the code defined in Suleimanov et al. (2012), where the equation of
transfer has part of the stimulated Compton scattering terms artificially
transferred to the denominator (see Eq. 12 their paper). As a result
expression for the source function was depleted by one of the Compton-induced
scattering terms (see Eq. 14). The impact of such an arbitrary change in their
paper on the atmosphere models near the Eddington limit is unknown. In the ATM
24 code, such an artificial manipulation does not exists.

In our model we applied fast rotation approximation; therefore, we included
the quadrupole moment of the mass distribution in the neutron star interior.
We calculated very carefully specific intensities $I_{\nu}(T_{\rm
eff}(\theta),g_{\rm eff}(\theta))$ at each patch on the surface. These
specific intensities are next integrated over the surface of the neutron star
to obtain spectra that could be seen by the distant observer (Różańska et al.,
2017, 2018). Spectra presented in this paper are not corrected for the
relativistic effects, which will be the subject of our future work. In our
paper, we did not include effects of the special relativity (relativistic
Doppler broadening), which change the shape of the spectrum. This has been
investigated by Bauböck et al. (2015) for a perfect black body spectrum.
Authors concluded that changes due to Doppler broadening are small compared to
the gravitational redshift.

In this paper, we clearly showed that the gravitational darkening effect
should be included in realistic calculations of the rotating neutron star
spectra. This effect modifies the shape of the emerging spectrum and therefore
affects the determination of neutron star parameters obtained by spectral
fitting of real data.

This work was supported by grant No. 2015/17/B/ST9/03422 from the Polish
National Science Center.

## References

  * AlGendy & Morsink (2014) AlGendy, M. & Morsink, S. M. 2014, ApJ, 791, 78
  * Bauböck et al. (2015) Bauböck, M., Özel, F., Psaltis, D., & Morsink, S. M. 2015, ApJ, 799, 22
  * Belvedere et al. (2014) Belvedere, R., Boshkayev, K., Rueda, J. A., & Ruffini, R. 2014, Nuclear Physics A, 921, 33
  * Bhattacharyya (2010) Bhattacharyya, S. 2010, AdSpR, 45, 949
  * Claret (2015) Claret, A. 2015, A&A, 577, A87
  * Claret (2021) Claret, A. 2021, A&A, 648, A111
  * Djurašević et al. (2003) Djurašević, G., Rovithis-Livaniou, H., Rovithis, P., et al. 2003, A&A, 402, 667
  * Djurašević et al. (2006) Djurašević, G., Rovithis-Livaniou, H., Rovithis, P., et al. 2006, A&A, 445, 291
  * Domiciano de Souza et al. (2014) Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569, A10
  * Galloway et al. (2008) Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, ApJS, 179, 360
  * Galloway et al. (2006) Galloway, D. K., Psaltis, D., Muno, M. P., & Chakrabarty, D. 2006, ApJ, 639, 1033
  * Haensel et al. (2007) Haensel, P., Potekhin, A. Y., & Yakovlev, D. G. 2007, Neutron Stars I. Equation of State and Structure, Springer
  * Hessels et al. (2006) Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901
  * Jaroszynski (1986) Jaroszynski, M. 1986, Acta Astron., 36, 355
  * Lucy (1967) Lucy, L. B. 1967, Zeitschrift für Astrophysik, 65, 89
  * Madej (1991) Madej, J. 1991, ApJ, 376, 161
  * Madej et al. (2004) Madej, J., Joss, P. C., & Różańska, A. 2004, ApJ, 602, 904
  * Madej et al. (2017) Madej, J., Różańska, A., Majczyna, A., & Należyty, M. 2017, MNRAS, 469, 2032
  * Majczyna & Madej (2005) Majczyna, A. & Madej, J. 2005, Acta Astron., 55, 349
  * Majczyna et al. (2005) Majczyna, A., Madej, J., Joss, P. C., & Różańska, A. 2005, A&A, 430, 643
  * Muno et al. (2002) Muno, M. P., Chakrabarty, D., Galloway, D. K., & Psaltis, D. 2002, ApJ, 580, 1048
  * Nagirner & Poutanen (1993) Nagirner, D. I. & Poutanen, Y. J. 1993, Astronomy Letters, 19, 262
  * Nagirner & Poutanen (1994) Nagirner, D. I. & Poutanen, Y. J. 1994, Astrophysics and Space Physics Reviews, 9, 1
  * Pomraning (1973) Pomraning, G. C. 1973, The equations of radiation hydrodynamics
  * Różańska et al. (2017) Różańska, A., Bełdycki, B., Madej, J., Adhikari, T. P., & You, B. 2017, Acta Astron., 67, 51
  * Różańska et al. (2018) Różańska, A., Bresler, K., Bełdycki, B., Madej, J., & Adhikari, T. P. 2018, A&A, 612, L12
  * Sampson (1959) Sampson, D. H. 1959, ApJ, 129, 734
  * Silva et al. (2021) Silva, H. O., Pappas, G., Yunes, N., & Yagi, K. 2021, Phys. Rev. D, 103, 063038
  * Suleimanov et al. (2012) Suleimanov, V., Poutanen, J., & Werner, K. 2012, A&A, 545, 120
  * Suleimanov et al. (2020) Suleimanov, V. F., Poutanen, J., & Werner, K. 2020, A&A, 639, A33
  * Vincent et al. (2018) Vincent, F. H., Bejger, M., Różańska, A., et al. 2018, ApJ, 855, 116
  * von Zeipel (1924) von Zeipel, H. 1924, MNRAS, 84, 665
  * Webbink (1976) Webbink, R. F. 1976, ApJ, 209, 829


NEMO-3 Collaboration

# Search for Periodic Modulations of the Rate of Double-Beta Decay of 100Mo in
the NEMO-3 Detector

R. Arnold  IPHC, ULP, CNRS/IN2P3, F-67037 Strasbourg, France  C. Augier  LAL,
Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay,
France  A.S. Barabash  NRC “Kurchatov Institute”, ITEP, 117218 Moscow, Russia
A. Basharina-Freshville  UCL, London WC1E 6BT, United Kingdom  S. Blondel
LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay,
France  S. Blot  University of Manchester, Manchester M13 9PL, United Kingdom
M. Bongrand  LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay,
F-91405 Orsay, France  D. Boursette  LAL, Université Paris-Sud, CNRS/IN2P3,
Université Paris-Saclay, F-91405 Orsay, France  R. Breier  FMFI, Comenius
University, SK-842 48 Bratislava, Slovakia  V. Brudanin  JINR, 141980 Dubna,
Russia  National Research Nuclear University MEPhI, 115409 Moscow, Russia  J.
Busto  Aix Marseille Université, CNRS, CPPM, F-13288 Marseille, France  A.J.
Caffrey  Idaho National Laboratory, Idaho Falls, ID 83415, U.S.A.  S. Calvez
LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay,
France  C. Cerna  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170
Gradignan, France  J.P. Cesar  University of Texas at Austin, Austin, TX
78712, U.S.A.  M. Ceschia  UCL, London WC1E 6BT, United Kingdom  A. Chapon
LPC Caen, ENSICAEN, Université de Caen, CNRS/IN2P3, F-14050 Caen, France  E.
Chauveau  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170 Gradignan,
France  A. Chopra  UCL, London WC1E 6BT, United Kingdom  L. Dawson  UCL,
London WC1E 6BT, United Kingdom  D. Duchesneau  LAPP, Université de Savoie,
CNRS/IN2P3, F-74941 Annecy-le-Vieux, France  D. Durand  LPC Caen, ENSICAEN,
Université de Caen, CNRS/IN2P3, F-14050 Caen, France  G. Eurin  LAL,
Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay,
France  UCL, London WC1E 6BT, United Kingdom  J.J. Evans  University of
Manchester, Manchester M13 9PL, United Kingdom  L. Fajt  Institute of
Experimental and Applied Physics, Czech Technical University in Prague,
CZ-11000 Prague, Czech Republic  D. Filosofov  JINR, 141980 Dubna, Russia  R.
Flack  UCL, London WC1E 6BT, United Kingdom  P. Franchini  Imperial College
London, London SW7 2AZ, United Kingdom  X. Garrido  LAL, Université Paris-Sud,
CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay, France  C. Girard-Carillo
LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay,
France  H. Gómez  LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-
Saclay, F-91405 Orsay, France  B. Guillon  LPC Caen, ENSICAEN, Université de
Caen, CNRS/IN2P3, F-14050 Caen, France  P. Guzowski  University of Manchester,
Manchester M13 9PL, United Kingdom  R. Hodák  Institute of Experimental and
Applied Physics, Czech Technical University in Prague, CZ-11000 Prague, Czech
Republic  A. Huber  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170
Gradignan, France  P. Hubert  Université de Bordeaux, CNRS, CENBG, UMR 5797,
F-33170 Gradignan, France  C. Hugon  Université de Bordeaux, CNRS, CENBG, UMR
5797, F-33170 Gradignan, France  M. H. Hussain  UCL, London WC1E 6BT, United
Kingdom  S. Jullian  LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-
Saclay, F-91405 Orsay, France  A. Klimenko  JINR, 141980 Dubna, Russia  O.
Kochetov  JINR, 141980 Dubna, Russia  S.I. Konovalov  NRC “Kurchatov
Institute”, ITEP, 117218 Moscow, Russia  V. Kovalenko  JINR, 141980 Dubna,
Russia  D. Lalanne  LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-
Saclay, F-91405 Orsay, France  K. Lang  University of Texas at Austin, Austin,
TX 78712, U.S.A.  Y. Lemière  LPC Caen, ENSICAEN, Université de Caen,
CNRS/IN2P3, F-14050 Caen, France  T. Le Noblet  LAPP, Université de Savoie,
CNRS/IN2P3, F-74941 Annecy-le-Vieux, France  Z. Liptak  University of Texas at
Austin, Austin, TX 78712, U.S.A.  X. R. Liu  UCL, London WC1E 6BT, United
Kingdom  P. Loaiza  LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-
Saclay, F-91405 Orsay, France  G. Lutter  Université de Bordeaux, CNRS, CENBG,
UMR 5797, F-33170 Gradignan, France  M. Macko  Institute of Experimental and
Applied Physics, Czech Technical University in Prague, CZ-11000 Prague, Czech
Republic  C. Macolino  LAL, Université Paris-Sud, CNRS/IN2P3, Université
Paris-Saclay, F-91405 Orsay, France  F. Mamedov  Institute of Experimental and
Applied Physics, Czech Technical University in Prague, CZ-11000 Prague, Czech
Republic  C. Marquet  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170
Gradignan, France  F. Mauger  LPC Caen, ENSICAEN, Université de Caen,
CNRS/IN2P3, F-14050 Caen, France  A. Minotti  LAPP, Université de Savoie,
CNRS/IN2P3, F-74941 Annecy-le-Vieux, France  B. Morgan  University of Warwick,
Coventry CV4 7AL, United Kingdom  J. Mott  UCL, London WC1E 6BT, United
Kingdom  I. Nemchenok  JINR, 141980 Dubna, Russia  M. Nomachi  Osaka
University, 1-1 Machikaneyama Toyonaka, Osaka 560-0043, Japan  F. Nova
University of Texas at Austin, Austin, TX 78712, U.S.A.  F. Nowacki  IPHC,
ULP, CNRS/IN2P3, F-67037 Strasbourg, France  H. Ohsumi  Saga University, Saga
840-8502, Japan  G. Oliviéro  LPC Caen, ENSICAEN, Université de Caen,
CNRS/IN2P3, F-14050 Caen, France  R.B. Pahlka  University of Texas at Austin,
Austin, TX 78712, U.S.A.  V. Palusova  Université de Bordeaux, CNRS, CENBG,
UMR 5797, F-33170 Gradignan, France  FMFI, Comenius University, SK-842 48
Bratislava, Slovakia  C. Patrick  UCL, London WC1E 6BT, United Kingdom  F.
Perrot  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170 Gradignan,
France  A. Pin  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170
Gradignan, France  F. Piquemal  Université de Bordeaux, CNRS, CENBG, UMR 5797,
F-33170 Gradignan, France  Laboratoire Souterrain de Modane, F-73500 Modane,
France  P. Povinec  FMFI, Comenius University, SK-842 48 Bratislava, Slovakia
P. Přidal  Institute of Experimental and Applied Physics, Czech Technical
University in Prague, CZ-11000 Prague, Czech Republic  W. S. Quinn  UCL,
London WC1E 6BT, United Kingdom  Y.A. Ramachers  University of Warwick,
Coventry CV4 7AL, United Kingdom  A. Remoto  LAPP, Université de Savoie,
CNRS/IN2P3, F-74941 Annecy-le-Vieux, France  J.L. Reyss  LSCE, CNRS, F-91190
Gif-sur-Yvette, France  C.L. Riddle  Idaho National Laboratory, Idaho Falls,
ID 83415, U.S.A.  E. Rukhadze  Institute of Experimental and Applied Physics,
Czech Technical University in Prague, CZ-11000 Prague, Czech Republic  R.
Saakyan  UCL, London WC1E 6BT, United Kingdom  A. Salamatin  JINR, 141980
Dubna, Russia  R. Salazar  University of Texas at Austin, Austin, TX 78712,
U.S.A.  X. Sarazin  LAL, Université Paris-Sud, CNRS/IN2P3, Université Paris-
Saclay, F-91405 Orsay, France  J. Sedgbeer  Imperial College London, London
SW7 2AZ, United Kingdom  Yu. Shitov  JINR, 141980 Dubna, Russia  Imperial
College London, London SW7 2AZ, United Kingdom  L. Simard  LAL, Université
Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay, France
Institut Universitaire de France, F-75005 Paris, France  F. Šimkovic  FMFI,
Comenius University, SK-842 48 Bratislava, Slovakia  A. Smetana  Institute of
Experimental and Applied Physics, Czech Technical University in Prague,
CZ-11000 Prague, Czech Republic  A. Smolnikov  JINR, 141980 Dubna, Russia  S.
Söldner-Rembold  University of Manchester, Manchester M13 9PL, United Kingdom
B. Soulé  Université de Bordeaux, CNRS, CENBG, UMR 5797, F-33170 Gradignan,
France  I. Štekl  Institute of Experimental and Applied Physics, Czech
Technical University in Prague, CZ-11000 Prague, Czech Republic  J. Suhonen
Jyväskylä University, FIN-40351 Jyväskylä, Finland  C.S. Sutton  MHC, South
Hadley, Massachusetts 01075, U.S.A.  G. Szklarz  LAL, Université Paris-Sud,
CNRS/IN2P3, Université Paris-Saclay, F-91405 Orsay, France  H. Tedjditi  Aix
Marseille Université, CNRS, CPPM, F-13288 Marseille, France  J. Thomas  UCL,
London WC1E 6BT, United Kingdom  V. Timkin  JINR, 141980 Dubna, Russia  S.
Torre  UCL, London WC1E 6BT, United Kingdom  Vl.I. Tretyak  Institute for
Nuclear Research, 03028, Kyiv, Ukraine  V.I. Tretyak  JINR, 141980 Dubna,
Russia  V.I. Umatov  NRC “Kurchatov Institute”, ITEP, 117218 Moscow, Russia
I. Vanushin  NRC “Kurchatov Institute”, ITEP, 117218 Moscow, Russia  C. Vilela
UCL, London WC1E 6BT, United Kingdom  V. Vorobel  Charles University, Faculty
of Mathematics and Physics, CZ-12116 Prague, Czech Republic  D. Waters  UCL,
London WC1E 6BT, United Kingdom  F. Xie  UCL, London WC1E 6BT, United Kingdom

###### Abstract

Double-beta decays of 100Mo from the 6.0195-year exposure of a 6.914 kg high-
purity sample were recorded by the NEMO-3 experiment that searched for
neutrinoless double-beta decays. These ultra-rare transitions to 100Ru have a
half-life of approximately $7\times 10^{18}$ years, and have been used to
conduct the first ever search for periodic variations of this decay mode. The
Lomb-Scargle periodogram technique, and its error-weighted extension, were
employed to look for periodic modulations of the half-life. Monte Carlo
modeling was used to study the modulation sensitivity of the data over a broad
range of amplitudes and frequencies. Data show no evidence of modulations with
amplitude greater than 2.5% in the frequency range of $0.33225\,{\rm y^{-1}}$
to $365.25\,{\rm y^{-1}}$.

###### pacs:

23.40.-s; 14.60.Pq

##  I Introduction

The invariance of fundamental constants of nature has been scrutinized in a
broad range of physics contexts including considerations discussed by Milne
Milne , Walker Walker , Dirac Dirac-Nature ; Dirac:1938mt , Chandrasekhar
Chadrasekhar , and Kothari Kothari  during the initial ascent of data-based
cosmology of the expanding Universe. Modern cosmology and the observed
evolution of the Universe are closely related to properties of interactions of
elementary particles, and impose tight bounds on possible changes in constants
that include the strength of gauge couplings as a function of time elapsed
since the Big Bang. A number of authors studied these constraints and their
implications, e.g., Weinberg-1983a ; Weinberg:1983xy ; Kolb:1985sj ;
Preskill:1988na ; Olive:2002tz ; Olive:2011xvo ; Uzan:2015uba ;
Fritzsch:2016ewd ; Braconi:2018gxo ; Balcerzak:2019smn . Many theoretical
ideas and implied phenomena are discussed in literature and have been reported
in reviews Uzan:2010pm ; Martins:2017yxk  that also contain exhaustive lists
of references on this subject.

Closely related to the time invariance of constants is their periodicity,
given the ubiquitous presence of cyclic processes in nature at almost all
distance and time scales. There are ongoing searches and tests of such
phenomena, mostly but not only connected to dark matter and dark energy
Lewin:1995rx ; Freese:2012xd ; Mayet:2016zxu . Some experiments have produced
controversial results suggesting yearly modulation of WIMP interactions, e.g.,
Bernabei:2010mq ; Bernabei:2008yi ; Bernabei:2018yyw ; Aalseth:2010vx ;
Aalseth:2011wp  as an expected ‘smoking gun’ observable in direct searches of
dark matter Drukier:1986tm . However, these claims of low-mass WIMP dark
matter are not supported by other experiments Aprile:2018dbl ; Aprile:2019jmx
; Akerib:2016vxi ; Agnese:2018gze ; Agnes:2018oej . Also controversial are
results of a posteriori data analyses of measurements of nuclear decay half-
lives which have yielded unexpected periodicities, including annual
modulations that authors have linked to the periodicity of the Earth–Sun
distance and solar activity Jenkins:2008tt ; Jenkins:2008vn ; Fischbach:2009zz
; Javorsek:2010sr ; Sturrock:2010 ; Sturrock:2010bu ; Sturrock:2014caa ;
Sturrock:2012gs ; Jenkins:2012jc ; Sturrock:2012re ; Sturrock:2019dfx ;
parkhomov . However, other analyses of the same data do not reveal any
significant modulations Pomme:2017dfi ; Pomme:2019mhq ; Dhaygude:2019hka .

Difficulties with precision testing the time variation of half-lives of long-
lived radioisotopes can be attributed, in part, to two important factors: the
time duration of measurements, and the necessity of long-term control of
background phenomena, which often exhibit time and/or seasonal periodicities.
Examples of such phenomena include the average environmental temperature;
radon levels in the ground, buildings, and caverns; the seasonal cosmic-ray
flux modulation; the solar wind intensity related to the Earth–Sun distance;
the phase of the lunar cycle and tides; solar activity, etc.

In this Letter, we are reporting results of a search for periodicity of
double-beta ($\rm\beta\beta$) decays of 100Mo $\rightarrow$100Ru. Data were
collected over a period of approximately eight years by the NEMO-3 experiment
NEMO-3:2019gwo ; Arnold:2015wpy . The experiment, designed to search for
neutrinoless double beta decay ($\rm 0\nu\beta\beta$), had an exquisite
capability of background identification and suppression. NEMO-3 collected an
unprecedentedly large data set of 2-neutrino double-beta decays ($\rm
2\nu\beta\beta$). These extremely rare events with two electrons in the final
state of the 100Mo decay, whose half-life is about $7\times 10^{18}$ years,
served as a unique testing ground for the first ever search for periodicity of
a second-order weak transition on time scales shorter than or comparable to
the measuring period.

##  II The NEMO-3 Experiment

The NEMO-3 detector Arnold:2004xq  was designed to detect two electrons in the
final state of neutrinoless double-beta decays. Thin foils of the source
isotopes were surrounded by a tracking chamber and plastic calorimeter blocks
that reconstructed the full kinematics of various decays and interactions
within the detector. The source foils were strips about 65 ${\rm mm}$ wide,
2,480 ${\rm mm}$ long, and $40-60$ ${\rm mg/cm{{}^{2}}}$ thick and were made
of various double-beta decay isotopes (100Mo, 82Se, 130Te, 116Cd, 150Nd, 96Zr,
and 48Ca) totaling about 10 kg.

These foil strips were arranged vertically to form a cylinder such that the
tracking and calorimetric volumes on either side formed a toroidal geometry
for the whole detector. For the purposes of this analysis, only events from
the 100Mo foils were considered. This particular isotope constituted the
majority, 6.914 kg, of the total source mass in NEMO-3. The tracking volume
that surrounded either side of the foils was comprised of 6,180 drift wire
cells, operating in Geiger mode, within a gas mixture of helium-argon
(95%-1%), and ethanol with water vapor (4%). The tracking volume was then
further enclosed by the calorimeter walls composed of 1,940 plastic
scintillator blocks coupled to photomultiplier tubes. Finally, a large
solenoid encircled the detector to produce a 25 ${\rm G}$ magnetic field to
help with e+/e- discrimination. The entire detector was shielded from external
backgrounds by a combination of iron, wood, paraffin, and borated water. The
detector was placed in the Modane Underground Laboratory in the Fréjus tunnel
in the Alps which provided 4800 m.w.e. overburden to shield from cosmic rays.

The detector was operated from early 2003 until early 2011 with data taking
split into two run periods known as Phase 1 (February 2003 - September 2004)
and Phase 2 (October 2004 - January 2011). During Phase 2, an additional
enclosure was installed, surrounding the detector. Filled with radon-filtered
air, this enclosure greatly reduced radon permeation into the detector. The
result was a significant increase in the purity of the signal channel for
Phase 2. Only these lower background runs were considered for the final
analysis presented herein. The total span of the runs from this period amounts
to $6.0195$ ${\rm y}$.

##  III Event Selection

The purpose of this analysis was to search for periodic trends in the rate of
double-beta decays (with no distinction between $\rm 2\nu\beta\beta$ and
theoretical $\rm 0\nu\beta\beta$ decays) originating in the 100Mo source
foils, which required counting the number of such decays per unit time
recorded in the NEMO-3 detector. Data taking was divided into specific run
periods with durations ranging from tens of minutes to just over two days. The
100Mo activity yielding double-beta events was about $0.1$ ${\rm Bq}$ and,
since most runs lasted longer than 20 minutes, most of them accumulated over
100 double-beta events. A discrete time series of the decay rate was
constructed by taking the number of selected events in a run divided by the
duration of the run to yield a value for the observed rate. The average
double-beta event rate for each run was associated with a timestamp,
corresponding to the midpoint of the run.

Events that were selected for inclusion into this calculation were chosen
based on a wide range of criteria to minimize the contribution from background
processes which could mimic topologies of a double-beta decay. The primary
characteristics of a double-beta decay event are the identification of two
tracks with curvatures consistent with negatively-charged particles,
originating from the 100Mo foil, and with associated energy deposits in
scintillator blocks; with no coincident alpha particles (short straight
tracks); and no gamma particles (unassociated scintillator hits) with energy
deposits greater than or equal to $150$ ${\rm keV}$.

The electron track lengths were required to sum to at least $60$ ${\rm cm}$
and trace back to a common vertex in one of the 100Mo source foils (with no
more than $4$ ${\rm cm}$ of difference between the two tracks either
transverse to, or in the plane of the foils). Vertices were required to fall
outside of “hot spots” – regions that were identified as contaminated with
radio-impurities. The electrons were required to fire Geiger cells within $50$
${\rm cm}$ of the vertex on the foil. Electron energy deposits were required
to exceed $200$ ${\rm keV}$ in each of two separate scintillator blocks, each
of which had no neighboring blocks with energy deposits. The extrapolated
electron tracks were required to hit the front face of each block, and each
block had to be coupled to a PMT that had not been flagged with potential
issues such as excess noise or a lack of calibration data. Electrons that hit
the blocks nearest to the foil on the endcaps (scintillator blocks at the top
and bottom of the detector) were also rejected due to a higher chance of
incorrectly identifying the track curvature and thus particle charge.

Time-of-flight measurements obtained from calorimeter hit times were also used
to calculate the probability of the event originating within the foil,
compared to being an external event in which a particle passes through the
foil. The event had to satisfy the condition that the probability of
originating within the foil was greater than or equal to 4% while the
probability of being an external crossing event was less than or equal to 1%.
Only runs that were deemed to be of good quality were used. An event display
for a candidate double-beta decay which passed all such selection criteria is
shown in Figure 1. This event selection resulted in a very high sample purity
(a very similar event selection resulted in a signal to background ratio of 76
Arnold:2013dha ) and no background subtraction was performed in calculating
per-run event rates. The possibility of time varying backgrounds was also
considered, in particular due to the potential seasonal variation of radon
levels. No modulations appeared in the analysis, of either the Phase 1 or
Phase 2 data, at frequencies corresponding to such processes.

Figure 1: An example event display from NEMO-3 showing hit Geiger cells (blue
circles) and calorimeter modules (red boxes). A helix has been fit to the hit
cells to show the curvature of the two electron tracks (due to the magnetic
field) and to show that they share a common origin on the source foil.

The decay rate based on the raw number of events that passed these selection
requirements had to be corrected by the efficiency of the detector during each
run period. This correction scaling could vary from one run to the next and
was accounted for using precise and comprehensive Monte Carlo simulations of
the NEMO-3 experiment.

The final corrected 100Mo double-beta decay rate time series as measured by
the NEMO-3 detector is shown in the upper plot of Figure 2. The data yielded
an average per-run rate of $0.085$ ${\rm Hz}$ with a standard deviation of
11.8%. The average run was 8.3 hours in duration (standard deviation of 4.7
hours) and saw approximately 115.9 events in that time (standard deviation of
66.7 events). In total, 449,733 events were collected, with a mean efficiency
of 4.6%, across the 3,869 run periods. A study of the most significant
systematic errors showed that their contributions were small compared to
statistical fluctuations from one run to the next.

##  IV The Search for Periodicities

A common approach to searching for periodic trends in data involves
decomposing the time series of interest into its spectral frequency
components. In one such approach, a periodogram (also called a power spectrum)
can be constructed by calculating, for any number of desired sample
frequencies, a quantity known as the power. Frequencies which produce a large
power are those which have a stronger presence in the underlying data set. The
basic periodogram technique developed by Lomb lomb  was eventually modified to
allow for unevenly sampled data to be analyzed and became known as the Lomb-
Scargle (LS) periodogram scargle . Further developments eventually allowed the
technique to also account for weighted data and an overall offset term, known
as the Generalized Lomb-Scargle (GLS) periodogram zechmeister .

Periodogram analyses, which are commonplace in astronomy and astrophysics,
have also seen use in the fields of nuclear and particle physics. Some of the
previously mentioned dark matter searches, e.g., Bernabei:2010mq ;
Bernabei:2008yi ; Bernabei:2018yyw  and radioactive decay analyses, e.g.,
Javorsek:2010sr ; Sturrock:2010 ; Sturrock:2010bu ; Sturrock:2014caa  and
parkhomov  have used LS periodogram analyses similar or identical to the
approach used herein. Other examples also include searches for periodic
variations in neutrino fluxes across various experiments, e.g., Milsztajn ;
yoo ; sturrock ; Sturrock:1999 ; aharmim ; Altmann , while uses of the GLS
periodogram appear in other recent analyses Tejas ; Dhaygude:2019hka ; Nakano
.

For a discrete time series $X(t_{j})$ consisting of $N$ entries, the basic
Lomb-Scargle power, $P_{\text{LS}}(\omega)$ can be calculated at a given
sample frequency, $f$ (where $\omega=2\pi f$), by

$\begin{split}P_{\text{LS}}(\omega)=\frac{1}{2\sigma^{2}}&\left\\{\frac{\left(\displaystyle\sum_{j=1}^{N}[X(t_{j})-\bar{X}]\cos[\omega(t_{j}-\tau)]\right)^{2}}{\displaystyle\sum_{j=1}^{N}\cos^{2}[\omega(t_{j}-\tau)]}\right.\\\
&+\left.\frac{\left(\displaystyle\sum_{j=1}^{N}[X(t_{j})-\bar{X}]\sin[\omega(t_{j}-\tau)]\right)^{2}}{\displaystyle\sum_{j=1}^{N}\sin^{2}[\omega(t_{j}-\tau)]}\right\\},\end{split}$
(1)

where $\bar{X}$ is the mean of the data points, $\sigma$ is their standard
deviation, and $\tau$ is defined by the relation
$\tan(2\omega\tau)=\sum_{j=1}^{N}\sin(2\omega
t_{j})/\sum_{j=1}^{N}\cos(2\omega t_{j})$. The periodogram is constructed by
calculating this power over a range of frequencies of interest. While powerful
in its handling of unevenly spaced data, the LS technique doesn’t weight each
data point by its uncertainty.

A well known property of the LS technique is its equivalence to least-squares
fitting of sine waves scargle . By taking into account an offset term and
weights, the GLS technique extends this equivalence to a full ${\chi}^{2}$
fitting approach. The end result is a new expression for the Generalized Lomb-
Scargle power, $P_{\text{GLS}}(\omega)$, given by

$\begin{split}P_{\text{GLS}}(\omega)=\frac{1}{XX\cdot
D}&\left[SS\cdot(XC)^{2}+CC\cdot(XS)^{2}\right.\\\ &\left.-2CS\cdot XC\cdot
XS\right]\end{split}$ (2)

where $D=CC\cdot SS-(CS)^{2}$ and the following abbreviations are used (with
summations running over the same indices as in Eq. (1)):
$XX=\sum w_{j}[X(t_{j})-\bar{X}]^{2},\\\ XC=\sum
w_{j}[X(t_{j})-\bar{X}]\cos(\omega t_{j}),\\\ XS=\sum
w_{j}[X(t_{j})-\bar{X}]\sin(\omega t_{j}),\\\ CC=\sum w_{j}\cos^{2}(\omega
t_{j})-[\sum w_{j}\cos(\omega t_{j})]^{2},\\\ SS=\sum w_{j}\sin^{2}(\omega
t_{j})-[\sum w_{j}\sin(\omega t_{j})]^{2},\\\ CS=\sum w_{j}\cos(\omega
t_{j})\sin(\omega t_{j})-[\sum w_{j}\cos(\omega t_{j})\times\sum
w_{j}\sin(\omega t_{j})]$.
Here the $w_{j}$ are the weights for each $X(t_{j})$, given by
$w_{j}=\frac{1}{W}\frac{1}{\sigma_{j}}$ for $W=\sum\frac{1}{\sigma_{j}}$ where
the $\sigma_{j}$ are the errors, and so the mean is now $\bar{X}=\sum
w_{j}X(t_{j})$.

Another known feature of the LS periodogram is that if the time series
$X(t_{j})$ are made up of Gaussian random values with no underlying modulation
then the resultant periodogram powers should be exponentially distributed with
unit mean scargle . This allows one to estimate the false-alarm probability
(F.A.P.), which gives the probability of finding a power larger than $P$, via
the expression

$\text{F.A.P.}(P)=1-(1-e^{-P})^{M},$ (3)

where $M$ is the number of frequencies sampled. This lets one calculate a
percentage confidence level (C.L.) value for a given power as

$\text{C.L.}(P)=(1-e^{-P})^{M}\times 100\%.$ (4)

This approach to estimating the significance of periodogram peaks is very
straightforward but depends critically on the initial assumption of the
Gaussian distribution of the data points. Furthermore, the GLS periodogram
expression given in Eq. (2) requires a normalization factor, for which the
authors of zechmeister  offer multiple approaches, in order for Eq. (4) to
apply. Here, the normalization scheme proposed by Baluev baluev  was used so
that LS and GLS powers could be shown on the same scale. However, due to these
added complications, the most reliable way to assess the significance of
periodogram peaks is via Monte Carlo methods which are described in section V.
The simple but rough significance approximation obtained from Eq. (4) was thus
employed for sensitivity studies which were optimized to reduce computational
loads.

Both the LS and GLS techniques were employed in this analysis. The two were
used as both a cross-check to one another as well as for consistency when
comparing to other analyses which predominantly relied on the more commonly-
used LS technique. To prevent analysis bias, a blinded approach was taken by
first applying both techniques to time-shuffled Phase 1 data and then to
unshuffled Phase 1 data before finally applying them to the Phase 2 time
series. Here, the term shuffling is used to describe a randomized re-ordering
of data points so as to destroy any potential underlying temporal trends in
the data.

The final time series contained 3,869 runs and an oversampling factor of two
was used. This meant that 7,738 frequencies, twice the number of data points,
were sampled. These frequencies were evenly distributed in the range
$[0.33225,365.25]\,\text{y${}^{-1}$}$ which was chosen based on the maximum
and minimum modulation periods that could be detectable within the given
duration of data taking. The minimum period was limited by the average run
spacing ($\Delta T_{avg}$) such that $f_{\text{max}}=1/(2\times\Delta
T_{avg})$. A conservative value of 12 hours was used for the average run
spacing which translates to a minimum period of one day or
$f_{\text{max}}=365.25\,\text{y${}^{-1}$}$. Similarly, the maximum period was
determined by the total span of the data. In this case, the minimum sample
frequency was chosen such that at least two full periods of a modulation would
be contained in the data which implies that
$f_{\text{min}}=2/(\text{span})=0.33225\,\text{y${}^{-1}$}$. The resultant LS
and GLS periodogram are shown in the lower plot of Figure 2.

Upon constructing the LS and GLS periodograms for the data, the largest LS
power of 8.78 was found at the frequency of $76.26\,\text{y${}^{-1}$}$,
corresponding to a periodicity of approximately 4.8 days, while the largest
GLS power of 6.16 was found at the frequency of $0.47\,\text{y${}^{-1}$}$,
corresponding to a periodicity of approximately 777 days or 2.1 years. The
fact that the locations of these peaks disagreed between the two techniques,
and that neither had any correlation with periodicities found in any of the
previously mentioned references, indicates they are very likely truly random
fluctuations. As further evidence, the GLS power found at the same frequency
as the largest LS power is smaller than any of the top three largest GLS
powers, and vice versa, wich emphasizes how uncorrelated the maximal peaks are
between the two techniques. The sizes of the largest periodogram peaks are
analyzed in more depth in the next section to further corroborate this
assertion.

Figure 2: Results of applying both of the periodogram techniques to the
double-beta decay rate time series of Phase 2 100Mo data. The upper plot shows
the efficiency-corrected time series that was analyzed. The lower plot shows
the periodograms obtained from applying the LS and GLS techniques. The
horizontal lines denote the threshold powers, in each technique, for reaching
a 90% confidence level based on the “shuffle test” method described in Section
V.

##  V Significance and Sensitivity Studies

A common approach to determining the significance of periodogram peaks is to
test how often a specific power is exceeded in a similar time series comprised
of randomized pseudo-data with no modulation. The same shuffling procedure,
described in the previous section to blind the Phase 1 data, was applied to
the Phase 2 data 10,000 times, to create a collection of randomized, null-
hypothesis time series. This preserved the structure of the original data in
terms of the actual values and their associated errors as well as their
temporal spacing. The LS and GLS techniques were then applied to each of these
time series and the largest powers were recorded for each resultant
periodogram. The largest powers from the true data could then be compared
against these maximal powers from the pseudo-data to estimate their
significance.

This so called “shuffle test” method has been used in a similar decay rate
analysis by Sturrock Sturrock:2010  and has also been attributed to an
analysis searching for modulations in the Homestake solar neutrino experiment
bahcall . The results of applying the shuffle test for the LS and GLS
techniques are shown in Figure 3. The division in each distribution shows what
percentage of maximal powers lie above the largest power found in the data.
For the LS (GLS) technique where the largest power was 8.78 (6.16), a larger
power was found in the shuffled pseudo-data sets 54.1% (44.8%) of the time.

Figure 3: A distribution of the maximal LS (top) and GLS (bottom) powers in
10,000 randomly shuffled null-hypothesis data sets (no applied modulations).
The fraction of data sets with maximal power greater than the maximum LS (GLS)
power observed in the data constitute 54.1% (44.8%). This estimates the
significance of the observed maximum LS (GLS) power in the data to be at the
45.9% (55.2%) confidence level.

Further studies with these Monte Carlo data sets were also undertaken to
estimate the sensitivity of the data to detecting different modulations. To do
this, a modulation to the pseudo-data sets was applied in the form of

$R(t_{i})=N[1+A\times\sin(2\pi ft_{i}+\phi)],$ (5)

where $N$ is a normalization constant to match the mean of the un-modulated
data; $A$ represents a fractional or relative amplitude (which will be denoted
as a percentage relative to the mean rate) for the applied modulation; $f$ is
the frequency of the applied modulation; and $\phi$ is the modulation phase.

Although the final analysis only sampled frequencies corresponding to
modulations in which at least two full periods are present in the data, the
sensitivity studies sampled down to even smaller frequencies. This helped
determine an apt lower bound and allowed for the study of the behavior of the
two techniques in the low frequency domain. Due to the computationally
intensive nature of these studies, the sample frequency range was broken down
into three regimes (low, mid, and high) that included different spacing
between frequencies to limit the total number of trials that were needed. A
summary of the various amplitudes and frequencies that were used for the
injected modulations are shown in Table 1. Early studies showed that the
modulation phase only affected the sensitivity contours in the low frequency
regime (which was irrelevant to the final analysis parameters). Thus the phase
was set to zero in generating the results shown below (more will be discussed
about the effects of the phase value further on).

| Amplitude ($A$) [%] | Step Size ($\Delta A$) [%]
---|---|---
| 0.5 - 4 | 0.1
Range | Freq. ($f$) [y-1] | Step Size ($\Delta f$) [y-1]
Low | 0.03 - 0.1 | 0.005
Mid | 0.15 - 2 | 0.05
High | 15 - 360 | 15
Table 1: Summary of the amplitude and frequency ranges and spacings used for
the injected modulation signals to test for detection sensitivity.

For each point in this amplitude-frequency phase space 100 different Monte
Carlo pseudo-data sets were analyzed to average out random variations. At each
combination of modulation amplitude and frequency, the LS and GLS periodograms
were constructed for the different pseudo-data sets and the average (across
all 100 sets) of the largest power was recorded, as well as its estimated C.L.
value derived from Eq. (4). This was used to create a contour plot of C.L.
values in the modulation amplitude-frequency space. These are shown, for the
three frequency regimes, in Figure 4 for both techniques.

Figure 4: An illustration of how the LS and GLS detection sensitivity changes
as the input modulation parameters $A$ and $f$ are varied (for $\phi=0$ and
three different frequency regimes - low, mid, and high). The z-axis, color, at
each point denotes the estimated significance (in terms of C.L. value) of the
largest peak, averaged over 100 periodograms.

For mid and high ranges of modulation frequencies, the detection sensitivity
only depended on the modulation amplitude. In these regimes the threshold for
detection at approximately a 95% C.L. was about 2.0% (2.5%) relative amplitude
for the LS (GLS) technique. The difference between the two values is due to
the GLS periodogram effectively using fewer data points (which correlates with
less sensitivity) by down-weighting some of the runs with larger errors on the
measured rate. Once the frequency dropped low enough, wherein the sample
period exceeded the time span of the data, this detection threshold amplitude
began to vary with frequency (as seen in the left-most plots of Figure 4). In
these regimes, only partial modulation waveforms are being captured in the
time span of the data and the modulation phase began to have an affect on the
results. If the phase was such that a mostly-linear portion of the modulation
was captured (well between a peak and a trough in the underlying oscillation)
then the sensitivity was reduced compared to other phase values.

##  VI Summary and Conclusions

The large, highly-pure sample of double-beta decays observed by the NEMO-3
collaboration provided a unique opportunity to probe the variability of a
second-order, weak nuclear process. A power spectrum analysis was used to
search for periodicities in the double-beta decay rate of 100Mo based on
$6.0195\,{\rm y}$ of data. Periodograms were generated using both the Lomb-
Scargle technique (to be consistent with various previous searches for
periodically varying decay rates), as well as its error-weighted extension,
the Generalized Lomb-Scargle technique, with both resulting in power spectra
consistent with the null hypothesis of no underlying modulation.

This conclusion was reached by noting that the largest LS (GLS) power in the
data-generated periodogram was exceeded 54.1% (44.8%) of the time by the
largest LS (GLS) power from completely randomized, time-shuffled data sets.
Furthermore, the frequencies at which these powers were found did not
correlate with each other nor with any previously claimed periodicities at, or
around, $1\,{\rm y^{-1}}$ for those relating to Earth’s orbital period or
those in the range $(10-15)\,{\rm y^{-1}}$ for those relating to solar synodic
rotation rates Jenkins:2008tt ; Jenkins:2008vn ; Fischbach:2009zz ;
Javorsek:2010sr ; Sturrock:2010 ; Sturrock:2010bu ; Sturrock:2014caa ;
Sturrock:2012gs ; Jenkins:2012jc ; Sturrock:2012re ; Sturrock:2019dfx ;
parkhomov . We estimate that the analyzed data set was sensitive to
modulations with periods between one day and three years, if the relative
amplitude of such modulations had exceeded approximately 2.0% (2.5%) when
using the LS (GLS) technique. Although these constraints are an order of
magnitude weaker than modulation searches in single-beta decay, nevertheless
they represent the first ever limits for the second order process of double-
beta decay.

###### Acknowledgements.

We thank the staff of the Modane Underground Laboratory for their technical
assistance in running the experiment. We acknowledge support by the grant
agencies of the Czech Republic, CNRS/IN2P3 in France, RFBR in Russia (NCNIL
No19-52-16002), APVV in Slovakia, the Science and Technology Facilities
Council, part of U.K. Research and Innovation, and the NSF in the U.S.

## References

  * (1) E. A. Milne, Proc. Roy. Soc. London A 158, 324 (1937).
  * (2) A. G. Walker, Proc. Lond. Math. Soc. 42, 90 (1937).
  * (3) P. A. M. Dirac, Nature (London) 139, 323 (1937).
  * (4) P. A. M. Dirac, Proc. Roy. Soc. Lond. A 165, 199 (1938).
  * (5) S. Chandrasekhar, Nature (London) 139, 757 (1937).
  * (6) D. S. Kothari, Nature (London) 142, 354 (1938).
  * (7) S. Weinberg, Phil. Trans. R. Soc. London A 310, 249 (1983).
  * (8) S. Weinberg, Phys. Lett. B 125, 265 (1983).
  * (9) E. W. Kolb, M. J. Perry and T. P. Walker, Phys. Rev. D 33, 869 (1986).
  * (10) J. Preskill, Nucl. Phys. B 323, 141 (1989).
  * (11) K. A. Olive et al., Phys. Rev. D 66, 045022 (2002).
  * (12) K. A. Olive, Can. J. Phys. 89, no. 4, 361 (2011).
  * (13) J. P. Uzan, Comptes Rendus Physique 16, 576 (2015).
  * (14) H. Fritzsch, J. Solà and R. C. Nunes, Eur. Phys. J. C 77, no. 3, 193 (2017).
  * (15) A. Braconi et al., Phys. Rev. D 100, 015032 (2019).
  * (16) A. Balcerzak and K. Marosek, Eur. Phys. J. C 79, no. 7, 563 (2019).
  * (17) J. P. Uzan, Living Rev. Rel. 14, 2 (2011).
  * (18) C. J. A. P. Martins, Rep. Prog. Phys. 80, 126902 (2017).
  * (19) J. D. Lewin and P. F. Smith, Astropart. Phys. 6, 87 (1996).
  * (20) K. Freese, M. Lisanti and C. Savage, Rev. Mod. Phys. 85, 1561 (2013).
  * (21) F. Mayet et al., Phys. Rept. 627, 1 (2016).
  * (22) R. Bernabei et al. [DAMA Collaboration], Eur. Phys. J. C 56, 333 (2008).
  * (23) R. Bernabei et al. [DAMA Collaboration], Eur. Phys. J. C 67, 39 (2010).
  * (24) R. Bernabei et al. [DAMA Collaboration], Nucl. Phys. Atom. Energy 19, no.4, 307-325 (2018).
  * (25) C. E. Aalseth et al. [CoGeNT Collaboration], Phys. Rev. Lett. 106, 131301 (2011).
  * (26) C. E. Aalseth et al., Phys. Rev. Lett. 107, 141301 (2011).
  * (27) A. K. Drukier, K. Freese and D. N. Spergel, Phys. Rev. D 33, 3495 (1986).
  * (28) E. Aprile et al. [XENON Collaboration], Phys. Rev. Lett. 121, 111302 (2018).
  * (29) E. Aprile et al. [XENON Collaboration], Phys. Rev. Lett. 123, 241803 (2019).
  * (30) D. S. Akerib et al. [LUX Collaboration], Phys. Rev. Lett. 118, 021303 (2017).
  * (31) R. Agnese et al. [SuperCDMS Collaboration], Phys. Rev. D 99, 062001 (2019).
  * (32) P. Agnes et al. [DarkSide Collaboration], Phys. Rev. Lett. 121, 111303 (2018).
  * (33) J. H. Jenkins and E. Fischbach, Astropart. Phys. 31, 407 (2009).
  * (34) J. H. Jenkins et al., Astropart. Phys. 32, 42 (2009).
  * (35) E. Fischbach et al., Space Sci. Rev. 145, 285 (2009).
  * (36) D. Javorsek II et al., Astropart. Phys. 34, 173 (2010).
  * (37) P. A. Sturrock et al., Astropart. Phys. 34, 121 (2010).
  * (38) P. A. Sturrock et al., Solar Phys. 267, 251 (2010).
  * (39) P. A. Sturrock, E. Fischbach and J. Jenkins, Astrophys. J. 794, no. 1, 42 (2014) [Erratum: Astrophys. J. 796, no. 2, 149 (2014)].
  * (40) P. A. Sturrock et al., Astropart. Phys. 36, 18 (2012).
  * (41) J. H. Jenkins et al., Astropart. Phys. 37, 81 (2012).
  * (42) P. A. Sturrock et al., Astropart. Phys. 42, 62 (2013).
  * (43) P. A. Sturrock et al., arXiv:1907.11749 [hep-ph].
  * (44) A. G. Parkhomov, J. Mod. Phys. 2(11), 1310 (2011).
  * (45) S. Pomme, K. Kossert and O. Nahle, Solar Phys. 292, no.11, 162 (2017).
  * (46) S. Pomme, Eur. Phys. J. C 79, no.1, 73 (2019).
  * (47) A. Dhaygude and S. Desai, Eur. Phys. J. C 80, no.2, 96 (2020).
  * (48) R. Arnold et al. [NEMO-3 Collaboration], Eur. Phys. J. C 79, no. 5, 440 (2019).
  * (49) R. Arnold et al. [NEMO-3 Collaboration], Phys. Rev. D 92, 072011 (2015).
  * (50) R. Arnold et al. [NEMO-3 Collaboration], Nucl. Instrum. Meth. A 536, 79 (2005).
  * (51) R. Arnold et al. [NEMO-3 Collaboration], Phys. Rev. D 89, no.11, 111101 (2014)
  * (52) N. R. Lomb, Astrophys. Space Sci. 39 447 (1976).
  * (53) J. D. Scargle et al., Phys. Rev. Lett. 77, 5186 (1996).
  * (54) M. Zechmeister and M. Kürster, Astron. Astrophys. 496, 577 (2009).
  * (55) A. Milsztajn, arXiv:hep-ph/0301252v1
  * (56) J. Yoo et al. [Super-Kamiokande Collaboration], Phys. Rev. D 68, 092002 (2003).
  * (57) P. A. Sturrock et al., Phys. Rev. D 72, 113004 (2005).
  * (58) P. A. Sturrock et al., Astrophys. J. 523, 177 (1999).
  * (59) B. Aharmim et al. [SNO Collaboration], Phys. Rev. D 72, 052010 (2005).
  * (60) M. Altmann et al. [GNO Collaboration], Phys. Lett. B 616, 174 (2005).
  * (61) P. Tejas and S. Desai, Eur. Phys. J. C 78, 554 (2018).
  * (62) Y. Nakano et al. [Super-Kamiokande Collaboration], J. Phys. Conf. Ser. 1342, 012037 (2020).
  * (63) R. V. Baluev, M.N.R.A.S. 385, 1279 (2008).
  * (64) J. N. Bahcall and W. H. Press, Astrophys. J. 370, 730 (1991).