Validity of path thermodynamic description of

reactive systems: Microscopic simulations

F. Baras${}^{\,a}$, Alejandro L. Garcia${}^{\,b}$, and M. Malek
Mansour${}^{\,c}$

(a) Laboratoire Interdisciplinaire Carnot de Bourgogne,
UMR 6303 CNRS-Université Bourgogne Franche-Comté,
9 Avenue A. Savary, BP 47 870,
F-21078 Dijon Cedex, France

(b) Dept. Physics and Astronomy,
San Jose State University,
San Jose, California, 95192 USA

(c) Université Libre de Bruxelles CP 231, Campus Plaine,
B-1050 Brussels, Belgium

###### Abstract

Traditional stochastic modeling of reactive systems limits the domain of
applicability of the associated path thermodynamics to systems involving a
single elementary reaction at the origin of each observed change in
composition. An alternative stochastic modeling has recently been proposed to
overcome this limitation. These two ways of modeling reactive systems are in
principle incompatible. The question thus arises about choosing the
appropriate type of modeling to be used in practical situations. In the
absence of sufficiently accurate experimental results, one way to address this
issue is through the microscopic simulation of reactive fluids, usually based
on hard-sphere dynamics in the Boltzmann limit. In this paper, we show that
results obtained through such simulations unambiguously confirm the
predictions of traditional stochastic modeling, invalidating a recently
proposed alternative.

##  1 Introduction

Thermodynamic description of non-equilibrium systems is based on the concept
of entropy production, the correct evaluation of which necessarily requires
full knowledge of each elementary process therein [1, 2]. This prerequisite
plays an important role in the corresponding statistical formulation, commonly
modeled by means of an appropriate Markovian stochastic process. In previous
papers, we proved that the thermodynamic formulation of homogeneous reactive
systems based on their sample path, widely known as “path thermodynamics”,
will lead to erroneous results whenever the system involves more than one
elementary reaction leading to the same change in composition [3, 4]. Similar
observations were reported in [5] and [6] (see also [7]).

This result can also be appreciated from a physical perspective. Let us
consider a perfectly homogeneous (no local fluctuations) isothermal reactive
system, such as can be produced experimentally in a "continuously stirred tank
reactor" (CSTR). The state of such a system is entirely characterized by its
composition, which is indeed the only pertinent quantity accessible to
experimental investigations. Suppose now that we have at our disposal an ideal
experimental device allowing us to measure the precise number of each chemical
species over an arbitrarily long interval of time. Using this device, we may
gain access to the exact state trajectory of the system, traditionally
referred to as its "sample path" within the framework of stochastic processes.
The question is whether such a state trajectory encompasses sufficient
information to allow us to determine the thermodynamic properties of that
system. The key issue here is obviously the fact that entropy production in
reactive systems is the sum of the entropy production that is associated with
each individual elementary reaction (cf. Section 9.5 in [1]). The answer to
this question will thus be positive if and only if each of the measured
changes in composition can be attributed to a specific elementary reaction.
Such an attribution becomes impossible if the reactive system includes more
than one elementary reaction leading to the same change in composition. No
matter what method we use, the path thermodynamic properties of such a
reactive system can never be determined from its state trajectory.

Despite the undeniability of this result, from both the mathematical and the
physical point of view, its validity was recently contested by Gaspard in a
Comment article [8]. Yet this very issue had already been raised by that same
author over a decade ago [9]. In short, extending the work of Lebowitz and
Spohn [10], Gaspard initiated the path thermodynamic formulation of reactive
systems in 2004 [11, 12]. The theory is sound but, like all scientific
theories, its domain of applicability has limits; specific situations will
exist outside the scope of the theory. Indeed, three years later that theory
was found to encounter certain inconsistencies when applied to the
Schnakenberg graph formulation of "current" fluctuations in reactive systems
involving more than one elementary reaction leading to the same change in
composition [9]. A remedy was proposed in the specific case of "current"
fluctuations but, as we have shown in [3], this remedy is inapplicable to
traditional modeling of reactive systems based on jump Markov processes.

In short, a pure jump process, $\chi(t)$, is entirely determined by the
concept of "transition rates", $W(X\,|\,X^{\prime})$, defined by Kolmogorov as
[13]

$P(X,t+\Delta t\,|\,X^{\prime},\,t)=W(X\,|\,X^{\prime})\,\Delta t\,+\,{\rm
o}(\Delta t)\,,\,\,\forall\,X\,\neq\,X^{\prime}$ (1)

The function $P(X,t+\Delta t\,|\,X^{\prime},t)$ represents the conditional
probability to have $\chi(t+\Delta t)=X$, given that $\chi(t)=X^{\prime}$. The
description of jump Markov processes in terms of their sample path is based on
this fundamental Kolmogorov equality [14]. Consider now a reactive system that
involves $n$ elementary reactions ${\cal R}_{1}\cdots{\cal R}_{n}$ leading the
same change in composition $X\rightarrow X^{\prime}$, and denote by
$W_{\rho}(X\,|\,X^{\prime})$ the transition rate associated with the reaction
$\mathcal{R}_{\rho}$. It was claimed in [9, 8] that we can write the sample
path of $\chi(t)$ in terms of any individual transition rate
$W_{\rho}(X\,|\,X^{\prime})$. This claim, however, contradicts the fundamental
property of jump processes.

Recall that the probability associated with a random event is unique. This
property results from the basic definition of the concept of probability,
first stated by Pascal, and refined over the years by various mathematicians
to Kolmogorov’s axiomatic formulation [15, 16]. A transition
$X^{\prime}\rightarrow X$ is clearly a random event. The probability that this
event occurs in the time interval $[t,\,t+\Delta t]$ is precisely
$P(X,t+\Delta t\,|\,X^{\prime},t)$. Accordingly, this probability cannot
assume several values simultaneously, that is, for any given $X,X^{\prime},t$,
and $\Delta t$, the conditional probability distribution $P(X,t+\Delta
t\,|\,X^{\prime},t)$ is unique. The Kolmogorov equality (1) then implies that
this is also the case for the transition rate $W(X\,|\,X^{\prime})$. In
particular, if a reactive system involves several elementary reactions leading
to the same change in composition, then the resulting transition rate is
necessarily the sum of the transition rates associated with each of them,
i.e., $W(X\,|\,X^{\prime})=\sum_{\rho}W_{\rho}(X\,|\,X^{\prime})$.
Consequently, expressing the sample path of $\chi(t)$ in terms of individual
transition rates $W_{\rho}(X\,|\,X^{\prime})$ does not apply to jump
stochastic processes.

Given this limitation, Gaspard proposed in his Comment a brand new type of
stochastic modeling of reactive systems that extends the domain of validity of
the associated path thermodynamics to the controversial situation of reactive
systems involving more than one elementary reaction leading to the same change
in composition [8]. The resulting stochastic process proves to be quite
different from that associated with the traditional stochastic modeling of
reactive systems. They don’t even share the same state space. The state space
of the stochastic process associated with the traditional modeling of an $n$
component isothermal homogeneous (CSTR) reactive system is simply
$\mathbb{Z}^{n}$, in one-to-one correspondence with quantities that can
actually be measured in laboratory experiments ($\mathbb{Z}$ represents the
set of non-negative integers). First proposed by McQuarrie in 1967 [17],
traditional stochastic modeling was then refined over the years by Van Kampen
[15], Haken [18], Nicolis and Prigogine [19], Kurtz [20], Gardiner [16], and
many others. Everything we currently know about the statistical properties of
non-equilibrium reactive systems was established in this way, including the
stochastic formulation of path thermodynamics by Seifert [21], Lebowitz and
Spohn [10], and Gaspard [11, 12].

The state space of the stochastic process associated with the new type of
modeling proposed in [8] is different. In addition to the number of particles
of chemically active molecules, it also includes an extra variable designed to
select the precise elementary reaction at the origin of an observed change in
composition. As clearly stated by the author, the main consequence of this
modeling is that it allows us to define several Markov processes associated
with a given reactive system (Section III, in [8]). Among them, there will
undoubtedly be found a process guaranteeing the validity of the associated
path thermodynamics in controversial situations.

One may argue that in this new formulation the model has been changed simply
to shoehorn it to fit the theory. As they have different state spaces, it is
impossible to use the framework of one of these stochastic processes to
confirm or deny the validity of the other. Conversely, it is impossible to
discuss the validity of the proposed new modeling within the framework of the
traditional stochastic modeling.

Nevertheless, the main purpose of the two approaches is to provide a
theoretical description of a "real-world" system: a system with properties
that can be observed in laboratory experiments. In this respect, we know that
reactive processes result from local interactions between chemically active
molecules (reactive collisions). If the reactive system is well stirred
(CSTR), it will be impossible to determine experimentally which reaction led
to an observed change in composition, unless that reaction is unique. In other
words, determining the precise state trajectory of a well-stirred reactive
system, as required by the new modeling [8], is beyond the reach of real-life
laboratory experiments.

On the other hand, the new stochastic modeling of reactive systems, however
strange, may nevertheless represent reality. The only way to investigate this
possibility is through laboratory experiments. Unfortunately, the accuracy of
the available experimental data are insufficient to address this issue. The
alternative option is to perform microscopic simulations of reactive fluids,
usually based on Newtonian hard sphere dynamics. Introduced in the mid-
seventies [22, 23], this technique provides quite useful information on the
relevance and accuracy of theoretical developments in non-equilibrium reactive
systems [24, 25, 26, 27, 28] (see [29] for a review).

In Section 2, we consider the microscopic simulation of Schlögl-like reactive
systems [30], often used to illustrate some peculiar aspects of path
thermodynamic properties of reactive systems [11, 12]. The results obtained
are in perfect agreement with the traditional stochastic modeling of reactive
systems, thus calling into question the alternative modeling proposed by
Gaspard. Conclusions and perspectives are presented in Section 3; algorithmic
details of the simulations are in an appendix.

##  2 Microscopic simulation

Let us consider the reactive system

$A\,\,+\,\,X\,\,\,\mathop{\rightleftharpoons}^{k_{1}}_{k_{-1}}\,\,\,2\,X\quad\quad
B\,\,+\,\,C\,\,\,\mathop{\rightleftharpoons}^{k_{2}}_{k_{-2}}\,\,\,B\,+\,X$
(2)

where the mole fractions of the reactants $A$, $B$ and $C$ are supposed to
remain constant. Following the traditional practice, we shall use the same
symbols (here $A,B,C,X$) to denote the number of particles of the
corresponding chemical species. Upon setting $k_{-1}=k_{1}$, $k_{-2}=k_{2}$,
$C=A/2$ and $A/B=k_{2}/k_{1}=5/6$, the number of $X$ particles at the
stationary state reads $X_{s}=A$. This choice of parameters guaranties that
the system operates under non-equilibrium conditions. For instance, we can
check that in dilute (ideal) systems, the thermodynamic entropy production,
$\sigma_{s}$, is strictly positive at the stationary state

$\sigma_{s}\,=\,\frac{6}{5}\,N\,k_{B}\,k_{2}\,a^{2}\,\ln(2)\,\,>\,\,0$ (3)

where $k_{B}$ is the Boltzmann constant, $a$ the mole fraction of $A$, and $N$
the total number of particles present in the system, including solvent or
other non-reactive particles (extensivity parameter) [1].

The model (2) belongs to the class of one variable reactive systems where the
reactions lead either to the change in composition $X\rightarrow X+1$
(forward) or $X\rightarrow X-1$ (backward). A well known example is the
Schlögl model [30]. As such, the associated state trajectory (sample path)
contains no information allowing one to distinguish these reactions from each
other. We showed that this property implies necessarily that the state
trajectory is time-reversible at the stationary regime, so that the resulting
path entropy production is zero [3, 4]. It is worthwhile to recall that the
time-reversibility of Schlögl type of reactive systems at the stationary
regime was first established by Graham and Haken in 1971 for diffusion
processes [31]. It was then generalized by several authors in order to include
jump processes as well (see for example [15]). An exhaustive demonstration is
given in Section 6.3 of Gardiner’s textbook [16].

