Big Bang Nucleosynthesis Hunts Chameleon Dark Matter

We study the chameleon field dark matter, dubbed scalaron, in F (R) gravity in the Big Bang Nucleosynthesis (BBN) epoch. With an R2-correction term required to solve the singularity problem for F (R) gravity, we first find that the scalaron dynamics is governed by the R2 term and the chameleon mechanism in the early universe, which makes the scalaron physics model-independent regarding the low-energy scale modification. In viable F (R) dark energy models including the R2 correction, our analysis suggests the scalaron universally evolves in a way with a bouncing oscillation irrespective of the low-energy modification for the late-time cosmic acceleration. Consequently, we find a universal bound on the scalaron mass in the BBN epoch, to be reflected on the constraint for the coupling strength of the R2 term, which turns out to be more stringent than the one coming from the fifth force experiments. It is then shown that the scalaron naturally develops a small enough fluctuation in the BBN epoch, hence can avoid the current BBN constraint placed by the latest Planck 2018 data, and can also have a large enough sensitivity to be hunted by the BBN, with more accurate measurements for light element abundances as well as the baryon number density fraction. ∗ huachen@mails.ccnu.edu.cn † taishi@mail.ccnu.edu.cn ‡ synya@jlu.edu.cn § qiutt@mail.ccnu.edu.cn 1 ar X iv :1 90 8. 04 14 6v 1 [ he pph ] 1 2 A ug 2 01 9


I. INTRODUCTION
While the observations [1,2], again and again, suggest the existence of dark energy and dark matter, their nature remains unknown: what are they and how they arise? Such a dark side of the universe strongly indicates extensions of general relativity (GR) or standard model (SM) of particle physics. In cosmology, the ΛCDM model fits the observation very well, although it relies on the fine-tuning of the cosmological constant. As an alternative, a large number of works have attempted to modify Einstein's gravity to dynamically explain the dark energy problem (see [3] for a thorough review). Among them, F (R) gravity is one of the straightforward extensions for Einstein's gravity.
In order to be realistic, the modified gravities have to recover the GR at small scales or dense regions. For example, F (R) gravity fails to meet such a requirement without the chameleon mechanism [4]. With the chameleon mechanism, F (R) gravity introduces an additional field degree of freedom in Einstein frame, which could serve as dynamical dark energy. There are several viable models applied to the late-time cosmology; for instance, Hu-Sawicki model [5], Starobinsky model [6], Tsujikawa model [7], and exponential model [8].
However, as was revealed in [9,10], dark energy models in F (R) gravity generally suffer from the so-called singularity problems. Thus, we need further modifications to cure the singularity problems, and a resolution by adding a second-order correction αR 2 was proposed in [11,12].
As for the dark matter, although it is possible to regard the dark matter as a pure gravitational effect in the modified gravity [13], the cosmic microwave background (CMB) data fits well with the ΛCDM model, and also the evidence of rotation curve [14] and weak gravitational lensing [15] favor the presence of dark matter as particles. Therefore, it is widely believed that the existence of dark matter is associated with extensions of the SM which could also have some correlation with other mysteries in the universe; for instance, the axion-dark matter that is responsible for the violation of charge and parity (CP) symmetry in the quantum chromodynamics (see [16][17][18], for a recent review on the cosmological impact of axion-dark matter).
On the contrary to the folklore, it has been suggested that dark matters may also arise from a gravitational origin [19], for instance through an R 2 gravity [20]. Inspired by those works, two of authors recently proposed a scenario to make both the dark energy and dark matter played by a chameleon field, dubbed a scalaron, in the so-called Starobinsky model of F (R) gravity [21,22]. There the dark side originates from the chameleon field, whose potential minimum corresponds to the dark energy and the excitation mode around the potential minimum behaves like the dark matter. As a consequence, we can address the coincidence problem [21]. A similar scenario has also been applied to the logarithmic F (R) model recently [23].
The scalaron can have other significant characteristics in the sense of astrophysical and phenomenological consequences: The scalaron mass becomes massive in the dense region while gets light in the cosmic background. The scalaron can thus avoid being detected in the current dark matter searches but affects the cosmological (possibly, astrophysical) phenomena as the dark matter, which is compared to other scalar-dark matter scenarios, such as axion(-like) particles [16], light scalar dark matter [24,25], and ultra-light scalar dark matter [26]. In further comparison with some dark matter models which also predict a variable mass [27,28], the scalaron could provide a characteristic signal.
In this paper, we first show that with an R 2 -correction term required to solve the singularity problem for F (R) gravity, the scalaron dynamics is governed by the R 2 term and the chameleon mechanism in the early universe. It makes the scalaron physics modelindependent regarding the low-energy scale modification in the early universe. In viable F (R) dark energy models including the R 2 correction, our analysis suggests the scalaron universally evolves in a way with a bouncing oscillation irrespective of the low-energy modification for the late-time cosmic acceleration.
In relation to the dark matter physics in the early universe, the Big Bang Nucleosynthesis (BBN) would give us the stringent observational constraint on the new particles beyond the SM of particle physics. Therefore, one may naively expect that the BBN constraint should be applied to the scalaron dark matter, which, in contrast to the other dark matter candidate, has the environment-dependence reflecting the chameleon mechanism.
We find a generic bound on the scalaron mass in the BBN epoch, to feedback to the constraint on the coupling strength for the R 2 term, which is universally applicable to the viable modified gravity. It turns out that this bound is more stringent than the one placed by the fifth force experiments. We then demonstrate that the scalaron naturally develops a small enough fluctuation in the BBN epoch, hence can avoid the current BBN constraint placed by the latest Planck 2018 data. The size of the detection sensitivity depends on the initial condition for the scalaron fluctuation in the early universe, so the BBN data with more accurate measurements for light element abundances as well as the baryon number density fraction can exclude, or probe some parts of the scalaron scenario.
Throughout the main body of the present paper, we will work on the flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric, diag (−1, a 2 , a 2 , a 2 ), while in the appended discussion given in Appendix A, we will take the Minkowski metric with diag (−1, 1, 1, 1) for convenience.

