Nonlinear resonant oscillation of gravitational potential induced by ultralight axion in $f(R)$ gravity

TL;DR

This paper investigates how ultralight axion dark matter causes oscillations in gravitational potential within $f(R)$ gravity, revealing that nonlinear resonance effects could enhance detection prospects via gravitational wave detectors.

Contribution

It demonstrates that nonlinear resonance between ultralight axions and scalarons in $f(R)$ gravity amplifies gravitational potential oscillations, improving detectability.

Findings

Oscillations of gravitational potential are sensitive to gravity models.

Nonlinear resonance enhances detection signals.

Potential for gravitational wave detection of axion-induced oscillations.

Abstract

We study the ultralight axion dark matter with mass around $10^{-22}$ eV in $f(R)$ gravity which might resolve the dark energy problem. In particular, we focus on the fact that the pressure of the axion field oscillating in time produces oscillations of gravitational potentials. We show that the oscillation of the gravitational potential is sensitive to the model of gravity. Remarkably, we find that the detectability of the oscillation through the gravitational wave detectors can be significantly enhanced due to the nonlinear resonance between the ultralight axion and the scalaron.

Equations

$$\varphi (t,\vec{x})={\varphi }_0(\vec{x})\cos [mt+\alpha (\vec{x})],$$

$$\rho$$

$$p$$

$$T_{\mu \nu }=\left(\begin{array}{cc}{\rho }_0& 0\\ 0& -{\rho }_0\cos (2mt){\delta }_{ij}\end{array}\right).$$

$$T=-{\rho }_0[1+3\cos (2mt)].$$

$$S=\frac{1}{2}\int d^{4}x\sqrt{-g}[R+f(R)]+S_{\text{m}},$$

$$G_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }f+(R_{\mu \nu }+g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu })f_R=T_{\mu \nu },$$

$$3\square f_R-R+R{f}_R-2f=T.$$

$$3{\ddot{f}}_R+R=-T,$$

$$g_{\mu \nu }=\left(\begin{array}{cc}-1-2\Psi & 0\\ 0& (1-2\Phi ){\delta }_{ij}\end{array}\right).$$

$$R=-6\ddot{\Phi }+2{\nabla }^2(2\Phi -\Psi ).$$

$$R=R_0+\delta R,$$

$$2{\nabla }^2{\Phi }_0={\rho }_0.$$

$$\delta R=-6\delta \ddot{\Phi }.$$

$$\delta \Phi =\frac{{\rho }_0}{8m^2}\cos (2mt)+\frac{1}{2}(f_R-⟨f_R⟩),$$

$$\delta \Phi =\delta {\Phi }_{\text{E}}\cos (2mt),$$

$$\delta {\Phi }_{\text{E}}\equiv \frac{{\rho }_0}{8m^2}=5\times 10^{-18}{\left(\frac{10^{-22}\text{eV}}{m}\right)}^2.$$

$$f=\frac{2m}{2\pi }=5\times 10^{-8}\text{Hz}\left(\frac{m}{10^{-22}\text{eV}}\right).$$

$$S=\frac{1}{2}\int d^{4}x\sqrt{-g}\,\big{[}(1+\varphi)R-2U(\varphi)\big{]}+S_{\text{m}}\ ,$$

$$U(\varphi)\equiv\frac{1}{2}\left[\varphi A(\varphi)-f\big{(}A(\varphi)\big{)}\right]\ ,$$

$$f^{\prime }(A)=\phi .$$

$$R=2U_{,\phi }=A(\phi ).$$

$$(1+\phi )G_{\mu \nu }+g_{\mu \nu }U(\phi )+\left(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu }\right)\phi =T_{\mu \nu }.$$

$$3\Box\varphi-A(\varphi)+\varphi A(\varphi)-2f\big{(}A(\varphi)\big{)}=T\ ,$$

$$V_{,\varphi}\equiv\frac{1}{3}\big{[}A(\varphi)-\varphi A(\varphi)+2f\big{(}A(\varphi)\big{)}\big{]}\ .$$

