# Domain walls in light of Cosmic Microwave Background and Pulsar Timing Array
data

D. Grüber<EMAIL_ADDRESS>Departamento de Física e Astronomia,
Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre 687,
PT4169-007 Porto, Portugal  Instituto de Astrofísica e Ciências do Espaço,
Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal
P.P. Avelino<EMAIL_ADDRESS>Departamento de Física e Astronomia,
Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre 687,
PT4169-007 Porto, Portugal  Instituto de Astrofísica e Ciências do Espaço,
Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire
Lagrange, France  L. Sousa<EMAIL_ADDRESS>Instituto de Astrofísica e
Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762
Porto, Portugal  Departamento de Física e Astronomia, Faculdade de Ciências,
Universidade do Porto, Rua do Campo Alegre 687, PT4169-007 Porto, Portugal

###### Abstract

In this paper, we study the compatibility of biased domain wall scenarios with
current gravitational wave data. We show that the Cosmic Microwave Background
bounds on the fractional density of gravitational waves at the time of
decoupling may only slightly improve on the constraints that result from
requiring that domain walls never dominate the cosmic energy budget. We show
that, despite this, the range of energy scales of the domain-wall forming
phase transitions are already quite constricted, even if the networks decay
early in cosmological history. We also show that, if domain walls are to
provide an explanation to the stochastic gravitational wave background that
was recently detected by pulsar timing arrays, they not only have to decay
early in the radiation dominated era but also their energy density would have
to be close to dominating the energy density of the universe, which would
require some fine tuning of the parameters of the models.

##  I Introduction

The recent confirmation of the detection of a stochastic gravitational wave
background (SGWB) by Pulsar Timing Arrays (PTAs) [1, 2, 3, 4] has ignited
speculations about its origin. While astrophysical phenomena are considered a
plausible explanation for the detected signals [1], the possibility of a
cosmological source remains open. Among those, Gravitational Waves (GWs) from
inflation, phase transitions in the early universe, or topological defect
networks, such as cosmic strings, are frequently explored in the literature
[5].

Like other cosmic defects, domain walls may arise as a consequence of the
spontaneous breaking of a symmetry associated with a scalar field [6]. Domain
walls, however, if they are to survive until the present day, are tightly
constrained by the Zel’dovich bound, which essentially limits the energy scale
of the wall-forming phase transition to be below approximately 1 MeV [7].
Domain wall networks, however, may not be inherently stable. Certain models,
such as those linked to axion-like particles [8, 9, 10, 11, 12] — which are
receiving increasing interest in the literature — suggest that domain wall
networks may have decayed in the past due to a bias in their scalar field
potential. In this case, domain wall networks would evade this stringent bound
and could be formed much earlier in cosmological history. Currently,
constraints on stable domain wall network scenarios are mainly based on the
Cosmic Microwave Background (CMB) [13], but they should also emit
gravitational radiation throughout their evolution. Their GW signatures should
then allow us to constrain particle physics scenarios.
A framework was recently developed in [14], where the authors developed a
semi-analytical model to characterize the SGWB produced by domain walls. In
the present work, we will build upon this work to constraint the energy scale
at domain wall formation as a function of the time of their decay, by
confronting this model with observational constraints on the SGWB at the time
of decoupling from the CMB. We also explore the possibility of domain walls
providing a major contribution to the SGWB that was recently detected by PTAs
and derive the region of parameter space that is consistent with these
observations.
In section II we will briefly recap the framework presented in [14] that
describes the SGWB by domain walls. In section III we will then show how the
constraint on the number of relativistic particle species at the time of
decoupling can be used to constraint the SGWB from domain wall networks and
introduce some further, more general constraints. We will apply the
aforementioned framework in section IV to the recent PTA signal to investigate
the possibility of the latter stemming from biased domain walls. Then we
conclude in section V.
Throughout this paper, we will work in natural units, where $c=\hbar=1$, with
$c$ being the speed of light in vacuum and $\hbar$ being the reduced Planck
constant. In these units Newton’s gravitational constant is given by
$G=6.70711\cdot 10^{-57}\,\text{eV}^{-2}$. Moreover, we will use the
cosmological parameters measured by the Planck mission [15], where the values
of the density parameters for radiation, matter and dark energy at the present
time are respectively given by $\Omega_{\rm r}=9.1476\cdot 10^{-5}$,
$\Omega_{\rm m}=0.308$, $\Omega_{\Lambda}=1-\Omega_{\rm r}-\Omega_{\rm m}$ and
the Hubble parameter is $H_{0}=2.13\cdot h\cdot 10^{-33}\,\rm eV$, with
$h=0.678$.

##  II SGWB generated by domain wall networks: an overview

In this section, we provide a brief overview of the stochastic gravitational
wave background generated by domain walls, which was studied in detail in
[14]. This spectrum is sourced by the continuous collapse of domain walls that
detach from the network as a result of the expansion of the background.
However, since the gravitational wave emission in the late stages of domain
wall collapse is not yet well understood (as numerical simulations are unable
to resolve the domain walls in the ultrarelativistic limit), some details of
the power spectrum are uncertain. In fact, in [14], the SGWB spectrum depends
on two unknown parameters. The total energy density of GWs emitted, however,
may be computed with minimal assumptions and, for this reason, it will serve
as the basis for our analysis in this paper.

###  II.1 GW energy density

In [14] the fractional contribution of the SGWB produced by standard and
biased domain wall networks to the cosmic energy budget has been estimated,
assuming that a constant fraction $\mathcal{F}$ of the energy of collapsing
domain walls is emitted in the form of GWs. This has been achieved by
resorting to the Velocity-dependent One-Scale (VOS) model [16, 17, 18] to
describe the cosmological evolution of their characteristic scale
$L\equiv\sigma/\rho=\eta^{3}/\rho$ — where $\sigma$, $\eta\equiv\sigma^{1/3}$
and $\rho$ are, respectively, the domain wall tension, their characteristic
energy scale and the average energy density of the domain wall network — and
root-mean-squared velocity $\bar{v}$.

For a power-law evolution of the scale factor — so that $a\propto t^{\lambda}$
with $\lambda=1/2$ and $\lambda=2/3$ in the radiation and matter eras,
respectively — the (biased) domain wall network evolves in a linear scaling
regime of the form $L=\xi t$ and $\bar{v}=\rm const$ until its eventual decay
at a time $t_{\star}$ (corresponding to a scale factor $a_{\star}$ and a
redshift $z_{\star}=1/a_{\star}-1$). In this case, one should have that [14]

$\Omega_{\rm gw}\equiv\left.\frac{8\pi G\rho_{\rm
gw}}{3H^{2}}\right|_{t=t_{0}}=\sigma\mathcal{F}C_{1}(1+z_{\star})^{-2}\mathcal{G}_{\star}^{1/2}\,,$
(1)

for $z_{\star}\leq z_{\rm eq}$, and