II. DARK ENERGY AND DARK MATTER IN F (R) GRAVITY
This section provides the brief introduction of the F (R) gravity motivated by the latetime cosmic acceleration. F (R) models for the dark energy, in general, suffer from the singularity problem arising at a finite value of the field, which can be improved with the higher curvature term, and we introduce the model of our interest in this paper. We also mention the concept and advantage to unify the dark energy and dark matter in F (R) gravity.

A. Dynamical Dark Energy in F (R) Gravity
In this subsection, we review the F (R) gravity and related dark energy models. The action of F (R) gravity reads Here we can recover the Einstein's gravity with an additional scalar field ϕ, so-called the scalaron, in the Einstein frame: where the scalaron potential V (ϕ) is defined as The variation of the action Eq. (3) with respect to g µν gives whereT The variation of the action Eq. (3) with respect to the metricg µν gives where˜ =∇ µ∇ µ , and∇ µ denotes the covariant derivative composed of the Einstein-frame metricg µν . Then we define the effective potential, Note that the trace of the energy momentum tensor T µ µ includes a nontrivial dependence of ϕ via the Weyl transformation of the metric in the matter Lagrangian. When one ignores such a dependence, Eq. (8) is reduced to be which has been used in the earlier works [29,30] with a fluid approximation of the matter sector.
The square of the scalaron mass is computed as the second derivative of the effective potential at the potential minimum ϕ = ϕ min : where ϕ min satisfies V eff ,ϕ (ϕ min ) = 0. Due to the exponential coupling between ϕ and T µ µ , the effective potential and the mass are environment-dependent quantities. Then we can find that there is a positive correlation between the scalaron mass and energy density: the scalaron mass becomes heavy in the high-density region, for instance in the solar system, while it becomes light in the low-energy region, for instance in the cosmic background. Such an environment dependence is called the chameleon mechanism, which enables F (R) gravity to satisfy the observational constraints.
The additional field degree of freedom has been broadly studied as the source of dark energy. However, it has been shown that dark energy models in F (R) gravity may suffer from the curvature singularity problem [9]. To review it in detail, we consider the following dark energy models, 1. Hu-Sawicki model [5]: 2. Starobinsky model [6]: 3. Tsujikawa model [7]: 4. Exponential model [8]: where R c ∼ Λ 4 × 10 −84 [GeV 2 ] corresponds to the dark energy scale, and c i , n, and β are model parameters. It is easy to find that R → ∞ when κϕ → 0 in Eq. (2). Looking at the Starobinsky model as an example (Fig. 1), we observe that the original potential is lifted due to the matter contribution. Eventually, the minimum of the effective potential gets closer to κϕ = 0, and the scalaron field can smoothly reach κϕ = 0 when one considers the perturbation around the minimum, because of the quite shallow potential form. This implies the curvature singularity at R = ∞ appears at a finite field value and energy level. To cure the curvature singularity problem above, one could extend the models listed in Eqs. (11), (12), (13), and (14) by adding a second-order correction, αR 2 [11,12]. Note that those models can explain the dark energy problem in the late-time universe, although there are some differences quantitatively. As we will see later, however, if we add a secondorder correction to cure the singularity problem, the differences among those models become negligible in the large-curvature region, and so the models would be indistinguishable in the early universe when the R 2 term dominates. So, we shall momentarily work only on the Starobinsky model of our interest.
We should note that R 2 correction also plays a significant role in the finite-time future singularity problem [10]. The earlier works have shown that R 2 term can remove several types of the above singularities although the so-called type IV singularity remains (see Refs. [31,32] for the classification of these singularities). As was suggested in [33], the finitetime future singularity problem has the same origin as the curvature singularity problem stems from, and thus the curvature singularity can be rephrased as one of the finite-time future singularities. Furthermore, one can find that not only R 2 term but also R n , where 1 < n < 2, can work to remove those singularities from F (R) gravity.

