# Probing a Dark Sector with Collider Physics, Direct Detection, and
Gravitational Waves

Giorgio Arcadi<EMAIL_ADDRESS>Dipartimento di Scienze Matematiche e
Informatiche, Scienze Fisiche e Scienze della Terra, Universita degli Studi di
Messina, Via Ferdinando Stagno d’Alcontres 31, I-98166 Messina, Italy  INFN
Sezione di Catania, Via Santa Sofia 64, I-95123 Catania, Italy  Glauber C.
Dorsch<EMAIL_ADDRESS>Departamento de Física, Universidade Federal de
Minas Gerais, 31270-901, Belo Horizonte, MG, Brasil  Jacinto P. Neto
<EMAIL_ADDRESS>Departamento de Física, Universidade Federal do
Rio Grande do Norte, 59078-970, Natal, RN, Brasil  International Institute of
Physics, Universidade Federal do Rio Grande do Norte, 59078-970, Natal, RN,
Brasil  Farinaldo S. Queiroz<EMAIL_ADDRESS>Departamento de
Física, Universidade Federal do Rio Grande do Norte, 59078-970, Natal, RN,
Brasil  International Institute of Physics, Universidade Federal do Rio Grande
do Norte, 59078-970, Natal, RN, Brasil  Millennium Institute for Subatomic
Physics at High-Energy Frontier (SAPHIR), Fernandez Concha 700, Santiago,
Chile  Y. M. Oviedo-Torres<EMAIL_ADDRESS>International Institute of
Physics, Universidade Federal do Rio Grande do Norte, 59078-970, Natal, RN,
Brasil  Departamento de Física, Universidade Federal da Paraíba, 58051-970,
João Pessoa, PB, Brasil

###### Abstract

We assess the complementarity between colliders, direct detection searches,
and gravitational wave interferometry in probing a scenario of dark matter in
the early universe. The model under consideration contains a $B-L$ gauge
symmetry and a vector-like fermion which acts as the dark matter candidate.
The fermion induces significant a large dark matter-nucleon scattering rate,
and the $Z^{\prime}$ field produces clear dilepton events at colliders. Thus,
direct detection experiments and colliders severely constrain the parameter
space in which the correct relic density is found in agreement with the data.
Nevertheless, little is known about the new scalar responsible for breaking
the B-L symmetry. If this breaking occurs via a first-order phase transition
at a TeV scale, it could lead to gravitational waves in the mHz frequency
range detectable by LISA, DECIGO, and BBO instruments. The spectrum is highly
sensitive to properties of the scalar sector and gauge coupling. We show that
a possible GW detection, together with information from colliders and direct
detection experiments, can simultaneously pinpoint the scalar self-coupling,
and narrow down the dark matter mass where a thermal relic is viable.

##  I  Introduction

The presence of dark matter (DM) in the universe is compelling, and among the
possible explanations for it, thermal relics stand out. They typically
experience interactions that could be probed by current near-future
experiments and easily yield the correct relic density, thus attracting much
attention in the community. The nature of these particles allowed us to
explore the so-called dark matter complementarity, which refers to the use of
data from various sources such as direct and indirect detection experiments,
as well as colliders, as a way to narrow down the viable properties of a dark
matter particle in a given model (See Arcadi  _et al._ (2018) for a review).
With the detection of gravitation waves in 2015 Abbott  _et al._ (2016), a new
era surfaced. No property of the dark matter particle is imprinted on the
gravitation wave signal. Although, if the scalar in the dark sector was
responsible for a first-order phase transition that would lead to the
detection of a gravitation wave, we could indirectly constrain the dark
sector. In other words, gravitation waves have become an interesting
laboratory for dark sectors that feature a scalar particle inducing a first-
order phase transition. This is precisely what we exploit here.

The detection of gravitational waves (GWs) by LIGO/Virgo/KAGRA Abbott  _et
al._ (2019, 2021a, 2021b), and the near-future launch of new-generation space-
based interferometers Amaro-Seoane  _et al._ (2017), open up other detection
venues that could be relevant for probing these scenarios, which are often
complementary to collider searches and direct detection methods. Indeed,
because our dark matter particle is associated to physics at the electroweak
scale, a first-order cosmological phase transition predicted by such model
would typically result in a GW spectrum peaked at mHz frequencies, which is
precisely the range of optimal detectability at LISA/DECIGO/BBO Crowder and
Cornish (2005); Kawamura  _et al._ (2006).

In this work, we investigate the complementarity between GW detectors,
collider searches and direct detection experiments for probing the
phenomenology of a model with a gauged $B-L$ symmetry. The model contains an
additional gauge boson $Z^{\prime}$, together with a DM candidate which is a
Dirac vector-like fermion coupled to the visible sector via a $Z^{\prime}$
portal. By computing the spin-independent nucleon-DM scattering cross-section,
we can impose constraints on the model due to direct detection limits from
XENONnT and LUX-ZEPLIN (LZ) Aprile  _et al._ (2023); Aalbers  _et al._ (2022).
Moreover, collider searches by ATLAS and LEP II impose bounds on the mass and
coupling of the $Z^{\prime}$ Aad  _et al._ (2019, 2020); Aaboud  _et al._
(2019). But these experiments can say little about the scalar sector
responsible for breaking the $U(1)_{B-L}$ symmetry: the scalar is not relevant
for the DM phenomenology and may decay only into invisible particles, making
it difficult to detect at colliders. However, the $U(1)_{B-L}$ symmetry
breaking process may be a first-order phase transition in the early Universe,
resulting in GWs testable at the aforementioned detectors, as shown in Ref.
Hasegawa  _et al._ (2019); Bosch  _et al._ (2023). Here, we show that the GW
spectrum is highly sensitive to the scalar sector and could be used to
constrain it experimentally. In fact, a detection of cosmological GWs, allied
to collider/direct detection measurements on the gauge and DM sector, could
even lead to a measurement of the scalar self-coupling at a precision of $\sim
50\%$ using GWs. This illustrates how a GW detection could be complementary to
collider and direct detection searches as probes for dark sector and beyond
the Standard Model (SM) physics in general.

This work is organized as follows: In Section II, we present the minimal $B-L$
model describing the particle content and the properties of the scalar singlet
and DM candidate that we are interested in. The generation of GWs via first-
order phase transition in the context of the minimal $B-L$ model is discussed
in Section III, and the DM production in Section IV. We explain the DM model
constraints in Section V. Finally, we discuss the results and present the
conclusions in Sections VI and VIII, respectively.

##  II  The Minimal B-L Model

The minimal $B-L$ model is a simple extension of the Standard Model (SM). An
additional $U(1)_{B-L}$ Abelian symmetry enlarges the SM gauge structure to
$SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}\otimes U(1)_{B-L}$, where $B$
stands for baryon number and $L$ for lepton number. Hence, the gauge content
increases by one new gauge boson, say, $Z^{\prime}$. We highlight that we are
assuming a negligible kinetic mixing. In the fermion sector, we add three
Majorana right-handed neutrinos $N_{iR}$ ($i=1,2,3)$ to cancel out the gauge
anomalies, and the DM candidate is a vector-like Dirac fermion $\chi$, which
also carries $B-L$ charge as the other fermions. Consequently, the DM-SM
interactions happen through a $Z^{\prime}$ portal. Moreover, we include a
scalar singlet $\Phi_{s}$, which breaks the $B-L$ symmetry and gives rise to a
Majorana mass term for right-handed neutrinos that is essential to realize the
seesaw mechanism type I to turn active neutrinos to massive particles. In
Table 1, we present the matter particle content and their respective charge
assignments under each symmetry, including the $Z_{2}$, which arises after
$B-L$ breaking and ensures DM stability Klasen  _et al._ (2017).

| $SU(3)_{C}$ | $SU(2)_{L}$ | $U(1)_{Y}$ | $U(1)_{B-L}$ | $Z_{2}$
---|---|---|---|---|---
$q_{iL}$ | $\boldsymbol{3}$ | $\boldsymbol{2}$ | $1/6$ | $1/3$ | $+$
$u_{iR}$ | $\boldsymbol{3}$ | $\boldsymbol{1}$ | $2/3$ | $1/3$ | $+$
$d_{iR}$ | $\boldsymbol{3}$ | $\boldsymbol{1}$ | $-1/3$ | $1/3$ | $+$
$\ell_{iL}$ | $\boldsymbol{1}$ | $\boldsymbol{2}$ | $-1/2$ | $-1$ | $+$
$e_{iR}$ | $\boldsymbol{1}$ | $\boldsymbol{1}$ | $-1$ | $-1$ | $+$
$N_{iR}$ | $\boldsymbol{1}$ | $\boldsymbol{1}$ | $0$ | $-1$ | $+$
$\chi$ | $\boldsymbol{1}$ | $\boldsymbol{1}$ | $0$ | $1/3$ | $-$
$H$ | $\boldsymbol{1}$ | $\boldsymbol{2}$ | $-1/2$ | $0$ | $+$
$\Phi_{s}$ | $\boldsymbol{1}$ | $\boldsymbol{1}$ | $0$ | $2$ | $+$
Table 1:  The matter particle content of the minimal B-L model and their
respective charge assignments under each symmetry.

The Lagrangian that describes the DM phenomenology and encodes the process of
$B-L$ symmetry breaking can be written as

$\displaystyle\mathcal{L}$
$\displaystyle\supset\overline{\chi}\gamma^{\mu}\partial_{\mu}\chi-
m_{\chi}\overline{\chi}\chi$
$\displaystyle-\frac{1}{4}F^{\prime\,\mu\nu}F^{\prime}_{\mu\nu}+g_{BL}n_{\chi}\overline{\chi}\gamma^{\mu}\chi
Z^{\prime}_{\mu}$
$\displaystyle+g_{BL}n_{\ell}\sum_{\psi=\ell,\nu_{\ell},N_{iR}}\overline{\psi}\gamma^{\mu}\psi
Z^{\prime}_{\mu}+g_{BL}n_{q}\sum_{i=1}^{6}\overline{q}_{i}\gamma^{\mu}q_{i}Z^{\prime}_{\mu}$
$\displaystyle+y_{ij}^{D}\overline{L}_{i}\tilde{H}N_{jR}+y_{ij}^{M}\overline{(N^{C}_{iR})}\Phi_{s}N_{jR}$
$\displaystyle+(D_{\mu}\Phi_{s})^{\dagger}(D^{\mu}\Phi_{s})-V_{0}(\Phi_{s}),$
(1)

