# Recent Neutrino Parameters Impact on the Effective Majorana Neutrino Mass in
0$\nu\beta\beta$ Decay

Dongming Mei<EMAIL_ADDRESS>Physics Department, University of South
Dakota, Vermillion, SD, 57069  Kunming Dong  Physics Department, University of
South Dakota, Vermillion, SD, 57069  Austin Warren  Physics Department,
University of South Dakota, Vermillion, SD, 57069  Sanjay Bhattarai  Physics
Department, University of South Dakota, Vermillion, SD, 57069

###### Abstract

We investigate how recent updates to neutrino oscillation parameters and the
sum of neutrino masses influence the sensitivity of neutrinoless double-beta
(0$\nu\beta\beta$) decay experiments. Incorporating the latest cosmological
constraints on the sum of neutrino masses and laboratory measurements on
oscillations, we determine the sum of neutrino masses for both the normal
hierarchy (NH) and the inverted hierarchy (IH). Our analysis reveals a narrow
range for the sum of neutrino masses, approximately 0.06 eV/c2 for NH and
0.102 eV/c2 for IH. Utilizing these constraints, we calculate the effective
Majorana masses for both NH and IH scenarios, establishing the corresponding
allowed regions. Importantly, we find that the minimum neutrino mass is non-
zero, as constrained by the current oscillation parameters. Additionally, we
estimate the half-life of 0$\nu\beta\beta$ decay using these effective
Majorana masses for both NH and IH. Our results suggest that upcoming ton-
scale experiments will comprehensively explore the IH scenario, while 100-ton-
scale experiments will effectively probe the parameter space for the NH
scenario, provided the background index can achieve 1 event/kton-year in the
region of interest.

††preprint: AIP/123-QED

##  I Introduction

The revelation of neutrino oscillations in solar and atmospheric neutrino
experiments SNO ; Kamland ; SuperK  has marked a significant leap in our
comprehension of neutrinos, uncovering their mass and mixing within the lepton
sector. This breakthrough holds profound implications for extending the
standard model of particle physics and pursuing a grand unified theory of
nature. Yet, despite this progress, numerous facets of neutrinos remain
enigmatic, such as their absolute mass scale, mass hierarchy, and the
determination of whether they exhibit a Dirac or Majorana nature Maj .

Various experiments, spanning accelerator-based and reactor-based neutrino
oscillation studies, beta decay investigations, and neutrinoless double-beta
($0\nu\beta\beta$) decay research, contribute to our understanding of the
absolute mass scale and mass hierarchy of neutrinos. However, it’s noteworthy
that $0\nu\beta\beta$ decay experiments uniquely offer a practical means to
explore the Majorana nature of neutrinos. These experiments aim to address
fundamental questions by probing the effective Majorana mass of electron
neutrinos and determining the absolute neutrino mass scale, hierarchy, and
character Ell02 ; Ell04 ; Avi04 ; Bar04 .

In a $0\nu\beta\beta$ decay experiment, the connection between the measured
half-life (T${}_{\beta\beta}^{0\nu}$) and the effective Majorana mass of the
electron neutrino, $m_{\beta\beta}$, is approximated by the following
expression:

$(T_{\beta\beta}^{0\nu})^{-1}=G^{0\nu}(E_{0},Z)|(\frac{m_{\beta\beta}}{m_{e}})^{2}||M_{f}^{0\nu}-(\frac{g_{A}}{g_{V}})^{2}M_{GT}^{0\nu}|^{2}.$
(1)

In Equation 1, $G^{0\nu}(E_{0},Z)$ includes couplings and a phase space
factor, where $g_{A}$ and $g_{V}$ represent the axial vector and vector
coupling constants, and $M_{f}^{0\nu}$ and $M_{GT}^{0\nu}$ denote the Fermi
and Gamow-Teller nuclear matrix elements, respectively. Consequently, the
precision with which $m_{\beta\beta}$ can be determined from the measured
half-life relies on the accuracy of the theoretical nuclear matrix elements.

If the $0\nu\beta\beta$ decay process is mediated by light Majorana neutrinos,
the effective neutrino mass is determined by the following coherent sum:

$|m_{\beta\beta}|=\left|\sum_{i}m_{i}U_{\beta i}^{2}\right|,$ (2)

where $m_{i}$ represents the mass of the $i$th neutrino mass eigenstate,
$U_{\beta i}$ are elements of the leptonic matrix $U$, and the sum extends
over all light neutrino mass eigenstates bpo . In the standard three-flavor
scheme, $U$ is a unitary matrix and can therefore be parameterized in terms of
three flavor mixing angles and three CP-violating phases:

$U=\left(\begin{array}[]{ccc}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta}\\\
-s_{12}c_{23}-c_{12}s_{13}s_{23}e^{i\delta}&c_{12}c_{23}-s_{12}s_{13}s_{23}e^{i\delta}&c_{13}s_{23}\\\
s_{12}s_{23}-c_{12}s_{13}c_{23}e^{i\delta}&-c_{12}s_{23}-s_{12}s_{13}c_{23}e^{i\delta}&c_{13}c_{23}\end{array}\right)$
(3)

where $c_{ij}=cos\theta_{ij}$, $s_{ij}=sin\theta_{ij}$ (for ij = 12; 13; 23),
$\delta$ represents the Dirac charge-parity (CP) phase.

The effective Majorana mass is determined by the absolute value of the element
$m_{\beta\beta}$ within the mass matrix in the charged lepton flavor basis.
This mass matrix characterizes the charged leptons in their diagonal basis.
Consequently, the Majorana neutrino mass matrix in the charged lepton flavor
basis can be represented as:

$|m_{\beta\beta}|=|U\cdot D\cdot U^{T}|.$ (4)

In this context, the diagonal matrix can be denoted as
$D=\text{diag}[m_{1}e^{i\alpha},m_{2},m_{3}e^{i\beta}]$, where both $\alpha$
and $\beta$ represent the Majorana CP-violating phases.

The results derived from Equation 4 can be expressed as:

$|m_{\beta\beta}|=|m_{\beta\beta}^{(1)}e^{i\alpha}+m_{\beta\beta}^{(2)}+m_{\beta\beta}^{(3)}e^{i\beta}|,$
(5)

and

$\begin{array}[]{ll}|m_{\beta\beta}^{(1)}|=m_{1}c_{12}^{2}c_{13}^{2},\\\
|m_{\beta\beta}^{(2)}|=m_{2}s_{12}^{2}c_{13}^{2},\\\
|m_{\beta\beta}^{(3)}|=m_{3}s_{13}^{2}.\end{array}$ (6)

After years of experimental exploration, significant progress has been made in
determining the values of $\theta_{12}$, $\theta_{13}$, and $\theta_{23}$ with
a high degree of accuracy from current neutrino oscillation data. However, the
three phase parameters ($\delta$, $\alpha$, and $\beta$) remain elusive pdg .
Additionally, while the value of $\Delta m_{21}^{2}\equiv m_{2}^{2}-m_{1}^{2}$
and the absolute value of $\Delta m_{31}^{2}\equiv m_{3}^{2}-m_{1}^{2}$ have
been measured, the sign of $\Delta m_{31}^{2}$ and the absolute neutrino mass
scale remain unknown.

Therefore, the magnitude of $|m_{\beta\beta}|$ is subject to three types of
uncertainties, namely the unknown absolute neutrino mass scale ($m_{1}$ or
$m_{3}$), the unknown neutrino mass hierarchy ($\Delta m_{31}^{2}>0$ or
$\Delta m_{31}^{2}<0$), and the unknown Majorana phases $\alpha$ and $\beta$
appearing in $|m_{\beta\beta}|$, even in the absence of new physics
contamination. Given that the magnitude of $|m_{\beta\beta}|$ is closely tied
to the decay half-life of the $0\nu\beta\beta$ decay process, it is of
significant interest to explore how these three uncertainties will influence
its value. Evaluating this impact will provide valuable insights into the
practical feasibility of conducting an experiment.

Until now, extensive phenomenological endeavors have been dedicated to
exploring the parameter space of $|m_{\beta\beta}|$ and evaluating its
sensitivity to potential new physics phenomena ago . Current experiments have
achieved remarkable sensitivities, with half-life measurements reaching levels
of $\sim 10^{26}$ years for isotopes such as 136Xe and 76Ge, consequently
setting upper limits on the effective Majorana neutrino mass within the range
of 36-156 meV KamLandZen ; Gerda .

In this study, we aim to investigate the parameter space of $|m_{\beta\beta}|$
by incorporating the sum of neutrino masses constrained by recent cosmological
data cos  and utilizing the best-fitted oscillation parameters derived from
global particle physics datasets pdg . We will explore several intriguing
scenarios, including the allowed regions for the sum of neutrino masses and
the minimum neutrino mass $m_{1}$ or $m_{3}$, and analyze the behavior of
$|m_{\beta\beta}|$ within these regions. Given an allowed sum of neutrino
masses, we determine the value of the minimum neutrino mass. We point out that
if the minimum neutrino mass falls outside the range of [2$\times 10^{-3}$
eV/c2, 7$\times 10^{-3}$ eV/c2], $|m_{\beta\beta}|$ will exceed 1 meV/c2 in
magnitude. This level of sensitivity can be achieved with a 100-ton experiment
over 10 years. This observation elucidates the relationship between neutrino
mass constraints and the effective Majorana mass.

##  II the Sum of Neutrino Masses versus the minimum Neutrino Mass

The most recent cosmological constraints on the sum of neutrino masses cos ,
coupled with the latest laboratory measurements on oscillations pdg , offer
valuable insights into the constraints on the effective Majorana neutrino
mass. In this work, we utilize the squared mass splitting constraints obtained
from a comprehensive global fit to observations of neutrino oscillations pdg .

In the context of the normal neutrino mass hierarchy, often abbreviated as the
normal hierarchy (NH), characterized by $m_{1}<m_{2}<m_{3}$, the parameters
employed in this study are:

$\left\\{\begin{array}[]{c}\Delta m_{21}^{2}=(7.41_{-0.20}^{+0.21})\times
10^{-5}eV^{2}/c^{4}\\\ \\\ \Delta m_{32}^{2}=(2.437_{-0.027}^{+0.028})\times
10^{-3}eV^{2}/c^{4}\\\ \\\ \theta_{13}=(8.54_{-0.12}^{+0.11})^{o}\\\ \\\
sin^{2}\theta_{13}=(2.203_{-0.059}^{+0.056})\times 10^{-2}\\\ \\\
\theta_{23}=(49.1_{-1.3}^{+1.0})^{o}\\\ \\\
sin^{2}\theta_{23}=(5.71_{-0.23}^{+0.18})\times 10^{-1}\\\ \\\
\theta_{12}=(33.41_{-0.72}^{+0.75})^{o}\\\ \\\
sin^{2}\theta_{12}=(3.03_{-0.11}^{+0.12})\times 10^{-1}\end{array}\right\\}$
(7)

Similarly, in the case of the inverted neutrino mass hierarchy, commonly
referred to as the inverted hierarchy (IH), characterized by
$m_{3}<m_{1}<m_{2}$, the parameters utilized in this study are:

$\left\\{\begin{array}[]{c}\Delta m_{21}^{2}=(7.41_{-0.20}^{+0.21})\times
10^{-5}eV^{2}/c^{4}\\\ \\\ \Delta m_{32}^{2}=(-2.498_{-0.025}^{+0.032})\times
10^{-3}eV^{2}/c^{4}\\\ \\\ \theta_{13}=(8.57_{-0.11}^{+0.12})^{o}\\\ \\\
sin^{2}\theta_{13}=(2.219_{-0.057}^{+0.060})\times 10^{-2}\\\ \\\
\theta_{23}=(49.5_{-1.2}^{+0.9})^{o}\\\ \\\
sin^{2}\theta_{23}=(5.78_{-0.21}^{+0.16})\times 10^{-1}\\\ \\\
\theta_{12}=(33.41_{-0.72}^{+0.75})^{o}\\\ \\\
sin^{2}\theta_{12}=(3.03_{-0.11}^{+0.12})\times 10^{-1}\end{array}\right\\}$
(8)

Using the measured values of $\Delta m_{12}^{2}$ and $\Delta m_{32}^{2}$ in
both the NH and IH cases, we can derive the sum $\Sigma=m_{1}+m_{2}+m_{3}$ as
follows:

$\left\\{\begin{array}[]{c}\Sigma=m_{1}+\sqrt{(7.41_{-0.20}^{+0.21})\times
10^{-5}+m_{1}^{2}}+\\\ \sqrt{(7.41_{-0.20}^{+0.21})\times
10^{-5}+(2.437_{-0.027}^{+0.028})\times 10^{-3}+m_{1}^{2}},\\\
\end{array}\right\\}$ (9)

and

$\left\\{\begin{array}[]{c}\Sigma=m_{3}+\sqrt{(2.498_{-0.025}^{+0.032})\times
10^{-3}+m_{3}^{2}}+\\\ \sqrt{(7.41_{-0.20}^{+0.021})\times
10^{-5}+(2.498_{-0.025}^{+0.032})\times 10^{-3}+m_{3}^{2}}\\\
\end{array}\right\\}$ (10)

The oscillation constraints outlined in Equations 9 and 10 can be visualized
in Figure 1, illustrating the permitted region in the $m_{L}-\Sigma$ plane.
Here, $m_{L}$ denotes the minimum neutrino mass, representing $m_{1}$ for NH
and $m_{3}$ for IH.

Figure 1: The sum of masses is represented by the symbol $\Sigma$, defined as
$\Sigma=m_{1}+m_{2}+m_{3}$, and is depicted as a function of the minimum mass
of neutrino mass eigenstates. The estimated sensitivity of KATRIN is
$\Sigma<0.8$ eV (95% C.L.) according to reference kat . Cosmology 2016 data is
sourced from reference ajc , while Cosmology 2021 data is cited from reference
cos .

As depicted in Figure 1, the values of $\Sigma$ fall within two narrow bands:
approximately $0.06$ eV/c2 for NH and $\sim$0.102 eV/c2 for IH. This suggests
a minor influence of the value of the minimum neutrino mass when
$m_{L}<10^{-2}eV/c^{2}$. It is noteworthy that the Particle Data Group pdg  in
2022 provided three additional sets of best-fit parameters. Despite slight
differences among these parameters, they do not alter the conclusion regarding
the narrow allowed bands for $\Sigma$ in both NH and IH cases.

