VARYING ALPHA FROM N-BODY SIMULATIONS

The Lagrangian density for the BSBM model could be written as

$$\begin{equation} \mathcal {L} = \frac{1}{2}\left[ \frac{R}{\kappa }-\nabla ^{a}\varphi \nabla _{a}\varphi \right]-\mathcal {L}_{m}- e^{-2\sqrt{\kappa }\varphi } \mathcal {L}_{\mathrm{EM}}-\mathcal {L}_{r}, \end{equation} \tag{ 1 }$$

where R is the Ricci scalar, κ = 8πG with G being the gravitational constant, φ is the scalar field; $\mathcal {L}_{m}, \mathcal {L}_{\mathrm{EM}}, \mbox{and } \mathcal {L}_{r}$ represent respectively the Lagrangian densities for dust, electromagnetic field (including photons), and other radiation (such as neutrinos). The coupling function between the scalar field and the electromagnetic field in the BSBM model is $e^{-2\sqrt{\kappa }\varphi }$, where $\sqrt{\kappa }$ is added so that $\sqrt{ \kappa }\varphi \equiv \psi$ is dimensionless. In the simplest version of the model, there is no potential for the scalar field.

The dust Lagrangian for a point particle with mass m₀ is

$$\begin{equation} \mathcal {L}_{m}(\mathbf {y}) = {-}\frac{m_{0}}{\sqrt{-g}}\delta (\mathbf {y}- \mathbf {x}_{0})\sqrt{g_{ab}\dot{x}_{0}^{a}\dot{x}_{0}^{b}}, \end{equation} \tag{ 2 }$$

where y is the general coordinate and x₀ is the coordinate of the center of the particle. From this equation, we derive the corresponding energy–momentum tensor

$$\begin{equation} T_{m}^{ab} = \frac{m_{0}}{\sqrt{-g}}\delta (\mathbf {y}-\mathbf {x}_{0})\dot{ x}_{0}^{a}\dot{x}_{0}^{b}. \end{equation} \tag{ 3 }$$

Also, because $g_{ab}\dot{x}_{0}^{a}\dot{x}_{0}^{b}\equiv g_{ab}u^{a}u^{b}=1,$ in which u^a is the 4-velocity of the dark matter particle centered at x₀, the Lagrangian can be rewritten as

$$\begin{equation} \mathcal {L}_{m}(\mathbf {y}) = {-}\frac{m_{0}}{\sqrt{-g}}\delta (\mathbf {y}- \mathbf {x}_{0}). \end{equation} \tag{ 4 }$$

This result will be used below.

Equation (3) is just the energy–momentum tensor for a single matter particle. For a fluid consisting of many particles, the energy–momentum tensor will be

$$\begin{eqnarray} T_{m}^{ab} &=& \frac{1}{\mathcal {V}}\int _{\mathcal {V}}d^{4}y\sqrt{-g}\frac{ m_{0}}{\sqrt{-g}}\delta (y-x_{0})\dot{x}_{0}^{a}\dot{x}_{0}^{b}\nonumber\\ &=&\rho _{\mathrm{CDM}}u^{a}u^{b}, \end{eqnarray} \tag{ 5 }$$

where $\mathcal {V}$ denotes a volume that is microscopically large but macroscopically small, and we have extended the three-dimensional δ function to a four-dimensional one by adding a time component. Here, u^a is the averaged 4-velocity of the matter fluid.

Using

$$\begin{equation} T^{ab} = {-}\frac{2}{\sqrt{-g}}\frac{\delta (\sqrt{-g}\mathcal {L}) }{\delta g_{ab}}, \end{equation} \tag{ 6 }$$

it is straightforward to show that the energy–momentum tensor for the scalar field is

$$\begin{equation} T^{\varphi ab} = \nabla ^{a}\varphi \nabla ^{b}\varphi -\case{1}{2} g^{ab}\nabla _{c}\varphi \nabla ^{c}\varphi. \end{equation} \tag{ 7 }$$

Therefore, the total energy–momentum tensor is

$$\begin{eqnarray} T_{ab} &=& \nabla _{a}\varphi \nabla _{b}\varphi - \frac{1}{2} g_{ab}\nabla _{c}\varphi \nabla ^{c}\varphi\nonumber\\ && +T^{m}_{ab}+T^{r}_{ab} +e^{-2\sqrt{\kappa }\varphi }T_{ab}^{\mathrm{EM}}, \end{eqnarray} \tag{ 8 }$$