where $F^{\prime}_{\mu\nu}$ and $g_{BL}$ are the strength tensor and coupling
of the $B-L$ symmetry. The $n_{j}$ stands for the $B-L$ quantum number of the
particles $j=\chi,\ell,\nu_{\ell},N_{iR},$ and $q$, with $\ell=e,\mu,\tau$,
and $q=u,d,c,s,t,$ and $d$. We remark that the DM charge, $n_{\chi}$, must be
other than $\pm 1$ to avoid DM decay via an additional Yukawa term involving
$\chi_{R}$. Notice that $\tilde{H}=i\sigma_{2}H$ is the isospin transformation
of the SM Higgs doublet $H=\left(\phi^{+},\phi^{0}\right)^{T}$. Moreover, in
the scalar sector, we have the kinetic term of $\Phi_{s}$ which will give mass
to $Z^{\prime}$, and the scalar singlet potential at tree level
$V_{0}(\Phi)=\mu_{s}^{2}\Phi^{\dagger}_{s}\Phi_{s}+\lambda_{s}\left(\Phi^{\dagger}_{s}\Phi_{s}\right)/2$.
We leave further details about the effective potential
$V_{\textrm{eff}}(\Phi_{s})$ that leads to the first-order phase transition
for the next section.

The parametrization of the scalar field singlet is given by,

$\Phi_{s}=\frac{1}{\sqrt{2}}\left(v_{s}+\phi_{s}+i\rho_{s}\right),$ (2)

where $v_{s}=\sqrt{2}\langle\Phi_{s}\rangle$ is the $U(1)_{B-L}$ vacuum
expectation value (VEV), and $\rho_{s}$ is the Goldstone boson which will be
eaten by the $Z^{\prime}$ field after spontaneous symmetry breaking. Its mass
arises from the kinetic term of scalar singlet, and the right-handed neutrinos
$N_{iR}$ via the Majorana mass term in Eq. (1),

$\displaystyle m_{N_{iR}}$ $\displaystyle=\frac{y^{M}_{i}}{\sqrt{2}}v_{s},$
(3) $\displaystyle m_{Z^{\prime}}$ $\displaystyle=2g_{BL}v_{s}.$ (4)

We highlight that it happens at high-energy scales, say, $v_{s}\gg v$, where
$v$ is the VEV of the SM Higgs field, $H$. In the same way, the tree-level
mass of the scalar singlet is given by

$m_{\phi_{s}}=\sqrt{\lambda_{s}}v_{s}.$ (5)

The active neutrinos become massive via the popular type I seesaw mechanism,
which nicely reproduces the neutrino data Minkowski (1977); Mohapatra and
Senjanovic (1980); Brdar  _et al._ (2019).

Figure 1: The main channels for DM-SM interactions in this $B-L$ DM model: (a)
DM annihilation into SM fermions $f$; (b) DM annihilation into on-shell pair
of $Z^{\prime}$, which is only relevant in the limit
$m_{\chi}>m_{Z^{\prime}}$; and (c) DM-nucleon scattering process, where $q$
stands for quarks.

In summary, the free parameters that govern the DM phenomenology are the DM
mass $m_{\chi}$, the $Z^{\prime}$ mass $m_{Z^{\prime}}$, and the $B-L$
coupling $g_{BL}$. In Fig. 1, we show the relevant Feynman diagrams. In the
next section, we assess the GWs production in our model due to a first-order
phase transition, where $g_{BL}$ and $v_{s}$ play a crucial role in the
spectra.

##  III  Stochastic GW spectrum from first-order phase transition

If the $U(1)_{B-L}$ breaking is a first-order phase transition, there are
regions of broken phase in the plasma where the symmetry remains unbroken.
These so-called bubbles expand and induce plasmatic motion in the form of
sound waves and turbulence. At the end of the transition, when the bubbles
collide and fill the entire space, GWs are produced due to a time-varying
quadrupole moment in the kinetic energy-momentum of the plasma Caprini and
Figueroa (2018); Caprini  _et al._ (2020).

To estimate the shape of this spectrum we need to study the dynamics of the
phase transition. In the presence of a thermal plasma, the dynamics of the
scalar field $\phi_{s}$ is described by an effective potential which, at
1-loop order, takes the form

$\begin{split}V_{\text{eff}}&(T,\phi_{s})=\dfrac{\lambda_{s}}{8}(\phi_{s}^{2}-v_{s}^{2})^{2}+V_{\text{th}}(T,\phi_{s})\\\
&+\dfrac{3\times(2g_{BL})^{4}}{64\pi^{2}}\left[\phi_{s}^{4}\left(\log\dfrac{\phi_{s}^{2}}{v_{s}^{2}}-\dfrac{3}{2}\right)+2\phi_{s}^{2}v_{s}^{2}\right].\end{split}$
(6)

The second line corresponds to the Coleman-Weinberg potential plus counter-
terms added to ensure that the minimum of the 1-loop zero-temperature
potential remains at $v_{s}$, and that the mass of the scalar is still given
by Eq. (5). We take into account the $Z^{\prime}$ running in the loop, and
neglect the right-handed neutrinos (assuming small Yukawas111Including the
right-handed neutrinos would amount to shifting the prefactor in the second
line of Eq. (6) to
$3\\!\times\\!(2g_{BL})^{4}-2\\!\times\\!\sum_{i}(y_{i}^{M})^{4}/4$. Here we
assume $\sum_{i}(y_{i}^{M})^{4}\ll 96g_{BL}^{4}$.) and the scalar (since
typically $\lambda_{s}\ll g_{BL}$ and it has only one degree of freedom, so
its effect is subdominant against the $Z^{\prime}$). The thermal part
$V_{\text{th}}$ is computed via standard methods of thermal field theory Dolan
and Jackiw (1974); Hindmarsh  _et al._ (2021).

Notice that the value of the zero-temperature potential at the unbroken state
$\phi_{s}=0$ is parametrized by $\lambda_{s}$, whereas at the broken vacuum
$\phi_{s}=v_{s}$ it goes with $g_{BL}$. Hence, for too small values of
$\lambda_{s}$ the symmetric minimum may become the lowest energy state and
$U(1)_{B-L}$ would remain unbroken, which is unphysical and should be avoided.
Explicitly, defining the energy difference between the two vacua $\Delta
V(T,\phi)\equiv V_{\text{eff}}(T,\phi)-V_{\text{eff}}(T,0)$, one finds that at
zero-temperature

$\Delta
V(0,v_{s})=\dfrac{3}{8\pi^{2}}g_{BL}^{4}v_{s}^{4}-\dfrac{\lambda_{s}}{8}v_{s}^{4}<0,$
(7)

implying $\lambda_{s}>3g_{BL}^{4}/\pi^{2}$.

The bubble nucleation rate per unit volume is
$\Gamma_{\text{nuc}}/\mathcal{V}\sim T^{4}e^{-S_{3}/T}$, where

$S_{3}=4\pi\int
dr\,r^{2}\left[\frac{1}{2}\left(\dfrac{d\phi_{c}}{dr}\right)^{2}+\Delta
V(T,\phi_{c})\right]$ (8)

is the Euclidean action of the critical bubble, corresponding to the saddle
point configuration that connects the two vacua in field space Hindmarsh  _et
al._ (2021); Mukhanov (2005). At some temperature $T$, this configuration can
be found by solving the bounce equation
$-\nabla^{2}\phi+dV_{\text{eff}}/d\phi=0$ using a shooting method. The
nucleation temperature is found by imposing that
$\Gamma_{\text{nuc}}/\mathcal{V}$ equals the Hubble expansion rate. This
corresponds roughly to the temperature at which $S_{3}/T\approx 140$ Caprini
_et al._ (2020). We can then define the fractional amount of energy released
by the transition Espinosa  _et al._ (2010) to be,

$\alpha\equiv\dfrac{Q_{\text{lat}}-3\Delta V}{4\rho_{\text{rad}}},$ (9)

where $Q_{\text{lat}}=Td\Delta V/dT-\Delta T$ is the latent heat of the phase
transition, the numerator $Q_{\text{lat}}-3\Delta V$ is the difference of the
trace of the energy-momentum tensor in both phases Espinosa  _et al._ (2010);
Giese  _et al._ (2020), $\rho_{\text{rad}}=\pi^{2}g_{\text{eff}}T^{4}/90$ is
the radiation energy density at the time of the transition and
$g_{\text{eff}}=116.5$ is the effective number of relativistic degrees of
freedom in the plasma. The amplitude of the GW spectrum will depend on how
efficiently this energy is converted into fluid motion of the plasma. For that
we construct efficiency factors $\kappa_{\text{sw}}$ and
$\kappa_{\text{turb}}$ such that the fractional energy in sound waves
(respectively in plasmatic turbulence) is proportional to
$\kappa_{\text{sw}}\alpha$ (resp. $\kappa_{\text{turb}}\alpha$). An
approximate formula for $\kappa_{\text{sw}}$ can be found in ref. Espinosa
_et al._ (2010), but for turbulence this conversion factor is unknown.
Sometimes it is estimated to be $1-10\%$ of $\kappa_{\text{sw}}$ Caprini  _et
al._ (2020), and its contribution is subdominant Caprini  _et al._ (2016). Due
to this indeterminacy and subdominance, we neglect here the contribution from
turbulence and consider only sound waves. Note that, by doing this, we
underestimate the spectrum, so our results are conservative.

Another important parameter for estimating the GW spectrum is the (inverse)
duration of the transition, estimated as222For the phase transitions
considered here, $\beta/H\sim 100-1000\gg 1$. The spectrum of GWs from sound
waves scales as $(\beta/H)^{-1}$ whereas bubble collisions is damped by
$(\beta/H)^{-2}$, which is comparatively negligible Caprini  _et al._ (2016).

$\dfrac{\beta}{H}\equiv T\dfrac{d(S_{3}/T)}{dT}.$ (10)

Given that the bubble expands at a velocity $v_{w}$, its radius at collision
will be proportional to $v_{w}(H/\beta)$, and the larger the bubble, the
larger the GW amplitude. Calculating the wall velocity $v_{w}$ is a daunting
task, since it involves non-equilibrium phenomena and depends on how we model
the plasma away from equilibrium. There have been recent discussions in the
literature on the appropriate way to achieve this description, but the debate
is still unsettled Cline and Kainulainen (2020); Dorsch  _et al._ (2022);
Laurent and Cline (2022). Here we approximate $v_{w}=1$ for simplicity, which
should suffice for an adequate estimate of the spectra at the correct order of
magnitude.

More specifically, for the amplitude of the GW spectrum from sound waves we
find Caprini  _et al._ (2016),

$\begin{split}\Omega_{\text{sw}}h^{2}=&\leavevmode\nobreak\ 2.65\times
10^{-6}\,v_{w}\left(\dfrac{H}{\beta}\right)\left(\dfrac{\kappa_{\text{sw}}\alpha}{1+\alpha}\right)^{2}\left(\dfrac{100}{g_{\text{eff}}}\right)^{1/3}\times\\\
&\quad\times\left(\dfrac{f}{f_{\text{peak}}}\right)\left(\dfrac{7}{4+3(f/f_{\text{peak}})}\right)^{7/2}.\end{split}$
(11)

This spectrum has a peak at

