© The Authors 2025

1. Introduction

The standard model of cosmology, referred to as ΛCDM (Planck Collaboration VI 2018), includes two dark components: cold dark matter (CDM) and dark energy in the form of a cosmological constant (Λ). Despite immense investigatory efforts, it remains unclear what the physical nature of these components may be. In the case of dark matter, many candidate particles have been proposed over the years (see Arbey & Mahmoudi 2021 for a recent review), as well as the possibility that the phenomenology of dark matter is in fact explained by a modification of gravity in the weak-field limit (for a recent review of this approach, see Banik & Zhao 2022). In the case of dark energy, the standard model assumes this to be represented by a cosmological constant; however, there are some very recent indications that dark energy may be dynamical and not in fact a cosmological constant (DES Collaboration 2025; DESI Collaboration 2025). Among the multitude of dynamical dark energy models, perhaps the most straightforward is that of quintessence (Wetterich 1995; Caldwell et al. 1998), where a scalar field drives the late-time accelerated expansion.

An interesting extension of these models introduces a coupling between the two dark sector components (Amendola 2000; Wang et al. 2016). There is also a multitude of possible forms that this coupling could take, with an appealing possibility being a momentum coupling between the dark matter component and the quintessence scalar field (Pourtsidou et al. 2013). Several studies have investigated momentum-coupled models, sometimes referred to as dark scattering models (Simpson 2010), through the lens of numerical simulations (Baldi & Simpson 2015, 2017). In our previous work reported in Palma & Candlish (2023), a specific class of quintessence scalar field dark energy models with a momentum coupling to dark matter was investigated for the case of a positive coupling. In that study it was found that the coupling enhanced structure formation at larger scales, as measured by the power spectrum, but significantly suppressed structure on small scales, and reduced the slope of the inner density profiles of large mass halos. Furthermore, the velocities of the substructure and the particle content of the halos was found to be significantly enhanced by the coupling. The evolution of the linear structure in these models has been studied in Pourtsidou & Tram (2016), Chamings et al. (2020), Spurio Mancini & Pourtsidou (2022) with a view to resolving the σ₈ tension (Ade et al. 2014).

In this work we extend our previous study in Palma & Candlish (2023) by considering higher resolution simulations, allowing us to analyse a much larger population of dark matter halos. In addition, we study the effect of a negative coupling between the dark sector components, something that was absent in our previous study. The particular focus of this study is to explore statistically the effects of the coupling on halos and their substructure, motivated specifically by an interesting effect that was found in our earlier work whereby the velocity distributions of the most massive halos exhibited a bimodal behaviour, suggestive of a dynamically disturbed state, as is expected in cluster mergers. Given that the comparison simulations without the presence of the coupling showed that the same structure had fully virialised by redshift zero, our result suggested that the dynamical states of clusters may be significantly different in coupled models. However, the strength of this conclusion was severely limited by the low number statistics available due to the low resolution of our previous simulations. For this work we use higher resolution simulations to undertake a statistical study of this phenomenon.

2. Theory

We refer to our previous work Palma & Candlish (2023) and the original theoretical papers (Pourtsidou et al. 2013; Skordis et al. 2015) for the full details of the theoretical background of the model we consider in this work. For a self-contained work, we here summarise the theoretical background to our study.

The energy-momentum tensor of the scalar field for Type 3 models is written as

$$\begin{array}{c}\hfill T_{\mu \nu }^{(\varphi )}=F_Y{\varphi }_{\mu }{\varphi }_{\nu }-F{g}_{\mu \nu }-Z{F}_Z{u}_{\mu }{u}_{\nu },\end{array}$$(1)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Y& =\frac{1}{2}{\varphi }_{\mu }{\varphi }^{\mu }\hfill \\ \hfill Z& =u^{\mu }{\varphi }_{\mu },\hfill \end{array}\end{array}$$(2)

and F = F(Y, Z, ϕ) is some function. F_Y and F_Z denote derivatives of this function with respect to Y and Z, respectively. The four-velocity of the dark matter perfect fluid is given by u^μ and ϕ_μ ≡ ∂_μϕ. The equation of motion for the scalar field is given by

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\mu }(F_Y{\varphi }^{\mu }+F_Z{u}^{\mu })-F_{\varphi }=0,\end{array}$$(3)

while the dark matter fluid satisfies the usual equation for the conservation of density. The momentum transfer equation is given by

$$\begin{array}{c}\hfill (\rho -Z{F}_Z)u^{\beta }{\mathrm{\nabla }}_{\beta }{u}_{\mu }={\mathrm{\nabla }}_{\beta }(F_Z{u}^{\beta }){\stackrel{\sim }{\varphi }}_{\mu }+F_Z{D}_{\mu }Z,\end{array}$$(4)

where $D_{\mu }=q_{\mu }^{\nu }{\mathrm{\nabla }}_{\nu }$ is the spatial derivative operator given in terms of the projection operator $q_{\mu }^{\nu }\equiv u_{\mu }{u}^{\nu }+{\delta }_{\mu }^{\nu }$, and ${\stackrel{\sim }{\varphi }}_{\mu }=q_{\mu }^{\nu }{\mathrm{\nabla }}_{\nu }\varphi =D_{\mu }\varphi ={\partial }_{\mu }\varphi +u^{\nu }{u}_{\mu }{\partial }_{\nu }\varphi$ is the spatial projection of the derivative of the scalar field. We now specialise to the case of coupled quintessence by choosing F = Y + V(ϕ)+γ(Z). We work in the Newtonian gauge to facilitate passing to the non-relativistic limit. The line element in this gauge is given by

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=a^2(\tau )[-(1+2\mathrm{\Psi })\mathrm{d}{\tau }^2+(1+2\mathrm{\Phi }){\delta }_{\mathit{ij}}\mathrm{d}{x}^i\mathrm{d}{x}^j],\end{array}$$(5)

where Φ and Ψ are spatial scalars and δ_ij is the three-dimensional Kronecker delta. The perturbed fluid four-velocities (to linear order) in this gauge are

$$\begin{array}{cc}\hfill u_0& =-a(1+\mathrm{\Psi }),\hfill \\ \hfill u_i& =a{v}_i.\hfill \end{array}$$

