Cosmological consequences of a scalar field with oscillating equation of state. IV. Primordial nucleosynthesis and the deuterium problem

TL;DR

This paper investigates how a scalar field with an oscillating equation of state affects primordial nucleosynthesis, demonstrating the model's compatibility with observed light element abundances and proposing a solution to the deuterium problem.

Contribution

It introduces a stepwise scalar field model compatible with BBN constraints and provides a public code for its analysis, offering a new perspective on early Universe evolution.

Findings

The model survives BBN constraints.

Early dark energy can address the deuterium problem.

Universe's radiation era may evolve chaotically rather than oscillating.

Abstract

We study the primordial nucleosynthesis (BBN) in the stepwise scalar field model proposed by Ti\'an [arXiv:1912.13208, Phys. Rev. D 101, 063531 (2020)], which provides a multiaccelerating Universe solution to the cosmological coincidence problem and predicts that the scalar field may be non-negligible even in the early Universe. The observed abundances of the light elements can be used to constrain the energy density of the scalar field during the BBN era. We present a public \texttt{Matlab} code to implement the BBN calculation in the stepwise scalar field model. We show that the model can survive the BBN constraints. In particular, this model incorporates a new solution to the possible deuterium problem: very early dark energy that appears at the end of BBN. In addition, the BBN constraints, along with constraints from the cosmic late-time acceleration, suggest that the Universe in the radiation era evolves in a chaotic accelerating manner, rather than an oscillating scaling manner.

Tables

Table 1. Table 1: Nuclides considered in the nuclear network, together with mass excess Δ ​ m i Δ subscript 𝑚 𝑖 \Delta m_{i} , mass number A i subscript 𝐴 𝑖 A_{i} , proton number Z i subscript 𝑍 𝑖 Z_{i} and spin s i subscript 𝑠 𝑖 s_{i} . The data comes from [ 54 ] . Mass excess is given in MeV MeV {\rm MeV} . Note that the nuclide mass m i = A i ​ m u + Δ ​ m i − Z i ​ m e subscript 𝑚 𝑖 subscript 𝐴 𝑖 subscript 𝑚 𝑢 Δ subscript 𝑚 𝑖 subscript 𝑍 𝑖 subscript 𝑚 𝑒 m_{i}=A_{i}m_{u}+\Delta m_{i}-Z_{i}m_{e} , where m u subscript 𝑚 𝑢 m_{u} is the atomic mass constant.

No.	1	2	3	4	5	6	7	8	9
Nuclide	$n$	$p$	$D$	$T$	${}_{}{}^3\mathrm{He}$	${}_{}{}^4\mathrm{He}$	${}_{}{}^6\mathrm{Li}$	${}_{}{}^7\mathrm{Li}$	${}_{}{}^7\mathrm{Be}$
$\Delta m_i$	$8.0713$	$7.2890$	$13.1357$	$14.9498$	$14.9312$	$2.4249$	14.0869	$14.9071$	$15.7690$
$A_i$	1	1	2	3	3	4	6	7	7
$Z_i$	0	1	1	1	2	2	3	3	4
$s_i$	$1/2$	$1/2$	$1$	$1/2$	$1/2$	$0$	1	$3/2$	$3/2$

Table 2. Table 2: The nuclear reaction network implemented in BBNLab . The n ↔ p ↔ 𝑛 𝑝 n\leftrightarrow p reaction includes three weak reactions (see Eq. (37) in [ 55 ] ).

No.	Reaction	No.	Reaction	No.	Reaction	No.	Reaction
01	$n\leftrightarrow p$	09	$T+D\leftrightarrow n+{}_{}{}^4\mathrm{He}$	17	${}_{}{}^3\mathrm{He}+{}_{}{}^3\mathrm{He}\leftrightarrow p+p+{}_{}{}^4\mathrm{He}$	25	${}_{}{}^6\mathrm{Li}+D\leftrightarrow n+{}_{}{}^7\mathrm{Be}$
02	$p+n\leftrightarrow \gamma +D$	10	$T+T\leftrightarrow n+n+{}_{}{}^4\mathrm{He}$	18	${}_{}{}^4\mathrm{He}+D\leftrightarrow \gamma +{}_{}{}^6\mathrm{Li}$	26	${}_{}{}^6\mathrm{Li}+D\leftrightarrow p+{}_{}{}^7\mathrm{Li}$
03	$D+n\leftrightarrow \gamma +T$	11	${}_{}{}^3\mathrm{He}+n\leftrightarrow \gamma +{}_{}{}^4\mathrm{He}$	19	${}_{}{}^4\mathrm{He}+T\leftrightarrow \gamma +{}_{}{}^7\mathrm{Li}$	27	${}_{}{}^7\mathrm{Li}+p\leftrightarrow {}_{}{}^4\mathrm{He}+{}_{}{}^4\mathrm{He}$
04	$D+p\leftrightarrow \gamma +{}_{}{}^3\mathrm{He}$	12	${}_{}{}^3\mathrm{He}+n\leftrightarrow p+T$	20	${}_{}{}^4\mathrm{He}+{}_{}{}^3\mathrm{He}\leftrightarrow \gamma +{}_{}{}^7\mathrm{Be}$	28	${}_{}{}^7\mathrm{Li}+p\leftrightarrow \gamma +{}_{}{}^4\mathrm{He}+{}_{}{}^4\mathrm{He}$
05	$D+D\leftrightarrow n+{}_{}{}^3\mathrm{He}$	13	${}_{}{}^3\mathrm{He}+D\leftrightarrow p+{}_{}{}^4\mathrm{He}$	21	${}_{}{}^6\mathrm{Li}+n\leftrightarrow \gamma +{}_{}{}^7\mathrm{Li}$	29	${}_{}{}^7\mathrm{Li}+D\leftrightarrow n+{}_{}{}^4\mathrm{He}+{}_{}{}^4\mathrm{He}$
06	$D+D\leftrightarrow p+T$	14	${}_{}{}^3\mathrm{He}+T\leftrightarrow \gamma +{}_{}{}^6\mathrm{Li}$	22	${}_{}{}^6\mathrm{Li}+n\leftrightarrow T+{}_{}{}^4\mathrm{He}$	30	${}_{}{}^7\mathrm{Be}+n\leftrightarrow p+{}_{}{}^7\mathrm{Li}$
07	$T\to {\overline{\nu }}_e+e^{-}+{}_{}{}^3\mathrm{He}$	15	${}_{}{}^3\mathrm{He}+T\leftrightarrow D+{}_{}{}^4\mathrm{He}$	23	${}_{}{}^6\mathrm{Li}+p\leftrightarrow \gamma +{}_{}{}^7\mathrm{Be}$	31	${}_{}{}^7\mathrm{Be}+n\leftrightarrow {}_{}{}^4\mathrm{He}+{}_{}{}^4\mathrm{He}$
08	$T+p\leftrightarrow \gamma +{}_{}{}^4\mathrm{He}$	16	${}_{}{}^3\mathrm{He}+T\leftrightarrow n+p+{}_{}{}^4\mathrm{He}$	24	${}_{}{}^6\mathrm{Li}+p\leftrightarrow {}_{}{}^3\mathrm{He}+{}_{}{}^4\mathrm{He}$	32	${}_{}{}^7\mathrm{Be}+D\leftrightarrow p+{}_{}{}^4\mathrm{He}+{}_{}{}^4\mathrm{He}$

