© The Authors 2024

1. Introduction

The asevolution code (Christiansen et al. 2023) extends the relativistic N-body particle-mesh code gevolution (Adamek et al. 2016; Daverio et al. 2016) by adding the (a)symmetron (Perivolaropoulos & Skara 2022b; Hinterbichler & Khoury 2010), which is a dynamical dark energy component. The theory is constructed so that the scalar field remains dormant during the early universe and undergoes a late-time phase transition that partitions the universe into domains of either the positive or negative field value minima. As the field takes on non-zero values, it starts mediating a fifth force that enhances the clustering of matter. The phenomenology of the (a)symmetron has been thought about in relation to cosmological tensions such as the Hubble tension (Perivolaropoulos & Skara 2022b,a; Hägås & Mörtsell 2023) and cosmic dipole tension (Perivolaropoulos 2023), and it is plausibly connected to many more (Peebles 2022, 2023; Naik & Burrage 2022; Afzal et al. 2023). The model has also been thought about as a candidate for dark matter (Burrage et al. 2017, 2019); because of the environmental dependence of the mass of the field, the dark matter mass would correlate with the halo and galaxy masses, which has recently been shown to alleviate the core-radius scaling tension (van Dissel et al. 2024). In the companion paper (Christiansen et al. 2024), we simulate the production of stochastic gravitational waves by the (a)symmetron scalar field and indicate a region in the parameter space that may account for the observed spectral index of the stochastic gravitational wave signal as observed by the NANOGrav collaboration (Afzal et al. 2023). The symmetron screening mechanism, similarly to chameleon or Vainshtein screening, makes the fifth force inactive in local environments where the model would otherwise be subject to strict constraints from Solar System experiments (Bertotti et al. 2003; Esposito-Farese 2004; Tsujikawa et al. 2008). Requiring the local Solar System or galaxy to be screened imposes some constraints on the allowed parameter space (Hinterbichler & Khoury 2010) that have recently been suggested, which were overestimated in the past due to the non-trivial interaction of the field with the environment (Burrage et al. 2024). We make asevolution publicly available to download¹.

While we previously found convergence in results such as the background evolution and power spectra in Christiansen et al. (2023), we did not explore the convergence and resolution criteria for the domain wall network thoroughly. For the parameter choice and simulation settings we investigated, the agreement in the power spectrum between the dynamic and quasi-static solvers was indeed good, indicating that time-convergence of dynamics must be established using other quantities. Furthermore, the sizes of the simulations that formed stable domain walls were large, so that a systematic exploration of the convergence of the dynamic and non-linear part of the model would have been too expensive. Here, we set out to more rigorously establish the physical validity and convergence of the solutions found by our code. We do this in Sect. 4 by considering temporal and spatial resolution separately, and by constructing ideal situations that we take as representative of subsystems of the cosmological setup. We then move on to low-resolution cosmological simulations to find results on the convergence that agree with those of the ideal setups. Finally, in Sect. 5, we extrapolate our results to a suite of five high-resolution cosmological simulations. We visualise the behaviour of the scalar field in animations available online and on our website².

2. Conventions

We set the speed of light c equal to 1 and use the metric signature ( − + +   + ). The conformal factor A has an implicit dependence on the scalar field ϕ. All symbols with a tilde over are defined in the Jordan frame, such as the trace of matter stress-energy tensor in the Jordan frame ${\stackrel{\sim }{T}}_m$; in the Einstein frame, it is written T_m. By dot ẋ, we mean that the derivative is made with respect to cosmic time; primed quantities (x′) have their derivatives taken with respect to conformal time. We refer to the Compton wavelength corresponding to the mass μ as L_C, while the general Compton wavelength valid for all T_m and spatial wavenumbers k is referred to as L. We refer to a_* as the scale factor at spontaneous symmetry breaking in a homogeneous universe, and a_SSB as the scale factor of actual symmetry breaking. Finally, we refer to the ‘true vacuum’ as the deepest minimum of the (a)symmetron potential when T_m = 0, while we just write ‘vacuum’ for the minimum otherwise, or ‘false vacuum’.

3. Theory

The asymmetron (Perivolaropoulos & Skara 2022b) is a generalisation of the symmetron (Hinterbichler & Khoury 2010), both of which are defined with the action

$$\begin{array}{c}\hfill S={\displaystyle \int {\mathrm{d}}^4}x\sqrt{-g}[\frac{1}{2}{M}_{\mathrm{pl}}^{2}R-\frac{1}{2}{\mathrm{\nabla }}_{\mu }\varphi {\mathrm{\nabla }}^{\mu }\varphi -V(\varphi )]+S_m({\psi }_i,{\stackrel{\sim }{g}}_{\mu \nu }),\end{array}$$(1)

where the integral is taken over spacetime, g is the determinant of the Einstein frame metric g_μν, R is the Ricci scalar corresponding to g_μν, ∇ is the covariant derivative, V(ϕ) is the potential of the (a)symmetron scalar field ϕ, and $S_m({\psi }_i,{\stackrel{\sim }{g}}_{\mu \nu })$ is the matter action. S_m is defined with the Jordan frame metric ${\stackrel{\sim }{g}}_{\mu \nu }$, but is otherwise as we know from the standard model of cosmology, consisting of cold dark matter (CDM), the cosmological constant, radiation, neutrinos, and baryons. The relation between the Jordan and Einstein frame metrics is given by the conformal transformation

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_{\mu \nu }=A^2(\varphi )g_{\mu \nu },\end{array}$$(2)

in which the conformal factor A(ϕ) is defined

$$\begin{array}{c}\hfill A(\varphi )\equiv 1+\frac{1}{2}{\left(\frac{\varphi }{M}\right)}^2\equiv 1+\mathrm{\Delta }A,\end{array}$$(3)

where M is the conformal coupling. Equation (3) may be viewed as the attenuated Taylor expansion of some more generic function that is symmetric in ϕ as long as ΔA is small, which it is for the parameters that we consider in this paper. We can define the effective potential V_eff that enters into the equations of motion of the (a)symmetron,

$$\begin{array}{c}\hfill V_{\mathrm{eff}}\equiv V-lnA(\varphi )T_m=V_0-\frac{1}{2}{\mu }^2{\varphi }^2-\frac{1}{3}\kappa {\varphi }^3+\frac{1}{4}\lambda {\varphi }^4-lnA(\varphi )T_m,\end{array}$$(4)

which has the additional term lnA(ϕ)T_m that comes from the conformal transformation of the metric in the matter action (see e.g. Christiansen et al. 2023). V₀ is a constant, μ is the mass, λ is the coupling constant, and T_m is the trace of the Einstein frame stress-energy tensor. κ is the coupling constant in the cubic term, and we recover the symmetron for the choice κ = 0. We can define a set of phenomenological parameters that map to the Lagrangian parameters as in Brax et al. (2012), Davis et al. (2012), which was extended to the case of the asymmetron in Christiansen et al. (2023)

$$\begin{array}{cc}\hfill L_C& =\frac{1}{\sqrt{2}\mu }\equiv \frac{{\xi }_{\ast }}{H_0}\approx {\xi }_{\ast }·2998{\mathrm{Mpc}\mathrm{h}}^{-1},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill M& ={\xi }_{\ast }\sqrt{6{\mathrm{\Omega }}_{m,0}/a_{\ast }^3},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill \lambda & ={\left(\frac{M_{\mathrm{pl}}}{M^2}\right)}^2\frac{{\mu }^2}{{\overline{\beta }}^2-{(\mathrm{\Delta }\beta /2)}^2},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \kappa & =\frac{\lambda M^2\mathrm{\Delta }\beta }{M_{\mathrm{pl}}},\hfill \end{array}$$(8)

where L_C is the Compton wavelength corresponding to the mass μ of the field, H₀ is the present time Hubble parameter, Ω_m is the present time background energy density parameter for matter, $M_{\mathrm{pl}}=1/\sqrt{8\pi G_N}$ is the Planck mass, and G_N is Newton’s gravitational constant. The Lagrangian parameters are then defined in terms of (ξ_*, z_*, β₊, β₋), where ξ_* = L_CH₀ is the ratio of the Compton wavelength to the present time Hubble horizon, z_* is the cosmological redshift at symmetry breaking for a homogeneous universe, and β_± are the force strengths relative to the Newtonian force in the positive and negative symmetry-broken minima. We also defined Δβ = β₊ − β₋ and $\overline{\beta }=({\beta }_{+}+{\beta }_{-})/2$. The field equations can be found using the Euler-Lagrange equation, and they can be set in the form (Llinares & Mota 2013) when neglecting relativistic corrections,

$$\begin{array}{cc}& \dot{q}=a{\mathrm{\nabla }}^2\chi -a^3{\mu }^2\{{\chi }^3-\overline{\kappa }{\chi }^2+(\eta -1)\chi \},\hfill \end{array}$$(9)

$$\begin{array}{cc}& q=a^3\dot{\chi },\chi =\varphi /v,\hfill \end{array}$$(10)

where $v=\mu /\sqrt{\lambda }$ is the true vacuum expectation value of the scalar field in the symmetron scenario. The extension of the equations of motion to include relativistic corrections is provided in Christiansen et al. (2023), but is not considered for the purposes of the present consideration.

4. Convergence