The evolution of the cold dark matter fluid at the background level and of the CDM fluid perturbations at linear level are given by the standard equations. However, the background evolution of the scalar field depends on the coupling and is given by Equation (3) as

$$\begin{array}{c}\hfill \ddot{\varphi }-{\gamma }_{\mathit{ZZ}}\ddot{\varphi }+2\mathcal{H}\dot{\varphi }+{\gamma }_{\mathit{ZZ}}\mathcal{H}\dot{\varphi }-3a\mathcal{H}{\gamma }_Z+a^2{V}_{\varphi }=0.\end{array}$$(6)

From Equation (3) at first order in the perturbations we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill V_{\varphi \varphi }\phi a^2& +3{\gamma }_Z\mathrm{\Psi }a\mathcal{H}+3{\gamma }_{Z}a\dot{\mathrm{\Phi }}-{\gamma }_{Z}a{\mathrm{\nabla }}^2\theta -{\gamma }_{\mathit{ZZZ}}\frac{\ddot{\overline{\varphi }}}{a}\dot{\phi }\hfill \\ \hfill & +{\gamma }_{\mathit{ZZZ}}\frac{\dot{\overline{\varphi }}}{a}\dot{\phi }\mathcal{H}+2{\gamma }_{\mathit{ZZ}}\mathrm{\Psi }\ddot{\overline{\varphi }}-2{\gamma }_{\mathit{ZZ}}\mathrm{\Psi }\dot{\overline{\varphi }}\mathcal{H}+{\gamma }_{\mathit{ZZ}}\dot{\mathrm{\Psi }}\dot{\overline{\varphi }}\hfill \\ \hfill & -{\gamma }_{\mathit{ZZ}}\ddot{\phi }-2{\gamma }_{\mathit{ZZ}}\dot{\phi }\mathcal{H}-2\mathrm{\Psi }\ddot{\overline{\varphi }}-4\mathrm{\Psi }\dot{\overline{\varphi }}\mathcal{H}-3\dot{\mathrm{\Phi }}\dot{\overline{\varphi }}-\dot{\mathrm{\Psi }}\dot{\overline{\varphi }}+\ddot{\phi }\hfill \\ \hfill & -{\mathrm{\nabla }}^2\phi +2\dot{\phi }\mathcal{H}=0.\hfill \end{array}\end{array}$$(7)

Passing to Fourier space and taking the Newtonian limit (non-relativistic velocities are implicit in the gauge choice as the DM fluid velocity perturbation is assumed to satisfy |v|≪1) of k ≫ ℋ, this equation simplifies drastically to

$$\begin{array}{c}\hfill \phi =a{\gamma }_z\theta .\end{array}$$(8)

With this we write the momentum transfer Equation (4) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{\theta }+\mathcal{H}\theta +\mathrm{\Psi }=& \frac{1}{a\overline{\rho }-{\gamma }_Z\dot{\overline{\varphi }}}[2{\gamma }_Z\dot{\overline{\varphi }}\theta \mathcal{H}+3a{\gamma }_Z^2\theta \mathcal{H}-{\gamma }_Z\mathrm{\Psi }\dot{\overline{\varphi }}\hfill \\ \hfill & +{\gamma }_Z\ddot{\overline{\varphi }}\theta +a{\gamma }_Z^2\mathcal{H}\theta +a{\gamma }_Z{\gamma }_{\mathit{ZZ}}\dot{\overline{Z}}\theta +a{\gamma }_Z^2\dot{\theta }]\hfill \\ \hfill -& \frac{1}{a^2\overline{\rho }-a{\gamma }_Z\dot{\overline{\varphi }}}[{\gamma }_{\mathit{ZZ}}{\dot{\overline{\varphi }}}^2\theta \mathcal{H}-\mathcal{H}{\gamma }_Z{\gamma }_{\mathit{ZZ}}\dot{\overline{\varphi }}\theta \hfill \\ \hfill & +{\gamma }_{\mathit{ZZ}}\dot{\overline{\varphi }}\ddot{\overline{\varphi }}\theta +{\gamma }_Z{\gamma }_{\mathit{ZZ}}\dot{\overline{\varphi }}\theta ],\hfill \end{array}\end{array}$$(9)

where $\dot{\overline{Z}}=1/a(\ddot{\varphi }+\mathcal{H}\dot{\varphi })$. We now must choose a form for the coupling, and so we follow Pourtsidou et al. (2013) and choose

$$\begin{array}{c}\hfill \gamma (Z)={\gamma }_0{Z}^2,\end{array}$$(10)

where γ₀ is a constant whose value is assumed to be in the range 0 ≤ γ₀ < 1/2. The equation (9) becomes

$$\begin{array}{c}\hfill (1+h_1){\dot{v}}_i+(1+h_2)\mathcal{H}{v}_i+(1+h_3){\mathrm{\nabla }}_i\mathrm{\Psi }=0\end{array}$$(11)

where the coefficients h₁, h₂, and h₃ are

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill h_1& =\frac{4{\gamma }_0^2{\dot{\varphi }}^2}{a^2\rho -2{\gamma }_0\dot{{\varphi }^2}},\hfill \\ \hfill h_2& =\frac{(8{\gamma }_0^2-2{\gamma }_0){\dot{\varphi }}^2+(8{\gamma }_0^2-4{\gamma }_0)\dot{\varphi }\ddot{\varphi }\frac{1}{\mathcal{H}}}{a^2\rho -2{\gamma }_0{\dot{\varphi }}^2},\hfill \\ \hfill h_3& =\frac{2{\gamma }_0{\dot{\varphi }}^2}{a^2\rho -2{\gamma }_0{\dot{\varphi }}^2}.\hfill \end{array}\end{array}$$(12)

In the h₂ term, we can replace $\ddot{\varphi }$ using the evolution equation of the background field, given by Equation (6), with

$$\begin{array}{c}\hfill \ddot{\varphi }(1-2{\gamma }_0)+2\mathcal{H}\dot{\varphi }(1-2{\gamma }_0)+a^2{V}_{\varphi }=0,\end{array}$$(13)

where we have used Equation (10). Thus, h₂ can be written as

