Gravitational Baryogenesis: Problems and Possible Resolution

Abstract

The coupling of baryonic current to the derivative of the curvature scalar, R, inherent to gravitational baryogenesis (GBG), leads to a fourth-order differential equation of motion for R instead of the algebraic one of general relativity (GR). The fourth-order differential equation is generically unstable. We consider a possible mechanism of stabilization of GBG by the modification of gravity, introducing an $R^2$ term into the canonical action of GR. It is shown that this mechanism allows for the stabilization of GBG with bosonic and fermionic baryon currents. We establish the region of the model parameters leading to the stabilization of R. Still, the standard cosmology would be noticeably modified.

1. Introduction

An excess of matter over antimatter in our Universe is crucial for our very existence and is well-supported by various observation. The local Universe is clearly matter-dominated. The amount of antimatter is very small, and it can be explained as the result of high-energy collisions in space. On the other hand, matter and antimatter seem to have similar properties; therefore, we could expect a matter–antimatter symmetric universe. The existence of large regions of antimatter in our neighbourhood would produce high-energy radiation created by matter–antimatter annihilation on the boundaries between matter and antimatter domains, which is not observed. A satisfactory model of our Universe should be able to explain the origin of the matter–antimatter asymmetry. Any initial asymmetry at inflation could not solve the problem of observed excess of matter over antimatter, because the energy density associated with the observed nonzero baryonic number density would not allow for a sufficiently long inflation.

The term baryogenesis is used to indicate the generation of the excess of matter (baryons) over antimatter (antibaryons) or vice versa.

In 1967, Andrey Sakharov formulated three conditions, today known as Sakharov’s principles [1], necessary to produce a matter–antimatter asymmetry in the initially symmetric universe. These conditions include:

1.

Nonconservation of baryonic number;

2.

Breaking of the symmetry between particles and antiparticles;

3.

Deviation from thermal equilibrium.

However, not all of the three Sakharov principles are strictly necessary. For example, spontaneous baryogenesis (SBG) and gravitational baryogenesis (GBG) do not demand an explicit C and CP violation and can proceed in thermal equilibrium. Moreover, these mechanisms are usually most efficient in thermal equilibrium.

The statement that the cosmological baryon asymmetry can be created by spontaneous baryogenesis in thermal equilibrium was mentioned in the original paper by A. Cohen and D. Kaplan in 1987 [2] and in the subsequent papers by A. Cohen, D. Kaplan, and A. Nelson [3,4] (for a review, see [5,6,7,8]).

The term “spontaneous” is related to a spontaneous breaking of the underlying symmetry of the theory, which ensures the conservation of the total baryonic number in the unbroken phase. This symmetry is supposed to be spontaneously broken, and in the broken phase, the Lagrangian density acquires the term

$$\begin{array}{c}\hfill {\mathcal{L}}_{SBG}=\left({\partial }_{\mu }\theta \right)J_B^{\mu },\end{array}$$

(1)

where $\theta$ is a (pseudo)Goldstone field, and $J_B^{\mu }$ is the baryonic current of matter fields, which becomes nonconserved as a result of the symmetry breaking.

For a spatially homogeneous field, $\theta =\theta \left(t\right)$, the Lagrangian is reduced to a simple form

$$\begin{array}{c}\hfill {\mathcal{L}}_{SBG}=\dot{\theta }{n}_B,n_B\equiv J_B^0.\end{array}$$

(2)

Here, $n_B$ is the baryonic number density, so it is tempting to identify $\dot{\theta }$ with the chemical potential, ${\mu }_B$, of the corresponding system. However, such an identification is questionable [9,10]. It depends upon the representation chosen for the fermionic fields and is heavily based on the assumption $\dot{\theta }\approx const$. In Ref. [9], the assumption $\dot{\theta }\approx const$ was relaxed.

Stimulated by spontaneous baryogenesis, the idea of gravitational baryogenesis was put forward [11]. The scenario of SBG was modified by the introduction of the coupling of the baryonic current to the derivative of the curvature scalar R:

$$\begin{array}{c}\hfill {𝒮}_{GBG}=-\frac{1}{M^2}\int d^{4}x\sqrt{-g}\left({\partial }_{\mu }R\right)J_B^{\mu },\end{array}$$

(3)