where T^m_ab = ρ_mu_au_b, T^r_ab is the energy–momentum tensor for radiation fields except photons, and T^EM_ab for photons. The Einstein equations are

$$\begin{equation} G_{ab} = \kappa T_{ab} \end{equation} \tag{ 9 }$$

in which $G_{ab}=R_{ab}-\frac{1}{2}g_{ab}R$ is the Einstein tensor. Note that due to the extra coupling between the scalar field, φ, and the electromagnetic field, the energy–momentum tensors for either will no longer be separately conserved, but instead we have

$$\begin{equation} \nabla _{b}T_{\mathrm{EM}}^{ab}= 2\sqrt{\kappa }\big(g^{ab}\mathcal {L}_{ \mathrm{EM}}+T_{\mathrm{EM}}^{ab}\big) \nabla _{b}\varphi. \end{equation} \tag{ 10 }$$

However, the total energy–momentum tensor is certainly conserved.

Meanwhile, the scalar-field equation of motion is

$$\begin{equation} \square \varphi = {-}2\sqrt{\kappa }e^{-2\sqrt{\kappa }\varphi }\mathcal {L}_{ \mathrm{EM}}, \end{equation} \tag{ 11 }$$

where □ ≡ ∇^a∇_a. This equation governs the time evolution and spatial configuration of the scalar field.

Equations (8)–(11) summarize all the physics needed for the following analysis. However, when making use of them we should also specify the form of the electromagnetic matter. For example, if it is a pure electromagnetic field (photons), then we have $\mathcal {L}_{\mathrm{EM}}=\frac{1}{2} (E^2-B^2)=0$ in which EandB stand for the electric and magnetic fields, respectively. Thus, from the time component of Equation (10) we obtain the (background) evolution equation for photon density as

$$\begin{equation} \dot{\rho }_r + 4H\rho _r = 2\dot{\psi }\rho _{r}, \end{equation} \tag{ 12 }$$

where H is the Hubble expansion rate and remember that $\psi =\sqrt{\kappa }\varphi$.

It then seems that, by the same reason, the right-hand side of Equation (11) also vanishes, leaving the scalar field unsourced. This may not be true however, as non-relativistic matter could also contribute to $\mathcal {L}_{\mathrm{EM}}$ and thus T^EM_ab. For example, in baryonic matter $\mathcal {L}_{\mathrm{EM}}\approx E^{2}/2$, and for neutrons and protons this electromagnetic contribution to the total mass can be of order 10⁻⁴; in superconducting cosmic strings $\mathcal {L}_{\mathrm{EM} }\approx -B^{2}/2,$ where ρ_EM ≈ B²/2 so that $| \mathcal {L}_{\mathrm{EM}}/\rho _{\mathrm{EM}}|\sim 1$. In the BSBM model, in order to simplify the situation, it is assumed that $\mathcal {L}_{ \mathrm{EM}}/\rho _{m}=\zeta$, where ζ is a constant with a modulus between 0 and ≈1, either positive or negative, and ρ_m is the density for non-relativistic matter.

Thus, the scalar-field equation gets sourced by a term proportional to ζ:

$$\begin{equation} \square \varphi = {-}2\sqrt{\kappa }\zeta e^{-2\sqrt{\kappa }\varphi }\rho _{m}. \end{equation} \tag{ 13 }$$

Since the part of $\mathcal {L}_{\mathrm{EM}}$ that affects the scalar field is a constituent of the non-relativistic matter and is presumably moving with the matter particles, we could combine Equation (10) and the conservation equation for the dust matter (no including electromagnetic contribution) to write a new conservation equation for the particle:

$$\begin{equation} \nabla _{b}T_{m+\mathrm{EM}}^{ab} = 2\sqrt{\kappa }\big(g^{ab}\mathcal {L}_{ \mathrm{EM}}+T_{\mathrm{EM}}^{ab}\big) \nabla _{b}\varphi. \end{equation} \tag{ 14 }$$

