Bumblebee cosmology: The FLRW solution and the CMB temperature anisotropy

Abstract

We put into test the idea of replacing dark energy by a vector field against the cosmic microwave background (CMB) observation using the simplest vector-tensor theory, where a massive vector field couples to the Ricci scalar and the Ricci tensor quadratically. First, a remarkable Friedmann−Lemaître−Robertson−Walker (FLRW) metric solution that is completely independent of the matter-energy compositions of the universe is found. Second, based on the FLRW solution as well as the perturbation equations, a numerical code calculating the CMB temperature power spectrum is built. We find that though the FLRW solution can mimic the evolution of the universe in the standard \mathrm{\Lambda }CDM model, the calculated CMB temperature power spectrum shows unavoidable discrepancies from the CMB power spectrum measurements.

Graphical abstract

Keywords

cosmological perturbation theory / Bumblebee gravity / Friedmann−Lemaître−Robertson−Walker metric / cosmic microwave background anisotropy

1 Introduction

The observation of the cosmic microwave background (CMB) anisotropy has promoted cosmology to a high-precision scientific field of research [1, 2]. Based on the standard \mathrm{\Lambda }CDM model in general relativity (GR), CMB data from the Planck mission updated our understanding of the geometry and the compositions of the universe [3]. Moreover, the precisely measured CMB anisotropy provides stringent tests for alternative theories of gravity, helping us to judge between dark energy and its substitutions (e.g., see Refs. [4-8]). Motivated by examing the idea of replacing dark energy by an auxiliary vector field rigorously, we build a code to calculate the power spectrum of the CMB temperature anisotropy in an appealing vector−tensor theory, called the bumblebee gravity [9], to compare with the standard \mathrm{\Lambda }CDM result.

The bumblebee gravity theory has the action

(1)\begin{array}{rl}S=& \frac{1}{2\kappa }\int {\mathrm{d}}^{4}x\sqrt{-g}\left(R+{\xi }_1{B}^{\mu }{B}^{\nu }{R}_{\mu \nu }+{\xi }_2{B}^{\mu }{B}_{\mu }R\right)\\ & -\int {\mathrm{d}}^{4}x\sqrt{-g}\left(\frac{1}{4}{B}^{\mu \nu }{B}_{\mu \nu }+V\right),\end{array}

where B^{\mu } is the auxiliary vector field called the bumblebee field, and B_{\mu \nu }:=D_{\mu }{B}_{\nu }-D_{\nu }{B}_{\mu } with D_{\mu } being the covariant derivative. The theory possesses the two simplest nonminimal coupling terms between the vector field and the curvature quantities, namely the Ricci tensor R_{\mu \nu } and the Ricci scalar R, with {\xi }_1 and {\xi }_2 being the coupling constants. The potential V is generally a function of B^{\mu }. In this work, it takes the simplest form

(2)V=V_1{B}^{\mu }{B}_{\mu },

where V_1 is a positive constant and represents the square of the effective mass of the bumblebee field. We will mostly use the geometrized units where the gravitational constant G and the speed of light c are set to unity (\kappa =8\pi), though physical quantities are sometimes given in conventional units for a better perception of their sizes.

A theory similar to the one in Eq. (1) but without the potential V was originally proposed and studied by Hellings and Nordtvedt [10]. They calculated coefficients for the parametrized post-Newtonian (PPN) formula and the Friedmann−Lemaître−Robertson−Walker (FLRW) metric solution in the theory. Later, Kostelecký started to use this theory to introduce spontaneous breaking of local Lorentz symmetry in gravity by adding in the potential V and requiring its minimum to be realized with a nonzero configuration of the bumblebee field [9]. The theory then received intensive studies in seeking for features and effects of Lorentz-symmetry violation in gravity (e.g., see Refs. [11-17]). Recently, interests in this theory as the simplest vector−tensor theory have returned since two branches of solutions of spherical black holes (BHs) were found for the special case of V=0 and {\xi }_2=0 [18, 19]. One branch of the solutions generalizes the Reissner−Nordström (RN) BH while the other is a deformation of the Schwarzschild BH due to the accompanied nontrival bumblebee field. Based on the generalized RN BH solutions, the first law of BH thermodynamics has been extended and numerically checked [20], and interesting theoretic bounds on the charge of the BHs have been discovered by investigating the dynamic instabilities of the solutions [21]. We point out that the potential which we use in Eq. (2) simply has its minimum when the bumblebee field vanishes. So it generates no vacuum expectation value of the bumblebee field and cannot cause spontaneous spacetime symmetry breaking. This is usually different in the literature where the potential takes other forms and has its minimum when the bumblebee field is shifted from zero (e.g., see Ref. [9]); the nonzero vacuum expectation value of the bumblebee field causes spontaneous spacetime symmetry breaking.

Appealed by the intriguing BH solutions in the bumblebee theory, we continue the work of Hellings and Nordtvedt by extending their FLRW solution to the case of a massive vector field and carrying out the calculation of the power spectrum of the CMB temperature anisotropy. We find that although the FLRW solution can mimic the expansion of the universe in the standard \mathrm{\Lambda }CDM scenario, the CMB power spectrum deviates radically from the standard \mathrm{\Lambda }CDM result at very large scales, disfavoring substituting dark energy with the bumblebee field. It demonstrates how effective the modern CMB observation can be in ruling out the alternative gravity theory that has escaped tests from weak-field observations [10, 18, 19], BH images [19, 22], and redshift-distance measurements in cosmology [23]. In building our CMB code, we use CAMB [24] and CLASS [25] as important references. Compared to modifications based on them, such as MGCAMB [26, 27], we have to derive the approximate solutions at the early universe according to the field equations in the bumblebee theory to set up proper initial conditions for the numerical integration. It is the first public CMB code working for an action-based modified gravity to our knowledge. The code is available at github.com/ryxxastroat/bumblebeecmb.

2 FLRW solution in the bumblebee theory

We start with the FLRW metric ansatz

(3)\mathrm{d}{s}^2=a^2\left(-\mathrm{d}{\eta }^2+{\displaystyle \frac{\mathrm{d}{r}^2}{1-{\mathcal{K}}_0{r}^2}}+r^2\mathrm{d}{\mathrm{\Omega }}^2\right),

where a is the time-depending scale factor and the constant {\mathcal{K}}_0 represents the current spatial curvature of the universe. The homogeneous and isotropic FLRW metric requires the bumblebee field to take the form of

(4)B_{\mu }=\left(b_{\eta },0,0,0\right),

where b_{\eta } depends only on the time coordinate \eta. The energy−momentum tensors for matter and radiation take the usual perfect-fluid form,

(5)\begin{array}{rl}& {\left(T_m\right)}_{\mu \nu }=\left({ϵ}_m+p_m\right){\left(u_m\right)}_{\mu }{\left(u_m\right)}_{\nu }+p_m{g}_{\mu \nu },\\ & {\left(T_r\right)}_{\mu \nu }=\left({ϵ}_r+p_r\right){\left(u_r\right)}_{\mu }{\left(u_r\right)}_{\nu }+p_r{g}_{\mu \nu },\end{array}

where {ϵ}_A,p_A,{\left(u_A\right)}_{\mu }, A=m,r, are energy densities, pressures, and four-velocities for matter and radiation. The equations of state for matter and radiation are

(6)\begin{array}{rl}& p_m=0,\\ & p_r=\frac{1}{3}{ϵ}_r.\end{array}

The four-velocities are

(7){\left(u_m\right)}_{\mu }={\left(u_r\right)}_{\mu }=\left(-a,0,0,0\right),

to be consistent with the FLRW metric.

With the above setup and using the field equations presented in Appendix A, we find two ordinary differential equations (ODEs) for a and b_{\eta },

(8)\begin{array}{rl}0=& -\frac{3({\xi }_1+2{\xi }_2)a^{\prime }{b}_{\eta }^{\prime }{b}_{\eta }}{a}+\left(\frac{3{\xi }_2{a}^{\mathrm{\prime }2}}{a^2}+\kappa V_1{a}^2-3{\mathcal{K}}_0{\xi }_2\right)b_{\eta }^2\\ & +3a^{\mathrm{\prime }2}-\kappa \left({ϵ}_m+{ϵ}_r\right)a^4+3{\mathcal{K}}_0{a}^2,\\ 0=& b_{\eta }\left[3({\xi }_1+2{\xi }_2)a^{″}-3{\xi }_1\frac{a^{\mathrm{\prime }2}}{a}+6{\xi }_2{\mathcal{K}}_{0}a-2\kappa V_1{a}^3\right],\end{array}

where the primes denote derivatives with respect to the conformal time \eta. The first equation is derived from the Einstein field equations in Eq. (A.1), and the second equation is from the temporal component of the vector field equation in Eq. (A.1). It is also possible to derive the second equation in Eq. (8) from the Einstein field equations, but it is less straightforward than using the temporal component of the bumblebee field equation.

Equation (8) interestingly admits two solutions. The first one, with b_{\eta }=0, simply reduces to the GR solution without dark energy. The second one, with b_{\eta }\ne 0, is a nontrivial new FLRW solution and remarkably has an elegant expression for the expansion rate,

(9)\mathcal{H}:={\displaystyle \frac{a^{\prime }}{a}}=H_0\sqrt{{\mathrm{\Omega }}_{V_1}{a}^2+\left(1-{\mathrm{\Omega }}_{V_1}-{\mathrm{\Omega }}_{\mathcal{K}0}\right)a^{\alpha }+{\mathrm{\Omega }}_{\mathcal{K}0}},

where H_0 is the Hubble constant, and