where g is the determinant of the space–time metric tensor and the mass parameter M determines the energy scale of baryogenesis. There are a lot of articles on the subject, and a partial list of references is included in Refs. [12,13,14,15,16]. According to these papers, the GBG mechanism can successfully explain the magnitude of the cosmological baryon asymmetry of the universe.

However, it was argued in Refs. [17,18] that the backreaction of the created nonzero baryonic density on the space–time curvature led to a strong instability of the cosmological evolution. In this paper, we show that the problem of stability can be solved by adding to the Hilbert–Einstein action the quadratic in curvature term.The underlying gravitational action takes the form:

$$\begin{array}{c}\hfill S_{Grav}=-\frac{M_{Pl}^2}{16\pi }\int d^{4}x\sqrt{-g}\left(R-\frac{R^2}{6M_R^2}\right),\end{array}$$

(4)

where $M_{Pl}=1.22·{10}^{19}$ GeV is the Planck mass, and we use the metric signature $(+,-,-,-)$.

The $R^2$ term (among many others) in the canonical action of GR appears as a result of one-loop corrections to the energy–momentum tensor of matter in curved spacetime as was first found in [19]. Subsequently, this approach was developed by V. T. Gurovich and A. A. Starobinsky [20]. Afterwards, it was noticed by A. A. Starobinsky [21] that the $R^2$ term would lead to a successful cosmological inflationary model, albeit with a large coefficient ${(M_{Pl}/M_R)}^2$. This is now a very popular inflationary scenario, since it predicts a small amplitude of the relic gravitational waves at low frequencies in agreement with CMB restrictions.

As is known, the $R^2$ term leads to the excitation of the scalar degree of freedom, named scalaron, and $M_R$ is the scalaron mass. In the very early universe, the $R^2$ term could generate inflation [21] and density perturbations. The amplitude of the observed density perturbations demands that $M_R=3·{10}^{13}$ GeV [22] if the scalaron is the inflaton. Otherwise, $M_R>3·{10}^{13}$ GeV is allowed. Below, we presume that the scalaron is the inflaton.

2. Instability Problem of Gravitational Baryogenesis

The essential ingredient of gravitational baryogenesis is the coupling of the baryonic current to the derivative of the curvature scalar ${\partial }_{\mu }R$ (3). Taken over the canonical cosmological Friedmann–Lemaitre–Robertson–Walker background, this interaction can successfully fulfil the task of generating the proper value of the baryon asymmetry of the universe.

However, any curvature-dependent term in the Lagrangian of the theory would modify the equations of general relativity (GR). The modified GR equations have been analysed in Refs. [17,18]. Since interaction (3) is not just linear in the curvature term multiplied by a constant, it leads to higher-order equations describing the evolution of gravitational fields. Higher-order equations of motion are typically unstable with respect to small perturbations. According to the results of Refs. [17,18], it indeed happens in the frameworks of the GBG scenario and the characteristic time of the exponential instability is much shorter than the cosmological time. It creates a serious problem for the realisation of the GBG mechanism.

In this work, we suggest to consider the possible stabilisation of GBG and prove that it can be realised, although the resulting cosmological model suffers from too large a value of R, much larger than that in the classical Friedmann cosmology. Possible ways to cure this shortcoming are mentioned.

3. Stabilisation of Gravitational Baryogenesis in Modified Gravity

3.1. Bosonic Case

Let us first consider the case when a baryonic number is carried by a complex scalar field $\varphi$ [17]. The total action has the form:

$$\begin{array}{c}\hfill S_{tot}\left[\varphi \right]=-\int d^{4}x\sqrt{-g}\left[\frac{M_{Pl}^2}{16\pi }\left(R-\frac{R^2}{6M_R^2}\right)+\frac{1}{M^2}\left({\partial }_{\mu }R\right)J_{\left(\varphi \right)}^{\mu }-g^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }{\varphi }^{*}+U(\varphi ,{\varphi }^{*})\right]\\ \hfill +S_{matt},\end{array}$$

(5)

where $U(\varphi ,{\varphi }^{*})$ is the potential of field $\varphi$ and $S_{matt}$ is the matter action which does not include the field $\varphi$. In Equation (5), $R\left(t\right)$ is the classical curvature field, while $\varphi (\overrightarrow{x},t)$ is the quantum operator of light scalar particles.