Equations

$$H^2=\frac{8\pi G}{3}({\rho }_{\gamma }+{\rho }_{\nu }+{\rho }_e+{\rho }_b+{\rho }_{\varphi }),$$

$${\rho }_{\gamma }$$

$${\rho }_{\nu }$$

$${\rho }_e$$

$$\left(\frac{1}{e^{x-{\varphi }_e}+1}+\frac{1}{e^{x+{\varphi }_e}+1}\right)\mathrm{d}x,$$

$$=\frac{2m_e^4{c}^3}{{\pi }^2{\hslash }^3}\sum _{i=1}^{\infty }(-1{)}^{i+1}\cosh (i{\varphi }_e)M(i{z}_{\gamma }),$$

$${\rho }_b$$

$${\rho }_{\varphi }$$

$${\dot{\rho }}_{\mathrm{plasma}}+3H\left({\rho }_{\mathrm{plasma}}+p_{\mathrm{plasma}}/c^2\right)=0,$$

$$p_{\gamma }$$

$$p_e$$

$$\left(\frac{1}{e^{x-{\varphi }_e}+1}+\frac{1}{e^{x+{\varphi }_e}+1}\right)\mathrm{d}x,$$

$$=\frac{2m_e^4{c}^5}{{\pi }^2{\hslash }^3{z}_{\gamma }}\sum _{i=1}^{\infty }\frac{(-1{)}^{i+1}}{i}\cosh (i{\varphi }_e)L(i{z}_{\gamma }),$$

$$p_b$$

$${\dot{T}}_{\nu }+H{T}_{\nu }=0.$$

$$\ddot{\varphi }+3H\dot{\varphi }+c^2{V}^{\prime }=0,$$

$${\dot{n}}_b+3H{n}_b=0.$$

$$\frac{n_{e^{-}}-n_{e^{+}}}{n_{\gamma }}=\eta \sum _{i=1}^{N_{\mathrm{nuclei}}}{Z}_i{Y}_i,$$

$$\frac{n_{e^{-}}-n_{e^{+}}}{n_{\gamma }}$$

$$=\frac{z_{\gamma }^3}{\zeta (3)}\sum _{i=1}^{\infty }(-1{)}^{i+1}\sinh (i{\varphi }_e)L(i{z}_{\gamma }),$$

$$N_{i_1}{(}^{A_{i_1}}{Z}_{i_1})+N_{i_2}{(}^{A_{i_2}}{Z}_{i_2})+\cdots +N_{i_p}{(}^{A_{i_p}}{Z}_{i_p})⟷N_{j_1}{(}^{A_{j_1}}{Z}_{j_1})+N_{j_2}{(}^{A_{j_2}}{Z}_{j_2})+\cdots +N_{j_q}{(}^{A_{j_q}}{Z}_{j_q}),$$

$${\dot{Y}}_{i_1}=\sum _{\text{reactions with }{i}_1}{N}_{i_1}\left({\Gamma }_{j_1\cdots j_q\to i_1\cdots i_p}\frac{Y_{j_1}^{N_{j_1}}\cdots Y_{j_q}^{N_{j_q}}}{N_{j_1}!\cdots N_{j_q}!}-{\Gamma }_{i_1\cdots i_p\to j_1\cdots j_q}\frac{Y_{i_1}^{N_{i_1}}\cdots Y_{i_p}^{N_{i_p}}}{N_{i_1}!\cdots N_{i_p}!}\right)\equiv {\Gamma }_{i_1},$$

$${\kappa }_1{\dot{\varphi }}_e+{\kappa }_2{\dot{z}}_{\gamma }={\kappa }_3,$$

$${\kappa }_1$$

$${\kappa }_2$$

$${\kappa }_3$$

$${\kappa }_4{\dot{\varphi }}_e+{\kappa }_5{\dot{z}}_{\gamma }={\kappa }_6,$$

$${\kappa }_4$$

$${\kappa }_5$$

$$+\sum _{i=1}^{\infty }(-1{)}^{i+1}i\cosh (i{\varphi }_e)M^{\prime }(i{z}_{\gamma }),$$

$${\kappa }_6$$

$$-n_b\sum _{i=1}^{N_{\mathrm{nuclei}}}{\Gamma }_i\left(m_i+\frac{3m_e}{2z_{\gamma }}\right)],$$

$${\dot{\varphi }}_e$$

$${\dot{z}}_{\gamma }$$

$$V(\varphi )=V_0\exp (-\lambda \varphi ),$$

$$V(\varphi )=V_0\exp \left[-\frac{{\lambda }_1+{\lambda }_2}{2}\varphi -\frac{\alpha ({\lambda }_1-{\lambda }_2)}{2}\sin \frac{\varphi }{\alpha }\right],$$

$${\Gamma }_{n\to p}$$

$$+K{\int }_1^{\infty }\frac{x(x-q{)}^2\sqrt{x^2-1}}{(1+e^{-x{z}_{\gamma }})\left[1+e^{(x-q)z_{\nu }}\right]}\mathrm{d}x,$$

$${\Gamma }_{p\to n}$$

$$+K{\int }_1^{\infty }\frac{x(x+q{)}^2\sqrt{x^2-1}}{(1+e^{-x{z}_{\gamma }})\left[1+e^{(x+q)z_{\nu }}\right]}\mathrm{d}x,$$

$${\Gamma }_{i_1\cdots i_p\to j_1\cdots j_q}={\left[\frac{n_b{m}_u}{\mathrm{g}/{\mathrm{cm}}^3}\right]}^{(N_{i_1}+\cdots +N_{i_p})-1}\cdot R_{\mathrm{t}}\cdot [{\mathrm{s}}^{-1}],$$

Full text

I Introduction

The current standard cosmological model holds that there are two stages of accelerating expansion in the Universe [1, 2]. One is the early inflation, which is proposed to solve the horizon and flatness problems [3, 4, 5]. The other is the cosmic late-time acceleration, which was confirmed by observations of supernovae [6, 7]. Theoretically, there are many dark energy models to explain the late-time acceleration, such as the cosmological constant (the standard model), quintessence [8, 9, 10, 11], phantom [12], $k$-essence [13, 14], quintom [15, 16, 17]. However, these models share a common problem — the cosmological coincidence problem, which states why the energy density of dark energy is now on the same order of magnitude as that of matter [18, 19]. To solve this problem, Dodelson et al. [20] proposed that the Universe may have experienced periodic accelerating expansion in the past and the current acceleration is just a natural continuation. Meanwhile they implemented this scenario with a canonical scalar field model. This multiacceleration scenario inspired many follow-up studies, including model extensions [21, 22, 23, 24, 25] and confronting observations [26, 27, 28, 29]. In particular, inflation can be naturally unified into this scenario [21, 23, 24], which makes the model more attractive. We proposed a new canonical scalar field model to realize the multiacceleration scenario [30]. Hereafter we call this type of model a stepwise scalar field model since the potential is similar to the staircase (see Fig. 1 in [30] for an illustration). Furthermore, we found chaos and period-doubling bifurcations in our model [31], which were not identified in previous models. Our model also naturally unifies inflation and predicts an extremely small tensor-to-scalar ratio [32], which is consistent with current observations [33].

