3.1. Convergence Is Not Sufficient

Confidence in predictions about the abundance of DM subhalos extracted is usually explained by the fact that they seem unaffected by an increase of resolution above 50–100 particles per subhalo (see, e.g., Springel et al. Reference Springel, Wang, Vogelsberger, Ludlow, Jenkins, Helmi, Navarro, Frenk and White2008). As mentioned above, however, this prediction is very sensitive to modeling assumptions. Up to the close of the last century, simulations were suffering from an ‘overmerging’ problem, in that not enough substructure was predicted to match the observations. Several culprits had been proposed back then, with no consensus on the cause of the subhalos disruption but with an agreement that it was a numerical problem. Moore, Katz, and Lake (Reference Moore, Katz and Lake1996) blamed inadequate force softening.Footnote ⁶ Carlberg (Reference Carlberg1994), however, argued that a low mass resolution could cause two-body heating and artificially enhance matter disruption. Since this problem was superseded by a ‘missing satellite’ problem as the resolution of simulations increased, no definitive conclusion was drawn about the cause of the overmerging. Yet, as shown by van den Bosch et al., subhalos disruption is still “extremely prevalent in modern simulations, with … ∼ 65 percent of all subhaloes accreted around $z=1$ … disrupted by $z=0$” (Reference van den Bosch, Ogiya, Hahn and Burkert2017, 3044). The question thus arises: Are N-body simulations still suffering from overmerging, or can they be considered reliable, on the basis of convergence studies? Is the subhalos disruption a physical mechanism or the result of numerical artifacts?

This question is addressed by van den Bosch and Ogiya (Reference van den Bosch and Ogiya2018), in a paper also aiming at gaining a better understanding of the nonlinear effects of tidal stripping on subhalos.Footnote ⁷ Tidal processes are very difficult to describe, as the stripping of matter causes the subhalo remnant to fall out of virial equilibrium and then to revirialize by expanding, thereby provoking more stripping of matter, and to fall out of equilibrium once again.Footnote ⁸ No analytical theory is available to describe such complicated processes of de- and revirialization: only simulations can tell how the density distribution of the subhalo is affected by tidal stripping. It is thus crucial to disentangle what pertains to this stripping mechanism and what artificially results from an inadequate parametrization in simulations.

In order to do so, the authors came up with an idealized scenario in which the physical hypothesis of tidal stripping could be tested against the numerical hypothesis that matter disruption is mostly caused by inadequate force softening. They focused on the tidal evolution of an isolated subhalo on a circular orbit, given that, in such a scenario, only the host halo’s tidal field and the numerical parameters can affect the subhalo disruption. Its isolation forbids high-speed encounters with other subhalos causing more matter to disrupt, and the circular orbit excludes a fast pericentric passage that would deform the orbiting body and cause internal heating and thus more mass loss. The goal is to isolate two possible causes of the disruption and find the culprit by looking at the evolution of the bound mass of subhalos when varying the strength of the tidal field and the force softening.

Figure 4 in van den Bosch and Ogiya (Reference van den Bosch and Ogiya2018, 4071) suggests that increasing the strength of the tidal field by decreasing the orbital radius does not lead to the disruption of subhalos, except in the extreme case in which the orbital radius is chosen so that $r_{\mathrm{orb}}/r_{\mathrm{vir},h}=0.1$.Footnote ⁹ However, at $r_{\mathrm{orb}}/r_{\mathrm{vir},h}=0.1$, the subhalo totally disrupts after ∼ 13 Gyr. This results concords with the results of the Millennium and the Bolshoi simulations according to which $r_{\mathrm{orb}}/r_{\mathrm{vir},h}=0.1$ is the typical radius at which subhalos undergo disruption. If this prediction is a physical prediction, as it is taken to be by most astrophysicists, it should not be significantly affected by variations in numerical parameters.