We assume that the potential $U(\varphi ,{\varphi }^{*})$ is not invariant with respect to phase transformation $\varphi \to exp\left(iq\beta \right)\varphi$ and thus the corresponding current

$$\begin{array}{c}\hfill J_{\left(\varphi \right)}^{\mu }=iq{g}^{\mu \nu }({\varphi }^{*}{\partial }_{\nu }\varphi -\varphi {\partial }_{\nu }{\varphi }^{*})\end{array}$$

(6)

is not conserved. Here, q is the baryonic number of field $\varphi$. The nonconservation of the current is necessary for the proper performance of the model, otherwise $S_{GBG}$ in Equation (3) can be integrated away by parts.

Varying action (5) over $g^{\mu \nu }$, we come to the following equations:

$$\begin{array}{ccc}\hfill & & \frac{M_{Pl}^2}{16\pi }\left[R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R-\frac{1}{3M_R^2}\left(R_{\mu \nu }-\frac{1}{4}{g}_{\mu \nu }R+g_{\mu \nu }{D}^2-D_{\mu }{D}_{\nu }\right)R\right]\hfill \\ \hfill & & -\frac{1}{M^2}\left[\left(R_{\mu \nu }-(D_{\mu }{D}_{\nu }-g_{\mu \nu }{D}^2)\right)D_{\alpha }{J}_{\left(\varphi \right)}^{\alpha }+\frac{1}{2}{g}_{\mu \nu }{J}_{\left(\varphi \right)}^{\alpha }{D}_{\alpha }R-\frac{1}{2}\left(J_{\left(\varphi \right)\nu }{D}_{\mu }R+J_{\left(\varphi \right)\mu }{D}_{\nu }R\right)\right]\hfill \\ \hfill & & -\frac{1}{2}\left(D_{\mu }\varphi D_{\nu }{\varphi }^{*}+D_{\nu }\varphi D_{\mu }{\varphi }^{*}\right)+\frac{1}{2}{g}_{\mu \nu }\left[D_{\alpha }\varphi D^{\alpha }{\varphi }^{*}-U\left(\varphi \right)\right]-\left(D_{\mu }\varphi \right)\left(D_{\nu }{\varphi }^{*}\right)\hfill \\ & & =\frac{1}{2}{T}_{\mu \nu }^{\left(matt\right)},\hfill \end{array}$$

(7)

where $D_{\mu }$ is the covariant derivative in metric $g_{\mu \nu }$ (of course, for scalars $D_{\mu }={\partial }_{\mu }$) and $T_{\mu \nu }^{\left(matt\right)}$ is the energy–momentum tensor of matter obtained from action $S_{matt}$.

Taking the trace of Equation (7) with respect to $\mu$ and $\nu$ and changing the sign, we obtain:

$$\begin{array}{c}\hfill \frac{M_{Pl}^2}{16\pi }\left(R+\frac{1}{M_R^2}{D}^{2}R\right)+\frac{1}{M^2}\left[(R+3D^2)D_{\alpha }{J}_{\left(\varphi \right)}^{\alpha }+J_{\left(\varphi \right)}^{\alpha }{D}_{\alpha }R\right]-D_{\alpha }\varphi D^{\alpha }{\varphi }^{*}+2U\left(\varphi \right)\\ \hfill =-\frac{1}{2}{T}^{\left(matt\right)}=0,\end{array}$$

(8)

where $T^{\left(matt\right)}=g^{\mu \nu }{T}_{\mu \nu }^{\left(matt\right)}$ is the trace of the energy–momentum tensor of matter. For the usual relativistic matter, $T^{\left(matt\right)}=0$, while for scalar field $\varphi$, the trace of the energy–momentum tensor is nonzero:

$$\begin{array}{c}\hfill T_{\mu }^{\mu }\left(\varphi \right)=-2D_{\alpha }\varphi D^{\alpha }{\varphi }^{*}+4U\left(\varphi \right).\end{array}$$

(9)

The equation of motion for field $\varphi$ is:

$$\begin{array}{c}\hfill D^2\varphi +\frac{\partial U}{\partial {\varphi }^{*}}=-\frac{iq}{M^2}\left(2D_{\mu }R{D}^{\mu }\varphi +\varphi D^{2}R\right).\end{array}$$