B. Chameleon Field as Dark Matter in F (R) Gravity
As an illustration, we consider the Starobinsky model with R 2 correction.
where α 10 22 [GeV −2 ] to be consistent with the constraint from the fifth force experiment [34]. Note that this R 2 does not have to be related to an inflaton dynamics in the current analysis. Since the R 2 term will make the curvature sharp as the fluctuating ϕ goes to zero, the ϕ cannot perturbatively access the singular point; hence, the singularity problem is to be gone if the potential value at R → ∞ is high enough. Now, we can safely consider the fluctuation of the scalaron field around the potential minimum. If we quantize such an excitation mode, we obtain the particle picture, which literally brings us a new particle beyond the SM of the particle physics. Two of the authors have proposed that this new particle would be a dark matter candidate [22]. With this scenario at our hand, we may resolve the dark matter problem in the framework of the modified gravity in addition to the dark energy problem, which not only gives a unification of the dark side of the universe but also allows us to make the coincidence problem understandable [21].
Since the R 2 term enable us to link the original Starobinsky model to the early universe, it is natural to speculate that dark matter is produced in the early universe. In relation to the dark matter physics in the early universe, the BBN gives us the stringent observational constraint on the new particles beyond the SM of particle physics. Therefore, one may naively expect that the BBN constraints should be applied to the scalaron dark matter.
In contrast to the other dark matter candidate, the scalaron and, of course, observational constraints on it have the environment-dependence, reflecting the chameleon mechanism. In order to address such a chameleon dark matter, it is inevitable to introduce an appropriate external matter. In the following, we shall look into more details on the environmentdependence in the early universe and the BBN constraints on the scalaron.

III. CHAMELEON DARK MATTER IN EARLY UNIVERSE
Cosmological studies for the early universe have been investigated with focus on the scalaron dynamics and the chameleon mechanism [34] (also see [29] in the scalar-tensor theory). In the following analysis, we study the time evolution of the scalaron field in the early universe during the epochs with T ∼ 100 [GeV] (after the electroweak phase transition) down to ∼ 1 [MeV] (around the BBN), by explicitly solving the equation of motion for the scalaron. We will adopt a toy-model to evaluate the chameleon mechanism to get a rough prospect for the chameleon in the overall environmental effect through the cosmic history.

A. Analytical Approximation
To help our analytic understanding of the scalaron potential structure, by using the fluid approximation, we shall model the environment surrounding the scalaron overall in the target-early universe. In this case, the trace of energy-momentum tensor, which shows up in Eqs. (9) and (10), is given by where ρ and p denote the energy density and pressure of the matters, respectively. For further convenience, we define where ξ(T ) = (ρ − 3p)/ρ and ξ(T ) is dimensionless function of temperature T . We consider the FLRW background, which suffices to study the scalaron evolution in the early universe.
Then, the (0, 0)-component of Eq. (6) gives Substituting Eq. (18) into the equation of motion, Eq. (7), we obtain the time evolution equation for the scalaron, Eq. (19) is a highly nonlinear differential equation, which is practically so hard to solve.
To grab a rough insight into the scalaron evolution, in this subsection we shall first make several approximations to the model, Eq. (15), so that Eq. (19) is written in an analytic and manageable form, which can include the essential part in the early universe. We will perform the numerical analysis and discuss the validity of those approximations in the next subsection.
For convenience, we define r ≡ R/R c and γ ≡ e 2 Eqs. (2) and (15), we have where the second line can be achieved if r 1 as in our targeted-early universe. When where consistently with the constraint from the fifth force experiment [34].
As r is large for our targeted-early universe, we can safely neglect the second term in Eq. (15), which corresponds to the present dark energy scale, so that the model Eq. (15) is approximated as Here the second line is achieved by using the Weyl transformation Eq. (2). It is crucial to note that the approximated expression above is also valid for other models quoted in Eqs. (11), (13) and (14) with the R 2 correction, because it is free from the present dark energy scale.
Thus we obtain analytical forms of potential terms as follows: V ,ϕϕ (ϕ) = 1 6α The stationary condition for the effective potential V eff ,ϕ (ϕ min ) = 0 reads where γ min = e 2 √ 1/6κϕ min − 1. The second derivative of the effective potential is written in terms of γ as We derive the analytically approximated formula for the scalaron mass m(ϕ min ) in Eq. (10) (28)