$$V_{,\phi }\simeq \frac{A(\phi )}{3}.$$

$$\ddot{\varphi}+V_{,\varphi}=\frac{\rho_{0}}{3}\big{[}1+3\cos(2mt)\big{]}\ ,$$

$$V_{\text{eff},\phi }\equiv V_{,\phi }-\frac{{\rho }_0}{3}=\frac{1}{3}(A(\phi )-{\rho }_0).$$

$$\ddot{\phi }+V_{\text{eff},\phi }={\rho }_0\cos (2mt).$$

$$\delta \Phi =\delta {\Phi }_{\text{E}}\cos (2mt)+\frac{1}{2}(\phi -⟨\phi ⟩),$$

$$f(R)=\frac{R^2}{6M^2},$$

Full text

††preprint: KOBE-COSMO-17-03

I Introduction

The cold dark matter (CDM) with a cosmological constant ($\Lambda$CDM model) is known as the standard cosmological model. So far, supersymmetric particles, the so-called neutralinos, have been regarded as the most promising candidate for the CDM. While the CDM works quite well on large scales, there exist problems on small scales. In fact, the CDM model predicts an overabundance of structure on galactic scales, which is not consistent with observations. Some common problems with CDM are the cusp-core problem 94:Flores ; 94:Moore ; 95:Burkert ; 96:Navarro ; 98:Moore ; 04:Gentile , the missing satellites problem 93:Kauffmann ; 99:Klypin ; 99:Moore , and the too-big-to-fail problem 11:Boylan-Kolchin ; 12:Boylan-Kolchin . Moreover, there was no signature of supersymmetry at the LHC. Thus, it is worth investigating another possibility, namely the axion dark matter.

Originally, the axion was invented to resolve the strong CP problem in QCD 77:Peccei ; 77:Peccei-2 ; 78:Weinberg ; 78:Wilczek . Nowadays, however, it is known that the string theory also predicts axions with a wide range of mass scales 06:Svrcek ; 10:Arvanitaki . Remarkably, the axion interacting very weakly with standard model particles is regarded as a candidate for the dark matter. In particular, the ultralight axion with mass around $10^{-22}$ eV can naturally resolve the problems on galactic scales because of its wave nature 00:Hu .111 Actually, there is a long history for the scalar dark matter. Some early works are Refs. 83:Baldeschi ; 94:Sin ; 96:Lee ; 00:Sahni . For complete references, see Ref. 16:Lee .

Indeed, according to numerical simulations of dark matter halo density profiles with the ultralight scalar field, the mass around $10^{-22}$ eV is favored by the data of dwarf spheroidal galaxies 12:Lora ; 14:Schive ; 16:Calabrese ; 16:Gonzales-Morales . For this reason, the ultralight axion has recently attracted much attention. Here, it would be fair to mention that the “small scale problems of CDM” may be resolved even in the CDM framework by taking into account astrophysical processes including baryonic matter.

Remarkably, the pressure of the ultralight axion is oscillating in time with angular frequency at twice the axion mass, $\omega=2m$. Therefore, in order to find the axion dark matter, we should detect the oscillation of gravitational potential induced by this oscillating pressure. Since the oscillation of the gravitational potential can be seen as a fluctuation of spacetime like gravitational waves, we would be able to detect the oscillation by means of gravitational-wave detectors. Indeed, Khmelnitsky and Rubakov pointed out that the effect of oscillating pressure might be detected with pulsar timing array experiments 14:Khmelnitsky . We also pointed out that laser interferometer detectors can be used for this purpose 17:Aoki . One may think interferometers have no sensitivity to the isotropic oscillation of the pressure since its two arm lengths seem to change exactly the same amount. However, the Solar System moves through the dark matter halo at a velocity about $v\sim 300\,\text{km}/\text{s}=10^{-3}$, and thus we feel the wind of the axion. We see the axion wind as scalar gravitational waves, and the gravitational-wave interferometer detector does have sensitivity to the axion oscillation.