(10)

According to definition (6) and Equation (10), the current divergence is:

$$\begin{array}{c}\hfill D_{\mu }{J}^{\mu }=\frac{2q^2}{M^2}\left[D_{\mu }R({\varphi }^{*}{D}^{\mu }\varphi +\varphi D^{\mu }{\varphi }^{*})+{\left|\varphi \right|}^2{D}^{2}R\right]+iq\left(\varphi \frac{\partial U}{\partial \varphi }-{\varphi }^{*}\frac{\partial U}{\partial {\varphi }^{*}}\right).\end{array}$$

(11)

For a homogeneous curvature scalar $R\left(t\right)$ in a spatially flat FLRW-metric

$$\begin{array}{c}\hfill d{s}^2=d{t}^2-a^2\left(t\right)d{\mathbf{r}}^2\end{array}$$

(12)

Equation (8) is reduced to:

$$\begin{array}{c}\hfill \frac{M_{Pl}^2}{16\pi }\left[R+\frac{1}{M_R^2}({\partial }_t^2+3H{\partial }_t)R\right]+\frac{1}{M^2}\left[(R+3{\partial }_t^2+9H{\partial }_t)D_{\alpha }{J}_{\left(\varphi \right)}^{\alpha }+\dot{R}{J}_{\left(\varphi \right)}^0\right]\\ \hfill +2U\left(\varphi \right)-\left(D_{\alpha }\varphi \right)\left(D^{\alpha }{\varphi }^{*}\right)=0.\end{array}$$

(13)

where $J_{\left(\varphi \right)}^0$ is the baryonic number density of the $\varphi$-field, $H=\dot{a}/a$ is the Hubble parameter, and the divergence of the current is given by the expression:

$$\begin{array}{c}\hfill D_{\alpha }{J}_{\left(\varphi \right)}^{\alpha }=\frac{2q^2}{M^2}\left[\dot{R}({\varphi }^{*}\dot{\varphi }+\varphi {\dot{\varphi }}^{*})+(\ddot{R}+3H\dot{R}){\varphi }^{*}\varphi \right]+iq\left(\varphi \frac{\partial U}{\partial \varphi }-{\varphi }^{*}\frac{\partial U}{\partial {\varphi }^{*}}\right).\end{array}$$

(14)

As we see in what follows, the last two terms in Equation (13) do not have an essential impact on the cosmological instability found in Ref. [17] and are disregarded below. Indeed, as shown in Ref. [17], the field $\varphi$ does not exponentially rise together with R and thus can be neglected in comparison with the terms containing R. In the case considered here of a modified $R^2$-gravity, the curvature also initially strongly rises before the $R^2$ term starts to operate and, in this sense, the situation is the same as that studied in Ref. [18]. In fact, Equation (15) is a good argument in favour of the subdominant nature of the terms containing $\varphi$ above.

Let us note that the statement of exponential instability of $R\left(t\right)$ [17] does not depend on the conservation or nonconservation of the current from the potential term $\left(\varphi \partial U/\partial \varphi -{\varphi }^{*}\partial U/\partial {\varphi }^{*}\right)$ in Equation (14). However, if the current from this term is conserved, then the baryon asymmetry is not generated. On the other hand, the term in square brackets in Equation (14) does not lead to the generation of the baryon asymmetry but leads to the exponential instability of $R\left(t\right)$. Below, we ignore the last term of Equation (14).

Performing thermal averaging of the normal-ordered bilinear products of field $\varphi$ in the high temperature limit (see Appendix of Ref. [17]) in accordance with equations:

$$\begin{array}{c}\hfill \langle {\varphi }^{*}\varphi \rangle =\frac{T^2}{12},\langle {\varphi }^{*}\dot{\varphi }+{\dot{\varphi }}^{*}\varphi \rangle =0,\end{array}$$

(15)

and using Equation (14), we obtain the fourth-order differential equation:

$$\begin{array}{c}\hfill \frac{M_{Pl}^2}{16\pi }\left(R+\frac{1}{M_R^2}{D}^{2}R\right)+\frac{q^2}{6M^4}\left(R+3{\partial }_t^2+9H{\partial }_t\right)\left[\left(\ddot{R}+3H\dot{R}\right)T^2\right]+\frac{1}{M^2}\dot{R}\langle J_{\left(\varphi \right)}^0\rangle \\ \hfill =-2U\left(\varphi \right)+\left(D_{\alpha }\varphi \right)\left(D^{\alpha }{\varphi }^{*}\right).\end{array}$$