As the literal meaning of multiacceleration expresses, there should also be stages of accelerating expansion in the early Universe, that is, stages dominated by dark energy. How to detect dark energy in the early Universe? One approach is the global cosmological parameter constraints as did in [26, 27, 29]. This approach has strong model dependencies and thus is indirect. Another approach is the primordial nucleosynthesis (Big Bang nucleosynthesis, abbreviated BBN, see [34, 35] for reviews), which could probe the Universe with temperatures between $10^{12}\,{\rm K}$ and $10^{8}\,{\rm K}$. Given the great success (mainly about the helium-4 and deuterium abundances) of the standard BBN theory at the time, Dodelson et al. [20] attempted to completely hide the scalar field during the BBN era in their model. Here standard means there is no dark energy. There are other similar works using BBN to constrain the upper limit of the early dark energy density in specific models [36, 37]. However, with the improvement of the measurement uncertainties about nuclear reaction rates and cosmological parameters, the standard BBN theory now seems to be broken. In addition to the well-known primordial lithium problem [38], there are inconsistencies between the theoretical prediction and observational abundance about deuterium [39, 40]. The essence of this deuterium problem remains unclear. On the one hand, it is still debatable whether the deuterium problem really exists [41, 42]. On the other hand, if the deuterium problem does exist, its solution may hint at new physics, such as time-varying fundamental constants [43]. The deuterium and lithium problems may also shed light on the possible existence of dark energy in the early Universe.

In this paper, we use BBN to test the stepwise scalar field model proposed by us [30]. We focus on the following two issues:

Is it possible to completely hide the scalar field in our model as did in [20]? 2. 2.

Is there some kind of early dark energy profiles in our model to solve the deuterium problem, or even the lithium problem?

Meanwhile we require that the model under the same parameter settings can explain the observed cosmic late-time acceleration. In this paper, we intend to provide case studies rather than global parameter constraints.

This paper is organized as follows. Section II summarizes the BBN theory, including how to incorporate the stepwise scalar field into the BBN calculation. Section III describes our numerical results, including discussions about the BBN results of the standard model, and models with exponential scalar field and stepwise scalar field. Conclusions are presented in Sec. IV. Throughout this paper, we adopt the SI Units and retain all physical constants unless otherwise stated. All reported uncertainties represent $68\%$ confidence intervals.

II Basic theory

BBN calculation requires one to solve a nuclear reaction network in the expanding Universe at temperatures between $10^{12}\,{\rm K}$ and $10^{8}\,{\rm K}$ (equivalently energy scales between $100\,{\rm MeV}$ and $0.01\,{\rm MeV}$). The main BBN physics was summarized in [44, 45], and the results are widely used in current BBN codes, e.g., NUC123 [46, 47], PArthENoPE [48, 49, 50], AlterBBN [51, 52] and PRIMAT [35, 39]. Here we summarize the BBN theory implemented in this paper, which include all the key physics. The results are integrated into a public Matlab code, which is named BBNLab and is available at GitHub111https://github.com/tshuxun/BBNLab. Note that no physics in this section is original. However, this summary, especially the explicit expressions given in Eqs. (13) and (15), may be useful for beginners and presents a clear basis for our future BBN work.

II.1 Evolution equations

The Universe is assumed to be described by the flat Friedmann-Lemaître-Robertson-Walker metric ${\rm d}s^{2}=-c^{2}{\rm d}t^{2}+a^{2}({\rm d}x^{2}+{\rm d}y^{2}+{\rm d}z^{2})$, where $a=a(t)$ is the cosmic scale factor and $c$ is the speed of light, and contains photons, neutrinos, electrons, positrons, baryons, and a canonical scalar field. Hereafter we use the subscripts $\{\gamma,\nu,e^{-},e^{+},b,\phi\}$ to denote these ingredients respectively, use plasma to refer {photons + electrons + positrons + baryons}, use the subscript matter to refer {plasma + neutrinos}, and use the subscript total to refer {matter + $\phi$}. During BBN era, we assume that the plasma is in thermal equilibrium, while the neutrino gas is thermally decoupled from the plasma222Note that neutrinos are in thermal equilibrium with the plasma through weak interactions when the energy scale is higher than $2\,{\rm MeV}$ [35]. As the temperature drops, neutrino decoupling and $e^{\pm}$ annihilation occur one after the other. Neither process is instantaneous and there is a temporal overlap between them. This overlap results in a neutrino heating process from the $e^{\pm}$ annihilation [53]. For simplicity, we assume that all the influences of this neutrino heating process on BBN are reflected in a constant $N_{\nu}^{\rm eff}$ [see Eq. (2b)]. After considering $N_{\nu}^{\rm eff}$, the neutrino gas can be regarded as expanding adiabatically. In addition, the scalar field that appears during the $e^{\pm}$ annihilation may slightly change the value of $N_{\nu}^{\rm eff}$ as it affects the Hubble expansion rate. In this paper, we ignore this possible modification and would like to leave this issue to the future.. The scalar field interacts with other species only through gravity. Nuclides considered in the nuclear network and their basic properties are listed in Table 1. Substituting the metric into the Einstein field equations gives the Friedmann equation

[TABLE]

where $H\equiv\dot{a}/a$ is the Hubble parameter, $\dot{}\equiv{\rm d}/{\rm d}t$, $G$ is the Newtonian constant of gravitation, $\rho_{e}=\rho_{e^{-}}+\rho_{e^{+}}$. The mass densities (${\rm energy\ densities}/c^{2}$) are given by

[TABLE]

where $k_{\rm B}$ is the Boltzmann constant, $\hbar$ is the reduced Planck constant, $T_{\gamma}$ is the plasma temperature, $T_{\nu}$ is the neutrino temperature, $N_{\nu}^{\rm eff}=3.046$ is the effective number of neutrinos [53], $\infty$ means positive infinity, $z_{\gamma}\equiv m_{e}c^{2}/(k_{\rm B}T_{\gamma})$, $\phi_{e}\equiv\mu_{e}/(k_{\rm B}T_{\gamma})$, $m_{e}$ is the mass of electron, $\mu_{e}$ is the chemical potential of the electron gas, $M(x)$ is a special function defined in Appendix A, $m_{i}$ is the mass of nuclide $i$ (see Table 1), $Y_{i}\equiv n_{i}/n_{b}$, $n_{i}$ is the number density of nuclide $i$, $n_{b}$ is the total number density of baryons, $N_{\rm nuclei}$ is the total species of nuclides, $V(\phi)$ is the potential of the scalar field (see Sec. II.1.1). The conventions about the scalar field and also $\rho_{\phi}$ are consistent with [30]. The black-body spectrum is used to derive Eq. (2a). The Fermi-Dirac distribution with zero neutrino mass and chemical potential is used to derive Eq. (2b). This assumption is reasonable as the neutrino mass is much smaller than ${\rm MeV}/c^{2}$ and the zero chemical potential is consistent with current observational constraints [35]. The Fermi-Dirac distribution with $\mu_{e}+\mu_{e^{+}}=0$, where $\mu_{e^{+}}$ is the positron chemical potential, is used to derive the first equality of Eq. (2c). Here $\mu_{e}+\mu_{e^{+}}=0$ corresponds to the chemical equilibrium of the $e^{-}+e^{+}\leftrightarrow\gamma+\gamma$ reaction. The series expansion in the second equality of Eq. (2c) facilitates numerical calculations. The ideal gas model is used to derive Eq. (2d). The internal degeneracy of particles has been taken into account in the above results.