It’s notable that recent searches in cosmology and accelerator experiments
lean towards favoring the NH scenario for neutrino mass raul ; nova ; t2k .
This suggests that the dominant components of $|m_{\beta\beta}|$ would be the
lighter neutrino mass eigenstates, leading to longer decay half-lives and
increased experimental challenges. While the reliability of the cosmological
claim has faced scrutiny sch , recent results from experiments such as NOVA
nova  and T2K t2k  lack statistical significance.

To address the potential preference for the normal hierarchy, Agostini,
Benato, and Detwiler ago  found that the likelihood of detecting
$0\nu\beta\beta$ decay could exceed 50% in the most promising experiments in
the absence of neutrino mass mechanisms that drive the mass of the lightest
state or the effective Majorana mass to zero. Furthermore, there is
consideration of decay facilitated by the exchange of heavy particles petr ;
shao , which involves Majorana neutrinos according to the Black-Box Theorem
jsch , thereby violating lepton number conservation. However, determining the
mass scale of the light neutrinos directly from the data in this scenario
remains challenging.

Therefore, while the adverse implications of NOVA and T2K results may be
overstated, the possibility of exotic heavy particle exchange physics could
diminish their significance. Nevertheless, the primary motivation for
investigating decay—testing lepton number conservation and distinguishing
between the Dirac or Majorana character of neutrinos—remains compelling.

Given the narrow allowed regions of $\Sigma$ in both NH and IH cases, we aim
to explore the parameter space of $|m_{\beta\beta}|$ versus the minimum
neutrino mass $m_{L}$ for both scenarios, particularly when accounting for the
Majorana CP-violating phases.

##  III The parameter Space of $|m_{\beta\beta}|$

In the early 2000s, pioneers developed formulations to calculate
$|m_{\beta\beta}|$ as a function of the minimum neutrino mass fps . Recently,
the coupling-rod diagram of $|m_{\beta\beta}|$, introduced by Zhi-Zhong Xing
and Ye-Ling Zhou xing , illustrates the key characteristics of
$|m_{\beta\beta}|$ in the complex plane. This diagram provides an intuitive
understanding of how the neutrino mass ordering and CP-violating phases impact
$|m_{\beta\beta}|$. Below, we present a summary of the maximum and minimum
values of $|m_{\beta\beta}|$ using Equation 5 for both NH and IH cases.

For the case of NH where $m_{1}<m_{2}<m_{3}$, the equations are:

$\left\\{\begin{array}[]{c}|m_{\beta\beta}|_{max}=\sqrt{\Delta
m_{21}^{2}+m_{1}^{2}}sin^{2}\theta_{12}cos^{2}\theta_{13}[1+\\\
\sqrt{1-\frac{\Delta m_{21}^{2}}{\Delta
m_{21}^{2}+m_{1}^{2}}}cot^{2}\theta_{12}+\\\ \sqrt{1-\frac{\Delta
m_{21}^{2}}{\Delta m_{21}^{2}+m_{1}^{2}}+\frac{\Delta m_{21}^{2}+\Delta
m_{32}^{2}}{\Delta
m_{21}^{2}+m_{1}^{2}}}\frac{tan^{2}\theta_{13}}{sin^{2}\theta_{12}}].\\\ \\\
|m_{\beta\beta}|_{min}^{(1)}=\sqrt{\Delta
m_{21}^{2}+m_{1}^{2}}sin^{2}\theta_{12}cos^{2}\theta_{13}[1-\\\
\sqrt{1-\frac{\Delta m_{21}^{2}}{\Delta
m_{21}^{2}+m_{1}^{2}}}cot^{2}\theta_{12}-\\\ \sqrt{1-\frac{\Delta
m_{21}^{2}}{\Delta m_{21}^{2}+m_{1}^{2}}+\frac{\Delta m_{21}^{2}+\Delta
m_{32}^{2}}{\Delta
m_{21}^{2}+m_{1}^{2}}}\frac{tan^{2}\theta_{13}}{sin^{2}\theta_{12}}].\\\ \\\
|m_{\beta\beta}|_{min}^{(2)}=\sqrt{\Delta
m_{21}^{2}+m_{1}^{2}}sin^{2}\theta_{12}cos^{2}\theta_{13}[\sqrt{1-\frac{\Delta
m_{21}^{2}}{\Delta m_{21}^{2}+m_{1}^{2}}}\times\\\ cot^{2}\theta_{12}-1-\\\
\sqrt{1-\frac{\Delta m_{21}^{2}}{\Delta m_{21}^{2}+m_{1}^{2}}+\frac{\Delta
m_{21}^{2}+\Delta m_{32}^{2}}{\Delta
m_{21}^{2}+m_{1}^{2}}}\frac{tan^{2}\theta_{13}}{sin^{2}\theta_{12}}].\\\
\end{array}\right\\}$ (11)

For the IH case, where $m_{3}<m_{1}<m_{2}$, the equations are as follows:

$\left\\{\begin{array}[]{c}|m_{\beta\beta}|_{max}=\sqrt{m_{3}^{2}-\Delta
m_{32}^{2}}sin^{2}\theta_{12}cos^{2}\theta_{13}[1+\\\ \sqrt{1-\frac{\Delta
m_{21}^{2}}{m_{3}^{2}-\Delta m_{32}^{2}}}cot^{2}\theta_{12}+\\\
\sqrt{1-\frac{\Delta m_{21}^{2}}{m_{3}^{2}-\Delta m_{32}^{2}}+\frac{\Delta
m_{21}^{2}-\Delta m_{32}^{2}}{m_{3}^{2}-\Delta
m_{32}^{2}}}\frac{tan^{2}\theta_{13}}{sin^{2}\theta_{12}}].\\\ \\\
|m_{\beta\beta}|_{min}=\sqrt{m_{3}^{2}-\Delta
m_{21}^{2}}sin^{2}\theta_{12}cos^{2}\theta_{13}[\sqrt{1-\frac{\Delta
m_{21}^{2}}{m_{3}^{2}-\Delta m_{32}^{2}}}\times\\\ cot^{2}\theta_{12}-1-\\\
\sqrt{1-\frac{\Delta m_{21}^{2}}{m_{3}^{2}-\Delta m_{32}^{2}}+\frac{\Delta
m_{21}^{2}-\Delta m_{32}^{2}}{m_{3}^{2}-\Delta
m_{32}^{2}}}\frac{tan^{2}\theta_{13}}{sin^{2}\theta_{12}}].\\\
\end{array}\right\\}$ (12)

We illustrate the relationship between $|m_{\beta\beta}|$ and the minimum
neutrino mass $m_{L}$ in Figure 2, utilizing the parameters specified in
equations 7 and 8. We vary the relevant CP-violating phases within the range
of 0 to $2\pi$.

Figure 2: The variation of $|m_{\beta\beta}|$ with $m_{L}$ in the NH or IH of
three neutrino masses is depicted. The input parameters from equations 7 and 8
are utilized, with the relevant CP-violating phases allowed to range from 0 to
$2\pi$. Additionally, the current sensitivities from KamLAND-Zen and GERDA are
included, along with the planned ton-scale experiments such as LEGEND-1000
legend , nEXO nexo , and CUPID cupid .

In Figure 2, the gray area represents the allowed regions for
$|m_{\beta\beta}|$ concerning $m_{L}$. In the IH scenario, when $m_{L}$ is
below 1 meV/c2, the allowed region for $|m_{\beta\beta}|$ appears nearly flat.
This region is targeted by upcoming ton-scale experiments like LEGEND-1000,
nEXO, and CUPID. Conversely, in the NH scenario, the allowed region exhibits a
pronounced dependence on $m_{L}$. For $m_{L}$ below $\sim$5$\times$10-4 eV/c2,
the $|m_{\beta\beta}|$ region remains almost flat, surpassing 1 meV/c2.
However, within the range [$\sim 2\times 10^{-3}$, $\sim 7\times 10^{-3}$]
eV/c2, $|m_{\beta\beta}|$ tends towards zero.

It’s evident from Figure 2 that the planned ton-scale experiments will
comprehensively explore the IH region. Therefore, our focus will be on
accessing sensitivity for the NH case. We aim to investigate the feasibility
of constructing a sensitive experiment capable of detecting the
0$\nu\beta\beta$ decay process if neutrino masses adhere to normal ordering.

##  IV Discussion on $m_{L}$ and Prospective Experiments Sensitive to the NH
Region

To investigate the NH region for the $0\nu\beta\beta$ decay process, it’s
imperative to avoid $m_{L}$ falling within the range of [$\sim 2\times
10^{-3}$, $\sim 7\times 10^{-3}$] eV/c2, as indicated in Figure 2. With this
consideration, our first step is to examine the value of $m_{L}$ = $m_{1}$ for
a given allowed $\Sigma$. The solution must satisfy the following condition:

$m_{L}+\sqrt{\Delta m_{21}^{2}+m_{L}^{2}}+\sqrt{\Delta m_{32}^{2}+\Delta
m_{21}^{2}+m_{L}^{2}}-\Sigma=0.$ (13)

Equation 13 must be solved numerically for a given value of $\Sigma$. Figure 3
illustrates the relationship between $m_{L}$ and a given $\Sigma$.

Figure 3: The plot illustrates the solution for $m_{1}+m_{2}+m_{3}$-$\Sigma$ =
0 when a $\Sigma$ value is specified. The gray area represents where the value
of $|m_{\beta\beta}|$ approaches zero, while the light-red region indicates
excluded values of $m_{L}$ due to the absence of a solution.

Numerical solutions indicate that when $\Sigma$ ranges between 0.059 and 0.061
eV/c2, Equation 13 is satisfied. This corresponds to $m_{L}$ falling within
the region [3$\times 10^{-4}$, 2$\times 10^{-3}$] eV/c2. However, when
$\Sigma$ = 0.058 eV/c2, Equation 13 is violated, indicating no suitable values
of $m_{L}$ exist and thus excluding that range. Furthermore, for $\Sigma$ in
the range [0.061, 0.068] eV/c2, $m_{L}$ drives $|m_{\beta\beta}|$ to zero, as
depicted in Figure 2.

It’s evident that only a narrow range of $m_{L}$ is viable for determining
$|m_{\beta\beta}|$ in future experiments beyond the ton-scale. This viable
region spans $m_{L}$ = [3$\times 10^{-4}$, 2$\times 10^{-3}$] eV/c2.
Additionally, the distribution of $m_{1}$, $m_{2}$, and $m_{3}$ for a given
$\Sigma$ is visualized in Figure 4.

Figure 4: Displayed is the solution for the neutrino masses $m_{1},m_{2}$, and
$m_{3}$ for a given allowed $\Sigma$. The box outlined by dashed blue lines
indicates the potential non-zero values of $m_{1}$.

In the NH scenario, it’s noteworthy that the sum of masses falls within the
range of 0.059 to 0.061 eV/c2, consequently leading to a non-zero $m_{1}$, as
depicted in Figure 4.

All the results depicted in Figures 3 and 4 indicate that the minimum neutrino
mass is non-zero. This suggests that only the Majorana-violating phases can
drive $|m_{\beta\beta}|$ to zero when $m_{L}$ falls within the range of
$[2\times 10^{-3},7\times 10^{-3}]$ eV/c2. There exists a parameter space for
$|m_{\beta\beta}|$ where its value is in the region of $\sim$ 1 meV/c2,
corresponding to the sensitivity achievable with a 100-ton-scale experiment.

Using 76Ge as an example, the sensitivity to the effective neutrino mass is
linked to the upper limit of the half-life, determined as:

$T_{\frac{1}{2}}=4.17\times 10^{26}(\frac{\varepsilon\cdot a}{W})\frac{M\cdot
t}{\mu},$ (14)

where $T_{\frac{1}{2}}$ is the decay half-life in years, $\varepsilon$ is the
detection efficiency, $a$ is the isotopic abundance (92% for ${}^{76}Ge$), $W$
is the molar mass of the source material, $M$ is the total mass of the source
in kg, $t$ is the run time in years, and $\mu$ is the Poisson signal mean
constrained by the number of background events ($B$) in the region of interest
(ROI). We express $B$ as:

$B=M\cdot t\cdot b\cdot\Delta E,$ (15)

where $b$ is the number of background events per kg per keV per year, and
$\Delta E=3.0$ keV is the width of the ROI.

A 90% confidence level interval for the Poisson signal mean ($\mu$) is
provided in TABLE XII of Ref. Gary  when the total number of events observed
($n_{0}$) and the mean background ($b$) are known.

Using equation 14 and parameters similar to those of LEGEND-1000 legend , we
demonstrate the sensitivity of such an experiment in Figure 5. It’s apparent
that achieving a half-life of $10^{30}$ years with a 10-year exposure requires
a background index of 0.001 events/t-y, which presents a considerable
challenge. This underscores the necessity for significant research and
development (R&D) efforts to be initiated as early as possible.

Figure 5:  The sensitivity for measuring the half-life of 76Ge decay is
illustrated for a future 100-ton scale experiment utilizing similar technology
to LEGEND-1000.

Using equation 1 and defining
$F_{N}=G^{0\nu}(E_{0},Z)|M_{f}^{0\nu}-(\frac{g_{A}}{g_{V}})^{2}M_{GT}^{0\nu}|^{2}$,
if a nuclear structure parameter $F_{N}=7.01\times 10^{14}/y$ is employed 34 ,
the effective neutrino mass $|m_{\beta\beta}|$ can be expressed in terms of
the calculated $F_{N}$ and the measured half-life as:

$|m_{\beta\beta}|=\frac{m_{e}}{\sqrt{F_{N}T_{1/2}^{0\nu}}},$ (16)

where $m_{e}$ is the electron mass. The measurable half-life for a 100-ton
scale 76Ge-based experiment is depicted in Figure 6.

Figure 6:  Displayed is the sensitivity for a prospective 100-ton scale 76Ge
experiment, leveraging technology akin to that of LEGEND-1000.

