Time evolution of the local gravitational parameters and gravitational wave polarizations in a relativistic MOND theory

TL;DR

This paper revisits a relativistic MOND theory, demonstrating time-varying gravitational parameters and confirming the existence of only two tensor polarizations in gravitational waves, aligning with general relativity predictions.

Contribution

It shows how to obtain time-varying Newtonian and MOND parameters by relaxing static assumptions and reproduces observed redshift dependence of MOND acceleration scale.

Findings

Time-varying effective Newtonian constant $G_N$ and MOND acceleration scale $a_{MOND}$.

Reproduction of redshift dependence of $a_{MOND}$ from simulations.

Confirmation of two tensor polarizations and gravitational wave speed equal to light speed.

Abstract

The recently proposed Skordis-Z\l{}o\'{s}nik theory is the first relativistic MOND theory that can recover the success of the standard $\Lambda$CDM model at matching observations of the cosmic microwave background. This paper aims to revisit the Newtonian and MOND approximations and the gravitational wave analysis of the theory. For the local gravitational parameters, we show that one could obtain both time-varying effective Newtonian gravitational \textit{constant} $G_\textrm{N}$ and time-varying characteristic MOND acceleration scale $a_\textrm{MOND}$, by relaxing the static assumption extensively adopted in the literature. Specially, we successfully demonstrate how to reproduce the redshift dependence of $a_\textrm{MOND}$ observed in the \textit{Magneticum} cold dark matter simulations. For the gravitational waves, we show that there are only two tensor polarizations, and reconfirm that its speed is equal to the speed of light.

Tables

Table 1. Table 1: Models discussing the possible cosmological evolution of a mond subscript 𝑎 mond a_{\textsc{mond}} . Note that Eqs. ( 12 ) and ( 13 ) are extensions of the original Skordis-Złośnik theory.

00footnotetext: Most of the existing relativistic MOND theories give constant $a_{\text{mond}}$. Here we only consider the models that the constant $a_{\text{mond}}$ can still be obtained after analyzing the relevant cosmic background evolutions. 00footnotetext: Here $l_{\text{de}}$ is the characteristic length scale of dark energy, which is of the order of ${\Lambda }^{-1/2}$, i.e., $c/H_0$, at today and could be time-varying in the dynamical models.
$a_{\text{mond}}\left(z\right)$	Skordis-Złośnik theory	Other theories & Phenomenological motivations
$\sim \mathrm{const}.$11footnotemark: 1	$𝒥=c_1{𝒴}^{3/2}/𝒬$	TeVeS theory [58, 33]                                               a subclass of nonlocal MOND models [59]
$\propto {\left(1+z\right)}^3$	$𝒥=c_1{𝒴}^{3/2}/\left(𝒬-{𝒬}_c\right)$	—
$\propto cH\left(z\right)$	$𝒥=c_1{𝒴}^{3/2}/{\nabla }_{\mu }{A}^{\mu }$	the numerical coincidence between $a_{\text{mond},0}$ and $c{H}_0$ [1, 60] a subclass of nonlocal MOND models [59]
$\propto c^2/l_{\text{de}}$22footnotemark: 2	$𝒥=c_1{𝒴}^{3/2}/\sqrt{V_{\text{de}}}$	the numerical coincidence between $a_{\text{mond},0}$ and $c^2/l_{\text{de},0}$ [60] relativistic theories linking MOND to dark energy [61, 62]
Magneticum	Eq. (11)	—

Equations

$${\mathcal{L}}_{\text{textsc}mond}=-\frac{K_{\text{textsc}b}}{2}{F}^{\mu \nu }{F}_{\mu \nu }+2(2-K_{\text{textsc}b})J^{\mu }{\nabla }_{\mu }\varphi$$

$$-(2-K_{\text{textsc}b})\mathcal{Y}-\mathcal{F}(\mathcal{Y},\mathcal{Q})-\lambda (A^{\mu }{A}_{\mu }+1),$$

$$G_{\mu \nu }-K_{\text{textsc}b}{F}_{\mu }^{\alpha }{F}_{\nu \alpha }+(2-K_{\text{textsc}b})\{2J_{(\mu }{\nabla }_{\nu )}\varphi -A_{\mu }{A}_{\nu }\square \varphi$$

$$+2[A_{(\mu }{\nabla }_{\nu )}{A}_{\alpha }-A_{(\mu }{\nabla }_{|\alpha |}{A}_{\nu )}]{\nabla }^{\alpha }\varphi \}-{\mathcal{F}}_{\mathcal{Q}}{A}_{(\mu }{\nabla }_{\nu )}\varphi$$