The validity of this result was contested by Gaspard [8], who proposed a new
type of stochastic modeling of reactive systems where the resulting path
entropy production satisfies, in average, the macroscopic thermodynamic
prediction (3). In the absence of sufficiently accurate experimental results,
an alternative way to clarify this issue is through the microscopic simulation
of the reactive system (2). Measuring the number of $X$ particles over a
sufficiently long time interval allows us to estimate the joint probability
distribution $P(X,t\,;X^{\prime},t+\tau)$, for a given $\tau>0$, and compare
it with the reverse probability distribution $P(X^{\prime},t\,;X,t+\tau)$. The
main challenge with this approach is that obtaining statistically reliable
results requires an extremely large number of data points.

The simplest way to deal with this problem is to resort to the simulation of
the reactive Boltzmann equation for which there exists a well-established
algorithm introduced half a century ago by G.A. Bird [32]. Not only is this
technique much faster than traditional hard spheres molecular dynamics, but it
can also be adapted easily to simulate a perfectly homogeneous systems while
conserving the main characteristics of microscopic dynamics. Even though the
latter procedure is a well-established and widely known algorithm [33], for
the sake of completeness we present in Appendix A a short review of its basic
features.

For the microscopic simulation we consider an assembly of $N=5,000$ hard
spheres of diameter $d$ confined in a box of volume $V$ with a number density
of $3\times 10^{-3}$ particles per $d^{3}$. This choice guaranties that the
system is well within the range of validity of Boltzmann equation (see ref.
[32] for more details). The other parameters are set as follows: $A=1000$
particles (thus $B=1200$ and $C=500$ particles) and $k_{2}=0.75\times\nu$,
where $\nu$ is the collision frequency. In other words, $75\%$ of collisions
between $A$ and $B$ or $B$ and $X$ particles (second reaction) are assumed to
be reactive (recall that $k_{1}=6/5\,k_{2}=0.9\,\nu$). Note that, by
definition, the number of solvent particles is changing over time since their
role is to maintain constant the number of $A$, $B$ and $C$ particles (cf.
Appendix A for details). We verified that the solvent mole fraction remains
always above $20\%$ during the simulation.

In addition, we considered a slightly different reactive system:

$A\,\,+\,\,A\,\,\,\mathop{\rightleftharpoons}^{k_{1}}_{k_{-1}}\,\,\,X\,\,+\,\,X\quad;\quad
B\,\,+\,\,C\,\,\,\mathop{\rightleftharpoons}^{k_{2}}_{k_{-2}}\,\,\,B\,+\,X$
(4)

Unlike the reactive system (2), here each elementary reaction leads to a
different change in composition so that the corresponding state trajectory is
expected to be time-irreversible in a non-equilibrium stationary regime. The
comparison of the statistical properties of these two systems proves to be
quite helpful for the correct interpretation of the results. Setting
$k_{1}=k_{2}=k_{-1}=k_{-2}$, $B=5\,A/3$ and $C=A/3$, the macroscopic number of
$X$ particles at the stationary state takes the same value as in system (2),
that is $X_{s}=1000$. Furthermore, the stationary state entropy production
$\sigma_{s}$ is strictly positive. For instance, in ideal systems, we find

$\sigma_{s}\,=\,\frac{5}{9}\,N\,k_{B}\,k_{1}\,a^{2}\,\ln(9/2)\,\,>\,\,0$ (5)

so that, as before, the system operates under non-equilibrium conditions.

For the microscopic simulation of (4), the total number of hard spheres is set
to $N=7000$ with the same number density as in system (2), i.e., $3\times
10^{-3}$ per $d^{3}$. The other parameters are as follows: $A=1500$ particles
(thus $B=2500$ and $C=500$ particles) and $k_{1}=0.9\times\nu$ (recall that
$k_{1}=k_{2}=k_{-1}=k_{-2}$).

As stated above, our purpose in performing these simulations is to estimate
the probability distributions $P(X,t\,;X^{\prime},t+\tau)$ and
$P(X^{\prime},t\,;X,t+\tau)$. To this end, successive values of $X$, separated
by a relatively small time interval $\tau$, were recorded for statistical
analysis. For each reactive system, (2) and (4), we collected a total of
$10^{9}$ sample points (SP), after the stationary regime has been reached (at
about $10^{4}$ SP).

The results are presented in Figure 1. This figure shows
$P(X,t\,\,;\,X_{ref},t+\tau)$ and $P(X_{ref},t\,\,;\,X,t+\tau)$ as a function
of $X$, for both systems (2) and (4), where $X_{ref}=1000$ (i.e., macroscopic
stationary state value in both systems) and $\tau=10$ mean reactive collision
time (MRCT). The statistical error, estimated from 10 successive runs of
$10^{8}$ SP, does not exceed $1\%$ (less than the size of the marker symbols
in Figure 1). While the sample path associated with model (4) is distinctly
time-irreversible, as expected, this is not the case of that associated with
model (2) where $P(X,t\,\,;\,X_{ref},t+\tau)\approx
P(X_{ref},t\,\,;\,X,t+\tau)$. This observation is further confirmed in Figure
2 where the ratio $P(X,t\,;\,X_{ref},t+\tau)/P(X_{ref},t\,;\,X,t+\tau)$ is
depicted for both models.

The time-reversible symmetry of a sample path is the key signature of systems
at thermodynamic equilibrium. But the system (2) operates under strict non-
equilibrium conditions (cf. eq. (3)). Consequently, the thermodynamic
properties of this reactive system based on its sample paths will necessarily
lead to a zero entropy production at the stationary-state, in obvious
contradiction with assertions presented in the Comment [8]. The situation is
different for the case of the second model (4) since in this reactive system
each elementary reaction leads to a different change in composition. The
microscopic simulation results thus confirm perfectly the validity of our
results [3, 4].

Figure 1: Joint probability distributions $P(X,t\,\,;\,X_{ref},t+\tau)$
(crosses) and $P(X_{ref},t\,\,;\,X,t+\tau)$ (bullets), as a function of $X$,
with $X_{ref}=1000$ (macroscopic stationary state) and $\tau=10$ MRCT. Top:
model (2). Bottom: model (4). Note that in both cases the entropy production
$\sigma_{s}>0$.

Figure 2: Probability ratio
$P(X,t\,\,;\,X_{ref},t+\tau)\,/\,P(X_{ref},t\,\,;\,X,t+\tau)$ (see the caption
of Fig. 1)

##  3 Concluding remarks and perspectives

The study of the statistical properties of reactive systems is traditionally
based on jump Markov process type of modeling, introduced in the mid-sixties
[17]. As shown in [3, 4], this traditional stochastic modeling limits the
domain of applicability of the associated path thermodynamics to reactive
systems involving only a single elementary reaction at the origin of each
observed change in composition. An alternative modeling that straightened this
limitation was proposed recently [8]. Using a microscopic simulation of dilute
reactive systems, in the Boltzmann limit, we showed in the present work that
results obtained through this procedure are in excellent agreement with the
predictions of the traditional modeling.

A peculiar consequence of this result concerns the class of one variable
reactive systems in which all reactions lead either to the change in
composition $X\rightarrow X+1$ (forward) or $X\rightarrow X-1$ (backward). A
well-known example is the Schlögl model [30]. The model (2) considered in our
microscopic simulations is precisely of this type. With the chosen parameter
values for this model, the entropy production is strictly positive, which
guarantees that the system operates under non-equilibrium conditions (cf. eq.
(3)). However, the associated observed state trajectory proves to be time-
reversible, in the sense that a sample path joining an arbitrary state
$\Gamma_{1}$ to another arbitrary state $\Gamma_{\\!2}$ will occur with the
same probability as the corresponding reverse path joining $\Gamma_{\\!2}$ to
$\Gamma_{1}$. But such a time-reversal symmetry is the key signature of
thermodynamic equilibrium state where, on average, each forward reaction is
exactly balanced by its reverse. That is not the case here. In addition, as
shown previously [4, 3], upon restricting ourselves to a well stirred system
(no diffusion) and adopting a traditional Markovian modeling, the results
observed in the simulation can be proved rigorously (see Section 6.3 in [16]
for more details). We are thus faced with a strange paradox.

Meticulous readers may object to our conclusions by pointing out that they may
just result from an over simplified theoretical modeling of the system. And
they are right. In particular, a problematic feature concerns the "perfect
homogeneity" assumption. Keeping a reactive system out of equilibrium requires
fixing the concentration of some chemically active components to prescribed
values. Not only do these chemical intermediates fluctuate locally because of
reactive collisions, even in well stirred systems, but they also diffuse. Even
though we are only interested in global (space averaged) quantities, there is
a priori no guaranty that the effect of these local fluctuations cancels out
through space averaging. In other words, local fluctuations could compromise
the "perfect homogeneity" assumption.

Here again, in the absence of sufficiently accurate experimental results, the
only way to address this issue is through microscopic simulations of reactive
fluids. For evident efficiency reasons, in the present article we used Bird’s
algorithm for microscopic simulations. Not only this algorithm is up to 3
orders of magnitude faster than the traditional hard spheres molecular
dynamics, but in addition it allows the simulation of perfectly homogeneous
(zero dimensional) Boltzmann equation. This choice was motivated by our main
purpose: investigate the validity of different types of stochastic modeling of
a perfectly homogeneous reactive systems. The next step is to perform
microscopic simulations of reactive fluids using the exact hard spheres
molecular dynamics procedure. Work in this direction is in progress.

Finally, while numerical simulations are useful one should not conflate a
computational algorithm with a mathematical process. For example, it is
misleading (and incorrect) to define a jump Markov process as being equivalent
to the Gillespie algorithm [34, 35] (see eqs. (2) and (3) in [8]). Such a
process is, in fact, defined by the "transition rate" concept, introduced by
Kolmogorov [13], while the algorithm is simply a numerical procedure to
generate the so-called "minimal process" associated with jump Markov processes
[36].

## Acknowledgments

The use of computational facilities at the Computing Center of the University
of Bourgogne, DNUM-CCUB, is gratefully acknowledged. The authors thank Carmela
Chateau-Smith for the careful reading of the manuscript. One author (AG)
acknowledges support by the U.S. Department of Energy, Office of Science,
Office of Advanced Scientific Computing Research, Applied Mathematics Program
under contract No. DE-AC02-05CH11231.

## Appendix A: Direct Simulation Monte Carlo for Reactive Systems

As in molecular dynamic simulation (MD), the state of the system in Direct
Simulation Monte Carlo (DSMC) is the set of particle positions and velocities,
$\\{{\bf r}_{i},{\bf v}_{i}\\}$. The evolution equations are integrated over
successive time steps $\Delta t$, typically a fraction of the mean collision
time for a particle. Within a time step, the free flight motion and the
particle interactions (collisions) are assumed to be decoupled. The free
flight motion for a particle $i$ is trivially computed as ${\bf
r}_{i}(t+\Delta t)={\bf r}_{i}(t)+{\bf v}_{i}(t)\Delta t$. After all particles
have been moved, they are sorted into "collisional" cells, typically a
fraction of mean free path in length. The main hypothesis in Bird’s DSMC
algorithm is that the cells are assumed to be perfectly homogeneous, i.e., all
particles within a cell are considered to be potential collision partners,
regardless of their exact positions. This basic hypothesis simplifies
considerably the dynamics and allows the algorithm to be up to three orders of
magnitude faster than the corresponding exact hard sphere MD.

A set of representative collisions, for the time step $\Delta t$, are then
chosen in each cell. A collision probability is assigned to each selected pair
based on their relative speed; a random impact parameter is selected and the
collision is performed. After the collision process has been completed in each
cell, the particles are moved according to their updated velocities and the
procedure is repeated. At this point we may recall that our main purpose here
is the study of the statistical properties of perfectly homogeneous reactive
systems. The DSMC algorithm is particularly well adapted for this case since
it allows the simulation of a homogeneous Boltzmann gas simply by associating
the entire system volume to a single collisional cell.

Reactions are modeled using "hard sphere chemistry" which was introduced in
the mid-1970’s [22, 24]. The basic idea is quite elegant and simple. We first
assign to each species an attribute, say a "color". A reactive collision
occurs if the colliding particles have "enough" energy, i.e., if the relative
kinetic energy of the colliding particles exceeds some threshold related to
the activation energy of the reaction [37]. If this is the case, then the
colors of the particles are changed, according to the chemical step under
consideration.