$\displaystyle\Omega_{\rm gw}$ $\displaystyle=$
$\displaystyle\sigma\mathcal{F}\mathcal{G}_{\star}^{1/2}$ (2)
$\displaystyle\times$
$\displaystyle\left[C_{2}\left((1+z_{\star})^{-\frac{5}{2}}-(1+z_{\rm
eq})^{-\frac{5}{2}}\right)+C_{3}\right]\,,$

for $z_{\star}>z_{\rm eq}$. Here, $z_{\rm eq}=\Omega_{m}/\Omega_{r}-1$ is the
redshift of radiation-matter equality,

$\displaystyle C_{1}$ $\displaystyle=$ $\displaystyle\frac{16\pi
G\Omega_{r}^{\frac{1}{2}}\mathcal{A}_{r}}{3H_{0}}\,,$ (3) $\displaystyle
C_{2}$ $\displaystyle=$ $\displaystyle\frac{12\pi
G\Omega_{m}^{\frac{1}{2}}\mathcal{A}_{m}}{5H_{0}}\,,$ (4) $\displaystyle
C_{3}$ $\displaystyle=$ $\displaystyle
C_{1}\left(\frac{\Omega_{r}}{\Omega_{m}}\right)^{2}\,,$ (5)

and $\mathcal{A}_{\rm i}\equiv\tilde{c}\bar{v}_{\rm i}/\xi^{2}_{\rm i}$. The
values of the dimensionless linear scaling parameters $\xi$ and $\bar{v}$ are
respectively given by $\xi_{\rm r}=1.532$, ${\bar{v}}_{\rm r}=0.479$, and
$\xi_{\rm m}=1.625$, ${\bar{v}}_{\rm m}=0.339$ in the radiation and matter
eras for $\tilde{c}=0.5$ (an excellent approximation in both eras [19]). The
factor of $\mathcal{G}(z)$ is a small correction to the result derived in [14]
that accounts for the impact of the change in the effective number of massless
degrees of freedom on the expansion rate. Here, we follow the approach in Ref.
[20] and approximate the Hubble parameter by

$H(z)=H_{0}\left[\Omega_{\Lambda}+\Omega_{\rm
m}(1+z)^{3}+\mathcal{G}(z)\Omega_{\rm r}(1+z)^{4}\right]^{1/2}\,,$ (6)

where

$\mathcal{G}(z)=\begin{cases}1&\text{for }z<z_{\rm e^{-}e^{+}}\\\
0.83&\text{for }z_{\rm e^{-}e^{+}}<z<z_{\rm QCD}\\\ 0.39&\text{for }z>z_{\rm
QCD}\end{cases}$ (7)

describes the decreases in the number of degrees of freedom expected to occur,
in the Standard Model of Particle Physics, at the QCD phase transition and as
a result of the electron-positron annihilation (throughout this paper we shall
take $z_{\rm e^{-}e^{+}}=10^{9}$ and $z_{\rm QCD}=10^{12}$ [20]). This factor
is responsible for a suppression of the fractional contribution of the
gravitational waves emitted by domain wall networks to the cosmic energy
budget by a factor of
$\mathcal{G}_{\star}^{1/2}\equiv\mathcal{G}(z_{\star})^{1/2}$.

Notice that, throughout this paper we shall assume that the network is in a
linear scaling regime before it is rendered unstable and instantaneously
decays without producing further GWs. Including the potential contribution of
the decay of the network to the SGWB could increase its amplitude by up to an
order of magnitude [14] (but there is still some uncertainty regarding this
contribution). We will briefly comment on the impact of this assumption in the
discussion of this paper.

###  II.2 Power spectrum

The power spectrum of the SGWB generated by domain walls, characterized by

$\Omega_{\rm gw}(f)=\frac{8\pi G}{3H_{0}^{2}}\frac{d\rho_{\rm gw}}{d\log
f}\,,$ (8)

can be roughly approximated by a broken power-law shape [14, 21], with the
power being a monotonically increasing (decreasing) function of the frequency
for frequencies smaller (larger) than the peak frequency. The frequency of the
peak of the spectrum is expected to be inversely proportional to the
characteristic length of the network at the time of decay of the network, and
is thus roughly given by111Aside from $\mathcal{O}(1)$ factors associated to
the changes in the number of relativistic degrees of freedom and to the
uncertainties regarding the power spectrum [14].

$f_{\rm p}=\frac{2}{\alpha}\frac{a_{\star}}{L_{\star}}\,,$ (9)

where $L_{\star}=\xi_{\star}t_{\star}=\xi\lambda_{\star}H_{\star}^{-1}$ and
$\alpha<1$ is a parameter introduced to characterize the typical scale of GW
in terms of $L_{\star}$ (see [14] for more details). Most of the GW energy
density is, in general, concentrated within a narrow frequency range around
the peak of the spectrum and, in fact, the peak amplitude is, aside from a
factor that is generally of order unity, approximately given by Eqs. (1)-(2).

##  III Cosmic Microwave Background constraints

The domain wall network is expected to have been formed at a temperature
$T_{\rm birth}\sim\eta$. Since the network (obviously) cannot decay before it
is formed, the characteristic energy scale of a domain wall network that
survives until $z_{\star}$ necessarily has to fulfill the condition
$T_{\star}<T_{\rm birth}$. We should then have that:

$\eta>\left(\frac{45}{4\pi^{3}}\frac{H_{0}^{2}\mathcal{G}_{\star}\Omega_{\rm
r}}{Gg_{\star}}\right)^{1/4}(1+z_{\star})\,,$ (10)

where $g_{\star}=g(z_{\star})$ and [20]

$g(z)=\begin{cases}3.36&\text{for }z<z_{\rm e^{+}e^{-}}\\\ 10.75&\text{for
}z_{\rm QCD}<z<z_{\rm e^{+}e^{-}}\\\ 106.75&\text{for }z>z_{\rm
QCD}\end{cases}\,,$ (11)

describes the effective number of degrees of freedom at a redshift $z$.

Moreover, domain walls should always provide a subdominant contribution to the
cosmic energy budget. This condition results in an additional constraint on
their characteristic energy scale of the form:

$\eta<\left(\frac{3H_{\star}\lambda_{\star}\xi_{\star}}{8\pi
G}\right)^{1/3}\,.$ (12)

Here, $\lambda_{\star}=\lambda_{\rm r}$ and $\xi_{\star}=\xi_{\rm r}$ (or
$\lambda_{\star}=\lambda_{\rm m}$ and $\xi_{\star}=\xi_{\rm m}$) for
$z_{\star}$ greater (or smaller) than $z_{\rm eq}$.

Finally, the power spectrum of temperature and polarization anisotropies of
the CMB is extremely sensitive to the expansion rate at the time of
decoupling, which, for a flat universe, is fully determined by the energy
density at that time. As a result, CMB observations impose stringent
constraints on the effective number of relativistic degrees of freedom at the
time of decoupling, which provides a tight limit on the total energy density
in gravitational waves generated before that time. One recently derived
constraint that takes into account the Planck 2018 data [15] places this limit
at