$$-(2-K_{\text{textsc}b}+{\mathcal{F}}_{\mathcal{Y}})[{\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi +2\mathcal{Q}{A}_{(\mu }{\nabla }_{\nu )}\varphi ]$$

$$-\lambda A_{\mu }{A}_{\nu }-g_{\mu \nu }{\mathcal{L}}_{\text{textsc}mond}/2=\kappa T_{\mu \nu },$$

$${\nabla }_{\mu }{\mathcal{I}}^{\mu }=0,$$

$$K_{\text{textsc}b}{\nabla }_{\nu }{F}^{\nu \mu }+(2-K_{\text{textsc}b})[({\nabla }^{\mu }{A}_{\nu }){\nabla }^{\nu }\varphi -{\nabla }_{\nu }(A^{\nu }{\nabla }^{\mu }\varphi )]$$

$$-\left[(2-K_{\text{textsc}b}+{\mathcal{F}}_{\mathcal{Y}})\mathcal{Q}+{\mathcal{F}}_{\mathcal{Q}}/2\right]{\nabla }^{\mu }\varphi -\lambda A^{\mu }=0.$$

$$A^{\mu }{A}_{\mu }+1=0.$$

$$H^2=\frac{8\pi \tilde{G}}{3}{\rho }_{\mathrm{m}}+\frac{c^2}{6}(\mathcal{F}-\mathcal{Q}{\mathcal{F}}_{\mathcal{Q}}),$$

$$H^2+2\frac{\ddot{a}}{a}=-\frac{8\pi \tilde{G}}{c^2}{p}_{\mathrm{m}}+\frac{1}{2}\mathcal{F}{c}^2,$$

$$\frac{\mathrm{d}{\mathcal{F}}_{\mathcal{Q}}}{\mathrm{d}t}+3H{\mathcal{F}}_{\mathcal{Q}}=0.$$

$${\rho }_{\text{textsc}mond}=\frac{c^2}{16\pi \tilde{G}}(\mathcal{F}-\mathcal{Q}{\mathcal{F}}_{\mathcal{Q}}),p_{\text{textsc}mond}=-\frac{\mathcal{F}{c}^4}{16\pi \tilde{G}}.$$

$$\mathcal{F}(\mathcal{Y},\mathcal{Q})=(2-K_{\text{textsc}b})\mathcal{J}(\mathcal{Y},\mathcal{Q})-2\mathcal{K}(\mathcal{Q}).$$

$$\frac{\mathcal{Q}-{\mathcal{Q}}_c}{{\mathcal{Q}}_c}=\frac{1}{2}\beta -\frac{3}{8}{\beta }^2+\mathcal{O}({\beta }^3),$$

$${\rho }_{\text{textsc}mond}=\frac{{\mathcal{K}}_2{c}^2{\mathcal{Q}}_c^2}{8\pi \tilde{G}}\left[\beta +\frac{1}{4}{\beta }^2+\mathcal{O}({\beta }^3)\right],$$

$$p_{\text{textsc}mond}=\frac{{\mathcal{K}}_2{c}^4{\mathcal{Q}}_c^2}{32\pi \tilde{G}}{\beta }^2+\mathcal{O}({\beta }^3),$$

$$w_{\text{textsc}mond}=\frac{\beta }{4}+\mathcal{O}({\beta }^2),$$

$$\beta =\frac{3{\Omega }_{\text{textsc}mond,0}{H}_0^2}{{\mathcal{K}}_2{c}^2{\mathcal{Q}}_c^2}{\left(\frac{a_0}{a}\right)}^3\ll 1.$$

$$S^{(2)}$$

$$+c^4\mathcal{J}(\mathcal{Y},\bar{\mathcal{Q}})+\mathrm{mass}\mathrm{terms}]+\rho \Phi \},$$

$$\mathcal{Y}=|\nabla \phi {|}^2+\mathrm{mass}\mathrm{terms},$$

$${\nabla }^2\hat{\Phi }=\frac{8\pi \tilde{G}}{2-K_{\text{textsc}b}}\rho ,$$

$$\nabla [{\mathcal{J}}_{\mathcal{Y}}\nabla \phi ]=\frac{8\pi \tilde{G}}{(2-K_{\text{textsc}b})c^2}\rho ,$$

$$G_{\text{textsc}n}=\frac{2\tilde{G}}{2-K_{\text{textsc}b}}(1+\frac{1}{{\lambda }_s}).$$

$$\frac{{\dot{G}}_{\text{textsc}n}}{G_{\text{textsc}n}}\approx -\frac{{\dot{\lambda }}_s}{{\lambda }_s^2}\approx \frac{3p{H}_0}{{\lambda }_{s,0}}.$$

$$\mathcal{J}(\mathcal{Y},\mathcal{Q})=\frac{{\mathcal{Y}}^{3/2}}{c_1{\mathcal{Q}}_c+c_2{\mathcal{Q}}_c^{1/3}(\mathcal{Q}-{\mathcal{Q}}_c{)}^{2/3}}.$$

$$\mathcal{J}(\mathcal{Y},{\nabla }_{\mu }{A}^{\mu })=c_1{\mathcal{Y}}^{3/2}/{\nabla }_{\mu }{A}^{\mu }$$

$$\mathcal{J}=c_1{\mathcal{Y}}^{3/2}/\sqrt{V_{\text{textsc}de}},$$

$$h_{tt}=2\varphi ,$$

$$h_{tj}={\beta }_j+{\partial }_j\gamma ,$$

$$h_{jk}=h_{jk}^{\mathrm{TT}}+\frac{1}{3}H{\delta }_{jk}+{\partial }_{(j}{ϵ}_{k)}+\left({\partial }_j{\partial }_k-\frac{1}{3}{\delta }_{jk}{\nabla }^2\right)\rho ,$$