In addition to the dark matter, the dark energy is also a big issue in current physics. In fact, the main energy component of the Universe is the dark energy. Hence, it would be necessary to consider the detectability of the dark matter in the context of dark energy models. Often, the cosmological constant is assumed when we discuss the detectability of dark matter. However, the cosmological constant has several problems, e.g., the fine-tuning problem and coincidence problem. One possibility to resolve these issues is to consider unknown matter such as the quintessence. Unfortunately, there is no natural candidate for the quintessence in particle physics. Therefore, it is worth investigating the possibility that the theory of gravity is different from Einstein gravity on cosmological scales. Thus, we study the detectability of the ultralight axion dark matter in the context of modified gravity.

In our previous paper 16:Aoki , we have discussed the detectability of the ultralight axion dark matter in the framework of $f(R)$ gravity, which is the simplest modified gravity. We derived the gravitational potential sourced by the axion oscillation in the $f(R)\propto R^{2}$ model. Remarkably, we found the resonance between the axion field and the scalaron field, which is the dynamical degree of freedom of $f(R)$ gravity, and the gravitational potential could be amplified dramatically. However, it is not obvious if the resonance behavior can be seen in more general models. In the $R^{2}$ model, the equation determining the gravitational potential is a linear equation under appropriate assumptions, while it should be nonlinear in more realistic models. Hence, it is interesting to investigate whether general models can have resonance behavior. In this paper, we will show that such a model does exist by constructing one specific model.

The paper is organized as follows: In Sec. II, we introduce a notion of axion dark matter and derive its energy-momentum tensor, which is a source term of the field equation for the metric. In Sec. III, we briefly summarize our previous work 16:Aoki . Then, we move on to the scalar-tensor formulation of $f(R)$ gravity and derive the formula (30) for calculating the oscillating part of the gravitational potential. In Sec. IV, we discuss two specific models. We show that the nonlinear resonance significantly enhances the amplitude of oscillation of the gravitational potential. The final section is devoted to our conclusion.

II Axion Oscillation

In this section, we introduce the axion dark matter. In particular, we derive the energy-momentum tensor of the axion field, a source term of the metric field equation.

Let us assume the situation where the dark matter halo is composed of the ultralight axion. Since the occupation number of the axion in the halo is huge, we can treat it as a classical scalar field. The axion field satisfies the Klein-Gordon equation in the flat spacetime at the leading order, and the solution is given by the superposition of plane waves with different wave numbers and frequencies. The wave number has a certain cutoff $k_{\text{max}}\sim mv$ due to the uncertainty principle, roughly set by the inverse of the de Broglie wavelength of an axion particle. Since a typical velocity in the galaxy is $v\sim 10^{-3}$, we can assume that the axion field oscillates monochromatically with the angular frequency corresponding to its mass. Under these assumptions, we can write the axion field as

[TABLE]

where $\phi_{0}(\vec{x})$ is the amplitude and $\alpha(\vec{x})$ is the phase of the oscillation. We can neglect the space dependence of $\phi_{0}(\vec{x})$ and $\alpha(\vec{x})$ at the leading order, and hereafter we will omit the phase $\alpha(\vec{x})$ for simplicity.

The energy density $\rho$ and the pressure $p$ of the axion field are given by

[TABLE]

The pressure oscillates in time with the angular frequency $\omega=2m$. Its amplitude is fixed by the local dark matter density $\rho_{0}$. The energy-momentum tensor of the axion field is then given by

[TABLE]

Finally, we obtain the trace of the energy-momentum tensor of the axion field as

[TABLE]

This is the source term of the field equation. We can use Eq. (5) as long as the axion field minimally couples to gravity.