$$\begin{array}{c}\hfill h_2=\frac{4{\gamma }_0(\frac{3}{2}-2{\gamma }_0){\dot{\varphi }}^2+4{\gamma }_0\dot{\varphi }{a}^2{V}_{\varphi }/\mathcal{H}}{a^2\rho -2{\gamma }_0{\dot{\varphi }}^2}.\end{array}$$(14)

Now that we have the modified Euler equation, we must now implement this in our cosmological N-body code. We use the well-known RAMSES code (Teyssier 2002), which is a grid-based adaptive mesh refinement (AMR) code, using a particle-mesh (PM) scheme to evolve the dark matter particle distribution. Our implementation begins with revisiting the supercomoving coordinates (Martel & Shapiro 1998) used in RAMSES, which are defined as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathit{v}& =H_{0}L\frac{1}{a}\stackrel{\sim }{\mathit{u}},\hfill \\ \hfill \mathit{x}& =\frac{1}{a}\frac{\stackrel{\sim }{\mathit{x}}}{L},\hfill \\ \hfill \mathrm{d}t& =a^2\frac{d\stackrel{\sim }{t}}{H_0},\hfill \\ \hfill \mathrm{\Psi }& =\frac{L^2{H}_0^2}{a^2}\stackrel{\sim }{\mathrm{\Phi }},\hfill \end{array}\end{array}$$(15)

where L is the length of the simulation box. The coordinates denoted with a tilde are the supercomoving coordinates. To simplify the notation, we apply the transformation and then remove the tildes. Thus, using Eq. (15) in Eq. (11) we get

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathit{u}}{\mathrm{d}t}=-\frac{h_2-h_1}{1+h_1}{a}^2\frac{H}{H_0}\mathit{u}-\frac{1+h_3}{1+h_1}\mathit{\nabla }\mathrm{\Phi }.\end{array}$$(16)

We now have the modified Euler Eq. (16) in a form in which it may be discretised and solved numerically. To connect with the implementation of the modified Euler equation in the code, we write the finite difference update of the velocity, as implemented in RAMSES using the Leapfrog scheme, in the following manner:

$$\begin{array}{c}\hfill \frac{v_p^{n+1/2}-v_p^n}{(1/2)\mathrm{\Delta }{t}^n}=F.\end{array}$$(17)

Here F is the force acting on the particle. We can now modify the velocity update as required to implement Equation (16) in the following way:

$$\begin{array}{c}\hfill v_p^{n+1/2}=v_p^n-\frac{h_2-h_1}{1+h_1}{a}^2\frac{H}{H_0}{v}_p^n\mathrm{\Delta }{t}^n/2+\frac{1+h_3}{1+h_1}F\mathrm{\Delta }{t}^n/2.\end{array}$$(18)

We define two new coefficients, ϵ₁ and ϵ₂, to simplify the expression:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {ϵ}_1& =1-\frac{h_2-h_1}{1+h_1}{a}^2\frac{H}{H_0}\mathrm{\Delta }{t}^n/2,\hfill \\ \hfill {ϵ}_2& =\frac{1+h_3}{1+h_1}.\hfill \end{array}\end{array}$$(19)

Thus, finally Eq. (18) becomes

$$\begin{array}{c}\hfill v_p^{n+1}={ϵ}_1{v}_p^n+{ϵ}_{2}F\mathrm{\Delta }{t}^n/2.\end{array}$$(20)

This is the equation we have implemented in RAMSES. The standard dynamics is recovered by setting ϵ₁ = ϵ₂ = 1, which is equivalent to having all the h_i equal to zero.

3. Simulations

We focus on the same scalar field model, referred to as model C, as used in our previous work. This model has a potential of the form of Albrecht & Skordis (2000), given by

$$\begin{array}{c}\hfill V(\varphi )=({(\varphi -\beta )}^{\alpha }+\mathrm{\Gamma })e^{-\lambda \varphi },\end{array}$$(21)

where the parameters are chosen to be β = 3.8, Γ = 20.0, α = 17.0, and λ = 0.934. The last value is fixed by demanding that the closure condition of the energy densities is satisfied (assuming spatial flatness). The initial values of the scalar field and its derivative are ϕ = 100, $\dot{\varphi }=10$. In Albrecht & Skordis (2000) the form of this potential is justified as providing a model where all parameters are roughly of order unity in Planck units, that inherits useful properties from the well-known tracker solutions of exponential models, and that deviates from a pure exponential in the correct manner to provide the appropriate late-time acceleration. Furthermore, potentials of this form are expected to emerge naturally in the low energy limit of string theory-inspired models. For our purposes, this model was chosen primarily based on the phenomenological constraint of closely matching to the Hubble constant evolution of ΛCDM in order to reproduce the observational successes of that model.

We generated the initial conditions for our simulations using the MonofonIC code (Hahn et al. 2020) (sometimes referred to as MUSIC2). This code incorporates the CLASS code for the generation of the initial transfer function. We modified the CLASS code used in order to incorporate the coupling between the dark matter and the scalar field in order to generate consistent initial conditions, although the effect of the coupling is minimal at high redshift. The cosmological parameters that we used are those of Planck Collaboration VI (2018), with the obvious exception that we did not use Ω_Λ but a quintessence field. The parameters are Ω_b = 0.0494, Ω_CDM = 0.264979, Ω_ϕ = 0.68412, h = 0.67321, σ₈ = 0.8102, and n_s = 0.9661.

We note that the presence of the coupling leads to an overall scaling of the background potential for the scalar field:

$$\begin{array}{c}\hfill \ddot{\varphi }+2\mathcal{H}\dot{\varphi }+\frac{1}{1-2{\gamma }_0}{a}^2{V}_{\varphi }=0.\end{array}$$(22)

Therefore, the background evolution is modified by changing the value of the coupling constant. In order to separate the effects of the coupling on structure formation due to the modified dark matter dynamics as compared to simply a modified Hubble parameter, we ran simulations (as in our previous study) where the background evolution is modified by the coupling, but the Euler equation is fixed to the standard form (i.e. the coupling is switched off for the particle dynamics).

The background evolution of the Hubble parameter for the coupled models is given in Fig. 1, as compared to ΛCDM. It is apparent that there is a more pronounced and longer deviation from ΛCDM seen in the case of positive coupling, while for small negative coupling the deviation is much reduced, and for large negative coupling it is very close to ΛCDM. It should be noted that the positive coupling model is affected by a strong coupling problem at γ₀ = 1/2, and so it is to be expected that values of γ₀ closer to that critical value exhibit more pronounced effects. In any case, the deviations from the standard model do not exceed ∼4%. In the case of a negative coupling there is no theoretical restriction in the allowed value of the coupling constant.