In this section, we review our results on the convergence of the field solver in finding the domain wall network. In Christiansen et al. (2023), we made use of the leapfrog solver to evolve the scalar field. During convergence testing, we initially found a high degree of sensitivity of the domain wall network to initial conditions as well as to changes in both temporal and spatial resolution. Therefore, we found it appropriate to first investigate a small selection of solvers to determine whether a different field solver might perform better. This made us choose the explicit fourth-order Runge-Kutta solver in the following. We refer to Appendix C for a review of the different solvers and their convergence. A useful property of the explicit Runge-Kutta solver that we make use of here is that it becomes unstable and diverges after accumulating a sufficient amount of error. Implicit schemes have larger stability regions and can provide (sometimes unconditionally) stable and reasonable solutions even when the physical system is not being resolved (Linge et al. 2017). We use the stability of the explicit scheme as an indicator for whether the physical solution is being resolved, which is motivated by our ideal system considerations in Sect. 4.1. We first review the relevant scales of the model. Next, we describe the convergence of physical solutions during the time after symmetry breaking. We use both idealised situations and realistic cosmological cases, where we also consider the initial conditions with care. Unless otherwise specified, we use in all instances below the symmetron parameters (L_C, z_*, β_*) = (1 Mpc h⁻¹, 2, 1) that were also studied in Llinares & Mota (2014), Christiansen et al. (2023). In this section, we refer to the Courant factor as the factor defined

$$\begin{array}{c}\hfill C_f=v_{\mathrm{vel}}\mathrm{d}\tau /\mathrm{d}x\equiv vC,\end{array}$$(11)

where for the scalar field C_f, ϕ = C_ϕ, we set $v_{\mathrm{vel}}^{(\varphi )}=1$, and dτ → dτ_ϕ; τ is the conformal time coordinate. On the other hand, for CDM, C_f, cdm, we use v_vel = v_cdm, max ∼ 0.02 and dτ → dτ_cdm. In other words, we vary the time stepping of CDM and the scalar field separately, requiring dτ_ϕ ≤ dτ_cdm. For convenience, we show the time convergence for the CDM equations in terms of C_cdm = dτ/dx.

4.1. Relevant scales

In order to reliably evolve the scalar field, the frequency scale of its oscillations should be smaller than the Nyquist wavenumber of our solver. In Christiansen et al. (2023), we found the linear solution of the perturbed symmetron in Fourier space, which is a damped oscillation with phase frequency ω, so that

$$\begin{array}{c}\hfill {\omega }^2=k^2+a^2{m}^2-{\mathcal{H}}^2\equiv \frac{1}{2L^2},\end{array}$$(12)

where the mass m is defined below, ℋ = aH is the comoving Hubble parameter, k is the comoving wave number, and we defined the general Compton scale L. For large scales kL ≪ 1, for a subhorizon general Compton wavelength LH ≪ 1, the relevant conformal timescale is

$$\begin{array}{c}\hfill L=\frac{1}{\sqrt{2}am}.\end{array}$$(13)

The mass m above makes use of the definition $m^2={\left|{\partial }_{\varphi }^{2}V\right|}_{\varphi =v}$, which gives

$$\begin{array}{cc}& \mathrm{before}\mathrm{SSB},\left|T_m\right|\ge {\rho }_{\ast }:m=\mu \sqrt{-\frac{T_m}{{\rho }_{\ast }}-1},\hfill \end{array}$$(14)

$$\begin{array}{cc}& \mathrm{after}\mathrm{SSB},\left|T_m\right|<{\rho }_{\ast }:m=\sqrt{2}\mu \sqrt{1+T_m/{\rho }_{\ast }},\hfill \end{array}$$(15)

where ρ_* ≡ μ²M² is the energy scale of the symmetry breaking. On the background, when Ω_ϕ ≪ Ω_m and ΔA ≪ 1, T_m/ρ_* ≃ −(a_*/a)³. For the fiducial symmetron model choice (L_C, z_*, β_*) = (1 Mpc h⁻¹, 2, 1), we see the background evolution of the conformal Compton scale shown in Fig. 1.

Fig. 1. General Compton wavelength of the symmetron field at the background as a function of the scale factor. The parameter choice is (LC, a*, β*) = (1 Mpc h−1, 0.33, 1), and the horizontal lines indicate the initial and final Compton scales.

We demonstrate more clearly in Sects. 4.3 and 4.4 that this scale relates to whether we are resolving the system. Taking this for granted for now, we show in Fig. 1 that the timescale of the oscillations is initially set very small due to the high mass gained from the background energy density. It is therefore very expensive to evolve the field reliably from z ∼ 100 until symmetry breaking, for instance. The mass m goes to zero around the time of the symmetry breaking. From expression (12), the Hubble horizon is then $k^2-{\mathcal{H}}^2\sim \frac{1}{2L^2}\to 0$, which is often thought to set the scale of percolation along with causality, as discussed for example in Llinares & Pogosian (2014).

4.2. Courant condition

The Sect. 4.1 indicates the resolution scale we expect to need in order to resolve the symmetron on the background. When the field is resolved, we still expect a condition on the Courant factor, which comes from the stability of the numerical scheme. We determine this condition next. We constructed an idealised system where the matter field is homogeneous and evolves as ρ_m = ρ_m, 0/a³. We initialised the scalar field with its expectation value for a homogeneous background v and with a sign difference ±v on the different halves of the box, so that there was a centre domain wall. At first, we ran the Gaussian relaxation iterations presented in Llinares et al. (2014), Christiansen et al. (2023) to thermalise the field configuration and smooth the steep gradient at the centre. The higher the tolerance we used for the relaxation, the closer the configuration of the field became to the analytic profile. If it had exactly the analytic profile, the field was unexcited and stayed dormant. If the box was too small for the analytic field profile to establish, the relaxation instead placed the field entirely at one minimum, decided by the numerical error. Since we are interested in considering a dynamic scenario, we set the relaxation tolerance to not relax the field entirely. In the right panel of Fig. 2, we show the initial profile of the field. Overlaid is the analytic profile corresponding to a smaller Compton wavelength of 0.5 Mpc h⁻¹, although the field we used here had a Compton wavelength of L_C = 1 Mpc h⁻¹. The field was then in an excited state, where the profile was initially squeezed and started to oscillate when it was evolved.

Fig. 2. Initial setup for the bare wall spatial and temporal convergence experiments. The right panel shows the overlaid analytic profile corresponding to a Compton wavelength of 0.5 Mpc h−1, although the model here is defined with a Compton wavelength of 1 Mpc h−1.

In Fig. 3 we show the evolution of a fraction of lattice points in the positive minimum for different choices of temporal and spatial resolutions, in addition to different Compton and initial scales. We take the convergence of this fraction to be a proxy for whether the field configuration is properly resolved. To gain more control of the initial conditions as we compare different spatial resolutions, we simply initialised the field in the analytical profile corresponding to a smaller Compton wavelength of L_IC = 0.5 Mpc h⁻¹. Because we used the explicit Runge-Kutta solver, the solutions diverged immediately when they accumulated a sufficient amount of error, and this indicated the critical Courant factor for stability of the numerical scheme, which we identify here to be C_ϕ, crit ∼ 1.4. This means that we require dτ ≤ 1.4 dx. As long as the Compton wavelength is resolved, the critical Courant seems to be independent of dx,  L_IC, and L_C.

Fig. 3. Volume fraction of the field that is in the positive minimum as a function of redshift for the oscillating bare wall system for variations of the Courant factor Cϕ. The four plots vary in spatial resolution dx, initial wall scale LIC, and Compton wavelength LC.

4.3. Spatial resolution

The smallest spatial scale introduced by the symmetron is the interface between two different domains. The profile of the field along an axis, x, orthogonal to the domain wall looks like what is shown in the right panel of Fig. 2. As was shown in Llinares & Pogosian (2014), for example, the analytic expression for a domain wall in a homogeneous matter distribution is given by

$$\begin{array}{c}\hfill \chi (x,a)=\sqrt{1-{\left(\frac{a_{\ast }}{a}\right)}^3}tanh\left(\frac{\mathit{ax}}{2L_{\mathrm{C}}}\sqrt{1-{\left(\frac{a_{\ast }}{a}\right)}^3}\right),\end{array}$$(16)

where L_C is the Compton wavelength corresponding to the Lagrangian mass parameter μ, as shown in Eq. (5). In Fig. 4, we show the oscillation in a fraction of lattice points in the positive minimum for different spatial resolutions, keeping the initial wall and Courant scales L_IC, L_C = 0.5, 1 constant. In each case, we fixed the Courant factor C_ϕ = 1 < 1.4 according to the numerical stability found in Sect. 4.2. For a spatial resolution of dx = 2 Mpc h⁻¹, which is higher than the general Compton resolution, we saw an aliasing of the correct oscillation. In Fig. 4 we instead used a Courant factor C_ϕ = 0.5 for the case when dx = 2 Mpc h⁻¹ to determine whether the spatial resolution is important by itself. In this case, some of the fluctuations were captured, but the second and third beats were missed while the fourth was slightly displaced. Our conclusion then is that the spatial resolution itself has to be smaller than L, which in turn is $L\ge L_{\mathrm{C}}/\sqrt{2}$ in the symmetry-broken phase. For still lower spatial resolutions, the gradients are increasingly well resolved, and the oscillation is smoother.

Fig. 4. Volume fraction of the field in the positive minimum as a function of redshift for the oscillating bare wall system for a selection of spatial resolutions. The coarsest simulation has the Courant factor choice Cϕ = 0.5 and would otherwise not resolve the oscillations. The rest are kept with Cϕ = 1.