We use the value $\rho_{0}=0.3\,\text{GeV}/\text{cm}^{3}$ as a typical energy density of the dark matter halo throughout the paper.222 While this value is traditionally used, slightly higher values are reported in some papers, e.g., Refs. 10:Catena ; 10:Salucci ; 10:Pato ; 11:McMillan ; 12:Garbari .

The period of the oscillation corresponds to about 1 year for $m=10^{-22}$ eV, and this time scale is much shorter than the cosmological time scale, i.e., $H_{0}^{-1}\sim 10^{10}$ years. Hence, after averaging the oscillating pressure over the cosmological time scale, the axion behaves as pressureless dust on cosmological scales. Thus, the axion can be a candidate for the dark matter.

III Formula for Gravitational Potential Oscillation in $f(R)$ Gravity

In this section, we will summarize our previous work 16:Aoki and derive the formula (15). Then, we move on to the scalar-tensor formulation of $f(R)$ gravity and derive the main formula (30) for calculating the time-dependent part of the gravitational potential in $f(R)$ gravity.

III.1 Gravitational potential oscillation in $f(R)$ gravity

The action for $f(R)$ gravity is given by

[TABLE]

where $f(R)$ is a function of the Ricci scalar $R$ and $S_{\text{m}}$ is the action for matter fields. We assume $f(R)\ll R$ and $f_{R}\equiv f^{\prime}(R)\ll 1$ so that the deviation from Einstein’s theory is small. Taking the variation of the action with respect to the metric, we obtain the metric field equation:

[TABLE]

where $G_{\mu\nu}\equiv R_{\mu\nu}-(1/2)g_{\mu\nu}R$ is the Einstein tensor, $T_{\mu\nu}$ is the energy-momentum tensor for the matter field, and $\Box\equiv\nabla_{\mu}\nabla^{\mu}$ is the d’Alembert operator. The trace of the field equation reads

[TABLE]

We assume that the spatial derivative of $f_{R}$ is much smaller than the time derivative of it, i.e., $\Box f_{R}\simeq-\ddot{f}_{R}$. This is because the length scale of the dark matter halo, about 10 kpc or larger, is much larger than the time scale of the oscillation, $m^{-1}\sim 0.1$ pc for $m=10^{-22}$ eV. We use this kind of approximation throughout the paper. Under the assumption, the field equation becomes

[TABLE]

where we used $f\ll R$ and $f_{R}\ll 1$. Although we assumed $f_{R}\ll 1$, the contribution of $f_{R}$ to the equation of motion is not so small. This is because the $f_{R}$ term appears as $\ddot{f}_{R}\sim(2m)^{2}f_{R}$ in Eq. (9). Thus, the $f_{R}$ term is enhanced by a factor of $m^{2}/R_{0}=m^{2}/\rho_{0}\sim 10^{17}(m/10^{-22}\,\text{eV})^{2}$ and can be comparable to the Ricci scalar $R$ in Eq. (9). Hereafter, we consider the axion as the matter field. Since the axion field minimally couples to gravity in $f(R)$ gravity, we can use Eq. (5) for the energy-momentum tensor of the axion field.

Since the gravitational potentials are small even in the dark matter halo, they can be treated as perturbations. Let us use the Newtonian gauge for the metric:

[TABLE]

Note that the expansion of the Universe is completely negligible on the scale of the dark matter halo. At the first order of the potentials, the Ricci scalar is calculated as

[TABLE]

Let us write the Ricci scalar as the sum of the time-independent part $R_{0}$ and the time-dependent part $\delta R$:

[TABLE]

where $R_{0}$ is defined as the long-term average of $R$, $R_{0}\equiv\braket{R}$. We also separate the gravitational potential $\Phi$ ($\Psi$) into the time-independent part $\Phi_{0}\equiv\braket{\Phi}$ ($\Psi_{0}\equiv\braket{\Psi}$) and the time-dependent part $\delta\Phi$ ($\delta\Psi$). We have the equation $\Phi_{0}=\Psi_{0}$ from the traceless part of the space-space component of the Einstein equation. Hence, $R_{0}$ can be written as $R_{0}=2\nabla^{2}\Phi_{0}$. The field equation (9) gives $R_{0}\equiv\braket{R}=\rho_{0}$, and this is nothing but the Poisson equation