Fig. 1. Hubble parameter normalised by that of ΛCDM for the three values of the coupling constant.

We also show the evolution of the parameters c₁ and c₂ in the modified Euler equation in Fig. 2. We can see that a positive coupling leads to a significant enhancement of the gravitational acceleration term c₂, which we would naively expect to lead to an increase in small-scale structure in this model. As we have shown before, this is in fact not the case. The reason is that the cosmological friction term c₁ for positive coupling leads to an effective push (which we call the cosmological push, as opposed to cosmological friction) where the particles are further accelerated beyond what is expected from the gravitational effects alone. For negative coupling, we have the opposite behaviour: the cosmological friction is now enhanced, while the gravitational acceleration is reduced. We note that the change in the effective gravitational force (the c₂ term) is not as pronounced as the change in the effective cosmological friction (the c₁ term) when comparing the difference between the two negatively coupled models. The c₂ coefficient is approximately 12% smaller for γ₀ = −20 as compared to γ₀ = −0.3, whereas the c₁ coefficient is approximately 35% larger for γ₀ = −20 as compared to γ₀ = −0.3. Again we might naively expect to see a reduction in structure due to the increased relevance of the cosmological friction term, but in fact we see, at small scales, an enhancement of structure.

Fig. 2. Modified Euler equation coefficients for the three values of the coupling constant.

All of our simulations use 512³ particles in a 128³ Mpc h⁻¹ volume, leading to a particle mass of order 1.3 × 10⁹ M_⊙ h⁻¹. We set the minimum particle number for a dark matter halo at 20, and thus we do not resolve any halos less massive than approximately 2.6 × 10¹⁰ M_⊙ h⁻¹. The minimum grid resolution is 250 kpc h⁻¹ and the maximum refined resolution is 7.8 kpc h⁻¹. The full list of simulations that we use is given in Table 1.

4. Results

4.1. Density projections

As a first simple check of our results, we confirmed that the overall density field of the structure formed in our simulations is visually identical. We did this using simple density projection plots where the particles are binned into a two-dimensional pixel mesh of 512 × 512 pixels. Each particle contributes in the same way to the total density within each pixel (i.e. the density calculation is unweighted). These results are shown in Fig. 3. It is clear, at the qualitative level of a visual inspection of these density plots, that all models generate structure in mostly the same locations within the computational box, although perhaps there is some indication of differing degrees of density contrast, as indicated by the colour associated with the value of the projected particle density within that pixel. This is to be expected as the coupling affects the degree of structure formation that takes place, as we see from the examination of the power spectra of our simulations in the next section.

Fig. 3. Density projection plots of all simulations considered in this work. The leftmost panel is our ΛCDM run, and the other three panels (from left to right) are coupled models with: γ0 = 0.3, γ0 = −0.3, γ0 = −20.

4.2. Power spectrum

Then we analysed the power spectra of our coupled models, normalised with respect to the associated uncoupled model. In other words, we divided the power spectrum of model p0.3 by that of p0.3_u, the power spectrum of model m0.3 by that of m0.3_u, and finally the power spectrum of m20.0 by that of m20.0_u. This was to extract the effect of the coupling from the modified background evolution in each model. In all cases the power spectra were calculated using the POWMES code (Colombi et al. 2009).

In Fig. 4 we plot the normalised power spectra. We can see that the positive coupling model behaves, as expected, in the same manner as we found in our earlier analysis described in Palma & Candlish (2023): on large linear scales we see a mostly scale-independent increase in power, while at small non-linear scales there is a pronounced (scale-dependent) reduction in power compared to the uncoupled model. In our previous study the loss of small-scale power was understood to be due to the significantly modified dynamics of the dark matter whereby high velocities (such as those found within a major overdensity) are further enhanced by what we refer to as the cosmological push. This push removes dark matter from the overdense regions, thus reducing structure. An interesting aspect of this phenomenon is that the modification of the velocity vector due to this effect is proportional to the Hubble parameter and the velocity vector itself, rather than being proportional to the gravitational potential gradient. Thus, the dynamics induced by the cosmological push are less directly connected to the density distribution, leading to late-time dynamical states of the substructure that exhibit less virialisation than is expected in a non-coupled model. We explore this much more fully in the following sections.

Fig. 4. Power spectra normalised by the corresponding uncoupled model.

In the case of a negatively coupled model we see opposite behaviour: the large-scale linear power spectrum is slightly suppressed compared to the coupled model, while the small-scale non-linear power spectrum shows a strong enhancement of power, particularly for the γ₀ = −20 model, although it is a less pronounced effect than the reduction in power shown in the positively coupled model. Inverting the argument given above for positive coupling, we expect that the enhanced cosmological friction acts to inhibit structure formation in a scale-independent manner, thus leading to large-scale suppression of structure formation. At small scales, within overdensities, the particle velocities are enhanced, leading to a more pronounced friction effect, which (despite the reduced gravitational acceleration) is sufficient to reduce velocities such that more bound virialised structures are formed. We note that the positive coupling model shows more pronounced deviations from the uncoupled model due to the coupling γ₀ = 0.3 being close to the strong coupling value of γ₀ = 1/2 where singular behaviour arises in the equations of motion.

4.3. Halo mass function and velocity dispersions

The halo catalogues from our simulation runs were generated by using the Amiga Halo Finder (AHF) (Knollmann & Knebe 2009). This uses an overdensity criterion to identify which particles belong to which halos. In addition, the escape velocity of each particle is calculated and multiplied by a factor of 1.5 to determine whether the particle is actually gravitationally bound to the halo. We note that in our analysis in our prior study (Palma & Candlish 2023), we modified the effective Newtonian gravitational constant used in the calculation of the escape velocity according to the value of the c₂ coefficient. We did this again for this study, where the effective values of the gravitational constant (at z = 0) are given in Table 2. We also considered how the analysis is affected if we assume the Newtonian value of the gravitational constant, supposing that in an observational context we would not be privy to the effective gravitational constant experienced by the dark matter. For the halo mass function, however, the difference between using these effective gravitational constants and using the standard Newtonian value are negligible. We used the conventional value of Δ₂₀₀ (i.e. 200 times the background density) to define the maximum radius of our halos.