(10)\begin{array}{rl}& {\mathrm{\Omega }}_{V_1}=\frac{2{\tilde{V}}_1}{{\xi }_1+4{\xi }_2},\\ & {\mathrm{\Omega }}_{\mathcal{K}0}=-\frac{{\mathcal{K}}_0}{H_0^2},\\ & \alpha =-\frac{4{\xi }_2}{{\xi }_1+2{\xi }_2},\end{array}

with {\tilde{V}}_1=V_1/{ϵ}_{\mathrm{c}\mathrm{r}\mathrm{i}0} and {ϵ}_{\mathrm{c}\mathrm{r}\mathrm{i}0}=3H_0^2/\kappa being the current critical energy density of the universe. At the background level, the energy−momentum tensors for matter and radiation are conserved separately so that

(11)\begin{array}{rl}& {ϵ}_m=\frac{{ϵ}_{m0}}{a^3},\\ & {ϵ}_r=\frac{{ϵ}_{r0}}{a^4},\end{array}

where {ϵ}_{m0} and {ϵ}_{r0} are the current energy densities of matter and radiation, respectively. The evolution of b_{\eta } can then be solved numerically using the first equation in Eq. (8).

The new bumblebee FLRW solution in Eq. (9) resembles the standard \mathrm{\Lambda }CDM solution

(12){\mathcal{H}}_{\mathrm{G}\mathrm{R}}=H_0\sqrt{{\mathrm{\Omega }}_{\mathrm{\Lambda }0}{a}^2+{\mathrm{\Omega }}_{r0}{a}^{-2}+{\mathrm{\Omega }}_{m0}{a}^{-1}+{\mathrm{\Omega }}_{\mathcal{K}0}},

where {\mathrm{\Omega }}_{\mathrm{\Lambda }0},{\mathrm{\Omega }}_{m0} and {\mathrm{\Omega }}_{r0} are the current energy fractions of dark energy, matter and radiation. But there are two important aspects where they differ.

　1) There is no cosmological constant in setting up the bumblebee theory; the effect of V_1 turns out to mimic the cosmological constant.

　2) Matter and radiation play no role in Eq. (9); they are replaced by a new term proportional to a^{\alpha }, where \alpha depends on the ratio between the two coupling constants.

The first point is appealing because it provides a possible origin for the cosmological constant or dark energy. The second point, though seems to contradict the doctrine that matter must influence spacetime, one finds that the nonminimal couplings between the bumblebee field and the curvature quantities contribute a negative effective energy density that cancel out the energy densities of matter and radiation. So the expasion rate of the universe depends on the coupling constants {\xi }_1,{\xi }_2, and the effective mass parameter V_1 of the bumblebee field, rather than the energy fractions of matter and radiation.

To make the point clear, the bumblebee cosmological model in Eq. (9) therefore provides an example where the expansion of the universe is a consequence of spacetime itself interacting with an auxiliary vector field, which possibly emerges from an underlying theory. Matter and radiation, presumably the source of spacetime curvature, become guests visiting the prefixed background universe and cannot influence it at the homogeneous and isotropic background level. The idea is so radically different from the conventional \mathrm{\Lambda }CDM scenario, yet there are no a priori reasons to exclude it. Predictions from this new bumblebee cosmological model need to be tested against observations to learn more about its validity.

Tests of the FLRW solution in Eq. (9) using distance-redshift data from selected standard type Ia supernovae and measurements of baryon acoustic oscillations have been done carefully in a separate work [23]. Sensible best-fit values for the 4 parameters, {\xi }_1/{\tilde{V}}_1,{\xi }_2/{\tilde{V}}_1,{\mathrm{\Omega }}_{\mathcal{K}0}, and H_0, have been obtained. In this work, we concentrate on testing the solution independently using the observation of CMB temperature power spectrum. We will take {\mathcal{K}}_0=0, corresponding to a spatially flat cosmological model, to simplify the calculation. For {\mathcal{K}}_0\ne 0, we expect the results to be qualitatively the same. This is because the most contribution to the CMB power spectrum happens at recombination when a\sim {10}^{-3}, so that the behavior of the expansion rate \mathcal{H} at a\to 0 is crucial. If \alpha \ge 0, \mathcal{H} is finite at a\to 0, which leads to pathological behaviors of the CMB perturbation variables (see Appendix C). Therefore, we restrict to , in which case a nonzero {\mathrm{\Omega }}_{\mathcal{K}0} does not affect the a^{\alpha /2} behavior of \mathcal{H} at a\to 0 and thus is not expected to qualitatively change the CMB results that we are going to show. To further reduce the parameter space, we require the solution to satisfy two brief observational conditions: (i) The deceleration parameter q_0=-\mathrm{d}\ln \mathcal{H}/\mathrm{d}\ln a{|}_{a=1} is negative [28]. (ii) The age of the universe is greater than 9.5\times {10}^9 years according to radioactive dating [29]. The resultant parameter space is shown in Fig. 1. In Fig. 2, a few examples of the bumblebee FLRW solution with {\mathrm{\Omega }}_{\mathcal{K}0}=0 are plotted. Note that the parameters q_0 and \alpha are used in place of the parameters {\xi }_1/{\tilde{V}}_1 and {\xi }_2/{\tilde{V}}_1, with the relation being

(13){\displaystyle \frac{{\xi }_1}{{\tilde{V}}_1}}=-{\displaystyle \frac{2\alpha +4}{\alpha +2q_0}},{\displaystyle \frac{{\xi }_2}{{\tilde{V}}_1}}={\displaystyle \frac{\alpha }{\alpha +2q_0}}.

3 CMB temperature anisotropy in the bumblebee theory

The CMB temperature anisotropy is represented by the correlation function C\left(\hat{\mathit{n}}\cdot {\hat{\mathit{n}}}^{\prime }\right):=⟨\mathrm{\Theta }(\hat{\mathit{n}})\mathrm{\Theta }({\hat{\mathit{n}}}^{\prime })⟩, with \mathrm{\Theta }:=\mathrm{\Delta }{T}_{\mathrm{C}\mathrm{M}\mathrm{B}}/T_{\mathrm{C}\mathrm{M}\mathrm{B}} being the fractional fluctuation of the CMB temperature and -\hat{\mathit{n}} representing the unit vector along the received CMB photons’ path. On the observational side, the multipoles of C\left(\hat{\mathit{n}}\cdot {\hat{\mathit{n}}}^{\prime }\right) at the current time \eta ={\eta }_0 and the location of the Earth \mathit{x}=0 are measured in Planck’s mission [30]. On the theoretical side, as presented here for the bumblebee theory and in the literature and textbooks of cosmology for GR (e.g., see Refs. [25, 31-34]), the multipoles of C\left(\hat{\mathit{n}}\cdot {\hat{\mathit{n}}}^{\prime }\right) at \eta ={\eta }_0,\mathit{x}=0 are calculated in the following three steps.

First, a truncated system of differential equations governing the evolution of the multipoles of \mathrm{\Theta } in the Fourier space needs to be solved. These equations involve perturbations in the spacetime metric and the velocity of baryon matter, on top of the background FLRW description of the universe, so they are solved together with the linear perturbation of the Einstein field equations. The complete set of the equations in the bumblebee theory is presented in Appendix B. The appropriate initial condition to solve the set of the equations is derived in Appendix C. For explaining our results, we list here the perturbation variables used.

\begin{array}{rl}& \mathrm{\Psi }:\mathrm{P}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{i}\mathrm{n}{g}_{\eta \eta };\\ & \mathrm{\Phi }:\mathrm{S}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{r}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{i}\mathrm{n}{g}_{ij};\\ & \delta b_{\eta }:\mathrm{P}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}{B}_{\eta };\\ & \delta b_S:\mathrm{S}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{r}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{i}\mathrm{n}{B}_i;\\ & {\delta }_b:\mathrm{F}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y};\\ & {\delta }_c:\mathrm{F}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{k}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y};\\ & v_b:\mathrm{S}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{r}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{b}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{t}\mathrm{y};\\ & v_c:\mathrm{S}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{a}\mathrm{r}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{k}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{i}\mathrm{t}\mathrm{y};\\ & {\mathrm{\Theta }}_l:\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{C}\mathrm{M}\mathrm{B}\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}.\end{array}

Second, the line-of-sight formula [34]

(14)\begin{array}{rl}{\mathrm{\Theta }}_l(k)=& {\int }_0^{{\eta }_0}\mathrm{d}\eta [\left(g\left({\mathrm{\Theta }}_0+\mathrm{\Psi }\right)+{\mathrm{e}}^{-\tau }\left({\mathrm{\Phi }}^{\prime }+{\mathrm{\Psi }}^{\prime }\right)\right)j_l(k\chi )\\ & -g{v}_b\frac{\mathrm{d}}{\mathrm{d}(k\chi )}{j}_l(k\chi )]\end{array}

is used to generate {\mathrm{\Theta }}_l in the Fourier space at \eta ={\eta }_0 for any l. In Eq. (14), j_l are the spherical Bessel functions, k is the magnitude of the wave vector for each Fourier mode, and \chi ={\eta }_0-\eta. The optical depth \tau and the visibility function g depend on the background number density of electrons, which can be calculated using the Peebles equation given the expansion history of the universe [35].

Third, the multipoles of C\left(\hat{\mathit{n}}\cdot {\hat{\mathit{n}}}^{\prime }\right) at \eta ={\eta }_0,\mathit{x}=0 are calculated via

(15)C_l=4\pi {\displaystyle \int {\displaystyle \frac{\mathrm{d}k}{k}}{\mathrm{\Theta }}_l^2(k){\mathrm{\Delta }}_{\mathcal{R}}^2(k),}