B. Numerical Evaluation for Time Evolution of Scalaron Field
For the numerical analysis, we should carefully check whether the above approximations can work or not. Firstly, we shall observe the form of the effective potential in Eq. (9), as depicted in Fig. 2, where we choose −T µ µ ∼ 10 −18 [GeV 4 ] and ρ Λ ∼ 10 −48 [GeV 4 ] following the previous work [21]. One could find that in the negative region, κϕ < 0, the matter contribution with −T µ µ (BBN) ∼ 10 30 ρ Λ is highly dominated. It allows us to safely set the first term to be zero in Eq. (9), so that only the matter contribution can be significant in the negative region.
Remarkably enough, the potential minimum emerges at around the origin ϕ = 0, under the matter contributions in the BBN epoch. It happens because of the balance created by the R 2 term and nonzero matter contributions, i.e., the chameleon mechanism. Thus this fact gives us great support to study the scalaron particle physics in the BBN epoch, as will be described later in details.
Secondly, we compare the approximated potential in Eq. (23) with the original potential in Eq. (4) (without matter contribution). See Fig. 3. We find that those two potential forms fit well in the positive region, while there is a large gap in the negative region. So we should set the approximated potential to be zero in the negative region. An easy way is to multiply it with Heaviside step function Θ(ϕ), V ,ϕ (ϕ) ≈ 1 6κα where we have Taylor expanded the exponential part for later convenience. We can find that the minimum locates at around κϕ ∼ 10 −16 in the BBN epoch as shown in Fig. 4. Note that the mass formula (28) holds because the potential minimum always locates in the region Now that we have confirmed reliably approximated analytical expressions for the effective potential (Eqs. (23) and (30)), we shall substitute them into Eq. (19) to study the scalaron evolution. By assuming the radiation dominance, H = 1 2t , the equation of motion Eq. (7) thus goes like where φ = κϕ and τ = t/ √ α, both of which are dimensionless.
Because of the successive decouplings of SM species during the time evolution in the early universe (scaling T ∼ 100[GeV] down to ∼ 1[MeV]), ξ(T ) is expected to roughly oscillate between 0.01 and 0.1 [21,30,35], as depicted in Fig. 5 which has been drawn by using the same evaluation of the ξ(T ) as in the literature [21]. Though the ξ depends on the time in a nontrivial way, to get a rough insight, we shall take it to be constant in time. According to the form of the time evolution Eq. (31) with the vicinity of the potential minimum (φ ∼ 0) focused on, this brute ansatz could be operative as long as the time scale is large enough compared to the contribution from the αR 2 term, i.e. τ 1. oscillate with damping for a larger time scale. Of importance is to notice that this is the charactereistic feature universally realized in a class of F (R) gravity, as listed in Eqs. (11) - (13).
This behavior can be understood by looking back Fig 2: the scalaron falls from a sharp and high potential barrier when it moves from the positive region, then turns to the negative region, where the matter contribution becomes dominant. Thus the scalaron decelerates till stop due to the matter contribution T µ µ and Hubble friction H(t). Then it will go back and will be blocked by the potential barrier, which is built very close to the origin, that is how it begins to bounce from the potential minimum.
The eminent turnaround takes place near τ ∼ 10 2 (i.e. t ∼ 10 −11 [s]), which correspond to a far past time even compared to the (expected) electroweak phase transition epoch around where the matter contribution can be highly model-dependent in a sense of possible models beyond the SM of particle physics. In the present study, we have simply assumed the order of ξ does not alter, and the ξ is not going to be zero. Note that the short periodicity and small amplitude of the oscillating scalaron field around the potential minimum at κϕ min ∼ 0 allows Eq. (18) to be almost reduced to the ordinary Friedman equation. The above result is consistent with our assumption H = 1 2t in solving the equation of motion Eq. (31). As we will see in the next section, the scalaron field affects observables in the early universe through the Hubble parameter. Hereafter, we shall assume the radiation-dominant Hubble parameter to be evaluated by the relativistic species of the SM particles. First of all, it should be noted that the scalaron has to be non-relativistic in the BBN epoch, otherwise the successful nucleation processes would be spoiled (for a review on light particle contributions to the BBN, see the review part in Particle Data Group [36]). It gives a universal bound for the modified gravity scenario to survive, including the class of models listed in Eqs. (11) - (14). As has been demonstrated in [21], the scalaron mass evolves as an almost constant in time (temperature) around the BBN where our approximated model Eq. (22) works. In fact, we can easily check it from Eq. (28) when α is small enough to satisfy the fifth-force constraints: Since the BBN took place at around T ∼ 1 [MeV], for the scalaron to be non-relativistic, the parameter α thus gets constrained to be which is more stringent than the bound that the fifth force has currently placed.

A. Chameleon Dark Matter in BBN Epoch
A typical time scale for the BBN (t ∼ 10 2 [s]) corresponds to τ ∼ 10 15 in our setup. Thus, we can study the scalaron dynamics, which connects to the results in the previous section, by taking the large τ limit (with φ 1) in Eq. (31). In such a limit, one can approximate Eq. (31) as follows:φ for φ > 0 withφ(τ ) 0 assumed. Eq. (34) implies that the matter contribution gives a perturbative effect on the harmonic oscillation of the scalaron as a damping factor (hence it is bounced), as long as the ξ is small enough (like ξ ∼ 0.01 − 0.1), but nonzero. Thus, we can deduce that in the BBN epoch, the scalaron is present as a stable field at the vacuum close to the origin of the potential (φ = κϕ ∼ O(10 −16 )) and keeps undergoing a slow damping oscillation, as seen in the case with orange dot-dashed line in Fig. 6 Note that the amplitude of the scalaron oscillation depends on the initial condition of the scalaron field. We have assumed that the scalaron is staying close to the potential minimum in the cosmic history. The above is a natural assumption and allows the scalaron to roll down the potential in the earlier universe if we incorporate the all cosmic history of the scalaron field. Moreover, if the amplitude of the oscillation is small enough in the earlier universe, as shown in Fig. (6), it is natural to expect following Eq. (34) that the amplitude of the oscillation should be smaller than the order of unity in the BBN epoch.
A similar argument was done in Ref. [29]. In the case that the initial value of φ, φ i , locates in the region φ < 0, the potential term induced by the R 2 term is negligible, and the matter contribution T µ µ and the Hubble friction 3Hφ dominate in the equation of motion, which forces the scalaron to get closer to the potential minimum. On the other hand, in the case that the φ i locates in the region φ > 0, the original potential term including the R 2 term dominates, and it makes the scalaron rolled down to the potential minimum. Those two effects may justify our assumption that the scalaron locates around the potential minimum in the early universe, which also does not allow us to set the unnaturally large initial value of the scalaron.
Keeping such an oscillation behavior of the scalaron in our mind, we shall next discuss the frame-difference between the Jordan and Einstein frame. Because the matter sector is introduced in the Jordan frame, all dimension-full observables receive the effect of the Weyl transformation. According to the oscillating solution of the scalaron as in Eq. (34), we can naturally expand the scalaron field as where ϕ min and δϕ denote the field value at the potential minimum and the fluctuation of the scalaron field. Then the exponential form of the scalaron field in the Weyl transformation can be rewritten as where C is an arbitrary coefficient.
The above equation allows having the perturbative picture for the scalaron couplings to the matter fields, and the overall factor e Cκϕ min is corresponding to the frame-difference associated with the Weyl transformation of the metric. We now note that the frame-difference is negligible because κϕ min is almost equal to zero. Consequently, there is almost no framedifference in the observables such as the S-matrix elements including the mass spectra.