[TABLE]

Assuming $|\delta\ddot{\Phi}|\gg|\nabla^{2}\delta\Phi|$ and $|\delta\ddot{\Phi}|\gg|\nabla^{2}\delta\Psi|$, $\delta R$ is approximately given by

[TABLE]

Integrating this twice and using the field equation (9), we obtain

[TABLE]

where $\braket{f_{R}}$ is the average value of $f_{R}$. This is the formula for calculating the time-dependent part of the gravitational potential. Therefore, in order to obtain the amplitude of the gravitational potential, we first solve the field equation (9) and then substitute the solution into Eq. (15).

For Einstein theory, we have $f_{R}=0$. Thus, Eq. (15) gives

[TABLE]

where

[TABLE]

The frequency of the gravitational potential is

[TABLE]

This is consistent with the result derived by Khmelnitsky and Rubakov 14:Khmelnitsky .

III.2 Formula for gravitational potential oscillation

While we have discussed the axion oscillation in $f(R)$ gravity, it is convenient to move on to the equivalent scalar-tensor theory for qualitative understanding of the physics 10:Sotiriou ; 10:DeFelice ; 11:Nojiri . In this section, we will reformulate $f(R)$ gravity in terms of the scalar-tensor theory.

The action for $f(R)$ gravity can be rewritten into that of the scalar-tensor theory as

[TABLE]

where $U(\varphi)$ is defined by

[TABLE]

and $A(\varphi)$ is the solution of

[TABLE]

Taking the variation of the action (19) with respect to $\varphi$, we obtain the constraint

[TABLE]

Substituting this into the action (19), we obtain the original $f(R)$ action (6). As long as $f^{\prime\prime}(A)\not=0$, the solution of Eq. (21) is uniquely determined and thus the two theories are completely equivalent. The variation of the action (19) with respect to the metric gives the field equation:

[TABLE]

The trace of this equation gives the equation of motion for the scalar field $\varphi$:

[TABLE]

where we used Eq. (20). Hence, we can interpret $\varphi$, which is often called a scalaron, as an extra dynamical degree of freedom in the theory. Let us define the potential of the scalar field $V(\varphi)$ by

[TABLE]

Using the approximations $f(A)\ll A$ and $\varphi=f^{\prime}(A)\ll 1$, the derivative of the potential $V_{,\varphi}$ is approximated by

[TABLE]

Hereafter we will use this approximate form for $V_{,\varphi}$. Under the assumption $\Box\varphi\simeq-\ddot{\varphi}$, the equation of motion is written as

[TABLE]

where we substituted the energy-momentum tensor (5). Let us introduce the effective potential $V_{\text{eff}}(\varphi)$ by

[TABLE]

The effective potential has a minimum at $R=A(\varphi)=\rho_{0}$ as desired. Using the effective potential, the equation of motion is written in a simple form:

[TABLE]

Since $\varphi=f_{R}$, Eq. (15) can be rewritten in terms of the scalar field $\varphi$ as

[TABLE]

where $\delta\Phi_{\text{E}}\equiv\rho_{0}/8m^{2}$. Therefore, in order to obtain the time-dependent part of the gravitational potential in the scalar-tensor formulation, we first calculate the effective potential from a function $f(R)$ by using Eq. (28), solve the equation of motion (29), and then substitute the solution into the formula (30).

IV Nonlinear Resonant Oscillation

In this section, we will study two specific models. First, we illustrate the resonance of oscillation using the $R^{2}$ model. Then, we show that the exponential model exhibits nonlinear resonance.

IV.1 $R^{2}$ model

We first consider the $f(R)\propto R^{2}$ model, which can be solved analytically. We have already studied this model in the previous paper in the framework of $f(R)$ gravity 16:Aoki . However, since the model is useful for qualitative understanding of general $f(R)$ models, here we again discuss the same model in the scalar-tensor theory.