The halo mass functions for our simulations are shown in Fig. 5, with a comparison to the Watson et al. (2013) fitting line extracted from a suite of very large particle number N-body simulations. We can see that the halo counts in the intermediate mass bins compare reasonably well with the reference line, except at the low and high mass extremes, as expected, where insufficient mass resolution and insufficient box size, respectively, limit our results. It is worth noting that the presence or otherwise of the coupling, and the sign of the coupling, have a negligible impact on the halo mass function.

Fig. 5. Halo mass functions for all simulations considered in this study. The comparison line given is from Watson et al. (2013), as defined in the Colossus Python package.

To begin our analysis of the velocities of the particles within the dark matter halos, we plot in Fig. 6 the velocity dispersions, as measured by AHF, for halos of all masses found in the simulations. Within each mass bin we determine the standard deviation of the measured velocity dispersions, which are indicated by the vertical error bars. We immediately see that the positive coupling model exhibits a clear increase in the velocity dispersions across all mass bins, as compared to all other models. The negatively coupled models show slightly suppressed velocity dispersions as compared to ΛCDM but within the errors they are identical. The uncoupled models (not plotted) are also identical to ΛCDM within the errors. We note that the velocity dispersions are also not significantly sensitive to the use of the effective gravitational constant as opposed to the Newtonian value.

Fig. 6. Velocity dispersion as a function of mass for all halos in all the coupled models considered in this study and the reference ΛCDM model.

4.4. Bimodality in the velocity distributions

4.4.1. Bimodality parameter

As discussed earlier, in previous work we found a strikingly bimodal velocity distribution for the most massive halo in the positively coupled model, indicating a halo that has yet to reach a fully virialised state. This was, however, an effect seen in only a single halo. For this work we studied this phenomenon in much more detail, and from a statistical perspective given the much higher number of halos available in our present simulations.

We analysed the degree of bimodality in the following manner. Firstly, we separate the halos into three categories: host, isolated, or subhalo. A host halo is simply any halo that contains at least one other halo, as determined by AHF. This category thus includes large clusters as well as small groups where a single subhalo is within a host. A subhalo is a halo that is contained within another halo. Finally, an isolated halo is a halo that is neither a host nor a subhalo.

After this categorisation, for each of the 500 most massive halos within each category in our simulations, we generate the histogram of 3D velocity magnitudes (normalised by each halo’s velocity dispersion) for all particles within the halo, using 60 evenly spaced bins. We then fit to that histogram the following bimodal Gaussian distribution:

$$\begin{array}{c}\hfill f(v)=A_{1}exp(-\frac{{(v-{\mu }_1)}^2}{2{\sigma }_1})+A_{2}exp(-\frac{{(v-{\mu }_2)}^2}{2{\sigma }_2}).\end{array}$$(23)

Here v denotes the 3D velocity magnitude, the μ_i values are the means of each component, the A_i values are the amplitudes, and the σ_i values are the standard deviations. We compare the values of A₁ and A₂ to determine which of the two Gaussian components has the larger amplitude. The parameters of the larger Gaussian is referred to as A_L, μ_L, and σ_L, while the parameters of the smaller Gaussian will be referred to as A_S, μ_S, and σ_S.

We then define the following bimodality parameter:

$$\begin{array}{c}\hfill B=\frac{A_S}{A_L}\frac{|{\mu }_S-{\mu }_L|}{{\sigma }_L}.\end{array}$$(24)

The motivation for this parameter is to quantify the degree of bimodality of the velocity magnitude histograms. The factor of the amplitude ratios will reduce the value of this parameter towards zero if the smaller Gaussian is significantly smaller than the larger, as in this case the smaller Gaussian would be essentially irrelevant to the distribution. The absolute value of the difference between the mean values is a straightforward measure of the bimodality of the double Gaussian distribution, which we measure in units of the standard deviation of the larger Gaussian. We find that this parameter is an effective measure of the bimodality in almost all cases, except for the smallest subhalos, where the relatively low particle number leads to noisy histograms, which causes some of the larger amplitude Gaussian components to be assigned very small standard deviations, thus leading to very high bimodality parameter values. For the subhalos we therefore limit the histograms to values of B < 5.

In Fig. 7 we plot the histograms of the bimodality parameters for the 500 most massive host halos in each simulation. Here we see a clear distinction between the positive coupling model (blue) and all other models, where the host halos exhibit a clear bimodality, and relatively infrequent unimodal distributions. This confirms the preliminary results found in our previous study (Palma & Candlish 2023). The bimodality found in a single massive host halo in that study has been found to be prevalent throughout the host halo population of the positively coupled model in the present work.

Fig. 7. Histogram of bimodality parameter values calculated from the bi-Gaussian distribution fitted to the velocity magnitude histograms of all particle velocities in the 500 most massive host halos of each simulation.

The uncoupled models and ΛCDM all show broadly similar behaviour, as does the small negative coupling model, with a preference towards unimodal distributions (small parameter values) and a tail of bimodality. Interestingly, however, in the model with a large negative coupling (green) the preference towards unimodal velocity distributions is considerably stronger, with bimodality far less frequent when compared to the uncoupled models (or ΛCDM). This suggests that, statistically, halos in this model are more virialised than either the uncoupled models or ΛCDM.

Now we turn our attention to the isolated halos in Fig. 8. We can see in this case that for most of the models, particularly ΛCDM and the small negatively coupled model, bimodality is relatively scarce in these halos; instead, we see a statistical preference towards unimodal distributions. This is to be expected given that these are isolated halos that either have not undergone any merger or have not undergone a recent merger such that any subhalo has fully dissolved within the host. Interestingly, this contrasts with the positive coupling model where there is still an increased probability of finding a bimodal distribution, even though these halos are isolated. This would suggest that the material stripped during some past merger is still undergoing virialisation. In addition, the large negative coupling model shows a slightly reduced probability for bimodality and an increased likelihood of unimodal distributions.

Fig. 8. As in Fig. 7, but for isolated halos.