4.4. Convergence in cosmological setups

In the above Sect. 4.2, we found the condition for numerical stability, and in Sects. 4.1 and 4.3, we indicated the relevant scale to resolve the physical system. We applied these criteria to a cosmological scenario. Some complications to the previous ideal situation arose: first, the full set of partial differential equations that we considered was extended by the dark matter and metric perturbations, which are no longer treated at the background. In principle, this may affect the numerical stability. Next, since the matter fields are not treated homogeneously, both regions of the field are in the symmetry-broken and -unbroken phases. In the unbroken symmetry, for large ρ_m, $L\sim M/\sqrt{{\rho }_m}$ may be smaller than L_C, and it may be important to resolve if we want convergence on the exact field configuration of the scalar is to be achieved. In the following, we ensured to resolve the Compton scale of the symmetry-broken phase by fixing dx = 0.8 Mpc h⁻¹ and C_ϕ ≲ 1.

We initialised the scalar field from a scale-invariant power spectrum at z = 100 for a demonstration: When the field is evolved at a time resolution that is too coarse, for example 0.7 Mpc h⁻¹, when the general Compton at the background at initialisation time is 0.53 Mpc h⁻¹, the realisations of the domain wall network are completely different in general when the time-stepping is varied slightly. The evolution in the pre-symmetry-broken phase is not very interesting for our parameter choices, however, and the power spectrum of the field remains close to scale-invariant right up until the symmetry breaking (see Appendix B). We demonstrate in Fig. 5 that we obtain a qualitatively similar domain decompositioning when we instead initialise the scale-invariant power spectrum at some short duration before the symmetry breaking, where we used a very coarse spatial resolution and the more stable leapfrog solver for demonstration. A late initialisation amounts to a large optimisation of the computational cost.

Fig. 5. Comparison of snapshots of the scalar field χ at redshifts z = 1.9 and 1.8 for a late initialisation at z = 2.5 and for an early initialisation at z = 100. The amplitude of the late-time initial conditions is chosen to be similar to the evolved early-time conditions. In the bottom left snapshot, a black square indicates the relative size of a (500 Mpc h−1)3 volume box.

By evolving the scalar field in a cosmological scenario together with the CDM and metric perturbations, using the resolution criteria developed in the ideal-case bare wall scenario, we find solutions that blow up and diverge. This indicates that not all relevant scales were resolved; in the current setup, we used a spatial resolution of dx = 0.8 Mpc h⁻¹ and a Courant factor of the scalar C_ϕ ∼ 1 < 1.4. We expected the stability to change both because new relevant scales were introduced by the density contrasts that maintained the ℤ₂ symmetry in high-density areas and because the numerical scheme changed from co-evolving CDM, metric perturbations, and the scalar field. Because we varied dτ_ϕ separately, we expected the numerical scheme to be similar to the ideal case when the other dynamical degrees of freedom evolved slowly in comparison, which we expected them to. The new minimum scale in the symmetric phase owing to the density contrasts is

$$\begin{array}{c}\hfill L^{\min }\sim \frac{L_{\mathrm{C}}}{\sqrt{a_{\ast }^3}}\sqrt{\frac{a}{{\delta }_m^{\max }}}.\end{array}$$(17)

For the spatial resolution that we considered in the high-resolution simulations, dx ∼ 0.4 Mpc h⁻¹, ${\delta }_m^{\text{max}}~3000$, meaning that conservatively, L^min ∼ 0.013 Mpc h⁻¹. In other words, we did not expect to resolve the false vacuum fluctuations of the field in the highest-density environments for computationally reasonable resolutions, and we therefore expect some small degree of phase randomisation. We were still able to avoid divergence of the solutions by adjusting the Courant factor to C_ϕ ∼ 0.15. In this event, we resolved the Compton scale in the symmetric phase for density contrasts δ_m ≲ 140 in the event of dx ∼ 0.4 Mpc h⁻¹, which was satisfied by almost the entire simulation box volume. For the choice dx ∼ 0.8 Mpc h⁻¹, it is instead sufficient to use C_ϕ ∼ 0.2. As in the ideal case, in Sect. 4.2, we expect the non-divergence of the solutions to indicate that the system is sufficiently resolved.

In the case of ΛCDM, the timescale of the evolution for the CDM is much longer than what we considered here for the symmetron. For large simulations, Courant factors can typically be used in gevolution so that C_cdm = dt/dx ∼ 50. Since the CDM particles are slow, for the type of resolution we considered, dx ∼ 0.4 Mpc h⁻¹, the maximum velocity is $v_{\mathrm{cdm}}^{\max }\sim 0.02$, which means C_f, cdm ∼ 1 in this instance. However, due to the interaction of the symmetron and the dark matter, additional dynamics of the CDM is introduced on the symmetron timescale. As an intuitive requirement for resolving this dynamic, we suggest $v_{\mathrm{cdm}}^{\max }\mathrm{d}{\tau }_{\mathrm{cdm}}\ll L(T=0)\le L_{\mathrm{C}}/\sqrt{2}$. The fifth force is screened when the symmetry is restored, so that the CDM time-stepping only needs to resolve the symmetry-broken Compton scale. For dx ∼ 0.4 Mpc h⁻¹, we then used a Courant factor of C_cdm = 1, while for dx ∼ 0.8 Mpc h⁻¹, we chose C_cdm = 0.4. In both instances, the maximum displacement of a dark matter particle during one iteration is $\mathrm{\Delta }{x}_{\mathrm{cdm}}^{\max }\sim 0.4·v_{\mathrm{cdm}}^{\max }{\mathrm{Mpc}\mathrm{h}}^{-1}\sim 0.01$ Mpc h⁻¹ ∼ 1% of the scale of the domain walls. In Fig. 6 we show convergence in the fraction of lattice points in the positive symmetry-broken minimum as a function of redshift for decreasing Courant factors and for a spatial resolution of dx ∼ 0.8 Mpc h⁻¹. There is approximate convergence for C_cdm ∼ 1.

Fig. 6. Fraction of nodes in positive minimum for simulations run with initial conditions found from a ΛCDM simulation with a Courant factor Ccdm = 1. The box size is B = 410 Mpc h−1, and to the right, we show the fraction of domains within a B = 205 Mpc h−1 subvolume. We used 5123 grids. We chose Cϕ = 0.15.

5. High-resolution simulations

In Sect. 4, we established resolution criteria, guided by physical intuition, ideal system simulations, and cosmological system convergence testing. In this section, we apply them to a suite of five 500 Mpc h⁻¹, 1280³ grid simulations in order to predict the non-linear and dynamically dependent observables that are not accessible by use of the non-converged, semi-analytic, or quasi-static approaches used in the past.

In Table 1 the (a)symmetron parameters are shown for each of the simulations I–V. The cosmological parameters are otherwise as found by the Planck satellite (Planck Collaboration VI 2020), with the Hubble parameter h = 0.674, and the cosmological energy density parameters at z = 0 are Ω_b, Ω_m, Ω_Λ = 0.049, 0.264, 0.687. We adjusted Ω_Λ to satisfy the flatness criterion ∑_iΩ_i = 1, although the symmetron parameter choices considered have Ω_ϕ ≲ 10⁻⁵, so that the adjustment is very small. We used the cosmic microwave background temperature T_cmb = 2.7255 K, the spectral amplitude and index A_s, n_s = 2.1 × 10⁻⁹, 0.965, and the spectral tilt k_p = 0.05 Mpc⁻¹. For simplicity, we did not consider massive neutrinos and set the effective number of relativistic degrees of freedom N_eff = 3.046. We chose the number of CDM particles equal to the number of grids.

5.1. Description of the data

In Fig. 7 we give an overview of the symmetry-breaking process for the different parameter choices, illustrated by the fraction of the simulation volume in the positive field-value minimum and the average of the square of the field as a function of redshift. We comment on some interesting features of the background:

1. The domain wall for the early symmetry-breaking scenario, model I, is only quasi-stable and collapses at redshift z ∼ 0.5.

2. The asymmetron scenario, model V, is only quasi-stable and undergoes a slow collapse from redshift z ∼ 0.25.

3. All three instances of symmetron Δβ = 0 late-time symmetry breaking z_* = 0.1, models II–IV, seem to be stable.

4. While the timescale from the symmetry breaking of the first voids to the last voids, that is, the width of the transition, is rather small for these parameter choices, it grows slightly longer for smaller L_C.

5. The timescale for the minima to reach the true vacuum, indicated by ⟨χ²⟩, is comparatively long, so that the vacua still evolve throughout the simulation time.

6. The bump in ⟨χ²⟩ at the onset of the symmetry breaking for model I bears witness of an initially non-adiabatic evolution where the field falls into the minimum. The phase transition seems to occur more smoothly in the late-time symmetry breaking, models II–V.

Fig. 7. Overview of background quantities for the suite of simulations. Top: fraction of the simulation volume where the scalar field occupies the positive minimum for the five different models I–V. Bottom: average over the simulation volume of the square of the scalar field χ2. Both quantities are plotted as a function of the redshift z.

We discuss the observational implications in Sect. 6. In Figs. 8 and 9 we give an overview of the evolution of the exact field configuration as viewed through snapshot slices through the simulation volume at four different redshifts.

1. In addition to the earlier symmetry breaking in the smaller Compton wavelength model, IV, a different domain partitioning occurs in this instance as a result of the different field phases at the time of symmetry breaking. We expect this both from the different oscillation timescale of the field ∼L and from the different symmetry-breaking time.