As depicted in Figures 5 and 6, constructing a 100-ton scale experiment allows
for comprehensive exploration of the $\sim$1 meV/c2 region in the NH scenario,
with the corresponding half-life approaching approximately $10^{30}$ years.
This endeavor entails the preparation of enriched materials over a few years,
alongside ongoing R&D efforts aimed at continuously reducing backgrounds.
Simultaneously, the scientific community is advancing towards ton-scale
experiments to fully investigate the 10 meV/c2 region in the IH scenario.
Notably, the JUNO Collaboration is investigating the feasibility of dissolving
130Te or 136Xe isotopes into liquid scintillator to create a 100-ton scale
experiment aimed at exploring the $|m_{\beta\beta}|\approx 1$ meV/c2 region
cao . Additionally, the DUNE experiment is contemplating the use of a Xe-Doped
liquid argon TPC as a platform for $0\nu\beta\beta$ decay, targeting a
sensitivity of $|m_{\beta\beta}|\approx 2$ meV dune . The critical aspect for
achieving sensitivity to $\sim$1 meV in NH with a 100-ton scale experiment
lies in reducing the background to a level of 1 event per kiloton per year in
the ROI, as illustrated in Figure 5 and 6. Additionally, exceptional energy
resolution is essential for discerning events from 2$\nu\beta\beta$ decays.
Both concepts proposed by JUNO and DUNE will face considerable challenges in
meeting these two critical requirements.

##  V Conclusion

Utilizing constraints from cosmology and neutrino oscillation experiments, we
discuss the permissible parameter space of $|m_{\beta\beta}|$ for both the
inverted hierarchy (IH) and normal hierarchy (NH) scenarios. Given the current
planned ton-scale experiments aimed at addressing the sensitivity across the
entire region of IH, we focus our discussion on the allowed regions of
$|m_{\beta\beta}|$ for NH. It becomes evident that the minimum neutrino mass
($m_{L}$) is non-zero in the NH case. To explore a non-zero
$|m_{\beta\beta}|$, it is imperative to circumvent the region where the
minimum mass lies between 3$\times 10^{-4}$ to 2$\times 10^{-3}$ eV/c2 for NH.
An achievable value of $|m_{\beta\beta}|$ is approximately 1 meV/c2. We
elucidate the importance for future experiments to attain this sensitivity.
With such heightened sensitivity, attainable through a high-capacity
experiment (100-ton scale), the likelihood of observing signal events of
0$\nu\beta\beta$ decays substantially increases.

##  VI Acknowledgement

This work was supported in part by NSF OISE 1743790, NSF PHYS 2310027, DOE DE-
SC0024519, DE-SC0004768, and a research center supported by the State of South
Dakota.

## References

  * (1) The SNO collaboration, Phys. Rev. Lett. 89 (2002) 011301 and nucl-ex/0204008.
  * (2) The KamLAND collaboration,PRL v90, 021802,(2003).
  * (3) The Super-Kamiokande Collaboration, Phys. Rev. Lett. 81 (1998) 1562-1567.
  * (4) Majorana, E. Theory of the symmetry of electrons and positrons [in Italian]. Nuovo Cim. 14, 171–184 (1937).
  * (5) Steve R. Elliott and Petr Vogel, Annu. Rev. Nucl. Part. Sci. 115 (2002)
  * (6) Steven R. Elliott and Jonathan Engel, J. Phys. G: Nucl. Part. Phys. 30, R183 (2004).
  * (7) F.T Avignone III, G.S. King III and Yuri Zdesenko, New Journal of Physics 7, 6 (2005).
  * (8) A.S. Barabash, Physics of Atomic Nuclei, 67, No. 3, 438 (2004).
  * (9) B. Pontecorvo, Zh. Eksp. Teor. Fiz. 33, 549 (1957); 34, 247 (1957); Z. Maki, M. Nakagawa, and Sakata, Prog. Thror. Phys. 28, 870 (1962).
  * (10) R.L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022) and 2023 update 1st December, 2023.
  * (11) Matteo Agostino, Giovanni Benato, and Jason Detwiler, Phys. Rev. D 96, 053001 (2017).
  * (12) S. Abe at el. (KamLAND-Zen Collaboration), Phys. Rev. Lett. 130, 051801.
  * (13) M. Agostini et al. (GERDA Collaboration), Eur. Phys. L. C 2023, 83(9):778.
  * (14) S. Alam, M. Aubert, S. Avila, C. Balland, J.E. Bautista, M.A. Bershady et al., Completed sdss-iv extended baryon oscillation spectroscopic survey: Cosmological implications from two decades of spectroscopic surveys at the apache point observatory, Physical Review D 103 (2021).
  * (15) Katrin Collaboration, M. Aker, A. Beglarian, J. Behrens, A. Berlev, U. Besserer et al., Direct neutrino-mass measurement with sub-electronvolt sensitivity, Nature Physics 18 (2022) 160.
  * (16) A.J. Cuesta, V. Niro and L. Verde, Neutrino mass limits: Robust information from the power spectrum of galaxy surveys, Physics of the Dark Universe 13 (2016) 77 [1511.05983].
  * (17) Raul Jimenez et al., JCAP09(2022)006.
  * (18) P. Adamson et al. (NOVA Collaboration), Phys. Rev. Lett. 118, 231801 (2017).
  * (19) K. Abe et al. (T2K Collaboration) Phys. Rev. Lett. 118, 151801 (2017).
  * (20) T. Schwetz et al., JCAP 03, 011 (2018).
  * (21) Gary Prézeau, Phys. Lett. B 633, 93 (2006).
  * (22) Shao-Feng Ge et al., Phys. Rev. D 96, 055019 (2017).
  * (23) J. Schechter and J.W. Valle, Phys. Rev. 25, 2951 (1982).
  * (24) Y. Farzan, O. L. G. Peres, A. Yu, Smirnov, Nucl.Phys.B612:59-97,2001.
  * (25) Xing Zhi-zhong and Zhou Ye-Ling 2015 Chinese Phys. C 39 011001.
  * (26) N. Abgrall et al. (LEGEND Collaboration), Cite as: arXiv:2107.11462.
  * (27) T. Stiegler et al. (nEXO Collaboration), arXiv: 2009.10231.
  * (28) K. Alfonso et al. (CUPID Collaboration), Journal of Low Temperature Physics, Vol. 211, Issue 5-6; ISSN 0022-2291.
  * (29) Gary J. Feldman and Robert D. Cousins, Phys. Rev. D 57 3873 (1998).
  * (30) Jun Cao et al., Chinese Physics C Vol. 44, No. 3 (2020) 031001.
  * (31) J. Suhonen and O. Civitarese, Phys. Rep. 300, 124 (1998).
  * (32) A. Mastbaum, F. Psihas, and J. Zennamo, Phys. Rev. D 106, 092002 (2022).


# CubismAMR - A C$++$ library for Distributed Block-Structured Adaptive Mesh
Refinement

Michail Chatzimanolakis, Pascal Weber,
Fabian Wermelinger, Petros Koumoutsakos

###### Abstract

We present CubismAMR, a C$++$ library for distributed simulations with block-
structured grids and Adaptive Mesh Refinement. A numerical method to solve the
incompressible Navier-Stokes equations is proposed, that comes with a novel
approach of solving the pressure Poisson equation on an adaptively refined
grid. Validation and verification results for the method are presented, for
the flow past an impulsively started cylinder.

##  1 Introduction

Fluid flow simulations involve a multitude of scales and quickly exhaust the
capabilities of simulations that employ uniform grids. Adaptive Mesh
Refinement (AMR) tackles the issue by dynamically adapting the grid to capture
the emerging structures in a flow field. Simulations using AMR that
efficiently exploit modern, massively parallel computer architectures are a
subject of ongoing investigations. Existing AMR software that targets large-
scale simulations include AMReX (Zhang et al., 2019), Chombo (Adams et al.,
2019), SAMRAI (Wissink et al., 2001), and Basilisk (Van Hooft et al., 2018).
Patch-based AMR dynamically groups cells to form rectangular patches in AMReX,
Chombo, and SAMRAI. In contrast, Basilisk uses cell-based AMR by organizing
the cells in an octree. The AMR framework presented here is based on the
Gordon-Bell-Prize winning Cubism library that partitions the simulation domain
into cubes and distributes them to multiple compute nodes. These cubes are
divided into blocks for cache-optimized thread-level parallelism
(Hejazialhosseini et al., 2012; Rossinelli et al., 2013; Hadjidoukas et al.,
2015). CubismAMR organizes this block-structure into an octree to perform
compression and refinement in different regions of space.

In the present paper we present a numerical method for solving the
incompressible Navier-Stokes equations through the use of CubismAMR. The paper
is structured as follows: section 2 describes the governing equations and
their discretization. Section 3 presents the AMR algorithm we employ, while
section 4 shows the validation and verification of our method.

##  2 Governing Equations

We solve the incompressible Navier-Stokes equations in the velocity-pressure
formulation

$\begin{split}&\mathbf{\nabla}\cdot\boldsymbol{u}=0\,,\\\
&\frac{\partial\boldsymbol{u}}{\partial
t}+(\boldsymbol{u}\cdot\mathbf{\nabla})\boldsymbol{u}=-\frac{1}{\rho}\nabla
p+\nu\mathbf{\nabla}^{2}\boldsymbol{u},\end{split}$ (1)

where $\boldsymbol{u},p$ and $\nu$ are the fluid velocity, pressure and
kinematic viscosity. The no-slip boundary condition is enforced on the surface
of moving bodies with a prescribed velocity $\boldsymbol{u}^{s}$ through the
penalisation approach (Angot et al., 1999; Coquerelle and Cottet, 2008;
Gazzola et al., 2011). The penalisation method augments the governing Navier-
Stokes equations by a penalty term
$\lambda\chi(\boldsymbol{u}^{s}-\boldsymbol{u})$ where $\chi$ is a
characteristic function with values $\chi=1$ inside the cylinder and $\chi=0$
outside. As the penalisation coefficient, denoted by $\lambda$, tends to
infinity, eq. 1 converge to the incompressible Navier-Stokes equations. Ueda
and Kida (2021) demonstrate that the tangential and slip velocities on the
surface of the solid body are in the order of $\lambda^{-1/2}$ and
$\lambda^{-1}$ respectively.

###  2.1 Discretisation

Equation 1 is solved using a pressure projection method (Chorin, 1968) and a
second-order time-stepping scheme. These methods are adapted to take into
account the penalisation force terms through a split-step algorithm.
Quantities measured at time $t=n\Delta t$ will have a superscript $n$. Given
the numerical solution at timestep $n$, an intermediate velocity
$\boldsymbol{u}^{*}$ is computed with the midpoint method

$\begin{split}&\boldsymbol{u}^{n+1/2}=\boldsymbol{u}^{n}+\frac{1}{2}\Delta
t\mathbf{F}(\boldsymbol{u}^{n})\,,\\\
&\boldsymbol{u}^{*}=\boldsymbol{u}^{n}+\Delta
t\mathbf{F}(\boldsymbol{u}^{n+1/2})\,,\\\ \end{split}$ (2)

where
$\mathbf{F}(\boldsymbol{u})=-(\boldsymbol{u}\cdot\mathbf{\nabla})\boldsymbol{u}+\nu\mathbf{\nabla}^{2}\boldsymbol{u}$.
Diffusion terms present in $\mathbf{F}(\boldsymbol{u})$ are discretised with
centered second-order finite differences and advection terms are handled with
an upwind fifth-order WENO scheme (Shu, 1999).

A second intermediate velocity is computed by performing an implicit
penalisation step

$\boldsymbol{u}^{**}=\boldsymbol{u}^{*}+\lambda\Delta
t\chi^{n+1}(\boldsymbol{u}^{s,n+1}-\boldsymbol{u}^{**})\,,$ (3)

where the discontinuous characteristic function $\chi^{n+1}$ is discretised by
using a second-order accurate mollification, proposed by Towers (2009).
Equation 3 is stable for any $\lambda>0$, allowing for arbitrarily large
values for the penalisation coefficient. The timestep is concluded with a
pressure correction step

$\boldsymbol{u}^{n+1}=\boldsymbol{u}^{**}-\Delta t\nabla p^{n+1}\,,$ (4)

where the pressure field is obtained by solving the Poisson equation

$\nabla^{2}p^{n+1}=\frac{1}{\Delta t}\nabla\cdot\mathbf{u}^{**}\,.$ (5)

The equations are solved using block-structured grids. The present work relies
on the CubismAMR software that is an adaptive version of the Cubism library,
which partitions the simulation domain into cubes of uniform resolution that
are distributed to multiple compute nodes. These cubes are further divided
into blocks for cache-optimised parallelism (Hejazialhosseini et al., 2012;
Rossinelli et al., 2013; Hadjidoukas et al., 2015). CubismAMR organizes these
blocks in an octree data structure (for three-dimensional simulations) or a
quadtree data structure (for two-dimensional simulations), which allows for
grid refinement or compression in different regions. In contrast to uniform
grids or body fitted meshes, this allows to dynamically adopting the grid to
capture the emerging structures in a flow field. CubismAMR fully supports two-
and three-dimensional distributed simulations. What follows in section 3 only
describes the two-dimensional case, which is relevant to the simulations of
this paper.

###  2.2 Force computation

The total force acting on an solid body is

$\boldsymbol{F}=\int\limits_{\partial\Omega_{s}}(2\mu\mathbf{D}\cdot\boldsymbol{n}-p\boldsymbol{n})\mathrm{d}S\,,$
(6)

where $\mathrm{d}S$ denotes the infinitesimal surface element with normal
$\boldsymbol{n}$, $\mu$ is the dynamic viscosity, $\Omega_{s}\in\Omega$ is the
subset of the simulation domain occupied by the body and
$\mathbf{D}=\frac{1}{2}\left(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{\top}\right)$
the strain-rate tensor. The first term corresponds to viscous forces, and the
second to pressure-induced forces.

Following Towers (2008) the surface integral in eq. 6 is expressed as a volume
integral

$\boldsymbol{F}=\int\limits_{\Omega}(2\mu\mathbf{D}\cdot\boldsymbol{n}-p\boldsymbol{n})\delta(\operatorname{S_{d}})\;\mathrm{d}\Omega\,,$
(7)

where $\delta$ is the Dirac delta and $\operatorname{S_{d}}$ is the signed
distance function from the surface of the body to any point in $\Omega$. Since
$\chi^{(s)}=H(\operatorname{S_{d}})$, where $H$ is the Heaviside function, we
find that
$\delta(\operatorname{S_{d}})=\frac{\partial\chi}{\partial\mathbf{n}}=\nabla\chi\cdot\boldsymbol{n}$
where the normal vector is computed from the signed distance function as
$\boldsymbol{n}=\frac{\nabla\operatorname{S_{d}}}{|\nabla\operatorname{S_{d}}|}\big{|}_{\partial\Omega_{s}}$.
Thus, the total force is computed from