As we discussed before, aside from the gravitational interaction, the thermodynamics of the plasma and neutrinos, and the dynamics of the scalar field are independent of each other. Therefore, energy conservation provides three other evolution equations. For plasma, we obtain

[TABLE]

where $\rho_{\rm plasma}=\rho_{\gamma}+\rho_{e}+\rho_{b}$, $p_{\rm plasma}=p_{\gamma}+p_{e}+p_{b}$, $p_{e}=p_{e^{-}}+p_{e^{+}}$, and the pressures are given by

[TABLE]

where $L(x)$ is a special function defined in Appendix A. The basic theory and assumptions used to derive the above results are discussed in the previous paragraph. For neutrinos, considering its pressure $p_{\nu}=\rho_{\nu}c^{2}/3$, we obtain

[TABLE]

Integrating the above equation gives $T_{\nu}\propto a^{-1}$. For the scalar field, the equation of energy conservation is exactly its equation of motion [30]

[TABLE]

where $V^{\prime}\equiv{\rm d}V/{\rm d}\phi$.

There are two other equations that are closely related to baryons. One is the equation of baryon number conservation

[TABLE]

The other is a constraint equation given by the electrical neutrality of the Universe, which can be written as

[TABLE]

where $\eta\equiv n_{b}/n_{\gamma}$, and $n_{e^{-}}$, $n_{e^{+}}$, $n_{\gamma}$ are the number densities of relevant ingredients, $n_{\gamma}=2\zeta(3)m_{e}^{3}c^{3}/(\pi^{2}\hbar^{3}z_{\gamma}^{3})$ and

[TABLE]

where $\zeta(x)$ is the Riemann zeta function and $\zeta(3)\approx 1.2021$. So far, roughly speaking, the above results determine the thermodynamics of the Universe: There are six equations $\{$Eqs. (1), (3) and (5–8)$\}$ and six variables $\{a,T_{\gamma},T_{\nu},\phi_{e},n_{b},\phi\}$. Note that it is unnecessary to find the solution of $a(t)$ because other equations depend directly on $H$ rather than $a$.

The rest is the nuclear reaction network which rules the evolution of $Y_{i}$. The general form of a nuclear reaction can be written as

[TABLE]

where ${}^{A_{i}}Z_{i}$ denotes nuclide $i$ (see Table 1), $N_{i}$ is the stoichiometric coefficient. The evolution of $Y_{i}$ takes the form [35]

[TABLE]

where the sum includes all reactions involving nuclide $i_{1}$, $\Gamma_{i_{1}\cdots i_{p}\rightarrow j_{1}\cdots j_{q}}$ and $\Gamma_{j_{1}\cdots j_{q}\rightarrow i_{1}\cdots i_{p}}$ are the reaction rates (see Sec. II.1.2), and the last equality defines $\Gamma_{i}$. The nuclear reactions considered in this work are listed in Table 2. These reactions are selected from Table 2 in [48] with the criterion that the reactions involve only the nuclides listed in Table 1. As we will see, this nuclear reaction network gives reasonable BBN predictions.

The above results form complete BBN evolution equations. However, Eqs. (3) and (8) are implicit and not suitable for direct numerical calculations. Here we convert these two equations into explicit form. Multiplying Eq. (8) by $\zeta(3)/z_{\gamma}^{3}$, differentiating the result with respect to $t$, and substituting Eq. (11) into the result to eliminate $\dot{Y}_{i}$, we obtain

[TABLE]

where the coefficients are given by

[TABLE]

where $L^{\prime}(x)\equiv{\rm d}L/{\rm d}x$ and its explicit expression is given in Appendix A. Expanding $\dot{\rho}_{\rm plasma}$ by $\{\dot{\phi}_{e},\dot{z}_{\gamma},\dot{n}_{b},\dot{Y}_{i}\}$, and substituting Eqs. (7) and (11) to eliminate $\{\dot{n}_{b},\dot{Y}_{i}\}$, then Eq. (3) can be rewritten as

[TABLE]

where the coefficients are given by

[TABLE]

where $M^{\prime}(x)\equiv{\rm d}M/{\rm d}x$ and its explicit expression is given in Appendix A. The combination of Eqs. (12) and (14) gives

[TABLE]

So far, we can obtain the explicit expressions of the first time-derivative of $\{z_{\gamma},z_{\nu},\phi_{e},\eta,\phi,\dot{\phi},Y_{i}\}$, where $z_{\nu}\equiv m_{e}c^{2}/(k_{\rm B}T_{\nu})$. Here we replace $\{T_{\nu},n_{b}\}$ with $\{z_{\nu},\eta\}$ to make the variables dimensionless. In addition, $\dot{z}_{\nu}=Hz_{\nu}$ and $\dot{\eta}=3\pi^{2}\hbar^{3}z_{\gamma}^{2}n_{b}(\dot{z}_{\gamma}-Hz_{\gamma})/[2\zeta(3)m_{e}^{3}c^{3}]$ can be easily derived from the previous equations.

BBNLab adopts the variables $\{z_{\gamma},z_{\nu},\phi_{e},\eta,\phi,\dot{\phi},Y_{i}\}$ in its numerical calculations. Furthermore, we change the “time” variable from $t$ to the $e$-folding number $N$, where $N\equiv\ln(a/a_{1})$, $a_{1}$ is the scale factor at one specific time and its absolute value is meaningless. The transformation of the evolution equations directly follows ${\rm d}X/{\rm d}N=\dot{X}/H$. Hereafter we use the subscript “ini” to denote the beginning of BBN. If we set the initial temperature to $10^{12}\,{\rm K}$, then $N\in[N_{\rm ini},N_{\rm ini}+10]$ roughly corresponds to $10^{12}\,{\rm K}\leq T_{\gamma}\lesssim 10^{8}\,{\rm K}$. The exact value of $N_{\rm ini}$ have no physical effect as the BBN evolution equations do not depend on $N$ explicitly. The reason for keeping $N_{\rm ini}$ here is that it is an auxiliary value used to calculate the initial conditions of the stepwise scalar field (see Sec. II.2). For clarity, we set $N_{\rm ini}=0$ for the standard model and the exponential scalar field model (see Sec. II.1.1). The BBN evolution equations are stiff. In BBNLab, we use the Matlab/ode23s solver with suitable option structure to integrate this stiff system. We set the options as odeset( ’RelTol’,1E-10,’AbsTol’,1E-10,’MaxStep’,1E-3) for Fig. 1, and odeset(’RelTol’,1E-10,’AbsTol’, 1E-10,’MaxStep’,2E-4) for Fig. 3 and Fig. 4. Details about ode23s can be found in https://www.mathworks.com/help/matlab/ref/ode23s.html. For all series expansions involved in our numerical calculations, e.g., Eq. (2c), we adopt the results up to the 20th order.