Although we have assumed above that $\mathcal {L}_{\mathrm{EM}}=\zeta \rho _m$, we still lack knowledge about T^ab_EM, whose relation to $\mathcal {L}_{\mathrm{EM}}$ could be complicated. Here, for simplicity, we assume that T^ab_EM = −ζρ_mu^au^b. Then it is easy to find that the time component of this equation reads

$$\begin{equation} \dot{\rho }_{m+\mathrm{EM}}+3H\rho _{m+\mathrm{EM}} = 0, \end{equation} \tag{ 15 }$$

while the ith spatial component of it gives the following (modified) geodesic equation:

$$\begin{equation} \ddot{x}^{i}_{0} + \Gamma ^{i}_{ab}\dot{x}^{a}_{0}\dot{x}^{b}_{0} = 2\zeta \sqrt{\kappa } (g^{ib}-u^{i}u^{b})\nabla _{b}\varphi, \end{equation} \tag{ 16 }$$

where x₀ is the coordinate of the center of a particle, and the right-hand side represents a fifth force on the particle (Li & Zhao 2009, 2010). The assumption T^ab_EM = −ζρ_mu^au^b might seem unappealing, but as we shall see below, because |ζ| ≪ 1, the fifth force is much weaker than gravity and thus has negligible effects in the clustering of matter in any case; ultimately, it is only the BSBM assumption $\mathcal {L}_{\mathrm{EM}}=\zeta \rho _m$ that is important in theoretical predictions of the spatial and temporal variations of α, given that $\alpha =e^{2\psi }\frac{e^{2}_0}{c\hbar }$, where e₀ is the present-day value for e.

The scalar-field equation of motion, which controls the dynamics of the scalar field φ, is generally complicated, because it depends nonlinearly on φ, which both evolves in time and fluctuates in space. Fortunately, for the majority of applications the scalar-field potential and coupling function are not nonlinear enough to give the scalar field a very heavy mass so as to make it fluctuate strongly. The nice thing about this is that in certain places of the equations one may then forget the scalar-field perturbation and simplify these equations accordingly. In Li & Barrow (2011), for example, it is shown that such a simplification is a very good approximation (see, however, Li & Zhao 2009, 2010 for an opposite extreme for which the scalar-field potential is very nonlinear so that such simplification does not work).

In the BSBM model, there is no scalar-field potential and the coupling function is close to linear if $\sqrt{\kappa }|\varphi |\ll 1$ (which is the case for parameter space that we are interested in; Barrow et al. 2002a). We therefore expect the fluctuation of φ to be very weak and assume that $\sqrt{\kappa }|\delta \varphi |\ll \sqrt{\kappa }|\varphi |\ll 1$ (which we shall confirm below using numerical simulations). Under such an assumption the Poisson equation could be written as

$$\begin{eqnarray} \nabla _{\mathbf {x}}^{2}\Phi &=&\! 4\pi Ga^{3}\rho _m[1+\zeta e^{-2\sqrt{ \kappa }(\bar{\varphi }+\delta \varphi)}]\,{-}\, 4\pi Ga^{3}\bar{\rho }_m[1+\zeta e^{-2\sqrt{\kappa }\bar{\varphi }} ]\nonumber\\ &\approx & 4\pi Ga^3(\rho _{m}-\bar{\rho }_m), \end{eqnarray} \tag{ 17 }$$

where $\bar{\varphi }$ is the background value of φ and δφis its perturbation; $\rho _m, \bar{\rho }_m$ are respectively the local and background matter density; Φ is the gravitational potential and a is the cosmic scale factor; ∇_x is the derivative with respect to the comoving coordinate x. To obtain Equation (17), we have used the fact that in the BSBM model $|\zeta e^{-2\bar{\varphi } }|\ll 1$. Note that Equation (17) clearly indicates that the gravitational potential is essentially not influenced by the scalar field φ.

For the scalar-field equation of motion Equation (13), because the background part (which has no spatial dependence) can be solved easily, we subtract that from the full equation to obtain an equation of motion for δφ only (remember that $\varphi =\bar{\varphi }+\delta \varphi$). Furthermore, we drop the time derivative terms of δφ, as they are small compared with the spatial gradients (i.e., work in the quasi-static limit). The final equation for δφ then becomes