B. Scalaron Couplings to BBN
Now we are ready to discuss the scalaron couplings to the BBN. The BBN constraints mainly arise from the precise measurement on the helium abundance (Y4 He ). When there exists a nontrivial field configuration not involved in the standard BBN model, the Y4 He can generically get modified through the corrections to the nucleon (neutron-proton) mass difference ∆m N ≡ m n − m p and the neutron lifetime τ n ∼ [G 2 F m 5 e ] −1 , with the electron mass m e and the Fermi constant G F (and possibly also to the effective degrees of freedom for relativistic particles at that time) because these are related to models of particle physics on which one works. The issue is, however, a bit complicated and somewhat nontrivial for the scalaron couplings to the nucleon, because the nucleon is a composite particle of quarks as a consequence of quantum chromodynamics, QCD. For preparation to study the interactions between the scalaron and such a composite particle, we begin by reviewing the scalaron couplings to the elementary fermions, such as quarks and leptons, following Ref. [22].
Let us first consider the following Lagrangian of the massive fermion field ψ: where m F denotes the fermion mass, and the generalized Dirac gamma matrix in the curved space-time satisfies {γ µ , γ ν } = −2g µν . After the Weyl transformation of the metric g µν = µν , the action of the fermion field is transformed as where the Dirac gamma matrix in the Einstein frame is transformed asγ µ = e − √ 1/6κϕ γ µ for the consistency with the Weyl transformation of the metric.
Here, we can redefine the field ψ to realize the canonical kinetic term: For the above field redefinition, the action in Eq. (39) is written in terms of ψ : This reads the minimal coupling form between the scalaron and massive Dirac fermion fields.
In a similar way, we may derive the scalaron coupling to a four-fermion interaction (with the coupling constant G 4f ) which would be relevant to the neutron decay process.
Thus, the interaction terms between the scalaron and massive (elementary) fermions arise like dilatonic scalar couplings reflecting the scale dimension for the corresponding interaction operators. Note also that the scalaron coupling does not distinguish any charges and flavor of the leptons and quarks, which gives us a naive expectation that the similar coupling forms show up for composite particles, such as the baryons and mesons after the Weyl transformation.
However, this is not the case for the nucleon: since the nucleon dynamics is highly nonperturbative in the sense of conventional perturbation theories, higher-dimensional nucleon interactions cannot simply be neglected. Among those non-minimal interaction terms, for instance, four-nucleon interaction terms of Lorentz vector(or scalar) type might be induced in QCD, which can arise along with the isospin breaking, such as (N γ µ σ 3 N ) 2 with N and σ 3 being a nucleon doublet (p, n) T and the third component of the Pauli matrices, respectively.
Generically, such interactions can generate nontrivial mass-split for the nucleon at quantum loop levels for the nucleon dynamics, and can also couple to the scalaron in a nontrivial way, due to the Weyl transformation and rescaling of the nucleon field, as demonstrated above. Thus, one naively suspects that the nonperturbative nucleon dynamics would drastically make the scalaron coupling to the nucleon mass modified from the minimal and the simplest form as in Eq. (41). What is worse, one cannot naively suppress all higher-order terms beyond the four-fermion type, because of the nonperturbativity, so that the form of the scalaron coupling to the nucleon mass would be out of control at the quantum level.
Note that this complexity is characteristic to the generic argument on the scalaron coupling to the nucleon, which does not come in the case of axion-like particle not associated with the Weyl transformation.
A way out of this complex issue can be hinted by noticing that the nucleon mass difference signals the explicit breaking of the chiral (isospin) symmetry for the nucleon.
Since the observed size of the nucleon mass difference is quite small compared to its mass, It is possible to demonstrate the above expansion based on the framework of the so-called chiral perturbation theory, where the multi-body nucleon interactions including the isospin breaking can be treated as perturbations compared to the leading-order correction set by the tree-level contribution from the quark mass difference to the nucleon mass.
As explicitly given in Appendix A, the scalaron coupling to the nucleon mass difference ∆m N can thus be evaluated based on the baryon chiral perturbation theory (BChPT). This BChPT is formulated in the Einstein frame but can be converted to the theory in the Jordan frame by taking into account a negligibly small discriminant exponential factor, as argued in the above. It will be ignorable compared to other uncertainties, e.g., in the BChPT. It turns out that at the nontrivial leading order (LO) of the perturbation, the scalaron affects the mass difference just like in a dilatonic coupling form, which coincides with the one read off from the minimal Dirac fermion-scalaron action in Eq. (41): Beyond the leading order corrections, the nucleon-mass splitting could get affected by interactions with pions and also by electromagnetic (EM) interactions with photon when the chiral symmetry is (in part) gauged. As in the literature [37], among those next LO corrections, the EM corrections will be most significant in magnitude. Indeed, the total nucleon -proton mass difference is shared by the EM correction part with a "pure QCD" part (the up-and down-quark mass difference) like where the pure QCD part has been extracted just by subtracting the observed value [36] by the EM term estimated assuming a one-photon exchange contribution based on the current algebra technique [37]. Eq. (44) implies that the EM correction pulls the mass difference down by about 37%, so it is indeed somewhat a significant destructive interference, which implies that the scalaron coupling to the EW correction part needs to be taken into account as well.
In Eq. (44), the ∆m QCD N is identified as the ∆m N in Eq. (43). Regarding the EM correction, the scalaron can actually modify it through the EM scale anomaly ∼ e 2 /(4π) 2 F µν F µν , as discussed in [22]. The scalaron thus makes the EM coupling e modified as e 2 → e 2 (ϕ) = e − √ 1/6κϕ · e 2 [22] along with the beta function coefficient arising from the charged fermion loops. Note that the scaling factor for the EM correction part is the same as that for the quark mass difference (∆m N ) part. Thus, the scalaron just gives an overall scaling for the neutron -proton mass difference in total, in such a way that We shall later on use this Eq. (45) to discuss a rough size of effects on the BBN.
As to the nucleon lifetime τ n ∼ [G 2 F m 5 e ] −1 which is saturated by the beta decay n → p + e +ν, it is not due to the strong QCD, but the weak interaction process. Hence we do not need to apply the BChPT but scale it by the scalaron field factor of e −κϕ/ √ 6 with respect to the scale dimension of minus one, following the four-fermion interaction and elementary fermion mass terms in Eqs. (41) and (42). Thus we have τ n (ϕ) = τ n · e +κϕ/ √ 6 . (46) Now that we have derived the relevant formulae (Eqs. (45) and (46)) for the scalaron effect on the BBN, in the next subsection we shall discuss the BBN constraint on the scalaron and the size of its sensitivity -possibility to hunt the scalaron -in the near future.