2. Variation in the coupling β shows no strong effects, but some small collapses in the weaker coupling case, model II, indicate that the domain walls are less strongly pinned to filaments in this instance.

3. Some of the domain wall collapses are especially clear in Fig. 9, where the velocity shows large amplitude oscillations around the site of the collapse. For example, the dark regions in the bottom right and top left corner of the z = 0.2 plot in the fourth column (labelled V) of Fig. 9 correspond to regions where the domain wall has moved between z = 0.4 and z = 0.2 in Fig. 8.

4. Otherwise, we note that the velocity field correlates on the scale of the Compton wavelength ∼L, being on smaller scales in the instance of model IV.

Fig. 8. Snapshots of the scalar field χ. Shown for models II–V from left to right, for four different redshifts from top to bottom.

Fig. 9. Snapshots of the velocity of the scalar field $q=a^3\dot{\chi }$. Shown for models I–V from left to right, for four different redshifts from top to bottom. The field is normalised to the conformal Hubble parameter.

Finally, in Figs. 10 and 11, we show the scalar configuration at an instance for models I and IV, respectively. We show the same simulation slice in four different fields: The scalar field χ; the isotropic gradient $\left|\mathbf{\nabla }\chi \right|/\sqrt{3}$, normalised to the bare wall analytic expectation at the centre, ${\left|\mathbf{\nabla }\chi \right|}_{\mathrm{wall}}=\frac{\sqrt{3}a}{2L_{\mathrm{C}}}(1-{(\frac{a_{\ast }}{a})}^3)$, in the case of model I and to the Hubble horizon ℋ in the case of model II–IV, where the analytic expectation is not available since a < a_*; the equation of state of the field, ω_ϕ in the rest frame of the simulation; and finally, the energy density parameter of the field Ω_ϕ = ρ_ϕ/ρ_c, where ${\rho }_c=3H^2{M}_{\text{pl}}^2$ is the critical density. Our comments are listed below.

1. The energy scale of the domain walls is set by $V\sim \frac{\lambda }{4}{v}^4\sim \frac{9{\beta }^2{\xi }^2{\mathrm{\Omega }}_m^2}{2a_{\ast }^6}{H}_0^2{v}_{\chi }^4(a)$, and the energy scale is ∼100 times smaller in model IV than in model I, whereas the above equation predicts ∼5%. We defined $v_{\chi }=\sqrt{1+T_m/{\rho }_{\ast }}$. However, it is complicated because v should be evaluated locally in models II–IV.

2. The amplitude of the isotropic gradient is ∼1 for both normalisations, indicating that the analytic homogeneous background solution correctly sets the scale in model I, whereas the Hubble scale is a good indication in model IV. The animations show that the gradients gradually grow larger with time in models II–V, while they are the largest for the domain formation and collapse in model I.

3. The equation-of-state parameter ω_ϕ correlates with the domain walls and with the structure, where it is ω_ϕ ∼ 0. On the other hand, in the underdense parts, it goes towards ω_ϕ ∼ −1. This agrees with the observation that was made in Christiansen et al. (2023).

Fig. 10. Snapshots of slices of the simulation volume. From top left to bottom right: the scalar field χ, the isotropic $\mathrm{\nabla }\chi \equiv \left|\mathrm{\nabla }\chi \right|/\sqrt{3}$ normalised to the one expected for the centre of the analytic bare wall, the equation-of-state parameter of the field ωϕ, and the energy density parameter Ωϕ = ρϕ/ρc. The snapshots are from model I and are shown at redshift z = 0.9 at a critical moment right before the two blue domains join.

Fig. 11. Snapshots of the simulation volume. From top left to bottom right: the scalar field χ, the isotropic gradient ∇χ≡$\left|\mathrm{\nabla }\chi \right|/\sqrt{3}$ normalised to the conformal Hubble parameter ℋ, the equation-of-state parameter of the field ωϕ, and the energy density parameter Ωϕ = ρϕ/ρc. The snapshots are from model IV.

We discuss some specific measurements/quantities in more detail below. We start with the matter power spectra.

5.2. Power spectra

In Figs. 12 and 13, we show the matter power spectrum enhancement relative to the ΛCDM of models I and II–V, respectively. Figure 12 also shows a comparison of the power spectra to those found by the linear prediction from CLASS (Lesgourgues 2011). The agreement is good in general, except at the very largest and smallest scales. This is due to cosmic variance and to the grid resolution, respectively. We expect the power spectra to be reliable to about k_max ∼ k_Nyq/7, where k_Nyq = π/dx is the Nyquist frequency. On account of the power spectra, we make the following observations:

1. Each of the relative enhancements in the spectra has a bump at a characteristic length scale; this scale seems to be related to the Compton wavelength L_C of the field, and it seems to move towards smaller scales when comparing model III and IV.

2. The matter power spectrum appears not to be particularly affected by the asymmetry of the potential in model V when comparing it to model III.

3. In models II and III, which have relative force strengths β_II = 1, β_III = 8, respectively, the shapes in the enhancement are similar, but the overall amplitude difference corresponds to ${({\overline{\beta }}_{\mathrm{II}}/{\overline{\beta }}_{\mathrm{III}})}^2$, as expected from Brax et al. (2012), for example.

4. Finally, comparing models I and II, we see the effect of varying the critical density ${\rho }_{*}~a_{*}^{-3}$. In this instance, shape and amplitude both appear to be affected. The amplitude difference of the bumps is roughly equal to ${\left(a_{*}^{\text{I}}/a_{*}^{\text{II}}\right)}^5$ and approaches ${\left(a_{*}^{\text{I}}/a_{*}^{\text{II}}\right)}^4$ for large k. Additionally, the power spectrum is less suppressed for larger modes k in model I, indicating a lower degree of screening, as expected from the higher critical energy ${\rho }_{*}\propto a_{*}^{-3}$. We expect a different scaling for large a_*, and this needs to be investigated separately.

Fig. 12. Matter power spectra for model I as a function of wavenumber k indicated for different redshifts z. Top: overlaid power spectra from the symmetron asevolution simulation plotted as a solid line, and the corresponding ΛCDM simulation is shown as dots. The CLASS power spectrum for the same ΛCDM parameters is shown as a dashed line. Bottom: relative difference with respect to ΛCDM of the power spectra from redshift z = 3 to z = 0 indicated in graded colours from blue to yellow. The gap in the spectra in the lower panel is not a physical feature and is caused by a typo in the settings file of the simulation. This caused two spectra to be skipped.

Fig. 13. Relative difference of matter power spectra with respect to the ΛCDM case, with redshift z = 1 to z = 0 indicated in graded colours from blue to yellow. The spectra are plotted as a function of wave number k. The top left to bottom right panels show models II–V.

In other words, the results agree with our physical intuition regarding the effect of the parameter choices and with the results obtained in previous approximate treatments of the symmetron in Brax et al. (2012), Llinares et al. (2014), for instance.

5.3. Halo mass functions

The halo mass functions (HMFs) for models I and II–V are displayed in Figs. 14 and 15, respectively. They indicate the abundances of dark matter haloes at different mass ranges, and are thought to roughly correspond to the galaxy abundances (Guo et al. 2010). We found the dark matter haloes using the Rockstar phase-space halo finder (Behroozi et al. 2013), and computed the Tinker10 (Tinker et al. 2008) ΛCDM expectation for the HMFs using Pylians (Villaescusa-Navarro 2018) and CLASS (Lesgourgues 2011). We comment on Figs. 14 and 15 below.

1. The top panel of Fig. 14 shows the agreement with analytic estimates from Tinker et al. (2008) constructed out of CLASS power spectra for the ΛCDM model. The agreement is better for later times z → 0 and higher masses. This is expected because higher-mass haloes are resolved better at later times, and in addition, more haloes provide improved statistics. We expect that the agreement would be better for all times for larger boxes and better mass resolutions.

2. The relative differences in the HMFs for model I are shown for five different redshifts in the lower panel of Fig. 14. A general enhancement of ∼10% evolves to ∼100% at redshift z = 1. That the enhancement appears larger at the intermediate redshift z = 1 than at z = 0 is a surprising feature of the derived HMFs. Compared with the relative matter power spectra in Fig. 12, this feature is not clearly visible. A small suppression at large k-scales at later time can be detected. It might be speculated whether the suppression in the spectra and the HMFs can both relate to the domain collapse at redshift z ∼ 0.5 by adding kinetic energy to dark matter particles and disturbing screening and bound structures. The enhancement already begins to shrink before the collapse, so that this might not be the best explanation. Another explanation may be that the initial period of enhanced clustering leaves underdensities in the near neighbourhood of the filaments with respect to the ΛCDM; this would throttle the structure growth in the continuing period. This observation can be made from the density fields in Fig. 19, which we present in Sect. 5.5. The bump in the enhancement in Fig. 14 may be related to the similar enhancement scale viewed in the power spectra of Fig. 12

3. In Fig. 15, the relative enhancement of the HMFs compared to ΛCDM for models I–V is shown at redshift z = 0. A similar overall amplitude to model I is visible in Fig. 14 for models III–V, indicating that the enhancement might scale with $\mathrm{\Delta }A\sim \frac{1}{2}{(\frac{v}{M})}^2$, which was kept approximately constant between the two initialisation points. That is, with the exception of model II, for which ΔA is much smaller, and the HMF is suppressed by a factor of ∼β². This is similar to the case with the power spectra in Fig. 13.