$$\begin{eqnarray} \nabla ^{2}_{\mathbf {x}}(a\sqrt{\kappa }\delta \varphi) &=& 2\zeta \kappa [\rho _me^{-2\sqrt{\kappa }(\bar{\varphi }+ \delta \varphi)}-\bar{\rho }_me^{-2\sqrt{\kappa }\bar{\varphi }}]a^3\nonumber\\ &\approx & 16\pi Ga^3\zeta (\rho _m-\bar{\rho }_m), \end{eqnarray} \tag{ 18 }$$

where we have used κ = 8πG and $\sqrt{\kappa }|\delta \varphi |\ll \sqrt{ \kappa }|\varphi |\ll 1$.

Comparing Equations (17) and (18), it is evident that the source terms are the same up to a constant coefficient 4ζ. Consequently, we shall have

$$\begin{equation} a\sqrt{\kappa }\delta \varphi (\mathbf {x}) \approx 4\zeta \Phi (\mathbf {x}). \end{equation} \tag{ 19 }$$

Note that this equation, together with the geodesic equation (16), implies that the magnitude of the fifth force (force due to exchange of scalar-field quanta between particles) |f| satisfies

$$\begin{equation} |\mathbf {f}| \sim \zeta |\vec{\nabla }(a\sqrt{\kappa }\delta \varphi)| \sim 4\zeta ^2|\vec{\nabla }\Phi |. \end{equation} \tag{ 20 }$$

Therefore, the ratio between the magnitudes of the fifth force and gravity is of order ζ² ≲ 10⁻¹² to 10⁻⁸ ≪ 1. This implies that the fifth force cannot significantly influence the structure formation.

To study in details the behavior of the scalar field and hence the fine structure constant α in the BSBM model, we have performed N-body simulations for four different models, with ζ = ±2 × 10⁻⁶ and ±5 × 10⁻⁶, respectively. The physical parameters that we adopt in all simulations are as follows: the present-day dark-energy fractional-energy density Ω_DE = 0.743 and Ω_m = Ω_CDM + Ω_B = 0.257, H₀ = 71.9 km s⁻¹ Mpc⁻¹, n_s = 0.963, and σ₈ = 0.761. These are in accordance with the concordance cosmological model preferred by current data sets. Our simulation box has a size of 64 h⁻¹ Mpc, in which h = H₀/(100 km s^-1 Mpc⁻¹). In all of those simulations, the mass resolution is 1.114 × 10⁹ h⁻¹ M_☉; the particle number is 256³; the domain grid (i.e., the coarsest grid which covers the whole simulation box) has 128³ equal-sized cubic cells; and the finest refined grids have 16,384 cells on each side, corresponding to a force resolution of ∼12 h⁻¹ kpc. Detailed description about the N-body simulation technique and code for the (coupled) scalar-field models could be found in Li & Zhao (2009, 2010) and will not be presented here.

Because the BSBM model (with the parameter ζ constrained by data) involves a very weak coupling between matter and the scalar field φ, the presence of the latter and its coupling have negligible influences on the background (ΛCDM) cosmology, although the opposite is not true since there is no scalar-field potential and thus the dynamics of φ is controlled entirely by the coupling. In our simulations, we compute the full background cosmology and evolution of φ on a predefined time grid using the MAPLE mathematical software, and interpolate to obtain the corresponding quantities which are needed in N-body simulations. Details of this procedure can be found in Appendix C of Li & Barrow (2011), and Figure 1 shows the background evolution of $\sqrt{\kappa } \varphi$ for the four models considered here. Obviously the condition $\sqrt{ \kappa }\bar{\varphi }\ll 1$ is satisfied, justifying the approximation we used above to derive Equation (19).

Zoom In Zoom Out Reset image size

One nice thing about the BSBM model is its simplicity, and it turns out that the background evolution of φ (and thus α) in different cosmic epochs can be well described by some analytical formulae (Barrow et al. 2002a). Therefore, in the present study, we shall mainly focus on the spatial variation of φ and α (especially in virialized halos).

As mentioned above, because there is no potential for the scalar field φ and because $\sqrt{\kappa }\varphi \ll 1$ for our choices of ζ, the scalar-field equation of motion (in the quasi-static limit) for δφ and the Poisson equation share the same source up to a constant coefficient, and therefore we expect aδφ ∝ Φ across the whole space. This is what we have found in Li & Barrow (2011) for a different coupled scalar-field model where the scalar-field potential is negligible. Indeed, this could serve as a test of the scalar-field solver in our N-body simulation code.