$\boldsymbol{F}=\int\limits_{\Omega}(2\mu\mathbf{D}\cdot\boldsymbol{n}-p\boldsymbol{n})(\nabla\chi^{(s)}\cdot\boldsymbol{n})\;\mathrm{d}\Omega\,.$
(8)

When computing viscous forces, Verma et al. (2017) observed that velocity
gradients near walls can be underestimated, when a penalisation method is
used. Therefore, they suggested computing the gradients on a “lifted” surface.
Here we employ a similar approach: the necessary gradients are computed two
grid points away from the surface and are then extrapolated back to it through
a second-order Taylor expansion

$\frac{\partial u}{\partial x_{j}}=\frac{\partial u_{L}}{\partial
x_{j}}+\Delta x_{i}\frac{\partial^{2}u_{L}}{\partial x_{j}\partial
x_{i}}~{},~{}j=1,2\,,$ (9)

where the subscript $L$ is used to denote a quantity on the lifted surface,
$\Delta x_{j}$ is distance in the $j$ direction of the grid point of the
actual surface from the grid point on the lifted surface. All derivatives are
computed with second-order one-sided differences, facing away from the wall.
From this the drag $F_{D}$ and drag coefficient $C_{D}$ can be computed as the
projection of the total force onto the direction of motion

$F_{D}=\boldsymbol{F}\cdot\frac{\mathbf{u}^{(s)}}{|\mathbf{u}^{(s)}|}\,,\qquad
C_{D}=\frac{2F_{D}}{\rho|\mathbf{u}^{(s)}|^{2}A}\,,$ (10)

where $A$ is a characteristic area and $\rho=1$ is the fluid density.

##  3 Block-Structured AMR

The grid is composed of square blocks, each with the same number of cells.
Each block has a locally uniform resolution that is defined by its level of
refinement $l=0,\dots,L-1$ as $h_{l}=2^{-l}h_{0}$, where $h_{0}$ is the
coarsest grid spacing possible. At mesh refinement from level $l$ to level
$l+1$ a block is divided into four blocks , whereas mesh compression from
level $l$ to $l-1$ is achieved by combining four blocks to one. Blocks that
are adjacent are not allowed to differ by more then one refinement level.
Therefore, the blocks of the grid are logically arranged in a quadtree data
structure.

Whether a block should be refined or compressed is determined every few
timesteps. Our simulations use the vicinity to solid body surfaces and the
magnitude of vorticity as criteria for mesh refinement or compression.
Additional examples of such criteria could be the magnitude of pressure
gradients, or the magnitude of Wavelet detail coefficients (Rossinelli et al.,
2011; Vasilyev and Kevlahan, 2005). Following the multiresolution framework by
Harten (1994, 1996), restriction and prolongation operators are used to map
values between blocks of different resolutions. Let

$f_{i,j}^{l}=f\left(\frac{h_{l}}{2}+ih_{l},\frac{h_{l}}{2}+jh_{l}\right)$ (11)

with indices $i,j\in\mathbb{Z}$ denote the value of
$f:\mathbb{R}^{2}\rightarrow\mathbb{R}$ at a cell center. The restriction
operator $\mathcal{R}$ is used during mesh compression to replace four blocks
at level $l+1$ by one block at level $l$. Grid point values at level $l$ are
computed by averaging (see fig. 1)

$f^{l}_{i,j,k}=\frac{1}{4}\sum_{I=0,1}\sum_{J=0,1}f^{l+1}_{2i+I,2j+J}+{\cal
O}(h_{l}^{2})\,.$ (12)

The prolongation operator $\mathcal{I}$ of one block at level $l$ to four
blocks at level $l+1$ is defined via a third-order Taylor expansion

$f^{l+1}_{2i+I,2j+J}=f^{l}_{i,j}+\frac{h_{l}}{4}\boldsymbol{d}\cdot\nabla
f^{l}_{i,j}+\frac{h_{l}^{2}}{32}(\boldsymbol{d}\cdot\nabla)^{2}f^{l}_{i,j}+{\cal
O}(h_{l}^{3})\,,$ (13)

where $\boldsymbol{d}=(2\delta_{I0}-1,2\delta_{J0}-1)\text{ with
}I,J\in\\{0,1\\}$ and $\delta$ is the Kronecker delta. Derivatives are
approximated with second-order central finite difference schemes.

###  3.1 Coarse-Fine Interfaces

Figure 1: Illustration of the AMR. Figure a) illustrates of how cells are
refined or compressed. b) shows the fluxes at coarse-fine-interface. In order
to ensure conservation, the sum of green fluxes is replacing the blue flux.
Figures c) illustrate load-balancing using diffusion of work. The black line
illustrates the Hilbert curve. The colour represents the process that owns the
block. The block in dark green got refined and thus the green process has more
work than the other three processes. After diffusion of work the load is again
equally distributed.

For the approximation of spatial derivatives, we use finite difference
schemes. This is done by creating a frame of uniform resolution around each
gridpoint (Rossinelli et al., 2015) through the interpolation of ghost cell
values. At the interface between different refinement levels $l$ and $l+1$,
this requires the interpolation of ghost cell values from the coarse to the
fine level and vice-versa. From a fine to a coarse level this is done by
averaging as described in eq. 12. For the ghost cells that need to be
interpolated from a coarse level to a finer one, we use third-order accurate
quadratic interpolation as proposed by Adams et al. (2019) and Almgren et al.
(1998). When used with a second-order finite difference scheme for
approximation of first and second derivatives, this guarantees second and
first order accuracy respectively.

Whenever needed we follow Berger and Oliger (1984) and Berger and Colella
(1989) and use a second-order accurate conservative discretisation of the
divergence operator $\nabla\cdot\boldsymbol{u}$, by expressing it as the sum
of fluxes passing through four faces of a cell

$\begin{split}&\nabla\cdot\boldsymbol{u}|_{i,j}^{l}=\frac{F_{i+1/2,j}^{l}-F_{i-1/2,j}^{l}}{h_{l}}+\frac{G_{i,j+1/2}^{l}-G_{i,j-1/2}^{l}}{h_{l}}+{\cal
O}(h_{l}^{2})\,,\end{split}$ (14)

where $F_{i+1/2,j}^{l}=\frac{u^{l}_{i+1,j}+u^{l}_{i,j}}{2}$. Similar
expressions are used for the other cell faces. Whenever a cell is next to
cells of different resolutions, missing values are interpolated. This results
in a non-conservative discretisation of the divergence operator at the
interface between two different refinement levels $l$ and $l+1$, as the flux
computed for level $l$ need not be equal to the sum of the two fluxes at level
$l+1$ that make up the same cell face. The situation is illustrated in fig. 1.
Conservation is achieved by replacing the flux at level $l$ by the sum of the
fluxes at level $l+1$. For the fluxes in the +$x$ direction this reads

$F_{i+1/2,j}^{l}=\frac{1}{2}\left(F_{2i,2j}^{l+1}+F_{2i,2j+1}^{l+1}\right)\,.$
(15)

The computation in the other directions is similar.

###  3.2 Parallelisation

CubismAMR is parallelised with the Message Passing Interface (MPI) programming
model. The quadtree data structure is traversed by a space-filling Hilbert
curve (Hilbert, 1935), that assigns each block a unique index. Initially,
blocks are distributed to different MPI processes based on their index along
the Hilbert curve, the locality property of which guarantees that each process
will own blocks that are spatially close to each other. Load-imbalance between
different processes, introduced because of mesh compression or refinement, is
handled by redistributing work through a one-dimensional diffusion-based
scheme. As was first proposed by Cybenko (1989), the number of blocks
$N_{p}^{t}$ of process $p$ at timestep $t$ is updated as

$N^{t+1}_{p}=N^{t}_{p}+c(N_{p+1}^{t}-2N^{t}_{p}+N_{p-1}^{t})\,.$ (16)

$c$ is a user-defined constant, set to $c=0.25$ in the present work. This
load-balancing scheme limits communication between consecutive processes along
the one-dimensional Hilbert curve while gradually redistributing workload. In
fig. 1 (bottom) we illustrate the process in two-dimensions. In addition to
this, we track the load-imbalance ratio (defined as the ratio between the
maximum and the minimum number of blocks any process has). When it exceeds a
user-defined threshold (set to 1.1 for all the presented results), all blocks
are evenly redistributed to all processes.

###  3.3 Solving the Pressure Equation on the Adaptively Refined Mesh

We write eq. 5 as

$\nabla^{2}\phi=\frac{1}{\Delta
t}\nabla\cdot\mathbf{u}^{**}-\nabla^{2}p^{n}\,,$ (17)

where $\phi=p^{n+1}-p^{n}$. Equation 17 is discretised with a conservative,
second-order accurate discretisation of the divergence operator (Martin et
al., 2008). All unknown values are concatenated in a vector
$\mathbf{\phi}=(\mathbf{\phi}_{1},\dots,\mathbf{\phi}_{B})\in\mathbb{R}^{Bm^{2}}$,
where $\mathbf{\phi}_{i}\in\mathbb{R}^{m^{2}}$ corresponds to values found in
the $m^{2}$ grid points owned by block $i$ (assuming a mesh that is composed
of $B$ square blocks of locally uniform resolution with $m^{2}$ grid points
each). The arising linear system $A\boldsymbol{\phi}=\boldsymbol{b}$ is solved
by using the preconditioned biconjugate gradient stabilised method (van der
Vorst, 1992), with a custom preconditioner.

The preconditioner is constructed by first defining the Poisson matrix

$\mathbf{M}=\begin{bmatrix}\mathbf{P}&-\mathbf{I}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}\\\
-\mathbf{I}&\mathbf{P}&-\mathbf{I}&\mathbf{0}&\mathbf{0}&\dots&\mathbf{0}\\\
\mathbf{0}&-\mathbf{I}&\mathbf{P}&-\mathbf{I}&\mathbf{0}&\dots&\mathbf{0}\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\mathbf{0}&\dots&\mathbf{0}&-\mathbf{I}&\mathbf{P}&-\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&\dots&\dots&\mathbf{0}&-\mathbf{I}&\mathbf{P}&-\mathbf{I}\\\
\mathbf{0}&\mathbf{0}&\dots&\dots&\mathbf{0}&-\mathbf{I}&\mathbf{P}\end{bmatrix},~{}~{}\mathbf{P}=\begin{bmatrix}4&-1&0&0&0&\dots&0\\\
-1&4&-1&0&0&\dots&0\\\ 0&-1&4&-1&0&\dots&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\ 0&\dots&0&-1&4&-1&0\\\
0&\dots&\dots&0&-1&4&-1\\\ 0&0&\dots&\dots&0&-1&4\end{bmatrix}\,,$

where $\mathbf{0}\in\mathbb{R}^{m\times m}$ is a zero matrix,
$\mathbf{I}\in\mathbb{R}^{m\times m}$ is the identity matrix and
$\mathbf{P}\in\mathbb{R}^{m\times m}$. This describes the discretisation of
the Poisson equation on an $m\times m$ uniform grid with zero Dirichlet
boundary conditions, using centered, second-order accurate finite differences.
Then, we define the block diagonal matrix

$\mathbf{A}=\text{diag}(\mathbf{M},\dots,\mathbf{M})\,.$ (18)

The linear system is preconditioned with $\mathbf{A}^{-1}$. Note that
$\mathbf{M}$ is the same for all blocks and it is therefore sufficient to
compute its Cholesky decomposition once, and use it whenever it is needed to
invert $\mathbf{A}$. The chosen preconditioner results in the solution of a
smaller Poisson equation with homogeneous Dirichlet boundary conditions per
block. The (per block) boundary condition $\phi=0$ is equivalent to assuming
$p^{n+1}=p^{n}$.

##  4 Verification and Validation

We simulate the impulsively started flow around a cylinder, for Reynolds
numbers in the range $[550,9~{}500]$. All the simulations presented in this
section use a reference time unit $T=\frac{D}{2U}$, where $D$ is the cylinder
diameter and $U$ its velocity. They were performed in a $[0,40D]\times[0,20D]$
domain and the cylinder was placed at $(10D,10D)$. The timestep was controlled
by a Courant number of $0.45$ , based on the minimal grid spacing present in
the domain. The coarsest grid possible had a resolution of $128\times 64$; the
finest resolution varied per case. The penalisation coefficient was set to
$\lambda=10^{7}$. Grid refinement and compression were based on vorticity
magnitude and the vicinity to the cylinder surface (regions within a distance
of $0.1D$ of non-zero values of $\nabla\chi$ were refined to the maximum
resolution allowed).

|  | Drag coefficient | Vorticity
---|---|---|---
$L$ | $\hat{N}^{L}$ | $\varepsilon_{C_{D}}^{L}$ | Rate | $\varepsilon_{\omega}^{L}$ | Rate
4 | 89 | 0.167 | 1.01 | 0.524 | 1.78
5 | 119 | 0.124 | 1.97 | 0.313 | 2.26
6 | 166 | 0.0469 | 2.94 | 0.125 | 2.63
7 | 242 | 0.0154 | 2.23 | 0.0486 | 2.04
8 | 381 | 0.00779 | 1.51 | 0.0238 | 1.59
9 | 637 | 0 | $-$ | 0 | $-$
Table 1: Verification study errors and convergence rates, impulsively started
cylinder at $\textit{Re}=1~{}000$.

Figure 2: Convergence of vorticity field with increasing resolution at
$\textit{Re}=1~{}000$. From top left to bottom right: 6,7,8 and 9 levels of
refinement.

Figure 3: Drag coefficient, validation results. Comparison of early drag
history with analytical solution by Collins and Dennis (1973) (left) and with
simulations by Koumoutsakos and Leonard (1995) (right), for various Reynolds
numbers.

To verify consistency and convergence, we fix the Reynolds number at
$\textit{Re}=1~{}000$ and perform simulations with increasing resolution, by
changing the maximum levels of refinement allowed, denoted by $L$. The
simulation with the finest resolution (with $L=L^{*}=9$ levels) is used as a
reference solution for two quantities of interest: the instantaneous cylinder
drag coefficient $C_{D}^{L}(t)$ and the vorticity field $\omega^{L}(\tau)$
with $\tau=10T$. For the drag, the error is defined as

$\varepsilon_{C_{D}}^{L}=\frac{1}{t_{\text{max}}-t_{\text{min}}}\int_{t_{\text{min}}}^{t_{\text{max}}}|C_{D}^{L}(t)-C_{D}^{L^{*}}(t)|\mathrm{d}t\,,$
(19)