(16)

Here, $\langle J_{\left(\varphi \right)}^0\rangle$ is the thermal average value of the baryonic number density of $\varphi$, which is supposed to vanish initially, but created through the process of gravitational baryogenesis. This term can be neglected because the baryon asymmetry is normally quite small. Even if it is not small, it does not have a considerable impact on the explosive rise of the curvature scalar. As we see in what follows, the evolution of $R\left(t\right)$ proceeds much faster than the cosmological evolution, that is $\ddot{R}/\dot{R}\gg H$. Consequently, we neglect the terms proportional to R with respect to the terms proportional to the second derivative of R, $\ddot{R}$. We also consider the terms of the type $HR$ as small with respect to $dR/dt$. We can check that this presumption is true a posteriori with the obtained solution for $R\left(t\right)$.

Keeping only the dominant terms we simplify the above equation to:

$$\begin{array}{c}\hfill \frac{d^{4}R}{d{t}^4}+\frac{{\kappa }^4}{M_R^2}\frac{d^{2}R}{d{t}^2}+{\kappa }^{4}R=-T_{\mu }^{\mu }\left(\varphi \right)\frac{M^4}{q^2{T}^2},\end{array}$$

(17)

where

$$\begin{array}{c}\hfill {\kappa }^4=\frac{M_{Pl}^2{M}^4}{8\pi q^2{T}^2}.\end{array}$$

(18)

While studying the instability of the solution, we do not take into account the r.h.s. of Equation (17) which does not depend upon R. Looking for the solution of Equation (17) in the form $R=R_{in}exp\left(\lambda t\right)$, we obtain the characteristic equation:

$$\begin{array}{c}\hfill {\lambda }^4+\frac{{\kappa }^4}{M_R^2}{\lambda }^2+{\kappa }^4=0\end{array}$$

(19)

with the eigenvalues $\lambda$ defined by the expression:

$$\begin{array}{c}\hfill {\lambda }^2=-\frac{{\kappa }^4}{2M_R^2}\pm {\kappa }^2\sqrt{\frac{{\kappa }^4}{4M_R^4}-1}.\end{array}$$

(20)

There is no instability if ${\lambda }^2<0$ and Equation (17) has only oscillating solutions. It is realised if ${\kappa }^4>4M_R^4$. Using the expression in Equation (18) for ${\kappa }^4$ and taking $M_R=3·{10}^{13}$ GeV, we find the stability condition

$$\begin{array}{c}\hfill M>3·{10}^4\mathrm{GeV}{\left(\frac{qT}{\mathrm{GeV}}\right)}^{1/2},\end{array}$$

(21)

which is fulfilled for all interesting values of M.

The value of $\lambda$ depends upon the relation between $\kappa$ and $M_R$. If $\kappa \sim M_R$, then the frequency of the oscillations of curvature is of the order of $M_R$ and $\left|\lambda \right|\sim M_R$. If $\kappa \gg M_R$, then there are two possible solutions $\left|\lambda \right|\sim M_R$ and $\left|\lambda \right|\sim \kappa (\kappa /M_R)\gg M_R$. High-frequency oscillations of R would lead to an efficient gravitational particle production and, as a result, to a damping of the oscillations.

In fact, both conditions $\ddot{R}/\dot{R}\gg H$ and $\dot{R}/R\gg H$ are essentially the same at the stage of exponential rise of $R\sim exp\left(\lambda t\right)$, since the r.h.s in both cases is just $\lambda$. Since H drops down with decreasing temperature and $\lambda \sim (M_{Pl}{M}^2/\left(T{M}_R\right))$ on the opposite rises up, these conditions should be true at sufficiently small temperatures.

3.2. Fermionic Case

In this section, we consider the case when a baryonic number is carried by fermions. The gravitational part of the action has the form as in Equation (4), while the fermionic part of the action is the same as in Refs. [10,18]:

$$\begin{array}{ccc}\hfill \mathcal{L}[Q,L]& =& \frac{i}{2}(\overline{Q}{\gamma }^{\mu }{\nabla }_{\mu }Q-{\nabla }_{\mu }\overline{Q}{\gamma }^{\mu }Q)-m_Q\overline{Q}Q\hfill \\ \hfill & +& \frac{i}{2}(\overline{L}{\gamma }^{\mu }{\nabla }_{\mu }L-{\nabla }_{\mu }\overline{L}{\gamma }^{\mu }L)-m_L\overline{L}L\hfill \\ & +& \frac{g}{m_X^2}\left[\left(\overline{Q}{Q}^c\right)\left(\overline{Q}L\right)+\left({\overline{Q}}^{c}Q\right)\left(\overline{L}Q\right)\right]+\frac{d}{M^2}\left({\partial }_{\mu }R\right)J^{\mu }+{\mathcal{L}}_{matt},\hfill \end{array}$$

(22)

where Q is the quarklike field with nonzero baryonic number $B_Q$, $Q^c$ is the charged conjugated quark operator, L is another fermionic field (lepton), and ${\nabla }_{\mu }$ is the covariant derivative of the Dirac fermions in tetrad formalism. The quark current is $J^{\mu }=B_Q\overline{Q}{\gamma }^{\mu }Q$ with ${\gamma }^{\mu }$ being the curved space gamma matrices, and ${\mathcal{L}}_{matt}$ describes all other forms of matter. The four-fermion interaction between quarks and leptons is introduced to ensure the necessary nonconservation of the baryon number with $m_X$ being a constant parameter with dimension of mass and g being a dimensionless coupling constant. In the term describing the interaction of the baryonic current of fermions with the derivative of the curvature scalar, M is a constant parameter with a dimension of mass and $d=\pm 1$ is a dimensionless coupling constant which is introduced to allow for an arbitrary sign of the above expression.

Gravitational equations of motion with an account of $R^2/M_R^2$-term in analogy with Equation (7) take the form:

$$\begin{array}{ccc}\hfill & & \frac{M_{Pl}^2}{8\pi }\left[R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R-\frac{1}{3M_R^2}\left(R_{\mu \nu }-\frac{1}{4}{g}_{\mu \nu }R+g_{\mu \nu }{D}^2-D_{\mu }{D}_{\nu }\right)R\right]\hfill \\ \hfill & =& \frac{g_{\mu \nu }}{2}\frac{g}{m_X^2}\left[\left(\overline{Q}{Q}^c\right)\left(\overline{Q}L\right)+\left({\overline{Q}}^{c}Q\right)\left(\overline{L}Q\right)\right]\hfill \\ \hfill & +& \frac{i}{4}\left[\overline{Q}({\gamma }_{\mu }{\nabla }_{\nu }+{\gamma }_{\nu }{\nabla }_{\mu })Q-({\nabla }_{\nu }\overline{Q}{\gamma }_{\mu }+{\nabla }_{\mu }\overline{Q}{\gamma }_{\nu })Q\right]\hfill \\ \hfill & +& \frac{i}{4}\left[\overline{L}({\gamma }_{\mu }{\nabla }_{\nu }+{\gamma }_{\nu }{\nabla }_{\mu })L-({\nabla }_{\nu }\overline{L}{\gamma }_{\mu }+{\nabla }_{\mu }\overline{L}{\gamma }_{\nu })L\right]\hfill \\ & -& \frac{2d}{M^2}\left[R_{\mu \nu }+g_{\mu \nu }{D}^2-D_{\mu }{D}_{\nu }\right]D_{\alpha }{J}^{\alpha }+\frac{d}{2M^2}\left(J_{\mu }{\partial }_{\nu }R+J_{\nu }{\partial }_{\mu }R\right)+T_{\mu \nu }^{matt}.\hfill \end{array}$$

(23)

Taking the trace of Equation (23) with an account of fermion equations of motion, we obtain:

$$\begin{array}{c}\hfill -\frac{M_{Pl}^2}{8\pi }\left(R+\frac{1}{M_R^2}{D}^{2}R\right)=m_Q\overline{Q}Q+m_L\overline{L}L+\frac{2g}{m_X^2}\left[\left(\overline{Q}{Q}^c\right)\left(\overline{Q}L\right)+\left({\overline{Q}}^{c}Q\right)\left(\overline{L}Q\right)\right]\\ \hfill -\frac{2d}{M^2}(R+3D^2)D_{\alpha }{J}^{\alpha }+T_{matt},\end{array}$$