$\Omega_{\rm gw}^{\rm(CMB)}h^{2}<1.7\cdot 10^{-6}$ (13)

for frequencies $f\gtrsim 10^{-15}\,\rm Hz$ [22]. There are also similar
constraints at the time of the primordial nucleosynthesis, but these apply
only at much larger redshifts ($z\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}10^{9}$).

Figure 1: Combined constraints on $\mathcal{F}^{1/3}\eta$ as a function of the
redshift $z_{\star}$ at the time of decay of the domain wall network. The
green line represents the upper limit on $\mathcal{F}^{1/3}\eta$ resulting
from the CMB bound in Eq. (13). The red lines represent the limit in which
domain walls dominate the universe’s energy budget (area above is excluded),
while the purple lines mark the boundary where the time of decay coincides
with the time of formation of the domain wall network (area below is
excluded). In the two latter cases, the dashed- and dashed dotted lines
correspond to values of $\mathcal{F}=1$ and $\mathcal{F}=10^{-6}$
respectively. The Zel’dovich bound is indicated by the black marker in the
plot. The shaded area represents the allowed region in parameter space if
$\mathcal{F}$ is treated as a free parameter.

In Fig. 1, we display the values of $\mathcal{F}^{1/3}\eta$ allowed by CMB
data (shaded area) as a function of the redshift of decay of the network. This
allowed region of the $(z_{\star},\mathcal{F}^{1/3}\eta)$-parameter space is
obtained by combining the bound obtained by requiring that $\Omega_{\rm
gw}(z_{\rm dec})\leq\Omega_{\rm gw}^{\rm(CMB)}$ (green line) — where $z_{\rm
dec}$ is the redshift at the time of decoupling — and that domain walls never
dominate the energy budget (red lines). Regarding the latter bound, we include
the cases in which $\mathcal{F}=1$ (dashed line) and $\mathcal{F}=10^{-6}$
(dot-dashed line). For $\mathcal{F}=1$, this actually corresponds to an
absolute constraint and it is only slightly less stringent than the upper
bound on $\mathcal{F}^{1/3}\eta$ obtained from the CMB bound on $\Omega_{\rm
gw}$ for $z<z_{\rm dec}$. However, since the CMB constrains only the total
energy density of GWs created until decoupling, the resulting bounds for walls
that survive past decoupling are independent of $z_{\star}$ and weaker than
those that result from imposing that domain walls are always subdominant for
$\mathcal{F}=1$. For values of $\mathcal{F}\lesssim 0.3$, Eq. (12) is always
more constraining than the CMB bound.

One may also derive lower bounds on $\mathcal{F}^{1/3}\eta$, using Eq. (10),
but these also depend on $\mathcal{F}$. These bound are included in Fig. 1, as
dashed and dot-dashed purple lines for $\mathcal{F}=1$ and
$\mathcal{F}=10^{-6}$, respectively. Notice that, for a given value of
$\mathcal{F}$, the allowed region of parameter space is also bounded from
below by this constraint and then is given by a narrow band.

##  IV Domain wall interpretation of PTA data

Recently, different PTAs — including the North American Nano-hertz Observatory
for Gravitational waves (NANOGrav) [1], the Parks PTA (PPTA) [2], the European
PTA (EPTA) [3] and the Chinese PTA (CPTA) [4] — have produced the first
convincing evidence for a SGWB in the nanohertz frequency band. These results
all indicate the existence of an excess signal with a characteristic strain
amplitude of the order $\mathcal{O}(10^{-15})$ at a reference frequency of
$1\,{\rm yr}^{-1}\simeq 31.7\,{\rm nHz}$. Assuming that the spectral index of
the spectrum $d$ — defined such that $\Omega_{\rm gw}\propto f^{d}$ locally —
is free, this data may be translated into the following value for the spectral
density of gravitational waves:

$\Omega_{\rm gw}(f)h^{2}\simeq 6.3\times 10^{-10}A^{2}\left(\frac{f}{{\rm
yr}^{-1}}\right)^{d}\,,$ (14)

with a range of allowed parameters of $A=6.4^{+4.2}_{-2.7}$ and $d=1.8\pm 0.6$
for NANOGrav [1], $A=2.9^{+2.6}_{-1.8}$ and $d=0.81^{+0.63}_{-0.73}$ for EPTA
[3] and $A=3.1^{+1.3}_{-0.9}$ and $d=1.1\pm 0.4$ for PPTA [2]. This spectral
shape is valid within a frequency range of roughly $[f_{\rm min},f_{\rm
max}]$, with $f_{\rm min}=3.2\cdot 10^{-8}\,\rm Hz$ and $f_{\rm
max}=10^{-7}\,\rm Hz$. Here, we will investigate whether such a signature
could be produced by a biased domain wall network as was proposed in the
literature [9, 5, 11, 23, 24].

Given the uncertainty associated to the shape of the SGWB power spectrum, here
we will again focus on the total gravitational wave energy density instead. A
conservative lower bound on the fractional contribution of the SGWB to the
cosmic energy budget in the PTA frequency range may be obtained by calculating

$\Omega_{\rm gw}^{\rm(PTA)}=\int_{f_{\rm min}}^{f_{\rm max}}\Omega_{\rm
gw}(f)\frac{df}{f}\,,$ (15)

while considering the lower bounds on the parameters $A$ and $d$. This gives
us the minimum GW energy density one would need to have in this range to
generate the observed signal in each PTA experiment: $\left(\Omega_{\rm
gw}^{\rm(NANO)},\,\Omega_{\rm gw}^{\rm(EPTA)},\,\Omega_{\rm
gw}^{\rm(PPTA)}\right)=\left(4.643\cdot 10^{-8},\,1.995\cdot
10^{-9},\,1.170\cdot 10^{-8}\right)$.