To check that our code does recover the analytical approximation, we have plotted in Figure 2 the comparison of the scalar-field perturbation $a\sqrt{\kappa }\delta \varphi$ and gravitational potential Φ from a slice of the simulation box. As indicated by this figure, the agreement between the numerical results (green dots) and analytical approximation (black solid line) is remarkably good, implying that the numerical code works well. Therefore to a high precision we can assume that $a\sqrt{\kappa }\delta \varphi \propto \Phi$ everywhere, a fact that will be used below to obtain an analytical expression of δφ in dark matter halos.

Zoom In Zoom Out Reset image size

In the standard picture, the galaxies where observers live generally locate inside the dark matter halos, which to the simplest approximation are just spherical clusters of matter with a universal NFW (Navarro et al. 1996) density profile.

We are certainly interested in the (possible) variation of α inside the halo in which we reside. For example, there has been a great deal of analytical work on how significantly the local value of α could deviate from its cosmological counterpart (Mota & Barrow 2004a, 2004b; Jacobson 1999; Wetterich 2003; Shaw & Barrow 2006a, 2006b, 2006c). Also, if the spatial variation of α is strong enough, then it might have an impact on our observation of the spectra for the stars from our Galaxy and other galaxies.

From our simulation output, it is easy to identify the dark matter halos (Knebe & Gibson 2004; Li & Barrow 2011). Now what we want to do here is to measure the quantity $\sqrt{\kappa }\delta \varphi$ as a function of distance R to the halo centers, assuming that the halos are exactly spherical. For this we have recorded the value for $\sqrt{\kappa }\delta \varphi$ at the position of each particle, and then presumably we could divide each halo into a number of spherical shells, determine the radius of each shell and compute the average value of $\sqrt{\kappa }\delta \varphi$ in the shells. However, there is some subtlety in the computation of the averaged $\sqrt{\kappa } \delta \varphi$.

The problem is that, since we have only recorded the information for $\sqrt{ \kappa }\delta \varphi$ at the positions of the particles, we do not have fair sampling points. Because $\sqrt{\kappa }\delta \varphi$ tends to be different in regions of different particle number densities, the high-density regions will be oversampled and low-density regions undersampled, resulting in a bias as we are trying to determine the spatially averaged, rather than the particle-averaged, value of $\sqrt{\kappa } \delta \varphi$. To test how large the bias could be, we use an approximation as follows: first, we divide each spherical shell into N equal-sized volumes which are small enough so that the particle number density does not change much inside each of them, then we compute the particle average of $\sqrt{\kappa }\delta \varphi$ in each of these volumes, call it $\langle \sqrt{\kappa }\delta \varphi \rangle _{i}$ for i = 1, ..., N, then the spatial (volume) average of $\sqrt{\kappa }\delta \varphi$ is given by

$$\begin{equation} \langle \sqrt{\kappa }\delta \varphi \rangle _{\mathrm{Vol}}\approx \frac{1}{N} \sum _{i=1}^{N}\langle \sqrt{\kappa }\delta \varphi \rangle _{i}. \end{equation} \tag{ 21 }$$

More precise treatments of the volume average could be obtained by using space tessellations, such as Delaunay triangulation, but this is too technical and thus beyond the scope of this work. Anyway, using our approximation Equation (21), we find that the bias caused by using particle-number average is at most 1%–2%, which is not unacceptable considering that the sphericity of halos is already an approximation.

Figure 3 shows the profile of $\sqrt{\kappa }\delta \varphi$ inside the dark matter halos. Instead of plotting this halo by halo, we have selected the 80 most massive halos from our simulation box, and divided them into six bins with M/(10¹⁴ M_☉) ∈ [2.0, ∞), [1.0, 2.0], [0.6, 1.0], [0.3, 0.6], [0.2, 0.3], and [0.1, 0.2]. Then, we compute the averaged profile of $\sqrt{\kappa }\delta \varphi$ in each of the six bins (the six curves in Figure 3). As expected, the larger the halo is, the deeper the gravitational potential is and, because $a \sqrt{\kappa }\delta \varphi \approx \Phi$, the larger $|\sqrt{\kappa }\varphi |$ is. Meanwhile, going from the halo center outward one finds that $|\sqrt{ \kappa }\varphi |$ gradually decreases (toward 0), which makes sense because $\sqrt{\kappa }\delta \varphi \rightarrow 0$ as R → ∞. Finally, we also note that the value of $\sqrt{\kappa }\delta \varphi$ depends sensitively on ζ and $|\sqrt{\kappa }\delta \varphi |$ increases with |ζ|.