where {\mathrm{\Delta }}_{\mathcal{R}}^2(k) is the initial power spectrum describing the size of the initial perturbations for each Fourier mode. Motivated by the inflation theory [36, 37], {\mathrm{\Delta }}_{\mathcal{R}}^2(k) takes the form

(16){\mathrm{\Delta }}_{\mathcal{R}}^2(k)=A_s{k}^{n_s-1},

where the constants A_s and n_s are determined by fitting the theoretically calculated C_l to the observational data.

Following the three steps, we have developed a code using Wolfram Mathematica (solving the equations in step one) and Python (calculating numerical integrals in steps two and three) to compute the multipoles C_l in the bumblebee theory. The code is available at github.com/ryxxastroat/bumblebeecmb. Here we present our numerical results at three levels accordingly: (i) evolution of the perturbation variables, (ii) {\mathrm{\Theta }}_l at \eta ={\eta }_0 as functions of k, and (iii) C_l changing with l.

We start with the final results of C_l in Fig. 3. The plot shows the CMB power spectrum in the bumblebee gravity calculated for representative values of the parameters q_0 and \alpha. An impressive observation is that the parameters q_0 and \alpha each determines one important feature in the power spectrum plot: the peaks of the power spectrum shift to right as q_0 increases while the distances between two consecutive peaks shrink as \alpha increases. This observation makes it possible to adjust q_0 and \alpha to more or less fit the standard \mathrm{\Lambda }CDM result at large l. But there is then the fatal discrepency at small l between the bumblebee results and the standard C_l curve: The unexpected raise of the power spectrum as l goes from about 100 to 1 in the bumblebee results rules out any possibility of a sensible match to the standard CMB power spectrum.

Before shedding light on the cause of the discrepency, let us point out that the results are found to be independent of the parameter V_1 as long as {\xi }_1/{\tilde{V}}_1 and {\xi }_2/{\tilde{V}}_1 are fixed. So q_0 and \alpha, equivalently {\xi }_1/{\tilde{V}}_1 and {\xi }_2/{\tilde{V}}_1, are the only two parameters to adjust in our model. There are the conventional parameters H_0,T_{\mathrm{C}\mathrm{M}\mathrm{B}}, and the fraction of baryonic matter {\mathrm{\Omega }}_{b0}, but they are more or less fixed by observations [38-40]. Small changes of them alter the results insignificantly. Then there is the fraction of dark matter {\mathrm{\Omega }}_{c0}, but it turns out to have no significant effect at all. In fact, we have set it to zero. The initial power spectrum index n_S changes the overall tilt of the calculated power spectrum in a predictable way, so the approximation n_s=1 has been used.

Now let us look at the middle-level results, {\mathrm{\Theta }}_l at \eta ={\eta }_0 as functions of k. Typical examples of {\mathrm{\Theta }}_l are presented in Fig. 4. They demonstrate that for large l, {\mathrm{\Theta }}_l in the bumblebee model and in the standard \mathrm{\Lambda }CDM model have comparable sizes, while for small l, {\mathrm{\Theta }}_l in the bumblebee model are visibly larger. Let us point out that the global maxima of {\mathrm{\Theta }}_l in Fig. 4 appear at k\sim l/{\eta }_0 with {\eta }_0\sim 3H_0^{-1}, because that is where the integrand in Eq. (14) contributes most. To be specific, the spherical Bessel functions j_l(x) and their derivatives have global maxima at x\sim l, and the visibility function g has a sharp peak at the epoch of recombination when the conformal time is \eta ={\eta }_{\ast }\approx 0. So they together pick up the integrand at k\sim l/({\eta }_0-{\eta }_{\ast })\sim l/{\eta }_0. Note that the term {\mathrm{e}}^{-\tau }\left({\mathrm{\Phi }}^{\prime }+{\mathrm{\Psi }}^{\prime }\right) is much smaller than the other terms, so it can be neglected when undertanding the qualitative features in the results of {\mathrm{\Theta }}_l.

Lastly, to explain why {\mathrm{\Theta }}_l at small l in the bumblebee model are larger than those in the standard \mathrm{\Lambda }CDM model, we come to the evolution of the perturbation variables. In Fig. 5, typical solutions of the perturbation variables used in Eq. (14) to calculate {\mathrm{\Theta }}_l are shown. The biggest difference between the bumblebee results and the standard \mathrm{\Lambda }CDM results is the evolution of the metric perturbation variable \mathrm{\Psi } (\mathrm{\Phi } has a similar behavior to \mathrm{\Psi }). In the bumblebee model, \mathrm{\Psi } is incredibly small. This is obtained from numerically solving the complete set of the equations governing the evolution of the perturbation variables. The key of the numerical integration of the system is finding out the initial conditions. The detailed derivation is presented in Appendix C. There we find that the modified Einstein field equations accompanied with the vector field equation in the bumblebee theory lead to \mathrm{\Psi }\sim \mathrm{\Phi }\sim a^{-\alpha /2} at a\to 0. It is vastly different from the GR case where \mathrm{\Psi } and \mathrm{\Phi } are approximately constant at a\to 0. As \alpha is negative under our consideration, \mathrm{\Psi } and \mathrm{\Phi } are condemned to be small.

So we can explain the size discrepancy between {\mathrm{\Theta }}_l in the bumblebee model and in the standard \mathrm{\Lambda }CDM model at small l now. Neglecting the relatively small term {\mathrm{e}}^{-\tau }\left({\mathrm{\Psi }}^{\prime }+{\mathrm{\Phi }}^{\prime }\right), the main contribution of the integrand in Eq. (14) comes from {\mathrm{\Theta }}_0+\mathrm{\Psi } and v_b around the recombination time \eta ={\eta }_{\ast } (\ln a\sim -7 at recombination) where the visibility function g peaks sharply. From Fig. 5, the magnitudes of {\mathrm{\Theta }}_0 and v_b in the bumblebee model and in the standard \mathrm{\Lambda }CDM model are comparable at recombination. But in the bumblebee model, \mathrm{\Psi } is extremely small so {\mathrm{\Theta }}_0+\mathrm{\Psi }\approx {\mathrm{\Theta }}_0. In the standard \mathrm{\Lambda }CDM model, the contribution from \mathrm{\Psi } at recombination is significant for small-k Fourier modes. It cancels {\mathrm{\Theta }}_0 so that the standard \mathrm{\Lambda }CDM model has smaller {\mathrm{\Theta }}_0+\mathrm{\Psi } at recombination compared to the bumblebee model. The smaller k is, the more significant for the cancellation. Multiplied by the spherical Bessel functions and their derivatives, the global maximum of the integrand in Eq. (14) satisfies l\sim k{\eta }_0. So small-k Fourier modes are correspondingly the dominant contribution for {\mathrm{\Theta }}_l with small l. Therefore, {\mathrm{\Theta }}_l in the standard \mathrm{\Lambda }CDM model have smaller magnitudes than those in the bumblebee model for small l.

4 Conclusions

First, an interesting FLRW solution in the bumblebee gravitational theory has been found. The solution depicts the expansion of the universe completely through the coupling constants and the mass of the bumblebee vector field, while the compositions of the universe, namely matter and radiation, play no role in the FLRW solution. This makes the bumblebee cosmology conceptually distinctive from the \mathrm{\Lambda }CDM cosmology where it is the compositions of the universe that determine the expansion rate. Despite the profound difference from the standard cosmological model, a decently large parameter space in the bumblebee theory has been found to survive from the two preliminary requirements from astrophysical observations: (i) currently the expansion of the universe is accelerating [28], and (ii) the age of the universe needs to be at least greater than 9.5\times {10}^9 years [29].

Then, to further explore the possibility for the bumblebee cosmology to replace the standard \mathrm{\Lambda }CDM model, we have calculated the power spectrum of the CMB temperature in the bumblebee cosmological model. We find that when l approaches zero, C_l calculated in the bumblebee cosmological model unavoidably increases as l decreases as shown in Fig. 3, which forbids any sensible match between the bumblebee C_l curves and the standard \mathrm{\Lambda }CDM one. Tracing this disastrous discrepancy between the bumblebee cosmological model and the standard \mathrm{\Lambda }CDM model for C_l at small l in each step of our calculation, we find that it is ultimately because the metric perturbation is extremely small in the bumblebee cosmological model. While the metric perturbation is sourced by the perturbations of the matter density and the CMB photon density in GR, the existence of the perturbation of the bumblebee field in the bumblebee theory seems to cause an effectively negative energy density that balances the perturbations of the matter density and the CMB photon density to an extraordinary accuracy, so that the total source for the metric perturbation almost vanishes in the bumblebee theory.

By building our own CMB code from scratch and analyzing the calculated results, we have demonstrated how effective the CMB observation can be in testing action-based gravitational theories. Our code contains only the minimal ingredients necessary for calculating the CMB power spectrum by sacrificing the precison (about 10\mathrm{\%} relative error in C_l, estimated by applying our code to the \mathrm{\Lambda }CDM model; see Appendix D) and the speed. So it is much more transparent to be adapted for other modified gravity theories, compared to CAMB [24] and CLASS [25]. In addition, existing CMB code dealing with modified gravity theories, such as MGCAMB [26, 27], requires the theories to admit parametrizations to fit in its framework, which is usually not the case for action-based gravity theories. Our code then serves as an adequate starting point to get a grasp on the CMB power spectrum in numerous such theories, providing quick tailored CMB tests against them complementary to other celebrated tests, such as the solar-system weak-field observations, pulsar timing and pulse profiles, gravitational waves from compact binary coalescences, and black hole images.