A major problem with this procedure is that it leads to the deformation of the
Maxwell-Boltzmann distribution since only the most energetic particles can
actually undergo a reactive transformation [38, 25]. To avoid this non-
equilibrium effect, the frequency of reactive collisions must be significantly
smaller than the frequency of elastic collisions, which results in an
significant waste of CPU time. One way to overcome this difficulty is to
further simplify the reactive collision rules by the following procedure. If
the intensive quantities are expressed through mole fractions, instead of
concentrations, then the kinetic constants are proportional to the collision
frequency, that is

$k_{i}\,=\,\nu_{i}\,\exp{\big{\\{}-\,E_{i}/k_{B}\,T\,\big{\\}}}\,\equiv\,\nu_{i}\,\widetilde{k}_{i}$
(6)

where $E_{i}$ is the activation energy of the reaction $i$ and $\nu_{i}$ is
the collision frequency between the corresponding reactant particles. After a
collision between two such reactive particles has occurred, we choose randomly
$\widetilde{k}_{i}\%$ of the collisions to be reactive, where
$\widetilde{k}_{i}$ stands for the Arrhenius factor defined in eq. (6). Note
that since the first reverse reaction in (2) involves a pair of the same
particles (i.e., $X$ particles), the relation (6) must be replaced by
$k_{-1}=\nu_{-1}\,\widetilde{k}_{-1}/2$ for that reaction. In any case, this
procedure avoids the deformation of the Maxwell-Boltzmann distribution since
it does not involve any systematic energy transfer between reactants and
products. It is, however, restricted to isothermal second-order (binary
collisions) reactions (see [29] for a review).

A final issue concerns the appropriate microscopic procedure to maintain
constant the mole fractions of some of the chemical species. To this end, in
addition to the reactive chemicals, we also consider "solvent" particles $S$.
Their role is precisely to maintain constant the number of $A$, $B$ and $C$
particles through the following strategy. Each time one of these particles is
created through a reactive collision, it is replaced by an $S$ particle.
Similarly, when one of these particles is destroyed through a reactive
collision, an $S$ particle is chosen randomly and transformed into that
species. Since the solvent particles don’t intervene directly in the reaction
scheme (2), they don’t modify the system’s dynamics while maintaining it out
of equilibrium [39]. The very same way as the presence of particle reservoirs
do, for example in CSTR (well stirred tank reactor). Finally, note that for
both reaction models (2) and (4) the number of $A$, $B$, and $C$ particles and
the sum of $X$ and solvent particles $X(t)+S(t)$ remain constant. As such,
knowledge of $X(t)$ determines entirely the state of the system at each
instant of time.

## References

  * [1] D. Kondepudi, Introduction to Modern Thermodynamics, Wiley (2008).
  * [2] A.L. Garcia, Essentials of Modern Thermodynamics, Amazon (2022).
  * [3] M. Malek Mansour and A. L. Garcia, Phys. Rev E, 101, 052135 (2020).
  * [4] M. Malek Mansour and F. Baras, Chaos 27, 104609 (2017).
  * [5] Y. De Decker, A. Garcia Cantu Ros and G. Nicolis, Euro. Phys. J. 224, 947 (2015) ; Y. De Decker, J-F. Derivaux and G. Nicolis, Phys. Rev. E 93, 042127 (2016).
  * [6] D. Seleznev, G. A. Zhernokleeva and L. M. Martyushe, JETP Letters 102, 557 (2015).
  * [7] J. Kurchan, J. Phys. A 31, 3719 (1998) ; J. Stat. Mech., 07, P07005 (2007).
  * [8] P. Gaspard, Phys. Rev E 103, 016101 (2021).
  * [9] D. Andrieux and P. Gaspard, J. Stat. Phys. 127, 107 (2007).
  * [10] J. L. Lebowitz and H. Spohn, J. Stat. Phys., 95, 333 (1999).
  * [11] P. Gaspard, J. Chem. Phys. 120, 8898 (2004).
  * [12] D. Andrieux and P. Gaspard, J. Chem. Phys. 121, 6167 (2004).
  * [13] See for example, B. Øksendal, Stochastic Differential Equations: An Introduction with Applications. Springer (2003).
  * [14] S. N. Ethier and T. G. Kurtz, Marko Processes: characterization and Convergence, Wiley- Interscience (2009).
  * [15] N. G. Van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam (1983).
  * [16] C.W. Gardiner, Handbook of Stochastic Methods, Springer-Verlag (2009).
  * [17] D. McQuarrie, Suppl. Rev. Ser: Appl. Prob., Methuen, London (1967).
  * [18] H. Haken, Synergetics: An Introduction, Springer-Verlag, Berlin, (1983).
  * [19] G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems, Wiley- Interscience (1977).
  * [20] T. G. Kurtz, Math. Progr. Stud. 5, 67 (1976) ; Stoch. Proc. Appl. 6, 223 (1978) ; M. Malek Mansour, C. Van Den Broeck, G. Nicolis and J.W. Turner, Ann. Phys. (USA) 131, 283 (1981).
  * [21] U. Seifert, Phys. Rev. Lett. 95, 040602 (2005) ; Eur. Phys. J. B 64, 423 (2008).
  * [22] J. Portnow, Phys. Lett. A, 51, 370 (1975).
  * [23] P. Ortoleva and S. Yip, J. Chem. Phys., 65, 2045 (1976).
  * [24] J. Boissonade, Phys. Lett. A 74, 285 (1979) ; Physica A 113, 607 (1982).
  * [25] F. Baras and M. Malek Mansour, Phys. Rev. Lett., 63, 2429 (1989); M. Malek Mansour and F. Baras, Physica A, 188, 253 (1992).
  * [26] A. Lemarchand and B. Nowakowski, J. Chem. Phys. 111, 6190 (1999) ; EuroPhys. Lett. 94, 48004 (2011).
  * [27] J. S. Hansen, B. Nowakowski and A. Lemarchand, J. Chem. Phys. 124, 034503 (2006) ; P. Dziekan, L. Signon, B. Nowakowski, and A. Lemarchand, J. Chem. Phys. 139, 114107 (2013).
  * [28] F. Baras, M. Salazar, E. Kestemont and M. Malek Mansour, Europhys. Lett., 67 (6), 900-906 (2004).
  * [29] F. Baras and M. Malek Mansour, Adv. Chem. Phys. 100, 393 (1997).
  * [30] F. Schlögl, Z. Phys. 248, 446 (1971) ; Z. Phys. 253, 147 (1972).
  * [31] R. Graham and H. Haken, Generalized thermodynamic potential for Markov system in detailed balance and far from thermal equilibrium, Z. Phys. 243 289 (1971).
  * [32] G.A. Bird, Molecular Gas Dynamics (Clarendon, Oxtord, 1976).
  * [33] A.L. Garcia, Numerical methods for physics, Englewood Cliffs, NJ: Prentice Hall, (2000).
  * [34] D. T. Gillespie, J. Comput. Phys. 22, 403 (1976) ; J. Phys. Chem. 81, 2340 (1977).
  * [35] D. T. Gillespie, Markov Processes: An Introduction of Physical Scientists (Academic, New York, 1992).
  * [36] S. Karlin and H. Taylor, A First Course in Stochastic Processes, (Academic Press, 1975).
  * [37] R. D. Present, J. Chem. Phys., 31, 747 (1959).
  * [38] I. Prigogine and E. Xhrouet, Physica, 15, 913 (1949).
  * [39] F. Baras, J. E. Pearson, and M. Malek Mansour, J. Chem. Phys., 93, 5747 (1990).


# Freezing-in dark matter through a heavy invisible $Z^{\prime}$

Gautam Bhattacharyyaa<EMAIL_ADDRESS>Maíra Dutrab
<EMAIL_ADDRESS>Yann Mambrinib<EMAIL_ADDRESS>Mathias
Pierreb<EMAIL_ADDRESS>a Saha Institute of Nuclear Physics, HBNI,
1/AF Bidhan Nagar, Kolkata 700064, India  b Laboratoire de Physique Théorique
(UMR8627), CNRS, Université Paris-Sud, Université Paris-Saclay, 91405 Orsay,
France

###### Abstract

We demonstrate that in a class of the $U(1)^{\prime}$ extension of the
Standard Model (SM), under which all the Standard Model matter fields are
uncharged and the additional neutral gauge boson $Z^{\prime}$ couples to a set
of heavy nonstandard fermions, dark matter (DM) production mediated by
$Z^{\prime}$ can proceed through the generation of generalized Chern-Simons
(GCS) couplings. The origin of the GCS terms is intimately connected to the
cancellation of gauge anomalies. We show that the DM production cross section
triggered by GCS couplings is sufficient even for an intermediate scale
$Z^{\prime}$ . A large range of DM and Z masses is then allowed for reasonably
high reheating temperature ($T_{\text{\tiny{RH}}}\gtrsim
10^{10}~{}\text{GeV}$). This type of scenario opens up a new paradigm for
unified models. We also study the UV completion of such effective field theory
constructions, augmenting it by a heavy fermionic spectrum. The latter, when
integrated out, generates the GCS-like terms and provides a new portal to the
dark sector. The presence of a number of derivative couplings in the GCS-like
operators induces a high temperature dependence to the DM production rate. The
mechanism has novel consequences and leads to a new reheating dependence of
the relic abundance.

††preprint: LPT–Orsay 18-14

##  I I. Introduction

In spite of a lot of speculations about its origin, dark matter (DM) still
remains an enigma, and the best we can do is to assume that it has a particle
physics origin in the domain of natural extensions of the Standard Model (SM).
However, the twin pressure of the clear existence of DM in the energy budget
of the Universe planck  and simultaneously the lack of any DM signal in direct
detection experiments XENON XENON , LUX LUX , and PANDAX PANDAX  pushes the
limits on weakly interacting massive particles (WIMPs) toward unnatural
corners of the parameter space. The simplest extensions as Higgs portal
Higgsportal , $Z$-portal Zportal , or even $Z^{\prime}$ portal Zpportal  are
now severely constrained (for a review on WIMP searches and supporting models,
see Ref. Arcadi:2017kky ). This scenario motivates the assumption that the
interactions between the dark and visible sectors are even weaker, leading to
an out-of-equilibrium production of feebly interacting massive particles
(FIMPs) fimp  (see Ref. Bernal:2017kxu  for a review). Alternative generation
mechanisms of a weakly interacting dark sector by direct thermal production at
the reheating temperatures are discussed in Ref. alternative .

On the other hand, theoretical considerations ranging from neutrino mass
generation mechanisms to grand unified theories (GUT), as well as inflation,
reheating, leptogenesis, or Higgs stability, all hint toward the existence of
an intermediate scale between $10^{10}$ and $10^{16}$ GeV. To interpret the
absence of DM signals, instead of invoking unnaturally weak DM-SM couplings,
one could explain its secluded nature by suppressions arising from the high-
mass scale of the mediators involved in its interaction with the thermal bath.
Concrete realizations for the high-energy physics origin of the freeze-in
mechanism involve Planck suppressed portals that could be embedded in quantum
theories of gravity mambrinilast ; Garny:2015sjg ; PIDM , the left-right
symmetric model Biswas:2018aib  and $Z^{\prime}$ mediators in the $SO(10)$
framework Mambrini:2013iaa ; 2loops .

In usual $Z^{\prime}$-mediated constructions, the SM particles are charged
under the new gauge group $U(1)^{\prime}$. An interesting question, therefore,
is to ask if the DM production processes are still efficient even if the SM is
uncharged under $U(1)^{\prime}$. To generate the effective interaction of the
associated gauge boson $Z^{\prime}$ with the SM fields, we would need a set of
nonstandard fermions charged under $U(1)^{\prime}$ as well as under the SM
gauge group(s). Such set-up is quite common in string constructions, or in
$E_{6}$ models. In this case, effective interactions of the type represented
by the Lagrangian

${\cal
L}\supset\lambda~{}\epsilon^{\mu\nu\rho\sigma}Z^{\prime}_{\mu}A_{\nu}F_{\rho\sigma}\,,$
(1)

where $F_{\rho\sigma}=\partial_{\rho}A_{\sigma}-\partial_{\sigma}A_{\rho}$,
arise from diagrams leading to anomalies (à la Green-Schwarz in string models
or Peccei-Quinn in the presence of axionic couplings). The gauge boson $A$ in
the above expression may be the SM hypercharge gauge boson or could even be
any other nonstandard $U(1)$ gauge boson.