with $t_{\text{min}}=0.01T$ and $t_{\text{max}}=10T$ and for the vorticity
field as

$\varepsilon_{\omega}^{L}=|\omega^{L}(\tau)-\omega^{L^{*}}(\tau)|~{}~{},~{}~{}\tau=10T\,.$
(20)

Following Bergdorf and Koumoutsakos (2006) we also define the average number
of gridpoints per dimension

$\hat{N}^{L}=\left(\frac{1}{t_{\text{max}}-t_{\text{min}}}\int_{t_{\text{min}}}^{t_{\text{max}}}N^{L}(t)\right)^{1/2}\,,$
(21)

as a measure of the computational cost, where $N^{L}(t)$ is the instantaneous
number of grid points when $L$ refinement levels are used. A summary of the
simulations performed for this verification study is shown in table 1.
Convergence of the vorticity field is visualised in fig. 2, where it can be
seen that the solutions for $8$ and $9$ refinement levels are
indistinguishable from each other.

To validate our method, we compare our results against the drag coefficient of
the impulsively started cylinder at early times as computed analytically by
Collins and Dennis (1973) and against simulations by Koumoutsakos and Leonard
(1995), for longer times. As shown in fig. 3, our results are in excellent
agreement with the references they are compared against.

##  5 Conclusion and Outlook

We presented CubismAMR, a high-performance C$++$ library for distributed
simulations with Adaptive Mesh Refinement in two and three dimensions. The
library was used with a two-dimensional incompressible flow solver that
includes a method of solving the pressure Poisson equation on an adaptively
refined mesh. Validation and verification results were shown for the flow past
an impulsively started cylinder at several Reynolds numbers that are in
excellent agreement with previous simulations and theoretical work. Future
work involves making CubismAMR open-source, to provide a reliable tool for
making distributed, large-scale AMR simulations.

## References

  * Adams et al. (2019) M. Adams, P. Colella, D. T. Graves, J.N. Johnson, N.D. Keen, T. J. Ligocki. D. F. Martin. P.W. McCorquodale, D. Modiano. P.O. Schwartz, T.D. Sternberg, and B. Van Straalen.  Chombo software package for amr applications - design document.  _Lawrence Berkeley National Laboratory Technical Report_ , LBNL-6616E, 2019.  URL https://commons.lbl.gov/download/attachments/73468344/chomboDesign.pdf?version=1&modificationDate=1554672305006&api=v2.
  * Almgren et al. (1998) Ann S. Almgren, John B. Bell, Phillip Colella, Louis H. Howell, and Michael L. Welcome.  A conservative adaptive projection method for the variable density incompressible navier–stokes equations.  _Journal of Computational Physics_ , 142(1):1–46, 1998.  ISSN 0021-9991.  doi: https://doi.org/10.1006/jcph.1998.5890.  URL https://www.sciencedirect.com/science/article/pii/S0021999198958909.
  * Angot et al. (1999) P Angot, C H Bruneau, and P Fabrie.  A penalization method to take into account obstacles in incompressible viscous flows.  _Numerische Mathematik_ , 81(4):497–520, 1999\.
  * Bergdorf and Koumoutsakos (2006) Michael Bergdorf and Petros Koumoutsakos.  A lagrangian particle-wavelet method.  _mms_ , 5:–, 09 2006.  doi: 10.1137/060652877.
  * Berger and Oliger (1984) Marsha J Berger and Joseph Oliger.  Adaptive mesh refinement for hyperbolic partial differential equations.  _Journal of Computational Physics_ , 53(3):484 – 512, 1984.  ISSN 0021-9991.  doi: https://doi.org/10.1016/0021-9991(84)90073-1.  URL http://www.sciencedirect.com/science/article/pii/0021999184900731.
  * Berger and Colella (1989) M.J. Berger and P. Colella.  Local adaptive mesh refinement for shock hydrodynamics.  _Journal of Computational Physics_ , 82(1):64 – 84, 1989.  ISSN 0021-9991.  doi: https://doi.org/10.1016/0021-9991(89)90035-1.  URL http://www.sciencedirect.com/science/article/pii/0021999189900351.
  * Chorin (1968) Alexandre Joel Chorin.  Numerical Solution of the Navier-Stokes Equations.  _Mathematics of Computation_ , 22(104):745–762, 1968.  URL https://www.jstor.org/stable/2004575.
  * Collins and Dennis (1973) W. M. Collins and S. C. R. Dennis.  Flow past an impulsively started circular cylinder.  _Journal of Fluid Mechanics_ , 60(1):105–127, 1973.  doi: 10.1017/S0022112073000078.
  * Coquerelle and Cottet (2008) M. Coquerelle and G.-H. Cottet.  A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies.  _Journal of Computational Physics_ , 227(21):9121 – 9137, 2008.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2008.03.041.  URL http://www.sciencedirect.com/science/article/pii/S0021999108001976.  Special Issue Celebrating Tony Leonard’s 70th Birthday.
  * Cybenko (1989) George Cybenko.  Dynamic load balancing for distributed memory multiprocessors.  _Journal of Parallel and Distributed Computing_ , 7(2):279 – 301, 1989.  ISSN 0743-7315.  doi: https://doi.org/10.1016/0743-7315(89)90021-X.  URL http://www.sciencedirect.com/science/article/pii/074373158990021X.
  * Gazzola et al. (2011) Mattia Gazzola, Philippe Chatelain, Wim M. van Rees, and Petros Koumoutsakos.  Simulations of single and multiple swimmers with non-divergence free deforming geometries.  _Journal of Computational Physics_ , 230(19):7093 – 7114, 2011.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2011.04.025.  URL http://www.sciencedirect.com/science/article/pii/S0021999111002737.
  * Hadjidoukas et al. (2015) Panagiotis E. Hadjidoukas, Diego Rossinelli, Babak Hejazialhosseini, and Petros Koumoutsakos.  From 11 to 14.4 pflops: Performance optimization for finite volume flow solver.  In _Proceedings of the 3rd International Conference on Exascale Applications and Software_ , page 7–12. University of Edinburgh, 2015.  ISBN 9780992661519.
  * Harten (1994) Ami Harten.  Adaptive multiresolution schemes for shock computations.  _Journal of Computational Physics_ , 115(2):319 – 338, 1994.  ISSN 0021-9991.  doi: https://doi.org/10.1006/jcph.1994.1199.  URL http://www.sciencedirect.com/science/article/pii/S0021999184711995.
  * Harten (1996) Ami Harten.  Multiresolution representation of data: A general framework.  _SIAM Journal on Numerical Analysis_ , 33(3):1205–1256, 1996.  ISSN 00361429.  URL http://www.jstor.org/stable/2158503.
  * Hejazialhosseini et al. (2012) B. Hejazialhosseini, D. Rossinelli, C. Conti, and P. Koumoutsakos.  High throughput software for direct numerical simulations of compressible two-phase flows.  In _SC ’12: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_ , pages 1–12, 2012.  doi: 10.1109/SC.2012.66.
  * Hilbert (1935) David Hilbert.  _Über die stetige Abbildung einer Linie auf ein Flächenstück_ , pages 1–2.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1935.  ISBN 978-3-662-38452-7.  doi: 10.1007/978-3-662-38452-7˙1.  URL https://doi.org/10.1007/978-3-662-38452-7_1.
  * Koumoutsakos and Leonard (1995) P. Koumoutsakos and A. Leonard.  High-resolution simulations of the flow around an impulsively started cylinder using vortex methods.  _Journal of Fluid Mechanics_ , 296:1–38, 1995.  doi: 10.1017/S0022112095002059.
  * Martin et al. (2008) Daniel Martin, Phillip Colella, and Daniel Graves.  A cell-centered adaptive projection method for the incompressiblenavier-stokes equations in three dimensions.  _Journal of Computational Physics_ , 227:1863–1886, 01 2008\.  doi: 10.1016/j.jcp.2007.09.032.
  * Rossinelli et al. (2011) Diego Rossinelli, Babak Hejazialhosseini, Daniele G. Spampinato, and Petros Koumoutsakos.  Multicore/multi-GPU accelerated simulations of multiphase compressible flows using wavelet adapted grids.  _SIAM J. Sci. Comput._ , 33(2):512–540, jan 2011.  doi: 10.1137/100795930.  URL http://www.cse-lab.ethz.ch/wp-content/papercite-data/pdf/rossinelli2011b.pdf.
  * Rossinelli et al. (2013) Diego Rossinelli, Babak Hejazi, Panagiotis Hadjidoukas, Costas Bekas, Alessandro Curioni, Adam Bertsch, Scott Futral, Steffen Schmidt, Nikolaus Adams, and Petros Koumoutsakos.  11 pflop/s simulations of cloud cavitation collapse.  11 2013.  doi: 10.1145/2503210.2504565.
  * Rossinelli et al. (2015) Diego Rossinelli, Babak Hejazialhosseini, Wim van Rees, Mattia Gazzola, Michael Bergdorf, and Petros Koumoutsakos.  MRAG-i2d: Multi-resolution adapted grids for remeshed vortex methods on multicore architectures.  _J. Comput. Phys._ , 288:1–18, may 2015.  doi: 10.1016/j.jcp.2015.01.035.  URL http://www.cse-lab.ethz.ch/wp-content/papercite-data/pdf/rossinelli2015a.pdf.
  * Shu (1999) Chi-Wang Shu.  _High Order ENO and WENO Schemes for Computational Fluid Dynamics_ , pages 439–582.  Springer Berlin Heidelberg, Berlin, Heidelberg, 1999.
  * Towers (2008) John D. Towers.  A convergence rate theorem for finite difference approximations to delta functions.  _Journal of Computational Physics_ , 227(13):6591–6597, 2008.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2008.03.019.  URL https://www.sciencedirect.com/science/article/pii/S0021999108001708.
  * Towers (2009) John D. Towers.  Finite difference methods for approximating Heaviside functions.  _Journal of Computational Physics_ , 228(9):3478–3489, 2009.  ISSN 00219991.  doi: 10.1016/j.jcp.2009.01.026.  URL http://dx.doi.org/10.1016/j.jcp.2009.01.026.
  * Ueda and Kida (2021) Y. Ueda and T. Kida.  Asymptotic analysis of initial flow around an impulsively started circular cylinder using a brinkman penalization method.  _Journal of Fluid Mechanics_ , 929:A31, 2021.  doi: 10.1017/jfm.2021.869.
  * van der Vorst (1992) H. A. van der Vorst.  Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems.  _SIAM Journal on Scientific and Statistical Computing_ , 13(2):631–644, 1992.  doi: 10.1137/0913035.
  * Van Hooft et al. (2018) J. Antoon Van Hooft, Stéphane Popinet, Chiel C. Van Heerwaarden, Steven J. A. Van Der Linden, Stephan R. De Roode, and Bas J. H. Van de Wiel.  Towards Adaptive Grids for Atmospheric Boundary-Layer Simulations.  _Boundary Layer Meteorology_ , pages 1–23, February 2018.  doi: 10.1007/s10546-018-0335-9.  URL https://hal.archives-ouvertes.fr/hal-01689036.
  * Vasilyev and Kevlahan (2005) Oleg V. Vasilyev and Nicholas K.-R. Kevlahan.  An adaptive multilevel wavelet collocation method for elliptic problems.  _Journal of Computational Physics_ , 206(2):412 – 431, 2005.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2004.12.013.  URL http://www.sciencedirect.com/science/article/pii/S0021999104005200.
  * Verma et al. (2017) Siddhartha Verma, Gabriele Abbati, Guido Novati, and Petros Koumoutsakos.  Computing the force distribution on the surface of complex, deforming geometries using vortex methods and brinkman penalization.  _International Journal for Numerical Methods in Fluids_ , 85(8):484–501, 2017.  doi: https://doi.org/10.1002/fld.4392.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.4392.
  * Wissink et al. (2001) Andrew Wissink, Richard Hornung, Scott Kohn, Steve Smith, and Noah Elliott.  Large scale parallel structured amr calculations using the samrai framework.  _Proceedings of Supercomputing 2001_ , pages 22– 22, 12 2001.  doi: 10.1109/SC.2001.10029.
  * Zhang et al. (2019) Weiqun Zhang, Ann Almgren, Vince Beckner, John Bell, Johannes Blaschke, Cy Chan, Marcus Day, Brian Friesen, Kevin Gott, Daniel Graves, Max Katz, Andrew Myers, Tan Nguyen, Andrew Nonaka, Michele Rosso, Samuel Williams, and Michael Zingale.  AMReX: a framework for block-structured adaptive mesh refinement.  _Journal of Open Source Software_ , 4(37):1370, May 2019.  doi: 10.21105/joss.01370.  URL https://doi.org/10.21105/joss.01370.


# Effects of gravity in extra dimensions in atomic phenomena

V. A. Dzuba, V. V. Flambaum, P. Munro-Laylim  School of Physics, University of
New South Wales, Sydney 2052, Australia

###### Abstract

We use the difference between theory and experiment for energy intervals in
simple atomic systems (hydrogen, muonium, positronium and deuteron) to find
limits on the size of extra space dimensions in the Arkani-Hamed - Dimopoulos
- Dvali model for gravitation potential on short distances. As an additional
experimental fact we use absence of the small size gravitational bound states
of elementary particles. We demonstrate that the perturbation theory approach
does not work and more reliable results are obtained by solving the Dirac
equations for an electron in Coulomb and gravitational fields. These results
probe smaller distances than distance between nuclei in molecules and the
limits are significantly stronger than the limits on the size of extra
dimensions obtained using spectra of hydrogen molecules.

##  I Introduction

Observations indicate that our universe has three non-compact spatial
dimensions, with all forces and particles operating inside these dimensions.
However, there are popular theoretical models with extra spatial dimensions
such as string theories - see e.g. Witten (1995). These models have been
motivated by search for theory unifying all interactions and producing finite
results which do not require hiding infinities using the renormalisation
procedures.