$f_{\text{peak}}=1.9\times
10^{-2}\text{mHz}\,\left(\dfrac{1}{v_{w}}\dfrac{\beta}{H}\right)\left(\dfrac{T}{100\leavevmode\nobreak\
\text{GeV}}\right)\left(\dfrac{g_{\text{eff}}}{100}\right)^{1/6}.$ (12)

We will see that, for typical benchmark values of the parameters in our model,
the peak frequency lies in the mHz band, hence possibly within reach of future
interferometers such as LISA, DECIGO and BBO. We will now discuss the dark
matter relic density and scattering rate to later put our findings into
perspective.

##  IV Dark Matter Relic Abundance

In this section, we assess the production of DM particles in the standard
thermal freeze-out paradigm. In such a narrative, after reheating, the DM
particles were in thermal equilibrium with the SM particles, which means that
they were pair-annihilated and pair-produced in equal rate in the early
universe. However, as the universe expands and cools down, the expansion rate
increases, whereas the interaction rate decreases. When the expansion rate
overcomes the interaction rates, freeze-out takes place and the DM abundance
becomes constant.

In this paradigm, the evolution of the DM number density $n_{DM}$ is described
by the Boltzmann equation,

$\frac{dY_{DM}(x)}{dx}=-\frac{s(x)\langle\sigma
v\rangle}{x\,H(x)}\left[Y^{2}_{DM}(x)-Y^{\textrm{eq}\,2}_{DM}(x)\right]$ (13)

where $Y_{DM}=n_{DM}/s$ is the comoving number density, with

$s(x)=\frac{2\pi^{2}}{45}g_{\star s}(x)m^{3}_{DM}x^{-3}$ (14)

representing the entropy of the primordial plasma, with $g_{\star s}$ being
the relativistic degrees of freedom that contribute to the entropy, and
$x=m_{DM}/T$ is a “time” variable that helps us to simplify the integration
and physical interpretations. In this parametrization, the Hubble expansion
rate and comoving abundance are written as,

$H(x)=\frac{\pi}{M_{pl}}\sqrt{\frac{g_{\star}(x)}{90}}m_{DM}^{2}x^{-2},$ (15)

$Y^{\textrm{eq}}_{DM}(x)=\frac{45}{4\pi^{4}}\frac{g_{DM}}{g_{\star
s}}x^{2}K_{2}(x),$ (16)

where $g_{\star}$ accounts for the relativistic degrees of freedom, with
$g_{\star s}\approx g_{\star}=106.75$ at the time of freeze-out, $g_{DM}$
corresponds to the DM degrees of freedom, and $K_{2}(x)$ is the modified
Bessel function. After solving Eq. (13), we obtain the DM abundance,

$\Omega_{DM}h^{2}\simeq 2.82\times 10^{8}m_{DM}Y_{DM}(x\rightarrow\infty),$
(17)

where we use the limit of $x\rightarrow\infty$ in the solution of the
Boltzmann equation. The most current value is given by Planck collaboration,
$\Omega_{DM}h^{2}=0.1200\pm 0.0012$ within 68% C.L. Aghanim  _et al._ (2020).

In Section II, we described how a vector-like Dirac fermion $\chi$ can be the
DM candidate in a minimal $B-L$ model. We computed the DM relic density using
micrOMEGAs Belanger  _et al._ (2007, 2021). In Fig. 2, we present the behavior
of the relic abundance as a function of the DM mass. We have set $v_{s}=7$
TeV. The red and green solid curves are the relic abundance for $g_{BL}=0.45$
and $g_{BL}=0.80$ (corresponding respectively to $m_{Z^{\prime}}=6.3$ TeV and
$m_{Z^{\prime}}=11.2$ TeV, according to Eq. (4)). The gray horizontal line
delimits the region that reproduces the Planck data Aghanim  _et al._ (2020).

It is remarkable that the largest part of the parameter space of the thermal
relic abundance is overabundant. It occurs because, in general, the model
provides small annihilation cross-sections. However, at the resonance regime,
when $m\chi\approx m_{Z^{\prime}}/2$, it reaches sizeable values via the
s-channel in Fig. 1 (a). Such an enhancement brings the relic density down to
the observed value. Because of the (inverted) resonance peak, we see that
there are typically two viable DM masses: one slightly smaller than
$m_{Z^{\prime}}/2$, the other slightly larger.

Thus far, we have addressed how to produce dark matter and GWs. In what
follows, we will show that our model is amenable to collider and direct
detection constraints.

Figure 2: DM relic abundance as a function of the DM mass. We set $v_{s}=7$
TeV. The red solid curve is the relic abundance for $g_{BL}=0.45$ and
$m_{Z^{\prime}}=6.3$ TeV, while the green one represents $g_{BL}=0.80$ and
$m_{Z^{\prime}}=11.2$ TeV. The gray line is the abundance in agreement with
Planck observations Aghanim  _et al._ (2020). Notice that there are two viable
DM masses.

##  V  The Constraints

Let us now turn to a discussion of the limits on the $Z^{\prime}$ mass from
the ATLAS and LEP-II results, and direct detection bounds on the spin
independent (SI) cross-section as reported by XENONnT and LUX-ZEPLIN (LZ).

###  V.1  Collider limits

Extra gauge bosons with unsuppressed coupling to the SM fermions, as the one
considered in this work, produce strong collider signals in the form of
resonances decaying into SM pairs as for example dileptons Aad  _et al._
(2019), dijets Aad  _et al._ (2020) and di-top Aaboud  _et al._ (2019). We
will adopt the most restrictive limits to our model, which stems from searches
for heavy dilepton events conducted by ATLAS collaboration during Run 2 of the
Large Hadron Collider, and corresponds to an integrated luminosity of
$139fb^{-1}$ Aad  _et al._ (2019). It is well-known that the presence of
$Z^{\prime}$ couplings to DM pairs can weaken the latter constraints Arcadi
_et al._ (2014), however, such invisible decay yields no effect on the
dilepton bound, because of its small contribution to the total $Z^{\prime}$
width Klasen  _et al._ (2017). ATLAS collaboration has not searched for a B-L
$Z^{\prime}$ boson. Therefore, we had to derive our own limit by comparing the
theoretical prediction with the upper limit derived at 95% CL on the fiducial
cross-section times branching ratio quoted in Aad  _et al._ (2019).

For each combination of gauge coupling and $Z^{\prime}$ mass, we computed the
$Z^{\prime}$ branching ratio using CalcHEP Belyaev  _et al._ (2013) and fed
this information into Madgraph5 Alwall  _et al._ (2014); Frederix  _et al._
(2018) where we performed the Monte Carlo simulation adopting a parton
distribution function (PDF) NNPDF23LO Carrazza  _et al._ (2013). We computed
the signal event $pp\to Z^{\prime}\to\ell\bar{\ell}$ at $\sqrt{s}=13$ TeV,
with $\ell=e,\mu$ and compared our result with the public result from ATLAS
Collaboration Aad  _et al._ (2019). Following the collaboration, we required
the signal events to feature two opposite charge leptons, with transverse
momentum $p_{\mathrm{T}}>30$ GeV, and pseudorapidity $|\eta|<2.5$. In doing
so, we obtained the collider limits shown in Table .2. We emphasize that these
findings represent a new result in the literature. We highlight that the limit
for $g_{BL}=0.8$, the production cross-section times branching ratio falls
outside the sensitivity of ATLAS data. This happens because when we increase
the gauge coupling, we also increase the production cross section of
$Z^{\prime}$ bosons. This shift upward in the cross-section makes the
theoretical prediction cross the exclusion limit from ATLAS at larger
$Z^{\prime}$ masses. As the highest pole mass probed by ATLAS was $6$ TeV,
when we significantly increase the gauge coupling, our model can no longer be
excluded by it. For this particular case, we applied the Collider Reach
$\beta$ tool Thamm  _et al._ (2015) which allows us to forecast the new bounds
on the $Z^{\prime}$ mass for a different collider configuration. Here we
maintained $\sqrt{s}=13$ TeV, and simply ramped up the luminosity and selected
to be $\mathcal{L}=300fb^{-1}$. Doing so, we project the ATLAS bound of
$m_{Z^{\prime}}>7.7$ TeV. This is a simple estimation of the ATLAS reach. For
$g_{BL}=0.8$, we could have adopted a more solid and weaker bound,
$m_{Z^{\prime}}>5.6$ TeV, from the Large Electron Positron (LEP) collider that
will address below.

Gauge coupling | Lower bound - ATLAS 13TeV
---|---
$g_{BL}=0.2$ |  $m_{Z^{\prime}}>4.94$ TeV
$g_{BL}=0.3$ |  $m_{Z^{\prime}}>5.35$ TeV
$g_{BL}=0.4$ |  $m_{Z^{\prime}}>5.62$ TeV
$g_{BL}=0.45$ |  $m_{Z^{\prime}}>5.75$ TeV
$g_{BL}=0.5$ |  $m_{Z^{\prime}}>5.8$ TeV
$g_{BL}=0.6$ |  $m_{Z^{\prime}}>5.97$ TeV
$g_{BL}=0.7$ |  $m_{Z^{\prime}}>6$ TeV
$g_{BL}=0.8$ |  $m_{Z^{\prime}}>7.7$ TeV
Table 2: Lower mass bounds derived on the $Z^{\prime}$ mass for different
gauge couplings using ATLAS results reported in Aad  _et al._ (2019).

An old and relevant collider bound comes from LEP-II data that reads
$m_{Z^{\prime}}/g_{BL}>7\,\mbox{TeV}$ Carena  _et al._ (2004); Cacciapaglia
_et al._ (2006). This limit is not affected by an eventual presence of an
invisible branching fraction. It stems from the comparison between the SM
prediction and new physics prediction for dilepton events, rather than from
resonance searches. As LEP featured fantastic precision due to the leptonic
nature of the process, they were able to obtain a stringent limit on the
$Z^{\prime}$ mass. LEP bound is more relevant for $g_{BL}\sim 1$. We explore a
setup with $g_{BL}=0.8$, and as we discussed previously, current data from
ATLAS cannot probe this case. For this reason, for $g_{BL}=0.8$ we also quote
the limit from LEP that reads $m_{Z^{\prime}}/g_{BL}>5.6\,\mbox{TeV}$.

###  V.2  Direct detection bounds

Since the $Z^{\prime}$ boson features vector coupling with both quark and DM
pairs, spin-independent interactions arise between the latter and the
nucleons. The corresponding cross-section can be written as

$\sigma_{\chi N=p,n}^{\rm SI}=\frac{\mu_{\chi
N}^{2}}{\pi}\frac{9\,n_{q}^{2}n_{\chi}^{2}g_{BL}^{4}}{m_{Z^{\prime}}^{4}}$
(18)