4. Similarly to the matter power spectra, asymmetricity in the potential, for model V, seems to have very little effect on the HMF. A reduced Compton wavelength moved the enhancement towards smaller scales and lower masses, which reduced the level of enhancement of larger haloes. The reduced enhancement of higher-mass haloes would correspond to a more efficient screening mechanism.

Fig. 14. HMFs as a function of halo mass m200b in units of solar masses Msun/h of the different simulations. Top: model I HMFs shown for three different redshifts z = 2, 1, 0. The HMFs are plotted as dashed lines, compared with the analytic Tinker10 model for ΛCDM using the CLASS spectra as a solid line, and the result of the ΛCDM simulation as a dotted lines. Bottom: relative differences of the HMFs compared to their ΛCDM counterparts.

Fig. 15. Relative difference between the HMFs of simulations I–V and ΛCDM at redshift z = 0, plotted as a function of the halo mass m200b, in units of solar masses Msun.

5.4. Oscillating fifth force

In Fig. 16 from model III (L_C = 1 Mpc h⁻¹, z_* = 0.1,  β₊ = 8,  β₋ = 8), the panels labelled A-C show zoomed-in snapshots corresponding to the regions marked A-C in the top left panel. At the centre location of each zoomed-in panel, marked by a cross, circle, and triangle, we sampled the scalar field value throughout the simulation runtime. Part of the resulting time series of the scalar field values is shown in Fig. 17. In Fig. 18 we show a similar time series of a randomly selected location in model I. We list our observations below.

1. All of the locations sampled in Fig. 17 show a smooth low-frequency evolution, while Fig. 18 shows an extreme low-frequency evolution with several switches in the minimum and additional large-amplitude high-frequency oscillation modes.

2. Although sampling locations B and C in Fig. 16 both seem to reside in screened regions at redshift z = 0.6, it appears that neither remains screened. The magnitudes of their background values increase steadily, but the oscillation amplitudes are comparable and more or less constant in time.

3. Site C in Fig. 16 shows an effective screening initially, with variations in the field value of ≲1%.

Fig. 16. Different slices of the simulation volume of model III at redshift z = 0.6. Top left: snapshot of the scalar field. The black squares labelled A-C indicate the zoomed snapshots, which are displayed from the top right to the bottom right panels. The centre mark indicates the point from which the field value was sampled in Fig. 17.

Fig. 17. Field value χ in model III sampled from a single lattice point as a function of redshift z for the sites labelled A-C, shown in Fig. 16. The rows labelled dχ are the instantaneous field value χ subtracted the average of the field over the interval ±dτ/2 = ±0.016 Mpc h−1.

Fig. 18. Top panel: field value χ sampled from a single lattice point as a function of redshift z in model I. χ is the instantaneous field value χ averaged over the interval ±dτ/2, where dτ is either 1.6 or 4.7 Mpc h−1, as indicated in the legend. Bottom panel: difference between the field values and the field smoothed over 4.7 Mpc h−1. Blue shows the difference for the field value sampled at a 0.4 Mpc h−1 rate, and orange shows the difference for the field value averaged over ±1.6 Mpc h−1.

5.5. Planar structures

In Fig. 19, we show the formation of planar structures on cosmological scales by plotting the difference in the overdensity field between simulations of the (a)symmetron and ΛCDM models. To indicate the structures more clearly, we saturated the colour scale and display it in a symmetric logarithmic plot. We comment on this figure below.

1. The largest differences in the overdensities are found in the cosmic filaments for each simulation, where the clustering is stronger and the colour scale is mostly saturated.

2. The correspondence between the location of the planar structures and the domain walls in Figs. 10 and 8 is strong for models I and II–V, respectively.

3. In the bottom right panel, we show the density field of matter in the ΛCDM simulation ${\rho }_m/{\overline{\rho }}_m={\delta }_m+1$ for comparison. Several of the filaments that make up the skeleton of the domain walls in the (a)symmetron simulation are visible, which may indicate that domain wall pinning has a strong effect on the field configuration.

Fig. 19. Differences in the overdensity fields of the symmetron models with respect to ΛCDM, $\mathrm{\Delta }{\delta }_m\equiv {\delta }_m^{\text{symmetron}}-{\delta }_m^{\mathrm{\Lambda }\text{CDM}}$. The top left to bottom left panels show models I–V. The bottom right panel shows the matter density field ${\rho }_m/{\overline{\rho }}_m={\delta }_m+1$ for ΛCDM, where ${\overline{\rho }}_m$ is the background value of the density.

We discuss some more consequences of the planar overdensities in Sect. 6 and comment on relevant observational signatures.

5.6. Light propagation

We expect light to be redshifted by the expansion of space. Additionally, secondary effects arise from the interaction of the light with inhomogeneities that will perturb its redshift and direction. These secondary effects on light from a source field at z = 0.0835 are shown in Figs. 20 and 21. The source field is limited to a distance d(source)≤B/2 = 250 Mpc h⁻¹, where B is the box size of the simulation, corresponding to the chosen source redshift. The light-cone was constructed within gevolution on the fly, and the integration of the null geodesics from the source field was made by a shooting method that was presented in Adamek et al. (2019). In Fig. 20 we show the angular power spectrum of the secondary effects due to the (ISW-RS) effect, lensing potential, convergence (κ), Shapiro time delay, and the gravitational redshift (or the ordinary Sachs-Wolfe effect (SW)). These effects their computation, and observation were discussed in Hassani et al. (2020). Because the distance to the source, r < 250 Mpc h⁻¹, is small, the features of the angular power spectra correspond to physical scales, d, which are small, approximately

$$\begin{array}{c}\hfill l\simeq \pi /\mathrm{\Delta }\varphi \simeq r\pi /d,\end{array}$$(18)

Fig. 20. Secondary light propagation effects in the full-sky light cone of an observer at redshift z = 0 located at the origin position x0 = (0, 0, 0) Mpc h−1. The source field is located at redshift z = 0.0835. Top left: angular power spectra of the ΛCDM simulation shown as thin lines. The thick dotted line shows the CLASS estimate, which is only available for the ISW, convergence, and lensing potential. Top right to bottom right: relative angular power spectra with respect to ΛCDM for the simulations of models I–V. The meaning of the labels in the legend is SW: Ordinary Sachs-Wolfe effect. ISW-RS: integrated Sachs-Wolfe and non-linear Rees-Sciama effect. κ: convergence field. Shapiro: Shapiro time delay. ψ: lensing potential.

Fig. 21. Secondary light propagation effects in the full-sky light cone of an observer at redshift z = 0 located at the centre position x0 = (250, 250, 250) Mpc h−1. The source field is located at redshift z = 0.0835. Top left to bottom right: relative angular power spectra with respect to ΛCDM for the simulations of models II–V. The meaning of the labels in the legend is SW: ordinary Sachs-Wolfe effect. ISW-RS: integrated Sachs-Wolfe and non-linear Rees-Sciama effect. κ: convergence field. Shapiro: shapiro time delay. ψ: lensing potential.

for Δϕ ≪ 1, where l is the moment of the spherical harmonic corresponding to the azimuth. We first consider the light cone of an observer at redshift z = 0, located at the edge of the simulation domain x = (0, 0, 0) Mpc h⁻¹ (Fig. 20): In all of the late-time symmetry-breaking simulations, models II–V, the signature in the ISW signal is clear, giving a ∼10% relative enhancement with respect to ΛCDM and peaking at a characteristic scale of l ∼ 20 − 40. The signal is suppressed in model II, where β_* is smaller, but has an identical shape. The remaining secondary effects in the first light cone share a similar shape and peak around l ∼ 100, although the convergence field κ is slightly larger and peaks at larger scales l ∼ 40 − 100. We find the same qualitative behaviour in the second light cone of an observer at redshift z = 0 at the centre of the box x₀ = (250, 250, 250) Mpc h⁻¹, although the convergence of the angular power spectrum is larger. It peaks more sharply at l ∼ 40 − 50. The secondary effects on the light propagation in model I, for which we produced the first light cone alone, have a very different appearance with far weaker peaked angular spectra. The enhancements on the ISW-RS and convergence are still largest and reach ∼20% at small angular scales l ∼ 200. We place these observations into context from the literature and discuss some observational signatures in the discussions Sect. 6.

6. Discussions and conclusion

In this section, we discuss some of the observations made in the simulation data in Sect. 5 and their implications for the features of the models and signatures in observational surveys.

6.1. Observations made in simulation data

For domain partitioning and domain wall stability, considering the simulations listed in Table 1 and displayed in Fig. 7, and those shown in Figs. A.1 and A.2, we list our observations below.

1. The domain wall stability seems to be improved when L_C and z_* are decreased.

2. The domain wall stability and partitioning seems to be almost unaffected when $\overline{\beta }$ is varied. However, for larger $\overline{\beta }$s, we expect the secondary effect of enhanced clustering along the walls to which in turn the domain walls are pinned, which enhances stability.

3. Although we were unable to find stable domain walls for asymmetric potentials Δβ ≠ 0, we were able to make them quasi-stable on presumably arbitrarily long timescales by sufficiently decreasing L_C and z_*.

4. The domain wall instability and the effect on partitioning caused by the asymmetric potential seem to depend only on the relative level of asymmetry $\mathrm{\Delta }\beta /\overline{\beta }$. The stability was not affected when this was kept constant, but the βs were changed.