II.1.1 Models of the scalar field

The goal of this paper is to study BBN in the framework of the stepwise scalar field model proposed in [30]. As a comparison, we will also present the BBN results for the standard model, in which no scalar fields exist, and the exponential scalar field model [56, 57]. Note that the stepwise scalar field model is inspired by the exponential scalar field model and these two models share some similar features, e.g., the scaling property [31]. Discussions on the exponential scalar field model may shed light on what kind of the stepwise scalar field model can survive from the BBN constraints.

In our conventions [30], the exponential potential reads

[TABLE]

where $\lambda$ is a dimensionless constant, $V_{0}$ is a constant in unit of ${\rm length}^{-2}$, $\phi$ is dimensionless. Scaling solution is a main mathematical property of this model [57]. In the radiation-dominated Universe, the scaling solution requires $\lambda>2$. Furthermore, previous discussions on BBN require $\lambda\gtrsim 9$ [36]. The potential of the stepwise scalar field model is written as [30]

[TABLE]

where $\lambda_{1}$, $\lambda_{2}$ and $\alpha$ are dimensionless constants, $V_{0}$ and $\phi$ follow the previous conventions. Period-doubling bifurcation together with the oscillating scaling and chaotic accelerating solutions are the main mathematical properties of this model [31]. In order to explain the cosmic late-time acceleration and solve the accompanying coincidence problem, we require $\{\lambda_{1}+\lambda_{2}>4,\,0<\lambda_{2}\lesssim 0.4,\,\alpha=\mathcal{O}(1)\}$ [30]. Furthermore, we require $\lambda_{2}=\mathcal{O}(10^{-4})$ to unify inflation [32].

II.1.2 Nuclear reaction rates

For clarity, we adopt the conventions that the forward direction of the reactions listed in Table 2 and Eq. (10) is taken as left to right. For the $n\leftrightarrow p$ reaction, weak interaction theory gives [55, 58]

[TABLE]

where $q=(m_{n}-m_{p})/m_{e}\approx 2.5310$, $m_{n}$ and $m_{p}$ are the neutron and proton mass respectively, $K=6.9503\times 10^{-4}\,{\rm s}^{-1}$ is the normalization constant computed by requiring $\Gamma_{n\rightarrow p}=1/\tau_{n}$ at low temperatures, $\tau_{n}=879.4\,{\rm s}$ is the neutron lifetime [59]. The Fermi correction discussed in [58] is ignored here. At high temperatures, we have $z_{\nu}=z_{\gamma}$ so that $\Gamma_{n\rightarrow p}=\Gamma_{p\rightarrow n}e^{qz_{\gamma}}$, which is consistent with the equilibrium abundances (see Eq. (25)). For other reactions, generally, the forward thermonuclear reaction rates in the CGS units are given in the literature [35]. The forward reaction rates read [35]

[TABLE]

where $R_{\rm t}$ is the numerical value of the thermonuclear reaction rate in CGS units. The $R_{\rm t}$ actually used in the present calculation and relevant references are presented in the code. Most of the reaction rates we use are consistent with those used in PRIMAT (Version 0.2.0) [39]. Especially, as did in [39], we adopt the most recent measurements about the $D+p\rightarrow\gamma+{{}^{3}{\rm He}}$ reaction [60]. The deuterium abundance is sensitive to this reaction [39, 40] and is the core in our discussions. The last term in Eq. (21) is the physical unit, i.e., the unit of $\Gamma_{i_{1}\cdots i_{p}\rightarrow j_{1}\cdots j_{q}}$. The inverse reaction rates read [35]

[TABLE]

where the dimensionless ratio

[TABLE]

and the constant $\beta_{i}$ is given in our code. Note that each nuclear reaction corresponds to a set of $\beta_{i}$. Uncertainties [61] of the reaction rates are not included in the present calculation. More nuclear reactions and the uncertainties will be considered in the future.

II.2 Initial conditions

Considering the neutrino-plasma thermal equilibrium, we set the initial temperatures $T_{\gamma,{\rm ini}}=T_{\nu,{\rm ini}}=10^{12}\,{\rm K}$ in BBNLab. For the baryon-to-photon number ratio, we set $\eta_{\rm ini}\in 11/4\times[10^{-10},10^{-9}]$ so that $\eta\in[10^{-10},10^{-9}]$ at the ending of BBN. Here the coefficient $11/4$ is due to the entropy transfer from $e^{\pm}$ pairs to photons [35].

The initial high temperature makes all nuclides in statistical equilibrium, which includes thermal and chemical equilibrium. The initial number density of nuclide $i$ reads

[TABLE]

where $g_{i}=2s_{i}+1$ is the spin degeneracy, $\mu_{i,{\rm ini}}$ is the initial chemical potential. The chemical potential of leptons is negligible compared to that of baryons. Therefore we obtain $\mu_{n,{\rm ini}}=\mu_{p,{\rm ini}}$ based on the $n\leftrightarrow p$ reaction. Applying Eq. (24) to neutrons and protons, we obtain

[TABLE]

Initially, neutrons and protons are much more numerous than other nucleons. Therefore, it is reasonable to adopt $n_{n,{\rm ini}}+n_{p,{\rm ini}}=n_{b,{\rm ini}}$, which in turn gives

[TABLE]

For other nuclides, we set $Y_{i,{\rm ini}}=0$ in BBNLab. Note that we turn off all nuclear reactions except the $n\leftrightarrow p$ reaction until $T_{\gamma}=10^{10}\,{\rm K}$. Once the full nuclear reaction network is turned on, the heavier nuclides reach their equilibrium abundances very quickly (see Fig. 1 for an illustration). Therefore, it is reasonable to set $Y_{i,{\rm ini}}=0$ for the heavier nuclides and the final results is independent of this setting. This strategy is also adopted by PRIMAT [35].

The initial value of $\phi_{e}$ can be obtained by numerically solving Eq. (8). Here we present a high temperature approximate solution. The right side of Eq. (8) is much less than $1$ as $\eta\ll 1$. For its left side, i.e., Eq. (II.1), one can verify that $\lim_{z_{\gamma}\rightarrow 0}z_{\gamma}^{3}L(iz_{\gamma})=2/i^{3}=\mathcal{O}(1)$, where $z_{\gamma}\rightarrow 0$ corresponds to the high temperature approximation. Therefore, Eq. (8) indicates $\sinh(i\phi_{e})=\mathcal{O}(\eta)\ll 1$ at the beginning of BBN. Substituting the Taylor expansion $\sinh x=x$ into Eq. (8), we obtain