The model is defined by

[TABLE]

where $M$ is the mass scale of the model, which is independent of the Ricci scalar in this model. The effective potential is calculated as

[TABLE]

where $\varphi_{0}\equiv f^{\prime}(\rho_{0})=\rho_{0}/3M^{2}$. Hence, the mass scale $M$ is nothing but the mass of the scalaron field $\varphi$. We expect this model to capture some features of general models since any analytic function can be approximated by a quadratic function near its minimum. The equation of motion for $\varphi$ now is a linear equation:

[TABLE]

The solution is given by

[TABLE]

Here we focus on the induced solution by the axion oscillation and omitted the homogeneous solutions. This is in part because the homogeneous solutions decay in the expanding Universe by the Hubble friction, and we expect that only the induced solution remains in the present Universe 07:Starobinsky . Since the long-time average of $\varphi$ is $\braket{\varphi}=\varphi_{0}$, Eq. (30) gives

[TABLE]

Let us introduce a parameter $\mu$ as

[TABLE]

When $\mu\ll 1$, the amplitude of $\delta\Phi$ becomes $\delta\Phi_{\text{amp}}\simeq\delta\Phi_{\text{E}}$, and the result (16) is recovered. In the opposite case $\mu\gg 1$, the amplitude is suppressed as $\delta\Phi_{\text{amp}}\simeq(1/\mu^{2})\delta\Phi_{\text{E}}$. Remarkably, when $\mu\simeq 1$, the resonance occurs and the gravitational potential is amplified. The behavior near the resonance point should be strongly dependent on the details of models. To see this, we study another model in the next subsection.

IV.2 Exponential model

In the $R^{2}$ model, the equation of motion for $\varphi$ is linear. In general $f(R)$ models, however, the equation of motion for $\varphi$ should be nonlinear. In the nonlinear cases, it is not obvious whether the resonance occurs. Indeed, in the previous paper 16:Aoki , we saw that there are no stable resonance solutions near $\mu=1$ in the Hu-Sawicki model 07:Hu and the Starobinsky model 07:Starobinsky . In this subsection, we will show that there exist nonlinear models with resonance behavior by constructing one specific model. We will also see that the resonance phenomena acquires new interesting properties when nonlinearity is taken into account.

Let us consider the following model:

[TABLE]

where $M$ is the mass scale of the model at $R=R_{0}$, i.e., $M^{2}=1/3f^{\prime\prime}(R_{0})$, and $\lambda$ is a positive parameter. The exponential-type model was first introduced by Ref. 08:Cognola . The effective potential is calculated as

[TABLE]

where $\varphi_{0}$ is the field value at the minimum of the effective potential:

[TABLE]

The model satisfies the solar system constraint if the following condition is held 10:DeFelice :

[TABLE]

Hence, the model can pass the solar system test even when $M=2\times 10^{-22}\,\text{eV}$ if $\lambda>10^{-3}$. Of course, the larger the mass scale $M$ is, the easier the model passes the solar system test.

We plot the functional form of the effective potential in Fig. 1. The effective potential has a singularity at $\varphi=0$, where $R=A(\varphi)=3V_{,\varphi}$ diverges. This kind of singularity is known as a curvature singularity 07:Briscese ; 08:Nojiri ; 08:Frolov ; 08:Kobayashi , which often appears in $f(R)$ models. We can always remove the singularity without affecting the dynamics on the halo scale by adding a regularization term, e.g., $R^{2}/6\mathcal{M}^{2}$ with a large $\mathcal{M}$ 08:Nojiri ; 08:Dev ; 09:Kobayashi ; 10:Appleby . However, we have no need of such a modification for our purpose to construct a concrete model with the resonance behavior. Therefore, we will not discuss the regularization and focus on what the model (37) predicts. In this model the scalar field can practically move only within the range

[TABLE]