Such terms, characterized by the presence of three gauge bosons and one
derivative, are dubbed generalized Chern-Simons (GCS) terms GCS  (see also
Ref. anomaly-zpr  for a general discussion of anomaly-free $Z^{\prime}$
models). This operator can be generated at the dimension-4 level as a low-
energy effective term by integrating out a set of heavy fermions. The
effective coupling $\lambda$ is independent of the heavy fermion masses, and
hence its effect does not decouple by increasing those fermion masses. The
underlying dynamics behind this apparent nondecoupling is that this term is
gauge noninvariant, but its gauge variation cancels some triangle anomalies
from some lighter fermions still persisting in the theory. On the other hand,
when we integrate out a complete anomaly-free set of heavy fermions having
couplings with $Z^{\prime}$ and $A$, similar but higher-dimensional terms with
three gauge bosons and more derivatives are generated, albeit suppressed by
the heavy fermion mass scale. We will indistinctly refer to those also as GCS
terms as in Eq. (1). Such terms have already been studied as a connection
between the SM field content and dark sectors in thermal CSDM  and nonthermal
Farzan:2014foo  dark matter production mechanisms.

Now, let us suppose that $Z^{\prime}$ couples to our DM candidate, leaving
open the possibilities that this candidate could be a fermion or a vectorial
boson (Abelian or non-Abelian). If the DM is a vectorial boson, we also assume
that it couples to $Z^{\prime}$ through GCS interactions generated possibly by
a different set of heavy fermions. Notice that, even if the SM matter fields
are not charged under $U(1)^{\prime}$, $Z^{\prime}$ would still be produced in
the early Universe through the freeze-in mechanism fimp , thanks to the GCS
terms, which would then decay to DM particle $\chi$ : (SM) (SM) $\rightarrow
Z^{\prime}\rightarrow\chi\chi$. However, the DM production rate is doubly
suppressed, both by the mass scale of heavy fermions in loop generating the
GCS couplings as well as by the mass of the virtual $Z^{\prime}$ exchanged in
the process. In fact, through a large temperature dependence, this rate is
highly sensitive to the reheating phase, especially when considered
noninstantaneous, as shown in Refs. Giudice:2000ex ; Garcia:2017tuj  (and Ref.
gravitino  for gravitino dark matter). Moreover, allowing for vectorial dark
matter (Abelian or non-Abelian) brings in higher-derivative couplings in the
GCS operators, thereby inducing a significantly high temperature sensitivity
to the DM production rate.

The paper is organized as follows. We first describe the model under
consideration. Then, we compute the relic abundance from the freeze-in process
in the early thermal bath, taking into account noninstantaneous reheating. We
then construct models for UV completion that would naturally lead to our
framework. We finally conclude, highlighting the new aspects that emerged from
our analysis.

##  II II. Our model

The two portals that connect the DM to the SM sector are $(i)$ a $Z^{\prime}$
gauge boson of the $U(1)^{\prime}$ group and $(ii)$ a set of heavy fermions
charged under $U(1)^{\prime}$ as well as the SM gauge group(s). Depending on
the nature of the couplings of those heavy fermions with the SM gauge bosons
and the $Z^{\prime}$, three possibilities emerge

  * •

The heavy fermions are vectorlike with respect to SM gauge bosons as well as
$U(1)^{\prime}$. In this case, no GCS terms are generated, as in the absence
of chiral couplings, any potential operator that contains three gauge bosons
(and derivatives) vanishes identically Dudas:2013sia .

  * •

The heavy fermions are chiral under $U(1)^{\prime}$ but do not form an
anomaly-free set. Then, GCS terms at dimension-4 level are generated in the
effective action of which the gauge variations exactly cancel the anomaly GCS
, yielding

${\cal
L}=\alpha~{}\theta^{\prime}\text{Tr}[F^{\text{\tiny{SM}}}F^{\text{\tiny{SM}}}]+\beta~{}\epsilon^{\mu\nu\rho\sigma}Z^{\prime}_{\mu}A^{\text{\tiny{SM}}}_{\nu}F^{\text{\tiny{SM}}}_{\rho\sigma}\,,$
(2)

with $A^{\text{\tiny{SM}}}$ and $F^{\text{\tiny{SM}}}$ being the hypercharge
SM gauge boson and associated field strength, respectively, and
$\theta^{\prime}$ being the would-be Goldstone boson of $U(1)^{\prime}$.
Notice the absence of any suppression coming from the masses of heavy fermions
(generically denoted by $M$) that are integrated out in generating the
effective couplings $\alpha$ and $\beta$. This apparent nondecoupling
signifies that the set-up is not anomaly free. These constructions are common
in string-inspired models (Green-Schwarz mechanism) and lead to interesting
phenomenological consequences (e.g., advocated to justify monochromatic
spectral lines originating from the dark sector line ).

  * •

The heavy fermions form an anomaly-free set. They are chiral under
$U(1)^{\prime}$ but vectorlike with respect to the SM gauge group. In this
case, several gauge-invariant combinations can be written. The complete list
can be found in Ref. Dudas:2013sia . However, a lot of them111For instance,
terms of the type
$(1/M^{2})\epsilon^{\mu\nu\rho\sigma}Z^{\prime}_{\mu}(H^{\dagger}D_{\nu}H)F^{\text{SM}}_{\rho\sigma}$
do not contribute to CS-like couplings above the electroweak phase transition.
either are not generated by triangle loops or vanish when the SM gauge bosons
are on shell. In the limit where the field $\Phi$ that breaks the extra
$U(1)^{\prime}$ is much heavier than its vacuum expectation value (VEV), $V$,
which controls the $Z^{\prime}$ mass, the effective theory exhibits only the
axionic (longitudinal) component ($a$) of the field as
$\Phi=\frac{V+\phi}{\sqrt{2}}\exp(ia/V)$. Defining the dimensionless axion
$\theta^{\prime}\equiv\frac{a}{V}$, the only relevant coupling can then be
extracted Dudas:2013sia  from a gauge invariant Lagrangian, given by

${\cal L}=\frac{1}{M^{2}}\partial^{\alpha}{\cal
D}_{\alpha}\theta^{\prime}~{}\epsilon^{\mu\nu\rho\sigma}\mathrm{Tr}[F^{\text{\tiny{SM}}}_{\mu\nu}F^{\text{\tiny{SM}}}_{\rho\sigma}]+V^{2}|{\cal
D}\theta^{\prime}|^{2}$ (3)

with ${\cal
D}_{\alpha}\theta^{\prime}\equiv\partial_{\alpha}\theta^{\prime}-\frac{\tilde{g}}{2}Z^{\prime}_{\alpha}$,
$\tilde{g}$ being the gauge coupling associated to the extra $U(1)^{\prime}$.
Throughout the present work, we shall consider this particular set-up, even if
our results can be applied to a general class of GCS couplings just by a
redefinition of parameters.

In the unitary gauge, the term related to the $Z^{\prime}$-SM-SM vertex can be
extracted from Eq. (3) as

${\cal
L}=\frac{\tilde{g}}{M^{2}}\partial^{\alpha}Z^{\prime}_{\alpha}\epsilon^{\mu\nu\rho\sigma}\partial_{\mu}A^{a}_{\nu}\partial_{\rho}A^{a}_{\sigma}\,,$
(4)

where $A^{a}$ are the SM gauge bosons. From now on, without any loss of
generality, we consider the gluons as the gauge bosons appearing in Eq. (4)
and define $\frac{1}{\Lambda^{2}}\equiv\frac{\tilde{g}}{M^{2}}$, as the
results would be exactly the same with electroweak gauge bosons, just by
rescaling the couplings. We consider the heavy fermions generating the GCS
couplings to be charged under $SU(3)_{C}$ so that they dominate the production
process222Considering fermions without hypercharge $Y$ leads to a
simplification as kinetic mixing of the type
$\delta~{}Z^{\prime}_{\mu\nu}B^{\mu\nu}$ can be avoided. Effects of such
mixing have been extensively studied in the literature kineticmix . With this
approach, the relevant Lagrangian would then read Dudas:2013sia

$\mathcal{L}_{\text{eff}}=\frac{1}{\Lambda^{2}}\partial^{\alpha}Z^{\prime}_{\alpha}\epsilon^{\mu\nu\rho\sigma}\text{Tr}[G^{a}_{\mu\nu}G^{a}_{\rho\sigma}]+\mathcal{L}_{\text{\tiny{DM}}}^{i}\,,$
(5)

where ${\cal L}^{i}_{\text{\tiny{DM}}}$ represents the interactions between
the $Z^{\prime}$ and the DM candidate, which can be fermionic ($\chi$), or
vectorial of Abelian ($X_{1}$) or non-Abelian ($X_{N}$) types. The respective
Lagrangians are given by333Only axial coupling is present for the fermionic
dark matter. The derivative $\partial^{\alpha}$ before $Z^{\prime}_{\alpha}$
in Eq. (4) ensures that the vector coupling does not contribute in a
$GG\rightarrow\chi\chi$ process.

$\mathcal{L}_{\text{\tiny{DM}}}^{\chi}=\alpha~{}\bar{\chi}\gamma^{\mu}\gamma_{5}\chi
Z^{\prime}_{\mu},$ (6)
$\mathcal{L}_{\text{\tiny{DM}}}^{\text{\tiny{$X_{1}$}}}=\beta~{}\epsilon_{\mu\nu\rho\sigma}Z^{\prime\mu}X_{1}^{\nu}X_{1}^{\rho\sigma},$
(7)

and

$\mathcal{L}_{\text{\tiny{DM}}}^{\text{\tiny{$X_{N}$}}}=\gamma~{}\partial^{\alpha}Z^{\prime}_{\alpha}\epsilon_{\mu\nu\rho\sigma}\text{Tr}[X_{N}^{\mu\nu}X_{N}^{\rho\sigma}].$
(8)

##  III III. Results

The evolution of dark matter number density $n_{\text{\tiny{DM}}}$ is governed
by the Boltzmann equation

$\frac{\text{d}n_{\text{\tiny{DM}}}}{\text{d}t}=-3H(T)n_{\text{\tiny{DM}}}+R(T),$
(9)

where $H\left(T\right)$ is the Hubble expansion rate and $R(T)$ is the
temperature-dependent interaction rate. In the regime in which the abundance
of dark matter is much smaller than the abundance of particles in the thermal
bath, the backreaction term in the rate (dark matter producing standard
particles) may be neglected, which is usually the case in the freeze-in
mechanism444The correct amount of dark matter is generated in a regime in
which $n_{\text{\tiny{DM}}}\ll n_{\text{\tiny{SM}}}$, since
$\frac{\Omega_{\text{\tiny{DM}}}^{0}h^{2}}{0.12}\sim\frac{Y_{\text{\tiny{DM}}}}{10^{-10}}\frac{m_{\text{\tiny{DM}}}}{\,{\rm
GeV}}$ and $Y_{\text{\tiny{DM}}}\propto
n_{\text{\tiny{DM}}}/n_{\text{\tiny{SM}}}$..

Figure 1: Production of dark matter through gluon fusion in the early Universe

In our set-up, the freeze-in occurs through the process depicted in Fig. 1.
For a 1 + 2 $\rightarrow$ 3 + 4 process, the rate can be written as
$R(T)=n_{\text{sm}}^{2}\langle\sigma v\rangle$, where $n_{\text{sm}}$ is the
number density of the SM species and $\langle\sigma v\rangle$ is the thermal
averaged production cross section. For a dark matter particle $i$, the rate
reads

$R(T)=\int
f_{1}f_{2}\frac{E_{1}E_{2}\text{d}E_{1}\text{d}E_{2}~{}\text{d}\cos\theta_{12}}{1024\pi^{6}}\int|{\cal
M}|_{i}^{2}\text{d}\Omega_{13}\,,$ (10)

with $E_{1,2}$ and $f_{1,2}$ as the energy and the distribution function of
the initial SM particles, $\theta_{12}$ as the angle between them, and
$d\Omega_{13}$ as the solid angle between particles 1 and 3.

For the fermionic dark matter case,