5. By decreasing the Compton wavelength, we find less variability in the fraction of positive minimum, meaning that the domain size has decreased so that we resolved more minima. This means that the fraction of positive minima is ∼0.5, as expected from the initial ℤ₂ symmetry.

These observations suggest that for this part of the parameter space, the pinning of domain walls to the filaments of large-scale structure is enhanced by the higher density contrasts in the filaments (later-time symmetry breaking) and by the interaction of the field with smaller-scale structures (smaller L_C). Intuitively, we would expect there to be a feedback process in which the enhanced clustering in the domain walls by itself would increase their stability by either increasing the filament density contrasts or creating wall structures where no filaments were originally. We see no strong indications that this were the case, but this effect might become more important at higher values of β. Some subtle differences between model II and III are visible in Figs. 8 and 9, such as the z = 0.2 collapse of the small centre-left positive-minimum island (approximately located at coordinates (0.1, 0.4)B, where B is the box size), indicating that this effect is not entirely absent.

For the explored late SSB parameter space of the symmetron, the mismatch n all instances of z_* ≲ 0.7 is increasingly large between the selected time of symmetry breaking at the background z_* and the actual time z_SSB, Δz = z_SSB − z_* for increasingly smaller z_*. This is very interesting because it means that 0 ≲ z_SSB ≲ 1 for a large part of the parameter space (z_*, L_C). When this part of the symmetron parameter space were linked to the dark energy phenomenon of late-time cosmological acceleration, the coincidence problem would be explained because z ≲ 1 marks the onset of high density contrasts. We will explore this possibility in more detail in a future work.

6.2. Signatures

In Sect. 5, we listed several results. We discuss their significance below.

The peaked ISW-RS signal in Figs. 20 and 21 is especially interesting in the context of the anamolous ISW signals that were claimed to be detected in large voids (see Kovács et al. 2017, 2019; Kovács 2018; Vielzeuf et al. 2020). Intuitively, we expect voids to be repulsive due to the fifth force F₅ ∼ −∂_rϕ², and we therefore expect an enhanced clustering towards the filaments that enclose voids compared to ΛCDM. Figure 19 shows that this indeed occurs: the filaments and domain walls all have enhanced densities (red), while regions immediately next to filaments or within voids have lower densities (light and dark blue). The effect on the evolution of the potentials from the different clustering is imprinted on the ISW signal. The important scale for this for models II–V based on Eq. (18) is ∼18 Mpc h⁻¹ for the mean distance of $R_{\mathrm{lightcone}}/\sqrt{3}\sim 144$ Mpc h⁻¹. Since this scale does not appear to depend significantly on ξ_* because there is no clear difference of it in model IV, we speculate that it may be set by the parameter a_* instead, or equivalently, by ρ_*. This might be related to the scale on which ${\delta }_u{\overline{\rho }}_m<{\rho }_{\ast }$, where δ_u is the underdensity of voids. This should be examined quantitatively. The width of the peak of the ISW signal instead seems to depend on the Compton wavelength and possibly also on a_*, although this needs to be verified in detail. In a future work, we will investigate the parameter dependence of this signal and study the ISW signal of voids in particular.

The planar structures that we identified in Fig. 19 are interesting in the context of similar features that were claimed to be observed in the local universe (see e.g. Peebles 2023). It is especially interesting for structures that were claimed to be detected on scales larger than what the standard ΛCDM should be able to account for, such as the detected dipoles in large-scale structure (Secrest et al. 2021, 2022; Panwar & Jain 2024) or other claimed structures such as the Hercules-Corona Borealis (HCB) Great Wall (Horvath et al. 2015, 2020). The structures that appear in our simulations (Fig. 19) extend across the box of ∼500 Mpc h⁻¹ and remain approximately flat on scales up to 100–250 Mpc h⁻¹ in some cases. Their overdensity amplitude can be adjusted by a change in $\overline{\beta }$, and their scale seems surprisingly not to depend strongly on the parameter choices within our explored parameter space. The scale in the late-time simulations seems to come from the cosmic filaments and the scale at which they encapsulate domains. We expect variations in this scale for more extreme variations in the parameter space. We leave this study to a future work. The planar scales we found in the simulations are comparable to those reported in Peebles (2022), who reported flat sheet-like structures with a thickness of ∼5 Mpc and lengths of ∼150 Mpc. An interesting future work would be to study the ability of the (a)symmetron parameter space to account for such observations in detail, over a larger parameter space, and using CDM haloes instead of particles. The haloes would correlate more stronger with the galaxy number count statistics made use of in Peebles (2023). Furthermore, it would be interesting to analyse the full three-dimensional data to learn more about the sheet structure morphologies.

The observation of local fluctations in the scalar field, and thus also the fifth force, F₅ ∼ ∂_r(ϕ²), in Figs. 17 and 18 may have distinct observational signatures. The effect of oscillating ultralight scalar fields has been considered in the literature for coupling to the dynamics of pulsars (Ünal et al. 2024) and to the oscillations of stars (Sakstein & Saltas 2024), which would help to impose strict constraints on the allowed energy density parameter of it when considered as a dark matter candidate. In addition to the effect of an oscillating gravitational potential owing to the energy density of the scalar field, we have the oscillating fifth force in this instance, which might work towards imposing stricter constraints on the model parameter space. However, screening must be taken into account, which was shown in Fig. 17 to dramatically suppress the oscillation and the overall scalar field value for the duration of the screening activity. In the past, the possibility of disrupting the screening mechanism from so-called cosmic tsunamis was reported (Hagala et al. 2017). These are waves in the scalar field. A study of such disruption events in the simulation output will be made in the future. In addition to pulsars and stars, periodic systems on different length scales may help to constrain dark matter masses and the fifth force couplings that operate on these scales. In our simulations, the oscillation has a timescale of ∼Mpc h⁻¹, in which case it might affect the dynamics of galaxy clusters. Likewise, we expect the oscillations in models with L_C ∼ kpc h⁻¹ to be able to affect the dynamics of galaxies. For example, when the attractive potential of a cluster changes, the virial radius oscillates and prevents the particles from having stable orbits; for one instance, the force increases, giving the particles a higher centripetal acceleration and injecting kinetic energy; and in the next moment, the force decreases, and the particles become unbound and spread out. This was observed in Dutta Chowdhury et al. (2021) for an oscillating potential, and we expect the effect to be enhanced when a fifth force is included. It would be an interesting project to consider whether this dynamic would be able to account for the velocity dispersion of stars in galaxies for choices of ∼kpc h⁻¹ scale general Compton wavelengths. When the screening is included, the effects on systems with AU timescales might be considered for shorter Compton wavelengths to constrain lunar ranger experiments (Xue et al. 2020).

Enhanced clustering owing to the fifth force is visualised in the matter power spectra in Figs. 12 and 13 and in the HMFs in Figs. 14 and 15 for models I and I–V, respectively. Overall, the effect of the scalar field is to enhance clustering, but this is not uniform over the scales k; the relative spectral enhancements in Fig. 13 are clearly peaked and then suppressed for larger and smaller scales. At larger scales, we expect the suppression to be caused by the horizon H⁻¹, and on smaller scales by screening δ_m(k)ρ_m > ρ_*. The enhanced clustering presents an obstacle to the model, as the effect is likely to worsen the σ₈ tension reported in the literature (McCarthy et al. 2023), which would require a stronger mechanism to suppress structure formation in the late-time universe than what is needed for ΛCDM. The σ₈ tension is only ∼2 − 3σ in significance so far, which means that it is not clear whether it will persist as more data are included, however. Observations that are more consistent with enhanced clustering are also indicated in the literature (Mazurenko et al. 2023). Speculatively, thinking about suppression of clustering, a mechanism on cosmological scales might be expected for a model where the (a)symmetron makes up dark matter and has a conformal coupling so that it has a ∼Mpc h⁻¹ scale Compton wavelength scale in screened regions during the time between the last scattering surface and now. In this instance, it would operate similarly to fuzzy dark matter (Niemeyer 2020) cosmologically. With growing inhomogeneities at late times, higher-mass dark matter particles might be accounted for and might simultaneously have a Compton wavelength of ∼kpc h⁻¹ in galaxies today. Furthermore, the correlation ranges of dark matter particles in unscreened regions will correspond to the choice of the Compton wavelength, which can be chosen at cosmological scales ∼Mpc h⁻¹. In this instance, suppression and enhancement can both be derived depending on the environment. These suggestions require a more careful analysis and will be explored in a future work. For the simulations considered here, the energy budget of the symmetron was chosen to be almost negligible, and reliably simulating fuzzy symmetron dark matter would present an additional computational challenge that would be an interesting endeavour to consider. The shrinking relative enhancement observed in the HMF of model I presents another opportunity for not worsening the σ₈ tensions as much. Structure formation may continue into a period of weaker clustering due to a depletion of matter in neighbouring environments.

6.3. Conclusion