(24)

where $T_{matt}$ is the trace of the energy momentum tensor of all other fields. In the early universe when various species are relativistic, we can take $T_{matt}=0$. The average expectation value of the quark–lepton interaction term proportional to g is also small, so the contribution of all matter fields may be neglected and hence the only term which remains in the r.h.s. of Equation (24) is that proportional to $D_{\alpha }{J}^{\alpha }$.

A higher-order differential equation for R is obtained after we substitute the current divergence, $D_{\alpha }{J}^{\alpha }$, calculated from the kinetic equation in the external field R [18], into Equation (24). For the spatially homogeneous case,

$$\begin{array}{c}\hfill D_{\alpha }{J}^{\alpha }=({\partial }_t+3H)n_B=I_B^{coll},\end{array}$$

(25)

where the collision integral, $I_B^{coll}$, in the lowest order of perturbation theory is equal to:

$$\begin{array}{ccc}& & I_B^{coll}=-3B_q{\left(2\pi \right)}^4\int d{\nu }_{q_1,q_2}d{\nu }_{{\overline{q}}_3,l_4}{\delta }^4(q_1+q_2-q_3-l_4)\hfill \\ & & \left[|A(q_1+q_2\to {\overline{q}}_3+l_4){|}^2{f}_{q_1}{f}_{q_2}-{\left|A({\overline{q}}_3+l_4\to q_1+q_2)\right|}^2{f}_{{\overline{q}}_3}{f}_{l_4}\right].\hfill \end{array}$$

(26)

Here, $A(a\to b)$ is the amplitude of the transition from state a to state b, $B_Q$ is the baryonic number of quark, $f_a$ is the phase-space distribution (the occupation number), and

$$\begin{array}{c}\hfill d{\nu }_{q_1,q_2}=\frac{d^3{q}_1}{2E_{q_1}{\left(2\pi \right)}^3}\frac{d^3{q}_2}{2E_{q_2}{\left(2\pi \right)}^3},\end{array}$$

(27)

where $E_q=\sqrt{q^2+m^2}$ is the energy of a particle with three-momentum q and mass m. The element of the phase space of the final particles, $d{\nu }_{{\overline{q}}_3,l_4}$, is defined analogously.

We choose such representation of the quark operator, Q, for which the interaction of the baryonic current with the derivative of the curvature scalar in Equation (22) vanishes but reappears in the quark–lepton interaction term:

$$\begin{array}{c}\hfill \frac{2g}{m_X^2}\left[e^{-3id{B}_{Q}R/M^2}\left(\overline{Q}{Q}^c\right)\left(\overline{Q}L\right)+e^{3id{B}_{Q}R/M^2}\left({\overline{Q}}^{c}Q\right)\left(\overline{L}Q\right)\right].\end{array}$$

(28)

We make the simplifying assumption that the evolution of R can be approximately described by the law

$$\begin{array}{c}\hfill R\left(t\right)\approx R\left(t_0\right)+(t-t_0)\dot{R}.\end{array}$$

(29)

We assume that $\dot{R}\left(t\right)$ slowly changes at the characteristic time scale of the reactions, which contribute to the collision integral (26), and so we can approximately take $\dot{R}\approx const$.

According to the rules of quantum field theory, the reaction probability is given by the square of the integral over space and time of the amplitude of the corresponding process. In the case of a time-independent interaction, it leads to the energy conservation, $\Sigma E_{in}=\Sigma E_{fin}$. If the interaction depends upon time, the energy evidently is nonconserved and in our case, e.g., for the reaction $q_1+q_2\to {\overline{q}}_3+l_4$, the energy balance has the form:

$$\begin{array}{c}\hfill E\left(q_1\right)+E\left(q_2\right)=E\left(q_3\right)+E\left(l_4\right)+3d{B}_Q\dot{R}/M^2.\end{array}$$

(30)

In kinetic equilibrium, the phase-space distribution of fermions has the form

$$\begin{array}{c}\hfill f=\frac{1}{e^{(E/T-\xi )}+1}\approx e^{-E/T+\xi },\end{array}$$

(31)