In conclusion, the bumblebee cosmological model studied in this work is excluded by the CMB observation, though it has an interesting FLRW solution that does not depend on the compositions of the universe and needs neither dark energy nor dark matter. We point out that a generalization of our model by adding in spatial components of the bumblebee field and replacing the FLRW metric with the Bianchi I model might cure the issue for the small multipoles encountered in our work. Maluf and Neves [41] have studied the case for a slightly different bumblebee theory. For our theory, we recently notice that De Felice and Hell [42] have carried out a detailed mode analysis for cosmological perturbations. A comprehensive study on the behavior of the CMB power spectrum in the Bianchi I spacetime and the development of the CMB code in this generalized spacetime are interesting but challenging directions for future investigation.

5 Appendix A: Field equations in the bumblebee theory

The field equations given by the action in Eq. (1) are

(A.1)\begin{array}{c}\begin{array}{rl}& G_{\mu \nu }=\kappa T_{\mu \nu }+\kappa {\left(T_B\right)}_{\mu \nu }+{\xi }_1{\left(T_{B1}\right)}_{\mu \nu }+{\xi }_2{\left(T_{B2}\right)}_{\mu \nu },\\ & D^{\mu }{B}_{\mu \nu }-\frac{\mathrm{d}V}{\mathrm{d}{B}^{\nu }}+\frac{{\xi }_1}{\kappa }{B}^{\mu }{R}_{\mu \nu }+\frac{{\xi }_2}{\kappa }{B}_{\nu }R=0,\end{array}\end{array}

where T_{\mu \nu } consists of usual energy-momentum tensors for matter and radiation, the energy-momentum tensor for the bumblebee vector field is

(A.2)\begin{array}{c}\begin{array}{rl}{\left(T_B\right)}_{\mu \nu }=& B_{\mu \lambda }{B}_{\nu }^{\phantom{\nu }\lambda }-g_{\mu \nu }\left(\frac{1}{4}{B}^{\alpha \beta }{B}_{\alpha \beta }+V\right)\\ & +2B_{\mu }{B}_{\nu }\frac{\mathrm{d}V}{\mathrm{d}(B^{\lambda }{B}_{\lambda })},\end{array}\end{array}

and the contributions due to the couplings between the bumblebee field and the spacetime curvature are

(A.3)\begin{array}{c}\begin{array}{rl}{\left(T_{B1}\right)}_{\mu \nu }=& \frac{1}{2}{g}_{\mu \nu }{B}^{\alpha }{B}^{\beta }{R}_{\alpha \beta }-B_{\mu }{B}_{\lambda }{R}_{\nu }^{\phantom{\nu }\lambda }-B_{\nu }{B}_{\lambda }{R}_{\mu }^{\phantom{\mu }\lambda }\\ & +\frac{1}{2}[D_{\kappa }{D}_{\mu }\left(B^{\kappa }{B}_{\nu }\right)+D_{\kappa }{D}_{\nu }\left(B_{\mu }{B}^{\kappa }\right)\\ & -{◻}_g\left(B_{\mu }{B}_{\nu }\right)-g_{\mu \nu }{D}_{\alpha }{D}_{\beta }\left(B^{\alpha }{B}^{\beta }\right)],\\ {\left(T_{B2}\right)}_{\mu \nu }=& -B^{\lambda }{B}_{\lambda }{G}_{\mu \nu }-B_{\mu }{B}_{\nu }R+D_{\mu }{D}_{\nu }\left(B^{\lambda }{B}_{\lambda }\right)\\ & -g_{\mu \nu }{◻}_g\left(B^{\lambda }{B}_{\lambda }\right),\end{array}\end{array}

with {◻}_g=D_{\mu }{D}^{\mu }.

Besides the field equations, the energy−momentum conservation equation for matter and radiation

(A.4)\begin{array}{c}{D}^{\mu }{T}_{\mu \nu }=0\end{array}

is useful in calculations. Note that Eq. (A.4) can be derived from the field equations. In fact, the covariant divergences of {\left(T_B\right)}_{\mu \nu },{\left(T_{B1}\right)}_{\mu \nu } and {\left(T_{B2}\right)}_{\mu \nu } are found to be

(A.5)\begin{array}{c}\begin{array}{rl}& D^{\mu }{\left(T_B\right)}_{\mu \nu }=\left(B_{\nu }^{\phantom{\nu }\lambda }-B_{\nu }{D}^{\lambda }\right)\left(D^{\kappa }{B}_{\kappa \lambda }-\frac{\mathrm{d}V}{\mathrm{d}{B}^{\lambda }}\right),\\ & D^{\mu }{\left(T_{B1}\right)}_{\mu \nu }=\left(B_{\nu }^{\phantom{\nu }\lambda }-B_{\nu }{D}^{\lambda }\right)\left(B^{\kappa }{R}_{\lambda \kappa }\right),\\ & D^{\mu }{\left(T_{B2}\right)}_{\mu \nu }=\left(B_{\nu }^{\phantom{\nu }\lambda }-B_{\nu }{D}^{\lambda }\right)\left(B_{\lambda }R\right),\end{array}\end{array}

so that the conservation of the total of them is guaranteed by the vector field equation in Eq. (A.1). Together with the identity D^{\mu }{G}_{\mu \nu }=0, Eq. (A.4) is just a consequence of the Einstein field equations.

6 Appendix B: Linear perturbation equations in the bumblebee theory

The setup and equations of the linear perturbation theory used for the standard \mathrm{\Lambda }CDM cosmology can be found in the literature (e.g., see Refs. [31, 34]). For using with the bumblebee theory, we need to add in the bumblebee field and derive the perturbation equations from Eq. (A.1). To begin with, the following provides notations of the scalar-vector-tensor (SVT) decompositions of the perturbation variables used in our CMB code. Note that we restrict to the spatially flat FLRW background solution.

First, the metric in the conformal coordinates (\eta ,x,y,z) is

(B.1)\begin{array}{c}{g}_{\mu \nu }=a^2({\eta }_{\mu \nu }+h_{\mu \nu }),\end{array}

where h_{\mu \nu } is the perturbation on top of a Minkowski metric {\eta }_{\mu \nu }. The SVT decompositions of h_{\eta i} and h_{ij} are

(B.2)\begin{array}{c}\begin{array}{rl}{h}_{\eta i}=& {\mathrm{\partial }}_i{h}_S+{\left(h_V\right)}_i,\\ h_{ij}=& {\left(h_{TT}\right)}_{ij}+{\mathrm{\partial }}_i{\left(h_T\right)}_j+{\mathrm{\partial }}_j{\left(h_T\right)}_i\\ & +\frac{1}{2}\left({\delta }_{ij}{\mathrm{\nabla }}^2-{\mathrm{\partial }}_i{\mathrm{\partial }}_j\right)h_1\\ & +\frac{1}{2}\left(3{\mathrm{\partial }}_i{\mathrm{\partial }}_j-{\delta }_{ij}{\mathrm{\nabla }}^2\right)h_2,\end{array}\end{array}

where {\mathrm{\nabla }}^2={\delta }^{kl}{\mathrm{\partial }}_k{\mathrm{\partial }}_l. The vector parts {\left(h_V\right)}_i,{\left(h_T\right)}_i and the tensor part {\left(h_{TT}\right)}_{ij} satisfy the constraints

(B.3)\begin{array}{c}\begin{array}{rl}& {\mathrm{\partial }}_i{\left(h_V\right)}_i={\mathrm{\partial }}_i{\left(h_T\right)}_i={\mathrm{\partial }}_j{\left(h_{TT}\right)}_{ij}=0,\\ & {\delta }^{ij}{\left(h_{TT}\right)}_{ij}=0.\end{array}\end{array}

Second, the conventional energy-momentum tensor T_{\mu \nu } in our bumblebee cosmology model includes three ingredients,

(B.4)\begin{array}{c}\begin{array}{rl}{\left(T_b\right)}_{\mu \nu }=& \left({ϵ}_b+\delta {ϵ}_b\right){\left(u_b\right)}_{\mu }{\left(u_b\right)}_{\nu },\\ {\left(T_c\right)}_{\mu \nu }=& \left({ϵ}_c+\delta {ϵ}_c\right){\left(u_c\right)}_{\mu }{\left(u_c\right)}_{\nu },\\ {\left(T_{\gamma }\right)}_{\mu \nu }=& \frac{4}{3}\left({ϵ}_{\gamma }+\delta {ϵ}_{\gamma }\right){\left(u_{\gamma }\right)}_{\mu }{\left(u_{\gamma }\right)}_{\nu }\\ & +\frac{1}{3}\left({ϵ}_{\gamma }+\delta {ϵ}_{\gamma }\right)g_{\mu \nu }+{\mathrm{\Pi }}_{\mu \nu },\end{array}\end{array}

where the subscripts b,c and \gamma represent baryon, cold dark matter and photon. The background energy densities of the ingredients are denoted as {ϵ}_b,{ϵ}_c and {ϵ}_{\gamma } while their perturbations are \delta {ϵ}_b,\delta {ϵ}_c and \delta {ϵ}_{\gamma }. We have used the approximations that matter is pressureless and that photons’ pressure is a third of their energy density. The four-velocities take the form

(B.5)\begin{array}{c}\begin{array}{rl}& {\left(u_b\right)}_{\eta }={\left(u_c\right)}_{\eta }={\left(u_{\gamma }\right)}_{\eta }=-a\left(1-\frac{1}{2}{h}_{\eta \eta }\right),\\ & {\left(u_b\right)}_i={\mathrm{\partial }}_i\delta u_{bS}+{\left(\delta u_{bV}\right)}_i,\\ & {\left(u_c\right)}_i={\mathrm{\partial }}_i\delta u_{cS}+{\left(\delta u_{cV}\right)}_i,\\ & {\left(u_{\gamma }\right)}_i={\mathrm{\partial }}_i\delta u_{\gamma S}+{\left(\delta u_{\gamma V}\right)}_i,\end{array}\end{array}