$$a^t=\frac{1}{2}{h}_{tt}=\varphi ,$$

$$a^j={\mu }^j+{\partial }^j\omega .$$

$$h_{\mu \nu }\to h_{\mu \nu }-{\partial }_{\mu }{\xi }_{\nu }-{\partial }_{\nu }{\xi }_{\mu },$$

$$a^{\mu }\to a^{\mu }+{\bar{A}}^{\nu }{\partial }_{\nu }{\xi }^{\mu },$$

$$\phi \to \phi .$$

$$\phi ,h_{jk}^{\mathrm{TT}},$$

$$\Phi =-\varphi +\dot{\gamma }-\frac{1}{2}\ddot{\rho },$$

$$\Theta =\frac{1}{3}(H-{\nabla }^2\rho ),$$

$${\Xi }_j={\beta }_j-\frac{1}{2}{\dot{ϵ}}_j,$$

$${\Sigma }_j={\beta }_j+{\mu }_j,$$

$$\Omega =\omega +\frac{1}{2}\dot{\rho }.$$

Full text

I Introduction

Modified Newtonian dynamics (MOND) is an alternative to the dark matter paradigm, through the modification of Newton’s law of universal gravitation or Newton’s second law of motion [1, 2, 3]. The former belongs to the traditional modified gravity, and construction of its relativistic counterpart has been extensively discussed. Bekenstein and Milgrom [4] proposed the first one. However, there are two major problems: the acausal problem [4] and the gravitational lensing problem [5]. Further modifications of the theory have been proposed to address these issues, such as the phase coupling [6] and disformal transformations [7, 8, 9, 10]. These attempts made great progress in shaping the relativistic MOND theory and explaining the local gravitational phenomena [11]. However, for the cosmological linear perturbations, no such theory has been shown to successfully fit all the current data about the cosmic microwave background anisotropies and matter power spectra [12, 13, 14, 15, 16]. Recently, Skordis and Złośnik [17] proposed a new MOND theory to address this observational fitting problem. Analysis of this theory is the topic of this paper. In addition, we note that modification of Newton’s second law still requires further development to arrive at a complete and observationally accepted theory [18, 19, 20, 21].

MOND theories generally predict an universal radial acceleration relation (RAR) — correlation between the observed radial acceleration and that predicted by baryons with Newtonian gravity. McGaugh et al. [22] first observed the RAR in the SPARC database, and further data confirmed the conclusion [23, 24]. This may be regarded as an observational evidence supporting MOND. However, after McGaugh et al. [22], the same relation was also observed in the $N$-body simulations of cold dark matter (CDM) [25, 26, 27, 28]. The mass discrepancy acceleration relation, which is similar to RAR, was also predicted by MOND, and observed in both observations [29] and CDM simulations [30, 31]. In particular, Keller and Wadsley [26] found that the CDM simulated RAR depends on the cosmological redshift. This result implies that, in the framework of CDM, rotating galaxies still satisfy the universal RAR at high redshifts. However, the parameter characterizing the acceleration scale in RAR is redshift-dependent. Recently, Mayer et al. [32] presented an explicit redshift evolution of this characteristic MOND acceleration scale $a_{\textsc{mond}}$ in the Magneticum CDM simulations.

In the relativistic MOND theories, $a_{\textsc{mond}}$ is a parameter and could be time-varying. Milgrom [1] conjectured that $a_{\textsc{mond}}\propto cH$ based on the numerical coincidence of their values at today. In the framework of TeVeS theory (a relativistic MOND theory) [10], Bekenstein and Sagi [33] analyzed this issue after considering the cosmological background evolution of the relevant fields. They found that $a_{\textsc{mond}}$ changes on time scales much longer than the Hubble timescale. In this paper, we present the first analysis of the possible time evolution of the local Newtonian and MOND parameters in the Skordis-Złośnik theory [17]. The method is principally the same as that in [33]. Our result demonstrates how to reproduce the Magneticum redshift dependence [32] in such relativistic MOND theory.

The first direct detection of the gravitational wave signal GW150914 has marked the new era of gravitational wave astronomy [34]. In general relativity, there exist two well-known gravitational wave polarizations (plus and cross), traveling at the speed of light. GW170814 and GW170817 observations confirmed these predictions [35, 36, 37]. In this paper, we present a gauge-invariant gravitational wave analysis for the Skordis-Złośnik theory, in which the polarization content and the propagation speed are determined.

This paper is organized as follows. Section II introduces the Skordis-Złośnik MOND theory and summarizes the cosmic background evolutions. Note that, in principle, most of the results given in this section were obtained by [17]. This section is retained to provide a clear basis for our subsequent discussions. Section III analyzes the Newtonian and MOND approximations. Section IV discusses gravitational waves in the theory. Conclusions are presented in Sec. V.

Throughout this paper, we adopt the Hubble constant $H_{0}=67.4\,{\rm km}/{\rm s}/{\rm Mpc}$ and denote $h$ as its reduced value [38]. The subscript [math] indicates the cosmological redshift $z=0$. In order to compare with observations, we adopt the SI Units and retain all physical constants in Sec. II and Sec. III. We set the speed of light $c=1$ in Sec. IV for simplicity.