The observed PTA spectrum, fitted by the function given in Eq. (14), appears
to be an increasing function of frequency, as all experiments predict that
$d>0$. Hence, if the main contribution to this spectrum was to be provided by
the gravitational radiation produced by biased domain walls, it should mostly
correspond to the part of the spectrum that lies to the left of the peak (on
the lower frequency side). Moreover, since for $f\ll f_{\rm p}$, one should
have $\Omega_{\rm gw}(f)\propto f^{3}$ due to causality constraints [25], the
peak of the spectrum cannot be located at a frequency that is much larger than
$f_{\rm max}$. On one hand, PTA data seems to be consistent with smaller
values of $d$ — which would indicate that, if domain walls are to be the main
source, the detected signal should correspond to the transitional region
between the peak and the regime characterized by $d=3$. On the other hand, if
$f_{\rm p}\gg f_{\rm max}$, the amplitude of the peak would have to be very
large in order to explain the signal within the PTA window and $\Omega_{\rm
gw}$ would then exceed the stringent bounds discussed in the previous section.
In fact, if the spectrum between $f_{\rm max}$ and $f_{\rm p}$ has a spectral
index of $d=3$, the total energy contained in the frequency interval $[f_{\rm
max},\,f_{\rm p}]$ would exceed the total energy limits dictated by the CMB
bound for $f_{\rm p}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}5\cdot 10^{-7}\,\rm Hz=5\cdot f_{\rm
max}$. The upper limit we compute for the redshift $z_{\star}$ at the time of
decay of the network, under the assumption that $f_{\rm p}$ is contained in
$[f_{\rm min},f_{\rm max}]$, may then increase at most by this factor of 5. We
will then assume that the peak of the spectrum lies within the frequency
interval $[f_{\rm min},f_{\rm max}]$.

Figure 2: Allowed values for the characteristic energy scale $\eta$ of domain
walls as a function of the redshift of decay $z_{\star}$. The red and purple
lines represent the constraints resulting from imposing that walls cannot
dominate the energy budget of the universe and that their time of decay cannot
be before the time of their formation. The area shaded in blue represents the
region of parameter space compatible with these constraints. The color shaded
band represents the region of parameter space compatible with the NANOGrav
data, assuming domain wall networks constitute the main contribution to this
signal. The dark- and light green rectangles represent the allowed regions
assuming a fixed value of $\alpha=1$ and $\alpha=10^{-6}$ respectively. More
details in the text.

If $f_{\rm p}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$>$}}f_{\rm min}$ and assuming that
$\alpha\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise
1.0pt\hbox{$>$}}10^{-8}$, then the domain wall network must have decayed
during the radiation dominated era. Eq. (9) then implies

$f_{\rm
p}=\frac{8H_{0}\Omega_{r}^{\frac{1}{2}}}{\alpha\xi_{r}}(1+z_{\star})\,.$ (16)

The redshift of decay of the network may then be inferred from the location of
the peak for each value of $\alpha$. Then, using Eq. (1) and the lower bounds
we have derived for the value of $\Omega_{\rm gw}$ in different PTA
experiments, we may derive constraints on $\mathcal{F}^{1/3}\eta$. Note that,
in deriving these constraints, we will neglect the contribution to
$\Omega_{\rm gw}$ outside of the $[f_{\rm min},f_{\rm max}]$ range and
therefore the bounds derived here should be conservative. Moreover, since the
domain wall SGWB should be simultaneously consistent with the data of all PTA
experiments, we consider $\Omega_{\rm gw}^{\rm(NANO)}=4.64\cdot 10^{-8}$,
which corresponds to the strongest signal detected across the various PTA
experiments.