[TABLE]

This result is consistent with Eq. (20) in [51].

The initial conditions of the scalar field, i.e., in principle, the values of $\{\phi_{\rm ini},\dot{\phi}_{\rm ini},V_{0}\}$, are calculated as follows. For the exponential scalar field model, we only consider the scaling case and assume that the scalar field reaches the scaling attractor at the beginning of BBN. Without loss of generality, we can set $\phi_{\rm ini}=0$ and use $V_{0}$ to adjust the potential energy. Then the initial conditions can be calculated as follows

$x_{1,{\rm ini}}=4/(\sqrt{6}\lambda)$,

$x_{2,{\rm ini}}=2/(\sqrt{3}\lambda)$, 2. 2.

$\Omega_{\phi,{\rm ini}}=4/\lambda^{2}$,

$\rho_{\rm matter,ini}=\rho_{\rm plasma,ini}+\rho_{\nu,{\rm ini}}$,

$\rho_{\rm total,ini}=\rho_{\rm matter,ini}/(1-\Omega_{\phi,{\rm ini}})$,

$H_{\rm ini}=(8\pi G\rho_{\rm total,ini}/3)^{1/2}$, 3. 3.

$\dot{\phi}_{\rm ini}=\sqrt{6}H_{\rm ini}x_{1,{\rm ini}}$,

$V_{0}=3(H_{\rm ini}x_{2,{\rm ini}}/c)^{2}$.

Here $x_{1}$ and $x_{2}$ are the dimensionless variables defined in the dynamical system analysis (see Eq. (7a) in [30]) and their initial values are given by the scaling solution, $\Omega_{\phi}$ is the relative energy density of the scalar field. For the stepwise scalar field model, instead of the scaling solution, we can use the numerical solution of the dynamical system to give the values of $x_{1}$ and $x_{2}$ at the beginning of BBN. The dynamical system is given by Eq. (8) with $w_{\rm m}=1/3$ in [30] and we do not repeat the equations and conventions here. Other subsequent steps are similar to the exponential case. However, one difference is that the value of $\phi_{\rm ini}$ is no longer fixed. Without loss of generality, we assume $0\leq\phi_{\rm ini}\leq 2\alpha\pi$ and use $V_{0}$ to adjust the potential energy. Details are as follows

Integrating the dynamical system in $N\in[0,N_{\rm ini}]$ gives the values of $\{x_{1,{\rm ini}},\,x_{2,{\rm ini}},\,\lambda_{\rm ini},\,\nu_{\rm ini}\}$, i.e., the values of $\{x_{1},x_{2},\lambda,\nu\}$ at $N=N_{\rm ini}$ (see Fig. 2 for an illustration). In this step, we set $\{x_{1}=0.75,x_{2}=0.5,\lambda=\lambda_{2}+2,\nu=\nu_{+}(\lambda)\}$ at $N=0$ for all following calculations. Compared to setting the values at $N=N_{\rm ini}$ directly, this step allows the system to enter the possible attractor solution and makes the desired values more natural. Note that $\lambda$ is an variable in the stepwise scalar field model, while it is a constant in the exponential scalar field model; $N_{\rm ini}$ is an auxiliary parameter and its absolute value have no physical meaning (see also discussions in Sec. II.1). 2. 2.

$\Omega_{\phi,{\rm ini}}=x_{1,{\rm ini}}^{2}+x_{2,{\rm ini}}^{2}$,

$\rho_{\rm matter,ini}=\rho_{\rm plasma,ini}+\rho_{\nu,{\rm ini}}$,

$\rho_{\rm total,ini}=\rho_{\rm matter,ini}/(1-\Omega_{\phi,{\rm ini}})$,

$H_{\rm ini}=(8\pi G\rho_{\rm total,ini}/3)^{1/2}$, 3. 3.

$\phi_{\rm ini}=\alpha\arccos\displaystyle\left[\frac{2}{\lambda_{1}-\lambda_{2}}\cdot\left(\lambda_{\rm ini}-\frac{\lambda_{1}+\lambda_{2}}{2}\right)\right]$,

if $\nu_{\rm ini}>0$ then $\phi_{\rm ini}=2\alpha\pi-\phi_{\rm ini}$, 4. 4.

$\dot{\phi}_{\rm ini}=\sqrt{6}H_{\rm ini}x_{1,{\rm ini}}$,

$V_{\rm ini}=3(H_{\rm ini}x_{2,{\rm ini}}/c)^{2}$,

$V_{0}=V_{\rm ini}\exp\displaystyle\left[\frac{\lambda_{1}+\lambda_{2}}{2}\phi_{\rm ini}+\frac{\alpha(\lambda_{1}-\lambda_{2})}{2}\sin\frac{\phi_{\rm ini}}{\alpha}\right]$.

Here the third step gives the solution of Eqs. (7b) and (7c) in [30]. Note that the function arccos returns values in the interval $[0,\pi]$ and thus the final result satisfies our initial assumption $0\leq\phi_{\rm ini}\leq 2\alpha\pi$. Considering the initial temperature, the above procedure gives $V_{0}\sim G(100\,{\rm MeV})^{4}/\hbar^{3}c^{7}\sim 10^{-80}l_{\rm P}^{-2}$, where $l_{\rm P}$ is the Planck length. This result is much smaller than that given in [32]. However, there is no inconsistency in the model. The reason is that there is a degeneracy between the absolute values of $V_{0}$ and $\phi_{\rm ini}$. Because the function $\sin$ is periodic, increasing $\phi_{\rm ini}$ by several periods together with a corresponding rescaled $V_{0}$ does not affect the potential and dynamics of the scalar field. In [32], we assumed $\phi\approx\alpha\pi$ at the beginning of inflation and obtained $V_{0}\sim 10^{-14}l_{\rm P}^{-2}$. If we want to restore these results, what we need to do is increase the above obtained $\phi_{\rm ini}$ by several periods so that $\exp[-(\lambda_{1}+\lambda_{2})\phi_{\rm ini}/2]\sim 10^{-66}$. The BBN results presented in this paper are independent of these settings.

III Numerical results

To compare the theoretical predictions and observations, we consider $Y_{\rm P}\equiv 4\,Y_{{}^{4}{\rm He}}$ (the ratio of the abundance by mass of helium-4), ${\rm D}/{\rm H}$, ${{}^{3}{\rm He}}/{\rm H}$ and ${{}^{7}{\rm Li}}/{\rm H}$ (the ratio of the abundance of the isotope to hydrogen, and $i/{\rm H}\equiv Y_{i}/Y_{p}$). Observational results are $Y_{\rm P}=0.2453\pm 0.0034$ [62], ${\rm D}/{\rm H}=(2.527\pm 0.030)\times 10^{-5}$ [63], ${{}^{3}{\rm He}}/{\rm H}<(1.1\pm 0.2)\times 10^{-5}$ [64], and ${{}^{7}{\rm Li}}/{\rm H}=(1.58_{-0.28}^{+0.35})\times 10^{-10}$ [65]. The ${{}^{3}{\rm He}}$ result is given by its abundance in our Galaxy and the result is adopted as an upper limit of the primordial abundance because the post-BBN evolution of ${{}^{3}{\rm He}}$ is unclear [64, 35]. As stated in Sec. II.2, we compute the abundances for a given scalar field model over the range $\eta_{\rm today}\in[10^{-10},10^{-9}]$, which covers the Planck 2018 result $\eta_{\rm today}=(6.1374\pm 0.0383)\times 10^{-10}$ [66, 39]. Here the subscript “today” means the present epoch or equivalently the ending of BBN.