This is because if we start with $V_{\text{eff}}(\varphi)>0$, i.e., $\varphi/\varphi_{0}>e$, the scalar field easily hits the curvature singularity. Hence possible amplitudes of the scalar field are roughly limited to $(e/2)|\varphi_{0}|\sim|\varphi_{0}|$, and the maximum amplitude of the gravitational potential is approximated as

[TABLE]

where we assumed the second term in Eq. (30) is dominant. The amplitude could become large for sufficiently small $\lambda$ even when we fix $\mu\sim 1$, i.e., $M\sim 2m$.

Now, we will show that the solutions with the maximum amplitude (42) do exist. All we have to do is to solve the nonlinear equation of motion

[TABLE]

where we used the dimensionless quantities $\tau=Mt$, $\mu=2m/M$, $\bar{\varphi}=\varphi/\varphi_{0}$, and the prime denotes a derivative with respect to $\tau$. In the case of a linear forced oscillator, the general solution is given by the superposition of an induced solution and homogeneous solutions. However, in nonlinear cases, the superposition of two solutions no longer gives a solution. Therefore, it is hard to solve nonlinear systems analytically. Indeed, we have to rely on a perturbative method or numerical calculations.

First, in order to get an intuition, we use the following perturbative method. The deviation from the minimum $\chi$ is defined by $\bar{\varphi}=1+\chi$. When $\chi\ll 1$, we can approximate the equation by

[TABLE]

This can be analyzed perturbatively. First, we seek a resonant solution around $\mu\sim 1$. To this aim, we rewrite Eq. (44) as

[TABLE]

We regard the terms in the right-hand side as small perturbations. Then, we can solve Eq. (45) with the expansion

[TABLE]

Substituting the series into Eq. (45), we obtain the lowest order solution

[TABLE]

At the next order, we have secular sources

[TABLE]

If this secular term remains, we will have a secular solution with the growing amplitude. In other words, we have to renormalize this secular evolution into the frequency. This can be achieved by imposing the condition

[TABLE]

Apparently, this allows the order one solution $A\simeq 2\sqrt{\mu^{2}-1}$ for $\mu\geq 1$. Thus, we have shown that there exists the resonant oscillation with the amplitude close to the maximum one. In principle, we can repeat the same analysis to find other resonant solutions. We will soon show explicit examples with the numerical method.

Next, we investigate the same system numerically. We seek periodic solutions satisfying $(\bar{\varphi}(T),\bar{\varphi}^{\prime}(T))=(\bar{\varphi}(0),\bar{\varphi}^{\prime}(0))$, where $T\equiv 2\pi/\mu$ is the period of the oscillation. In other words, we look for the solutions with closed orbit in the phase space $(\bar{\varphi},\bar{\varphi}^{\prime})$. Note that in the $R^{2}$ model, this condition removes homogeneous solutions. In order to find such solutions, in general, we should study the map $(\bar{\varphi}(0),\bar{\varphi}^{\prime}(0))\mapsto(\bar{\varphi}(T),\bar{\varphi}^{\prime}(T))$ and find its fixed points. However, since the equation of motion (27) now has the time reflection symmetry $t\to-t$ thanks to the absence of a friction term in addition to the time translation symmetry $t\to t+T$, the orbit is closed if $\bar{\varphi}^{\prime}(T/2)=0$ when starting with the initial condition $(\bar{\varphi}(0),\bar{\varphi}^{\prime}(0))=(\bar{\varphi}_{i},0)$. Note that the reverse statement is not always true; there exist periodic solutions that do not satisfy $\bar{\varphi}^{\prime}(T/2)=0$ as will be seen later.