In Fig. 2, we plot the allowed region of the $(z_{\star},\eta)$-parameter
space (shaded blue area) for domain wall networks. This figure includes the
constraints in Eqs. (10) and (12) — again plotted as red and purple lines
respectively — but here they become fixed boundaries of the allowed parameter
space, since we display this only in terms of the energy scale $\eta$ (i.e.,
independently of $\mathcal{F}$). Additionally, the green dashed line marks the
minimal characteristic energy scale $\eta$ required for domain wall networks
to generate the PTA signals for $\mathcal{F}=1$. Since a larger contribution
of the fractional energy density in gravitational waves as well as a smaller
value of $\mathcal{F}$ would both imply a stronger bound on $\eta$, only the
region above this line is compatible with the observed PTA signal. There is
then only a tight range of about one order of magnitude in $\eta$ that is
compatible with PTA experiments and avoids domain wall domination. Notice
that, if one considers networks with $\mathcal{F}<1$ — as one would
realistically expect, since domain walls may also emit scalar radiation in
their decay — this allowed range of $\eta$ becomes even narrower. This means,
in particular, that values of $\mathcal{F}\mathrel{\hbox to0.0pt{\lower
4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}0.004$ can be ruled
out if domain walls are to be the main contributors to the PTA signal.

Additionally, under the assumption that the peak of the domain wall
gravitational wave spectrum lies within the frequency interval [$f_{\rm min}$,
$f_{\rm max}$], one may see that the allowed range of $z_{\star}$ is also very
narrow for each value of $\alpha$. Fig. 2 shows the region of the
($z_{\star}$, $\eta$) parameter space compatible with PTA data for $\alpha=1$
and $\alpha=10^{-6}$ (as the opaque dark and light green areas respectively).
For an $\alpha$ of order unity, this would imply a decay time of the domain
wall network shortly after the QCD phase transitions, while the redshift of
decay (roughly) linearly shifts with smaller $\alpha$ to later times as
indicated by the color bar on the right side of the plot. Considering the
constraints on the time of decay and on the energy scale $\eta$ at domain wall
formation, for a given value of $\alpha$, the region of parameter space
consistent with the signal detected by PTAs becomes extremely narrow,
requiring a substantial fine tuning of the properties of the domain wall
networks in order to explain the SGWB detected by the PTAs.

##  V Discussion and Conclusions

In this article we have investigated what are the allowed domain wall
scenarios in light of current CMB and PTA data. Since there are some
uncertainties regarding the spectrum (especially about the shape of its peak),
here we opted to perform our analysis using the total SGWB, instead of its
power spectrum, as this may be characterized with minimal assumptions. The
main uncertainty in this case is related to the efficiency of GW emission in
domain wall collapse, quantified here using the parameter $\mathcal{F}$, but
our main conclusions are mostly independent on this parameter.

Our results show that indirect CMB bounds on $\Omega_{\rm gw}$ at the time of
decoupling have a limited constraining power on domain wall scenarios and
generally provide constraints that are less stringent than those that result
from imposing that domain walls provide only a subdominant contribution to the
cosmic energy budget. Only if the emission of gravitational radiation, and not
scalar radiation, is the dominant decay mechanism of collapsing domain walls —
more precisely if $\mathcal{F}\sim\mathcal{O}(1)$ — does this CMB bound lead
to a decrease of the allowed region of the $(z_{\star},\eta$) parameter space
(and, even in this case, this decrease is not significant). Notice however
that the CMB constraints derived here are conservative because the computation
of $\Omega_{\rm gw}$ in Eqs. (1) and (2) only includes the GWs emitted before
the whole domain wall networks decays at $z_{\star}$. There may also be a
significant emission of GWs during the final collapse of the network, which
should enhance $\Omega_{\rm gw}$ by up to one order of magnitude 222In [9],
the authors find that this decay may be slightly delayed, which may further
enhance the amplitude of this background.. The inclusion of this contribution
may then result in a further reduction of the allowed range of parameter
space, but this — unlike the constraints presented here – is model dependent.

Despite this, we found that even without considering CMB bounds, the range of
allowed energy scales for the domain-wall forming phase transition is already
quite limited, especially for networks that decay early in cosmic history. For
instance, for walls decaying at $z_{\star}\sim 10^{18}$ (at a temperature of
roughly $T_{\star}\sim 100\,\rm{TeV}$), we should have that
$10^{2}\lesssim\eta\lesssim 10^{6}\,\rm{TeV}$ (roughly 4 order of magnitude in
range), while for $z_{\star}\sim 10^{8}$ ($T_{\star}\sim 10^{-1}\,\rm{MeV}$),
the allowed parameter range is much wider: $10^{-1}\lesssim\eta\lesssim
10^{6}\,\rm{MeV}$.

We have also found that, if domain walls are to provide an explanation to the
observed PTA signal, this would require a significant fine tuning of the
parameters of the model. Domain walls would not only need to have a
significant GW emission efficiency — we found that $\mathcal{F}\gtrsim
10^{-3}$, but, as previously explained, we were quite conservative in deriving
this constraint — but their energy density would have to be very close to
dominating the energy density of the universe. As a matter of fact, the
allowed range of parameters is very narrow: one should have, roughly, that
$\eta\sim\mathcal{O}(10^{2})\,\rm{TeV}$ and, for $\alpha=1$,
$T_{\star}\sim\mathcal{O}(10^{2})\,\rm{MeV}$ 333$T_{\star}$ would roughly
decrease proportionally to $\alpha$.. Like in the previous case, this result
was derived based on the contribution to the SGWB before the actual decay of
the network. Here, in fact taking an enhancement of the SGWB due to the
network’s decay into account could potentially increase the range of allowed
parameter space by lowering the minimal domain wall tension required to
reproduce the signal. As discussed before, this would be model dependent and
at most lower the minimal $\sigma$ by one order of magnitude, which would
translate into less than a factor of $1/2$ in the characteristic energy scale.
These results are in very good agreement with those presented in the
literature [9, 5, 11, 23, 24], but unlike those, ours are based on minimal
assumptions and as model-independent as possible. The results in the
literature all stem from numerical simulations of particular domain wall
models [21], in which the energy density of GWs, $\rho_{\rm gw}$, was
measured. These simulations seem to indicate that $\rho_{\rm gw}$ remains
constant and is roughly given by $\sim G\sigma^{2}$ (aside from an efficiency
parameter of order unity), which would correspond to having $\mathcal{F}$
growing linearly with time in this framework. Although such a regime may
indeed exist (as discussed in detail in [14]), it would necessarily have to be
transient as $\mathcal{F}$ should necessarily be smaller than unity at all
times — after all, the energy emitted in the form of GWs cannot ever exceed
the energy lost by the network. As a matter of fact, in such a regime, as
$\mathcal{F}$ increases and $\rho_{\rm gw}$ becomes comparable to the domain
wall energy density, gravitational backreaction should start having a
significant impact on GW emission. Such a regime cannot be captured by
numerical simulations and the last stages of collapse — in which walls becomes
ultra-relativistic and the emission GWs may be quite significant [14] — cannot
as well. This casts some doubt as to the ability of numerical simulations to
describe realistic domain wall scenarios. Assuming, as in simulation-inferred
models, that $\mathcal{F}\propto t$ results in an additional factor of
$G\sigma/H_{\star}$ in $\Omega_{\rm gw}$ (note that in our model $\rho_{\rm
gw}$ scales as $\sigma H_{\star}$), but despite this our interpretations of
PTA data roughly coincide. This coincidence may be easily explained: since, in
both cases, domain wall density would have to be comparable to that of the
background density at the time of decay, we should have, in both cases, that
$\sigma H_{\star}\sim H_{\star}^{2}/G$, which implies that
$G\sigma/H_{\star}\sim 1$. This is precisely the limit in which both
approaches agree — they would not if domain walls were to be subdominant — but
it is also the limit in which the potential impact of gravitational
backreaction is most severe.

###### Acknowledgements.

The authors would like to thank Dr. Paul Tol for developing the colorblind-
friendly color schemes [26] used in this work. D. G. is supported by FCT –
Fundação para a Ciência e a Tecnologia through the PhD fellowship
2020.07632.BD. L. S. is supported by FCT through contract No. DL
57/2016/CP1364/CT0001. Funding for this work has also been provided by FCT
through the research grants UIDB/04434/2020 and UIDP/04434/2020 and through
the R & D project 2022.03495.PTDC – Uncovering the nature of cosmic strings.

## References

  * Agazie  _et al._ [2023] G. Agazie _et al._ (NANOGrav), The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background, Astrophys. J. Lett. 951, L8 (2023), arXiv:2306.16213 [astro-ph.HE] .
  * Reardon  _et al._ [2023] D. J. Reardon _et al._ , Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array, Astrophys. J. Lett. 951, L6 (2023), arXiv:2306.16215 [astro-ph.HE] .
  * Antoniadis  _et al._ [2023] J. Antoniadis _et al._ (EPTA, InPTA:), The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals, Astron. Astrophys. 678, A50 (2023), arXiv:2306.16214 [astro-ph.HE] .
  * Xu  _et al._ [2023] H. Xu _et al._ , Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I, Res. Astron. Astrophys. 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Afzal  _et al._ [2023] A. Afzal _et al._ (NANOGrav), The NANOGrav 15 yr Data Set: Search for Signals from New Physics, Astrophys. J. Lett. 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Vilenkin and Shellard [2000] A. Vilenkin and E. P. S. Shellard,  _Cosmic Strings and Other Topological Defects_ (Cambridge University Press, 2000).
  * Zeldovich  _et al._ [1974] Ya. B. Zeldovich, I. Yu. Kobzarev, and L. B. Okun, Cosmological consequences of the spontaneous breakdown of discrete symmetry, Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki 67, 3 (1974).
  * Huang  _et al._ [2024] L. Huang, Y. Wang, H.-X. Zhang, S. Matsuzaki, H. Ishida, M. Kawaguchi, and A. Tomiya, Impact of local CP-odd domain in hot QCD on axionic domain-wall interpretation for NANOGrav 15-year Data, (2024), arXiv:2403.11444 [hep-ph] .
  * Ferreira  _et al._ [2024] R. Z. Ferreira, A. Notari, O. Pujolàs, and F. Rompineve, Collapsing Domain Wall Networks: Impact on Pulsar Timing Arrays and Primordial Black Holes, (2024), arXiv:2401.14331 [astro-ph.CO] .
  * Dunsky and Kongsore [2024] D. I. Dunsky and M. Kongsore, Primordial Black Holes from Axion Domain Wall Collapse, (2024), arXiv:2402.03426 [hep-ph] .
  * Kitajima  _et al._ [2024] N. Kitajima, J. Lee, K. Murai, F. Takahashi, and W. Yin, Gravitational waves from domain wall collapse, and application to nanohertz signals with QCD-coupled axions, Phys. Lett. B 851, 138586 (2024), arXiv:2306.17146 [hep-ph] .
  * Blasi  _et al._ [2023] S. Blasi, A. Mariotti, A. Rase, A. Sevrin, and K. Turbang, Friction on ALP domain walls and gravitational waves, JCAP 04, 008, arXiv:2210.14246 [hep-ph] .
  * Sousa and Avelino [2015] L. Sousa and P. P. Avelino, Cosmic Microwave Background anisotropies generated by domain wall networks, Phys. Rev. D 92, 083520 (2015), arXiv:1507.01064 [astro-ph.CO] .
  * Grüber  _et al._ [2024] D. Grüber, L. Sousa, and P. P. Avelino, Stochastic gravitational wave background generated by domain wall networks, (2024), arXiv:2403.09816 [gr-qc] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Avelino and Sousa [2011] P. P. Avelino and L. Sousa, Domain wall network evolution in (N+1)-dimensional FRW universes, Physical Review D 83, 043530 (2011), arxiv:1101.3360 .
  * Sousa and Avelino [2011a] L. Sousa and P. P. Avelino, p-brane dynamics in (N+1)-dimensional FRW universes: a unified framework, Phys. Rev. D 83, 103507 (2011a), arXiv:1103.1381 [hep-th] .
  * Sousa and Avelino [2011b] L. Sousa and P. P. Avelino, The cosmological evolution of p-brane networks, Phys. Rev. D 84, 063502 (2011b), arXiv:1107.4582 [hep-th] .
  * Leite and Martins [2011] A. M. M. Leite and C. J. A. P. Martins, Scaling Properties of Domain Wall Networks, Physical Review D 84, 103523 (2011), arxiv:1110.3486 [astro-ph, physics:hep-ph, physics:hep-th] .
  * Binetruy  _et al._ [2012] P. Binetruy, A. Bohe, C. Caprini, and J.-F. Dufaux, Cosmological Backgrounds of Gravitational Waves and eLISA/NGO: Phase Transitions, Cosmic Strings and Other Sources, JCAP 06, 027, arXiv:1201.0983 [gr-qc] .
  * Hiramatsu  _et al._ [2014] T. Hiramatsu, M. Kawasaki, and K. Saikawa, On the estimation of gravitational wave spectrum from cosmic domain walls, Journal of Cosmology and Astroparticle Physics 2014 (02), 031, arxiv:1309.5001 .
  * Clarke  _et al._ [2020] T. J. Clarke, E. J. Copeland, and A. Moss, Constraints on primordial gravitational waves from the Cosmic Microwave Background, JCAP 10, 002, arXiv:2004.11396 [astro-ph.CO] .
  * Zhang  _et al._ [2023] Z. Zhang, C. Cai, Y.-H. Su, S. Wang, Z.-H. Yu, and H.-H. Zhang, Nano-Hertz gravitational waves from collapsing domain walls associated with freeze-in dark matter in light of pulsar timing array observations, Phys. Rev. D 108, 095037 (2023), arXiv:2307.11495 [hep-ph] .
  * Gouttenoire and Vitagliano [2023] Y. Gouttenoire and E. Vitagliano, Domain wall interpretation of the PTA signal confronting black hole overproduction, (2023), arXiv:2306.17841 [gr-qc] .
  * Caprini  _et al._ [2009] C. Caprini, R. Durrer, T. Konstandin, and G. Servant, General Properties of the Gravitational Wave Spectrum from Phase Transitions, Phys. Rev. D 79, 083519 (2009), arXiv:0901.1661 [astro-ph.CO] .
  * [26] Paul tol’s colour schemes and templates, https://personal.sron.nl/~pault/, accessed: 22.05.2024.


# Discrete calculus of variations and Boltzmann distribution without
Stirling’s approximation

Q. H. Liu<EMAIL_ADDRESS>School for Theoretical Physics, School of
Physics and Electronics, Hunan University, Changsha 410082, China

###### Abstract

A _double extrema form_ of the calculus of variations is put forward in which
only the smallest one of the finite differences is physically meaningful to
represent the variational derivatives defined on the discrete points. The most
probable distribution for the Boltzmann system is then reproduced without the
Stirling’s approximation, and free from other theoretical problems.

Boltzmann distribution, most probable distribution, finite difference,
Discrete calculus of variations

##  I Introduction

Boltzmann distribution lies in the heart of statistics physics with
applications in mathematics, chemistry and engineering, and it was first given
more than 100 years ago Boltzmann . However, seldom feels comfortable with its
derivation with use of the Stirling’s approximation 61 ; 62 ; 02 ; 1989 ; 2009
. A more puzzling point is that, if one uses more accurate expression of the
Stirling’s approximation, or attempts to exact the derivation of the Boltzmann
distribution without invoking the Stirling’s approximation, he always runs
into the a surplus $\sim-1/2$ term 1989 ; 2009  in the distribution. One can
find the relevant derivation in most standard statistical physics or physical
chemistry textbooks, for instance, 1997 ; greiner . The key steps are outlined
in the following.

Considering a gaseous system of $N$ noninteracting, indistinguishable
particles confined to a space of volume $V$ and sharing a given energy $E$.
Let $\varepsilon_{i}$ denote the energy of $i$-th level and
$\varepsilon_{1}\prec\varepsilon_{2}\prec\varepsilon_{3}\prec...$, and $g_{i}$
denote the degeneracy of the level. In a particular situation, we may have
$n_{1}$ particles in the first level $\varepsilon_{1}$, $n_{2}$ particles in
the second level $\varepsilon_{2}$, and so on, defining a distribution set
$\left\\{n_{i}\right\\}$ 1997 . The number of the distinct microstates in set
$\left\\{n_{i}\right\\}$ is then given by,

$\Omega\left\\{n_{i}\right\\}=\prod_{i}\frac{\left(g_{i}\right)^{n_{i}}}{n_{i}!}.$
(1)

The distribution set $\left\\{n_{i}\right\\}$ must conform to two macroscopic
conditions,

$\sum_{i}n_{i}=N,\sum_{i}n_{i}\varepsilon_{i}=E.$ (2)

The entropy of the system is
$S=S\left(N,V,E\right)=k_{B}\ln\Omega\left\\{n_{i}^{\ast}\right\\}$ where
$k_{B}$ is the Boltzmann constant, and $\left\\{n_{i}^{\ast}\right\\}$ is the
most probable distribution, determined by maximizing the following
variational,

$f=k_{B}\ln\Omega\left\\{n_{i}\right\\}-\alpha\left(\sum_{i}n_{i}-N\right)-\beta\left(\sum_{i}n_{i}\varepsilon_{i}-E\right),$
(3)

where $\alpha$ and $\beta$ are two Lagrangian multipliers. The variational
calculation gives with the discrete indices $i=1,2,3,...$1989 ; 2009 ; 1997 ;
greiner ,

$\delta f=\sum_{i}\left(\delta n_{i}\left(\ln
g_{i}-\alpha-\beta\varepsilon_{i}\right)-\delta\ln n_{i}!\right).$ (4)

The most probable distribution $\left\\{n_{i}^{\ast}\right\\}$ satisfies 1989
; 2009 ; 1997 ; greiner ,

$\frac{\delta f}{\delta n_{i}^{\ast}}=\left(\ln
g_{i}-\alpha-\beta\varepsilon_{i}\right)-\frac{\delta\ln n_{i}^{\ast}!}{\delta
n_{i}^{\ast}}=0,\text{ i.e., }\psi(1+n_{i}^{\ast})=\ln
g_{i}-\alpha-\beta\varepsilon_{i},$ (5)

where

$\psi(1+x)\equiv\frac{d\ln x!}{dx}$ (6)

is the so-called Digamma function of variable $x$, and $\Gamma(1+x)\equiv x!$
is the Gamma function.

When $x\gg 1$, utilizing the usual form of the Stirling’s approximation $\ln
x!\approx x\ln x-x$, we have,

$\psi(1+x)\approx\ln x.$ (7)

In the limit $x\gg 1$, the standard expression for the Boltzmann statistics is
recovered 1989 ; 2009 ; 1997 ; greiner ,

$n_{i}^{\ast}\approx g_{i}e^{-\alpha-\beta\varepsilon_{i}}.$ (8)

However, it is not satisfactory. We can use grand ensemble statistics to prove
that the approximately-equal-to symbol ”$\approx$” is in fact the identical-to
one ”$=$” greiner , and a more convenient and detailed derivation is available
from a website lecture note lectnote .

The Stirling’s approximation $\ln x!\approx x\ln x-x$ holds only when $x$ is
large. When $x$ is small, a slight deviation from the distribution
$n_{i}^{\ast}\approx g_{i}e^{-\alpha-\beta\varepsilon_{i}}$ as
$i\rightarrow\infty$ may cause a ”long tail” contribution. How about a more
accurate approximation of $\ln x!$ is used? In the following section II, we
will show the continuous calculus of variations in step (4) is questionable.
In section III, with the proper use of the discrete calculus of variations,
the relevant difficulty is removed, and the approximate relation (7) turns out
to be an exact result. In section IV, some further comments on the derivatives
in the discrete calculus of variations are added. The final section V is a
brief conclusion.

##  II Appearance of additional $-1/2$ term with the continuous calculus of
variations

With more accurate expression of $\ln x!$ as $\ln
x!\approx\left(x+1/2\right)\ln x-x+1/2\ln 2\pi$, we have $d\ln x!/dx\approx\ln
x+1/\left(2x\right)$, and $\exp(d\ln x!/dx)\approx x+1/2$. Then we come up
against a worse result,

$n_{i}^{\ast}\approx-\frac{1}{2}+g_{i}e^{-\alpha-\beta\varepsilon_{i}}.$ (9)

When the energy levels are lower to the lowest, i.e., the number
$n_{i}^{\ast}$ of particles is large, we have approximately the Boltzmann
distribution $n_{i}^{\ast}\simeq g_{i}e^{-\alpha-\beta\varepsilon_{i}}$.
Meanwhile, the high energy levels are less likely to be populated as
$g_{i}e^{-\alpha-\beta\varepsilon_{i}}\rightarrow 0$ when
$\varepsilon_{i}\rightarrow\infty$, leading to a ”long tail” distribution from
(9),

$n_{i}^{\ast}\approx-\frac{1}{2},i\rightarrow\infty.$ (10)

This additional $-1/2$ term presents though the ”rigorous treatment” of Eq.
(4) is utilized. To see it, let us consider an asymptotic expression of the
Digamma function $\psi(1+x)$ in limit $x\rightarrow\infty$ which is, with
$O\left(x^{-2}\right)$ denoting a quantity of order $x^{-2}$ psi ,

$\psi(1+x)\rightarrow\ln x+\frac{1}{2x}+O\left(x^{-2}\right),\text{
}x\rightarrow\infty.$ (11)

A surprising property of the Digamma function $\psi(1+x)$ is that we can
numerically verify following equation which holds true for all
$x\in\left[0,\infty\right)$,

$e^{\psi(1+x)}\approx x+1/2.$ (12)

The largest deviation from $1/2$ occurs at the ending point $x=0$ and
$e^{\psi(1)}=e^{-\gamma}\approx 0.561$, with $\gamma\simeq 0.577$ is the
Euler–Mascheroni constant. As $x$ increases from zero, the difference
$e^{\psi(1+x)}-x$ monotonically converges to $1/2$; and at $x=10$ and
$10^{2}$, it gives $0.503$ and $0.500$, respectively. Clearly, Eq. (9) turns
out to be approximately valid _in whole range_ of $n_{i}^{\ast}\in[0,\infty)$.
When $x$ is large, $x+1/2\approx x$; however, when $x=0$,
$e^{\psi(1+x)}\approx 0.561\approx 1/2$ leads to the spurious ”long tail”
distribution again.

The additional $-1/2$ term has to be removed, because the energy levels can
never be _negatively_ occupied. However, one can not simply discard the $-1/2$
term in the limit $i\rightarrow\infty$, because this term is associated with a
divergence as the number of particles is summed up over all energy levels,

$\sum_{i}^{\infty}n_{i}^{\ast}\approx-\frac{1}{2}\sum_{i=1}^{\infty}1+\sum_{i=1}^{\infty}g_{i}e^{-\alpha-\beta\varepsilon_{i}},\text{
and }\sum_{i=1}^{\infty}1\rightarrow\infty.$ (13)

This situation bears a resemblance to the famous divergence in quantum
electrodynamics. Some feels comfortable with it but other are strongly against
it. Remembering that the lattice quantum fields are helpful to eliminate the
undesired divergences, we are confident that the similar and proper treatment
of the discrete calculus can be used to get rid of the difficulty. To note
that the discrete calculus of variations is a well-established discipline in
mathematics with much wider applications, and one of the fundamental
principles is clear for the derivatives should be replaced by finite
differences 1968 ; guo ; 2006 ; 2014 .

##  III Removal of additional $-1/2$ term with discrete calculus of variations

To note that the change of the number of the particle should be integers. The
application of the variational calculation to the Boltzmann system in Eq. (4)
must be with care. Since the number of particle is discrete, we must thus
define the Digamma function on discrete lattices $x=0,1,2,3,...$. Furthermore,
the smallest possible change of the particle can only be $+1$ or $-1$ and can
never be an infinitesimal.

The minimum finite difference for $\Delta\ln n!/\Delta n$ may take two forms.
One is the so-called the forward one $\psi^{+}(1+n)$, and another the backward
one $\psi^{+}(1+n)$,

$\displaystyle\frac{\Delta\ln n!}{\Delta n}$ $\displaystyle=$
$\displaystyle\frac{\ln\left(n+1\right)!-\ln\left(n\right)!}{1}=\psi^{+}(1+n),\text{
and, }$ (14a) $\displaystyle\frac{\Delta\ln x!}{\Delta n}$ $\displaystyle=$
$\displaystyle\frac{\ln\left(n\right)!-\ln\left(n-1\right)!}{1}=\psi^{-}(1+n).$
(14b) In simple forms, we have, respectively,  $\displaystyle\psi^{+}(1+n)$
$\displaystyle\equiv$ $\displaystyle\ln(n+1)!-\ln n!=\ln(n+1)\text{,
}\left(n\succeq 0\right),\text{ and}$ (15a) $\displaystyle\psi^{-}(1+n)$
$\displaystyle\equiv$ $\displaystyle\ln n!-\ln(n-1)!=\ln n.$ (15b) Two
important and somewhat trivial relations are $\ln(n+1)\succ\ln n$,
$\left(n\succeq 1\right)$, and when $n\rightarrow 0$, $\ln(n+1)\succ\ln n$ can
also be understood in the sense of $1\approx e^{\ln(n+1)}\succ e^{\ln
n}\rightarrow 0$.

The principles of the calculus of variations indicate that a variation of $\ln
n!$ is to find _a minimum possible amount_ of the changes of $\ln n!$ when $n$
is changed by a minimum possible value. Thus, $\psi^{-}(1+n)=$ $\ln n$ is
singled out. I.e., the approximate result (7) happens to be an exact one. Now
we have an exact relation,

$\text{ }n_{i}^{\ast}=g_{i}e^{-\alpha-\beta\varepsilon_{i}}.$ (16)

It is right the standard form of the Boltzmann distribution applicable in the
whole meaningful interval $n_{i}^{\ast}\in(0,n_{0}^{\ast})$.

##  IV Why the additional $-1/2$ term occcurs from the point of the discrete
calculus of variations

In above section, we stress that only a special choose of discrete derivatives
among $\Delta\ln n!/\Delta n$ suffices to represent the discrete variational
derivative $\delta\ln n!/\delta n$, and we single out the minimum one
$\psi^{-}(1+n)=\ln n$. In contrast, if including both $\psi^{+}(1+n)=\ln(n+1)$
and $\psi^{-}(1+n)=\ln n$, and using the arithmetic mean of $\psi^{+}(1+n)$
and $\psi^{-}(1+n)$ to present the discrete variational derivative $\delta\ln
n!/\delta n$ instead, we have,

$\psi_{\min}(1+n)\equiv\frac{1}{2}\left(\psi^{+}(1+n)+\psi^{-}(1+n)\right)=\left(\ln
n+\ln(n+1)\right)/2.$ (17)

we run across the unnecessary $-1/2$ term again, for we have,

$\exp\left(\psi_{\min}(1+n)\right)=\sqrt{n(n+1)}\approx n+1/2,(n\gg 1).$ (18)

Such the discrete variational derivative is applicable for the definition of
thermodynamic quantities for the _few-particle systems in the microcanonical
ensemble_. Since in the ensemble both the number of particle and the energy
are fixed, the temperature, for instance, can not be defined by usual
derivative, but the variational derivative instead. Following the same idea
presented in above section, we must choose a minimum possible value among all
possible differences. Recent treatments few  of the definition of the
temperature and other thermodynamic quantities in the microcanonical ensemble
for small systems demonstrate that we must pick up the backward difference of
$\Delta\ln n!/\Delta n$, based on the completely different arguments.

##  V Conclusions and discussions

The calculus of variations is to use a small changes in functions and
functionals, to find maxima and minima of them. Once the functions and
functionals are defined on the continuous intervals, the usual derivatives are
sufficient. However, once they are defined on the discrete points, the
variational derivatives correspond to the many finite differences. We in fact
put forward a _double extrema form_ of the calculus of variations in which
only the smallest one of the finite differences is physically meaningful to
represent the variational derivatives. By the _double extrema form_ , we mean
that the usual form of the calculus of variations deals with one extrema only.

The _double extrema form_ of the variational (4) means not only $\delta f=0$,
but also $\delta\ln n!/\delta n=\min\left\\{\Delta\ln n!/\Delta n\right\\}$.
The exact form the most probable distribution for the Boltzmann system is then
reproduced without the Stirling’s approximation, and free from other
theoretical problems. The _double extrema form_ of the calculus of variations
has in fact been used in literature on statistical mechanics for finite size
systems, but based on the physical arguments.

###### Acknowledgements.

This work is financially supported by National Natural Science Foundation of
China under Grant No. 11675051.

## References

  * (1) L. Boltzmann, Wissenschaftliche Abhandlungen, Vol. I, F. Hasenöhrl (ed.), (Leipzig: Barth, 1909; and New York: Chelsea, 1969 (reissued)).
  * (2) R. W. Hakala, A new derivation of the Boltzmann distribution law, J. Chem. Educ. 38(1), 33 (1961).
  * (3) P. A. H. Wyatt, Elementary statistical mechanics without Stirling’s approximation, J. Chem. Educ. 39(1), 27(1962) .
  * (4) R. B. Shirts, and M. R. Shirts, EDeviations from the Boltzmann distribution in small microcanonical quantum systems: Two approximate one-particle energy distributions, J. Chem. Phys. 117, 5564(2002).
  * (5) V. J. Menon and D. C. Agrawal, Method of most probable distribution: new solutions and results, Pramana Journal of Physics 33, 455(1989).
  * (6) S. Kakorin, Revision of Boltzmann statistics for a finite number of particles, Am. J. Phys. 77, 48(2009).
  * (7) R. K. Pathria, Statistical Mechanics, 2nd Ed., (Oxford: Butterworth-Heinemann, 1997).
  * (8) W. Greiner, L. Neise, H. Stöker, Thermodynamics and statistical mechanics, (New York: Springer-Verlag, 1995) pp.297.
  * (9) Clare C. Yu, https://ps.uci.edu/~cyu/p115A/LectureNotes/Lecture13/lecture13.pdf.
  * (10) J. A. Cadzow, Discrete calculus of variations, International Journal of Control 11, 393(1970).
  * (11) Han-Ying Guo, On variations in discrete mechanics and field theory, Journal of Mathematical Physics 44, 5978(2003).
  * (12) J. N. Reddy, An Introduction to the Finite Element Method (3rd ed.). (New York: McGraw-Hill, 2006).
  * (13) K. J. Bathe, Finite Element Procedures, (2nd ed.), (Watertown: MA. 2014).
  * (14) M. Abramowitz, and C. A. Stegun, (eds.). “Psi (Digamma) Function.” in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. (New York: Dover, 1972) pp. 258-259.
  * (15) E. N. Miranda, Statistical mechanics of few-particle systems: exact results for two useful models, Eur. J. Phys. 38, 065101(2017).