Now, look at the impact of varying force softening and mass resolution (van den Bosch and Ogiya Reference van den Bosch and Ogiya2018, fig. 7, 4074). According to Power et al. (Reference Power, Navarro, Jenkins, Frenk, White, Springel, Stadel and Quinn2003), for instance, the optimal value for the force softening ε_opt ranges between 0.02 and 0.06. Since a force softening ε that is too large results in smaller density distribution at small radii, and since less dense systems are more exposed to tidal stripping, for $\epsilon >{\epsilon }_{\mathrm{opt}}$ one expects enhanced stripping and matter disruption. Likewise, a small softening is more exposed to two-body relaxation effects that flatten the central density profile and enhance disruption. Thus, force softening should have no other consequence than increasing the mass loss. Surprisingly however, the simulations show the opposite: for $\epsilon <{\epsilon }_{\mathrm{opt}}$, the bound remnants are larger and survive longer. One can suspect then that the physical disruption is very sensitive to the value assigned to ε.

Consider now simulations at the orbital radius $r_{\mathrm{orb}}/r_{\mathrm{vir},h}=0.1$ when the mass resolution and the force softening are increased together (see fig. 1). If you look for a region in the parameter space where the instability caused by mass resolution is kept under control and a ‘converged’ fraction of bound mass f _bound(t) can be found, the upper-left corner of figure 1 gives you the closest result to such an ideal.Footnote ¹⁰ But there, subhalos do not disrupt after 13 Gyrs. On the contrary: a large fraction of bound mass survives.

Figure 1.

Bound fraction of mass as function of time for simulations with different mass resolution N_p and different force softening ε. The black line shows the ‘converged’ results of a simulation with $N_p={10}^7$ and $\epsilon =0.003$; the blue line, the results from 10 simulations; the red line, their average. (Source: van den Bosch and Ogiya [Reference van den Bosch and Ogiya2018], fig. 10, 4077.)

Note the contrast between this result at large N_p and small ε and those observed along the yellow-shaded band, which corresponds to the scaling between mass resolution and force softening defined by Power et al. (Reference Power, Navarro, Jenkins, Frenk, White, Springel, Stadel and Quinn2003). This scaling is obeyed by most of the state-of-the-art simulations since 2003. If one focuses on the red lines (i.e., on the averaged results of simulations), the bound fraction of mass appears converged: the red line prediction stays more or less the same despite increasing the mass and the force resolution, which is precisely why it is considered as robust by Power et al.’s standards. Yet, this red line indicates a full disruption of the subhalo between 5 and 8 Gyrs. Thus, simulations converge on predicting that subhalos fully disrupt after 8 Gyrs but also on predicting that they survive after 13 Gyrs.

What should alert the philosopher here is that convergence alone will not tell which of these results is the correct one or whether one of them is. Given the range of the optimal softening defined by Power et al. (Reference Power, Navarro, Jenkins, Frenk, White, Springel, Stadel and Quinn2003), increasing the mass resolution and the force softening gives good agreement among simulations and confidence in the prediction that subhalos at $r_{\mathrm{orb}}/r_{\mathrm{vir},h}=0.1$ will not survive after one Hubble time. But if another region of the parameter space is scrutinized, convergence is found on a prediction that contradicts the former. As summarized by van den Bosch and Ogiya, convergence is “not a sufficient condition to guarantee that the results are reliable” (Reference van den Bosch and Ogiya2018, 4067).

3.2. Pseudoconvergence or Convergence?

In the previous subsection, I argued that robustness analysis in the form of convergence studies is not sufficient to exclude numerical artifacts and that robust predictions cannot be considered reliable on such grounds. Here, I push this thought further and argue that convergence can result from artifacts, following the work of Anton Baushev on the cusp-core problem.

In N-body simulations, real DM is represented by a limited number of heavy test bodies, to make the computational task tractable while preserving the averaged density of the system. However, whereas DM is collisionless, heavy test bodies undergo collisional effects, which in turn affect the density profile of DM halos. These collisional effects are characterized by the relaxation time ${\tau }_r=N\left(r\right)/(8\ln \Lambda .{\tau }_d)$, with N(r) the number of test bodies inside a sphere of radius r, lnΛ the Coulomb logarithm, and τ_d the characteristic dynamical time of the system at radius r. This is why Power et al. (Reference Power, Navarro, Jenkins, Frenk, White, Springel, Stadel and Quinn2003) recommend, as one of their convergence criteria, that convergence be reached in a time $t\le 1.7{\tau }_r$ (see n. 5).