with $\mu_{\chi N}=\frac{m_{\chi}m_{N}}{m_{\chi}+m_{N}}$ being the DM-nucleon
reduced mass. Such interactions are strongly constrained by xenon-based direct
detection experiments. For our study, we will consider the most recent bounds
as given by the LZ Aalbers  _et al._ (2022) and XENONnT Aprile  _et al._
(2023). Notice also that the vector coupling with SM fermions implies an
s-wave dominated annihilation cross-section into SM fermions, hence
potentially testable via indirect detection. One could then use searches of
$\gamma$-ray signals, see e.g. Ackermann  _et al._ (2015); Acharyya  _et al._
(2021); Hooper  _et al._ (2013); Acharya  _et al._ (2018), to further
constrain the parameter of the model. Indirect detection constraints are,
however, not competitive with direct detection and collider for the model
under scrutiny. Consequently, we will not show them explicitly. That said, we
will now put our findings into perspective.

##  VI  Results

Our main results concerning the DM phenomenology are displayed in Fig. 3 for
$g_{BL}=0.45$ (top panel) and $g_{BL}=0.80$ (bottom panel). The red (top) and
green (bottom) curves yield the correct relic abundance, in agreement with
Planck’s observations. The region in between the curves gives an underabundant
relic, whereas the outside region is an overabundant one. Furthermore, notice
that one value of $Z^{\prime}$ mass is associated with two values of DM mass
that yield the correct relic density. This happens because the model reaches
$\Omega_{\chi}h^{2}\approx 0.11$ in the $Z^{\prime}$ resonance regime.

We also exhibit the direct detection constraint from the LZ collaboration that
is slightly stronger than the XENONnT one. The shape of the experiment curve
can be understood looking at Eq. (18). When the dark matter mass is much
larger than the nucleus mass, the scattering cross section is independent of
the dark matter mass, but decreases with $m_{Z^{\prime 4}}$. As the
experimental limit linearly weakens with the dark matter, the viable parameter
space in the $m_{\chi}$ vs $m_{Z^{\prime}}$ plane is mostly sensitive to the
$Z^{\prime}$ mass. We are plotting our findings in a Log-Log scale, thus what
we see is a line representing the direct detection bound weakening as the
$Z^{\prime}$ mass increases.

As for the collider constraint from ATLAS, $m_{Z^{\prime}}>6$ TeV. It is
simply a vertical line on the $Z^{\prime}$ mass once we fix the gauge
coupling. This has to do with the fact that the signal scales with
$g_{BL}^{2}\times BR(Z^{\prime}\rightarrow l\bar{l})$. Hence, once we fix
$g_{BL}$ the bound on the $Z^{\prime}$ can be directly extracted from the
collaboration report. The branching of the $Z^{\prime}$ into dileptons will
not change when $m_{Z^{\prime}}>2m_{\chi}$, because the partial width of
$Z^{\prime}\rightarrow\chi\bar{\chi}$ is also controlled by $g_{BL}$.
Therefore, when we open the invisible decay into dark matter, no meaningful
change in the branching into lepton is expected. In models, where the
$Z^{\prime}$ coupling to dark matter is different from the $Z^{\prime}$
coupling to leptons, a large decay into invisible can be more pronounced,
however, Alves  _et al._ (2014, 2015a, 2015b).

Figure 3: Combined results for DM phenomenology. The solid curves correspond
to the correct DM relic abundance $\Omega_{\chi}h^{2}\approx 0.12$, for
$g_{BL}=0.45$ (top) and $g_{BL}=0.80$ (bottom). The region between the two
curves corresponds to underabundance and is in principle allowed, while
outside the curves the Universe would overclose. The shaded blue region
represents the parameter space excluded by direct detection. In orange, we
have the corresponding collider limits, which substantially constrain the
model.

Notice that in Fig. 3, $v_{s}$ varies because the $Z^{\prime}$ mass changes
linearly with it. We are bringing this to the reader’s attention because the
value of $v_{s}$ will be rather relevant to the GW spectra. For this reason,
we drew a dashed vertical line representing $v_{s}=7$ TeV, which realizes the
GW spectra to be explored later in Fig. 5. In the same way, we benchmark with
blue and violet stars DM masses associated to those GW spectra. In the top
panel, $m_{\chi}=2.9$ or $3.4$ TeV, while in the bottom one $m_{\chi}$ has to
be $5.0$ or $6.15$ TeV. The shift upward in the relic density curve has to do
with the $Z^{\prime}$ width, which is larger as it grows with $g_{BL}$.

In Fig. 4 we fix $v_{s}=7$ TeV and obtain the relic density curve in the
$m_{\chi}$ vs $g_{BL}$ plane. The shaded region is ruled out by the LZ
collaboration. From Eq. (18) we notice that the DM-nucleon SI scattering
cross-section will depend only on $m_{\chi}$, since the dependence on $g_{BL}$
enters only in the ratio $\left(g_{BL}/m_{Z^{\prime}}\right)^{4}$, which is
constant for fixed $v_{s}$. As we are considering $m_{\chi}\gg m_{N}$, the
scattering cross-section no longer depends on the dark matter mass, and is
consequently constant in the plane of Fig. 4. Although, the experimental limit
from direct detection experiments linearly weakens with the dark matter mass.
In particular, we found a dark matter-nucleon scattering cross section of
$3.1\times 10^{-10}$ pb, which is consistent with the experiment limit that
reads $3.2\times 10^{-10}$ pb for $m_{\chi}\simeq 1.35$ TeV. Therefore, for
larger dark matter masses, the direct detection bound becomes too weak to
constrain the parameter space in Fig. 4. This explains why the exclusion
region from direct detection represents a horizontal blue curve.

As for collider bound, we show the collider bound following Table 2.

Figure 4: The observed DM relic for $v_{S}=7$ TeV is represented by the black
curve. The shaded orange and blue regions are the ATLAS and LZ-ZEPLIN limits,
respectively. The ATLAS bound reproduces Table 2. Direct detection imposes
$m_{\chi}>1.35$ TeV.

The above discussion highlights how DM searches can shed light on the
properties of the DM particle $\chi$, the $Z^{\prime}$ mediator, and the gauge
coupling $g_{BL}$. But these searches are basically insensitive to the
properties of the scalar $\phi_{s}$, such as its mass or self-coupling
$\lambda_{s}$, which do not affect the DM density. Moreover, if
$m_{\phi_{s}}\lesssim 2m_{Z^{\prime}}$ (i.e. $\lambda_{s}\lesssim
16g_{BL}^{2}$), the scalar will only decay to the invisible right-handed
neutrinos and will be hardly detectable at colliders. Fortunately, in this
case, one might constrain the scalar sector using gravitational waves, and an
eventual detection of GWs would allow us to determine the scalar self-
coupling.

This is shown in Fig. 5, where we display the GW spectra for the benchmark
value of $v_{s}=7$ TeV with $g_{BL}=0.45$ (top) and $g_{BL}=0.80$ (bottom).
The different spectra in each figure correspond to varying $\lambda_{s}$,
whose corresponding values are shown along each curve. Notice that decreasing
$\lambda_{s}$ generates larger spectra. This is because smaller $\lambda_{s}$
will lead to a smaller energy difference in the zero-temperature potential
between the broken and unbroken minima (cf. Eq. (7)), which leads to stronger
phase transitions Dorsch  _et al._ (2017). Eventually, the energy gap is so
small that the transition to the true minimum never occurs: the field remains
trapped in the metastable false vacuum and the symmetry is not broken, which
is obviously non-physical. Hence, there is an upper bound on the spectrum that
could be achieved for each pair $(v_{s},g_{BL})$.

Varying $v_{s}$ will shift the peak frequency, but since the plot is
logarithmic in $f$, one would need large differences in $v_{s}$ for this shift
to be noticeable. The differences in the spectra will not be significant for
the range of $v_{s}$ values considered here (cf. the allowed region of DM
production in Fig. 3).

Figure 5: Spectra of GWs for $v_{s}=7$ TeV. The figure at the top
(respectively bottom) refers to $g_{BL}=0.45$ (resp. $g_{BL}=0.80$). The
various spectra displayed are obtained by varying $\lambda_{s}$ at values
shown in each curve. Sensitivity curves for LISA, DECIGO and BBO have been
obtained from Ref. Schmitz (2021).

Fig. 6 shows the maximum values of $\lambda_{s}$, for given $g_{BL}$ and fixed
$v_{s}=7$ TeV, that would yield a detectable spectrum at BBO, DECIGO and LISA,
respectively. In other words, the parameter space below the curves could be
probed using gravitational waves. We also show the regions where the unbroken
state is metastable (red region) and stable (gray region). For concreteness,
considering BBO sensitivity, the region that lies between the dotted curve and
the metastability yields a detectable gravitation wave signal. Thus, the
detectability at LISA is achieved very close to the metastable situation. This
means that, for the perturbative range of $g_{BL}$ shown in the plot, one can
find spectra within the LISA sensitivity, but not by a huge margin. More
sensitive detectors, such as DECIGO and BBO, would be able to probe a larger
range of $\lambda_{s}$ values.

##  VII Discussions

We have shown that the combination $\lambda_{s}-g_{BL}$ governs the
gravitation wave signal, while $g_{BL}$ and $m_{Z^{\prime}}$, $m_{\chi}$ the
dark matter phenomenology. Therefore, if colliders happen to observe a
dilepton resonance consistent with a B-L model, both the $Z^{\prime}$ mass and
gauge couplings will be extracted. In the minimal B-L model, it means that
$v_{s}$ is also known. The benchmark points highlighted in Fig. 3 demonstrate
the power of exploring the correlation between colliders, direct detection,
and cosmology because we know the possible values for the dark matter masses
that satisfy those conditions. Furthermore, using this information we can
predict the spectrum of gravitational waves to be observed in future probes
depending on the value of quartic coupling $\lambda_{s}$. Ideally, if a
gravitational signal is observed, we can determine which combination of the
gauge coupling and $\lambda_{s}$ that could reproduce the data. After that,
one can easily check whether this value obtained for the gauge couplings would
produce a thermal dark matter particle with the correct relic density in
agreement with the data. In other words, gravitational waves have the power to
constrain the scalar self-coupling as well as the gauge coupling of the
theory, and through that indirectly probe dark sectors. Excitingly,
gravitational waves offer a rather orthogonal way to probe the dark sectors.

Figure 6: Maximum values of $\lambda_{s}$, for given $g_{BL}$ and $v_{s}=7$
TeV, that would yield a GW spectrum detectable at BBO (dotted curve), DECIGO
(dashed) and LISA (solid). The red region shows the values of $\lambda_{s}$
that lead to a metastable false vacuum, whereas the gray region corresponds to
stability of the unbroken state, i.e. the broken phase not being the global
minimum of the potential, as per Eq. (7).

##  VIII  Conclusions

In this work, we studied how a minimal $B-L$ model can provide both a viable
DM candidate and a stochastic spectrum of GWs that might be detected in future
experiments. The DM candidate is a vector-like fermion and the $U(1)_{B-L}$
symmetry is broken by a scalar singlet and induces a first-order cosmological
phase transition.