II The theory and cosmic evolutions

The Skordis-Złośnik MOND theory is constructed based on a scalar field $\phi$ and a vector field $A_{\mu}$ [17]. Its action is of the form $S=\int{\rm d}^{4}x\sqrt{-g}\left[R+\mathcal{L}_{\textsc{mond}}\right]/2\kappa+S_{\rm m}$, where $\kappa=8\pi\tilde{G}/c^{4}$ and $\tilde{G}$ is a constant with the same dimension of the Newtonian gravitational constant $G_{\textsc{n}}$. The MOND Lagrangian reads

[TABLE]

where $F_{\mu\nu}=2\nabla_{[\mu}A_{\nu]}$, $J_{\mu}=A^{\alpha}\nabla_{\alpha}A_{\mu}$, $\mathcal{Y}=q^{\mu\nu}\nabla_{\mu}\phi\nabla_{\nu}\phi$, $q^{\mu\nu}=g^{\mu\nu}+A^{\mu}A^{\nu}$, $\mathcal{Q}=A^{\mu}\nabla_{\mu}\phi$, $\mathcal{F}(\mathcal{Y},\mathcal{Q})$ is an arbitrary function, $\lambda$ is the Lagrange multiplier (a scalar), $K_{\textsc{b}}$ is a dimensionless constant. In our conventions, the dimension of $A_{\mu}$ relates to the metric ($[A^{\mu}A_{\mu}]=1$), $\phi$ is dimensionless, and $[\mathcal{Y}]=[\mathcal{Q}^{2}]=[\mathcal{F}]=[\lambda]={\rm length}^{-2}$.

The field equations can be derived from the variational principle. Variation of the action with respect to the metric gives the gravitational field equations

[TABLE]

Energy and momentum conservation $\nabla_{\nu}T^{\mu\nu}=0$ can be directly derived from Eq. (2).

As we discussed in Sec. I, one goal of this paper is to study the time evolution of the local gravitational parameters in the Skordis-Złośnik MOND theory. In a relativistic theory, parameters describing the local gravitational system can be time-varying due to the cosmic evolution of the relevant fields. For example, the effective Newtonian gravitational constant is time-varying in scalar-tensor gravity [39, 40, 41, 42, 43] and nonlocal gravity [44, 45, 46]. Here we summarize the cosmic background evolution for the Skordis-Złośnik MOND theory. We assume the Universe is described by the flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric ${\rm d}s^{2}=-c^{2}{\rm d}t^{2}+a^{2}{\rm d}\mathbf{x}^{2}$, where $a=a(t)$. To be consistent with Eq. (2d), we assume $A_{\mu}=[-c,0,0,0]$ for the vector field. The scalar field is assumed to be $\phi=\phi(t)$. For the normal matters, we adopt $T^{\mu}_{\ \nu}={\rm diag}\{-\rho_{\rm m}c^{2},p_{\rm m},p_{\rm m},p_{\rm m}\}$ [47]. Substituting Eq. (2) into Eq. (2) eliminates $\lambda$. Then, substituting the above assumptions into the result, we obtain

[TABLE]

To test self-consistency, we confirm that Eq. (2) gives only trivial results except one constraint equation on $\lambda$. Energy conservation $\dot{\rho}_{\rm m}+3H(\rho_{\rm m}+p_{\rm m}/c^{2})=0$ can be derived from Eq. (3) for arbitrary $\mathcal{F}$ function. In other words, Eqs. (3a), (3c) and the matter energy conservation equation form a complete and self-consistent set. Based on Eq. (3), we can define the effective MOND (dark matter) mass density and pressure as

[TABLE]

Then Eq. (3) can be rewritten as the two conventional Friedmann equations and one effective MOND energy conservation equation. The MOND relative mass density is defined as $\Omega_{\textsc{mond}}\equiv 8\pi\tilde{G}\rho_{\textsc{mond}}/(3H^{2})$.

In order to reveal the key properties of the cosmic background evolution, and to be consistent with the conventions adopted in [17], we rewrite the function

[TABLE]