Baushev, however, has shown that cuspy profiles are generated by an intensive energy relaxation. This fact raises important questions. DM halos undergo violent relaxation when they collapse, as density inhomogeneities create small-scale gravitational fields mediating the exchange of energy among DM particles. But they do so only at the moment of the collapse. The halo, once formed, has a stationary gravitational field. Why is it, then, that N-body simulations based on the CDM model predict cuspy density profiles for formed halos? If the cuspy profile stems from collisionality, where does it come from? Is it artificial or physical?

Baushev et al. (Reference Baushev, Valle, Campusano, Escala, Muñoz and Palma2017) address this question with a methodology very similar to that of van den Bosch and Ogiya (Reference van den Bosch and Ogiya2018). First, they propose an idealized scenario in which all sources of collisionality are turned off except the one under scrutiny. Thus, they simulate an isolated halo, to avoid the tidal influence of nearby halos and the gravitational capture of more mass from the surrounding, potentially leading to a secondary violent relaxation. They chose a halo with a Hernquist density profile that behaves like the Navarro et al. (Reference Navarro, Frenk and White1997) profile in the central region but is fully stable. Given that such a halo has a constant gravitational potential field ϕ(r), the density and velocity profiles should remain the same and all (the implicit functions of) the integrals of motion such as the specific energy $\omega =\varphi \left(r\right)+(v^2/2)$, the specific angular momentum $\overrightarrow{K}$, and the apocenter distance r ₀ should be conserved.Footnote ¹¹ Thus, any variation of them must be due to numerical effects. Such symptoms of artificial collisions can be tracked in their simplified scenario. Figure 2 of Baushev et al. (Reference Baushev, Valle, Campusano, Escala, Muñoz and Palma2017, 6) shows the variations of the integrals of motion as a function of the radius over 200 snapshots and confirms that they vary significantly even in a single time step. Thus, the convergence criterion of Power et al. (Reference Power, Navarro, Jenkins, Frenk, White, Springel, Stadel and Quinn2003) by no means guarantees that all sources of numerical artifacts have been excluded.

This preliminary conclusion raises another interesting question: Why does the density profile remain stable if the integrals of motion vary? In other words, is the stability of the cuspy profile a mere coincidence or the result of the numerical effects observed above? If the system is collisional, then it is better modeled by the Fokker-Planck equation, which models dynamical friction and diffusion, than by the collisionless Boltzmann equation usually appealed to. In that case, the diffusion streams created by the particle interactions could contribute in forming a stable density profile. And, as it turns out, the Fokker-Planck equation does have a stationary solution close to Navarro et al.’s (Reference Navarro, Frenk and White1997). If the cusp is a product of this Fokker-Planck diffusion, then artificial collisions should form a downward and an upward stream of particles, with increasing and decreasing r ₀ respectively, that complement each other and thereby explain the stability of the profile. This prediction can be tested by taking two adjacent snapshots S ₁ and S ₂ and calculating, for an array of radii r, the number of particles ΔN ₊(r) that had $r_0<r$ at S ₁ and $r_0>r$ at S ₂ and the number of particles ΔN ₋(r) that had $r_0>r$ at S ₁ and $r_0<r$ at S ₂. Figure 4 in Baushev et al. (Reference Baushev, Valle, Campusano, Escala, Muñoz and Palma2017, 9) confirms that the Fokker-Planck streams are present, complement each other very well outside a given radius, and are important enough to shape the density profile.

In sum, the variations of the integrals of motion unambiguously demonstrate that converged N-body simulations still suffer from artificial collisionality. Convergence fails to identify results independent of artifacts. Moreover, confidence in the density profile is not warranted by its stability since this stability is itself produced by virtue of numerical artifacts.