The dynamics of the cosmological phase transition and its gravitational wave
imprints are highly sensitive to the scalar effective potential, so
gravitational wave interferometers could be used to extract information on the
self-coupling $\lambda_{s}$ and $g_{BL}$. Interestingly, the expected LISA
sensitivity is barely enough to probe the most extreme cases, when the phase
transition is extremely strong and the unbroken vacuum is on the verge of
becoming metastable. On the other hand, more sensitive interferometers such as
DECIGO and BBO might be able to probe a larger range of self-couplings. There
is arguably some degree of tuning for gravitational wave detection, in the
sense that an exiguous change in the value of $\lambda_{s}$ might move the
spectrum from within the detectability range of LISA down to non-detectability
even at BBO. But this also means that a gravitational wave detection, allied
with measurements of $v_{s}$ and $g_{BL}$ from other experiments, would lead
to a measurement of the scalar self-coupling with good precision.

Colliders can probe the $Z^{\prime}$ mass and the gauge coupling $g_{BL}$, but
are insensitive to details of the scalar singlet. Direct detection experiments
rule out a larger fraction of the parameter space, which help us constrain the
dark sector. Knowing that, if we impose a gravitational wave detection in
future probes and a thermal production of dark matter, we can predict which
dark matter masses reproduce the correct relic density in agreement with the
data. In summary, our main result is that gravitational wave detectors offer a
complementary and orthogonal probe to dark sectors, allowing us to further
narrow down the parameter space of dark matter models.

###### Acknowledgements.

The authors thank Diego Cogollo, Carlos Pires, Carlos Yaguna for discussions.
JPN acknowledges support from Coordenação de Aperfeiçoamento de Pessoal de
Nível Superior (CAPES) under Grant No. 88887.712383/2022-00. Y.M.O.T.
acknowledges financial support from CAPES under grants 88887.485509 / 2020-00.
This work was financially supported by Simons Foundation (Award
Number:1023171-RC), FAPESP Grant 2021/01089-1, ICTP-SAIFR FAPESP Grants
2021/14335-0, CNPq Grant 307130/2021-5, FONDECYT Grant 1191103 (Chile) and
ANID-Programa Milenio-code ICN2019_044.

## References

  * Arcadi  _et al._ (2018) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, Eur. Phys. J. C 78, 203 (2018), arXiv:1703.07364 [hep-ph] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, VIRGO), (2021a), arXiv:2108.01045 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03606 [gr-qc] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Crowder and Cornish (2005) J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005), arXiv:gr-qc/0506015 .
  * Kawamura  _et al._ (2006) S. Kawamura _et al._ , Class. Quant. Grav. 23, S125 (2006).
  * Aprile  _et al._ (2023) E. Aprile _et al._ (XENON), (2023), arXiv:2303.14729 [hep-ex] .
  * Aalbers  _et al._ (2022) J. Aalbers _et al._ (LZ), (2022), arXiv:2207.03764 [hep-ex] .
  * Aad  _et al._ (2019) G. Aad _et al._ (ATLAS), Phys. Lett. B 796, 68 (2019), arXiv:1903.06248 [hep-ex] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), JHEP 03, 145 (2020), arXiv:1910.08447 [hep-ex] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), Phys. Rev. D 99, 092004 (2019), arXiv:1902.10077 [hep-ex] .
  * Hasegawa  _et al._ (2019) T. Hasegawa, N. Okada, and O. Seto, Phys. Rev. D 99, 095039 (2019), arXiv:1904.03020 [hep-ph] .
  * Bosch  _et al._ (2023) J. Bosch, Z. Delgado, B. Fornal, and A. Leon, (2023), arXiv:2306.00332 [hep-ph] .
  * Klasen  _et al._ (2017) M. Klasen, F. Lyonnet, and F. S. Queiroz, Eur. Phys. J. C 77, 348 (2017), arXiv:1607.06468 [hep-ph] .
  * Minkowski (1977) P. Minkowski, Phys. Lett. B 67, 421 (1977).
  * Mohapatra and Senjanovic (1980) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * Brdar  _et al._ (2019) V. Brdar, A. J. Helmboldt, S. Iwamoto, and K. Schmitz, Phys. Rev. D 100, 075029 (2019), arXiv:1905.12634 [hep-ph] .
  * Caprini and Figueroa (2018) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Caprini  _et al._ (2020) C. Caprini _et al._ , JCAP 03, 024 (2020), arXiv:1910.13125 [astro-ph.CO] .
  * Dolan and Jackiw (1974) L. Dolan and R. Jackiw, Phys. Rev. D 9, 3320 (1974).
  * Hindmarsh  _et al._ (2021) M. B. Hindmarsh, M. Lüben, J. Lumma, and M. Pauly, SciPost Phys. Lect. Notes 24, 1 (2021), arXiv:2008.09136 [astro-ph.CO] .
  * Mukhanov (2005) V. Mukhanov,  _Physical Foundations of Cosmology_ (Cambridge University Press, Oxford, 2005).
  * Espinosa  _et al._ (2010) J. R. Espinosa, T. Konstandin, J. M. No, and G. Servant, JCAP 06, 028 (2010), arXiv:1004.4187 [hep-ph] .
  * Giese  _et al._ (2020) F. Giese, T. Konstandin, and J. van de Vis, JCAP 07, 057 (2020), arXiv:2004.06995 [astro-ph.CO] .
  * Caprini  _et al._ (2016) C. Caprini _et al._ , JCAP 04, 001 (2016), arXiv:1512.06239 [astro-ph.CO] .
  * Cline and Kainulainen (2020) J. M. Cline and K. Kainulainen, Phys. Rev. D 101, 063525 (2020), arXiv:2001.00568 [hep-ph] .
  * Dorsch  _et al._ (2022) G. C. Dorsch, S. J. Huber, and T. Konstandin, JCAP 04, 010 (2022), arXiv:2112.12548 [hep-ph] .
  * Laurent and Cline (2022) B. Laurent and J. M. Cline, Phys. Rev. D 106, 023501 (2022), arXiv:2204.13120 [hep-ph] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Belanger  _et al._ (2007) G. Belanger, F. Boudjema, A. Pukhov, and A. Semenov, Comput. Phys. Commun. 176, 367 (2007), arXiv:hep-ph/0607059 .
  * Belanger  _et al._ (2021) G. Belanger, A. Mjallal, and A. Pukhov, Eur. Phys. J. C 81, 239 (2021), arXiv:2003.08621 [hep-ph] .
  * Arcadi  _et al._ (2014) G. Arcadi, Y. Mambrini, M. H. G. Tytgat, and B. Zaldivar, JHEP 03, 134 (2014), arXiv:1401.0221 [hep-ph] .
  * Belyaev  _et al._ (2013) A. Belyaev, N. D. Christensen, and A. Pukhov, Comput. Phys. Commun. 184, 1729 (2013), arXiv:1207.6082 [hep-ph] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Frederix  _et al._ (2018) R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H. S. Shao, and M. Zaro, JHEP 07, 185 (2018), [Erratum: JHEP 11, 085 (2021)], arXiv:1804.10017 [hep-ph] .
  * Carrazza  _et al._ (2013) S. Carrazza, S. Forte, and J. Rojo, in  _43rd International Symposium on Multiparticle Dynamics_ (2013) pp. 89–96, arXiv:1311.5887 [hep-ph] .
  * Thamm  _et al._ (2015) A. Thamm, R. Torre, and A. Wulzer, JHEP 07, 100 (2015), arXiv:1502.01701 [hep-ph] .
  * Carena  _et al._ (2004) M. Carena, A. Daleo, B. A. Dobrescu, and T. M. P. Tait, Phys. Rev. D 70, 093009 (2004), arXiv:hep-ph/0408098 .
  * Cacciapaglia  _et al._ (2006) G. Cacciapaglia, C. Csaki, G. Marandella, and A. Strumia, Phys. Rev. D 74, 033011 (2006), arXiv:hep-ph/0604111 .
  * Ackermann  _et al._ (2015) M. Ackermann _et al._ (Fermi-LAT), Phys. Rev. Lett. 115, 231301 (2015), arXiv:1503.02641 [astro-ph.HE] .
  * Acharyya  _et al._ (2021) A. Acharyya _et al._ (CTA), JCAP 01, 057 (2021), arXiv:2007.16129 [astro-ph.HE] .
  * Hooper  _et al._ (2013) D. Hooper, C. Kelso, and F. S. Queiroz, Astropart. Phys. 46, 55 (2013), arXiv:1209.3015 [astro-ph.HE] .
  * Acharya  _et al._ (2018) B. S. Acharya _et al._ (CTA Consortium),  _Science with the Cherenkov Telescope Array_ (WSP, 2018) arXiv:1709.07997 [astro-ph.IM] .
  * Alves  _et al._ (2014) A. Alves, S. Profumo, and F. S. Queiroz, JHEP 04, 063 (2014), arXiv:1312.5281 [hep-ph] .
  * Alves  _et al._ (2015a) A. Alves, A. Berlin, S. Profumo, and F. S. Queiroz, Phys. Rev. D 92, 083004 (2015a), arXiv:1501.03490 [hep-ph] .
  * Alves  _et al._ (2015b) A. Alves, A. Berlin, S. Profumo, and F. S. Queiroz, JHEP 10, 076 (2015b), arXiv:1506.06767 [hep-ph] .
  * Dorsch  _et al._ (2017) G. C. Dorsch, S. J. Huber, K. Mimasu, and J. M. No, JHEP 12, 086 (2017), arXiv:1705.09186 [hep-ph] .
  * Schmitz (2021) K. Schmitz, JHEP 01, 097 (2021), arXiv:2002.04615 [hep-ph] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Chanda  _et al._ (2020) P. Chanda, S. Hamdan, and J. Unwin, Journal of Cosmology and Astroparticle Physics 2020, 034–034 (2020).
  * Kolb and Turner (1990) E. W. Kolb and M. S. Turner,  _The Early Universe_, Vol. 69 (1990).
  * Bernal  _et al._ (2017) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, Int. J. Mod. Phys. A 32, 1730023 (2017), arXiv:1706.07442 [hep-ph] .
  * Arakawa and Tait (2021) J. Arakawa and T. M. P. Tait, SciPost Phys. 11, 019 (2021), arXiv:2101.11031 [hep-ph] .
  * Roszkowski  _et al._ (2018) L. Roszkowski, E. M. Sessolo, and S. Trojanowski, Rept. Prog. Phys. 81, 066201 (2018), arXiv:1707.06277 [hep-ph] .
  * Laurent and Cline (2020) B. Laurent and J. M. Cline, Phys. Rev. D 102, 063516 (2020), arXiv:2007.10935 [hep-ph] .


###### Abstract