up to the linear order of the perturbations, with the constraints

(B.6)\begin{array}{c}{\mathrm{\partial }}_i{\left(\delta u_{bV}\right)}_i={\mathrm{\partial }}_i{\left(\delta u_{cV}\right)}_i={\mathrm{\partial }}_i{\left(\delta u_{\gamma V}\right)}_i=0.\end{array}

The last term in (T_{\gamma }{)}_{\mu \nu } is a higher-order contribution to the energy−momentum tensor of the photons beyond the perfect-fluid description due to the CMB temperature fluctuation. It takes the form

(B.7)\begin{array}{c}\begin{array}{rl}{\mathrm{\Pi }}_{\eta i}=& 0,\\ {\mathrm{\Pi }}_{ij}=& {\left({\mathrm{\Pi }}_{TT}\right)}_{ij}+{\mathrm{\partial }}_i{\left({\mathrm{\Pi }}_T\right)}_j+{\mathrm{\partial }}_j{\left({\mathrm{\Pi }}_T\right)}_i\\ & +\frac{1}{2}\left(3{\mathrm{\partial }}_i{\mathrm{\partial }}_j-{\delta }_{ij}{\mathrm{\nabla }}^2\right)\mathrm{\Pi },\end{array}\end{array}

with the constraints

(B.8)\begin{array}{c}\begin{array}{rl}& {\mathrm{\partial }}_i{\left({\mathrm{\Pi }}_T\right)}_i={\mathrm{\partial }}_j{\left({\mathrm{\Pi }}_{TT}\right)}_{ij}=0,\\ & {\delta }^{ij}{\left({\mathrm{\Pi }}_{TT}\right)}_{ij}=0.\end{array}\end{array}

Third, the bumblebee field consists of the background part B_{\mu }^{(0)}=\left(b_{\eta },0,0,0\right) and the perturbation part B_{\mu }^{(1)} that has an SVT decomposition

(B.9)\begin{array}{c}\begin{array}{rl}& B_{\eta }^{(1)}=\delta b_{\eta },\\ & B_i^{(1)}={\mathrm{\partial }}_i\delta b_S+{\left(\delta b_V\right)}_i.\end{array}\end{array}

The background bumblebee field has only the temporal component to be consistent with the homogeneous and isotropic FLRW spacetime. The perturbation of the bumblebee field has two scalars \delta b_{\eta } and \delta b_S, and a vector {\left(\delta b_V\right)}_i subject to the constraint

(B.10)\begin{array}{c}{\mathrm{\partial }}_i{\left(\delta b_V\right)}_i=0.\end{array}

Now with the above setup, the field equations in Eq. (A.1) simplify to Eq. (8) at the zeroth order of perturbation with {ϵ}_m={ϵ}_b+{ϵ}_c and {ϵ}_r={ϵ}_{\gamma }; at the linear order of perturbation, they produce three sets of equations for the scalar, the vector, and the tensor variables. Each set of equations does not mix with others. For the calculation of the CMB temperature anisotropy, we only need the set of equations for the scalar variables. Fixing the freedom of the coordinates at the perturbation level by choosing the Newtonian gauge,

(B.11)\begin{array}{c}\begin{array}{rl}& h_{\eta \eta }=-2\mathrm{\Psi },\\ & h_S=0,\\ & \frac{1}{3}{\mathrm{\nabla }}^2{h}_1={\mathrm{\nabla }}^2{h}_2=-2\mathrm{\Phi },\end{array}\end{array}

so that the metric perturbation takes the simple form

(B.12)\begin{array}{c}\mathrm{d}{s}^2=a^2\left[-(1+2\mathrm{\Psi })d{\eta }^2+\left(1-2\mathrm{\Phi }\right){\delta }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j\right].\end{array}

We find 4 equations for the 4 scalar variables \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S,

(B.13)\begin{array}{c}\begin{array}{rl}0=& \left(1-\frac{({\xi }_1+{\xi }_2)b_{\eta }^2}{a^2}\right){\mathrm{\Phi }}^{\prime }-\frac{({\xi }_1+2{\xi }_2)b_{\eta }^2}{2a^2}{\mathrm{\Psi }}^{\prime }\\ & +\left(\mathcal{H}+\frac{3{\xi }_2{b}_{\eta }^2\mathcal{H}}{a^2}-\frac{3({\xi }_1+2{\xi }_2)b_{\eta }{b}_{\eta }^{\prime }}{2a^2}\right)\mathrm{\Psi }\\ & +\frac{({\xi }_1+2{\xi }_2)b_{\eta }}{2a^2}\delta b_{\eta }^{\prime }+\frac{-\mathcal{H}({\xi }_1+6{\xi }_2)b_{\eta }+({\xi }_1+2{\xi }_2)b_{\eta }^{\prime }}{2a^2}\delta b_{\eta }\\ & +\frac{\kappa }{2}a\left({ϵ}_b{u}_{bS}+{ϵ}_c{u}_{cS}+\frac{4}{3}{ϵ}_{\gamma }{u}_{\gamma S}\right),\\ 0=& \mathrm{\Phi }-\mathrm{\Psi }-\frac{{\xi }_2{b}_{\eta }^2}{a^2}(\mathrm{\Phi }+\mathrm{\Psi })+\frac{2{\xi }_2{b}_{\eta }}{a^2}\delta b_{\eta }+\frac{{\xi }_1{b}_{\eta }^{\prime }}{a^2}\delta b_S\\ & +\frac{{\xi }_1{b}_{\eta }}{a^2}\delta b_S^{\prime }-\frac{3}{2}\kappa \mathrm{\Pi },\\ 0=& \left(3\mathcal{H}-\frac{3\mathcal{H}{\xi }_1{b}_{\eta }^2}{2a^2}-\frac{3({\xi }_1+2{\xi }_2)b_{\eta }^{\prime }{b}_{\eta }}{2a^2}\right){\mathrm{\Phi }}^{\prime }\\ & -\frac{3\mathcal{H}({\xi }_1+2{\xi }_2)b_{\eta }^2}{2a^2}{\mathrm{\Psi }}^{\prime }-\left(1-\frac{{\xi }_2{b}_{\eta }^2}{a^2}\right){\mathrm{\nabla }}^2\mathrm{\Phi }\\ & +\frac{3\mathcal{H}({\xi }_1+2{\xi }_2)b_{\eta }}{2a^2}\delta b_{\eta }^{\prime }-\frac{\kappa b_{\eta }}{2a^2}{\mathrm{\nabla }}^2\delta b_S^{\prime }\\ & +\frac{b_{\eta }\left(-6\mathcal{H}({\xi }_1+2{\xi }_2)b_{\eta }^{\prime }+6{\mathcal{H}}^2{\xi }_2{b}_{\eta }+b_{\eta }({\xi }_1+2{\xi }_2){\mathrm{\nabla }}^2\right)}{2a^2}\mathrm{\Psi }\\ & +a^2\kappa \left({ϵ}_b+{ϵ}_c+{ϵ}_{\gamma }\right)\mathrm{\Psi }+\frac{\kappa }{2}{a}^2\left(\delta {ϵ}_b+\delta {ϵ}_c+\delta {ϵ}_{\gamma }\right)\\ & +\frac{b_{\eta }\left(-2\kappa a^2{V}_1-6{\mathcal{H}}^2{\xi }_2+\left(\kappa -{\xi }_1-2{\xi }_2\right){\mathrm{\nabla }}^2\right)}{2a^2}\delta b_{\eta }\\ & +\frac{3\mathcal{H}({\xi }_1+2{\xi }_2)b_{\eta }^{\prime }}{2a^2}\delta b_{\eta }-\frac{\mathcal{H}{\xi }_1{b}_{\eta }{\mathrm{\nabla }}^2}{a^2}\delta b_S,\\ 0=& -\frac{2{\xi }_1\left(3{\mathcal{H}}^2({\xi }_1+4{\xi }_2)-2a^2\kappa V_1\right)\delta b_S}{3\kappa ({\xi }_1+2{\xi }_2)}\\ & +\frac{2{\xi }_1{b}_{\eta }\left(\mathcal{H}\mathrm{\Psi }+{\mathrm{\Phi }}^{\prime }\right)}{\kappa }-\delta b_{\eta }^{\prime }+\delta b_S^{″},\end{array}\end{array}

where the primes denote derivatives with respect to \eta and \mathcal{H}=a^{\prime }/a. The quadratic potential V=V_1{B}^{\mu }{B}_{\mu } has been used.

To complete the set of equations for all the scalar variables, we also need the energy-momentum conservation equations and the Boltzmann equation to describe the evolution of the scalar variables in the matter and photon sector. Because the bumblebee field has no direct interaction with matter and photons, these equations for the perturbation variables of matter and photon are the same as those in GR, which are [34]

\begin{array}{rl}{\delta }_b^{\prime }=& 3{\mathrm{\Phi }}^{\prime }+k{v}_b,\\ v_b^{\prime }=& -k\mathrm{\Psi }-\mathcal{H}{v}_b-\frac{4{ϵ}_{\gamma }}{3{ϵ}_b}\mathrm{\Gamma }\left(3{\mathrm{\Theta }}_1+v_b\right),\end{array}