In this paper, we examine the Arkani-Hamed - Dimopoulos - Dvali (ADD) model
Arkani–Hamed et al. (1998) which aimed at solving the hierarchy problem by
proposing that gravity can propagate through $n$ extra spatial compact
dimensions as well as the regular three non-compact dimensions. Ordinary
Newtonian gravitational interaction between elementary particles is many
orders of magnitude smaller than other interactions. In ADD model the observed
Newton gravitational law in our three dimensions will be significantly
strengthened at distance smaller than the size of extra dimensions $R$. The
Newtonian gravitational potential will change from $1/r$ to a more singular
$1/r^{n+1}$ dependence at distances smaller than the size of the extra
dimensions, $R$. This may be easily explained by the Gauss integrated flux
formula for the gravitational force since size of the ”surface” in $(2+n)D$ is
proportional to $r^{2+n}$ Arkani–Hamed et al. (1998) (see also the Randall-
Sundrum multi-dimensional models Randall and Sundrum (1999a, b) which have
similar small-distance gravity). Search for macroscopic effects has not found
any deviation from Newton law. Gravity force has been observed to obey the
inverse-square law down to the $\mu$m scale Tan et al. (2020), see also Refs.
Kapner et al. (2007); Adelberger et al. (2007); Chen et al. (2016); Vasilakis
et al. (2009); Terrano et al. (2015). Smaller distances may be probed in
atomic phenomena. Several papers Feng and Hong-Ya (2006); Luo and Liu (2007);
Dahia and Lemos (2016); Zhi-Gang et al. (2008); Wang and Ni (2013) have
attempted to calculate this effect of extra dimensions on spectra of simple
atomic systems and placed constraints on $R$ through first order perturbation
theory. The problem is that the results strongly depend on an unknown (and
practically arbitrary) cut-off parameter which has to be introduced to obtain
finite energy shift produced by the gravitational potential $g/r^{n+1}$ for
$n>1$. Another problem is that perturabation theory is not applicable if the
cut-off parameter is too small and the gravitational potential near the cut-
off exceeds the Coulomb potential. These problems have been solved in Ref.
Salumbides et al. (2015) where the authors studied effects of potential
$g/r^{n+1}$ between nuclei in H2, D2 and HD+ molecules. Since nuclei in
molecules are separated by distance exceeding Bohr radius $a_{B}$, there is no
need in cut-off parameter here. However, the effects of the potential
$g/r^{n+1}$ are much bigger for subatomic distances which give the main
contribution to the energy shift in atoms.

In this paper we want to solve the cut-off problem using absence of small size
bound states between two elementary particles which could be produced by the
singular potential $g/r^{n+1}$ if the cut-off parameter is too small. Indeed,
the highly singular nature of the ADD gravitational potential introduces a
fall-to-centre problem, where must exist some new physics mechanism that cuts
off the gravitational potential below some cut-off radius $r_{c}$. To solve
the cut-off problem empirically, we will use simultaneously two experimental
facts: absence of the collapsed gravitational bound states of elementary
particles and difference between the measured and calculated (using QED)
transition energy. In this way we can get conservative limits on the size of
extra dimensions by excluding the area of unknown physics at distances smaller
than the electroweak scale.

##  II Overview of The ADD Model

The ADD model was first proposed in Ref. Arkani–Hamed et al. (1998) and
introduced $n$ extra spatial compactified dimensions of size $R$. There is
also a change to the Planck mass, where a higher dimensional Planck mass $M$
is defined with respect to $R$ and replaces the observed three-dimensional
Planck mass $M_{Pl}$,

$\left(\frac{cM_{Pl}}{\hbar}\right)^{2}=R^{n}\left(\frac{cM}{\hbar}\right)^{n+2}.$
(1)

This means that Gauss’s Law can be applied to the gravitational potential for
two masses $m_{1}$ and $m_{2}$ separated by distance $r$:

$\displaystyle V(r)=-\frac{Gm_{1}m_{2}}{r}$
$\displaystyle\mathrm{for}\,\,r>>R,$
$\displaystyle-\frac{Gm_{1}m_{2}R^{n}}{r^{n+1}}=-\frac{\hbar^{n+1}}{c^{n-1}M^{n+2}}\frac{m_{1}m_{2}}{r^{n+1}}$
$\displaystyle\mathrm{for}\,\,r<<R,$ (2)

where $G=6.674^{-11}$ m3 kg-1 s-2 is the regular gravitational constant. The
long-range component of the potential matches the regular Newtonian
gravitational potential that is known to be accurate from $\mu$m scales to
astronomical scales. However, the short-range component has a more singular
form, therefore this is the region where gravitational potential can become
significantly stronger with extra dimensions.

For atomic calculations it is convenient to present relation between Plank
mass $M$ for $4+n$ dimensions and radius of extra dimensions $R$ in the
following form:

$Mc^{2}=(1.22\cdot
10^{19})^{\frac{2}{n+2}}10^{-\frac{6n}{n+2}}\left(\frac{3.73a_{B}}{R}\right)^{\frac{n}{n+2}}\,\,{\rm
GeV}.$ (3)

We introduce the dimensionless variable $S$ to describe the ratio of the
three-dimensional gravitational and Coulomb potential

$S\equiv\frac{Gm_{1}m_{2}}{e^{2}}=\frac{m_{1}m_{2}}{\alpha M_{Pl}^{2}},$ (4)

and therefore the short-range ADD potential can be written as

$V(r)=-\frac{\hbar c\alpha SR^{n}}{r^{n+1}}.$ (5)

We can also consider the gravity-Coulomb boundary

$\frac{\hbar c\alpha SR^{n}}{r_{gc}^{n+1}}=\frac{e^{2}}{r_{gc}},$ (6)

which provides the cut-off radius $r_{gc}=S^{1/n}R$ for the applicability of
the perturbation theory for $n>2$. Indeed, for $n>2$ the matrix elements of
$V(r)=g/r^{n+1}$ are dominated by the small distance integral near the cut off
radius $r_{c}$, $\int_{b}d^{3}r/r^{n+1}\sim 1/r_{c}^{n-2}$. We can present the
correction to energy as

$\displaystyle\delta E=<\psi|H|\psi>-<\psi_{o}|H_{0}|\psi_{0}>=$
$\displaystyle<\psi_{0}|H-H_{0}|\psi_{0}>+<\psi|H|\psi>-<\psi_{0}|H|\psi_{0}>,$

where $H-H_{0}=V(r)$ is the gravitational potential. Since in the area
$r<r_{gc}$ the correction to the wave function is not small, i.e. condition
$\psi-\psi_{0}\ll\psi_{0}$ is violated for $r_{gc}>r_{c}$, the first order
perturbation theory result is incorrect even if $\delta E$ is small.

One may say that distance $r_{gc}$ may be associated with the scale where all
interactions have a comparable strength.

##  III The Perturbation Approach

As mentioned, we expect that perturbation theory for the potential
$V(r)=g/r^{n+1}$ is applicable for $n<2$ and not applicable for $n>2$. For
intermediate case $n=2$ the perturbation theory may provide an estimate with a
logarithmic accuracy. To make a link to ADD theory, we should find estimates
for the size of extra dimensions $R_{n}$, which determines area of highly
singular gravitational potential, $V(r)=g/r^{n+1}$ for $r<R_{n}$. We will do
these first estimates using experimental data from muonium, postronium and
hydrogen spectra and perturbation theory.

For $n=1$, we use the known expectation values of $1/r^{2}$ potentials for
$s$-waves Bethe and Salpeter (2014) and obtain the energy shift produced by
potential Eq. (5):

$\delta E_{1}=\frac{2\hbar c\alpha SRZ^{2}}{n_{p}^{3}a^{2}},$ (7)

where $a$ is the reduced atomic Bohr radius, $Z$ is the nuclear charge (in
units of proton electric charge $e$), and $n_{p}$ is the principal quantum
number.

For $n=2$ and $n\geq 3$, we estimate leading terms enhanced by the small cut
off parameter $r_{c}$ in the first order perturbation energy shift:

$\delta E_{2}=\frac{4\hbar c\alpha
SR^{2}Z^{3}}{n_{p}^{3}a^{3}}\ln\left(\frac{a}{r_{c}}\right),$ (8) $\delta
E_{n\geq 3}=\frac{4\hbar c\alpha
SR^{n}Z^{3}}{(n-2)n_{p}^{3}a^{3}r_{c}^{n-2}}.$ (9)

We consider two different choices for the cut-off parameter $r_{c}$. Firstly,
we consider the point $r_{gc}$ where the gravitational and Coulomb potentials
are equal. This choice corresponds to the boundary of applicability of the
perturbation theory. Defining the variable
$\lambda\equiv\frac{1}{n_{i}^{3}}-\frac{1}{n_{j}^{3}}$ as the inverse cubed
principal quantum number difference of the two energy states, we determine the
size for any extra dimensions $R_{n}$ for a given energy shift from gravity
for a transition $n_{i}\rightarrow n_{f}$, assuming cut-off radius
$r_{c}=r_{gc}=S^{1/n}R$ from equality of the gravitational and Coulomb
potentials:

$R_{1}=\frac{a^{2}\delta E}{2\hbar c\alpha SZ^{2}\lambda},$ (10)
$R_{2}=\sqrt{\frac{a^{3}\delta E}{4\hbar c\alpha
SZ^{3}\lambda\ln(a/S^{1/n}R_{2})}},$ (11) $R_{n\geq
3}=\sqrt{\frac{a^{3}(n-2)\delta E}{4\hbar c\alpha S^{2/n}Z^{3}\lambda}}.$ (12)

Note for $n=2$, $R_{2}$ must be obtained using iterations; all other $R_{n}$
can be trivially calculated. To have inputs for $\delta E$, we consider the
maximal energy shift between experimental and theoretical energy levels in
simple systems. The systems were chosen due to their small deviation between
theoretical and experimental results, which is presented in Table 1.
Spectroscopy results for hydrogen are new, while muonium and positronium
measurements are older; see Ref. Stadnik (2022); Adkins et al. (2022) for
recent reviews in muonium and positronium spectroscopy. Significant
improvements for muonium and positronium $1s-2s$ spectroscopy are expected in
the near future Crivelli (2018); Heiss et al. (2018).

Table 1:  Input parameters for simple systems. All systems have atomic charge $Z=1$ and reduced Bohr radius $a_{B}$ (with the exception of positronium with $2a_{B}$). System | $\Delta E=E_{exp}-E_{thr}$ | Max. $\Delta E$ (eV) | Ref.
---|---|---|---
H $1s-2s$ | -0.9(5.4) kHz | $2.2\times 10^{-11}$ | Ahmadi et al. (2018)
H $1s-3s$ | 1.2(4.1) kHz | $2.2\times 10^{-11}$ | Fleurbaey (2017)
Mu $1s-2s$ | 5.6(9.9) MHz | $6.4\times 10^{-8}$ | Meyer et al. (2000)
Ps $1s-2s$ | -5.78(3.5) MHz | $1.3\times 10^{-8}$ | Czarnecki et al. (1999); Manohar and Stewart (2000); Fee et al. (1993)

In our calculations, we use the $S$ values for each system: $S_{ee}=2.4\times
10^{-43}$ for positronium and $S_{e\mu}=4.9\times 10^{-41}$ for muonium. For
hydrogen stable gravitational bound states could be formed between electron
and point-like quark. However, when we calculate the gravitational energy
shift, the total proton mass density defines the strength of the gravitational
potential. For $n>2$ the quark-electron potential $V(r)=g/r^{n+1}$ with a cut-
off $r_{c}$, which is very small on atomic scale, is practically equivalent to
the contact potential $V(r)=C\delta(r)$, where $C=\int V(r)d^{3}r$. For a
finite size proton or nucleus this approximation gives potential proportional
to the mass density $m\rho(r)$, where $\int\rho(r)d^{3}r=1$. In light atoms,
s-wave electron wave function tends to constant at $r\to 0$ and finite size of
the nucleus has no effect, i.e. we may take $V(r)=C\delta(r)$ for $n>2$.

For $n=1$ the effect of the potential $V(r)=g/r^{2}$ is not sensitive to both
the cut-off radius and nuclear size. For $n=2$ there is only very weak
logarithmic sensitivity to these parameters. Thus we may conclude that
formulas presented above may be used to estimate energy shift for hydrogen
atom with $S_{ep}=4.4\times 10^{-40}$. Use of the proton mass (instead of the
quark mass) for the estimate of the cut-off parameter, overestimates $r_{c}$
and makes the limits on $R$ weaker but more reliable.

A smaller value of the cut-off distance $r_{c}$ comes from the condition of
absence of the gravitational bound states of elementary particles (since they
have not been observed experimentally). The estimate follows from the equality
of the relativistic kinetic energy $\hbar c/r$ and gravitational energy Eq.
(II) for weakly bound states. It leads to

$(r_{c}/R)^{n}\sim m_{1}m_{2}/M_{Pl}^{2},$ (13)

which gives

$r_{c}\sim\alpha r_{gc}.$ (14)

This choice leads to stronger constrains on the size of extra dimensions (see
Table 4). Note that we still use the non-relativistic expression for
gravitation interaction (II). The estimation of the role of the relativistic
effects is presented in the Appendix. Relativistic corrections increase
gravitational interaction and make constrain on the radius of extra dimensions
$R$ stronger.

###  III.1 Deuteron Binding Energy

Wave function of the deuteron may be found using the short range character of
the strong interaction and relatively small binding energy of the deuteron.
Outside the interaction range, we use solution to the Schrodinger equation for
zero potential. Within the interaction range $r_{0}=1.2$ fm, wave function has
a constant value for $s$ orbital.