Finally, the subhalos are given in Fig. 9. Given that these halos are undergoing tidal interactions within a host halo, and are likely to be far from virialised, we thus expect some degree of bimodality, and we do indeed see a larger tail of higher bimodality parameters in this case. This behaviour is, however, clearly independent of the coupling.

Fig. 9. As in Fig. 7, but for subhalos.

In the previous analysis of the bimodality, all halos are calculated using the effective gravitational constants as listed in Table 2. Clearly we know what these values must be from our simulations. However, it is worthwhile to explore how these results are modified if we instead assume a position of ignorance and work with the standard Newtonian value of G, as would be the case in a hypothetical observational analysis.

In Figs. 10, 11, and 12 we compare the coupled models with γ₀ = 0.3 and γ₀ = −20 using the modified gravitational constant and those models using the Newtonian value. We also include the ΛCDM results for reference. By using the Newtonian gravitational constant in the halo finding procedure for these models we are modifying which particles are considered bound to each halo. In the positive coupling case the effective gravitational constant is higher than the Newtonian value, thus using the lower Newtonian value in the halo finding process could lead to the unbinding of high velocity particles which are otherwise included. In the case of a negative coupling, the effective gravitational constant is lower, thus using the higher Newtonian value may lead to unbound high velocity particles as being considered bound.

Fig. 10. As in Fig. 7, but using the Newtonian gravitational constant in the halo finding process.

Fig. 11. As in Fig. 8, but using the Newtonian gravitational constant in the halo finding process.

Fig. 12. As in Fig. 9, but using the Newtonian gravitational constant in the halo finding process.

In Fig. 10 we compare only the host halos, which, as we have seen, exhibit the greatest degree of bimodality. We see that the histogram for the positively coupled model using the Newtonian value is identical to that seen by using the effective gravitational constant. Thus, any unbound particles have no significant effect on the degree of bimodality. In the negatively coupled case we see that the unimodality of the velocity distribution is even more pronounced when using the Newtonian gravitational constant (purple line) as opposed to the effective gravitational constant (orange line). For the isolated halos in Fig. 11 we see again that the bimodality of the positively coupled model is unaffected by the choice of gravitational constant. However, The negatively coupled model again shows some difference: the use of the Newtonian gravitational constant leads to a slight reduction in the unimodality of these halos. In either case, the isolated halos in the negatively coupled model exhibit more unimodality than bimodality. Finally, in Fig. 12 we see that the subhalo distributions are essentially unaffected by the gravitational constant used.

4.4.2. Probability of bimodality

We further quantify the degree of unimodality or bimodality by determining the probabilities of observing certain ranges of bimodality parameters for our models. We limit the bimodality distributions for all halo types to B = 3. We then calculate the areas under our histograms for certain ranges of the bimodality parameter. The first range we consider is values of B from 0 to 0.5, which we consider to be indicative of a unimodal distribution. The other range we consider is values of B from 1 to 2, which we consider to be indicative of a bimodal distribution. Furthermore, to connect our work more closely with observationally accessible quantities, we consider 1D velocity distributions by repeating our bimodal Gaussian fitting procedure with only the v_x component of the particle velocity vectors.

The results of the 3D analysis are shown in Figs. 13, 14, while the probabilities obtained from the 1D bimodality parameter histograms are given in Figs. 15 and 16. The coupled models are separated from the uncoupled models and ΛCDM by a vertical blue line to aid in reading these plots.

Fig. 13. Probabilities of unimodality in 3D velocity distributions for all models.

Fig. 14. Probabilities of bimodality in 3D velocity distributions for all models.

Fig. 15. Probabilities of unimodality in 1D velocity distributions for all models.

Fig. 16. Probabilities of bimodality in 1D velocity distributions for all models.

For the case of unimodality in the 3D velocity distributions (Fig. 13) we see for the host halos that there is a clear separation as we move from the positively coupled model (with either the effective or Newtonian gravitational constant) to the small negative coupling, large negative coupling, and then large negative coupling with the Newtonian gravitational constant. The probability of unimodal host distributions is below 15% for the positively coupled model, whereas it is around 35% for the large negatively coupled model (or 45% if using a Newtonian analysis). For the isolated halos, there is a higher probability for unimodality in negatively coupled models, around 40%, compared to 25% for the positively coupled model. The subhalo unimodality probability is roughly 25%−30% for all coupled models. In the uncoupled cases and ΛCDM the host halos have a unimodal probability of ∼25%, similarly for subhalos, while the isolated halos have a somewhat higher probability of unimodality, ∼35%.

In Fig. 14 we consider the probability of bimodality. First we look at host halos. For the positively coupled model this probability is high, at ∼50%, whereas for the large negative coupling model the probability is low, around 15% (or below 10% in a Newtonian analysis). Isolated halos in the positively coupled model have ∼35% probability of bimodality, while this is far lower for isolated halos in all the negatively coupled models, being ∼10%. For the uncoupled and ΛCDM models the hosts and subhalos have ∼30% probability of bimodality, whereas the isolated halos have a probability of ∼20% (with rather large variation across uncoupled models).

The trends clearly seen in the 3D velocity distributions become much less clear when using 1D velocity distributions, as expected. In the case of unimodality (Fig. 15) for the host halos, the positively coupled model, using the Newtonian analysis, has a probability of 65%, whereas it is ∼78% using the effective gravitational constant. For the negatively coupled models, the uncoupled models and ΛCDM the probability of unimodality is around 80%−85%. The isolated halos all show a similar probability, regardless of coupling. The subhalos for all models show a probability of a unimodal distribution of around 65%−70%. Thus, for the 1D velocity distributions only the positively coupled model using a Newtonian analysis is clearly discriminated.

For the case of bimodality, again only the positively coupled model in a Newtonian analysis shows a clear difference, with a bimodality probability in the host halos of ∼10%, whereas it lies at 6% or below for all other models. Isolated halos and subhalos appear to offer no means for model discrimination in the case of 1D velocities.

4.5. Halo properties

4.5.1. Mass, radius, and concentration

The halo masses, radii, and concentrations are compared across the entire halo catalogues in Figs. 17, 18, and 19. In these plots we sort the halo catalogues of each model according to each property and then plot these sorted values together, with the sorting index referring simply to the halo number according to the chosen ordering. Due to differences in the number of halos in each model, the curves do not necessarily overlap at the extreme right of the plots, where we see a shift in the horizontal position of the curves. If the distributions of properties are broadly similar, however, we see a significant overlap of the curves for all index values except the highest and, in particular, we do not see any noticeable shift in the vertical positions of the curves.