Here the first term satisfies $\mathcal{J}(0,\mathcal{Q})=0$, and is used to produce the MOND behavior (see Sec. III). For the second term, we adopt the Higgs-like function $\mathcal{K}(\mathcal{Q})=(\mathcal{K}_{2}/4\mathcal{Q}_{c}^{2})(\mathcal{Q}^{2}-\mathcal{Q}_{c}^{2})^{2}$, where $\mathcal{K}_{2}$ and $\mathcal{Q}_{c}$ are constant parameters [17, 48]. The $\mathcal{Q}$ approaches to $\mathcal{Q}_{c}$ in the infinite future. In the late-time Universe, we can adopt $\rho_{\rm m}=\rho_{\rm baryon}\propto a^{-3}$. Meanwhile we add the cosmological constant $\Lambda$ to Eq. (3a). Then Eq. (3c) completely determines the cosmological evolution of $\mathcal{Q}$. The solid lines in Fig. 1 show the numerical results for this ordinary differential system. The initial condition of $\mathcal{Q}$ and parameters of baryon and $\Lambda$ are set so that $\Omega_{\rm Baryon}h^{2}=0.0224$ and $\Omega_{\textsc{mond}}h^{2}=0.120$ at today [38]. The MOND parameters are $\mathcal{K}_{2}=8.5\times 10^{8}$ and $\mathcal{Q}_{c}=1\,{\rm Mpc}^{-1}$, which guarantee good fits to the Planck CMB measurements and SDSS matter power spectra results [17]. We set $\tilde{G}=G_{\textsc{n}}$, which corresponds to a special kind of $\mathcal{J}$ function (see Sec. III). Note that the cosmic background evolutions are independent of $\mathcal{J}$. We emphasize that high precision calculations are required to suppress numerical errors. The bottom part plots the relative energy density $\Omega_{i}$ of each component together with the result of dark matter in the standard $\Lambda$CDM model [38]. The coincidence between MOND and CDM indicates that such MOND can behave as cold as the CDM in the expanding Universe. The top part plots the evolution of $\mathcal{Q}$. Considering the extremely small value of the dimensionless $y$-axis, we conclude that, during the late-time era, $\mathcal{Q}$ keeps almost constant while the variation of $\mathcal{Q}-\mathcal{Q}_{c}$ is considerable. Generalizing the Taylor expansion discussed in [17], we obtain

[TABLE]

Note that $1+z=a_{0}/a$, where $z$ is the cosmological redshift. The $\Omega_{\textsc{mond},0}$ appearing in Eq. (6e) can be regarded as a boundary condition of the differential equation (3). The above result confirms that the MOND is cold for the previous parameter settings. In the top part of Fig. 1, we also plot the leading term of Eq. (6a), and the result shows it is a good approximation. Especially, the leading terms of Eq. (6) are valid for a general $\mathcal{K}(\mathcal{Q})$ once it satisfies $\mathcal{K}\approx\mathcal{K}_{2}(\mathcal{Q}-\mathcal{Q}_{c})^{2}$ when $\mathcal{Q}\rightarrow\mathcal{Q}_{c}$.

III Newtonian and MOND analysis

In the Skordis-Złośnik theory, the form of $\mathcal{J}(\mathcal{Y},\mathcal{Q})$ determines the local gravitational behaviors. Skordis and Złośnik [17] pointed out two key properties. For physically acceptable scenarios, in the strong field region [49], the scalar field is described by the tracking or screening solution, which corresponds to the strong asymptotic expression $\mathcal{J}\rightarrow\mathcal{Y}$ or $\mathcal{J}\rightarrow\mathcal{Y}^{p}$ with $p\geq 3/2$, respectively. The strong field solution determines the relation between $\tilde{G}$ and $G_{\textsc{n}}$. In the weak field region [49], MOND appears if $\mathcal{J}\rightarrow\mathcal{Y}^{3/2}$. These analyzes assumed that $\mathcal{Q}$ appearing in $\mathcal{J}$ reaches its cosmological minimum $\mathcal{Q}_{c}$. However, this may be invalid if variables such as $\mathcal{Q}-\mathcal{Q}_{c}$ appear in $\mathcal{J}$ (see Fig. 1 for the evolutions). Similar to the time-varying $G_{\textsc{n}}$ in the scalar-tensor theory [39, 40, 41, 42, 43], relaxing this static assumption may make the local MOND parameters time-varying.

Considering the great success of the $\Lambda$CDM model in both theories and observations [50, 51, 52], we wish to answer what kind of $\mathcal{J}(\mathcal{Y},\mathcal{Q})$ can reproduce the redshift dependence found in the Magneticum $\Lambda$CDM simulations [32]. However, we emphasize that the Magneticum trend has not been confirmed by observations. The difference between the MOND predictions and the Magneticum trend does not mean the failure of either theory. Instead, this possible difference provides an indicator to distinguish between MOND and $\Lambda$CDM observationally in the future. Besides MOND, dark sector models beyond $\Lambda$CDM, such as dynamical dark energy and ultralight dark matter, might also predict a different $a_{\textsc{mond}}$–$z$ relation. This is due to the fact that these models could affect galaxy formation [53, 54]. A complete model dictionary of $a_{\textsc{mond}}(z)$ is useful for future observational tests. The present paper only focuses on the part about the Skordis-Złośnik MOND theory.

Following [17], we adopt the perturbed metric ${\rm d}s^{2}=-c^{2}(1+2\Phi/c^{2}){\rm d}t^{2}+(1-2\Phi/c^{2}){\rm d}\mathbf{x}^{2}$, where the first-order infinitesimal $\Phi=\Phi(\mathbf{x})$. The vector field is assumed to be $A_{\mu}=[-c(1+\Phi/c^{2}),0,0,0]$, which is consistent with Eq. (2d). The scalar field is assumed to be $\phi=\bar{\phi}(t)+\varphi$, where the bar means cosmic background value and the first-order infinitesimal $\varphi=\varphi(\mathbf{x})$. The time derivative of the first-order infinitesimal is ignored because it is much smaller than the corresponding space derivative [17, 46]. The possible time dependence of the local MOND parameters is encoded in $\bar{\phi}(t)$, or strictly $\bar{\mathcal{Q}}(t)$. Note that we no longer assume $\bar{\mathcal{Q}}=\mathcal{Q}_{c}$. Calculating the quadratic terms in the action with the above perturbations, we obtain