where $\xi =\mu /T$ is the dimensionless chemical potential, different for quarks, ${\xi }_q$, and leptons, ${\xi }_l$. In the thermal equilibrium case, the condition of conservation of chemical potentials is fulfilled, that is, $\Sigma {\xi }_{in}=\Sigma {\xi }_{fin}$. In particular, it demands that chemical potentials of particles and antiparticles are equal by magnitude and have opposite signs: $\xi =-\overline{\xi }$, as follows, e.g., from the consideration of particle–antiparticle annihilation into different numbers of photons. If energy is not conserved, due to the time-dependent $R\left(t\right)$, the conservation of chemical potentials is also broken, as we see in what follows.

We assume that $\xi \ll 1$ and hence, distribution (31) turns into:

$$\begin{array}{c}\hfill f\approx e^{-E/T}(1+\xi ).\end{array}$$

(32)

We also assume that $3d{B}_Q\dot{R}/\left(M^{2}T\right)\ll 1$ and correspondingly, the balance of chemical potentials in equilibrium for the reactions $q_1+q_2\leftrightarrow {\overline{q}}_3+l_4$ leads to:

$$\begin{array}{c}\hfill 3{\xi }_q-{\xi }_l-\frac{3d{B}_Q\dot{R}\left(t\right)}{M^{2}T}=0.\end{array}$$

(33)

Following Ref. [18], we express

$$\begin{array}{c}\hfill n_B\approx \frac{g_s{B}_Q}{6}{\xi }_q{T}^3,\end{array}$$

(34)

where $g_s$ is the number of quark spin states. Since we are studying the instability of R whose timescale is presumed to be much smaller than the expansion rate of the Universe, we approximate

$$\begin{array}{c}\hfill D_{\alpha }{J}^{\alpha }\approx {\dot{n}}_B\approx \frac{g_s{B}_Q}{6}{\dot{\xi }}_q{T}^3\end{array}$$

(35)

$$\begin{array}{c}\hfill \approx \frac{g_s{B}_Q}{6}{\dot{\xi }}_q^{eq}{T}^3,\end{array}$$

(36)

${\xi }_q^{eq}$ is obtained from Equation (33), using the conservation of the sum of baryonic and leptonic numbers, which implies ${\xi }_l=-{\xi }_q/3$. Then,

$$\begin{array}{c}\hfill {\xi }_q^{eq}=\frac{9d{B}_Q\dot{R}\left(t\right)}{10M^{2}T}.\end{array}$$

(37)

Substituting Equation (37) in Equation (36) and neglecting the $\dot{T}$-term, Equation (24) gives the following fourth-order differential equation for the curvature scalar:

$$\begin{array}{c}\hfill \frac{d^{4}R}{d{t}^4}+\frac{{\kappa }_f^4}{M_R^2}\frac{d^{2}R}{d{t}^2}+{\kappa }_f^{4}R=0,\end{array}$$

(38)

where

$$\begin{array}{c}\hfill {\kappa }_f^4=\frac{5M_{Pl}^2{M}^4}{36\pi g_s{B}_Q^2{T}^2}.\end{array}$$

(39)

Once again, we consider terms containing R as small with respect to the terms containing $\ddot{R}$. The value of ${\kappa }_f$ is only slightly numerically different from $\kappa$ in Equation (18) and has the same dependence upon the essential parameters, so the solutions of Equations (17) and (38) practically coincide.

4. Discussion

We showed that as discovered in Refs. [17,18], the exponential instability of the curvature scalar inherent to the mechanism of gravitational baryogenesis could be successfully cured in modified gravity. There is an immense number of papers dedicated to gravitational baryogenesis but in none of them the problem of instability was considered. The length of our submission does not allow us the quote them all. A complete list of the literature is a task for a review paper.

The special form of gravity modification by the introduction of an $R^2$ term into the canonical Hilbert–Einstein action of general relativity was explored as a workable mechanism. However, the stabilized asymptotic value of R was extremely large and together with possibly successful baryogenesis would still strongly perturb canonical cosmology. At the present stage, a comparison with the data does not make much sense but further development indicated below may result in a model leading to a realistic cosmology and tested by astronomical observations. Possible ways out of this problem could either be a more complicated model of $F\left(R\right)$ gravity or a proper account of particle production created by the high-frequency oscillations of $R\left(t\right)$. Both options open interesting possibilities for future research.