Fig. 17. Sorted halo masses for all halos in each model.

Fig. 18. Sorted halo virial radii for all halos in each model.

For the masses of the halos we cannot distinguish the curves, except for the horizontal shift at high sorting index due to differing halo numbers, as discussed earlier. This is consistent with the similarity in the halo mass functions of all models. This tells us that halo mass is not an effective means of discriminating between models, and that modification of the gravitational constant does not break the degeneracy in this property.

For the radii, however, we see a vertical shift for the positively coupled model, indicating that the halos in this model are more extended than in the uncoupled models. The halo radii of the negatively coupled model, when using the appropriate effective gravitational constant, are slightly below those of the uncoupled models, indicating less radially extended halos. This difference is eliminated when working with the Newtonian gravitational coupling, indicating that gravitationally unbound particles at large radii are included in these halos. Given the overlap of the γ₀ = −20 curve (using the Newtonian value of G) and the γ₀ = 0.3 curve, we conclude that halo radius is not an effective means of discriminating between models if the gravitational constant is assumed Newtonian.

Finally, the halo concentrations are given in Fig. 19. In this case the uncoupled models are coincident, but all coupled models are discriminated, with a clear vertical separation in their curves. The negative coupling leads to more concentrated halos, while the positive coupling leads to less concentrated halos. The separation is even more stark when working with the Newtonian gravitational constant, leading to significantly higher values for the concentration parameter in the negatively coupled model.

Fig. 19. Sorted halo concentrations for all halos in each model.

4.5.2. Derivative of the density profile

In Fig. 20 we show the derivative of the averaged density contrast profile for the first 500 halos in each simulation, with the radii normalised by the virial radius of each halo, to ensure these values range from 0 to 1. We now explain the details of how we build this plot. Rather than work with the density itself, we consider the density contrast (i.e. the density in units of the background density value). Furthermore, we are interested in the slope of the density profile, and thus we calculate the derivative of the profile in a simple way: we calculate the difference in the density contrast values between each radial bin divided by the differences in the normalised radii between each bin. These derivatives are then averaged over all 500 halos. We note that the radii are plotted on a logarithmic scale in Fig. 20.

Fig. 20. Averaged density contrast derivative for the first 500 halos in each model.

We find that there is a difference in the inner slope, at around 1%−10% of the halo virial radius, depending on the presence or otherwise of coupling, and the sign and magnitude of that coupling. Firstly, we note that the lowest values are found for the positively coupled model (regardless of the gravitational constant used) indicating, on average, less steep inner density profiles for this model. The results for the other models are broadly similar, showing steeper inner profiles in the negatively coupled models, mostly consistent with those of the ΛCDM reference simulation. It is also very clear that the inner profile is drastically steepened if the Newtonian gravitational constant is used in the halo analysis, suggesting that the halo finding process in this case includes a substantial population of additional bound particles in the inner regions of these halos.

4.5.3. Velocity dispersion profile

We now stack the velocity dispersion profiles for the first 500 halos in each model, in the same manner as we did in the previous section for the density profile derivatives. The velocity dispersions are normalised by σ(R_halo), leading to a profile value of unity at R_halo. The results of this analysis are given in Fig. 21, where we see that the Newtonian velocity dispersions are either substantially larger (smaller) at all radii for negative (positive) coupling. A less pronounced radius-independent enhancement is seen for the negatively coupled model using the effective gravitational constant; the other models are mostly indistinguishable. It is worth noting, however, that the peak of the normalised velocity dispersion profiles shifts towards smaller radii, relative to R_halo, as the coupling changes from positive to negative, with this being most pronounced again for the Newtonian analysis.

Fig. 21. Averaged velocity dispersion profiles for the first 500 halos in each model.

4.5.4. Virial ratios

In Figs. 22, 23 and 24 we plot the halo virial ratios (determined using all halo particles) for hosts, isolated halos, and subhalos, respectively. These results are also binned in mass.

Fig. 22. Virial ratios for the host halos, binned in mass. The error bars show the standard deviation of the virial ratios within each bin.

Fig. 23. Virial ratios for the isolated halos, binned in mass. The error bars show the standard deviation of the virial ratios within each bin.

Fig. 24. Virial ratios for the subhalos, binned in mass. The error bars show the standard deviation of the virial ratios within each bin.

We find that the virial ratios for the hosts, Fig. 22, for all the uncoupled models, and for our reference ΛCDM model are broadly similar and indicate somewhat overvirialised halos for all masses except the least massive. A similar picture emerges for the isolated halos and the subhalos. The isolated halos are closer to being virialised, as expected, while the subhalos, especially those that are more massive, are substantially overvirialised, most likely due to tidal interactions within their hosts.

The positively coupled model, indicated with the blue line in all three plots, shows a clear positive vertical offset in the virial ratios of their halos, across all masses, especially for the hosts. The external tidal interactions are expected to be dominant for the subhalos, leading to a virialisation behaviour, which is mostly consistent with that seen for the other models.

When analysed using the effective gravitational constant arising from the coupling, the host halos in the negatively coupled model show a slight reduction in their virial ratios, with a behaviour consistent with the uncoupled models for the isolated halos and the subhalos. If we work with the Newtonian gravitational constant, however, almost all halos in the negatively coupled model are significantly undervirialised. This is especially the case for the isolated halos. The subhalos show a trend that is consistent with the other models, but with lower values across all masses of the virial ratio.

This analysis is consistent with the bimodalities seen earlier: the positive coupling leads to dynamically disturbed clusters and groups, and to isolated halos, as compared to ΛCDM, whereas the negative coupling (at least for γ₀ = −20) leads to these halos being more dynamically cold than expected in the standard model. This effect is even more pronounced if the gravitational constant assumed in the analysis is the Newtonian value.

5. Conclusions