$\begin{split}\int|\mathcal{M}|_{\chi}^{2}~{}\text{d}\Omega_{13}&=2^{10}\pi~{}\frac{\alpha^{2}}{\Lambda^{4}}\frac{m_{\chi}^{2}}{M_{\text{$Z^{\prime}$}}^{4}}\frac{s^{3}(s-M_{\text{$Z^{\prime}$}}^{2})^{2}}{(s-M_{\text{$Z^{\prime}$}}^{2})^{2}+M_{\text{$Z^{\prime}$}}^{2}\Gamma_{\text{$Z^{\prime}$}}^{2}}\\\
&\approx
2^{10}\pi~{}\frac{\alpha^{2}}{\Lambda^{4}}\frac{m_{\chi}^{2}}{M_{\text{$Z^{\prime}$}}^{4}}s^{3}\end{split}$
(11)

For the Abelian dark matter case,

$\begin{split}\int|\mathcal{M}|_{\text{\tiny{$X_{1}$}}}^{2}~{}\text{d}\Omega_{13}&=2^{10}\pi~{}\frac{\beta^{2}}{\Lambda^{4}}\frac{s^{3}}{M_{\text{$Z^{\prime}$}}^{4}}\frac{(s-4m_{\text{\tiny{$X_{1}$}}}^{2})(s-M_{\text{$Z^{\prime}$}}^{2})^{2}}{(s-M_{\text{$Z^{\prime}$}}^{2})^{2}+M_{\text{$Z^{\prime}$}}^{2}\Gamma_{\text{$Z^{\prime}$}}^{2}}\\\
&\approx
2^{10}\pi~{}\frac{\beta^{2}}{\Lambda^{4}}\frac{1}{M_{\text{$Z^{\prime}$}}^{4}}s^{4}\end{split}$
(12)

For the non-Abelian dark matter case,

$\begin{split}\int|\mathcal{M}|_{\text{\tiny{$X_{N}$}}}^{2}\text{d}\Omega_{13}&=2^{12}\pi~{}\frac{\gamma^{2}}{\Lambda^{4}}\frac{s^{5}}{M_{\text{$Z^{\prime}$}}^{4}}\frac{(s-4m_{\text{\tiny{$X_{N}$}}}^{2})(s-M_{\text{$Z^{\prime}$}}^{2})^{2}}{(s-M_{\text{$Z^{\prime}$}}^{2})^{2}+M_{\text{$Z^{\prime}$}}^{2}\Gamma_{\text{$Z^{\prime}$}}^{2}}\\\
&\approx
2^{12}\pi~{}\frac{\gamma^{2}}{\Lambda^{4}}\frac{1}{M_{\text{$Z^{\prime}$}}^{4}}s^{6}\end{split}$
(13)

Above, $\Gamma_{\text{$Z^{\prime}$}}$ is the total width of $Z^{\prime}$ (see
the Appendix for details); $s$ is the center-of-mass energy squared; and
$m_{\chi},m_{\text{\tiny{$X_{1}$}}},m_{\text{\tiny{$X_{N}$}}}$, and
$M_{\text{$Z^{\prime}$}}$ are the three types of dark matter and $Z^{\prime}$
masses, respectively. Note that we recover the Landau-Yang effect in the above
expressions, though the pole enhancement studied in Ref. mambrinilast  is not
present in our case. Note also that the vectorial nature of the mediator has
specific characteristics that we do not observe for other type of mediators.
Importantly, we notice that once the pole is reached
($s=M_{\text{$Z^{\prime}$}}^{2}$), the production rate vanishes exactly – see
Eqs. (11), (12), and (13). This is expected following the Landau-Yang theorem,
which states that a massive spin-1 particle cannot decay into two massless
spin-1 fields. This behavior is opposite of the traditional freeze-in scenario
in which, on the contrary, the majority of dark matter is produced when the
temperature of the thermal bath reaches $T\sim M_{\text{$Z^{\prime}$}}$
mambrinilast ; Blennow:2013jba .

We have integrated numerically the production rate, Eq. (10), considering the
Bose-Einstein distributions of the gluons and the exact squared amplitudes of
our three dark matter candidates. Our result is depicted in Fig.2.

We can obtain analytical approximations for the rates by assuming
$\Gamma_{\text{$Z^{\prime}$}}\ll M_{\text{$Z^{\prime}$}}$ and
$m_{\text{\tiny{DM}}}^{2}\ll s$:

$R(T)\approx\left\\{\begin{array}[]{rc}2\times
10^{2}~{}\dfrac{\alpha^{2}}{\Lambda^{4}}\dfrac{m_{\chi}^{2}}{M_{\text{$Z^{\prime}$}}^{4}}T^{10}~{}~{}~{}~{}(\text{fermionic
DM})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
10^{4}~{}\dfrac{\beta^{2}}{\Lambda^{4}M_{\text{$Z^{\prime}$}}^{4}}T^{12}~{}~{}~{}~{}~{}~{}(\text{Abelian
DM})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 2\times
10^{9}\dfrac{\gamma^{2}}{\Lambda^{4}M_{Z^{\prime}}^{4}}T^{16}~{}(\text{non-
Abelian DM})\end{array}\right.$ (14)

We also show in Fig. 2 our approximate solutions. In the inset of the figure,
we show when they depart from the exact solutions. We can distinguish two
regimes in which the approximations fail. First, let us consider when the
temperature of the thermal bath is close to the mediator mass. In this case,
the exact solutions are smaller than the approximate results as an effect of
the nonvanishing mediator decay width. Even though the departure from
approximations is small in this case, it carries a special feature of our set-
up, emerging from the consequence of the Landau-Yang theorem in a thermal bath
of gluons. The significant departure from approximations occurs at large
$x\equiv M_{Z^{\prime}}/T$, due to a threshold effect, as for $T\ll
m_{\text{\tiny{DM}}}$ the production rate is exponentially suppressed because
only the high-energy tail of the initial states distribution function have
sufficient energy to produce a DM pair, an effect that is not encapsulated in
the analytical approximations.

Another typical characteristic of a longitudinal (“would-be Goldstone”)
mediator appears in the generic expression for the rate. Indeed, the “light”
mediator regime ($M_{\text{$Z^{\prime}$}}\ll T_{\text{\tiny{RH}}}$), and the
“heavy” mediator regime ($M_{\text{$Z^{\prime}$}}\gg T_{\text{\tiny{RH}}}$)
give the same dependence of the rate $R(T)$ on $M_{\text{$Z^{\prime}$}}$, and
thus on temperature for a given nature of dark matter, as one can see from the
Eq. (14) and Fig. 2. In fact, there exists only one main regime, independent
of the mass of the $Z^{\prime}$ mediator555This is in contrast with what has
been observed in mambrinilast  for spin-2 mediator. for which the slope of the
rate is constant until $T\sim m_{\text{\tiny{DM}}}$.

This can be understood by noting that only the longitudinal mode of
$Z^{\prime}$ is exchanged, and hence it cannot feel any pole effect. The
longitudinal component has its origin in the Goldstone mode of a nonlinear
sigma model. The behavior of the amplitude squared is dominated by a term
proportional to powers of $1/M_{\text{$Z^{\prime}$}}$. This happens because
the Goldstone, which is the dominant mode exchanged in the DM production
process, carries the $1/V$ factor arising from $U(1)^{\prime}$ breaking. This
is similar to the gravitino production in supergravity in which the
longitudinal mode, carrying a factor $1/m_{3/2}$ ($m_{3/2}$ being the
gravitino mass), is generated in a high-scale supersymmetric scenario as was
shown in Ref. gravitino .

Figure 2: Rates for the fermionic, Abelian and non-Abelian DM (orange, green
and blue lines, respectively). Solid lines represent the exact numerical
computation while dashed lines represent the approximated results based on
Eqs. (14). We fix $N_{\psi}=Q_{\psi}=1$ and
$m_{\psi}=0.4~{}M_{\text{$Z^{\prime}$}}$ and
$T_{\text{\tiny{MAX}}}=100~{}T_{\text{\tiny{RH}}}$ for illustrative purposes.

###  III.1 Dark matter freeze-in

For instantaneous reheating, the Universe is dominated by radiation and
entropy is conserved. In this case,

$\frac{\text{d}}{\text{d}t}=-H(T)T\frac{\text{d}}{\text{d}T},\hskip
14.22636pt\textrm{with}~{}~{}H(T)=\sqrt{\frac{g_{e}}{90}}\pi\frac{T^{2}}{M_{P}}$
(15)

and Eq. (9) can be put in the familiar form

$\frac{\text{d}Y}{\text{d}T}=-\frac{R(T)}{H~{}T~{}s}\,,\hskip
14.22636pt\textrm{for}~{}~{}T<T_{\text{\tiny{RH}}}$ (16)

where $Y=\frac{n}{s}$ is the dark matter yield,
$s=\frac{2\pi^{2}}{45}g_{s}T^{3}$ is the entropy density of the thermal bath
with $g_{e,s}$ the energy and entropy density degrees of freedom, and
$M_{P}\approx 2.4\times 10^{18}$ GeV is the reduced Planck mass.

However, once we consider noninstantaneous reheating, dark matter can be
produced before the end of reheating. Indeed, as was shown in Refs.
Giudice:2000ex ; Garcia:2017tuj , a rate with a dependence $R(T)\propto T^{n}$
for $n\geq 12$ enhances drastically the dark matter production before
$T_{\text{\tiny{RH}}}$ if the reheating is considered as noninstantaneous. For
instance, for $n=12$, the ratio of the relic abundance computed with the
noninstantaneous reheating hypothesis ($\Omega h^{2}$) to the one with the
instantaneous reheating hypothesis ($\Omega h^{2}_{\text{\tiny{RH}}}$) is
$\Omega h^{2}/\Omega h^{2}_{\text{\tiny{RH}}}\simeq
0.4\times\dfrac{55}{6}\ln\left(\dfrac{T_{\text{\tiny{MAX}}}}{T_{\text{\tiny{RH}}}}\right)$
(where $T_{\text{\tiny{MAX}}}$ is the maximum temperature produced in the
reheating process), whereas $\Omega h^{2}/\Omega
h^{2}_{\text{\tiny{RH}}}\simeq
0.4\times\dfrac{8}{5}\left(\dfrac{n-5}{n-12}\right)\left(\dfrac{T_{\text{\tiny{MAX}}}}{T_{\text{\tiny{RH}}}}\right)^{n-12}$
for $n>12$ and $\Omega h^{2}/\Omega h^{2}_{\text{\tiny{RH}}}\sim 2$ for
$n=10$. These ratios can be seen as “boost factors” from the reheating process
and will be called $B^{n}_{F}\equiv\Omega h^{2}/\Omega
h^{2}_{\text{\tiny{RH}}}$ from now on. After integrating Eq. (16) from
$T_{\text{\tiny{MAX}}}$ until the present day, we deduce the parameter space
leading to the relic abundance

$\begin{split}&\frac{\Omega
h^{2}}{0.12}\approx\left\\{\begin{array}[]{lc}\left(\dfrac{B_{F}^{10}}{2}\right)\left(\dfrac{\alpha}{1}\right)^{2}\left(\dfrac{m_{\chi}}{6\times
10^{10}~{}{\rm GeV}}\right)^{3}\left(\dfrac{10^{14}~{}{\rm
GeV}}{M_{\text{$Z^{\prime}$}}}\right)^{4}\left(\dfrac{10^{16}~{}{\rm
GeV}}{\Lambda}\right)^{4}\left(\dfrac{T_{\text{\tiny{RH}}}}{10^{12}~{}{\rm
GeV}}\right)^{5}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\left(\dfrac{B_{F}^{12}}{21}\right)\left(\dfrac{\beta}{10^{-5}}\right)^{2}\left(\dfrac{m_{\text{\tiny{$X_{1}$}}}}{10^{9}~{}{\rm
GeV}}\right)\left(\dfrac{10^{14}~{}{\rm
GeV}}{M_{\text{$Z^{\prime}$}}}\right)^{4}\left(\dfrac{10^{16}~{}{\rm
GeV}}{\Lambda}\right)^{4}\left(\dfrac{T_{\text{\tiny{RH}}}}{10^{13}~{}{\rm
GeV}}\right)^{7}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\left(\dfrac{B_{F}^{16}}{1.76\times
10^{8}}\right)\left(\dfrac{m_{\text{\tiny{$X_{N}$}}}}{2\times 10^{9}~{}{\rm
GeV}}\right)\left(\dfrac{10^{14}~{}{\rm
GeV}}{M_{\text{$Z^{\prime}$}}}\right)^{4}\left(\dfrac{10^{16}~{}{\rm
GeV}}{\Lambda}\right)^{8}\left(\dfrac{T_{\text{\tiny{RH}}}}{10^{12}~{}{\rm
GeV}}\right)^{11}\end{array}\right.\end{split}$ (17)

where we set $\gamma=1/\Lambda^{2}$, and we consider $SU(2)$ as gauge group
for the non-Abelian dark matter case.

Our results are summarized in Fig. 3, in which we plotted the parameter space
allowed by the cosmological constraints in the plane
$(m_{\text{\tiny{DM}}},\Lambda)$ for the three natures of dark matter
considered in this work and for $T_{\text{\tiny{RH}}}=10^{10}~{}\text{GeV}$
and $T_{\text{\tiny{MAX}}}=100~{}T_{\text{\tiny{RH}}}$. To compute the relic
abundance numerically, we did not assume instantaneous reheating. Instead,
using the method developed in Ref. Garcia:2017tuj , we assumed that an
inflaton $\phi$ decays into radiation with a rate $\Gamma_{\phi}$ and that the
DM particles are created and annihilated into radiation of density
$\rho_{\text{R}}$ with a thermal-averaged cross section times velocity
$\langle\sigma v\rangle$. The corresponding energy and number densities
satisfy the differential equations Chung:1998rq ; Giudice:2000ex

$\displaystyle\frac{\text{d}n_{\text{\tiny{DM}}}}{\text{d}t}$ $\displaystyle=$
$\displaystyle-3H\,n_{\text{\tiny{DM}}}-\langle\sigma
v\rangle\left[n_{\text{\tiny{DM}}}^{2}-(n_{\text{\tiny{DM}}}^{\text{eq}})^{2}\right]\,,$
$\displaystyle\frac{\text{d}\rho_{\text{R}}}{\text{d}t}$ $\displaystyle=$
$\displaystyle-4H\,\rho_{\text{R}}+\Gamma_{\phi}\,\rho_{\phi}+2\langle\sigma
v\rangle\langle
E_{\text{\tiny{DM}}}\rangle\left[n_{\text{\tiny{DM}}}^{2}-(n_{\text{\tiny{DM}}}^{\text{eq}})^{2}\right]\,,$
$\displaystyle\frac{\text{d}\rho_{\phi}}{\text{d}t}$ $\displaystyle=$
$\displaystyle-3H\,\rho_{\phi}-\Gamma_{\phi}\,\rho_{\phi}\,,$ (18)

with $\langle E_{\text{\tiny{DM}}}\rangle$ the mean energy of the dark matter.
The Hubble expansion parameter $H$ is given by
$H^{2}=\frac{1}{3M_{P}^{2}}(\rho_{\phi}+\rho_{\text{R}}+\rho_{\text{\tiny{DM}}})$,
and its dependence in temperature is quite complex between
$T_{\text{\tiny{MAX}}}$ and $T_{\text{\tiny{RH}}}$ because of the composition
of the Universe (mixed between a decaying inflaton compensated by an
increasing radiation666We have defined $T_{\text{\tiny{RH}}}$ as the
temperature in a radiation-dominated Universe after the inflaton decay,
$\Gamma_{\phi}=H(T_{\text{\tiny{RH}}})$).

Figure 3: Relic density curves for the fermionic (orange dashed line), Abelian
(green dashed line) and non-Abelian (blue dashed line) dark matter.

We can see that our analytical expressions in Eq. (17) give an impressively
good approximation. Indeed, the right amount of relic abundance planck  is
obtained (in the fermionic case) for $m_{\chi}=10^{6}$ GeV,
$T_{\text{\tiny{RH}}}=10^{10}$ GeV, $M_{\text{$Z^{\prime}$}}=10^{11}$ GeV, and
$\Lambda\simeq 8\times 10^{12}$ GeV, in perfect agreement with Fig. 3. The
slopes of the curves depicted in Fig. 3 correspond also perfectly with the
ones predicted by our analytical solution in Eq. (17): it follows a line
$\Lambda\propto m_{\text{\tiny{DM}}}^{3/4}$ ($m_{\text{\tiny{DM}}}^{1/4}$,
$m_{\text{\tiny{DM}}}^{1/8}$) for fermionic (Abelian, non-Abelian) for
$m_{\text{\tiny{DM}}}<T_{\text{\tiny{MAX}}}$.

Without entering too much into detail, there is an interesting feature in the
change of slope between $T_{\text{\tiny{RH}}}$ and $T_{\text{\tiny{MAX}}}$ in
the fermionic dark matter case. This is a novel feature that was not treated
in Ref. Garcia:2017tuj  nor Ref. mambrinilast . Indeed, in the case in which
dark matter is heavier than $T_{\text{\tiny{RH}}}$, there is still a
possibility to produce it as long as $m_{\text{\tiny{DM}}}\lesssim
T_{\text{\tiny{MAX}}}$. If the temperature dependence of the rate is small
enough (fermionic case), most of the DM density is produced at the lowest
scale available, and we notice a change of slope in the curve giving the
correct relic density. It is worth commenting that, due to statistical
distribution, the production rate does not vanish completely when $T\lesssim
m_{\text{\tiny{DM}}}$, which explains why the DM production window is still
open when $m_{\text{\tiny{DM}}}>T_{\text{\tiny{MAX}}}$777The corresponding
region of parameter space as shown in Fig. 3 is quantitatively less precise as
the EFT approach becomes less reliable.. Therefore, in this regime, a small
effective scale $\Lambda$ is required to compensate the thermal suppression of
the rate, as one can see in Fig.3.

Moreover, a quick look at Fig.3 shows to what extent the allowed parameter
space is technically natural. Indeed, for a very large range of the DM mass,
from $\mathcal{O}$(TeV) to $T_{\text{\tiny{RH}}}$, values of the beyond SM
scale $\Lambda$ range from $T_{\text{\tiny{RH}}}$ to GUT/string scale and can
still populate the Universe with the correct relic abundance. This means that
the heavy spectrum of masses above the reheating temperature
$T_{\text{\tiny{RH}}}$ generates naturally small couplings of an invisible
$Z^{\prime}$ to the SM bath to satisfy the cosmological constraints through
the freeze-in process. This constitutes one of the most important observations
of our work.

##  IV IV. toward a microscopic approach

As mentioned earlier, we consider processes happening at a temperature below
the $U(1)^{\prime}$ phase transition scale. We have also assumed that the
radial component of the complex scalar that breaks $U(1)^{\prime}$ is way too
heavy compared to the corresponding VEV ($V$). Then, $Z^{\prime}$ is primarily
longitudinal absorbing the axion field ($a$), and the effective Lagrangian
containing $Z^{\prime}$ realizes the gauge symmetry nonlinearly à la
Stueckelberg. Now, we attempt to look deep inside the effective GCS vertices
searching for microscopic details. Importantly, the masses of the loop
fermions ($\Psi$) generating the GCS couplings, as shown in Fig. 4, must be
invariant both under the SM and the $U(1)^{\prime}$ gauge symmetries to ensure
that the induced low-energy GCS operators are gauge invariant. One can, in
fact, write the microscopic (gauge-invariant) Lagrangian introducing pairs of
heavy fermions ($\Psi$) that are vectorlike with respect to the SM group, but
necessarily chiral under $U(1)^{\prime}$. This generates the effective
Lagrangian (5) at energies below the $U(1)^{\prime}$ breaking scale,

$\begin{split}{\cal L}=&{\cal
L}_{\text{SM}}+\frac{1}{2}(\partial_{\mu}a-M_{Z^{\prime}}Z^{\prime}_{\mu})^{2}-M_{i}~{}\overline{\Psi}^{i}_{L}e^{i(q_{L}-q_{R})\frac{a}{V}}\Psi^{i}_{R}\\\
&+i\overline{\Psi}^{i}_{L}\gamma^{\mu}(\partial_{\mu}-i\frac{\tilde{g}}{2}q_{L}^{i}Z^{\prime}_{\mu})\Psi^{i}_{L}+i\overline{\Psi}^{i}_{R}\gamma^{\mu}(\partial_{\mu}-i\frac{\tilde{g}}{2}q_{R}^{i}Z^{\prime}_{\mu})\Psi^{i}_{R}\end{split}$
(19)

which is manifestly invariant under the (nonlinear) $U(1)^{\prime}$
transformation of parameter $\alpha$,

$\displaystyle\Psi^{i}_{R}\rightarrow\Psi^{i}_{R}e^{i\frac{\tilde{g}}{2}q_{R}\alpha}~{};~{}~{}~{}~{}\Psi^{i}_{L}\rightarrow\Psi^{i}_{L}e^{i\frac{\tilde{g}}{2}q_{L}\alpha}$
$\displaystyle Z^{\prime}_{\mu}\rightarrow
Z^{\prime}_{\mu}+\partial_{\mu}\alpha~{};~{}~{}~{}~{}a\rightarrow
a+\frac{\tilde{g}}{2}V~{}\alpha\equiv a+M_{Z^{\prime}}~{}\alpha$

Figure 4: Triangle diagram containing heavy chiral fermions $\Psi_{i}$ (left
panel) and the resulting effective vertex at low energy (right panel).

From the Lagrangian in Eq. (19), we compute the triangle loops shown in Fig. 4
and integrate out the heavy fermions. We then obtain the same effective
Lagrangian as in Eq. (5), but now we can express the effective coupling of the
dimension-6 Lagrangian in terms of the parameters of the microscopical theory.
In agreement with Ref. Dudas:2013sia , we obtain

${\cal
L}_{\text{loop}}=\frac{1}{\Lambda_{\text{loop}}^{2}}\partial^{\alpha}Z^{\prime}_{\alpha}\epsilon^{\mu\nu\rho\sigma}\mathrm{Tr}[G^{a}_{\mu\nu}G^{a}_{\rho\sigma}],$
(20)

with

$\frac{1}{\Lambda_{\text{loop}}^{2}}=\frac{g_{3}^{2}~{}\tilde{g}}{96\pi^{2}}\sum_{i}\frac{q_{L}^{i}-q_{R}^{i}}{M_{\Psi_{i}}^{2}}\mathrm{Tr}[T^{a}T^{a}].$
(21)

Defining for simplicity
$\sum_{i}\frac{q_{L}^{i}-q_{R}^{i}}{M_{\Psi_{i}}^{2}}\mathrm{Tr}[T^{a}T^{a}]=\frac{N_{\Psi}Q_{\Psi}}{M_{\Psi}^{2}}$
(which corresponds to a set of $N_{\Psi}$ fermions of effective charges
$Q_{\Psi}$ and masses $M_{\Psi}$) we obtain
$\Lambda_{\text{loop}}\simeq\frac{50}{\sqrt{N_{\Psi}Q_{\Psi}}}\frac{M_{\Psi}}{\sqrt{\tilde{g}}}$888We
take the SM expected value of $g_{3}$ at $10^{12}~{}\text{GeV}$.. We can now
reexpress the production rates in Eqs. (14) in terms of the fundamental
parameters of the microscopic theory. For the fermionic dark matter case, we
then have

$R(T)\simeq 5\times
10^{-4}\left(\frac{\alpha~{}N_{\Psi}Q_{\Psi}}{y_{\Psi}^{2}\tilde{g}}\right)^{2}\frac{m_{\chi}^{2}}{V^{8}}T^{10}$
(22)

where we defined $M_{\Psi}=y_{\Psi}V$ and
$M_{Z^{\prime}}=\frac{\tilde{g}}{2}V$. Solving Eq. (16) gives (with all mass
dimensional parameters in GeV units)

$\dfrac{\Omega
h^{2}}{0.12}\simeq\left(\dfrac{B_{F}^{10}}{2}\right)\left(\frac{\alpha
N_{\Psi}Q_{\Psi}}{\tilde{g}y_{\Psi}^{2}}\right)^{2}\left(\frac{m_{\chi}}{10^{10}}\right)^{3}\left(\frac{T_{\text{RH}}}{10^{12}}\right)^{5}\left(\frac{10^{14}}{V}\right)^{8}.$
(23)

We could keep $V$ as a free fundamental parameter of the model, which is
determined by the potential of the Higgs responsible for the extra
$U(1)^{\prime}$ breaking. However, to be more complete, we investigated UV
scenarios in which $V$ is determined as an intermediate scale by the
unification condition of the gauge coupling constants, in $SO(10)$ GUT
constructions (as an example).

Indeed, in such set-ups, the $SO(10)$ group is not directly broken into the SM
in one step but goes through an intermediate gauge group $G_{\text{int}}$ like
$SO(10)\rightarrow G_{\text{int}}\rightarrow SU(3)_{c}\times SU(2)_{L}\times
U(1)_{Y}$. The scale $M_{\text{int}}$ at which the intermediate gauge group is
broken is fixed by the unification condition $g_{1}=g_{2}=g_{3}$ at a higher
unified scale. It was shown in Refs. Mambrini:2013iaa  (at one loop) and
2loops  (at two loops) that $V=M_{\rm int}$ can range from $10^{9}$ to
$10^{15}$ GeV depending on $G_{\text{int}}$ and the representation in which
the Higgs field responsible for $G_{\text{int}}$ breaking lies. We show in
Fig. 5 the parameter space providing the correct relic density for our
fermionic dark matter candidate in several intermediate scenarios. Here, we
take $V=M_{\Psi}=M_{\text{int}}$, which is a reasonable approximation. The
numerical results, obtained by solving the complete set of Boltzmann
equations, Eq. (18), and numerical integration of the rate, Eq. (10), are in
perfect agreement with our analytical solution Eq. (23).

Figure 5: Relic density curves for the fermionic DM case in several $SO(10)$
breaking schemes and for different representations [$\boldmath{r}$] for the
Higgs field responsible for intermediate scale breaking. We take
$N_{\Psi}=Q_{\Psi}=y_{\Psi}=1$ for illustrative purposes.

We observe that in these unified scenarios DM density corresponding to the
Planck measurements planck  can be directly produced from annihilation of SM
particles even if the mediator $Z^{\prime}$ is extremely heavy with no SM
particles charged under the extra $U(1)^{\prime}$. The effective couplings of
the thermal bath to $Z^{\prime}$ are being generated through the GCS
interactions.

Figure 5 highlights the natural relation between the parameters of the theory
$M_{\text{int}}\sim T_{\text{\tiny{MAX}}}\sim
10^{2}~{}T_{\text{\tiny{RH}}}\sim 10^{2}~{}m_{\text{\tiny{DM}}}$ corresponding
to the correct DM relic abundance. Therefore, the intermediate scale in such
unified constructions could be closely related to the inflaton mass as one
expects it to be of the order of the maximum temperature reached by the SM
thermal bath. The large hierarchy between these scales and the SM electroweak
VEV naturally provides the suppressed DM-SM effective coupling required to
produce the correct DM density nonthermally via the freeze-in mechanism.

Though the heavy colored fermions $\Psi$ are cosmologically stable, their very
small number density owing to their heaviness ($\sim 10^{14}$ GeV) keeps them
hidden as benign999For a discussion of the possibility of the DM as a colored
composite object of mass $\sim 10$ TeV reproducing the relic abundance, see a
recent analysis DeLuca:2018mzn , and also an older one Smith:1982qu .. In
principle, these states contribute to the running of $g_{3}$, but they are
kicked into life so late that the nonstandard modification of running during
the remaining phase up to the GUT scale is inconsequential.

If, however, the $\Psi$ states additionally have SM hypercharge, there is a
nontrivial twist. Through loop contributions, they would induce a kinetic
mixing term like $\delta B^{\mu\nu}Z^{\prime}_{\mu\nu}$, which would ascribe
the kinetically diagonal $Z^{\prime}$ with a small, proportional to $\delta$,
coupling with the SM fermion $f$. This would allow the direct s-channel
production of DM from $f$ annihilation as follows: $\bar{f}+f\rightarrow
Z^{\prime}\rightarrow\text{DM}+\text{DM}$. Note that the natural size of
$\delta$ is a loop factor times a logarithm of the ratio of two scales, so
$\delta\sim 10^{-4}$ is a representative number. If $M_{Z^{\prime}}<T_{\rm
MAX}$, then through $Z^{\prime}$, the DM would obviously be in thermal
equilibrium with the SM. But this does not constitute the FIMP or freeze-in
scenario we are pursuing here. On the other hand, if $M_{Z^{\prime}}>T_{\rm
MAX}$, the above-mentioned tree-level $s$-channel annihilation, as shown in
Mambrini:2013iaa , would be the dominating process. The GCS coupling-induced
process we advocated in this paper would then be a subleading one. The
importance of our analysis lies in the fact that even if the heavy $\Psi$
states do not carry the SM hypercharge and the kinetic mixing is absent, the
novel mechanism triggered by the GCS interaction can still explain the relic
DM abundance albeit for an inaccessibly high-mass range of both the DM and its
portals.

##  V V. Conclusion

We have shown that a dark (very) massive $Z^{\prime}$, not charged under the
SM gauge group, can successfully play the role of a mediator between the
visible and the dark sectors even if the corresponding $U(1)^{\prime}$
breaking scale lies far above the maximum temperature of the Universe. Pair
annihilation of SM gauge bosons, proceeding through triangle loops containing
heavy fermions through this $Z^{\prime}$ portal, can produce a cosmologically
agreeable amount of DM. These types of effective couplings between the SM
gauge bosons and $Z^{\prime}$ find inherent justification in an anomaly-free
set-up in which Chern-Simons (more precisely, the GCS type discussed in the
beginning) terms are generated through the anomaly cancellation mechanism.

The large effective scale of the GCS operators is responsible for the weakness
of the DM-SM interaction strength without invoking unnaturally small couplings
as is often required in the context of freeze-in. The large dependence of the
production rate on temperature indicates that the majority of DM is produced
during the initial moments of reheating. Subsequently, the reheating process
itself lends important consequences to the computation of DM production.
Moreover, the assertion that the $U(1)^{\prime}$ breaking scale lies far above
the reheating temperature implies that only the longitudinal mode (axion à la
the Stueckelberg formalism) of the $Z^{\prime}$ contributes to the production
process, rendering its phenomenology very particular in comparison with other
type of mediators. The prominence of the longitudinal mode is also consistent
with how $Z^{\prime}$ is coupled via GCS interaction. Such a scenario can be
embedded in a unified $SO(10)$ framework in which the $Z^{\prime}$ mass scale
represents an intermediate breaking stage of $SO(10)$.

Acknowledgments: The authors want to thank especially E. Dudas for very
insightful discussions. This research has been supported by the (Indo-French)
CEFIPRA/IFCPAR Project No. 5404-2. Support from CNRS LIA-THEP and the INFRE-
HEPNET of CEFIPRA/IFCPAR is also acknowledged. G.B. acknowledges support of
the J.C. Bose National Fellowship from the Department of Science and
Technology, Government of India (SERB Grant No. SB/S2/JCB-062/2016). M.D.
acknowledges support from the Brazilian Ph.D. program “Ciências sem
Fronteiras” – CNPQ Process No. 202055/2015-9. This work was also supported by
the France-US PICS no. 06482, PICS MicroDark. M. D. and Y. M. acknowledge
partial support from the European Union Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie: RISE InvisiblesPlus (Grant No
690575) and the ITN Elusives (Grant No 674896).

## Appendix

##  VI noninstantaneous reheating: analytical estimations

As the reheating is not an instantaneous process, the inflaton dominates the
energy density of the Universe at some stage and we have a different relation
between time and temperature Giudice:2000ex  101010The numerical factors $c$
and $d$ are related to the convention chosen in the definition of reheating
temperature: $t_{\phi}=\frac{c}{\Gamma_{\phi}}$ is the time of inflaton decay
completion and $t_{H}=\frac{d}{H}$ is the Hubble time. The reheating
temperature is defined such that
$\Gamma_{\phi}=\frac{c}{d}H(T_{\text{\tiny{RH}}})$.

$\frac{\text{d}}{\text{d}t}=-\frac{3}{8}H(T)T\frac{\text{d}}{\text{d}T},~{}~{}\textrm{with}~{}~{}H(T)=\frac{d}{c}\sqrt{\frac{5g_{e}^{2}}{72g_{\text{\tiny{RH}}}}}\pi\frac{T^{4}}{T_{\text{\tiny{RH}}}^{2}M_{P}},$
(24)

with $g_{\text{\tiny{RH}}}=g_{e}(T_{\text{\tiny{RH}}})$.

We emphasize that while the inflaton dominates, we have $T\propto a^{-3/8}$
instead of $a^{-1}$ as in the radiation-dominated era, with $a$ the expansion
scale factor.

It is convenient to define the dimensionless quantity
$\mathcal{N}_{\phi}\equiv\rho_{\phi}/T_{\text{\tiny{RH}}}a^{3}$ Giudice:2000ex
, with $\rho_{\phi}$ the energy density of the inflaton field. It remains
constant before the end of reheating, with a value given by 111111We have
$\mathcal{N}_{\phi}^{I}\equiv
b_{2}^{-8}\frac{g_{\text{\tiny{MAX}}}^{2}}{g_{\text{\tiny{RH}}}}(\frac{T_{\text{\tiny{MAX}}}}{T_{\text{\tiny{RH}}}})^{8}$,
with
$b_{2}=\Big{(}\frac{3^{11/10}5^{1/2}}{2^{23/10}\pi}\frac{c}{d}\Big{)}^{1/4}$
and $g_{\text{\tiny{MAX}}}=g_{e}(T_{\text{\tiny{MAX}}})$.
$\mathcal{N}_{\phi}^{I}\propto(\frac{T_{\text{\tiny{MAX}}}}{T_{\text{\tiny{RH}}}})^{8}$.

We notice that, by using Eq. (24), Eq. (9) can be put in the form

$\frac{\text{d}Y_{\text{\tiny{ID}}}}{\text{d}T}=-\frac{8}{3}\frac{R(T)}{H~{}T~{}\epsilon},\hskip
14.22636pt\textrm{for}~{}~{}T_{\text{\tiny{RH}}}<T<T_{\text{\tiny{MAX}}}$ (25)

where the yield in the inflaton-dominated era is defined as the dimensionless
and comoving parameter $Y_{\text{\tiny{ID}}}\equiv n/\epsilon$ and $\epsilon$
is defined as

$\epsilon\equiv\mathcal{N}_{\phi}^{I}a^{-3}=\frac{5d^{2}g_{e}^{2}}{96c^{2}g_{\text{\tiny{RH}}}T_{\text{\tiny{RH}}}^{5}}T^{8}.$
(26)

By solving Eq. (25), we can find the contribution of the noninstantaneous
heating process to the relic density and define the boost factor discussed in
Sec. III.

##  VII Decay width

We present the decay width used in the computation of the production rate.
Allowing the $Z^{\prime}$ to decay into $N_{\psi}$ dark heavy fermions $\psi$
(with only vectorial coupling) of charge $Q_{\psi}$ and mass $m_{\psi}$, the
total decay width is given by

$\Gamma_{Z^{\prime}}=\frac{M_{\text{$Z^{\prime}$}}}{12\pi}N_{\psi}Q_{\psi}^{2}\sqrt{1-\frac{4m_{\psi}^{2}}{M_{\text{$Z^{\prime}$}}^{2}}}\Big{(}1+\frac{2m_{\psi}^{2}}{M_{\text{$Z^{\prime}$}}^{2}}\Big{)}+\Gamma_{\text{$Z^{\prime}$}\rightarrow\text{DM}},$
(27)

where

$\Gamma_{\text{$Z^{\prime}$}\rightarrow\text{DM}}=\frac{M_{\text{$Z^{\prime}$}}}{12\pi}\times\left\\{\begin{array}[]{rc}\alpha^{2}\Big{(}1-\frac{4m_{\chi}^{2}}{M_{\text{$Z^{\prime}$}}^{2}}\Big{)}^{3/2}~{}(\text{fermionic
DM})\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\frac{\beta^{2}}{2}\frac{M_{\text{$Z^{\prime}$}}^{2}}{m_{\text{\tiny{$X_{1}$}}}^{2}}\Big{(}1-\frac{4m_{\text{\tiny{$X_{1}$}}}^{2}}{M_{\text{$Z^{\prime}$}}^{2}}\Big{)}^{5/2}~{}~{}~{}(\text{Abelian
DM})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0~{}~{}~{}~{}~{}~{}~{}~{}~{}(\text{non-Abelian DM})\end{array}\right.$ (28)

The decay width of $Z^{\prime}$ into non-Abelian DM vanishes identically
because of the nature of the coupling.

## References

  * (1) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 594, A13 (2016) [arXiv:1502.01589 [astro-ph.CO]].
  * (2) E. Aprile et al. [XENON100 Collaboration], Phys. Rev. Lett. 109 (2012) 181301 [arXiv:1207.5988 [astro-ph.CO]].
  * (3) D. S. Akerib et al. [LUX Collaboration], Phys. Rev. Lett. 118 (2017) no.2, 021303 [arXiv:1608.07648 [astro-ph.CO]].
  * (4) C. Fu et al. [PandaX-II Collaboration], Phys. Rev. Lett. 118, no. 7, 071301 (2017) [arXiv:1611.06553 [hep-ex]].
  * (5) J. A. Casas, D. G. Cerdeño, J. M. Moreno and J. Quilis, JHEP 1705 (2017) 036 doi:10.1007/JHEP05(2017)036 [arXiv:1701.08134 [hep-ph]]. A. Djouadi, O. Lebedev, Y. Mambrini and J. Quevillon, Phys. Lett. B 709 (2012) 65 doi:10.1016/j.physletb.2012.01.062 [arXiv:1112.3299 [hep-ph]]; A. Djouadi, A. Falkowski, Y. Mambrini and J. Quevillon, Eur. Phys. J. C 73 (2013) no.6, 2455 doi:10.1140/epjc/s10052-013-2455-1 [arXiv:1205.3169 [hep-ph]]; O. Lebedev, H. M. Lee and Y. Mambrini, Phys. Lett. B 707 (2012) 570 doi:10.1016/j.physletb.2012.01.029 [arXiv:1111.4482 [hep-ph]]; Y. Mambrini, Phys. Rev. D 84 (2011) 115017 doi:10.1103/PhysRevD.84.115017 [arXiv:1108.0671 [hep-ph]].
  * (6) J. Ellis, A. Fowlie, L. Marzola and M. Raidal, arXiv:1711.09912 [hep-ph]. G. Arcadi, Y. Mambrini and F. Richard, JCAP 1503 (2015) 018 [arXiv:1411.2985 [hep-ph]]; J. Kearney, N. Orlofsky and A. Pierce, arXiv:1611.05048 [hep-ph]; M. Escudero, A. Berlin, D. Hooper and M. X. Lin, JCAP 1612 (2016) 029 doi:10.1088/1475-7516/2016/12/029 [arXiv:1609.09079 [hep-ph]].
  * (7) A. Alves, S. Profumo and F. S. Queiroz, JHEP 1404 (2014) 063 doi:10.1007/JHEP04(2014)063 [arXiv:1312.5281 [hep-ph]]. O. Lebedev and Y. Mambrini, Phys. Lett. B 734 (2014) 350 [arXiv:1403.4837 [hep-ph]]; G. Arcadi, Y. Mambrini, M. H. G. Tytgat and B. Zaldivar, JHEP 1403 (2014) 134 [arXiv:1401.0221 [hep-ph]]; O. Lebedev and Y. Mambrini, Phys. Lett. B 734 (2014) 350 doi:10.1016/j.physletb.2014.05.025 [arXiv:1403.4837 [hep-ph]].
  * (8) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo and F. S. Queiroz, Eur. Phys. J. C 78 (2018) no.3, 203 doi:10.1140/epjc/s10052-018-5662-y [arXiv:1703.07364 [hep-ph]].
  * (9) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 1003 (2010) 080 [arXiv:0911.1120 [hep-ph]]; X. Chu, T. Hambye and M. H. G. Tytgat, JCAP 1205 (2012) 034 [arXiv:1112.0493 [hep-ph]]; X. Chu, Y. Mambrini, J. Quevillon and B. Zaldivar, JCAP 1401 (2014) 034 [arXiv:1306.4677 [hep-ph]].
  * (10) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, Int. J. Mod. Phys. A 32 (2017) no.27, 1730023 doi:10.1142/S0217751X1730023X [arXiv:1706.07442 [hep-ph]].
  * (11) Z. G. Berezhiani, A. D. Dolgov and R. N. Mohapatra, Phys. Lett. B 375 (1996) 26 doi:10.1016/0370-2693(96)00219-5 [hep-ph/9511221]; P. Adshead, Y. Cui and J. Shelton, JHEP 1606 (2016) 016 doi:10.1007/JHEP06(2016)016 [arXiv:1604.02458 [hep-ph]]; E. Hardy and J. Unwin, JHEP 1709 (2017) 113 doi:10.1007/JHEP09(2017)113 [arXiv:1703.07642 [hep-ph]].
  * (12) N. Bernal, M. Dutra, Y. Mambrini, K. A. Olive, M. Peloso and M. Pierre, arXiv:1803.01866 [hep-ph].
  * (13) M. Garny, A. Palessandro, M. Sandora and M. S. Sloth, JCAP 1802, no. 02, 027 (2018) doi:10.1088/1475-7516/2018/02/027 [arXiv:1709.09688 [hep-ph]].
  * (14) M. Garny, M. Sandora and M. S. Sloth, Phys. Rev. Lett. 116 (2016) no.10, 101302 doi:10.1103/PhysRevLett.116.101302 [arXiv:1511.03278 [hep-ph]].
  * (15) A. Biswas, D. Borah and A. Dasgupta, arXiv:1805.06903 [hep-ph].
  * (16) Y. Mambrini, K. A. Olive, J. Quevillon and B. Zaldivar, Phys. Rev. Lett. 110 (2013) no.24, 241306 [arXiv:1302.4438 [hep-ph]];
  * (17) Y. Mambrini, N. Nagata, K. A. Olive, J. Quevillon and J. Zheng, Phys. Rev. D 91 (2015) no.9, 095010 doi:10.1103/PhysRevD.91.095010 [arXiv:1502.06929 [hep-ph]]; Y. Mambrini, N. Nagata, K. A. Olive and J. Zheng, Phys. Rev. D 93 (2016) no.11, 111703 doi:10.1103/PhysRevD.93.111703 [arXiv:1602.05583 [hep-ph]].
  * (18) P. Anastasopoulos, M. Bianchi, E. Dudas and E. Kiritsis, JHEP 0611 (2006) 057 doi:10.1088/1126-6708/2006/11/057 [hep-th/0605225].
  * (19) J. A. Dror, R. Lasenby and M. Pospelov, Phys. Rev. D 96 (2017) no.7, 075036 doi:10.1103/PhysRevD.96.075036 [arXiv:1707.01503 [hep-ph]]; F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz and S. Vogl, JHEP 1602 (2016) 016 doi:10.1007/JHEP02(2016)016 [arXiv:1510.02110 [hep-ph]]; Y. Cui and F. D’Eramo, Phys. Rev. D 96 (2017) no.9, 095006 doi:10.1103/PhysRevD.96.095006 [arXiv:1705.03897 [hep-ph]]; J. A. Dror, R. Lasenby and M. Pospelov, Phys. Rev. Lett. 119 (2017) no.14, 141803 doi:10.1103/PhysRevLett.119.141803 [arXiv:1705.06726 [hep-ph]].
  * (20) A. Ismail, A. Katz and D. Racco, JHEP 1710 (2017) 165 doi:10.1007/JHEP10(2017)165 [arXiv:1707.00709 [hep-ph]]; G. Arcadi, P. Ghosh, Y. Mambrini, M. Pierre and F. S. Queiroz, JCAP 1711 (2017) no.11, 020 doi:10.1088/1475-7516/2017/11/020 [arXiv:1706.04198 [hep-ph]]; S. M. Choi, Y. Hochberg, E. Kuflik, H. M. Lee, Y. Mambrini, H. Murayama and M. Pierre, JHEP 1710 (2017) 162 doi:10.1007/JHEP10(2017)162 [arXiv:1707.01434 [hep-ph]].
  * (21) Y. Farzan and A. R. Akbarieh, JCAP 1411 (2014) no.11, 015 doi:10.1088/1475-7516/2014/11/015 [arXiv:1408.2950 [hep-ph]].
  * (22) G. F. Giudice, E. W. Kolb and A. Riotto, Phys. Rev. D 64 (2001) 023508 doi:10.1103/PhysRevD.64.023508 [hep-ph/0005123].
  * (23) M. A. G. Garcia, Y. Mambrini, K. A. Olive and M. Peloso, Phys. Rev. D 96 (2017) no.10, 103510 doi:10.1103/PhysRevD.96.103510 [arXiv:1709.01549 [hep-ph]].
  * (24) K. Benakli, Y. Chen, E. Dudas and Y. Mambrini, Phys. Rev. D 95 (2017) no.9, 095002 doi:10.1103/PhysRevD.95.095002 [arXiv:1701.06574 [hep-ph]]; E. Dudas, Y. Mambrini and K. Olive, Phys. Rev. Lett. 119 (2017) no.5, 051801 doi:10.1103/PhysRevLett.119.051801 [arXiv:1704.03008 [hep-ph]]; E. Dudas, T. Gherghetta, Y. Mambrini and K. A. Olive, Phys. Rev. D 96 (2017) no.11, 115032 doi:10.1103/PhysRevD.96.115032 [arXiv:1710.07341 [hep-ph]]; E. Dudas, T. Gherghetta, K. Kaneta, Y. Mambrini and K. A. Olive, arXiv:1805.07342 [hep-ph].
  * (25) E. Dudas, L. Heurtier, Y. Mambrini and B. Zaldivar, JHEP 1311 (2013) 083 doi:10.1007/JHEP11(2013)083 [arXiv:1307.0005 [hep-ph]].
  * (26) E. Dudas, Y. Mambrini, S. Pokorski and A. Romagnoni, JHEP 1210 (2012) 123 doi:10.1007/JHEP10(2012)123 [arXiv:1205.1520 [hep-ph]]; Y. Mambrini, JCAP 0912 (2009) 005 doi:10.1088/1475-7516/2009/12/005 [arXiv:0907.2918 [hep-ph]]; E. Dudas, Y. Mambrini, S. Pokorski and A. Romagnoni, JHEP 0908 (2009) 014 doi:10.1088/1126-6708/2009/08/014 [arXiv:0904.1745 [hep-ph]].
  * (27) G. Bélanger, J. Da Silva and H. M. Tran, Phys. Rev. D 95 (2017) no.11, 115017 doi:10.1103/PhysRevD.95.115017 [arXiv:1703.03275 [hep-ph]]. J. L. Feng, J. Smolinsky and P. Tanedo, Phys. Rev. D 93 (2016) no.11, 115036 Erratum: [Phys. Rev. D 96 (2017) no.9, 099903] doi:10.1103/PhysRevD.93.115036, 10.1103/PhysRevD.96.099903 [arXiv:1602.01465 [hep-ph]]. Y. Mambrini, JCAP 1107 (2011) 009 doi:10.1088/1475-7516/2011/07/009 [arXiv:1104.4799 [hep-ph]].
  * (28) M. Blennow, E. Fernandez-Martinez and B. Zaldivar, JCAP 1401 (2014) 003 doi:10.1088/1475-7516/2014/01/003 [arXiv:1309.7348 [hep-ph]].
  * (29) D. J. H. Chung, E. W. Kolb and A. Riotto, Phys. Rev. D 60 (1999) 063504 doi:10.1103/PhysRevD.60.063504 [hep-ph/9809453].
  * (30) V. De Luca, A. Mitridate, M. Redi, J. Smirnov and A. Strumia, Phys. Rev. D 97 (2018) no.11, 115024 doi:10.1103/PhysRevD.97.115024 [arXiv:1801.01135 [hep-ph]].
  * (31) P. F. Smith, J. R. J. Bennett, G. J. Homer, J. D. Lewin, H. E. Walford and W. A. Smith, Nucl. Phys. B 206 (1982) 333. doi:10.1016/0550-3213(82)90271-1.