[TABLE]

and $\hat{\Phi}=\Phi-\varphi c^{2}$, $\rho$ is the local baryon mass density, and the mass terms indicate terms like ${\rm const.}\times\Phi^{n}$. The $\mathcal{K}(\mathcal{Q})$ only contributes to the mass terms because $\mathcal{Q}=\bar{\mathcal{Q}}\cdot(1-\Phi/c^{2}+2\Phi^{2}/c^{4})$. For the same reason, we can rewrite $\bar{\mathcal{Q}}$ as $\mathcal{Q}$ in Eq. (7), and only consider the perturbation of $\mathcal{Y}$ in the following discussions. Hereafter we ignore the mass terms. This is reasonable because suitable parameters can indeed suppress the corresponding influences on the Newtonian and MOND dynamics [17]. Integration by parts is used to eliminate the second derivative terms (e.g., $\Phi\nabla^{2}\Phi$) and obtain the above results. Equation (7) recovers Eq. (6) in [17] when $\bar{\mathcal{Q}}=\mathcal{Q}_{c}$. Hereafter we omit the bar in $\bar{\mathcal{Q}}$ and adopt $\mathcal{J}_{\mathcal{Y}}=\partial\mathcal{J}(\mathcal{Y},\mathcal{Q})/\partial\mathcal{Y}$. Variation of $S^{(2)}$ with respect to $\hat{\Phi}$ and $\varphi$, we obtain

[TABLE]

respectively. The Skordis-Złośnik theory is written in the Einstein frame with minimally coupling between matter and other fields. Therefore, $\Phi=\hat{\Phi}+\varphi c^{2}$ is the physical gravitational potential. In the weak field region, if $\mathcal{J}\propto\mathcal{Y}^{3/2}$, i.e., $\mathcal{J}_{\mathcal{Y}}\propto|\nabla\varphi|$, then $\varphi$ dominates $\Phi$ and produces the MOND behavior [4, 17].

Here we discuss the possible time evolution of $G_{\textsc{n}}$ and $a_{\textsc{mond}}$. Comparison of Eq. (8) and Poisson equation in strong field region determines $G_{\textsc{n}}$. In the scaling case [17], we assume $\mathcal{J}\rightarrow\lambda_{s}\mathcal{Y}$, where the dimensionless variable $\lambda_{s}=\lambda_{s}(\mathcal{Q})$. Then Eq. (8) gives $\varphi c^{2}\rightarrow\hat{\Phi}/\lambda_{s}$ and

[TABLE]

Note that $\tilde{G}$ is a constant introduced in the action, and $G_{\textsc{n}}$ could be time-varying because of its dependence on $\lambda_{s}$. Considering $(\mathcal{Q}-\mathcal{Q}_{c})/\mathcal{Q}_{c}\ll 1$ (see Fig. 1), if $\lambda_{s}\propto\mathcal{Q}^{p}$, where $p$ is a constant, then the time evolution of $G_{\textsc{n}}$ is unobservable. However, if $\lambda_{s}\propto(\mathcal{Q}-\mathcal{Q}_{c})^{p}\propto a^{-3p}$, then Eq. (9) gives

[TABLE]

in which we assumed $\lambda_{s}\gg 1$ and the last equality is valid at the low redshift Universe. Current observations give $|\dot{G}_{\textsc{n}}/G_{\textsc{n}}|\lesssim 10^{-12}\,{\rm yr}^{-1}\approx 0.01H_{0}$ [55, 56, 57]. Therefore we require $\lambda_{s,0}\gtrsim 100$ for this case. The screening case [17] corresponds to $\lambda_{s}=\infty$, and results in an exactly constant $G_{\textsc{n}}$. This requires $\mathcal{J}\propto\mathcal{Y}^{p}$, where $p\geq 3/2$ [17].

The parameter $a_{\textsc{mond}}$ is determined in the weak field region, in which $\varphi$ dominates $\Phi$. Hereafter, for simplicity, we adopt $\mathcal{J}\propto\mathcal{Y}^{3/2}$ throughout the MOND region to the Newtonian region. Considering Eqs. (7b), (8b), (9) with $\lambda_{s}=\infty$, and Eq. (3) in [4], we see that the coefficient of $\mathcal{J}\propto\mathcal{Y}^{3/2}$ equals to $2c^{2}/(3a_{\textsc{mond}})$. On the other hand, the coefficient can be written as a function of $\mathcal{Q}$. Theoretically, $a_{\textsc{mond}}$ could be redshift-dependent. Here we discuss several explicit cases. Considering the dimensions of the variables, one of the simplest cases is $\mathcal{J}=c_{1}\mathcal{Y}^{3/2}/\mathcal{Q}$, where $c_{i}$ ($i=1,2,3\cdots$) is dimensionless constants. This case gives nearly constant $a_{\textsc{mond}}$ in the late-time Universe (see Fig. 1). Replacing $\mathcal{Q}$ with $\mathcal{Q}-\mathcal{Q}_{c}$ in the denominator, we obtain $a_{\textsc{mond}}\propto(1+z)^{3}$ based on Eq. (6a). Figure 2 depicts this result together with the Magneticum result [32]. We see that this simple case fails to accurately describe the Magneticum trend. The best polynomial fit of the Magneticum result is $a_{\textsc{mond}}=[0.9+0.2(1+z)^{2}]\times 10^{-10}\,{\rm m}/{\rm s}^{2}$ [63]. Therefore, in the Skordis-Złośnik MOND theory, the Magneticum trend could be reproduced by