Complex structures can only form in a universe that allows for bound states.
While this is clearly observed in three-dimensions, added degrees of freedom
in a higher-dimensional space preclude the immediate assumption that binding
potentials can in fact exist. In this paper, we derive a constraint on the
dimensionality of a universe in the presence of an arbitrary set of forces. We
then apply this constraint to systems with several example potentials. In
doing so, we find that bound states in higher than 3 dimensions are in fact
possible under specific circumstances which we characterize.

Constraints on the Dimensionality of Space

C. A. Petway, R. D. Orlando, A. M. McNamara, E. A. Zweig,

B. C. Caminada, E. J. Kincaid, C. V. Landgraf,

C. M. Mohs, M. L. Schiff, E. Fischbach111Corresponding Author:
<EMAIL_ADDRESS>
Department of Physics and Astronomy, Purdue University

West Lafayette, Indiana, 47907, USA

Keywords: Dimensionality of Space, New Forces, Bound States

##  1 Introduction

In order to have a universe of atoms and stars, forces need to allow for bound
states to occur. Without the attraction supplied by fundamental forces, the
early universe would never have reached baryogenesis, and the universe would
remain a quickly cooling sea of unbound particles. However, attractive forces
do exist in our universe, allowing for quarks to form baryons, and eventually
for atoms, stars, and galaxies to form.

The local observable universe has three spatial dimensions. On a larger scale,
however, the larger universe could exhibit four or more spatial dimensions, at
least in some local regions. Mathematically, any curved surface in
$n$-dimensions on a small enough scale approaches a flat tangential surface in
$(n-1)$-dimensions; this concept is the foundation for taking derivatives.
Additionally, cosmologists have long hypothesized about the shape of the
universe, evidenced by the curvature term included in the Friedmann equation.

However, bound states must be possible in a higher-dimensional space for
complex structures to form. Since altering the dimensionality of space affects
many equations derived in our local universe, the instability provided by
increasing spatial degrees of freedom precludes the natural assumption that
bound states can also occur in regions with more than three spatial
dimensions. In Ref. [1], Rabinowitz demonstrated that the existence of
gravitational and electromagnetic forces demands that our observable universe
be only three-dimensional. Separately, L$\ddot{\text{a}}$mmerzahl and Macias
have presented general arguments why the dimensions of space-time should be
four [2]. In what follows we expand upon Rabinowitz’s work by deriving a
constraint on a higher dimensional universe from a potential with arbitrary
form. We show specifically that combinations of electromagnetic,
gravitational, and short range Yukawa forces can in fact produce bound states
in a Euclidean universe with greater than three spatial dimensions.

##  2 Deriving the Constraint

###  2.1 Hypersphere Surface Area

The $r$-dependence of typical potentials results from applying Gauss’ law over
some volume containing the source of a force field. An integral part of
generalizing this derivation to $n$ dimensions is determining the surface area
of an $n$-sphere. An easy way to do this is to first consider the volume
$\mathcal{V}_{n}$ of an $n$-dimensional Euclidean sphere ($n$-sphere). A
derivation of these volumes is rather involved but has been carried out by
Gipple in Ref. [3]. He finds that

$\mathcal{V}_{n}=\frac{2\pi^{n/2}(r_{n})^{n}}{n\Gamma(n/2)},$ (2.1.1)

where $r_{n}$ is the radius of the $n$-sphere, and $\Gamma$ is the usual gamma
function. To obtain the surface area $A_{n-1}$ of the $n$-sphere, we
differentiate the volume with respect to $r_{n}$, which yields,

$A_{n-1}=\frac{d\mathcal{V}_{n}}{dr_{n}}=\frac{2\pi^{n/2}(r_{n})^{n-1}}{\Gamma(n/2)}.$
(2.1.2)

###  2.2 General Potential

Rabinowitz provided a constraint on the number of spatial dimensions a
universe can have while maintaining the ability to have bound states [1]. In
his paper, he used the Virial Theorem for a Coulomb potential and found that
bound orbits are only possible in 3-dimensional universes. If, however, one
assumes that potentials can take on many different functional forms, other
solutions appear. Consider to start with a 3-dimensional Coulomb or
gravitational potential given by

$V_{3}=-\frac{K/q}{4\pi r_{3}},$ (2.2.1)

where $K=Ze^{2}/\epsilon_{0}$ for the Coulomb force, or $4\pi GMm$ for the
gravitational force, and $q$ is the generalized charge of the object the
potential is acting on. Specifically, $q=e$ for the electromagnetic force, or
$q=m$ for the gravitational force. The force is then given by the negative
gradient of such a potential when multiplied by the generalized charge $q$,

$\vec{F}_{3}=(-\vec{\nabla}V_{3})q=-\frac{K}{4\pi(r_{3})^{2}}\hat{r}_{3}.$
(2.2.2)

To generalize this result to an arbitrary functional form in higher Euclidean
spatial dimensions, we start by multiplying (2.2.1) by an appropriate function
of interest, $\phi(r_{n})$, and take the gradient of this new potential.
Noting that $A_{2}=4\pi(r_{3})^{2}$ from (2.1.2), we find

$V_{3}=-\frac{K/q}{4\pi r_{3}}\phi(r_{3}),$ (2.2.3)
$\vec{F}_{3}=(-\vec{\nabla}V_{3})q=-\frac{K}{A_{2}}\left[\phi(r_{3})-r_{3}\phi^{\prime}(r_{3})\right]\hat{r_{3}},$
(2.2.4)

where $\phi^{\prime}$ denotes $d\phi(r_{3})/dr_{3}$. To check that this
maintains the correct relationship between the force and the potential, we can
take $\phi(r_{3})=1$ in (2.2.4), which evidently reproduces (2.2.2) as
expected. By multiplying both sides of (2.2.4) by $A_{2}$, we recover Gauss’
law,

$\oint_{A_{2}}\vec{F}_{3}\cdot
d\vec{A}_{2}=-K\left[\phi(r_{3})-r_{3}\phi^{\prime}(r_{3})\right].$ (2.2.5)

###  2.3 $n$-Dimensional Forces

Since the expression for the surface area of a sphere in any dimension is
known, we can extend the relation in (2.2.5) to an arbitrary dimension,

$\oint_{A_{n-1}}\vec{F}_{n}\cdot
d\vec{A}_{n-1}=-K\left[\phi(r_{n})-r_{n}\phi^{\prime}(r_{n})\right].$ (2.3.1)

The force $\vec{F}_{n}$ will always be normal to the Gaussian surface in a
Euclidean space, and hence the right hand side of (2.3.1) simplifies to the
product of the force and the surface area. Solving for $\vec{F}_{n}$ leads to
the generalized force in an arbitrary dimension,

$\vec{F}_{n}=-\frac{K}{A_{n-1}}\left[\phi(r_{n})-r_{n}\phi^{\prime}(r_{n})\right]=-\frac{K\Gamma(n/2)}{2\pi^{n/2}r_{n}^{n-1}}[\phi(r_{n})-r_{n}\phi^{\prime}(r_{n})]\hat{r}_{n}.$
(2.3.2)

###  2.4 Existence of Bound States in Higher Dimensions

To simplify the problem, we focus on the case of circular orbits, since for
any elliptical orbit there will be a circular orbit of the same energy. For an
object to be bound while in orbit, its total energy must be negative,

$E_{n}=T_{n}+U_{n}<0,$ (2.4.1)

where $T_{n}$ is the kinetic energy and $U_{n}$ is the potential energy. For a
force $\vec{F}_{n}$ to maintain an object in a circular orbit, the magnitude
of the force $\absolutevalue{\vec{F}_{n}}$ must be given by
$m(v_{n})^{2}/r_{n}$. Combining this with the kinetic energy
$T_{n}=\frac{1}{2}m(v_{n})^{2}$, it follows that

$T_{n}=\frac{1}{2}\absolutevalue{F_{n}}r_{n}=\absolutevalue{\frac{K\Gamma(n/2)}{4\pi^{n/2}(r_{n})^{n-2}}\left[\phi(r_{n})-r_{n}\phi^{\prime}(r_{n})\right]}.$
(2.4.2)

The potential $U_{n}$ is given by the path integral of the force,

$U_{n}=-\int_{\infty}^{r_{n}}\vec{F}_{n}\cdot
d(\vec{r}_{n})^{\prime}=\frac{K\Gamma(n/2)}{2\pi^{n/2}}\left[\int_{\infty}^{r_{n}}\frac{\phi(x)}{x^{n-1}}dx-\int_{\infty}^{r_{n}}\frac{\phi^{\prime}(x)}{x^{n-2}}dx\right].$
(2.4.3)

Using (2.4.2) and (2.4.3), a constraint on $n$ and $\phi(r_{n})$ can be found
by looking for solutions which satisfy (2.4.1),

$\frac{1}{2}\absolutevalue{\frac{\phi(r_{n})}{(r_{n})^{n-2}}-\frac{\phi^{\prime}(r_{n})}{(r_{n})^{n-3}}}+\int_{\infty}^{r_{n}}\frac{\phi(x)}{x^{n-1}}dx-\int_{\infty}^{r_{n}}\frac{\phi^{\prime}(x)}{x^{n-2}}dx<0.$
(2.4.4)

For the force to be net attractive, the expression within the absolute value
must be positive, for if it were negative, the force would be repulsive and
bound states would not be possible. This leads to the following set of
constraints on the forces,

$\phi(r_{n})\geq r_{n}\phi^{\prime}(r_{n}),$ (2.4.5)
$\frac{1}{2}\left(\frac{\phi(r_{n})}{(r_{n})^{n-2}}-\frac{\phi^{\prime}(r_{n})}{(r_{n})^{n-2}}\right)+\int_{\infty}^{r_{n}}\frac{\phi(x)}{x^{n-1}}dx-\int_{\infty}^{r_{n}}\frac{\phi^{\prime}(x)}{x^{n-2}}dx<0.$
(2.4.6)

##  3 Applying the Constraint

Here we consider several different expressions for $\phi(r_{n})$ for which we
derive the specific constraint.

###  3.1 Coulomb Potential

We begin with a Coulomb potential which characterizes both the gravitational
and electromagnetic forces. Although these forces are typically referred to as
inverse square law forces, this name loses its descriptive meaning when
generalized to $n$ dimensions. This potential is a convenient starting point
because this case has already been analyzed by Rabinowitz in Ref. [1]. In the
notation of (2.2.3) and (2.2.4), $\phi(r_{n})$ and its derivative are then
given by

$\phi(r_{n})=1$ (3.1.1) $\phi^{\prime}(r_{n})=0.$ (3.1.2)

Clearly these functions satisfy (2.4.5). Combining these two values and
(2.4.6) gives the following inequality,

$0>\frac{1}{2}\left(\frac{1}{r^{n-2}_{n}}\right)+\int_{\infty}^{r_{n}}\frac{1}{x^{n-1}}dx.$
(3.1.3)