C. BBN Constraints
We begin by looking at the generic formula for the helium abundance at the BBN epoch Y4 He (t BBN ) [38]: in which for later concrete estimates, we will take t BBN 180[s], corresponding to 1 [MeV].
The overall yield Y4 He (t → 0) is given by integrating the balance equation involving a couple of n → p conversion processes (nν e → pe − , ne + → pν e , n → pe −ν e ) [38]: where y = ∆m N /T , and where g * (t BBN ) denotes the effective degrees of freedom for relativistic particles in the BBN era. Here we have assumed that the scalaron does not significantly affect the Hubble parameter H in Eq. (18) in the BBN epoch, and the H can be evaluated by the fluid picture for the relativistic particles at that time t BBN . The former assumption would be reasonable: in Fig. 6 with a small enough matter contribution (ξ O(10 −2 )), we have observed that the approximated time evolution Eq. (34) describes the damping oscillation very well for a large 24 . It implies that we can numerically neglect the implicit scalaron dependence in the Hubble parameter H as in Eq. (18), and replace the H as the standard BBN's.
As to the effective degrees of freedom g * (t BBN ), we will take g * (t BBN ) 3.36 that is a standard value after the e + e − annihilation [38], because the scalaron acts like a nonrelativistic particle to make no contribution to the background Hubble parameter, including the above argument on its numerically negligible time evolution effect.
Though the scalaron may not significantly affect the background Hubble evolution as argued above, the scalaron can be hunted by the BBN: Since the scalaron fluctuation in the BBN epoch can be controlled to be small enough, due to the presence of the balance between the R 2 term and the chameleon mechanism, the size of the deviations from the SM predictions in Eqs. (45) and (46) can be small as well. In that case, the helium abundance will get a small shift to be modified from the one in Eq. (47) by the small amount. Then we can easily see the correction arises in a perturbative series like with the higher order terms neglected. (Here we have used τ n 880[s] from Particle Data Group [36].) According to the current status on the standard BBN prediction estimated using the latest Planck 2018 data for the baryon number density fraction, we may quote Y4 He (t BBN )| SM 0.247 (at 95% C.L.) [39]. Taking the 1 sigma range for the observed helium abundance, Y4 He (t BBN )| obs = 0.243 +0.023 −0.024 (including gravitational lensing and baryon acoustic oscillation data and taking the best-fit number for the neutrino species, N eff = 3.046 at present) [39], we thus get the BBN bound on the scalaron field, −0.05 < (κϕ) < 0.03 .
In the previous section, we showed that thanks to the R 2 term and the chameleon mechanism, the scalaron always stays around the minimum of the effective potential, with the size of the amplitude |κϕ| = O(10 −3 − 10 −2 ), in the BBN epoch (see Figs. 2 and 6). We again note that the amplitude of the scalaron depends on the initial condition. Considering the damping behavior of the scalaron field (as seen in Fig. 6), we see that the scalaron fluctuation can naturally be small enough to survive the BBN constraint in Eq. (51), unless we impose the unnatural initial condition for the scalaron field. Note the current observational accuracy on the helium abundance is roughly on the order of the allowed size for the oscillation amplitude for the scalaron |κϕ| 10%. Therefore, future updated measurements on the light element abundance as well as the baryon number density fraction would exclude a class of scalaron scenarios with somewhat a larger initial value (∼ 0.1 − 1%) as in Fig. 6