In order to check the correctness of our code, we present the equilibrium abundances. Omitting the subscript “ini” in Eq. (26) gives the equilibrium abundances of neutron and proton. This result is available at $T_{\gamma}\sim 10^{12}\,{\rm K}$. For other nuclides, the equilibrium abundances are given by [35]

[TABLE]

where the binding energy $B_{i}\equiv[Z_{i}m_{p}+(A_{i}-Z_{i})m_{n}-m_{i}]\times c^{2}$, and the relevant values are given in Table 1. If the temperature is low enough that the neutrons and protons are not in thermal equilibrium, e.g., $T_{\gamma}\sim 10^{10}\,{\rm K}$, then $Y_{i,{\rm eq}}$ should be calculated based on the exact values of $Y_{n}$ and $Y_{p}$, rather than their equilibrium values.

III.1 Exponential scalar field

Figure 1 presents the main BBN results for the exponential scalar field model. The predicted abundances as a function of $\eta_{\rm today}$ and the observed results are plotted in the left part. The case of $\lambda=\infty$ corresponds to the standard BBN model, at which $\Omega_{\phi}=0$. Compared to PRIMAT, our code gives a slightly lower ${{}^{4}{\rm He}}$ abundance, and consistent abundances of ${\rm D}$, ${{}^{3}{\rm He}}$ and ${{}^{7}{\rm Li}}$. The lower $Y_{\rm P}$ should be due to the incomplete nuclear reaction network. However, the difference is small relative to the observational uncertainty and thus is ignorable. Importantly, our code recovers the deuterium problem, which is first pointed out in [39] and states that the predicted ${\rm D}/{\rm H}$ abundance is lower than the observational constraints. Using a scalar field to solve the deuterium problem is the core in our following discussions. Decreasing $\lambda$, that is, increasing $\Omega_{\phi}$ during BBN, can significantly increase the abundances of ${{}^{4}{\rm He}}$ and ${\rm D}$. In particular, we can solve the deuterium problem with $\lambda=10$. A potential disadvantage here is that the predicted abundance of ${{}^{4}{\rm He}}$ may exceed its observed value if we consider the full PRIMAT nuclear reaction network. The case of $\lambda=8$ gives too large $Y_{\rm P}$ and ${\rm D}/{\rm H}$ and will be ruled out by the observations. This is consistent with the result given by [36] ($\lambda>9$ at $2\sigma$ C.L.). The abundances of ${{}^{3}{\rm He}}$ and ${{}^{7}{\rm Li}}$ remain almost unchanged when $\lambda$ changes from infinity to $8$. Consistent with conventional cognition [37], our results show that the exponential scalar field cannot solve the primordial lithium problem [38].

The case of $\lambda=10$ deserves more discussion as it can solve the deuterium problem. The evolution of the isotopes as a function of $T_{\gamma}$ and their equilibrium abundance are depicted in the right part of Fig. 1. As the Universe expands (temperature decreases), the evolution of neutron and proton follow their equilibrium abundance until $T_{\gamma}\approx 2\times 10^{10}\,{\rm K}$. Heavier nuclei are taken into account when $T_{\gamma}=10^{10}\,{\rm K}$ and their evolution follow the equilibrium abundance until $T_{\gamma}\approx 6\times 10^{9}\,{\rm K}$. Therefore, it is reasonable for BBNLab to consider heavier nuclei starting from $T_{\gamma}=10^{10}\,{\rm K}$. We also plot the evolution of $\Omega_{\phi}$ and the equation of state (EOS) $w_{\rm matter}$. During BBN era, $\Omega_{\phi}$ remains almost constant ($4\%$), which corresponds to the scaling solution [57]. Small fluctuations in $\Omega_{\phi}$ are due to the fluctuations in $w_{\rm matter}$, which in turn results from the $e^{\pm}$ annihilation. Combining the left and right parts of Fig. 1, we conclude that if the relative scalar field energy density reaches $4\%$, then it can result in observable influences on the final BBN predictions. During BBN, such a scalar field may can affect the out-of-equilibrium process of neutron and proton around $T_{\gamma}\approx 2\times 10^{10}\,{\rm K}$, and the heavier nuclei evolution when $T_{\gamma}\lesssim 6\times 10^{9}\,{\rm K}$.

III.2 Stepwise scalar field

In the stepwise scalar field model, our previous work established the following unified cosmic evolution scenario: inflation [32], deflationary phase with very stiff EOS [32], radiation era with oscillating scaling solution (OSS) [31], and matter era with chaotic accelerating solution (CAS) [30, 31]. OSS during radiation era can attenuate the sensitive dependence of the late-time cosmic evolution on the early initial conditions and facilitate the cosmological parameter constraints [31]. This is a pragmatic choice rather than an observational confirmation. In principle, the Universe can also evolve as a CAS during the radiation era. The BBN discussions may shed light on the real type of cosmic evolution during the radiation era.

We start from the OSS. Roughly speaking, Sec. III.1 shows that $\Omega_{\phi}\approx 4\%$ is a critical point for the BBN constraints. Higher $\Omega_{\phi}$ would be excluded from observations, while lower $\Omega_{\phi}$ would have no observable effect. If we expect OSS in the radiation era, then we require $\lambda_{1}\gtrsim 20$ because $\Omega_{\phi}\approx 16/(\lambda_{1}+\lambda_{2})^{2}$ [30] and $\lambda_{2}\ll 1$ [32]. In particular, $\lambda_{1}\approx 20$ may can be used to solve the deuterium problem as the exponential scalar field did. In Fig. 2 (a), we plot the evolution of $\Omega_{\phi}(N)$ in the radiation and matter epochs for $\lambda_{1}=20$. Other parameter settings and calculation details can be found in the caption. The inserted Fourier plots show that, under this parameter settings, the Universe evolves as an OSS in the radiation era, and as a CAS in the matter era. In the radiation era, $\Omega_{\phi}\approx 4\%$ as we expected. However, this parameter setting is difficult to explain the cosmic late-time acceleration in the matter era. An apparent difficulty is that $\Omega_{\phi}$ cannot reach $70\%$, i.e., the current observed dark energy relative density [66]. Increasing $\lambda_{1}$ will exacerbate this problem. Increasing $\alpha$ can increase the maximum value of $\Omega_{\phi}$ in the matter era, but also make the cosmic evolution in the radiation era enter the CAS333In current parameter settings, $f_{\textsc{oss}}/2$ component already appears in the radiation era [see the left inset in Fig. 2 (a)]. Further parameter changes will cause the system to quickly enter the chaotic phase.. Therefore, if $\lambda_{1}\gtrsim 20$, then the OSS in the radiation era may be incompatible with the desired acceleration in the matter era.