[TABLE]

For the parameters used in Sec. II, we obtain $c_{1}=4.6\times 10^{-5}$ and $c_{2}=1.3\times 10^{6}$ for the best fit. Note that Eq. (11) gives an exactly constant $G_{\textsc{n}}$. Furthermore, if $\mathcal{J}\propto\mathcal{Y}^{3/2}$, then its specific form does not affect the cosmological linear perturbation analysis of the theory. The reason is that an equation similar to Eq. (7b) can be obtained in the case of expanding Universe.

Equation (11) could reproduce the Magneticum trend, but may not be the most natural way — the functional form and relevant parameters require slight fine-tuning. Figure 2 also plots the case of $a_{\textsc{mond}}\propto cH$, which is pretty close to the Magneticum result. Inspired by the numerical coincidence between $a_{\textsc{mond},0}$ and $cH_{0}$, Milgrom [1] firstly conjectured this relation. A theory that $a_{\textsc{mond}}$ is controlled by $cH$ multiplied by a $\mathcal{O}(1)$-valued redshift-dependent function seems more natural. The Skordis-Złośnik theory may not be able to realize $a_{\textsc{mond}}\propto cH$. The reason is that the explicit time-dependent variable that exists here is $\mathcal{Q}-\mathcal{Q}_{c}\propto(1+z)^{3}$, rather than an explicit $H$-dependent expression. Luckily, a minor extension of the Skordis-Złośnik theory can achieve the desired scenario. Considering $\nabla_{\mu}A^{\mu}=3H/c$ [59] (see also the Acknowledgements), we see that replacing $\mathcal{J}(\mathcal{Y},\mathcal{Q})$ with

[TABLE]

gives $a_{\textsc{mond}}=2cH/c_{1}$, where the dimensionless constant $c_{1}\approx 4\pi$ [60]. Note that, in the MOND analysis, we only need to consider the background value of $\nabla_{\mu}A^{\mu}$ because $\mathcal{Y}^{3/2}$ is a high-order infinitesimal. In addition, this extension does not destroy the success of the original Skordis-Złośnik theory in fitting the cosmological observations [17]. This is due to the facts that $\mathcal{Y}=0$ in the cosmological background and $\mathcal{Y}^{3/2}$ only contributes the high-order terms in the cosmological linear perturbation analysis. For the same reason, in the framework of the Skordis-Złośnik theory with minor extension, we can link MOND to dark energy with

[TABLE]

where $V_{\textsc{de}}$ is the field potential of dark energy. Following the conventions in [64], we know $c_{1}$ is dimensionless and $[V_{\textsc{de}}]={\rm length}^{-2}$. There are some work in the literature discussing possible links between MOND and dark energy [60, 61, 62]. Here Eq. (13) provides a new example. If dark energy is the cosmological constant, then this case gives constant $a_{\textsc{mond}}$. However, if dark energy is dynamical, then $a_{\textsc{mond}}$ could be redshift-dependent. In Fig. 2, the green line plots an illustration for the power-law potential with the index $p=1$ [65, 66]. Detailed evolution of dark energy can be found in Appendix A. We emphasize that current observations require $p<0.06$ [67], which in turn gives a flatter $a_{\textsc{mond}}(z)$. Therefore, relativistic theory linking MOND to dark energy is not a good option to reproduce the Magneticum trend.

Table 1 summarizes the models mentioned above. One thing is worth mentioning here. Generally, in the TeVeS theory [10], the $a_{\textsc{mond}}$ changes much more slowly than $cH$ [33]. To our knowledge, no specific TeVeS theory has been confirmed that can realize $a_{\textsc{mond}}\propto cH$.

IV Gravitational wave analysis

Gravitational wave properties in the Skordis-Złośnik theory can be determined by solving the linearized equations of motion about the Minkowski spacetime defined by $\bar{g}_{\mu\nu}=\eta_{\mu\nu}$, $\bar{A}^{\mu}=\delta^{\mu}_{0}$, and $\bar{\phi}={\rm const.}$. This background solution requires that $\bar{\mathcal{F}}=\bar{\mathcal{F}}(\bar{g}_{\mu\nu},\bar{A}^{\mu},\bar{\phi})=0$. The perturbed solutions are $g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$, $A^{\mu}=\delta^{\mu}_{0}+a^{\mu}$, and $\phi=\bar{\phi}+\varphi$. Since the linearized equations of motion are very complicated and coupled together, it is easier to use the gauge-invariant formalism to decouple the equations [68, 69]. Following Gong et al. [69], one can decompose the components of $h_{\mu\nu}$ and $a^{\mu}$ as