V. CONCLUSION
Dark energy and dark matter are two of the biggest mysteries in cosmology and particle physics. Although they are usually studied separately, it is still possible to postulate that they have the same origin if one tries to resolve the phenomenological coincidence problem.
In recent years, there are some attempts to incorporate the dark matter into the framework of modified gravity besides dark energy. In this paper, we have studied the unification of dark sector in the framework of F (R) gravity and investigate the evolution of chameleon dark matter in the BBN epoch.
We have introduced the second-order correction R 2 to cure the singularity problems which show up in generic models of the F (R) dark energy. We have confirmed that in the early universe, the model dependence of modifications for dark energy disappears, which is rea-sonable because the infrared dark energy has little effect on higher energy physics. Thus we have demonstrated that in the early universe, the model in F (R) gravity which corresponds to the dark sector could be approximated as R + αR 2 , to find that the scalaron universally evolves in a way with bouncing oscillation, irrespective of the low-energy modification for the late-time cosmic acceleration (Fig. 6). Such a universal feature shows up due to the R 2 term and the chameleon mechanism characteristic to the scalaron.
We At last, we will make several comments for future works. In our analysis, we have approximated ξ(T ) to be constant in time (ξ = 0.01, 0.05, 0.1 as in Fig. 6), which has given sufficient information about the scalaron field in the BBN, as far as the typical order of magnitude is concerned. Although the temperature dependence of ξ(T ) has computed in [21] and also in Fig. 5 of the present paper, we still have theoretical ambiguities for the earlier or later cosmic history; that is, we can include the various effects from models beyond the SM at the higher energy scale, and the decouplings and light elements at lower energy scales. Two of the authors have embedded the scalaron dynamics into the scale-invariant extension of the SM and shown the vanishing T µ µ at the electroweak phase transition. Thus, T µ µ goes from zero to nonzero, which would be the first kick to make the scalaron close to the potential minimum. In this case, the argument about the initial condition should be examined to relate the scalaron physics with the possible (classical) scale-invariance in the very early universe.
Moreover, towards the complete analysis of the scalaron physics in the BBN epoch, we have to address a nontrivially interacting system which comprises the scalaron and the other matters, instead of merely substituting the cosmic history of the other matters as external fields into the scalaron potential term. It is mandatory to reveal the scalaron couplings to the composite particles, that is, nucleons and light elements in the early universe, in which the BChPT approach we have employed in this paper shall play crucial roles. Notably, such research also gives us potential tools to study theoretical predictions for the DM direct search experiments in the current universe. Various experiments are designed to detect the dark matter candidate through the scattering with nucleus or electron in the atom. The detector in such experiments consists of, for instance, liquid xenon or argon, and the scalaron collides with the atomic nucleus. Thus, as a branch of this work, the study on the scalaron coupling to the composite particles has a potential impact for the future works to reveal the scalaron dynamics in both early and current universe.