We have presented the first cosmological simulations of the (a)symmetron for which spatial and temporal convergence in the field configuration was demonstrated, and we retrieved quantitative estimates of novel signatures such as secondary light-propagation effects, the formation of sheet-like structures imprinted on the dark matter distribution, and local oscillations in the fifth force. We analysed the evolution of the power spectra and the HMFs. In the latter, we observed a feature of greater enhancement at intermediate times, the cause of which will be interesting to determine in the future. The possibility of producing simulations that fully resolve the dynamic part of the symmetron parameter space opens several possible paths towards uncovering new phenomenological features of the model and further constraining its parameter space. In the future, we wish to study the ISW-RS signal in more detail, its correlation with void structures, and in particular, the evolution of the gravitational potential of large voids, taken in the context of the anomalous ISW signal found in the literature (Kovács et al. 2017). We wish to study the effect of oscillating potentials on the radial mass profiles of dark matter haloes in greater detail, and we would like to conduct a more extensive study of the imprint of the domain walls on the dark matter density field, where it creates planar structures, taken in the context of Peebles (2023). In the future, a larger project would be a simulation of the (a)symmetron in the part of its parameter space where it accounts for dark matter, which is well motivated in the literature (Burrage et al. 2017, 2019; Käding 2023), in addition to our findings presented here and in the companion paper Christiansen et al. (2024). The (a)symmetron model is an interesting and phenomenologically complex model that includes mechanisms for solving many of the observational tensions of modern day cosmology, along with a large array of new and unique predictions. With the improvement of computational resources and techniques, we expect to access predictions in a larger part of its parameter space in the near future.

Movies

Movie 1 (4wayplot_1280_sym1) Access here

Movie 1 (4wayplot_1280_sym2) Access here

Movie 1 (4wayplot_1280_sym2_asym) Access here

Movie 1 (4wayplot_1280_sym2_compton) Access here

Movie 1 (4wayplot_1280_sym2_time) Access here

Movie 1 (halfboxdomain_ngrid128_ideal_one) Access here

Acknowledgments

We would like to thank Julian Adamek and Jens Jasche for helpful discussions or comments. Ø.C. would like to thank Constantinos Skordis and the Central European Institute of Cosmology (CEICO) for hosting him during the final stages of the project. F.H. wishes to thank Martin Kunz and Jean-Pierre Eckmann for their insightful discussions and hosting him during the project’s early stages. We thank the Research Council of Norway for their support. The simulations were performed on resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway. This work is supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s1051.

References

Appendix A: Parameter space

In preparation of the more expensive simulations presented in section 5, we ran a large number of N_grid = 512³, (500 Mpc/h)³ volume simulations in order to identify an interesting region in the parameter space that permitted the existence of quasi-stable domain walls throughout the simulations. As a simple quantity to quantify the stability and behaviour of the domain walls, we compared the fraction of the simulation volume in the least occupied domain (labelled + here). The result is shown in figures A.1 and A.2. By choosing a resolution that barely resolved the Compton wavelength dx = 0.9 ⋅ L_C, but keeping the size of the simulation domain L_box = 500 Mpc/h, we expected the results to be representative for the dynamics in the higher-resolution simulations of section 5, according to section 4 and figure 4. Several of the observations listed in section 6.1 on the effect of variation of the symmetron parameters on the domain wall stability can be made from the figures A.1 and A.2.

Fig. A.1. Fraction of the simulation volume with the scalar field in the positive minimum as a function of redshift for a range of (a)symmetron parameters. The last three lines in pink all show the same parameter choice, but at an earlier initialisation point as the dotted line, and with a different seed number as the long dash line.

Fig. A.2. Fraction of the simulation volume with the scalar field in the positive minimum as a function of conformal time before collapse (normalised by the Hubble scale) for a range of (a)symmetron parameters. In the case of no collapse (densely dotted purple), it is instead plotted against the conformal time until the end of the simulation z = 0.

Appendix B: Initialisation

As was demonstrated in subsection 4.1, resolving the Klein-Gordon equation in the pre-symmetry-broken time at early redshifts z ∼ 100 would be computationally demanding, both because the Compton scale is smaller due to the high mass $m\sim \sqrt{-T}$ and because we would have to evolve the symmetron for a longer duration. In subsection 4.1, we argued that a late-time initialisation was sufficient and would save computational expenses. We substantiate this argument here.

The initialisation that was used in Christiansen et al. (2023) and that was also used here, initialises the scalar field as a draw from a Gaussian scale-invariant power spectrum of some small amplitude set by the user A_amp ∼ 10⁻²⁰. The following evolution of the field in the pre-symmetry-broken regime involves local oscillations on spatial scales of the Compton scale L. The energy scale made use of here for the symmetron is ∼10⁻⁶ ⋅ ρ_c, 0 at most post-symmetry breaking, and it is completely negligible before ∼10⁻⁴⁰ ⋅ ρ_c, 0, so that we neglect no relevant gravitational effects. The fifth force is ∝ϕ∂_rϕ, and it is therefore also suppressed because ϕ ∼ A_amp. The only information that might be lost by dismissing the early time field evolution is then the phase information of the field and the correlation structure that would be established throughout this time. In figure B.1 we demonstrate that the initially scale-invariant power spectrum of the scalar field remains very nearly scale-invariant up until the symmetry-breaking time, and the symmetry-breaking process otherwise takes place identically. Similar features for earlier-time initialisations are observed, although there is an amplitude modulation owing to the Hubble friction acting on the field. Since both the amplitude and initial phase of the scalar field were chosen arbitrarily, the information of its evolution in the pre-symmetry-broken time is not relevant as long as it leaves its statistics qualitatively similar. Small variations in the amplitude seems to affect the symmetry-breaking process only very little, and variations in the phase amount to changing the seed number of the simulation.

Fig. B.1. Pre-symmetry-breaking evolution of the scalar power spectrum for a simulation of 5123 grids and (500 Mpc h−1)3 volume. The left panel has an initialisation time at redshift zini = 1.5, and the right side panel has zini = 1.

For a late-time initialisation of the scalar field, a relevant question is how late the initialisation can be made. At this point, we discovered that there is some pre-symmetry-breaking evolution that needs to be resolved. Before the symmetry breaking occurs globally, the density contrast in the simulation causes |T|< ρ_* locally. The spatial Laplacian in the wave equation prevents the field from symmetry breaking in this region until the density has dropped below ρ_* for a sufficiently large scale that seems to be related to the Compton wavelength L_C. However, during the time before this occurs, the domain structure establishes itself, even though the field amplitude is still very small (see figure B.2). This evolution is relevant to resolve, and when we did not resolve it, the domain walls were more unstable in some instances. Additionally, because the field is expected to be correlated on the Compton scale L, which increases before the symmetry breaking (see figure 1), we expect to have to initialise the simulation so that the field thermalises on this scale before the symmetry breaking. We chose both initialisation times z_in = 3, 1 in section 5 by running 512³ grid (500 Mpc h⁻¹)³ volume simulations, ensuring that we started the simulations before the domain structure established, and finding convergence with yet later initialisation times. We expect the time required before symmetry breaking to depend on the density contrasts, among other parameters, and we leave the analytic prediction of this for the future.

Fig. B.2. Evolution in the scalar field, χ, of model I before the global symmetry breaking, which occurs at zSSB ∼ 1.9. The right panel shows the equation-of-state parameter of the field, which evolves until the global symmetry breaking and indicates the final domain walls.

Appendix C: Solvers

The solvers we chose were expensive in different ways. Using a method such as the fourth-order Runge-Kutta method required us to keep multiple auxiliary fields to avoid overwriting the field values before the neighbour lattice points had calculated their Laplacians. In addition, every partial step had to be done within its own loop over the lattice, giving more overhead. Finally, between each step, the auxiliary field that was used in the next loop Laplacian must have updated ghost layers of the computational domains before entering the loop, giving more communication overhead. Taking all of these obstacles into account, we compare the convergence of both solvers with respect to the wall-clock time and with the number of iterations in figure C.1. The comparison is made over a redshift interval of z ∈ (1.5, 1.46) for the symmetron model parameters (L_C, a_*, β_*) = (1Mpc/h, 0.33, 1). At each step i, the error was calculated as the average absolute difference between the solution obtained with a step size h and the previous with step size h/2 as

$$\begin{array}{c}\hfill \mathrm{err}\equiv ⟨|{\chi }_i-{\chi }_{i-1}|⟩.\end{array}$$(C.1)

Fig. C.1. Comparison of different solvers for the symmetron in a 3 Gpc/h box and 1283 grids. The convergence is checked in a redshift interval z ∈ (1.5, 1.46). The final slopes of the solver in logspace are 1.97 for the leapfrog, 4.3 for the fourth-order Runge-Kutta, 1.91 for PEC2k1, and 0.87 for the Euler solver.

Figure C.1 shows the expected convergence with smaller step sizes. The Runge Kutta method is of the fourth order, meaning err ∼ 𝒪(h⁴). This easily recognised by realising that if err ∼ 𝒪(hⁿ), then

$$\begin{array}{c}\hfill {\chi }_i(2h)-\chi (h)=\chi -\chi +[\mathrm{err}(2h)-\mathrm{err}(h)]\sim (2^n-1)h^n,\end{array}$$(C.2)

so that the function in logspace is log(err)∼n log h.

In figure C.2 we demonstrate this convergence in the relative error of the actual scalar field value for a variety of choices of scalar Courant factors for cosmological time evolutions using the fourth-order Runge Kutta solver. Because of the divergence for χ ∼ 0, we only averaged the points on the lattice for which χ > 0.1.

Fig. C.2. Comparison of different Courant factors for the scalar field in a 512 Mpc/h box and 1283 grids. The comparison is made with respect to the solution of the scalar field found with a Courant factor of 10−4. The percent-level agreement is indicated with the dashed horizontal line.

All Tables