[TABLE]

Here, $\partial^{k}h_{jk}^{\mathrm{TT}}=0$, $\eta^{jk}h_{jk}^{\mathrm{TT}}=0$, and $\partial_{j}\beta^{j}=\partial_{j}\epsilon^{j}=\partial_{j}\mu^{j}=0$. Under the infinitesimal coordinate transformation parameterized by $\xi^{\mu}=(\xi^{t},\xi^{j})=(A,B^{j}+\partial^{j}C)$ with $\partial_{j}B^{j}=0$, one knows that

[TABLE]

Therefore, one determines the following gauge-invariant variables,

[TABLE]

Then, one can try to reexpress the linearized equations of motion to conclude that

[TABLE]

assuming $\partial\bar{\mathcal{F}}/\partial\mathcal{Q}=0$, where barred quantities are to be evaluated at the flat spacetime background. The above equations show that the tensor and vector modes are propagating at the speed of light, while $\varphi$ generally travels at a different speed, which is smaller than the speed of light for the parameter values adopted in Sec. II. This result reconfirmed the conclusion presented in [17].

Provided that the ordinary matter couples with the metric minimally, one can calculate the geodesic deviation equation, $\ddot{x}^{j}=-R_{tjtk}x^{k}$, to determine the polarizations of gravitational waves [70]. It turns of that

[TABLE]

Therefore, there are only two polarizations (plus and cross), like in general relativity.

V Conclusions

In this paper, we discuss the Newtonian, MOND and gravitational wave analyses for the Skordis-Złośnik theory [17]. In the first two cases, after abandoning the static assumption adopted in [17], we find that whether $G_{\textsc{n}}$ and $a_{\textsc{mond}}$ are time-varying depends on the specific form of the $\mathcal{J}$ function of the theory. Screening the scalar field in strong field region is a sufficient condition to give a constant $G_{\textsc{n}}$, and this scenario may be a preferred choice in both theory and observations. For the $a_{\textsc{mond}}$, we highlight that the theory with Eq. (11) could reproduce the $a_{\textsc{mond}}$–$z$ dependence observed in the Magneticum simulations [32]. Minor extension of the original Skordis-Złośnik theory with Eq. (12) gives $a_{\textsc{mond}}\propto cH$. For the gravitational wave analysis, we show that there are only two tensor polarizations, which is preferred by the GW170814 observations [35, 36].

Acknowledgements

We especially thank the referee for pointing out $\nabla_{\mu}A^{\mu}\propto H(z)$ and suggesting that we discuss the model described by Eq. (12). This work was supported by the National Natural Science Foundation of China under Grants No. 11633001, No. 11920101003, No. 12021003 and No. 11690023, and the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000. S. T. was supported by the Initiative Postdocs Supporting Program under Grant No. BX20200065 and China Postdoctoral Science Foundation under Grant No. 2021M700481. S. H. was supported by the National Natural Science Foundation of China under Grant No. 12205222.

Appendix A Cosmic evolution of dark energy

Equation (13) describes a model linking MOND to dark energy. We adopt a quintessence model with field potential $V_{\textsc{de}}(\phi_{\textsc{de}})=V_{\textsc{de},0}\cdot(\phi_{\textsc{de},0}/\phi_{\textsc{de}})^{p}$ [65, 66], where the index $p\geq 0$. Following the conventions in [64], we have $[V_{\textsc{de}}]={\rm length}^{-2}$. For the flat FLRW Universe, the cosmic evolution equations are [64]

[TABLE]

where ${}^{\prime}\equiv{\rm d}/{\rm d}\phi_{\textsc{de}}$, $\rho_{\textsc{de}}=(c^{2}/8\pi\tilde{G})\cdot[\dot{\phi}_{\textsc{de}}^{2}/(2c^{2})+V_{\textsc{de}}]$ and the subscript f means fluid. Here we regard MOND as a pressureless dark matter and include its contribution in $\rho_{\textsc{f}}$. This is a good approximation as shown in Fig. 1. Then the equation of state $w_{\textsc{f}}$ is given by Eq. (4) in [71]. Introducing the dimensionless variables

[TABLE]

the above evolution equations can be rewritten as

[TABLE]

where $L=(1-w_{\textsc{f}})x_{1}^{2}+(1+w_{\textsc{f}})(1-x_{2}^{2})$ and $N=\ln(a/a_{0})$. The relative dark energy density $\Omega_{\textsc{de}}=x_{1}^{2}+x_{2}^{2}$. Figure 3 presents the numerical solutions of Eq. (21), and illustrates the frozen and tracker properties of this model [66]. We can use the tracker solution to calculate $a_{\textsc{mond}}(z)$ in the low-redshift Universe. Especially, we have $a_{\textsc{mond}}\propto\sqrt{V_{\textsc{de}}}\propto\lambda^{p/2}$. Considering the calibration at $z=0$ in Fig. 2, we can directly obtain $a_{\textsc{mond}}(z)$ from the solution of $\lambda$ without the value of $V_{\textsc{de},0}$. Parameters adopted in Fig. 3 are used to plot the green line in Fig. 2.