ACKNOWLEDGMENTS
We are grateful to Masato Yamanaka for useful discussions related to the BBN bounds. To formulate the mass difference coupled to the scalaron in the Einstien frame, in this Appendix, we introduce the techniques in the hadron physics, which allows one to incorporate isospin breaking effects in a systematic way, called the baryon chiral perturbation theory (BChPT).
Thus the chiral Lagrangian invariant under the G-symmetry can be written in terms of those 1-forms, and at the leading order of derivative terms (of O(p 2 )) we find where at this moment the kinetic terms of the pion fields have been canonically normalized.
To incorporate the explicit-chiral symmetry breaking effect reflecting the current quark mass terms in the underlying QCD (arising from the electroweak symmetry breaking via the Higgs), one may introduce a spurion field χ, which makes the current quark mass term chiral-invariant form, likeq L χq R +q R χ † q L , by allowing the χ to transform in the same way as the chiral field U does. The χ is then assumed to develop the "vacuum expectation value" χ = 2B 0 M with M = diag{m u , m d }, so that the chiral symmetry is explicitly broken consistently with the underlying QCD. Here the overall coefficient B 0 will be related to the chiral condensate later. Then one may have the following term, corresponding to the explicit breaking by the current quark mass, in a chiral invariant way: where in the second equality theχ readsχ = ξ L ·χ·ξ † R , which transforms likeχ → h(π, g L , g R )· χ · h † (π, g L , g R ). Expanding the L χ in powers of the π fields, one can readily find the pion mass, m 2 π = B 0 2 (m u + m d ), which can be compared to the Gell-Mann-Oakes-Renner relation based on the current algebra (m 2 π f 2 π = − qq (m u + m d )) to derive B 0 = −2 qq /f 2 π . The chiral perturbation theory is thus constructed from the 1-form α ⊥ µ and the spurion fieldχ, by assuming the pion mass m π to be much smaller than a typical chiral breaking scale Λ χ ∼ 4πf π (which can be estimated by equating tree-level and one-loop contributions for pion scattering amplitudes arising from the interaction terms in L (2) ), and respecting the low-energy theorem for the chiral symmetry which leads to vanishing amplitudes involving pions when the transfer momentum p (∼ soft pion mass m π ) is set to zero. Namely, the theory is built based on the derivative expansion with the order counting rule as ∂ µ ∼ m π ∼ O(p), so that the 1-form α ⊥ µ ∼ O(p) and the spurion fieldχ ∼ O(p 2 ). Thus one realizes that the leading-order chiral Lagrangian of O(p 2 ) is given by L (2) in Eq. (A1) and L χ in Eq. (A2) respecting the chiral symmetry, and also the charge (C) and parity (P ) conjugations. The C and P transformations can be set to the building blocks, following the corresponding quark current form coupled to the underlying QCD, so that α ⊥ µ ↔ −α µ⊥ andχ ↔χ † under P , α ⊥ µ ↔ (α ⊥ µ ) T andχ ↔χ T under C. The higher derivative terms can systematically be incorporated to renormalize loop corrections arising from the lower derivative interactions order by order [40,41].

Inclusion of Nucleon in Chiral Perturbation Theory
Nucleon is heavy enough to be comparable with the chiral breaking scale ∼ 4πf π ∼ 1 GeV, and does not have the soft mass limit (never goes to zero in the chiral limit) unlike the pion, so it involves the problematic issue to incorporate the nucleon into the chiral perturbation theory, which would naively lead to the loss of the systematic order counting as presented in the meson sector. Indeed, it has been founded so far that the higher loop corrections actually bring significant cancellations in the baryon mass spectra, although the estimated masses tend to be close to the result from lattice simulations. (For the BChPT up to the next-to-next-to-next-to leading order v.s. lattice simulations, e.g. see [42], and references therein). In the present study, we shall work on the nontrivial-leading order terms in the BChPT, and try to give a rough order of estimate on the effect from the chameleon in the nucleon physics.
We introduce the chiral-nucleon fieldsN L,R dressed by the "pion clouds" of the representatives ξ L,R asN L,R = ξ L,R · N L,R , where the N L,R = (p, n) T L,R denote the undressed ones transforming under the G-symmetry like N L,R → g L,R N L,R and henceN L,R → h(π, g L , g R )N L,R .
The chiral-(and also C-and P -) invariant nucleon Lagrangian based on the nonlinear realization is at the leading order of O(p) written to be where ellipses stand for higher order terms in the nucleon field (e.g. four-nucleon terms). The explicit-breaking effect on the nucleon-mass can be incorporated by coupling the nucleon field to the spurion field χ, in a chiral-invariant way, as in the meson sector. The leading order correction thus arises at O(p 2 ) to take the form Expanding the fields in the sum of L N (1) and L m N (2) leads to the nucleon masses and its splitting at the leading order of the chiral perturbation: (A5)

Coupling Scalaron to Baryon Chiral Perturbation Theory at Leading Order
Now we make the chameleon, the scalaron field ϕ, coupled to the BChPT, which is at the leading order of derivative expansion (up to O(p 2 )) described by L (2) + L χ + L N (1) + L m N (2) in Eqs. (A1), (A2), (A3) and (A4). To this end, only we need to do is just transform the Minkowski space-time metric g µν in Jordan frame into the one in Einstein frame, g µν → g µν = e 2σ g µν by the Weyl transformation, where σ = κϕ/ √ 6 with κ being inverse of reduced Planck mass (1/M pl = 1/ √ 8πG). Then we can find coupling terms between the scalaron field ϕ and the nucleon-pion system, in the same manner as done in [22] for the standard-model sector. After canonically normalizing the nucleon field by transformingN asN → e 3/2σN , we find the action corresponding to the nucleon part in Eqs. (A3) and (A4): Note that this scale transformation gives rise to the scale-anomaly at the quantum level of the baryon chiral perturbation theory, just like in the standard-model case coupled to the scalaron discussed in [22]. This anomaly can be matched to the one induced from the scalaron couplings to the quarks and leptons in the SM, to reproduce precisely the same as the scalaron couplings to diphoton and digluon, derived in the reference.
From the action in Eq. (A6) we can easily read off the nucleon mass formulae similar to those displayed in Eq. (A5), but modified by the presence of the scalaron profile as It is interesting to note that at the leading order of the BChPT, the scalaron effects on the nucleon mass by coincidence looks just like the Brown-Rho (BR) scaling [43], or the leading-order scale symmetry relation [44,45] as a sort of the extension of the BR scaling.