It follows that the only physically meaningful value of $n$ consistent with
(3.1.3) is $n=3$ as shown by Rabinowitz.

###  3.2 Yukawa Potential

The strong nuclear force can be described by a Yukawa potential, i.e.,
$V=-g^{2}e^{-r/\lambda}/r$ where $\mu$ is a constant. This can be analyzed in
our framework by assuming the following expression for $\phi(r_{n})$,

$\phi(r_{n})=e^{-r_{n}/\lambda},$ (3.2.1)
$\phi^{\prime}(r_{n})=-\frac{e^{-r_{n}/\lambda}}{\lambda}.$ (3.2.2)

Applying (2.4.5) to these functions gives $1\geq-r_{n}/\lambda$. Since $r_{n}$
and $\lambda$ are positive by definition, this condition is always satisfied.
Thus the constraint from (2.4.6) becomes

$0>\frac{1}{2}\left(\frac{e^{-r_{n}/\lambda}}{r_{n}^{n-2}}+\frac{e^{-r_{n}/\lambda}}{\lambda
r^{n-3}_{n}}\right)+\int_{\infty}^{r_{n}}\frac{e^{-x/\lambda}}{x^{n-1}}dx+\int_{\infty}^{r_{n}}\frac{e^{-x/\lambda}}{\lambda
x^{n-2}}dx.$ (3.2.3)

The contribution from $e^{-r_{n}/\lambda}/r_{n}^{n-2}$ can be factored out of
the kinetic energy term, and $-1/\lambda^{n-2}$ can be factored out of the
potential energy term, while inverting the bounds of the integral. Doing so
leads to

$0>\frac{e^{-r_{n}/\lambda}}{2r_{n}^{n-2}}\left(1+\frac{r_{n}}{\lambda}\right)-\frac{1}{\lambda^{n-2}}\left[\int_{r_{n}}^{\infty}\frac{\lambda^{n-2}e^{-x/\lambda}}{x^{n-1}}dx+\int_{r_{n}}^{\infty}\frac{\lambda^{n-3}e^{-x/\lambda}}{x^{n-2}}dx\right].$
(3.2.4)

We next substitute $u=x/\lambda$ and $du=dx/\lambda$, while multiplying all
terms by $2r_{n}^{n-2}$. This leads to the inequality,

$0>e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)-2\left(\frac{r_{n}}{\lambda}\right)^{n-2}\left[\int_{r_{n}/\lambda}^{\infty}u^{1-n}e^{-u}du+\int_{r_{n}/\lambda}^{\infty}u^{2-n}e^{-u}du\right].$
(3.2.5)

Each of the integrals can be expressed in terms of an upper incomplete gamma
function, $\Gamma(a,b)$, defined as

$\Gamma(a,b)=\int_{b}^{\infty}x^{a-1}e^{-x}dx.$ (3.2.6)

Combining (3.2.6) and (3.2.5) gives

$0>e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)-2\left(\frac{r_{n}}{\lambda}\right)^{n-2}\left[\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)+\Gamma\left(3-n,\frac{r_{n}}{\lambda}\right)\right].$
(3.2.7)

Upper incomplete gamma functions have the following property, which is useful
at this point,

$\Gamma(a+1,b)=a\Gamma(a,b)+b^{a}e^{-b}.$ (3.2.8)

Combining (3.2.7) and (3.2.8) leads to

$0>e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)-2\left(\frac{r_{n}}{\lambda}\right)^{n-2}\left[\left(\frac{r_{n}}{\lambda}\right)^{2-n}e^{-r_{n}/\lambda}+(3-n)\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)\right].$
(3.2.9)

Simplifying this further yields the following constraint,

$0>e^{-r_{n}/\lambda}\left(\frac{r_{n}}{\lambda}-1\right)-(6-2n)\left(\frac{r_{n}}{\lambda}\right)^{n-2}\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right).$
(3.2.10)

This constraint is graphically represented in Fig. 1. In this figure it is
clear that a Yukawa force can only bind particles in three or fewer
dimensions, and that there is a maximum range in which these states can occur
which is $\lambda$ in three dimensions and a multiple of lambda in one and two
dimensions.

Figure 1: Plot of the maximum dimension (shaded region) allowing bound states
to exist for a Yukawa potential. Bound states can only arise in 3 or fewer
dimensions, and then only for separations shorter than the range $\lambda$ of
the force in three dimensions and some multiple of $\lambda$ for one and two
dimensions. The region where bound states cannot arise and any dimensionality
is labeled ”Unbound.”

###  3.3 Coulomb$(1+\text{constant}/r)$ Potential

We can also consider potentials $\phi(r_{n})$ that have more than one term. In
this case, one of the forces decreases with distance faster than the other by
a factor of $1/(r_{n})$. There are simple expressions for $\phi(r_{n})$ and
$\phi^{\prime}(r_{n})$ which we can use to explore this possibility:

$\phi(r_{n})=1+\frac{\alpha}{r_{n}},$ (3.3.1)
$\phi^{\prime}(r_{n})=-\frac{\alpha}{r_{n}^{2}}.$ (3.3.2)

Here $\alpha$ denotes the ratio of the strength of the faster decaying force
to the slower one at $r_{n}=1$. Using (2.4.5), we can see that there is a
lower bound on the value of $\alpha$,

$1+\frac{\alpha}{r_{n}}\geq r_{n}(-\frac{\alpha}{r_{n}^{2}}).$ (3.3.3)

It follows from (3.3.3) that the lower bound on $\alpha$ is given by

$\alpha\geq-r_{n}/2.$ (3.3.4)

Inserting this value of $\alpha$ into $\phi(r_{n})$ and
$\phi^{\prime}(r_{n})$, and combining with (2.4.6) leads to the following
constraint on $n$,

$0>\frac{1}{2}\left(\frac{1+\alpha/r_{n}}{r_{n}^{n-2}}+\frac{\alpha/r_{n}^{2}}{r_{n}^{n-3}}\right)+\int_{\infty}^{r_{n}}\frac{1+\alpha/x}{x^{n-1}}dx+\int_{\infty}^{r_{n}}\frac{\alpha/x^{2}}{x^{n-2}}dx.$
(3.3.5)

By separating the first integral into two contributions, we note that the
second has the same functional form as the last integral. Combining the
various terms in (3.3.5) we find,

$0>\frac{1}{2}\left(\frac{1}{r_{n}^{n-2}}+\frac{\alpha}{r_{n}^{n-1}}+\frac{\alpha}{r_{n}^{n-1}}\right)+\int_{\infty}^{r_{n}}\frac{1}{x^{n-1}}dx+2\int_{\infty}^{r_{n}}\frac{\alpha}{x^{n}}dx.$
(3.3.6)

Assuming again that $n>1$ we can perform these integrations which gives,

$0>\frac{1}{2}\left(\frac{1}{r_{n}^{n-2}}+\frac{2\alpha}{r_{n}^{n-1}}\right)-\frac{1}{(n-2)r_{n}^{n-2}}-\frac{2\alpha}{(n-1)r_{n}^{n-1}}.$
(3.3.7)

After multiplying through by $2/r_{n}^{n-2}$, we find

$0>1+\frac{2\alpha}{r_{n}}-\frac{2}{n-2}-\frac{4\alpha}{(n-1)r_{n}},$ (3.3.8)

and hence,

$0>\frac{n-4}{n-2}+\frac{2\alpha}{r_{n}}\left(\frac{n-3}{n-1}\right).$ (3.3.9)

The application of both equations (3.3.4) and (3.3.9) is graphically
represented in Figure 2.

Figure 2: Plot of the maximum dimension allowing for bound states of a system
with both a Coulomb potential and a potential that decreases faster by a
factor of $r^{-1}$. When the forces oppose each other, there exists a small
band of distances in which there are bound states in higher spatial
dimensions. For computational reasons, six spatial dimensions were the highest
tested. However, bound states do exist in higher dimensions. The region where
bound states cannot arise and any dimensionality is labeled ”Unbound.”

It can be seen in this figure that the maximum dimensionality allowing for
bound states increases as the parameters approach the lower bound of $\alpha$.
If we take $\alpha$ to be at its minimum possible value for a given $r_{n}$,
that is $\alpha=-r_{n}/2$, and apply it to (3.3.9), we find

$\frac{n-3}{n-1}>\frac{n-4}{n-2}.$ (3.3.10)

Since (3.3.10) is equivalent to the statement that $6>4$, it follows that at
this minimum value of $\alpha$, any number of spatial dimensions are allowed.

###  3.4 Coulomb plus Yukawa Potential

We next consider a more physical form of a potential with two terms, namely
the sum of a Yukawa and a Coulomb potential. In fact, we have actual examples
of particles that interact via forces of both kinds, such as quarks. The
$\phi(r_{n})$ and $\phi^{\prime}(r_{n})$ are then as follows,

$\phi(r_{n})=1+\Omega e^{-r_{n}/\lambda},$ (3.4.1)
$\phi^{\prime}(r_{n})=-\frac{\Omega}{\lambda}e^{-r_{n}/\lambda}.$ (3.4.2)

Here $\Omega$ is the ratio of the strengths of the Yukawa and Coulomb
potentials, and $\lambda$ is the range of the Yukawa force. Again we have a
constraint on $\Omega$ from (2.4.5),

$1+\Omega e^{-r_{n}/\lambda}\geq
r_{n}\left(-\frac{\Omega}{\lambda}e^{-r_{n}/\lambda}\right),$ (3.4.3)

which is equivalent to

$\Omega e^{-r_{n}/\lambda}(1+\frac{r_{n}}{\lambda})\geq-1.$ (3.4.4)

It follows that (3.4.4) implies a lower bound the value of $\Omega$ given by

$\Omega\geq-\frac{e^{r_{n}/\lambda}}{1+r_{n}/\lambda}.$ (3.4.5)

A constraint on $n$ from (2.4.6) can be obtained by using the expression for
$\phi(r_{n})$ and $\phi^{\prime}(r_{n})$ for this potential:

$0>\frac{1}{2}\left(\frac{1+\Omega
e^{-r_{n}/\lambda}}{(r_{n})^{n-2}}+\frac{\Omega
e^{-r_{n}/\lambda}}{\lambda(r_{n})^{n-3}}\right)+\int_{\infty}^{r_{n}}\frac{1+\Omega
e^{-x/\lambda}}{x^{n-1}}dx+\int_{\infty}^{r_{n}}\frac{\Omega
e^{-x/\lambda}}{\lambda x^{n-2}}dx.$ (3.4.6)

The first of the integrals can be separated into two contributions, the first
of which is simple to evaluate. After doing so, $\Omega/\lambda^{n-2}$ can be
factored out of the remaining integrals and $r_{n}^{2-n}/2$ can be factored
out of the kinetic energy term. The inequality in (3.4.6) then assumes the
form