$\psi(r)=\left\\{\begin{array}[]{ll}\frac{Be^{-\kappa r}}{r}&\text{ for
}r>r_{0},\\\ \frac{B\,J(0)}{r_{0}}&\text{ for }r<r_{0},\end{array}\right.$
(15)

where the normalisation constant $B$ is given by $4\pi B^{2}=2\kappa$ for
$\kappa=\sqrt{2m|E|}=4.56\times 10^{7}$ eV (reduced mass $m=m_{p}/2$ and
binding energy $|E|=2.22$ MeV). The Jastrow factor, $J(0)=0.4$ Dmitriev et al.
(1983), is included to account for the nucleon repulsion at short distance.
For $n=1$ we obtain

$\delta E_{1}=-\frac{2\kappa R_{1}\alpha S_{pn}}{r_{0}},$ (16)

where $S_{pn}=0.81\cdot 10^{-36}$, $\alpha S_{pn}=0.59\cdot 10^{-38}$.
Following Ref.Stadnik (2018) we take difference between experimental Kessler,
Jr et al. (1999) and theoretical Ekström et al. (2015) results as
$E_{\rm{exp}}-E_{\rm{thr}}=-13.7$ eV. This gives $R_{1}<3\cdot 10^{16}$ m.

For $n>1$ the quark-quark potential $V(r)=-g/r^{n+1}$ with a cut-off $r_{c}$,
which is very small on nuclear scale, may be approximated by the contact
potential $V(r)=-C\delta(r)$, where $C=-\int V(r)d^{3}r$. This contact
interaction is similar to the weak interaction mediated by $Z$-boson. From the
results of the weak effects calculation we know that finite size of the
nucleons produces suppression of the effects described by the Jastrow factor
$J(0)$ Dmitriev et al. (1983). After accounting for this factor we may assume
that nucleons are point-like particles with interaction $V(r)=-C\delta(r)$.

Using wave function in Eq. (15) we obtain the following estimate for the
energy shift in the case of $n>1$:

$\delta
E_{n}=\langle\psi|-C_{n}\delta(\mathbf{r})|\psi\rangle=-\frac{C_{n}\kappa
J(0)^{2}}{2\pi r_{0}^{2}}=-C_{n}0.8\cdot 10^{36}\frac{eV}{m^{2}},$ (17)

where for $n=2$

$C_{2}=4\pi\alpha S_{pn}R_{2}^{2}\ln\frac{r_{0}}{r_{c}},$ (18)

for $n>2$

$C_{n}=\frac{4\pi\alpha
S_{pn}R_{n}^{2}}{(n-2)}\left(\frac{R_{n}}{r_{c}}\right)^{n-2}.$ (19)

We can estimate the cut-off radius $r_{c}$ from the condition of the absence
of the small size gravitational bound states Eq. (13). Note that gravitational
bound state may exist for two quarks while all quarks contribute to the
deuteron binding energy. To have a conservative estimate of the limits on $R$
and $M$, we use a constituent quark mass $m_{q}=300$ MeV in Eq. (13) (mass
$m_{q}\sim$ 5 MeV would give smaller $r_{c}$, bigger energy shifts and
stronger limits). Using $E_{\rm{exp}}-E_{\rm{thr}}=-13.7$ eV we obtain
estimates for $R_{n}$ and $M_{n}$ presented in Table 2. Formally, these
estimates look like the strongest constrains among two-body systems. However,
deuteron is a system with the strong interaction and these constrains are
probably less reliable than the constrains from hydrogen, muonium and
positronium.

Table 2:  Parameters of extra dimensions obtained from the deuteron data. $n$ is the number of extra dimensions, $R$ is the size of extra dimensions, $r_{c}$ is the cut-off parameter, $M$ is Planck mass for extra dimensions (3). Numbers in square brackets mean powers of ten. $n$ | $R$ | $r_{c}$ | $M$
---|---|---|---
| [m] | [m] | [GeV]
2 | 5.1 | 4.5[-19] | 22
3 | 4.8[-06] | 9.6[-19] | 25
4 | 3.8[-09] | 1.1[-18] | 32
5 | 5.2[-11] | 1.2[-18] | 38
6 | 3.0[-12] | 1.3[-18] | 43
7 | 4.0]-13] | 1.4[-18] | 47

##  IV Numerical calculations

Table 3:  Parameters of extra dimensions obtained from comparing experimental and theoretical energy shifts. Experimental shifts $\Delta E_{\rm expt}$ are taken from Table 1. Parameter $n$ is the number of extra dimensions; $\Delta E_{a}$ is the energy shift given by Eqs. (7,8,9); $\Delta E_{\rm PT0}=\langle\psi^{(0)}_{1s}|V|\psi^{(0)}_{ns}\rangle$, where $\psi^{(0)}$ are unperturbed wave functions; $\Delta E_{\rm PT}=\langle\psi_{1s}|V|\psi_{ns}\rangle$, where $\psi$ are wave functions obtained by solving Dirac equations including Coulomb and gravitation potentials. Energy shifts $\Delta E_{D}$, obtained using Dirac equations, are fitted to $\Delta E_{\rm expt}$ by varying size of extra dimensions $R$, i.e. $\Delta E_{D}=\Delta E_{\rm expt}$ . The values of $R_{a}$ are given by Eqs. (10,11,12), $R_{D}$ is chosen to fit experimental energy shift by solving Dirac equations; $r_{gca}=S^{1/n}R_{a}$, $r_{gc}=S^{1/n}R_{D}$. $M$ is Planck mass for extra dimensions (3). The calculation were done at the cut-off parameter $r_{c}=r_{gc}$. Numbers in square brackets mean powers of ten. $n$ | $\Delta E_{a}$ | $\Delta E_{\rm PT0}$ | $\Delta E_{\rm PT}$ | $R_{a}$ | $r_{gca}$ | $R_{D}$ | $r_{gc}$ | $M$
---|---|---|---|---|---|---|---|---
| [eV] | [eV] | [eV] | [m] | [m] | [m] | [m] | [GeV]
Hydrogen, $1s-2s$, $\Delta E_{D}=\Delta E_{\rm expt}=2.23\times 10^{-11}$ eV
1 | 2.23[-11] | 2.23[-11] | 2.23[-11] | 5.63[+16] | 2.48[-23] | 5.63[+16] | 2.48[-23] | 80.5
2 | 2.41[-11] | 2.27[-11] | 2.27[-11] | 3.06[+02] | 6.42[-18] | 3.18[+02] | 6.68[-18] | 2.75
3 | 1.69[-11] | 2.26[-11] | 2.26[-11] | 3.37[-04] | 2.56[-17] | 2.93[-04] | 2.23[-17] | 2.15
4 | 1.29[-11] | 2.15[-11] | 2.15[-11] | 2.50[-07] | 3.62[-17] | 1.90[-07] | 2.75[-17] | 2.36
5 | 1.09[-11] | 2.19[-11] | 2.19[-11] | 3.30[-09] | 4.44[-17] | 2.31[-09] | 3.11[-17] | 2.54
6 | 9.44[-12] | 2.20[-11] | 2.20[-11] | 1.86[-10] | 5.13[-17] | 1.21[-10] | 3.33[-17] | 2.70
7 | 8.04[-12] | 2.14[-11] | 2.14[-11] | 2.40[-11] | 5.73[-17] | 1.44[-11] | 3.44[-17] | 2.88
Hydrogen, $1s-3s$, $\Delta E_{D}=\Delta E_{\rm expt}=2.19\times 10^{-11}$ eV
1 | 2.22[-11] | 2.22[-11] | 2.22[-11] | 5.53[+16] | 2.44[-23] | 5.08[+16] | 2.24[-23] | 83.3
2 | 2.41[-11] | 2.27[-11] | 2.27[-11] | 3.03[+02] | 6.36[-18] | 3.03[+02] | 6.36[-18] | 2.82
3 | 1.62[-11] | 2.16[-11] | 2.16[-11] | 3.34[-04] | 2.54[-17] | 2.74[-04] | 2.08[-17] | 2.24
4 | 1.32[-11] | 2.20[-11] | 2.20[-11] | 2.48[-07] | 3.59[-17] | 1.83[-07] | 2.66[-17] | 2.42
5 | 1.12[-11] | 2.23[-11] | 2.23[-11] | 3.27[-09] | 4.40[-17] | 2.22[-09] | 2.99[-17] | 2.61
6 | 9.58[-12] | 2.24[-11] | 2.24[-11] | 1.84[-10] | 5.08[-17] | 1.16[-10] | 3.20[-17] | 2.79
7 | 8.54[-12] | 2.28[-11] | 2.28[-11] | 2.38[-11] | 5.68[-17] | 1.42[-11] | 3.38[-17] | 2.92
Muonium, $\Delta E_{D}=\Delta E_{\rm expt}=6.41\times 10^{-8}$ eV
1 | 6.41[-08] | 6.41[-08] | 6.41[-08] | 1.43[+21] | 7.12[-20] | 1.43[+21] | 7.12[-20] | 2.74
2 | 7.04[-08] | 6.48[-08] | 6.48[-08] | 5.67[+04] | 4.00[-16] | 5.95[+04] | 4.20[-16] | 0.201
3 | 4.74[-08] | 6.32[-08] | 6.32[-08] | 3.73[-02] | 1.37[-15] | 3.21[-02] | 1.18[-15] | 0.128
4 | 3.90[-08] | 6.50[-08] | 6.50[-08] | 2.31[-05] | 1.94[-15] | 1.80[-05] | 1.51[-15] | 0.114
5 | 3.32[-08] | 6.65[-08] | 6.65[-08] | 2.73[-07] | 2.38[-15] | 1.97[-07] | 1.71[-15] | 0.106
6 | 2.71[-08] | 6.32[-08] | 6.32[-08] | 1.43[-08] | 2.75[-15] | 9.31[-09] | 1.78[-15] | 0.104
7 | 2.54[-08] | 6.78[-08] | 6.79[-08] | 1.76[-09] | 3.07[-15] | 1.11[-09] | 1.93[-15] | 0.098
Ps, $1s-2s$, $\Delta E_{D}=\Delta E_{\rm expt}=1.34\times 10^{-8}$ eV
1 | 1.34[-8] | 1.34[-8] | 1.34[-8] | 2.49[+23] | 5.97[-20] | 2.49[+23] | 5.97[-20] | 4.9[-1]
2 | 1.50[-8] | 1.26[-8] | 1.26[-8] | 1.04[+06] | 5.08[-16] | 1.10[+06] | 5.39[-16] | 4.7[-2]
3 | 9.71[-9] | 9.71[-9] | 9.71[-9] | 2.86[-01] | 1.78[-15] | 2.43[-01] | 1.51[-15] | 3.8-2]
4 | 7.97[-9] | 7.97[-9] | 7.97[-9] | 1.14[-04] | 2.51[-15] | 8.75[-05] | 1.94[-15] | 4.0[-2]
5 | 6.58[-9] | 6.59[-9] | 6.59[-9] | 1.03[-06] | 3.08[-15] | 7.20[-07] | 2.16[-15] | 4.2[-2]
6 | 5.68[-9] | 5.68[-9] | 5.68[-9] | 4.51[-08] | 3.56[-15] | 2.93[-08] | 2.31[-15] | 4.4[-2]
7 | 5.16[-9] | 5.17[-9] | 5.17[-9] | 4.87[-09] | 3.97[-15] | 3.02[-09] | 2.46[-15] | 4.5[-2]
Table 4:  The same as Table 3 but the calculations were done at the smaller values of the cut-off parameter $r_{c}=r_{gca}/100$, chosen from the condition of the absence of the gravitational bound states of elementary particles.  $n$ | $\Delta E_{a}$ | $\Delta E_{\rm PT0}$ | $\Delta E_{\rm PT}$ | $R_{a}$ | $r_{c}$ | $R_{D}$ | $M$
---|---|---|---|---|---|---|---
| [eV] | [eV] | [eV] | [m] | [m] | [m] | [GeV]
Hydrogen, $1s-2s$, $\Delta E_{D}=\Delta E_{\rm expt}=2.23\times 10^{-11}$ eV
2 | 1.01[-12] | 9.49[-13] | 1.26[-08] | 3.06[+02] | 6.42[-20] | 5.73[+01] | 6.5
3 | 7.58[-13] | 7.59[-13] | 8.41[-10] | 3.37[-04] | 2.56[-19] | 2.35[-05] | 9.8
4 | 7.85[-13] | 7.86[-13] | 3.23[-10] | 2.50[-07] | 3.62[-19] | 1.08[-08] | 16
5 | 8.04[-13] | 8.05[-13] | 1.96[-10] | 3.30[-09] | 4.44[-19] | 1.07[-10] | 23
6 | 8.17[-13] | 8.17[-13] | 1.36[-10] | 1.86[-10] | 5.13[-19] | 4.97[-12] | 30
7 | 8.25[-13] | 8.26[-13] | 9.73[-11] | 2.40[-11] | 5.73[-19] | 5.59[-13] | 36
Hydrogen, $1s-3s$, $\Delta E_{D}=\Delta E_{\rm expt}=2.19\times 10^{-11}$ eV
2 | 1.10[-12] | 1.03[-12] | 1.33[-08] | 3.03[+02] | 6.36[-20] | 5.69[+01] | 6.5
3 | 8.18[-13] | 8.19[-13] | 7.21[-10] | 3.34[-04] | 2.54[-19] | 2.33[-05] | 9.8
4 | 8.48[-13] | 8.49[-13] | 3.49[-10] | 2.48[-07] | 3.59[-19] | 1.07[-08] | 16
5 | 8.61[-13] | 8.62[-13] | 1.09[-10] | 3.27[-09] | 4.40[-19] | 1.06[-10] | 23
6 | 8.74[-13] | 8.75[-13] | 7.92[-11] | 1.84[-10] | 5.08[-19] | 4.91[-12] | 30
7 | 8.85[-13] | 8.86[-13] | 7.33[-11] | 2.38[-11] | 5.68[-19] | 5.53[-13] | 36
Muonium, $\Delta E_{D}=\Delta E_{\rm expt}=6.41\times 10^{-8}$ eV
2 | 2.84[-09] | 2.62[-09] | 2.24[-05] | 5.67[+04] | 4.00[-18] | 1.01[+04] | 0.49
3 | 2.17[-09] | 2.17[-09] | 1.55[-06] | 3.73[-02] | 1.37[-17] | 2.60[-03] | 0.58
4 | 2.25[-09] | 2.25[-09] | 9.50[-07] | 2.31[-05] | 1.94[-17] | 1.00[-06] | 0.78
5 | 2.32[-09] | 2.33[-09] | 1.60[-06] | 2.73[-07] | 2.38[-17] | 8.89[-09] | 0.97
6 | 2.35[-09] | 2.35[-09] | 4.77[-07] | 1.43[-08] | 2.75[-17] | 3.83[-10] | 1.14
7 | 2.37[-09] | 2.37[-09] | 2.82[-07] | 1.76[-09] | 3.07[-17] | 4.09[-11] | 1.28
Ps, $1s-2s$,$\Delta E_{D}=\Delta E_{\rm expt}=1.34\times 10^{-8}$ eV
2 | 5.89[-10] | 5.21[-10] | 6.02[-06] | 1.04[+06] | 5.08[-18] | 1.85[+05] | 0.11
3 | 4.55[-10] | 4.55[-10] | 3.27[-07] | 2.86[-01] | 1.78[-17] | 1.99[-02] | 0.17
4 | 4.72[-10] | 4.73[-10] | 2.00[-07] | 1.14[-04] | 2.51[-17] | 4.92[-06] | 0.27
5 | 4.81[-10] | 4.82[-10] | 7.83[-08] | 1.03[-06] | 3.08[-17] | 3.34[-08] | 0.38
6 | 4.91[-10] | 4.92[-10] | 8.26[-08] | 4.51[-08] | 3.56[-17] | 1.21[-09] | 0.48
7 | 4.93[-10] | 4.94[-10] | 4.12[-08] | 4.87[-09] | 3.97[-17] | 1.13[-10] | 0.58