Zoom In Zoom Out Reset image size

Figure 3 also shows that for our choices of ζ the value of $|\sqrt{\kappa }\delta \varphi |$ is typically of order 10⁻¹⁰to10⁻⁹, which is quite small but gives us no idea how large the fluctuation in φ is. For the latter, we have instead plotted the profile of the quantity δφ/φ inside the 80 halos distributed in six bins. Interestingly, unlike δφ, the quantity δφ/φ does not depend on the sign of ζ, and indeed it almost does not depend on the magnitude of ζ either. This, together with the facts that (1) $a\sqrt{\kappa }\delta \varphi (R,a)\propto \Phi (R,a),$ and (2) the presence of the scalar field and its coupling to matter have negligible effect in the structure formation (so that the halo density profile remains NFW), indicates that the fluctuation of φ in the BSBM model only depends on the non-BSBM physical parameters, and once we have solved the background value of φ, we might gain some knowledge of the $\sqrt{ \kappa }\delta \varphi$ profile without solving it explicitly (see below for further details).

Because the BSBM model is designed for the variation of α, we are also interested in how large this variation, δα/α, is. Because $\alpha \propto e^{2\sqrt{\kappa }\varphi }$, we have $\delta \alpha /\alpha =2\sqrt{\kappa } \delta \varphi$. From Figure 3, we clearly see that this is of order 10⁻¹⁰to10⁻⁹ across galaxy clusters, where the exact value depends on both the size of the halo and the model parameters. Levshakov et al. (2010c) have given an upper limit of the spatial variation of α in Milky Way, |δα/α| < 2 × 10⁻⁷, and the BSBM prediction is well within this bound. Levshakov et al. (2010a, 2010b) also found that in the Milky Way the spatial variation of μ is $\delta \mu /\mu =26\pm 1(\rm {stat})\pm 3(\rm {sys})\times 10^{-9}$, which is comparable to the spatial variation of α predicted by the BSBM model.

Now, given that the NFW profile for dark matter halos is (expected to be) preserved in the varying-α simulations, we wonder whether it is possible to derive some analytical (approximate) formula for the profile of $\sqrt{\kappa }\varphi$. For this let us recall that the NFW profile is expressed as

$$\begin{equation} \frac{\rho (r)}{\rho _{c}} = \frac{\beta }{\frac{r}{R_s}\big(1+\frac{r}{R_s} \big)^2}, \end{equation} \tag{ 22 }$$

where ρ_c is the critical density for matter, β is a dimensionless fitting parameter, and R_s is a second fitting parameter with length dimension. β and R_s are generally different for different halos and should be fitted for individual halos.

We have checked the halos in our analysis and found that the majority of them are indeed very well described by Equation (22), confirming our earlier argument that the coupled scalar-field effect is too tiny to change the structure formation.⁴ However, we shall not use the fitting to Equation (22) in this work, mainly for two reasons: first, the dark matter density profile is in general difficult to measure directly or precisely, while in contrast the circular velocities V_c of the stars rotating about the halo center are easier to measure; second, V_c as a function of radius R is more closely related to $\sqrt{\kappa }\delta \varphi (R)$, which will become clear later. As a result, we shall use fittings to V_c from here on.

Assuming Equation (22) as the density profile and sphericity of halos, we derive V_c easily as

$$\begin{eqnarray} V_c^2(R) &=& \frac{G\!M(R)}{R}\nonumber\\ &=& 4\pi G\beta \rho _cR^3_s\left[\frac{1}{R}\ln \left(1+\frac{R}{R_s}\right) - \frac{1}{R_s+R}\right], \end{eqnarray} \tag{ 23 }$$