In Fig. 2, we plot the value $\bar{\varphi}^{\prime}(T/2)$ as a function of $(\mu,\bar{\varphi}_{i})$ with different values of $\lambda$. The red (blue) regions correspond to positive (negative) values of $\bar{\varphi}^{\prime}(T/2)$. Since $\bar{\varphi}^{\prime}(T/2)$ is a smooth function of $(\mu,\bar{\varphi}_{i})$, there are boundaries with $\bar{\varphi}^{\prime}(T/2)=0$ between two regions, which correspond to closed-orbit solutions. Hence, Fig. 2 shows that solutions exist with amplitude $|\bar{\varphi}|=|\varphi/\varphi_{0}|\sim\mathcal{O}(1)$. The figure also shows that there are three solutions around $\mu\gtrsim 1$. This is a general feature of nonlinear forced oscillators 60:Landau . Which solution is selected should be determined by initial conditions. We also show the corresponding resonance curves in Fig. 3. The resonance curves are bent by nonlinearity compared to the $R^{2}$ case. This plot shows that the resonance curves in the exponential model cannot be distinguished from those in the $R^{2}$ model except for $\mu\sim 1$. Namely, while the behavior near the resonance point is strongly dependent on models, the $R^{2}$ model extracts some features of general models as desired.

Remarkably, in the nonlinear case, there appear new resonances at $\mu=q/p$ for positive integers $p$ and $q$ 60:Landau . As an example, we show the same plot as Fig. 2 for the series with $p=5$ in Fig. 4, where we chose $\lambda=10^{-3}$. Note that the plotted function $\bar{\varphi}^{\prime}(T/2)$ is replaced by $\bar{\varphi}^{\prime}(5T/2)$ in Fig. 4. In these cases, the resonance of scalaron field oscillation occurs when the relation $M\sim 2m(p/q)$ holds. Hence, the observed frequency becomes higher or lower than $2m$. In the former cases, the frequencies could be accessible by ground-based gravitational-wave detectors such as LIGO.

V Conclusion

In this paper, we have studied the axion oscillation in $f(R)$ gravity. The point was that the pressure of the axion field oscillating in time produces the oscillation of the gravitational potential. We have derived the formula (30) for calculating the time-dependent part of the gravitational potential. It turns out that the amplitude of the oscillation of the gravitational potential is determined by the amplitude of the scalar field (scalaron) in the scalar-tensor formulation. Remarkably, the amplitude could be amplified dramatically by the resonance when the mass of the scalaron is sufficiently close to twice the axion mass. We have also shown that there appear subharmonics resonances at $\mu=q/p$ in the case of nonlinear models by studying one specific model (37) as an example.

In this paper, we have studied two $f(R)$ models on halo scales. We should mention that we did not attempt to construct complete models which also explain the accelerated expansion of the present Universe. In fact, in order to explain the dark energy, we need to modify the models on cosmological scales. Since the cosmological critical density is much lower than the dark matter density in the halo, $\rho_{\text{cr}}/\rho_{0}\sim 10^{-5}$, we expect that the modification of the functional form of $f(R)$ on cosmological scales does not affect the dynamics on halo scales.

The oscillating gravitational potential can, in principle, be observed by using gravitational-wave detectors. So far, two ideas for direct detection of the ultralight axion dark matter through the gravitational interaction have been proposed: pulsar timing array experiments 14:Khmelnitsky and laser interferometers 17:Aoki .333 Recently, a new detection method with binary pulsars was also proposed 16:Blas . The nonlinear resonance that we found in this paper further enhanced the detectability of the ultralight axion.

Although we have focussed on the present Universe, there may be other phenomena caused by the oscillating axion field. In particular, we should revisit the structure formation process in the presence of the ultralight axion field. Structure formation in modified gravity is often discussed assuming the validity of static or quasistatic approximation on galactic scales. However, if the ultralight axion is the dark matter, the situation is quite different. We should take into account the oscillating pressure of the axion field when we discuss structure formation in modified gravity.

Acknowledgements.

This work was in part supported by JSPS Grant-in-Aid for JSPS Research Fellow Grant No. 17J00568 (A.A.), JSPS KAKENHI Grant No. 17H02894, and MEXT KAKENHI Grant No. 15H05895 (J.S.).