3.3. Discussion

Although convergence failed in the two cases discussed above, in diagnosing artifacts, one could defend robustness by arguing that Baushev et al.’s (Reference Baushev, Valle, Campusano, Escala, Muñoz and Palma2017) and van den Bosch and Ogiya’s (Reference van den Bosch and Ogiya2018) studies exemplify the fourth step of Wimsatt’s procedure—that they demonstrated the unreliability of simulations based on the breaking down of their robustness. These two papers cannot, however, be considered instances of robustness analysis, for at least two reasons. Let us consider first the target of van den Bosch and Ogiya (Reference van den Bosch and Ogiya2018). What the authors examined is not the trustworthiness of the simulations outcomes but that of the convergence criterion itself. Their main conclusion is that “most, if not all, disruption of substructure in N-body simulations is numerical in origin,” a conclusion that “questions whether the fact that subhalo mass functions appear to be converged down to 50–100 particles per subhalo implies that results are reliable” (4084). In other words, the authors showed that, even though convergence was reached, numerical artifacts had not been excluded and that the results backed up through convergence were still not reliable. What is at stake here is not whether the amount of substructure predicted through N-body simulations is a reliable prediction but whether the reason it is taken to be a reliable prediction is sound.

Still, one could insist that this study makes a negative use of robustness analysis, that is, that it uses robustness to show that predictions made about the number of satellite galaxies that should be observed in halos are not reliable, because their apparent robustness breaks down when the resolution is increased further. However, this reconstruction of their argument is misleading: what we see in their paper is not that the appearance of convergence of the results based on Power et al. (Reference Power, Navarro, Jenkins, Frenk, White, Springel, Stadel and Quinn2003) breaks down while ‘true’ convergence is reached for higher resolution. What we see instead is a first instance of convergence for a given range of values and a second one for another range of values. Asserting that the robustness of the first breaks down when increasing the resolution is assuming that we have more reasons to trust the second set of results than the first. However, on the basis of convergence only, we have equally good reasons to accept one result or the other—and, conversely, equally good reasons to reject them. Convergence alone will not tell us which result to trust. But here is where the problem becomes particularly sharp: in the astrophysical context at least, no other criterion is available to back up or supplement the results of robustness analysis. Simulations are used to determine what rival DM models tell us about the structure of the universe precisely because nonlinear effects related to star formation and gas dynamics make it impossible to determine the distribution of matter in DM halos analytically. Thus, there is no analytical solution against which the results of simulations can be tested. Furthermore, a comparison between observations of ‘real’ DM halos and simulations is of no help, for at least two reasons. While the similarity between the model and the real world is often used to validate the former, what is under scrutiny here is precisely the pertinence of such a comparison given that we do not know whether the simulated outcome reliably tracks the predictions of the CDM model or crucially depends on the parametrization process. Second, parametrization—or ‘calibration’, or ‘tuning’—is also the process through which one improves the agreement with real world data when other processes of validation are not feasible (Oberkampf Reference Oberkampf2014, 33). How could one then appeal to the agreement between observations and simulations when trying to circumscribe the impact of calibration on the simulation outcome? But in the absence of any other criterion to assess reliability, how could one conclude anything from the fact that two mutually exclusive but converged results can be found in two different regions of the parameter space?

Even more undoubtedly, robustness analysis has no role to play in Baushev et al.’s result. Like that of van den Bosch and Ogiya (Reference van den Bosch and Ogiya2018), the aim of Baushev et al. (Reference Baushev, Valle, Campusano, Escala, Muñoz and Palma2017) is not to test the reliability of the Navarro et al. (Reference Navarro, Frenk and White1997) density profile per se but to test whether its robustness says anything about its reliability. This goal does not involve any attempt to make the stability of the profile break down. On the contrary, their argument shows that its robustness will not break down because it results from numerical artifacts. Simulations converge because the artificial diffusion streams generated by numerical artifacts complement each other. The lesson to be drawn from this article is that robustness by no means says anything about the physical nature of the prediction made and even less about its reliability.