where M(R) is the mass enclosed in radius R, and again this is parameterized by β and R_s. From a simulation point of view, it is straightforward to measure M(R) and then fit β and R_s; from an observational viewpoint, it is easy to measure V_c(R), which could again be used to fit β and R_s.

To show how good the fittings could be, we pick out nine halos with different masses and sizes from our simulation box, fit the corresponding β and R_s using the measured M(R), and plot these in Figure 5. As can be seen there, the fitting results (solid curves) agree with the simulation results (crosses) quite well (in particular for the halo with M = 3.895 × 10¹³ M_☉).

To see how this could be related to $\sqrt{\kappa }\varphi$, we remember that in the above we have shown that $a\sqrt{\kappa }\varphi =4\zeta \Phi$, and so all we need to do is to find an expression for Φ(R). For this, we use the fact that the potential inside a spherical halo is given as

$$\begin{equation} \Phi (R) = \int ^{R}_{0}\frac{G\!M(r)}{r^2}dr + C \end{equation} \tag{ 24 }$$

in which GM(r)/r² is the gravitational force and C is a constant to be fixed using the fact that Φ(R = ∞) = Φ_∞, where Φ_∞ is the value of the potential far from the halo.

Using the formula for GM(r)/r² given in Equation (23), it is not difficult to find that

$$\begin{equation*} \int ^{R}_{0}\frac{G\!M(r)}{r^2}dr = 4\pi G\beta \rho _cR^3_s\left[\frac{1}{R_s} -\frac{\ln \big(1+\frac{R}{R_s}\big)}{R}\right] \nonumber \end{equation*}$$

and so

$$\begin{equation} C = \Phi _\infty - 4\pi G\beta \rho _cR^2_s. \end{equation} \tag{ 25 }$$

Then it follows that

$$\begin{equation} \Phi (R) = \Phi _\infty - 4\pi G\beta \rho _c\frac{R^3_s}{R}\ln \left(1+\frac{R }{R_s}\right). \end{equation} \tag{ 26 }$$

If the halo is isolated, then Φ_∞ = 0 and we get

$$\begin{equation} \Phi (R) = {-} 4\pi G\beta \rho _c\frac{R^3_s}{R}\ln \left(1+\frac{R}{R_s} \right). \end{equation} \tag{ 27 }$$

However, in N-body simulations, we have a large number of dark matter halos and no halo is ideally isolated from the others. In such situations, Φ_∞ in Equation (26) should be replaced by

$$\begin{equation} \Phi _\ast \ \equiv \ \Phi (R=R_\ast \gg R_{\mathrm{vir}}) \ne 0, \end{equation} \tag{ 28 }$$

where R_* is some radius large compared with R_vir (the virialized halo radius) but small compared with inter-halo distances. Then we get

$$\begin{equation} a\sqrt{\kappa }\delta \varphi (R) = 4\zeta \left[\Phi _\ast - 4\pi G\beta \rho _c \frac{R^3_s}{R}\ln \left(1+\frac{R}{R_s}\right)\right].\ \ \ \end{equation} \tag{ 29 }$$

As an example to show how well Equation (29) works, we show in Figure 6 the results of $a\sqrt{\kappa } \delta \varphi (R)$ for the same halos used to fit β and R_s in Figure 5. Here, the crosses represent the values of $a\sqrt{\kappa }\delta \varphi$ measured from the N-body simulations and the curves represent our analytical approximations, in which the solid curve is obtained by setting Φ_* = Φ_∞ = 0, while the dashed curve is from tuning Φ_* appropriately. Obviously in most cases, there is a (nearly) constant shift of the approximation with respect to the numerical results, which accounts for Φ_* being non-zero.

Note that Equation (29) captures the shapes for $a\sqrt{\kappa } \delta \varphi$ in various halos, but there is still one free parameter, Φ_*, to be tuned to match the numerical results. This parameter summarizes our lack of knowledge about the environment in which the considered halo resides. As a result, formula (29) is most suitable for application in isolated halos, while for residential halos some extra work remains to be done to make it accurate.

Alternatively, one could also consider Equation (29) as a three-parameter parameterization of $a\sqrt{\kappa }\delta \varphi$, for which the three parameters β, R_s, and Φ_* could be fitted using the results from N-body simulations. Given all the above results, we expect that this could produce some nice fitting curves too, but we do not expand on this point here.