We performed the numerical calculations of the energy shifts caused by
gravitation potential in order to find limits on the size of extra dimensions
$R$. We do this in two ways. First, we use perturbation theory (PT) approach
described in previous section. Second, we solve Dirac equations for the $1s$,
$2s$, $3s$ states of hydrogen, muonium and positronium. The energy shift is
found as a difference of energies calculated in the Coulomb potential and in
the potential with added gravitational contribution. The size of extra
dimensions $R$ is used in this approach as a fitting parameter to fit input
energy shifts from Table 1 (we call these shifts ”experimental” shifts and use
notation $\Delta E_{\rm expt}$ for them). The calculations with Dirac
equations are done to check whether PT really works. Note that small value of
the correction to the energy does not necessary mean that PT is applicable.
The correction is small because it comes from a very small region in the
vicinity of the cut-off radius. However, in this region the change of the
potential is not small leading to large change in the wave function and
breaking of PT.

We perform the calculations using two ways of defining the cut-off parameter
$r_{c}$ ($V_{g}(r)=V_{g}(r_{c})$ for $r<r_{c}$). Using $r_{c}=r_{gc}$ (see Eq.
(6)) leads to conservative estimates of the size of extra dimensions which are
in agreement with the perturbation theory calculations. Corresponding results
are presented in Table 3. In another approach we use much smaller value of
$r_{c}$ which we find from the condition of the absence of the small size
gravitational bound states of elementary particles. In practice this means
that the wave function, found from solving the Dirac equation, does not
oscillate at small distances. We found that in all cases $r_{c}=r_{gca}/100$,
where $r_{gca}=S^{1/n}R_{a}$, in excellent agreement with an estimate
$r_{c}\sim\alpha r_{gc}$ from Eqs. (13,14). The results are presented in Table
4.

In Tables 3 and 4 we present energy shifts and corresponding values of
$r_{c}$, $R$ and $M$ obtained in several different ways. $\Delta E_{a}$ is
calculated using formulae (7,8,9) but with the radius of extra dimensions
found from the fitting of the energy shift $\Delta E_{D}$ obtained by solving
the Dirac equations. $\Delta E_{\rm
PT0}=\langle\psi^{(0)}_{1s}|V|\psi^{(0)}_{ns}\rangle$, where $\psi^{(0)}$ are
unperturbed wave functions.

To avoid a misunderstanding, we should note that we used slightly different
method of cut-off in the analytical approach and in the numerical solution of
the Dirac equation. $\Delta E_{\rm PT0}=\Delta E_{a}$ if integration for
$\Delta E_{\rm PT0}$ is done from $r_{c}$ to infinity, as in the analytical
approach. We did such calculations as a test and found good agreement between
the two. However, in the tables we present values of $\Delta E_{\rm PT0}$ and
$\Delta E_{\rm PT}$ found by integrating from zero to infinity. This is done
to make meaningful comparison with the energy shift obtained from solving the
Dirac equations. The same gravitation potential is used to calculate $\Delta
E_{\rm PT0}$, $\Delta E_{\rm PT}$ and $\Delta E_{D}$ in the Dirac equations.
If $\Delta E_{PT0}=\Delta E_{PT}=\Delta E_{\rm expt}$ ( we have $\Delta
E_{D}=\Delta E_{\rm expt}$) then one could say that PT works. We see that this
is the case only for $r_{c}=r_{gca}$ (Table 3).

The values of $R_{a}$ and $r_{gca}$ are given by Eqs. (10,11,12) and (6). The
values of $R_{D}$ and $r_{gc}$ are found from solving the Dirac equations. We
believe that the results obtained with the use of Dirac equations are more
reliable. The Planck mass for extra dimensions $M$ was found using the Dirac
equations value of $R$ ($R_{D}$).

##  V Conclusion

In the present paper we investigated possibilities to probe ADD gravitational
potential Eq. (II) at sub-atomic distances. This potential does not include
relativistic corrections. However, as it is demonstrated in Appendix,
relativistic corrections increase gravitational interaction and lead to much
stronger constrains. Therefore, our estimates of the limits on the size of
extra dimensions, based on the potential Eq. (II), are conservative. By
imposing the condition that there are no gravitational bound states between
stable elementary particles (which have not been observed), we cut-off area of
of unknown physics at distances smaller than the electroweak scale. Constrains
on the size of extra dimensions $R$ and Plank mass $M$, obtained using the
difference between the experimental and calculated value of the deuteron
binding energy, have been presented in the Table 2. However, deuteron is a
system with the strong interaction and these constrains are probably less
reliable than the constrains from hydrogen, muonium and positronium presented
in the Table 4. The perturbation theory is not applicable if the main
contribution to the energy shift comes from the area where the gravitational
potential exceeds the Coulomb potential. Therefore, we obtain the energy shift
using solution of the Dirac equation including both Coulomb and gravitational
potential. Corresponding limit on the size of extra dimensions is denoted by
$R_{D}$.

## Acknowledgements

This work was supported by the Australian Research Council Grants No.
DP190100974 and DP200100150.

##  VI Appendix: relativistic estimates

As known, all types of energy contribute to the gravitational potential. For
example, kinetic energy of nearly massless quarks ($m_{u}\approx$ 3 MeV,
$m_{d}\approx$ 5 MeV) give significant contribution to the nucleon mass.
Kinetic energy of ultra-relativistic particle may be estimated as
$E_{k}=pc\sim\hbar c/r$, so mass $m$ in the gravitational potential should be
replaced by $\hbar/cr$, and corresponding gravitational interaction energy of
two particles is

$U\sim-\frac{G}{r}\left(\frac{\hbar}{cr}\right)^{2}\left(\frac{R_{n}}{r}\right)^{n}$
(20)

For weakly bound states kinetic energy and potential energy compensate each
other, $E_{k}=|U|$, and this equality defines the cut-off radius excluding
formation of the bound states, $r_{c}=\hbar/M_{n}c$, where $M_{n}$ is the
Plank mass for extra $n$ dimensions -see Eq. (1). Note that one can find
similar estimates in Ref. Guenther et al. (2000).

Perturbation theory gives the following estimate for the energy shift produced
by the potential Eq. (20) with the cut-off $r_{c}=\hbar/M_{n}c$ (in natural
units $\hbar=c=1$):

$\delta E\sim\frac{4\pi\psi(0)^{2}}{nM_{n}^{2}}$ (21)

For hydrogen $\delta E<2.2\cdot 10^{-11}eV$ and we obtain the limit $M_{n}>$
100 GeV $/n^{1/2}$, which is much stronger than the limit obtain using the
non-relativistic gravitational potential. For deuterium the limit is even
stronger, $M_{n}>$ 170 GeV $/n^{1/2}$.

Thus, we conclude that use of the non-relativistic gravitational potential
underestimates the energy shift $\delta E$ and gives conservative limits on
the Plank mass $M_{n}$ and radius of extra dimensions $R_{n}$.

## References

  * Witten (1995) E. Witten, Nuclear Physics B 443, 85 (1995), URL https://doi.org/10.1016/0550-3213(95)00158-O.
  * Arkani–Hamed et al. (1998) N. Arkani–Hamed, S. Dimopoulos, and G. Dvali, Physics Letters B 429, 263–272 (1998), ISSN 0370-2693, URL http://dx.doi.org/10.1016/S0370-2693(98)00466-3.
  * Randall and Sundrum (1999a) L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 3370 (1999a), URL https://link.aps.org/doi/10.1103/PhysRevLett.83.3370.
  * Randall and Sundrum (1999b) L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 4690 (1999b), URL https://link.aps.org/doi/10.1103/PhysRevLett.83.4690.
  * Tan et al. (2020) W.-H. Tan, A.-B. Du, W.-C. Dong, S.-Q. Yang, C.-G. Shao, S.-G. Guan, Q.-L. Wang, B.-F. Zhan, P.-S. Luo, L.-C. Tu, et al., Phys. Rev. Lett. 124, 051301 (2020), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.051301.
  * Kapner et al. (2007) D. J. Kapner, T. S. Cook, E. G. Adelberger, J. H. Gundlach, B. R. Heckel, C. D. Hoyle, and H. E. Swanson, Phys. Rev. Lett. 98, 021101 (2007).
  * Adelberger et al. (2007) E. G. Adelberger, B. R. Heckel, S. Hoedl, C. D. Hoyle, D. J. Kapner, and A. Upadhye, Phys. Rev. Lett. 98, 131104 (2007).
  * Chen et al. (2016) Y.-J. Chen, W. K. Tham, D. E. Krause, D. Lopez, E. Fischbach, and R. S. Decca, Phys. Rev. Lett. 116, 221102 (2016).
  * Vasilakis et al. (2009) G. Vasilakis, J. M. Brown, T. W. Kornack, and M. V. Romalis, Phys. Rev. Lett. 103, 261801 (2009).
  * Terrano et al. (2015) W. A. Terrano, E. G. Adelberger, J. G. Lee, and B. R. Heckel, Phys. Rev. Lett. 115, 201801 (2015).
  * Feng and Hong-Ya (2006) L. Feng and L. Hong-Ya, Chinese Physics Letters 23, 2903 (2006), URL https://doi.org/10.1088/0256-307x/23/11/006.
  * Luo and Liu (2007) F. Luo and H. Liu, Int. J. Theor. Phys. 46, 606 (2007), eprint gr-qc/0602093.
  * Dahia and Lemos (2016) F. Dahia and A. S. Lemos, Phys. Rev. D 94, 084033 (2016), URL https://link.aps.org/doi/10.1103/PhysRevD.94.084033.
  * Zhi-Gang et al. (2008) L. Zhi-Gang, N. Wei-Tou, and A. P. Patón, Chinese Physics B 17, 70 (2008), URL https://doi.org/10.1088/1674-1056/17/1/013.
  * Wang and Ni (2013) L.-B. Wang and W.-T. Ni, Modern Physics Letters A 28, 1350094 (2013), ISSN 1793-6632, URL http://dx.doi.org/10.1142/S0217732313500946.
  * Salumbides et al. (2015) E. J. Salumbides, A. N. Schellekens, B. Gato-Rivera, and W. Ubachs, New Journal of Physics 17, 033015 (2015), URL https://doi.org/10.1088/1367-2630/17/3/033015.
  * Bethe and Salpeter (2014) H. A. Bethe and E. E. Salpeter, _Quantum mechanics of one- and two-electron atoms_ (Springer Berlin, 2014).
  * Stadnik (2022) Y. V. Stadnik, _Searching for ultralight scalar dark matter with muonium and muonic atoms_ (2022), URL https://arxiv.org/abs/2206.10808.
  * Adkins et al. (2022) G. Adkins, D. Cassidy, and J. Pérez-Ríos, Physics Reports 975, 1 (2022), ISSN 0370-1573, precision spectroscopy of positronium: testing bound-state QED theory and the search for physics beyond the Standard Model, URL https://www.sciencedirect.com/science/article/pii/S0370157322001508.
  * Crivelli (2018) P. Crivelli, Hyperfine Interactions 239 (2018), ISSN 1572-9540, URL http://dx.doi.org/10.1007/s10751-018-1525-z.
  * Heiss et al. (2018) M. Heiss, G. Wichmann, A. Rubbia, and P. Crivelli, _The positronium hyperfine structure: Progress towards a direct measurement of the $\text{2}^{\text{3}}\text{S}_{\text{1}}\rightarrow\text{2}^{\text{1}}\text{S}_{\text{0}}$ transition in vacuum_ (2018), URL https://arxiv.org/abs/1805.05886.
  * Ahmadi et al. (2018) M. Ahmadi, B. Alves, C. Baker, and et al., Nature 557, 71 (2018).
  * Fleurbaey (2017) H. Fleurbaey, Theses, Université Pierre et Marie Curie (UPMC) (2017), URL https://tel.archives-ouvertes.fr/tel-01633631.
  * Meyer et al. (2000) V. Meyer, S. N. Bagayev, P. E. G. Baird, P. Bakule, M. G. Boshier, A. Breitrück, S. L. Cornish, S. Dychkov, G. H. Eaton, A. Grossmann, et al., Physical Review Letters 84, 1136 (2000), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.84.1136.
  * Czarnecki et al. (1999) A. Czarnecki, K. Melnikov, and A. Yelkhovsky, Phys. Rev. A 59, 4316 (1999), URL https://link.aps.org/doi/10.1103/PhysRevA.59.4316.
  * Manohar and Stewart (2000) A. V. Manohar and I. W. Stewart, Phys. Rev. Lett. 85, 2248 (2000), URL https://link.aps.org/doi/10.1103/PhysRevLett.85.2248.
  * Fee et al. (1993) M. S. Fee, A. P. Mills, S. Chu, E. D. Shaw, K. Danzmann, R. J. Chichester, and D. M. Zuckerman, Phys. Rev. Lett. 70, 1397 (1993), URL https://link.aps.org/doi/10.1103/PhysRevLett.70.1397.
  * Dmitriev et al. (1983) V. F. Dmitriev, V. V. Flambaum, O. P. Sushkov, and V. B. Telitsin, Physics Letters B 125, 1 (1983).
  * Stadnik (2018) Y. V. Stadnik, Phys. Rev. Lett. 120, 223202 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.223202.
  * Kessler, Jr et al. (1999) E. Kessler, Jr, M. Dewey, R. Deslattes, A. Henins, H. Börner, M. Jentschel, C. Doll, and H. Lehmann, Physics Letters A 255, 221 (1999), ISSN 0375-9601, URL https://www.sciencedirect.com/science/article/pii/S037596019900078X.
  * Ekström et al. (2015) A. Ekström, G. R. Jansen, K. A. Wendt, G. Hagen, T. Papenbrock, B. D. Carlsson, C. Forssén, M. Hjorth-Jensen, P. Navrátil, and W. Nazarewicz, Phys. Rev. C 91, 051301(R) (2015), URL https://link.aps.org/doi/10.1103/PhysRevC.91.051301.
  * Guenther et al. (2000) U. Guenther, A. Kuklov, and A. Zhuk, _Comment on stability of atoms and nuclei in multidimensional tev scale gravity_ (2000), URL https://arxiv.org/abs/hep-ph/0012246.