All Figures

	Fig. 1. General Compton wavelength of the symmetron field at the background as a function of the scale factor. The parameter choice is (LC, a*, β*) = (1 Mpc h−1, 0.33, 1), and the horizontal lines indicate the initial and final Compton scales.
In the text

	Fig. 2. Initial setup for the bare wall spatial and temporal convergence experiments. The right panel shows the overlaid analytic profile corresponding to a Compton wavelength of 0.5 Mpc h−1, although the model here is defined with a Compton wavelength of 1 Mpc h−1.
In the text

	Fig. 3. Volume fraction of the field that is in the positive minimum as a function of redshift for the oscillating bare wall system for variations of the Courant factor Cϕ. The four plots vary in spatial resolution dx, initial wall scale LIC, and Compton wavelength LC.
In the text

	Fig. 4. Volume fraction of the field in the positive minimum as a function of redshift for the oscillating bare wall system for a selection of spatial resolutions. The coarsest simulation has the Courant factor choice Cϕ = 0.5 and would otherwise not resolve the oscillations. The rest are kept with Cϕ = 1.
In the text

	Fig. 5. Comparison of snapshots of the scalar field χ at redshifts z = 1.9 and 1.8 for a late initialisation at z = 2.5 and for an early initialisation at z = 100. The amplitude of the late-time initial conditions is chosen to be similar to the evolved early-time conditions. In the bottom left snapshot, a black square indicates the relative size of a (500 Mpc h−1)3 volume box.
In the text

	Fig. 6. Fraction of nodes in positive minimum for simulations run with initial conditions found from a ΛCDM simulation with a Courant factor Ccdm = 1. The box size is B = 410 Mpc h−1, and to the right, we show the fraction of domains within a B = 205 Mpc h−1 subvolume. We used 5123 grids. We chose Cϕ = 0.15.
In the text

	Fig. 7. Overview of background quantities for the suite of simulations. Top: fraction of the simulation volume where the scalar field occupies the positive minimum for the five different models I–V. Bottom: average over the simulation volume of the square of the scalar field χ2. Both quantities are plotted as a function of the redshift z.
In the text

	Fig. 8. Snapshots of the scalar field χ. Shown for models II–V from left to right, for four different redshifts from top to bottom.
In the text

	Fig. 9. Snapshots of the velocity of the scalar field $q=a^3\dot{\chi }$. Shown for models I–V from left to right, for four different redshifts from top to bottom. The field is normalised to the conformal Hubble parameter.
In the text

Fig. 10. Snapshots of slices of the simulation volume. From top left to bottom right: the scalar field χ, the isotropic $\mathrm{\nabla }\chi \equiv \left|\mathrm{\nabla }\chi \right|/\sqrt{3}$ normalised to the one expected for the centre of the analytic bare wall, the equation-of-state parameter of the field ωϕ, and the energy density parameter Ωϕ = ρϕ/ρc. The snapshots are from model I and are shown at redshift z = 0.9 at a critical moment right before the two blue domains join.

In the text

	Fig. 11. Snapshots of the simulation volume. From top left to bottom right: the scalar field χ, the isotropic gradient ∇χ≡$\left|\mathrm{\nabla }\chi \right|/\sqrt{3}$ normalised to the conformal Hubble parameter ℋ, the equation-of-state parameter of the field ωϕ, and the energy density parameter Ωϕ = ρϕ/ρc. The snapshots are from model IV.
In the text

Fig. 12. Matter power spectra for model I as a function of wavenumber k indicated for different redshifts z. Top: overlaid power spectra from the symmetron asevolution simulation plotted as a solid line, and the corresponding ΛCDM simulation is shown as dots. The CLASS power spectrum for the same ΛCDM parameters is shown as a dashed line. Bottom: relative difference with respect to ΛCDM of the power spectra from redshift z = 3 to z = 0 indicated in graded colours from blue to yellow. The gap in the spectra in the lower panel is not a physical feature and is caused by a typo in the settings file of the simulation. This caused two spectra to be skipped.

In the text

	Fig. 13. Relative difference of matter power spectra with respect to the ΛCDM case, with redshift z = 1 to z = 0 indicated in graded colours from blue to yellow. The spectra are plotted as a function of wave number k. The top left to bottom right panels show models II–V.
In the text

Fig. 14. HMFs as a function of halo mass m200b in units of solar masses Msun/h of the different simulations. Top: model I HMFs shown for three different redshifts z = 2, 1, 0. The HMFs are plotted as dashed lines, compared with the analytic Tinker10 model for ΛCDM using the CLASS spectra as a solid line, and the result of the ΛCDM simulation as a dotted lines. Bottom: relative differences of the HMFs compared to their ΛCDM counterparts.

In the text

	Fig. 15. Relative difference between the HMFs of simulations I–V and ΛCDM at redshift z = 0, plotted as a function of the halo mass m200b, in units of solar masses Msun.
In the text

	Fig. 16. Different slices of the simulation volume of model III at redshift z = 0.6. Top left: snapshot of the scalar field. The black squares labelled A-C indicate the zoomed snapshots, which are displayed from the top right to the bottom right panels. The centre mark indicates the point from which the field value was sampled in Fig. 17.
In the text

	Fig. 17. Field value χ in model III sampled from a single lattice point as a function of redshift z for the sites labelled A-C, shown in Fig. 16. The rows labelled dχ are the instantaneous field value χ subtracted the average of the field over the interval ±dτ/2 = ±0.016 Mpc h−1.
In the text

Fig. 18. Top panel: field value χ sampled from a single lattice point as a function of redshift z in model I. χ is the instantaneous field value χ averaged over the interval ±dτ/2, where dτ is either 1.6 or 4.7 Mpc h−1, as indicated in the legend. Bottom panel: difference between the field values and the field smoothed over 4.7 Mpc h−1. Blue shows the difference for the field value sampled at a 0.4 Mpc h−1 rate, and orange shows the difference for the field value averaged over ±1.6 Mpc h−1.

In the text

Fig. 19. Differences in the overdensity fields of the symmetron models with respect to ΛCDM, $\mathrm{\Delta }{\delta }_m\equiv {\delta }_m^{\text{symmetron}}-{\delta }_m^{\mathrm{\Lambda }\text{CDM}}$. The top left to bottom left panels show models I–V. The bottom right panel shows the matter density field ${\rho }_m/{\overline{\rho }}_m={\delta }_m+1$ for ΛCDM, where ${\overline{\rho }}_m$ is the background value of the density.

In the text

Fig. 20. Secondary light propagation effects in the full-sky light cone of an observer at redshift z = 0 located at the origin position x0 = (0, 0, 0) Mpc h−1. The source field is located at redshift z = 0.0835. Top left: angular power spectra of the ΛCDM simulation shown as thin lines. The thick dotted line shows the CLASS estimate, which is only available for the ISW, convergence, and lensing potential. Top right to bottom right: relative angular power spectra with respect to ΛCDM for the simulations of models I–V. The meaning of the labels in the legend is SW: Ordinary Sachs-Wolfe effect. ISW-RS: integrated Sachs-Wolfe and non-linear Rees-Sciama effect. κ: convergence field. Shapiro: Shapiro time delay. ψ: lensing potential.

In the text

Fig. 21. Secondary light propagation effects in the full-sky light cone of an observer at redshift z = 0 located at the centre position x0 = (250, 250, 250) Mpc h−1. The source field is located at redshift z = 0.0835. Top left to bottom right: relative angular power spectra with respect to ΛCDM for the simulations of models II–V. The meaning of the labels in the legend is SW: ordinary Sachs-Wolfe effect. ISW-RS: integrated Sachs-Wolfe and non-linear Rees-Sciama effect. κ: convergence field. Shapiro: shapiro time delay. ψ: lensing potential.

In the text

	Fig. A.1. Fraction of the simulation volume with the scalar field in the positive minimum as a function of redshift for a range of (a)symmetron parameters. The last three lines in pink all show the same parameter choice, but at an earlier initialisation point as the dotted line, and with a different seed number as the long dash line.
In the text

	Fig. A.2. Fraction of the simulation volume with the scalar field in the positive minimum as a function of conformal time before collapse (normalised by the Hubble scale) for a range of (a)symmetron parameters. In the case of no collapse (densely dotted purple), it is instead plotted against the conformal time until the end of the simulation z = 0.
In the text

	Fig. B.1. Pre-symmetry-breaking evolution of the scalar power spectrum for a simulation of 5123 grids and (500 Mpc h−1)3 volume. The left panel has an initialisation time at redshift zini = 1.5, and the right side panel has zini = 1.
In the text

	Fig. B.2. Evolution in the scalar field, χ, of model I before the global symmetry breaking, which occurs at zSSB ∼ 1.9. The right panel shows the equation-of-state parameter of the field, which evolves until the global symmetry breaking and indicates the final domain walls.
In the text

	Fig. C.1. Comparison of different solvers for the symmetron in a 3 Gpc/h box and 1283 grids. The convergence is checked in a redshift interval z ∈ (1.5, 1.46). The final slopes of the solver in logspace are 1.97 for the leapfrog, 4.3 for the fourth-order Runge-Kutta, 1.91 for PEC2k1, and 0.87 for the Euler solver.
In the text

	Fig. C.2. Comparison of different Courant factors for the scalar field in a 512 Mpc/h box and 1283 grids. The comparison is made with respect to the solution of the scalar field found with a Courant factor of 10−4. The percent-level agreement is indicated with the dashed horizontal line.
In the text