What we need is a small $\Omega_{\phi}$ during BBN. To achieve this, in addition to increasing $\lambda_{1}$, we can also adjust the evolution of the Universe in the radiation era to CAS. Note that $\Omega_{\phi}$ can be very close to [math] and $1$ in a CAS, but not in an OSS [31]. Technically, decreasing $\lambda_{1}$ or increasing $\alpha$ can bring the system into the CAS phase. In Fig. 2 (b), we plot such a scenario. A key setting here is $\lambda_{1}=10$ and other settings can be found in the caption. The model might survive the BBN constraints if $\Omega_{\phi}\lesssim 4\%$ during the main BBN era. In addition, this parameter setting should also be able to explain the cosmic late-time acceleration (mainly $\Omega_{\phi}\approx 70\%$ and the EOS $w_{\phi}\approx-1$) as shown in the right part of Fig. 2 (b) and the discussions in [30, 31]. Note that, in the matter era of Fig. 2 (b), we tend to emphasize the existence of the $\Omega_{\phi}=70\%$ points, rather than requiring it to occur exactly at today, i.e., $N=N_{\rm eq}+8.13$ [31]. Once the point exists, we can adjust the initial conditions so that it occurs right now, and this does not require fine-tuning [31].

For the same parameter settings as Fig. 2 (b), we perform detailed BBN calculations with BBNLab. The results are shown in Figs. 3 and 4. In order to set the initial conditions of the scalar field at $T_{\gamma}=10^{12}\,{\rm K}$ (see the discussions in Sec. II), we need an auxiliary parameter $N_{\rm ini}$, which is labeled in the figures. The colored shaded region in Fig. 2 (b) roughly corresponds to the BBN era of the $N_{\rm ini}=70.2$ case in Fig. 3. Note that this correspondence is not strict, since $w_{\rm matter}$ (i.e., $w_{\rm m}$ in [31]) remains almost constant in the former case, while $w_{\rm matter}$ is obviously time-varying in the latter case. Looking at the $Y_{\rm P}$ and ${\rm D}/{\rm H}$ results in Figs. 3 and 4, we conclude that there do exist viable parameter space in the stepwise scalar field model to explain the BBN observations. In particular, the case of $N_{\rm ini}=70.2$ in Fig. 3 or $N_{\rm ini}=78.5$ in Fig. 4 could be used to solve the possible deuterium problem [39].

For the cases of $N_{\rm ini}=70.2$ and $N_{\rm ini}=78.5$, the evolutions of $\{Y_{i},\Omega_{\phi},w_{\rm matter}\}$ as a function of $T_{\gamma}$ are depicted in the corresponding right parts. There is one difference between the two cases worth mentioning: the $Y_{\rm P}$ predictions are different. Compared with the standard model, the case of $N_{\rm ini}=70.2$ does not change $Y_{\rm P}$ while the case of $N_{\rm ini}=78.5$ increases $Y_{\rm P}$. This difference should be due to the different values of $\Omega_{\phi}$ when protons and neutrons are out of equilibrium (see the second vertical dash-dotted line in the right parts of the figures). If $\Omega_{\phi}$ drops below $4\%$ before protons and neutrons fall out of equilibrium, then the scalar field would have no observable effect on the isotope evolutions in this period. In the case of $N_{\rm ini}=70.2$, the scalar field affects the isotope evolutions only near the end of BBN, and only has an observable effect on the deuterium abundance. This discussion also provides a simple mechanism for solving the deuterium problem: very early dark energy appears near the end of BBN. Essentially, this mechanism is the core of solving the deuterium problem with the stepwise scalar field. Finally, we admit that no solution of the lithium problem was found in the stepwise scalar field model.

IV Conclusions

The stepwise scalar field model was constructed to solve the cosmological coincidence problem [20, 30]. A key idea of this model is that the cosmic expansion has been dominated by the scalar field many times, which indicates that a non-negligible scalar field could have arisen at any stage in the early Universe. If it appears during BBN era, then the abundances of light elements will be affected due to the increment of the Hubble expansion rate. In this paper, we analyze the BBN consequences of the stepwise scalar field proposed in [30]. We provide here a public Matlab program, BBNLab, which enables the BBN calculation in the stepwise scalar field model. We find that BBN provides a strong upper limit $(\sim 4\%)$ to the nearly constant $\Omega_{\phi}$ (see Fig. 1). However, in the stepwise scalar field model, $\Omega_{\phi}$ can vary over several orders of magnitude, and can remain less than $4\%$ over a wide period of time. Based on this, we show examples where the stepwise scalar field can survive the BBN constraints. The scalar field can be completely hidden in BBN era, which restores the standard BBN results (see the case of $N_{\rm ini}=70.0$ in Fig. 3 or $N_{\rm ini}=78.2$ in Fig. 4). A more interesting example is that the very early dark energy that appears near the end of BBN can be used to solve the possible deuterium problem [39] (see the case of $N_{\rm ini}=70.2$ in Fig. 3).

In the stepwise scalar field model, the evolution of the Universe can be divided into two categories: OSS and CAS [31]. In order to facilitate the cosmological parameter constraints, we argued that the Universe should evolve as OSS in the radiation era [31]. This is a subjective argument, not a result of observational constraints. However, the above $4\%$ upper limit, together with the constraints from the cosmic late-time acceleration, rule out the possibility of OSS in the radiation era. The examples provided in this paper that not only survive the BBN constraints, but also explain the cosmic late-time acceleration, are all belong to the category of CAS. Therefore, we conclude that the Universe should evolve as a CAS in the radiation era in the stepwise scalar field model. This does not rule out the model, but brings difficulties to the global cosmological parameter constraints: there may be many local maxima in the posterior distribution [30]. How to quantitatively describe such a posterior distribution is a major task for our future work.

What kinds of observations could verify or rule out the CAS scenario for the radiation era? Given the above discussed difficulty of cosmological parameter constraints, it is important to look for direct observations. Electromagnetic observations seem impossible due to the scattering of photons and free electrons. However, gravitational waves and dark matter observations might be possible. The stochastic gravitational wave background can be used to infer the cosmic expansion history before BBN [67, 68, 69]. The possible axion detection provides an auxiliary probe of the early Universe [70]. Future multimessenger astronomy can directly probe the early Universe and test the CAS scenario in our model.

Acknowledgements

This work was supported by the Initiative Postdocs Supporting Program under Grant No. BX20200065 and China Postdoctoral Science Foundation under Grant No. 2021M700481.

Appendix A Special functions

Here we summarize the special functions involved in the BBN calculation:

[TABLE]

where the second kind modified Bessel function [71]

[TABLE]

In Matlab, $K_{\nu}(x)$ is calculated by $\texttt{besselk}(\nu,x)$.