(B.14)\begin{array}{c}\begin{array}{rl}{\delta }_c^{\prime }=& 3{\mathrm{\Phi }}^{\prime }+k{v}_c,\\ v_c^{\prime }=& -k\mathrm{\Psi }-\mathcal{H}{v}_c,\\ {\mathrm{\Theta }}_0^{\prime }=& {\mathrm{\Phi }}^{\prime }-k{\mathrm{\Theta }}_1,\\ {\mathrm{\Theta }}_1^{\prime }=& \frac{k}{3}{\mathrm{\Theta }}_0+\frac{k}{3}\mathrm{\Psi }-\frac{2k}{3}{\mathrm{\Theta }}_2-\mathrm{\Gamma }\left({\mathrm{\Theta }}_1+\frac{1}{3}{v}_b\right),\\ {\mathrm{\Theta }}_l^{\prime }=& -\frac{k}{2l+1}\left[(l+1){\mathrm{\Theta }}_{l+1}-l{\mathrm{\Theta }}_{l-1}\right]\\ & -\mathrm{\Gamma }\left(1-\frac{1}{10}{\delta }_{l2}\right){\mathrm{\Theta }}_l,l\ge 2,\end{array}\end{array}

where {\mathrm{\Theta }}_l are the multipoles of the CMB temperature fluctuation \mathrm{\Theta }:=\mathrm{\Delta }{T}_{\mathrm{C}\mathrm{M}\mathrm{B}}/T_{\mathrm{C}\mathrm{M}\mathrm{B}}. Note that Eq. (B.14) has been written in the Fourier space with k being the magnitude of the wave vector of the Fourier mode. The scalar variables \delta {ϵ}_b,\delta {ϵ}_c,\delta {ϵ}_{\gamma },\delta u_{bS},\delta u_{cS},\delta u_{\gamma S} and \mathrm{\Pi } have been replaced by a set of variables more convenient to use,

(B.15)\begin{array}{c}\begin{array}{rl}& {\delta }_b=\frac{\delta {ϵ}_b}{{ϵ}_b},\\ & {\delta }_c=\frac{\delta {ϵ}_c}{{ϵ}_c},\\ & {\mathrm{\Theta }}_0=\frac{1}{4}{\delta }_{\gamma }=\frac{1}{4}\frac{\delta {ϵ}_{\gamma }}{{ϵ}_{\gamma }},\\ & v_b=\frac{k}{a}\delta u_{bS},\\ & v_c=\frac{k}{a}\delta u_{cS},\\ & {\mathrm{\Theta }}_1=-\frac{k}{3a}\delta u_{\gamma S},\\ & {\mathrm{\Theta }}_2=\frac{3k^2}{8a^2{ϵ}_{\gamma }}\mathrm{\Pi }.\end{array}\end{array}

The interaction rate between electrons and photons is

(B.16)\begin{array}{c}\mathrm{\Gamma }=a{n}_e{\sigma }_T=a{X}_e{n}_b{\sigma }_T,\end{array}

where n_e and n_b are the number densities of electrons and baryons, and {\sigma }_T is the Thomson scattering cross section. The free electron fraction X_e is calculated using the Saha equation before recombination and the Peebles equation during and after recombination [35].

7 Appendix C: Approximate solutions at a→0

For any given value of k, we want to solve Eq. (B.14) together with the Fourier transform of Eq. (B.13) numerically. The adiabatic initial condition, which is inherited from inflation, sets

(C.1)\begin{array}{c}{\delta }_b\approx {\delta }_c\approx 3{\mathrm{\Theta }}_0,\end{array}

at a\to 0. To start numerical integrations at a\to 0, we still need initial conditions for other variables \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S,v_b,v_c,{\mathrm{\Theta }}_1,{\mathrm{\Theta }}_2,..., namely that we need to find approximate solutions for those variables at a\to 0. To do this, we first simplify the equations by removing the velocity variables as well as the higher-order multipoles {\mathrm{\Theta }}_l. Thus, Eq. (B.14) gives

(C.2)\begin{array}{c}\begin{array}{rl}& {\delta }_b\approx 3\mathrm{\Phi }+A_b,\\ & {\delta }_c\approx 3\mathrm{\Phi }+A_c,\\ & {\mathrm{\Theta }}_0\approx \mathrm{\Phi }+A_0,\end{array}\end{array}

where the integral constants satisfy A_b=A_c=3A_0 to be consistent with Eq. (C.1), while Eq. (B.13) becomes a closed set of equations for \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S. Once approximate solutions for \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S are obtained, they can be used in Eq. (B.14) to find approximate solutions for the velocity variables as well as the higher-order multipoles. Then, we will be able to check if the velocity variables and the higher-order multipoles contribute relevant terms in solving {\delta }_b,{\delta }_c,{\mathrm{\Theta }}_0 and \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S approximately. If they do contribute relevantly and cannot be neglected in the first place, we can correct the obtained approximate solutions accordingly.

With the scheme stated, to find approximate solutions for \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S at a\to 0 from Eq. (B.13), the behaviors of \mathcal{H} and b_{\eta } at a\to 0 are necessary. Equation (9) gives

(C.3)

for the spatially flat solution. Note that \alpha =2 corresponds to {\xi }_1+4{\xi }_2=0, and it turns out in this case

(C.4)\begin{array}{c}\mathcal{H}\approx H_{0}a\sqrt{{\displaystyle \frac{8{\tilde{V}}_1}{{\xi }_1}}\ln a}.\end{array}

The \ln a factor spoils the power-law behavior of \mathcal{H} and subsequently power-law behaviors of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S. We will avoid considering the case of \alpha =2. The behavior of b_{\eta } at a\to 0 turns out to be more delicated than that of \mathcal{H}. By a careful analysis of the first equation in Eq. (8), we find

(C.5)

where D_1 and D_2 are integral constants with no presumed values. For the cases of \alpha =-4,0,2 and 3, we find that \ln a is encountered in the approximate solution of b_{\eta }, ruining the expected power-law behavior. We will not consider these pathological cases.

Neglecting the velocity variables and the higher-order multipoles, and using the approximate solutions in Eqs. (C.2), (C.3) and (C.5), Eq. (B.13) in the Fourier space becomes a closed set of ODEs for \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S. We expect to find approximate solutions for them also in the form of power-law relations with respect to a if the time evolution of these variables is physically sensible. This turns out to be true for . For \alpha >0, we find that the approximate solutions for \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S consist of terms like {\mathrm{e}}^{a^r} with constant values of r, which we take as a signal indicating the cases of \alpha >0 being unphysical.

We find power-law approximate solutions of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S for both the case of and the case of . But we will focus only on the case of , because the case of is excluded by the fact that the universe is currently expanding with an acceleration. This can be shown by using Eq. (13) to substitute {\xi }_1 in the approximate solution of b_{\eta } for in Eq. (C.5). The approximate solution exists when

(C.6)\begin{array}{c}-{\displaystyle \frac{\alpha +2q_0}{{\tilde{V}}_1\left(\alpha +4\right)}}>0.\end{array}

With {\tilde{V}}_1>0 and , the inequality cannot hold if the deceleration parameter q_0 is negative. Numerical integration also verifies that for and , b_{\eta } becomes imaginary when a approaches 0.

The approximate solutions of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta } and \delta b_S that we find for the case of take the form

(C.7)\begin{array}{c}\begin{array}{rl}& \mathrm{\Psi }\approx c_1{a}^r,\\ & \mathrm{\Phi }\approx c_2{a}^r,\\ & \delta b_{\eta }\approx c_3{a}^{r-\frac{\alpha }{4}},\\ & \delta b_S\approx c_4{a}^{r-\frac{3\alpha }{4}},\end{array}\end{array}

where c_1,c_2,c_3,c_4 and r are constants. Notice that (i) \mathrm{\Psi } can be eliminated using the second equation in Eq. (B.13), leading to two first-order ODEs for \mathrm{\Phi } and \delta b_{\eta } and one second-order ODE for \delta b_S in the Fourier space, and (ii) the ODEs are inhomogeneous as substituting \delta {ϵ}_b,\delta {ϵ}_c and \delta {ϵ}_{\gamma } using Eq. (C.2) introduces source terms proportional to A_0. Therefore, we find 4 homogeneous solutions and one inhomogeneous solution as follows.

1) Homogeneous solution I: r takes the value

(C.8)\begin{array}{c}{r}_1=0,\end{array}

and the coefficients are