$\begin{split}0>&\frac{r_{n}^{2-n}}{2}\bigg{[}1+\Omega
e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)\bigg{]}-\frac{r_{n}^{2-n}}{n-2}+\frac{\Omega}{\lambda^{n-2}}\bigg{[}\int_{\infty}^{r_{n}}\frac{\lambda^{n-2}e^{-x/\lambda}}{x^{n-1}}dx\\\
&+\int_{\infty}^{r_{n}}\frac{\lambda^{n-3}e^{-x/\lambda}}{x^{n-2}}dx\bigg{]}.\end{split}$
(3.4.7)

Substituting $u=x/\lambda$ and $du=dx/\lambda$, we find

$\begin{split}0>&1-\frac{2}{n-2}+\Omega
e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)+2\Omega\left(\frac{r_{n}}{\lambda}\right)^{n-2}\bigg{[}\int_{\infty}^{r_{n}/\lambda}u^{1-n}e^{-u}du\\\
&+\int_{\infty}^{r_{n}/\lambda}u^{2-n}e^{-u}du\bigg{]}.\end{split}$ (3.4.8)

The first two terms of the the previous equation simplify to the expression
for the constraint on the Coulomb potential given in (3.1.3). Again these
integrals are of the same form as the upper incomplete gamma function given in
(3.2.6), and hence,

$\begin{split}0>&\frac{n-4}{n-2}+\Omega
e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)-2\Omega\left(\frac{r_{n}}{\lambda}\right)^{n-2}\bigg{[}\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)\\\
&+\Gamma\left(3-n,\frac{r_{n}}{\lambda}\right)\bigg{]}.\end{split}$ (3.4.9)

We can also use the property of incomplete gamma functions given in (3.2.8) to
derive the following inequality:

$\begin{split}0>&\frac{n-4}{n-2}+\Omega
e^{-r_{n}/\lambda}\left(1+\frac{r_{n}}{\lambda}\right)-2\Omega\left(\frac{r_{n}}{\lambda}\right)^{n-2}\bigg{[}\left(\frac{r_{n}}{\lambda}\right)^{2-n}e^{-r_{n}/\lambda}\\\
&+(3-n)\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)\bigg{]}.\end{split}$
(3.4.10)

The inequality in (3.4.10) simplifies to the following constraint,

$0>\frac{n-4}{n-2}+\Omega
e^{-r_{n}/\lambda}\left(\frac{r_{n}}{\lambda}-1\right)-\Omega(6-2n)\left(\frac{r_{n}}{\lambda}\right)^{n-2}\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right).$
(3.4.11)

Unlike the case for the Coulomb potential, this constraint allows for bound
states in higher than 3 dimensions but only for specific values of $\Omega$
and $\lambda$. This is graphically represented in Figure 3.

Figure 3: Plot of the maximum dimension allowed for bound states of a system
in the presence of the sum of a Coulomb potential and a Yukawa potential. We
find two unbounded regions, where the upper one (labeled ”UB”) is due to the
the finite range of the Yukawa potential, and the lower (labeled ”Unbound”) is
due to the net force being repulsive. As the parameters approach this lower
unbounded region, the maximum dimension allowing for bound states
significantly increases. For computational reasons, six spatial dimensions
were the highest tested. However, bound states do exist in higher dimensions.

We see in this plot similar behavior to the case of the
Coulomb($1+\text{constant}/r)$ potential in that the maximum dimensionality
increases as the parameters approach the lower bound. It thus appears that we
might also see the possibility of near infinite dimensions at this minimum
$\Omega$. Noting that the kinetic energy term vanishes at the minimum
$\Omega$, we can derive the following constraint from (2.4.6), using
$\phi(r_{n})$, and $\phi^{\prime}(r_{n})$ for this potential:

$0>\int_{\infty}^{r_{n}}\frac{1+\Omega
e^{-x/\lambda}}{x^{n-1}}dx+\int_{\infty}^{r_{n}}\frac{\Omega
e^{-x/\lambda}}{\lambda x^{n-2}}dx.$ (3.4.12)

As before, we can separate the first integral into two, the first of which can
be easily evaluated. After removing a factor an $\Omega/\lambda^{n-2}$ from
the remaining integrals, we find

$0>-\frac{r_{n}^{2-n}}{n-2}+\frac{\Omega}{\lambda^{n-2}}\left[\int_{\infty}^{r_{n}}\frac{\lambda^{n-2}e^{-x/\lambda}}{x^{n-1}}dx+\int_{\infty}^{r_{n}}\frac{\lambda^{n-3}e^{-x/\lambda}}{x^{n-2}}dx\right].$
(3.4.13)

As we have done previously, we can express the the two remaining integrals as
an upper incomplete gamma function, leading to,

$0>-\frac{r_{n}^{2-n}}{n-2}-\frac{\Omega}{\lambda^{n-2}}\left[\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)+\Gamma\left(3-n,\frac{r_{n}}{\lambda}\right)\right].$
(3.4.14)

The expression in (3.4.14) can be simplified using (3.2.8) yielding,

$0>-\frac{1}{n-2}-\Omega\left(\frac{r_{n}}{\lambda}\right)^{n-2}\left[e^{-r_{n}/\lambda}\left(\frac{r_{n}}{\lambda}\right)^{2-n}+(3-n)\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)\right].$
(3.4.15)

Simplifying further gives an expression similar to the final constraint found
before,

$0>-\frac{1}{n-2}-\Omega
e^{-r_{n}/\lambda}-(3-n)\Omega\left(\frac{r_{n}}{\lambda}\right)^{n-2}\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right).$
(3.4.16)

Substituting in the minimum value of $\Omega$ given in (3.4.5) we find,

$0>-\frac{1}{n-2}+\frac{1}{1+r_{n}/\lambda}+\frac{3-n}{1+r_{n}/\lambda}\left(\frac{r_{n}}{\lambda}\right)^{n-2}e^{r_{n}/\lambda}\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right).$
(3.4.17)

Next we can multiply through by $1/(1+r_{n}/\lambda)$ which gives,

$0>-\frac{1+r_{n}/\lambda}{n-2}+1+(3-n)\left(\frac{r_{n}}{\lambda}\right)^{n-2}e^{r_{n}/\lambda}\Gamma\left(2-n,\frac{r_{n}}{\lambda}\right)\\\
.$ (3.4.18)

Let us now consider the limit of this expression as $n\rightarrow\infty$. The
first term vanishes for all finite $r_{n}$. For simplicity we substitute
$x=r_{n}/\lambda$, which gives the following constraint for bound states to
exist in any number of dimensions,

$0>\lim_{n\rightarrow\infty}\left[1+(3-n)x^{n-2}e^{x}\Gamma(2-n,x)\right].$
(3.4.19)

The expression in (3.4.19) can be simplified by expressing the inequality in
the form

$1<\lim_{n\rightarrow\infty}\left[nx^{n-2}e^{x}\Gamma(2-n,x)\right].$ (3.4.20)

If the gamma function is then re-expressed in an integral form this becomes,

$1<\lim_{n\rightarrow\infty}\left[nx^{n-2}e^{x}\int_{x}^{\infty}u^{1-n}e^{-u}du\right].$
(3.4.21)

Although we will now show that this constraint is not satisfied as
$n\rightarrow\infty$, the limit of this expression does yield an interesting
result. Integrating by parts the integral in (3.4.21) gives

$\begin{split}\int_{x}^{\infty}y^{1-n}e^{-y}dy&=\evaluated{\left(\frac{y^{2-n}}{2-n}e^{-y}\right)}_{x}^{\infty}+\frac{1}{2-n}\int_{x}^{\infty}{y^{2-n}e^{-y}dy}\\\
&=-\frac{x^{2-n}}{2-n}e^{-x}+\
\frac{1}{2-n}\int_{x}^{\infty}{y^{2-n}e^{-y}dy}.\end{split}$ (3.4.22)

Continuing to integrate by parts yields

$\int_{x}^{\infty}{y^{1-n}e^{-y}dy}=-\frac{x^{2-n}}{2-n}e^{-x}\left(1+\frac{x}{3-n}+\frac{x^{2}}{\left(3-n\right)\left(4-n\right)}+\ldots\right).$
(3.4.23)

Taking the limit as $n\rightarrow\infty$ in (3.4.20), and retaining only the
leading term in the expanded integral shown in (3.4.23), we find

$\begin{split}\lim_{n\rightarrow\infty}{\frac{\Gamma\left(2-n,x\right)}{\frac{1}{n}x^{2-n}e^{-x}}}&=\left.\lim_{n\rightarrow\infty}{-\frac{x^{2-n}}{2-n}e^{-x}\left(1+\frac{x}{3-n}+\frac{x^{2}}{\left(3-n\right)\left(4-n\right)}+\ldots\right)}\middle/{\left(\frac{1}{n}x^{2-n}e^{-x}\right)}\right.\\\
&=\left.\lim_{n\rightarrow\infty}{\left(\frac{x^{2-n}}{n-2}e^{-x}\right)}\middle/{\left(\frac{1}{n}x^{\left(2-n\right)}e^{-x}\right)}\right.=1.\end{split}$
(3.4.24)

It follows that the right hand side of (3.4.20) approaches 1 as
$n\rightarrow\infty$, at which point (3.4.20) is no longer satisfied. However,
we observe bound states in increasingly high dimensions just before this
constraint is violated, which implies that there should be bound states in all
finite dimensions as the parameters approach this lower bound. As a specific
example, we find that the constraint for $\Omega=-2002.405981$ and
$r_{n}/\lambda=10$ is $n\leq 324350$.

##  4 Conclusions

For the building blocks of our universe to do more than simply scatter off one
another, there must be situations in which it is energetically advantageous
for these to couple, and to thus form more complex structures. The constraints
derived in this paper provide a simplified tool for analyzing the Euclidean
spatial dimensions where such behavior can arise. The expected results for the
Coulomb and gravitational interactions, as well as the Yukawa potential, can
be derived from these constraints. When applied to situations where multiple
forces act on the same objects, these constraints allow for more spatial
dimensions than the 3 spatial dimensions one would normally expect. The
utility of this constraint in determining the allowed dimensionality of our
universe could have far ranging implications from elementary particle physics
to cosmology.

Although we have limited our discussion to regions of our universe which can
be characterized, at least approximately, as Euclidean, it appears possible to
extend these arguments to General Relativity as will will discuss elsewhere.

## Acknowledgements

The authors wish to thank Megan McDuffie and Ben Hayward for their
contributions in the early stages of this project.

## References

[1] Mario Rabinowitz. Why Observable Space is Solely Three Dimensional.
Advanced Studies in Theoretical Physics, 8 (2014), 689-700.

[2] Claus L$\ddot{\text{a}}$mmerzahl and Alfredo Macias. On the Dimensionality
of Space-time. Journal of Mathematical Physics, 34 (1993), 4540-4553.

[3] Jake Gipple. The Volume of n-Balls. Rose-Hulman Undergraduate Mathematics
Journal, 15 (2014), 238-248.