In this study we have extended previous work to analyse the effect of a coupling between dark matter and dark energy in the form of a quintessence scalar field, by focusing on the halo properties in a statistical manner. We have specifically concentrated on the dynamical state of the halos, particularly the hosts, confirming earlier results of a bimodality in the velocity distribution of the halo particles. This may be a general signal of positively coupled models, with an associated super-virial dynamical state. Furthermore, we have extended our work to negatively coupled models, showing a signature tendency towards unimodal velocity distributions and sub-virial halos. In addition, kinematical halo properties such as the steepness of the inner density profile and the halo concentration may help discriminate between negatively and positively coupled models. The halo mass-concentration relation has previously been shown to be a useful probe of the interaction strength in coupled models (Zhao et al. 2023). We now list our specific conclusions.

• There is a very significant reduction (40%) in structure at smaller scales for positive coupling that is sufficiently close to the value γ₀ = 1/2, as seen in our previous work. For negative coupling there is a significant enhancement (up to 30% for the larger coupling) in structure at smaller scales. On linear scales this behaviour is inverted. This is consistent with other studies of momentum-coupled (or dark scattering) models (Simpson 2010; Baldi & Simpson 2015, 2017).

• The dynamical state of the dark matter halos, particularly host halos, is more super-virial than seen in the uncoupled models and is more likely to exhibit bimodality in their particle velocities in the case of positive coupling. For (large) negative coupling the opposite is the case: host halos are somewhat less super-virial than seen in the uncoupled models and much less likely to exhibit bimodality in their particle velocities.

• The results regarding the halo dynamical states for negative coupling is even more pronounced if the analysis assumes that the dark matter dynamics are Newtonian. Then the host and isolated halos are significantly sub-virial, with very clear unimodality in the velocity distributions of most of the host halos.

• In the case of positive coupling the halo concentrations are reduced compared to uncoupled models and ΛCDM. For negative coupling the halo concentrations are increased compared to uncoupled models and ΛCDM, especially if the analysis assumes Newtonian gravitational accelerations.

• Consistent with the effect on the halo concentrations, we find that the slope of the inner density profiles is reduced (compared to the uncoupled models and ΛCDM) in the positive coupling case, whereas it is increased in the case of negative coupling, especially so in a Newtonian analysis.

The central conclusion of our work is that using averaged information from the full halo population, both kinematical information about the halos (halo concentration and slope of the inner density profile) and information about their dynamical states (bimodality and virialisation), allows us to discriminate the sign and magnitude of the coupling in these models. The discriminatory effectiveness of these properties is even stronger if the analysis is undertaken assuming Newtonian gravitational accelerations, as would be the case in an observational analysis, where a standard value of G would be taken as a minimal prior.

There are several possible extensions to our study that would be worthwhile to pursue. First, it would be beneficial to run hydrodynamical simulations in order to connect these results with the baryonic component, and thus open the door to a direct comparison with observations. In this vein, one possible line of investigation would be to consider very high resolution zoom simulations of individual galaxies or galaxy groups embedded in a cosmological context in the presence of these couplings. In addition, as always, it would be helpful to have higher resolution simulations to be able to explore these effects in lower mass halos. From a theoretical point of view, it would be interesting to further explore the parameter space of these models by considering larger negative coupling values, as well as other potentials for the scalar field, giving rise to alternative background evolutions. In this context, it would be of interest to explore models where the coupling has relevance at higher redshift, perhaps tending to the standard behaviour at late times. Generalisations of these models are also possible, such as considering other forms for the coupling term, or even mixing a momentum coupling with a density coupling.

The signatures of these models discussed in this study may eventually present an opportunity to test these ideas with observational data, allowing us to further constrain the vast space of possible cosmological models.

Acknowledgments

GNC wishes to thank Alkistis Pourtsidou for very useful discussions.

References

All Tables

All Figures

	Fig. 1. Hubble parameter normalised by that of ΛCDM for the three values of the coupling constant.
In the text

	Fig. 2. Modified Euler equation coefficients for the three values of the coupling constant.
In the text

	Fig. 3. Density projection plots of all simulations considered in this work. The leftmost panel is our ΛCDM run, and the other three panels (from left to right) are coupled models with: γ0 = 0.3, γ0 = −0.3, γ0 = −20.
In the text

	Fig. 4. Power spectra normalised by the corresponding uncoupled model.
In the text

	Fig. 5. Halo mass functions for all simulations considered in this study. The comparison line given is from Watson et al. (2013), as defined in the Colossus Python package.
In the text

	Fig. 6. Velocity dispersion as a function of mass for all halos in all the coupled models considered in this study and the reference ΛCDM model.
In the text

	Fig. 7. Histogram of bimodality parameter values calculated from the bi-Gaussian distribution fitted to the velocity magnitude histograms of all particle velocities in the 500 most massive host halos of each simulation.
In the text

	Fig. 8. As in Fig. 7, but for isolated halos.
In the text

	Fig. 9. As in Fig. 7, but for subhalos.
In the text

	Fig. 10. As in Fig. 7, but using the Newtonian gravitational constant in the halo finding process.
In the text

	Fig. 11. As in Fig. 8, but using the Newtonian gravitational constant in the halo finding process.
In the text

	Fig. 12. As in Fig. 9, but using the Newtonian gravitational constant in the halo finding process.
In the text

	Fig. 13. Probabilities of unimodality in 3D velocity distributions for all models.
In the text

	Fig. 14. Probabilities of bimodality in 3D velocity distributions for all models.
In the text

	Fig. 15. Probabilities of unimodality in 1D velocity distributions for all models.
In the text

	Fig. 16. Probabilities of bimodality in 1D velocity distributions for all models.
In the text

	Fig. 17. Sorted halo masses for all halos in each model.
In the text

	Fig. 18. Sorted halo virial radii for all halos in each model.
In the text

	Fig. 19. Sorted halo concentrations for all halos in each model.
In the text

	Fig. 20. Averaged density contrast derivative for the first 500 halos in each model.
In the text

	Fig. 21. Averaged velocity dispersion profiles for the first 500 halos in each model.
In the text

	Fig. 22. Virial ratios for the host halos, binned in mass. The error bars show the standard deviation of the virial ratios within each bin.
In the text

	Fig. 23. Virial ratios for the isolated halos, binned in mass. The error bars show the standard deviation of the virial ratios within each bin.
In the text

	Fig. 24. Virial ratios for the subhalos, binned in mass. The error bars show the standard deviation of the virial ratios within each bin.
In the text