(C.9)\begin{array}{c}\begin{array}{rl}& c_1=-\frac{c_4(\alpha -2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{2D_1},\\ & c_2=\frac{c_4(\alpha +3)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{D_1},\\ & c_3=-\frac{3}{4}{c}_4\alpha H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}.\end{array}\end{array}

2) Homogeneous solution II: r takes the value

(C.10)\begin{array}{c}{r}_2=\frac{\alpha }{2}+\sqrt{{\displaystyle \frac{{\alpha }^2}{4}}+{\displaystyle \frac{\left(2-\alpha \right){\xi }_1}{\kappa }}},\end{array}

and the coefficients are

(C.11)\begin{array}{c}\begin{array}{rl}{c}_1=& \frac{c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{2\alpha D_1\kappa }\\ & \times \left[\sqrt{\kappa \left({\alpha }^2\kappa -16(\alpha -2){\xi }_1\right)}+4\kappa \right],\\ c_2=& -\frac{2c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{\alpha D_1},\\ c_3=& \frac{c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{2\alpha \kappa }\\ & \times \left[\sqrt{\kappa \left({\alpha }^2\kappa -16(\alpha -2){\xi }_1\right)}-\alpha \kappa \right].\end{array}\end{array}

3) Homogeneous solution III: r takes the value

(C.12)\begin{array}{c}{r}_3={\displaystyle \frac{\alpha }{2}}-\sqrt{{\displaystyle \frac{{\alpha }^2}{4}}+{\displaystyle \frac{\left(2-\alpha \right){\xi }_1}{\kappa }}},\end{array}

and the coefficients are

(C.13)\begin{array}{c}\begin{array}{rl}{c}_1=& -\frac{c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{2\alpha D_1\kappa }\\ & \times \left[\sqrt{\kappa \left({\alpha }^2\kappa -16(\alpha -2){\xi }_1\right)}-4\kappa \right],\\ c_2=& -\frac{2c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{\alpha D_1},\\ c_3=& -\frac{c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{2\alpha \kappa }\\ & \times \left[\sqrt{\kappa \left({\alpha }^2\kappa -16(\alpha -2){\xi }_1\right)}+\alpha \kappa \right].\end{array}\end{array}

4) Homogeneous solution IV: r takes the value

(C.14)\begin{array}{c}{r}_4=\alpha ,\end{array}

and the coefficients are

(C.15)\begin{array}{c}\begin{array}{rl}& c_1=\frac{c_4(\alpha +2)H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{2D_1},\\ & c_2=-\frac{c_4{H}_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}{D_1},\\ & c_3=\frac{1}{4}{c}_4\alpha H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}.\end{array}\end{array}

5) Inhomogeneous solution: r takes the value

(C.16)\begin{array}{c}{r}_0=-{\displaystyle \frac{\alpha }{2}},\end{array}

and the coefficients are

\begin{array}{rl}{c}_1=& -\frac{16(\alpha -2)(\alpha -1)(\alpha +2){\mathrm{\Omega }}_{r0}{A}_0}{3\alpha (\alpha +4)D_1^2\left(15{\alpha }^2\kappa +16(\alpha -2){\xi }_1\right)}\\ & \times \frac{3\alpha (\alpha +12)\kappa -16(\alpha -2){\xi }_1}{2(\alpha +2){\tilde{V}}_1+(\alpha -2){\xi }_1},\end{array}

(C.17)\begin{array}{c}\begin{array}{rl}{c}_2=& -\frac{16(\alpha -2)(\alpha +2){\mathrm{\Omega }}_{r0}{A}_0}{3\alpha (\alpha +4)D_1^2\left(15{\alpha }^2\kappa +16(\alpha -2){\xi }_1\right)}\\ & \times \frac{3\alpha (\alpha +12)\kappa -16(\alpha -2){\xi }_1}{2(\alpha +2){\tilde{V}}_1+(\alpha -2){\xi }_1},\\ c_3=& \frac{8(\alpha -2)(\alpha +2){\mathrm{\Omega }}_{r0}{A}_0}{3\alpha (\alpha +4)D_1\left(15{\alpha }^2\kappa +16(\alpha -2){\xi }_1\right)}\\ & \times \frac{15{\alpha }^3\kappa +32(\alpha -2)(2\alpha +3){\xi }_1}{2(\alpha +2){\tilde{V}}_1+(\alpha -2){\xi }_1},\\ c_4=& -\frac{32(\alpha +2){\mathrm{\Omega }}_{r0}{A}_0}{3\alpha (\alpha +4)D_1\left(15{\alpha }^2\kappa +16(\alpha -2){\xi }_1\right)}\\ & \times \frac{3{\alpha }^2\kappa +8(\alpha -2){\xi }_1}{H_0{\xi }_1{\left(1-{\mathrm{\Omega }}_{V_1}\right)}^{3/2}}.\end{array}\end{array}

With the approximate solutions of \mathrm{\Psi } and \mathrm{\Phi }, approximate solutions of the velocity variables and the higher-order multipoles can be found from Eq. (B.14). First, the equation for v_c can be directly solved to get

(C.18)\begin{array}{c}{v}_c\approx -{\displaystyle \frac{k{c}_1}{H_0\left(1-{\displaystyle \frac{\alpha }{2}}+r\right)\sqrt{1-{\mathrm{\Omega }}_{V_1}}}}{a}^{-\frac{\alpha }{2}+r}.\end{array}

Then, for v_b and {\mathrm{\Theta }}_1, using the tight-coupling approximation

(C.19)\begin{array}{c}{v}_b\approx -3{\mathrm{\Theta }}_1,\end{array}

the equation for {\mathrm{\Theta }}_1 becomes

(C.20)\begin{array}{c}{\mathrm{\Theta }}_1^{\prime }\approx {\displaystyle \frac{k}{3}}\left(\mathrm{\Phi }+\mathrm{\Psi }\right)+{\displaystyle \frac{k{A}_0}{3}}.\end{array}

For the homogeneous solutions of \mathrm{\Phi } and \mathrm{\Psi }, we find

(C.21)\begin{array}{c}{\mathrm{\Theta }}_1\approx {\displaystyle \frac{k(c_1+c_2)}{3H_0\left(r-{\displaystyle \frac{\alpha }{2}}\right)\sqrt{1-{\mathrm{\Omega }}_{V_1}}}}{a}^{-\frac{\alpha }{2}+r},\end{array}

and for the inhomogeneous solution of \mathrm{\Phi } and \mathrm{\Psi }, we find

(C.22)\begin{array}{c}{\mathrm{\Theta }}_1\approx -{\displaystyle \frac{2k{A}_0}{3\alpha H_0\sqrt{1-{\mathrm{\Omega }}_{V_1}}}}{a}^{-\frac{\alpha }{2}}.\end{array}

Finally, the equations for the higher-order multipoles give the recurrence relations

(C.23)\begin{array}{c}\begin{array}{rl}& {\mathrm{\Theta }}_2\approx \frac{4k}{9}\frac{a^2}{{\mathrm{\Gamma }}_0}{\mathrm{\Theta }}_1,\\ & {\mathrm{\Theta }}_l\approx \frac{lk}{2l+1}\frac{a^2}{{\mathrm{\Gamma }}_0}{\mathrm{\Theta }}_{l-1},l\ge 3,\end{array}\end{array}

where we have used \mathrm{\Gamma }\approx {\mathrm{\Gamma }}_0/a^2 at a\to 0, with {\mathrm{\Gamma }}_0 being a constant.

Having the approximate solutions for the velocity variables and the higher-order multipoles, we need to check if neglecting them at the beginning when solving {\delta }_b,{\delta }_c,{\mathrm{\Theta }}_0 and \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S is justified or not. It is straightforward to verify that the velocity variables and the higher-order multipoles can be neglected in solving the approximate solutions of {\delta }_b,{\delta }_c,{\mathrm{\Theta }}_0 and the homogeneous solutions of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S, while the photon velocity variable \delta u_{\gamma S}=-3a{\mathrm{\Theta }}_1/k does contribute relevant terms in solving the inhomogeneous solutions of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S. Knowing that the inhomogeneous solutions of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S take the form of Eq. (C.7) with r=r_0=-\alpha /2, the coefficients c_1,c_2,c_3,c_4 can be recalculated taking into consideration of the contribution of \delta u_{\gamma S}. In fact, we have shown in Eq. (C.17) the correct coefficients for the inhomogeneous solutions of \mathrm{\Psi },\mathrm{\Phi },\delta b_{\eta },\delta b_S, as it is unnecessary to show the inaccurate result obtained with \delta u_{\gamma S} dropped.

Equipped with the above approximate solutions at a\to 0, numerical integrations for solving the ODE system consists of Eq. (B.14) and Eq. (B.13) in the Fourier space can be carried out from a small enough value of a (\ln a=-14 in our numerical code), with the initial values of the variables given by the sum of the 5 approximate solutions. In fact, we need to point out that out of the 5 approximate solutions, only the inhomogeneous solution is relevant for the CMB calculation, because it is the mode corresponding to an initial perturbation in the CMB photons’ energy density. The 4 homogeneous modes, including one mode corresponding to an initial metric perturbation and three modes corresponding to initial perturbations in the bumblebee field, turn out to contribute insignificantly in calculating the CMB anisotropy. This is similar to GR where there is no mode of bumblebee perturbation and therefore only one homogeneous mode corresponding to an initial metric perturbation.

The perturbation equations in Eq. (B.13) apply to GR by setting b_{\eta }=\delta b_{\eta }=\delta b_S=0 and adding {ϵ}_{\mathrm{\Lambda }}=\mathrm{\Lambda }/\kappa to the combination {ϵ}_b+{ϵ}_c+{ϵ}_{\gamma } in the third equation. Using the above discussed method, one homogeneous solution can be found by closing the equation of \mathrm{\Phi } using {\mathrm{\Theta }}_0\approx \mathrm{\Phi },

(C.24)\begin{array}{c}\begin{array}{rl}& \mathrm{\Psi }\approx \mathrm{\Phi }\approx {\mathrm{\Theta }}_0\approx \frac{c_1}{a^3},\\ & {\mathrm{\Theta }}_1\approx -\frac{k{c}_1}{3H_0\sqrt{{\mathrm{\Omega }}_{r0}}}\frac{1}{a^2},\\ & v_b\approx v_c\approx \frac{k{c}_1}{H_0\sqrt{{\mathrm{\Omega }}_{r0}}}\frac{1}{a^2},\end{array}\end{array}

where c_1 is an arbitrary constant, and one inhomogeneous solution can be found by closing the equation of \mathrm{\Phi } using {\mathrm{\Theta }}_0\approx \mathrm{\Phi }+A_0,

(C.25)\begin{array}{c}\begin{array}{rl}& \mathrm{\Psi }\approx \mathrm{\Phi }\approx -\frac{2A_0}{3},\\ & {\mathrm{\Theta }}_0\approx \frac{A_0}{3},\\ & {\mathrm{\Theta }}_1\approx -\frac{k{A}_0}{9H_0\sqrt{{\mathrm{\Omega }}_{r0}}}a,\\ & v_b\approx v_c\approx \frac{k{A}_0}{3H_0\sqrt{{\mathrm{\Omega }}_{r0}}}a.\end{array}\end{array}

For both solutions, the higher-order multipoles are given by Eq. (C.23). A set of general initial data is a combination of the homogeneous solution in Eq. (C.24) and the inhomogeneous solution in Eq. (C.25). But as the homogeneous solution contributes insignificantly, it is sufficient to use only the inhomogeneous solution to set up the initial condition for numerically integrating the full set of perturbation equations.

8 Appendix D: Checking our CMB code by applying it to the ΛCDM model

Using the conformal expansion rate for the \mathrm{\Lambda }CDM model in GR as shown in Eq. (12), switching off b_{\eta },\delta b_{\eta },\delta b_S, and adding {ϵ}_{\mathrm{\Lambda }}=\mathrm{\Lambda }/\kappa in Eq. (B.13), our CMB code applies to the \mathrm{\Lambda }CDM model with {\mathrm{\Omega }}_{\mathcal{K}0}=0. By calculating the CMB power spectrum in the \mathrm{\Lambda }CDM model and comparing with results from CAMB), we can have an estimate for the error in the results produced using our code.

As explained in Section 3, there are three steps in calculating the CMB temperature power spectrum and the results can be shown at three levels: (i) the solutions for the perturbation variables, (ii) the transfer function {\mathrm{\Theta }}_l(k), and (iii) the final power spectrum C_l. Figures A.1-A.3 show comparisons between results from our code and from CAMB. The parameters of the \mathrm{\Lambda }CDM model calculated are

(D.1)\begin{array}{c}\begin{array}{rl}& H_0=70\mathrm{k}\mathrm{m}\cdot {\mathrm{s}}^{-1}\cdot {\mathrm{M}\mathrm{p}\mathrm{c}}^{-1},\\ & {\mathrm{\Omega }}_{b0}=0.05,\\ & {\mathrm{\Omega }}_{c0}=0.25,\\ & T_{\mathrm{C}\mathrm{M}\mathrm{B}}=2.7\mathrm{K},\\ & n_s=1.\end{array}\end{array}

Neutrinos and reionization are switched off when using CAMB to match the simple cosmological model in our code.

References

[1] C. L. Bennett , A. Banday , K. M. Gorski , G. Hinshaw , P. Jackson , P. Keegstra , A. Kogut , G. F. Smoot , D. T. Wilkinson , and E. L. Wright , Four year COBE DMR cosmic microwave background observations: Maps and basic results, Astrophys. J. 464, L1 (1996)

[2] D. J. Fixsen , E. S. Cheng , J. M. Gales , J. C. Mather , R. A. Shafer , and E. L. Wright , The cosmic microwave background spectrum from the full COBE FIRAS data set, Astrophys. J. 473, 576 (1996)

[3] N. Aghanim, ., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020) [Erratum: Astron. Astrophys. 652, C4 (2021)], arXiv:

[4] T. Clifton , P. G. Ferreira , A. Padilla , and C. Skordis , Modified Gravity and Cosmology, Phys. Rep. 513, 1 (2012)

[5] K. Bamba , S. Capozziello , S. Nojiri , and S. D. Odintsov , Dark energy cosmology: the equivalent description via different theoretical models and cosmography tests, Astrophys. Space Sci. 342, 155 (2012)

[6] K. Koyama , Cosmological tests of modified gravity, Rep. Prog. Phys. 79(4), 046902 (2016)

[7] P. A. R. Ade , . Planck 2015 results. XIV. Dark energy and modified gravity, Astron. Astrophys. 594, A14 (2016)

[8] E. N. Saridakis,R. Lazkoz,V. Salzano,P. Vargas Moniz,S. Capozziello,J. Beltrań Jimeńez,M. De Laurentis,G. J. Olmo (Eds.), Modified Gravity and Cosmology: An Update by the CANTATA Network, Springer, 2021, arXiv:

[9] V. A. Kostelecky , Gravity , Lorentz violation, and the standard model, Phys. Rev. D 69, 105009 (2004)

[10] R. W. Hellings and K. Nordtvedt , Vector-metric theory of gravity, Phys. Rev. D 7, 3593 (1973)

[11] R. Bluhm , V. A. Kostelecky , and Spontaneous Lorentz violation , Nambu−Goldstone modes, and gravity, Phys. Rev. D 71, 065008 (2005)

[12] O. Bertolami and J. Paramos , The Flight of the bumblebee: Vacuum solutions of a gravity model with vector-induced spontaneous Lorentz symmetry breaking, Phys. Rev. D 72, 044001 (2005)

[13] Q. G. Bailey and V. A. Kostelecky , Signals for Lorentz violation in post-Newtonian gravity, Phys. Rev. D 74, 045001 (2006)

[14] R. Bluhm , S. H. Fung , and V. A. Kostelecky , Spontaneous Lorentz and diffeomorphism violation, massive modes, and gravity, Phys. Rev. D 77, 065020 (2008)

[15] R. Bluhm,N. L. Gagne,R. Potting,A. Vrublevskis, Constraints and stability in vector theories with spontaneous Lorentz violation, Phys. Rev. D 77, 125007 (2008) [Erratum: Phys. Rev. D 79, 029902(E)], arXiv:

[16] V. A. Kostelecky and J. D. Tasson , Matter-gravity couplings and Lorentz violation, Phys. Rev. D 83, 016013 (2011)

[17] A. Ovgün , K. Jusufi , and I. Sakallı , Exact traversable wormhole solution in bumblebee gravity, Phys. Rev. D 99(2), 024042 (2019)

[18] R. Casana , A. Cavalcante , F. P. Poulis , and E. B. Santos , Exact Schwarzschild-like solution in a bumblebee gravity model, Phys. Rev. D 97(10), 104001 (2018)

[19] R. Xu , D. Liang , and L. Shao , Static spherical vacuum solutions in the bumblebee gravity model, Phys. Rev. D 107(2), 024011 (2023)

[20] Z. F. Mai , R. Xu , D. Liang , and L. Shao , Extended thermodynamics of the bumblebee black holes, Phys. Rev. D 108(2), 024004 (2023)

[21] Z. F. Mai , R. Xu , D. Liang , and L. Shao , Dynamic instability analysis for bumblebee black holes: The odd parity, Phys. Rev. D 109(8), 084076 (2024)

[22] R. Xu , D. Liang , and L. Shao , Bumblebee black holes in light of event horizon telescope observations, Astrophys. J. 945(2), 148 (2023)

[23] X. Zhu,R. Xu,D. Xu, Bumblebee cosmology: Tests using distance and time-redshift probes, arXiv: 2024) (accepted for publication in Physics of the Dark Universe)

[24] A. Lewis , A. Challinor , and A. Lasenby , Efficient computation of CMB anisotropies in closed FRW models, Astrophys. J. 538, 473 (2000)

[25] D. Blas , J. Lesgourgues , and T. Tram , The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, J. Cosmol. Astropart. Phys. 07, 034 (2011)

[26] G. B. Zhao , L. Pogosian , A. Silvestri , and J. Zylberberg , Cosmological tests of general relativity with future tomographic surveys, Phys. Rev. Lett. 103, 241301 (2009)

[27] A. Hojjati , L. Pogosian , and G. B. Zhao , Testing gravity with CAMB and CosmoMC, J. Cosmol. Astropart. Phys. 08, 005 (2011)

[28] A. G. Riess , . Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998)

[29] R. Cayrel , . Measurement of stellar age from uranium decay, Nature 409, 691 (2001)

[30] N. Aghanim , . Planck 2018 results. I. Overview and the cosmological legacy of Planck, Astron. Astrophys. 641, A1 (2020)

[31] C. P. Ma and E. Bertschinger , Cosmological perturbation theory in the synchronous and conformal Newtonian gauges, Astrophys. J. 455, 7 (1995)

[32] P. Callin, How to calculate the CMB spectrum, arXiv: 2006)

[33] R. Durrer, The Cosmic Microwave Background, Cambridge University Press, Cambridge, 2008

[34] D. Baumann, Cosmology, Cambridge University Press, 2022

[35] P. J. E. Peebles , Recombination of the primeval plasma, Astrophys. J. 153, 1 (1968)

[36] A. H. Guth , The inflationary universe: A possible solution to the horizon and flatness problems, Phys. Rev. D 23, 347 (1981)

[37] A. D. Linde, A new inflationary universe scenario: A possible solution of the horizon, flatness, homogeneity, isotropy and primordial monopole problems, Phys. Lett. B 108, 389 (1982)

[38] W. L. Freedman , . Final results from the Hubble Space Telescope key project to measure the Hubble constant, Astrophys. J. 553, 47 (2001)

[39] D. J. Fixsen , The temperature of the cosmic microwave background, Astrophys. J. 707, 916 (2009)

[40] S. Cole , . The 2dF galaxy redshift survey: Powerspectrum analysis of the final dataset and cosmological implications, Mon. Not. Roy. Astron. Soc. 362, 505 (2005)

[41] R. V. Maluf and J. C. S. Neves , Bumblebee field as a source of cosmological anisotropies, J. Cosmol. Astropart. Phys. 10, 038 (2021)

[42] A. De Felice,A. Hell, On the cosmological degrees of freedom of Proca field with non-minimal coupling to gravity, arXiv: 2025)

RIGHTS & PERMISSIONS

Higher Education Press