Warm Dark Matter in Simulations

: In recent years, warm dark matter models have been studied as a viable alternative to the cold dark matter models. The warm dark matter particle properties are expected to imprint distinct signatures on the structure formation at both large and small scales and there have been many attempts to study these properties with numerical simulations. In this paper, we review and update on warm dark matter simulation studies from the past two decades and their most signiﬁcant results: structure formation mechanisms, halos evolution, sizes and distribution, and internal structure properties. We discuss the theoretical assumptions and the limitations of the methods employed. In this context, several controversial claims are scrutinized in the attempt to clarify these confusing and sometimes even contradictory conclusions in the numerical simulation literature. We address the circumstances in which a promising keV dark matter candidate should be properly treated in the simulations.


Introduction
Bridging theory with observations, numerical simulations represent an important tool in studying and testing the hypothesized dark matter models, and in making predictions that can be then tested against astronomical observations. Several warm dark matter (WDM) models have been theorized and proposed as viable dark matter candidates, e.g., [1][2][3][4][5][6]. Given the challenges that the cold dark matter (CDM) model faces, these models, especially those with particles in the keV range, are gaining a lot of traction as feasible alternatives [7][8][9][10][11][12][13][14][15]. The complex properties of the warm dark matter candidates are expected to manifest at various scales-from large-scale structure to internal structure of galaxies-into detectable features that will allow for a more comprehensive picture of dark matter and a better understanding of the observations. Unfortunately, many times these properties are inaccurately described in simulations, and the results appear confusing and even contradictory. These results, their limitations notwithstanding, have been misinterpreted as strong constraints against warm dark matter in general, giving way to confusion and unfounded criticism.
In this article, we review the research done in warm dark matter simulations spanning over two decades and what it teaches us about the warm dark matter candidates. We examine the theoretical models that are depicted in simulations, the numerical methods employed together with their approximations, and finally, the results and their interpretations. We will see how the results from the simulations are not accurate in giving a strong constraint on the mass of the warm dark matter particle in the absence of a universal particle physics model [16,17]. They can act, however, as an important tool for describing the qualitative differences between several models at both large and small scales, and they can provide hints of where these qualitative differences should be further investigated [16,17]. Several confusing claims found in the literature will be discussed and clarified.
In Section 2, we give a short overview of the theoretical warm dark matter models. In Section 3, we examine how these warm dark matter models are represented in simulations and discuss the rough approximations that are made. In Section 4, we address the simulation analysis methods and the numerical difficulties encountered. Sections 5 and 6

Simulating Warm Dark Matter
To describe WDM particles in simulations in a way in which these simulations can be easily compared with the CDM counterparts, several approximations have been made. These approximations, reviewed below, begin with a 'cutoff' in the power spectrum and continue with an optional representation of thermal velocity dispersions completely ignoring other properties of the WDM candidates such as quantum pressure.

The First Cut
The cosmological simulations of warm dark matter collisionless particles, be they pure N-body or Hydrodynamical simulations, begin with a cutoff in the power spectrum. And as far as the difference between cold dark matter and warm dark matter is concerned, most of the simulations also end with that.The streaming velocity of the warm dark matter particles suppresses the matter power spectrum P(k) on scales smaller than their free-streaming scale; A rough estimation of the free-streaming scale for a particle X is given by [18]: where k FS is the free-streaming wavenumber and λ FS is the free-streaming length.
To implement the cutoff in the power spectrum the most common method uses the fitting formula suggested in [19] to compute the transfer function: where α, the scale of the break, is a function of the WDM parameters, while the index ν is fixed. Using a Boltzmann code simulation it has been found that ν = 1.12 is the best fit for k < 5 h Mpc −1 [20,21], and the following formula for α has been used: 11 Ω X 0.25 0.11 h 0.7 1.22 Mpc h −1 , with the cosmological parameter Ω X representing the particle's contribution to the cosmic density. An example of how the linear power spectrum cutoff dampens the small scales power for several warm dark matter masses is given by [22] in their Figure 2 with respect to the cold dark matter power spectrum for a fully thermalized particle; The smaller the mass of the particle is, the more the power is damped at small scales, which results in a suppression of smaller mass halos. The damping of the power depends, however, on the properties of the particle and thus on the model, therefore different models can have different expressions for the cutoff, e.g., [23] with the linear power spectrum shown for different masses in their Figure 3 and Figure 4, this time for a nonthermal sterile neutrino.

Cutting Deeper
As we have seen, even for the cutoff in the power spectrum the consensus breaks down, as there is a strong dependence on the model, but most of the simulations are indeed using the aforementioned recipe. The approaches are about to become even more arbitrary when it comes to addressing the velocity dispersion in simulations. The discussion is, however, of vital importance because even the power spectrum recipe is dependent on the velocity of the particle and the model for the mass-velocity correspondence that is assumed.

Velocity Dispersion in Simulations-Divergent Approaches
In [16,17], we provided a lengthy analysis of why the velocities are important in describing the warm dark matter particle. To summarize, the arguments diverge in two different approaches: simulations that include a velocity dispersion in agreement with the cutoff in power spectrum and simulations that do not.
For the inclusion, the arguments can be synthesized as follows: • The finite phase space argument-since for a warm dark matter particle the distribution in velocity space is important, e.g., [24] taking a strict zero thermal velocity is inconsistent with a finite phase space collisionless fluid. • 2-body relaxation problem for a large homogenous medium [25,26] shows that the relaxation time is proportional to the velocity of the particle and is zero for v = 0. • Heisenberg's inequality sets a maximum phase space density that translates to a minimum particle velocity dispersion. • The non-negligible influence of a 'negligible' velocity dispersion-while for some studies the velocity dispersion for a particle in the keV range may have a negligible effect, for other studies, especially at small scales, the effect becomes important not just in phase space and for the 'physical' correctness of the problem, but also in real space.
For the exclusion, the arguments are: • The numerical resolution is very limited and thus the sampling of the phase space density distribution is poor in both 'real' space and velocity space. • Thermal velocities for warm dark matter particles, in the keV range for example, are comparable (or smaller) to the bulk Zel'dovich velocities and thus, negligible.
• Furthermore, so introducing in the simulations a velocity that is much larger, due to the resolution limitations, to describe a contribution that is much smaller in reality, may introduce artificial artifacts. This argument that emphasizes the important effects of adding the velocities is usually followed by an argument that downgrades the real effect of the 'real' velocities.
With no middle ground solution, the implementation of velocities becomes a personal choice, influenced often by the limitations of the computing resources. In reality, while the cutoff in the power spectrum that corresponds to a certain velocity dispersion for a particle with a mass in a few keV range is enough to reproduce the large scale structure to a good approximation, at small scales the addition of the thermal velocity dispersion becomes significant. Indeed, for phase space density studies, for example, the difference becomes apparent as we will show in Section 6.2.

No Escaping Velocities
Several studies have been evaluating this mass-velocity correspondence, but in most of the simulations the adopted relation is the one given in [19]. In [16,17], we inspected those results in detail, emphasizing some implications of the assumptions made.
To briefly review the frequently cited expression for the velocity of the WDM particle as a function of mass and redshift, we commence with the expression in [19] Appendix A that recalls the abundance with respect to photons of a particle X that decouples while relativistic: where g dec is the number of relativistic species present at decoupling, and g X is the number of spin states of the particle. Connecting the particle mass density with the cosmological parameters Ω X ≡ ρx/ρc and h, and using the critical universal density ρ c ≡ 3H 2 /8πG, and the Hubble constant H ≡ 100 h km s −1 kpc −1 , they derive Assuming that the distribution function scales as the non-thermal distribution (exp(v/v 0 ) + 1) −1 , the velocity v 0 at redshift z (in [19] Equation (A3)) is given, Eliminating Ω X h 2 in this previous equation, for a m X = 1 keV particle (rounding to 2 significant digits), the velocity-mass correspondence formula is obtained Thus, for a 1 keV particle, we find g dec = 1000 (g X /1.5) 1/3 . This is, perhaps, the first striking difference with respect to other models-the high number of degrees of freedom at decoupling. Indeed, the value here-much higher than the value for the standard model ∼107 and even the supersymmetric extension ∼229 [27]-is also much higher than in other warm dark matter studies that were more 'conservative' in their choice for the number of degrees of freedom. Ref. [27] cites a value of 150 for a massive neutrino; ref. [24] cites a value of 100 for right-handed neutrinos that would decouple before the electroweak phase transition. This high value can be linked to higher-dimensional gauge group representations, as recently shown in [28].
In a similar manner to [19], the velocity-mass correspondence is derived in [29] assuming the following freeze-out dimensionless distribution where = mc|v| is the kinetic energy and T is the (comoving) decoupling temperature. At redshift zero, the velocity dispersion σ 0 is given by with the following correction for the distribution 5 ζ(5) where ζ(x) is Riemann's ζ−function. This gives an expression for the velocity citing a value of g = 2 or g = 4 for a −1/2 spin sterile neutrino [1,2,30].
With yet another set of assumptions based on number conservation and non-entropy production, and taking into account the quantum pressure while assuming a thermalization caused by an exchange potential, we derived in [16,17] a mass-thermal velocity correspondence that is independent on the cosmological parameters. Considering the full Fermi-Dirac and Bose-Einstein distribution, e.g., [31].
where p is the particle momentum, g the spin-degeneracy factor (of order 1 or 2), µ the chemical potential, = p 2 c 2 + m 2 c 4 − mc 2 the particle energy, m the particle mass, and T the particle temperature, the thermodynamical quantities for fermions and bosons at all regimes can be calculated accurately by evaluating numerically the relativistic Fermi-Dirac or Bose-Einstein integrals for particle density n, energy density e, and pressure P as functions of temperature T and chemical potential µ [31]. Connecting the conserved particle density with the expansion via the scale factor a = 1/(1 + z), for a constant particle entropy, the solution of this system for Fermi-Dirac particles in the non-relativistic regime is: for n 0 = 115 cm −3 , m = 1 keV, g = 1. The approximated solution for velocities in the non-relativistic regime is thus: for the temperature The full discussion and derivation can be found in [17] following the considerations in [31,32]. We do note that for the Maxwell-Boltzmann distribution that is usually considered in simulations, the velocity coefficient is found to be 0.20592, a value between the values for Fermi-Dirac and Bose-Einstein distributions, respectively, but the correction remains negligible.
For studies of neutrinos with much smaller masses, below 1 eV, assuming a homogenous, isotropic distribution with no chemical potential, although mentioning that the chemical potential is not zero necessarily, ref. [33] estimated a mass-velocity relation of: with the present neutrino temperature (in eV units) of Distinguishing between thermal relics (TR) and neutrinos produced through a Non-Resonant Production mechanism (NRP)-non-resonant oscillations with active neutrinos [1,3,34]-the following expressions for the velocities are derived in [35]: and with the parameter T TR (for TR) inferred from (T TR /T ν ) 3 = Ω WDM h 2 (94 eV/m). To summarize, the theoretical models for warm dark matter particle candidates are flexible in their initial assumptions, and are thus far from converging towards a 'universal' WDM model. We can already see how depending on the actual model assumed for the production of a warm dark matter particle candidate, the parameters considered, and the distribution function, the correspondence between the mass of the particle and its velocity can vary by an important factor.
In Table 1, for particles with masses of 0.2, 1 and 3.5 keV, respectively, we give the different values for velocities corresponding to some of the various studies we reviewed. While these differences in value may seem small at redshift zero, their effect is cumulative for the higher redshift when the simulations are started, as discussed in [16]. Even a factor of 2 difference in the velocity dispersion may be of important consequence for the structure formation, mass function and internal structure of halos. This difference would translate into a difference between let us say 1 keV and 1.6 keV within the model in [19] and would reflect a significant difference in the halo mass function and the density/phase space density profiles, as will be discussed in the following sections.
We note that all the values for the velocities presented in this table that do depend on the cosmological parameters have been obtained under the assumption that there is only one species of dark matter particle that contributes to the dark matter content. If the dark matter content is made up of more than one species, the correspondence found in [19], for example, does not hold and simply calibrating the Ω X parameter does not provide an accurate relation between the mass of the particle and the velocity dispersion. Overlooking this aspect in the literature resulted in a mixed dark matter mix-up, as shown in [16].
Moreover, while we referred mainly to particles with velocities from hundreds of eV to a few keV that do not decay, the warm dark matter picture can contain, from the point of view of particle physics, a wide variety of particles, with different properties and behaviors than the cases we considered. Indeed, heavy sterile neutrinos with rest masses in the few hundred MeV that decay nonthermally [36], or particles with masses less than 10 MeV that remain in equilibrium with neutrinos until becoming non-relativistic, reheat the neutrinos with respect to the electromagnetic plasma and therefore lead to extra energy density in the early universe [37], have been theorized. Light sterile neutrinos, which decouple shortly after the QCD phase transition with a highly non-thermal distribution function, have been proposed to explain, for example, the 3.55 keV line and the cores in dwarf galaxies, as shown in [38]. Table 1. The mass-velocity dispersion relation found in the literature for 0.2, 1 and 3.5 keV. The first column shows the value originally given in [19]-note that while in the table the velocity is calculated using the 3.571 value for the rms velocity factor, as usually found in the literature, the rms value should be 3.5971 instead, as shown in [16]; the second column shows the value obtained using g dec = 150 [27] in Equation (6); the third one, the value given by Equation (14) ( [16]); and the fourth and fifth, the value given by [35] for thermal and non-resonant produced neutrinos, respectively. All the values are given at redshift z = 0. Since the simulations are, as we will see in the following sections, extremely sensitive to the velocity parameter even when the actual velocity dispersion is not included, in the absence of a universal particle physics model for the warm dark matter candidate, the constraints on the mass of the particle from the simulations are naive. Qualitative differences with respect to the structure formation, evolution and internal structure of galactic halos are, however, a possible constraint on the velocity of the particle and a rough constraint on the mass. Unfortunately, as we will see when we look at the analyses of the simulation results, even when assuming the same model, the conclusions may be conflicting and confusing.

The Deepest Cut
One parameter that is completely ignored in warm dark matter simulations is the quantum pressure of the fermionic particles. An intrinsic property of the WDM fermions is the non-zero quantum pressure that when balanced by the gravitational pressure results, in dwarf galaxies, in cores of sizes that are in agreement with the observations [11,12]. This implies that while dwarf galaxies act as quantum degenerate objects supported by the quantum pressure, in larger galaxies the phase space density in the classical dilute limit determines the size of the cores. Considering a lower limit of ∼2 keV, the quantum description is able to predict and explain the observed cores in a broad variety of galaxies, regardless of their mass [11,12]. This result is largely ignored and the quantum effects in general are omitted from the simulations and their analysis, even when the WDM model is compared with other models that imply quantum pressure and for which the quantum pressure is indeed included, as we will see in Section 7.

Simulation Methods and Analysis
The technical numerical description of dark matter in simulations is done using collisionless particles rather than a collisional fluid. Several parallel codes are used to describe the evolution of the system, either with pure N-body methods or, when the baryons are added, with hydrodynamics [39,40]. The main difference between the methods used is in general the choice between a grid method that uses a regular cubic lattice and a glass method where the sign of the peculiar gravitational acceleration is reversed [41,42]. For the grid lattice, when using the same initial conditions, simulations ran with different codes [39,40] show very little difference at very small scales [17]. While these differences are not significant for the structure formation analysis, they become significant for precision cosmology studies [43]. The initial conditions-extremely important for non-linear processes [44]-can be generated using one of the many elaborate codes designed for this purpose [45][46][47][48]. While discussing these methods is beyond the purpose of this work, we will address some significant differences between grid and glass-like methods in regards to the so-called fragmentation problem that is prevalent in all warm dark matter simulations.

The Fragmentation of WDM Filaments
The N-body method is proven successful in describing the behavior of CDM particles as far as the convergence tests are concerned [49][50][51][52]. Several studies, however, pointed out the effects of discreteness errors on the numerical accuracy at small scales, in particular, e.g., [42,[53][54][55][56][57][58]. These errors become more problematic for the WDM case [19,59]. Ever since the first simulations, e.g., [60][61][62][63], the fragmentation of filaments into small mass halos was observed, and for a long time this was placed on the discreteness effects that are influencing the growth of nonlinear structures below the streaming scale, although some studies argued its cause to be physical [19,64].
Indeed, N-body simulations of WDM particles showed filaments fragmenting into small similar mass halos equally spaced inside the filaments. Since most of these simulations have been performed with the grid method, which employs a cubic regular lattice, the numerical artifacts were explained by the preferred directions and the coherence of the lattice. The glass-like method, on the other hand, has the sign of the peculiar gravitational acceleration reversed so that each particle is repelled by all the others resulting in the total force to vanish on each particle when the system reaches a quasi-equilibrium state. Since there are no preferred directions and no memory of the initial grid, one would expect no artificial fragmentation in this case. The simulations painted a different picture, though, and while [65,66] found that in the glass simulations the number of spurious halos is reduced compared to the grid simulations, ref. [42] found regularly spaced artificial halos with a frequency similar to that present in the grid simulations. A further analysis showed that the fragmentation depends on and scales with the resolution, with fewer but more massive 'spurious' halos in the low-resolution runs [42]. While overall the fragmentation was more pronounced in the grid run, it still occurred in all glass simulations.
A similar resolution test has been performed for the grid method in [17] in a region where a halo forms top-down at the intersection of large filaments. The higher-resolution filament shows a gradient in density with fragmentation occurring in the highest density region in the central part of the filament, as opposed to the whole filament as in the low-resolution case. In the high-resolution run, the number of spurious halos increasesscaling with N 3 -while their mass decreases-scaling with N −1 -confirming thus the results in [42]. The same high-resolution simulation has been performed with an increased softening length and as expected, when the softening is comparable to or higher than the fragmentation scale, no spurious halos form [17].
Further studies showed [67] that the fragmentation due to large anisotropic force errors can be corrected using a more accurate force calculation that results in a reduced number of 'spurious' halos. To solve the same problem, ref. [68] proposed a modified force softening criterion, which minimizes the spurious two-body effects, while maintaining high force accuracy in collapsed regions. Other studies showed that infinite self-gravitating cylinders are unstable to fragmentation [69,70] and that for finite cosmological filaments the fragmentation may occur in some cases regardless of the resolution and the softening parameters [68]-so the question whether the fragmentation is physical in some filaments, still remains.
In practice, the fragmentation problem is overcome by eliminating the 'spurious' halos altogether from the analysis, since the mass scale at which these halos form can be computed [71]. The caveat with this approach is that in WDM simulations, some of these 'spurious' halos merge together into larger halos, well above the fragmentation mass scale [17], 'contaminating' thus the mass function at redshift zero. Furthermore, using popular 'friends-of-friends' (FoF) methods to identify halos in CDM simulations and compute the mass function proves to be insufficient in analyzing halos in WDM simulations [72]. Extreme caution is thus required when analyzing a halo simply by looking at the redshift zero simulation; tracing back its evolution in the simulations is therefore advised. Moreover, eliminating from the analysis a significant number of particles does question the 'physicality' of the simulation with respect to its mass content and so the jury is still out on what the best approach should be. How the baryonic processes influence and/or are influenced by the fragmentation, is another important question that has not been extensively studied in simulations.

Structure Formation in WDM Simulations
When comparing the structure formation in CDM simulations versus WDM ones, many times the analysis is done using the same tools, methods and assumptions. As we will see, however, the structure formation in WDM is more complex, with several phases depending on the mass of the particle on one hand, and the morphology of the region analyzed on the other. These qualitative differences can escape a minimalistic analysis, where only several simulations outputs are compared, but are obvious in a consistent sequential analysis employing a large number of simulation outputs and tracing the particles back to earlier redshifts. This is best conveyed by movies made of a large number of simulation outputs that are showing the structure formation and evolution with redshift, but since this method is expensive as far as the computational resources are concerned, this type of analysis is rarely used.

Top-Down or Bottom-Up?
Traditionally, the top-down structure formation refers to a structure formation scenario in which the biggest structures at redshift zero are also the earliest formed structures that then continue to grow through the accretion of matter from the filaments, causing the filaments to break into smaller halos. Conversely, in the bottom-up structure formation scenario, small halos form first from fluctuations at the smallest scales and via mergers, they form larger halos and clusters that show a higher concentration of particles at redshift zero. While technically both these mechanisms of structure formations are 'hierarchical'-with the hierarchy reversed in one model with respect to the other-in the literature, the term 'hierarchical' has been used to describe the bottom-up mechanism.
In earlier simulation studies, for dark matter particles with a mass around or below 100 eV, the structure formation was associated with the top-down scenario-the biggest structures forming first via gravitational collapse in the highest density regions, growing then through accretion. We note this mechanism differs from the mechanism originally proposed by [73] to describe how a galaxy forms through the collapse of a cloud of gas-the monolithic collapse mechanism [74,75]. While in both these cases structures are formed by collapse, in the monolithic collapse case the halo forms from a cloud of matter, while in the WDM case, it forms at the intersection of three-dimensional filaments where density reaches a peak. Later simulation studies, which employed colder WDM particles in the few keV regimes, claimed that the structure formation is rather 'hierarchical', probably because the redshift zero simulation outputs did not look much different than the CDM ones, with the exception, perhaps, of a lower number of satellites in the WDM case.
More comprehensive studies that performed, analyzed and compared simulations of particles with velocities in a broader range-from hot to cold dark matter-reported a much more complex structure formation mechanism for the WDM particles that depending on the velocity of the particle, go from top-down to bottom-up, with a lot of overlap in the middle, where different regions experience different structure formation scenarios depending on the morphology [16,17].
While these results are from pure N-body simulations, adding gas and baryonic processes may introduce, as we will see in the next section, yet another degree of complexity with respect to the structure formation in WDM scenarios.

A Hybrid Structure Formation Mechanism
When carefully analyzing WDM simulations in a broader range of velocities, particularly between what would correspond in the model cited in [19] to masses in the 50 eV-2 keV range, various structure formation scenarios have been reported [16,17]. Focusing on a transient intermediate regime of particles with dispersion velocities around 0.3 km/s (∼200 eV), where the complexity of the structure formation mechanism becomes obvious, the various layers of structure formation have been explored. We note that while otherwise ruled out by many observational constraints mainly because the structure formation is very much delayed in simulations with such light particles, it is this delay and the slow evolution of the system that allowed for a better examination of the structure formation differences. By extrapolation, these structure formation trends would be present even in keV scenarios, but at smaller scales and with a more rapid evolution [16,17], thus escaping superficial analysis.
The different phases of structure formation are summarized below and illustrated in Figure 1. A similar figure with the zoom-in evolution of these regions at different redshifts can be found in Appendix C of [16]. Several movies following the evolution of these simulations can be found on Youtube 1 .

•
The collapse starts with sheets of matter collapsing into filaments, which then collapse into halos. • The first halos are indeed formed top-down in the high density regions found at the intersection of well-contoured filaments. These halos begin accreting matter from the disrupted filaments. • Depending on the morphology of the region, some of these halos can become very massive very fast just by accretion, and they can survive without major mergers until redshift zero. • Later on, in medium density regions, depending on the spatial distribution of filaments in that region, halos merge into bigger ones, signaling the beginning of a bottom-up growth scenario. • In less dense regions usually situated in voids, the collapse is even more delayed, with filaments being formed and collapsing very late. This favors the merger-free survival of a halo formed top-down all the way to redshift zero. • Several of the early formed massive halos-depending once again on the region's morphology-are violently merging together at later times forming large clusters.
Therefore, we see that most of the first-formed halos, which are indeed formed through top-down gravitational collapse, are then merging into more massive ones. The statement that the most massive halos at redshift zero are the earliest ones that form is not valid for all the regions in the simulations, meaning the structure formation is not at all scales a top-down mechanism in the consecrated sense. The late-stage 'hierarchical' build-up that follows the formation of the first halos does not initiate from the small-scale fluctuations as would be the case in the classical bottom-up CDM scenario, but it is achieved through mergers. In the later stages, while in medium to high-density regions the structure evolves through mergers, the top-down trend still continues in low-density regions-voids, and thus we have the two distinct structure formation mechanisms operating at the same time.
The structure formation mechanism observed in WDM simulations is therefore a more complex and complicated hybrid mechanism, where long-range and short-range effects are both present.  [17]. The different regions corresponding to different mechanisms are encircled: in green-a region where a large halo forms top-down, in black-a region where several halos are approaching merging, in blue-a halo forming late in an isolated region.

The Baryon Component
When it comes to the evolution of structure in the universe, it is usually believed that dark matter dominates the large-scale structure formation, while baryonic processes come to dominate the small-scale internal structure of galaxies at their core. From the tapestry of phenomena that involve baryonic processes, predicting with high confidence what effect each individual process has on the structure formation at observable scales is problematic to begin with, because many of these processes take place at the same time in many of the cases observed. Simulations represent an important tool for testing the baryonic effects and in recent years there have been considerable developments related to the implementation of baryonic physics in cosmological simulations, but they are far from providing high-accuracy descriptions and/or predictions about the myriad of physical processes that may occur. Moreover, most of these studies have been done solely in the CDM structure formation scenario, with only a few attempts in WDM scenarios.
There is, however, one study that stands out, in which using unparalleled highresolution simulations of dark matter plus gas, the very early formation of stars has been studied in both WDM and CDM scenarios [76]. For the WDM case, the simulations employed a 3 keV gravitino-type particle and starting from a z ∼ 23 redshift, the gas behavior inside a 3 kpc filament was analyzed. While in most simulation papers in the literature a 3 keV particle is indistinguishable from a CDM one, the simulations in [76] showed a remarkable difference: in WDM stars form before the halo. Indeed, in the WDM case, the filament is preserved for a longer time due to the free streaming of the particles, allowing thus enough time for the gas to heat as it becomes compressed. In this context, the stars form inside the filament before the filament collapses and forms halos. Due to the formation of molecular hydrogen, the gas then cools, cooling also the center of the filament. This cooling can act back as a delaying factor in the collapse of the filament. This is very different from the CDM simulations in which the filament fragments first into small halos that then accrete the gas and reach the optimal conditions for star formation much later in the high-density regions at their centers. This is not a surprising result since the structure forms differently in WDM with respect to CDM, and it has been shown that even in 1.5 keV particle simulations stars form inside the filaments up to redshift z ∼ 2, resulting in stringy 'chain' galaxies that may be confirmed by observations [72].
The difficulty in observing and analyzing such effects lies in the very high resolution these simulations require, M dm = 272.6M and M gas = 41.9M in a commoving box with a Lagrangian radius of 600 kpc [76]. Other studies that do not benefit from the same high resolution claim that on the contrary, the structure formation is simply delayed in WDM [77]. In [16], it has been argued that the top-down trend is present also for particles corresponding to a few keV (in the estimation from [19]), but it is more difficult to see it with the commonly used resolution, because it is rapidly overtaken by the hierarchical growth and can be missed completely without a rigorous analysis. Therefore, using for WDM particles the same parameters and the same resolution that is commonly used for CDM, especially when including baryonic processes, does not justify strong claims, constraints and conclusions.

Mergers and Satellites
The gravitational interaction of two or more halos encountering each other is generally referred to as merging. Another distinction is required here for the case in which a large galaxy 'accretes' smaller halos-a minor merger event-when the orbits of small satellite halos around a large host galaxy are reduced due to dynamical friction such that the halos 'fall' towards the center of the more massive galactic halo. When two halos or galaxies of similar masses are interacting, their merging results in a new galaxy with a structure different than that of the preexisting objects. These major mergers are considered responsible for the formation of spheroidal galaxies. The structure of the merger remnant depends, however, on the structure of the major mergers. The merging of purely stellar disk halos will produce a spheroidal merger remnant with a central phase space density lower than that of the ellipticals observed [78].
Numerical simulations showed that mergers with a mass ratio M 2 /M 1 ∼ 0.25 destroy galactic disks and form spheroidal remnants, while mergers with a lower mass ratio (very minor mergers) may form a thickened disk system. In the pure N-body cold dark matter simulations, since mergers are ubiquitous, disk galaxies have not been found. The presence of baryons, however, may play an important role, and in some exceptional cases, where the major mergers are having a gas fraction larger than ∼50%, the formation of a disk may occur in the remnant galaxy [79][80][81].
In N-body simulations, merger trees can be extracted directly using several methods, e.g., [82][83][84]. In the case of substructures, however, a more complex problem arises when three-body encounters can result in the ejection of halos in the N-body simulation [85]. In WDM simulations, in order to compute an accurate mass function, one needs to first compute the merger trees, since the small halos formed due to the fragmentation may sometimes merge into larger halos that have a mass above the fragmentation scale, as shown in Section 4.1. Without a merger tree analysis, these halos are difficult to spot. This is an aspect usually overlooked in the halo-mass studies in the literature.

Halo Mass Distribution
Besides the significant differences in the structure formation that we mentioned in the previous section, the statistical distribution of halos at different times and with different masses can be yet another way to distinguish between CDM and WDM models and therefore, to constrain them with observations.
Based on the assumption that halos are correlated with peaks in the Gaussian random density field of dark matter in the early universe [86] and adding the halos that form in underdense regions, a rough estimation on the number of halos with a mass in a range M to M + δM per unit volume, δM(dn/dM) is given by [87][88][89] dn dM (M, t) = 2 where ρ 0 is the mean density in the universe, σ(M) is the fractional root variance in the density field smoothed using a filter that contains a mass M, on average, and δ c (t) is the critical overdensity for spherical collapse at time t [90]. This simple formula, however, does not accurately apply to N-body simulations [91] and several other methods have been proposed for fitting the mass function [92][93][94]. Dependent on simulation parameters, the estimation of the mass function is given by, with where A, a, b and c are parameters inferred from N-body simulations results. The mass variance inferred from the power spectrum of density fluctuations is where P(k) is the initial power spectrum, T(k) is the dark matter transfer function and W M (k) is the Fourier transform [95]. These fitting parameters depend on the simulation redshift and are therefore unable to predict a universal mass function and consequentially, a universal description of the halo mass distribution. In the case of warm dark matter, the mass function is also sensitive to the structure formation mechanism and it may become 'contaminated' by the presence of spurious halos formed by artificial fragmentation.
In a recent attempt to find a general parametrization of the WDM halo mass function, ref. [96] come to similar conclusions and caution that strong constraints will require properly tailored simulations and careful considerations of the halo mass estimations and methods of analysis.

Internal Structure of WDM Halos
In CDM simulations, halos form from the smallest scales fluctuations, merging into larger halos as the simulation evolves. In WDM, pure N-body simulations show a more complex hybrid mechanism [16], with the first halos forming top-down at the intersection of large filaments and evolving then in different ways, either by accretion or by merging, depending locally on the morphology of the region where they formed and globally on the velocity of the simulated particle. The influence of the warm dark matter properties is expected to extend from these structure formation effects to the internal structure of galaxies in several ways that will be discussed in this section.

Shells and Caustics
The dark matter shells or caustics are formed first through the gravitational collapse of dark matter and then enhanced by halo mergers and the tidal disruption of satellites. The formation of halos, as described by the Jeans-Vlasov-Poisson equation, results in the formation of shells or caustics where the dark matter streams meet at the surface of high density regions [97][98][99][100][101]. High-density caustics can also be formed by the tidal tails-the satellites falling inside the potential well of the host halos. Caustics and shells of stars formed during galaxy mergers have been observed in so-called shell galaxies [102][103][104][105].
As these shells and caustics cannot be destroyed once formed, they become permanent structures in both real and phase space. The density in shells and caustics exceeds the density of the neighboring regions by a factor that depends on the age and history of that particular galaxy. Their study is therefore of high importance for establishing the mechanism controlling the structure formation and evolution, and hence, the nature of dark matter.
In the CDM case caustics are being wrapped inside earlier generations of the merging history during hierarchical evolution, and are thus not visible even in high-resolution simulations. However, using CDM simulations, it was shown in [106] that accretion mechanisms of stars and dark matter clumps, and the disruption of the latter, can produce concentering shells that resemble those observed in NGC 7600.
In WDM high-resolution simulations shells and caustics are largely present and highly visible, as shown in [16]. Using a high-resolution refinement of a WDM halo that forms top-down at the intersection of filaments and then starts accreting matter from those filaments, reaching 18 million particles in its r 200 radius at redshift zero, [16] clearly show the formation of shells and caustics and the coherent infall of material from the smooth surrounding regions and the filaments. The tidal tails of small satellites that are accreted into the central region are also visible. Given the more 'extreme' choice for velocity in this case, the size of the halo is almost three time larger than a simulated CDM halo of the same mass and thus less concentrated; it does not achieve virialization by redshift zero and it is not spherically symmetric. It is, however, a good medium for displaying the details of forming shells and caustics.
The presence of caustic and shells has been signaled in other studies, e.g., [71,107,108] that use both grid and glass simulation methods, but without including thermal velocities.
Since the presence of shells is a direct effect of the structure formation mechanism it is natural they would occur in simulations with or without the thermal component, but their structure is different, with thicker, more apparent caustics in simulations with thermal velocities [17].
The magnitude of radial velocities in the high resolution region where the halo forms is very different than what one expects in a spherically symmetric CDM halo, strongly correlated in orientation with the shells rather than with the center of the halo, as shown in Figure 2 for the halo analyzed in [16,17].
Shells and caustics become important features for detection since the higher density of a caustic on the line of sight can translate into a high peak in the detection signal, such as the one reported from the XMM and Chandra data-the 3.55 keV emission line [109,110]. An increase in the signal can also appear in the detection alongside a filament since the density in its inner central region is higher than in the outer part. Using an analytical model for the secondary infall mechanism [111,112] with a finite velocity dispersion [113], and assuming spherical symmetry and smooth accretion, it has been shown in [114] that assuming the weakly interacting species strongly annihilate in such dense regions, a boost in the signal of the antiprotons and positrons produced is expected. The annihilation signal is not considerably enhanced in the case of cold dark matter caustics [115].

Density Profiles
One of the two main initial challenges of the CDM model was the so-called 'cusp-core' problem: the density profiles inferred from the observed galaxy rotation curves show a constant core, e.g., [116][117][118][119] rather than the cuspy density profiles of the simulated CDM halos [50,[120][121][122][123]. Given the properties of the WDM particles that result in a cored profilethe maximum phase space density and the velocity dispersion-WDM was sought out as a viable alternative.
If we consider the simple spherical collapse model, the overdensity of a collapsed halo with respect to the background density is approximately ∼200, depending slightly on the cosmology, and it corresponds roughly to the virialized region of the halo. The virialized radius is thus defined: In CDM simulations, the virial theorem is obeyed in most cases within this radius and can be fitted either by the NFW profile [124,125] where r s is a characteristic scale radius and ρ s is the density at r = r s , or by the Einasto profile [126], as shown in [127][128][129] where r −2 is a characteristic radius at which the logarithmic slope of the density profile is −2 and α is a parameter controlling the variation of the logarithmic slope with the radius. Using the fitting formula provided in [130], we have where ν(M, z) ≡ ρ crit /σ(M, z) is a dimensionless 'peak-height' parameter defined as the ratio of the linear density threshold for collapse at redshift z within spheres of mean enclosed mass M.
Since, as we saw in the previous section, the WDM halos form differently than the CDM ones, the spherical collapse does not accurately describe the WDM halo, but for the purposes of comparing and reviewing the results in the literature we will work within this framework.  [16,17].
Under the assumption that the maximum coarse-grained phase space density value is conserved for the WDM particle, the effects on the halo's internal structure can be estimated-the velocity dispersion becomes constant and isotropic in the halo center, the density becomes a constant, and a core radius can be thus derived.
From the maximum phase space density of a homogeneous neutrino background, we can derive the Tremaine-Gunn limit [131] on the mass of a neutrino, assuming once again that neutrinos form bound structures where central regions can be well-approximated by isothermal spheres. For a Maxwellian velocity distribution, the maximum phase space density is where ρ 0 is the central density and σ is the one-dimensional velocity dispersion. The core radius is then given by which gives a minimum mass for neutrinos of: For a halo in the spherical collapse model, the maximum phase space density is given by [132,133] and it can be converted into a minimum 'core' size following [134]: where σ halo is the velocity dispersion of the simulated dark matter halo. However, this pseudo phase space density Q ∼ ρ/σ 3 [134] that is often cited in simulation studies, is not a proper coarse-grained average of the fine-grained phase space, overestimating, in fact, the maximum phase space density by an order of magnitude, as shown in [135].
In reality, at the final epoch, the maximum coarse-grained phase space density could also be lowered significantly due to mixing. Moreover, as velocities may be anisotropic, the density may not necessary approach a constant value towards the center of the haloincrease or decrease in values towards the halo center will result in a rise or fall in the density slopes.

Cores in Simulated WDM Halos
Several studies in the literature focused on analyzing in simulations the density and phase space density profiles for dark matter halos of different masses in regimes from tens of eV to a few keV. In simulations with particles hotter than 50 eV (in the mass-velocity correspondence from [19]) halos do not form by redshift zero. Studying the Tremaine-Gunn limit, several simulations discussed in [136][137][138] have been performed with particles of different masses, using the cutoff in the power spectrum described in [19]. For a few keV, the core size radius from the analytical estimation in [134] is already too small when compared with cores observed in similar mass galaxies and also very close to the resolution limit, which is roughly 0.5% of the virial radius [55]. Therefore, to study the Tremaine-Gunn limit in a regime that is better represented by the resolution of the simulation, a choice was made to keep the 'colder' power spectrum cutoff, while artificially increasing the thermal velocity of the particles in much hotter regimes [136].

Fitting Profiles
In the scenarios mentioned above, after re-simulating with higher resolution particles halos with masses 7 × 10 11 M -7 × 10 12 M , the density profiles and phase space density profiles have been studied [136].
When fitting the density profiles with an α, β, γ law profile using a Monte-Carlo Markov Chain (MCMC) algorithm, the transition to the central density core was so sharp that the α, β, γ turned out to be a poor fit to the simulated profile. A different parametrization [51] where the density profile is linear down to a scale R λ beyond which it approaches the central maximum density ρ 0 as r → 0 gives a better fit, as shown in [136]. This fitting function is extremely flexible making it possible to reproduce both cuspy profiles such as the ones predicted by the CDM model and highly cored profiles. For the CDM halos it was shown [139] that for over two and a half decades in radius, the phase space density profiles follow a power law ρ/σ 3 ≡ r α with α ∼ −1.9. For larger radii outside the core, this fit holds also for the WDM halos. A 'by-eye' fit shows a constant density core where the density profile breaks from the power law, even though in the analytical approach from the spherical collapse model [134] the core is defined as the radius at which the density drops by a factor of 2. When compared to the pseudo phase space density profiles, it was seen that the turn of the density profile towards a constant value corresponds indeed to the turn in the phase space density towards a constant value.
Performing similar simulation studies, in [135] it has been found that the density profiles are best fitted within the core radius by a pseudo-isothermal profile with a core [140,141] where ρ 0 is the core density and r c the core radius. While the maximum coarse-grained phase space density of the simulated haloes was found to be close to the theoretical finegrained upper limit, the maximum pseudo phase space density was significantly larger, by up to an order of magnitude, as mentioned previously. A similar analysis has been performed with a much larger high-resolution halo for a 200 eV particle, both without the velocities and with velocities corresponding to the same cutoff [17,137,138]. It was found that while in both cases the density profile deviates from the cuspy CDM profile, the profile is shallower in the simulation with velocities, with a lower central density, resulting in a larger core. The same applies to the pseudo phase space density, which in the presence of velocities showed an almost flat behavior towards the center, as expected. It is safe to say at this point that while the cutoff in the power spectrum alone can soften the abrupt cuspy behavior of pure CDM halos, the constant-core shape of the profiles is a consequence of the maximum phase space density given by the velocity dispersion, and not of the suppression of power at small scales.
Given the results of all these studies, one conclusion was the same, that the phase space density of WDM particles cannot be solely responsible for the large cores observed in galaxies, mainly because the particle that would form such cores, both in simulations and the analytical predictions, is too hot to have formed the galaxy in the first place-a Catch 22 [136]. Unfortunately the conclusions drawn from these results have been misleadingly used as either a definitive proof against WDM or as evidence supporting the idea that the observed cores are ultimately the effect of baryonic processes. We will discuss these misconceptions in Section 8.

'Real Cores'
Several sets of observations have been pointing to the existence of large cores. Using the stellar velocity dispersion profile measured in the Fornax dwarf spheroidal galaxy, it has been shown [133] that the distribution of the globular clusters inside the galaxy is consistent with a large core of ∼1.5 kpc, and thus inconsistent with the expectations from the CDM models. Dwarf spheroidals are indeed an important testing bed for dark matter models-in the Milky Way, dwarf spheroidals with masses between 10 5 and 10 7 M seem to be dominated by dark matter at all radii [142]. However, the few hundred pc cores in Fornax, Sculptor and Ursa Minor cannot be explained by CDM alone.
To reconcile this discrepancy-and possibly to save the CDM model-several models have been proposed using dynamical considerations. In several studies, it was argued that that the transfer of energy and angular momentum from massive infalling objects, such as black holes, can explain the central density in dwarfs [143][144][145][146]. The dynamical coupling of dark matter to energetic baryonic outflows has also been explored as a possible solution [147][148][149][150]. These assumptions, however, were not able to explain or reproduce other properties of dwarf spheroidals, such as the luminosity functions, the star formation history, etc. [151][152][153].
Under the assumption that baryonic processes have a dominant effect in the inner regions of all galaxies, even though in dwarf galaxies the baryonic content is poor compared to the dark matter content, several baryonic processes have been incorporated in hydrodynamical simulations in the attempt to reproduce the observed cores-for example: • Molecular hydrogen cooling of the gas is crucial for the formation of stars and galaxies [154,155]. • Gas radiative cooling causes the gas to fall towards the center of the dark matter halo due to the loss of pressure support, making the center denser [156,157]. • Star formation feedback may explain the cores in dwarf galaxies [158], and is also important for the formation of late-type galaxies, e.g., [159][160][161][162]. • Supernovae feedback-repetitive outflows from supernova feedback-may induce the formation of a core overcoming the adiabatic contraction of the halo [163]; or it may not [164]. • AGN feedback, on the other hand, heats the gas with a much higher energy than the one available from supernova feedback [165].
Many of these studies, however, are not without their shortcomings and do not provide a comprehensive solution to the problem, i.e., the ones that claim to solve the cores fail to solve the missing satellite problem (and vice versa) or they give extreme structure formation predictions, as reviewed and discussed in [166]. Some implementations give contradictory results, e.g., when including repeated baryonic outflows large cored galaxies are not present in some simulations [164] but seem to be present in others [163]. Furthermore, since in reality many physical processes are happening on scales below the resolution limit of galaxy simulations, the implementation of these processes is done 'by hand' (using semi-analytic approaches) [167] and thus the results depend strongly on the fine-tuning of the parameters and the resolution, e.g., [168,169].
Coming back to the dark matter effects in the dark matter dominated galaxies, it has been shown [11,12] that in fact by properly modeling the quantum pressure of fermionic particles [170,171], which has been completely ignored in WDM simulations, the expected cores can be reproduced [11,12]. Indeed applying the Thomas-Fermi semiclassical approach to fermionic WDM for a particle in the 1-2 keV range [11,12], the mass halo-radius, phase space density and velocity dispersion are found to be fully consistent with observations of compact dwarf galaxies. This suggests that dwarf galaxies are "natural quantum objects for WDM" [11,12] and hints that dark matter is not as 'simple' as was represented in simulations.
The quantum effects have also been studied in relation to the cores in other scenarios, i.e., for a fully degenerate Fermi gas [172] found a good fit to the profiles of eight dwarf spheroidal galaxies of the Milky Way for a ∼200 eV particle, claiming that such low mass overcomes the constraints in a DM non-thermally produced via inflation decay model and [173] discussed a 'Flooded Dark Matter' model in which the mass limit for DM is relaxed by a late-entropy release and where assuming DM to be a quasi-degenerate Fermi gas surrounded by a thermal envelope may explain the cores in dwarf galaxies.

Interpreting the Results
When properly interpreting these results, one needs to keep in mind that the WDM simulations are first based on a set of unproven assumptions. First, a particular model of particle production is assumed and that gives a certain mass-velocity correspondence. Second, even within a particular model, corners are cut in the proper physical description of the warm dark matter candidate with the velocity implementation and the quantum pressure of the particle often ignored. It is within this very crude approximation that certain conclusions have been drawn-that in order to produce significant cores in the kpc scale in dwarf galaxies that are in agreement with observations [153,174,175] one would need primordial velocities corresponding to a particle mass below 0.1 keV, which in turn is excluded by large-scale independent constraints (Lyman-alpha and lensing). This would stand if and only if the dark matter velocities would be the sole cause responsible for the size of the cores observed in galaxies. It has been shown, however, that the Fermi quantum pressure of the WDM particle has a significant contribution to the size of the core. Several conclusions from these results are summarized below: • The finite initial fine grained PSD sets a maximum coarse grained PSD, resulting in PSD profiles of WDM haloes similar to CDM halos in the outer regions, but which in inner regions turn over to a constant value set by the initial conditions. • The pseudo phase space density is an overestimation of the six-dimensional phase space density and the coarse grained phase space density. • The turnover in PSD results in a constant density core with a characteristic size in agreement with the simplest expectations. • All WDM simulated halos show this turn in the PSD towards a constant value and have cored rather than cusped profiles. • While the cutoff in the power spectrum alone results in a shallower PSD and density profile when compared to the CDM counterpart, it is the velocity dispersion that results in a constant PSD and a significant core in the density profile. • For a WDM candidate with high enough velocities as to produce the observed large cores, the free streaming would erase all perturbations on that scale so that the halos would not be able to form in the first place-a catch 22-in the scenario in which velocities alone would be responsible for the production of the core. • This is not the case, however, as simulations ignore the Fermi pressure of the dark matter particles that together with the velocity dispersion was shown to explain the large cores observed in galaxies for particles in the keV regime [11,12].

Double Standards
In previous sections, we discussed how, depending on the actual model of particle production, the WDM properties vary, in particular the mass-velocity dispersion correspondence. Since the simulations can only constrain the velocity of the particle, in the absence of a universal mass-velocity correspondence, a constraint on the mass can only be placed within the particular model that has been assumed and engaged in the initial conditions. Even within that particular model, simulations that do not properly solve the physics at the scales where the inherent WDM properties are important and/or simulations that do not have enough resolution to explore the effects at these scales can only provide a very crude estimation of the effects of the WDM free streaming and thus a very loose constraint on the velocity. Many studies have claimed strong constraints, or even ruling out WDM, regardless.
One example is the constraint from the Lyman-alpha forest. First, Lyman-alpha constraints are based on the temperature of the intergalactic medium (IGM)-the hydrogen absorption lines are used to re-construct the matter power spectrum, under the assumption that the hydrogen distribution follows the matter distribution. This distribution is compared with a 'mock' distribution from dark matter simulations. Several such studies report constraints on the mass of the particle, with [176] reviewing previous methods and studies and posing a 3.3 keV constraint at 2 sigma, a 'stronger' bound than in other reports [177,178]. This constraint was confused in the literature with a definite argument against warm dark matter as a plausible dark matter model. The justification was that a warm dark matter model with such 'cold' velocities does not help with the other issues that were raised against CDM, such as the internal structure of halos, and in [179], originally titled "Warm Dark Matter: The End is Nigh" the failure of the WDM thermal relic model is inferred. Firstly, the hydrodynamical simulations that the data is compared against do not have high enough resolution to accurately describe the WDM effects, such as the ones reported in [76]. Secondly, the study completely ignores the quantum pressure that does 'help', in that it increases the core size, making WDM quite distinguishable from CDM. Finally, as an exercise in logical reasoning, even if we overlook these shortcomings and omissions and we accept the result as binding, from two models that produce similar results there is no logical reason to disfavor one because it does not solve the problems of the other.
Moreover, a closer look at the Lyman-alpha forest constraints shows that there are several astrophysical mechanisms that influence the IGM and may affect the shape of the flux power spectrum, as reviewed in [5]. Different considerations and assumptions that may change the interpretation of the data are also discussed in [180] when comparing the Lymanalpha forest with simulations that assume resonantly produced sterile neutrinos [2,181], recognizing the need for a higher resolution set of simulations and concluding that in the analyzed case, no model is strongly favored.
In the simulation community the trend remains, however, to ignore the physical properties of warm dark matter and a proper attempt of representing them in simulations. The most striking example of applying double standards is, perhaps, when comparing warm dark matter with fuzzy dark matter (FDM) [182]. The study, besides ignoring a significant number of results in the literature, omits completely the quantum pressure of the WDM particle, but includes it for FDM. Furthermore, when 'comparing' FDM (with quantum pressure) with WDM (without its quantum pressure) it is specifically the quantum pressure that is cited as an important factor in preventing the fragmentation and producing a flat core in FDM, claimed as unique to the model [182], even though it was previously shown that the WDM quantum pressure has a significant contribution to the cores [11,12]. WDM is once again stripped of one of its inherent properties and blatantly discriminated against.

How to Simulate Warm Dark Matter
Warm dark matter properties, even in the most 'simple' particle physics models, have a more complex influence on the structure's formation, evolution, and its final properties, than the CDM particle properties. To properly simulate these properties, one needs to first integrate them in the simulation by adding a better description of the physics. Simulating WDM by a crude approximation in reference to the CDM can lead to erroneous claims and conclusions.
Once the physical properties are introduced, they induce behaviors that are qualitatively different from the CDM counterparts at large scales, as well as at small scales. To properly analyze these effects, one needs a high resolution, much higher than in a CDM simulation. Because of the free streaming, the particles travel a larger distance in the WDM case, thus a larger simulation region will need to be refined when using high-resolution refinements. Attempting to compare a CDM halo formed in a 5 Mpc box, for example, with a WDM counterpart by simulating just the 5 Mpc WDM box would not allow for the structure formation to be properly resolved and evolved in the WDM case. Since the structure formation happens in phases and is sensitive to both the simulation parameters and the rapidly changing environmental conditions, several assumptions made in CDM simulations do not apply to WDM. For example, we saw that introducing gas and other physical processes in a pure N-body simulation, produces significant effects in the WDM case, starting from a higher redshift, even before the filament collapse. The assumption that the baryonic processes are important only at low redshift and for the internal structure of halos is thus wrong in the WDM case. Furthermore, since there are several regimes at different stages in which different properties become significant, using the same parameters throughout the simulation and with the value taken from CDM simulations, i.e., softening length, may not allow for an accurate description of the system at those stages.
Going into the simulation analysis, methods employed in the analysis of the CDM halos, such as 'FoF', may not be accurately describing the distribution of halos in WDM, even when simulations have been performed with high resolution. When analyzing the properties of a WDM halo at redshift zero, it is crucial to look at the halo's history by tracing the particles back to a much earlier redshift or the initial conditions, first to make sure that the halo is not formed through fragmentation or through merging of other 'spurious' halos, and second, to establish what structure formation mechanism was the halo formed through, since the structure formation in WDM depends on the morphology and the architecture of the environment. For this purpose, once again, merger tree methods that have been used in the analysis of CDM simulations may not be appropriate if the halo has formed through a top-down collapse and accretion, as opposed to mergers.
To summarize, the recipe to perform WDM simulations that properly solve the subtleties of structure formation at both large and small scales has to include the following: • Proper representation of the physics in the simulations-free streaming length, velocity dispersion, particle-particle interactions. • High-resolution simulations that will allow for the proper representation of the WDM effects at all scales-similar to the resolution in [76]. • Well-suited initial conditions-box size, starting redshift, simulation parameters. • Inclusion of gas and baryonic processes early on in the simulation. • Suitable analysis with suitable analysis tools-i.e., halo formation history, mass function. • Proper representation of the physics in the analysis and interpretation of the resultsi.e., quantum pressure, phase space analysis. While these steps may be self-evident, unfortunately in practice, due to the expensive numerical treatment required by WDM, many of them are skipped or crudely approximated, leading to unreliable results and imprecise, or even incorrect, conclusions. Add to that the fact that even when performing and analyzing the simulations in a proper manner, the constraints can only be applied within a certain theoretical frame that was originally assumed-simulations cannot strongly constrain the mass of the WDM particle in generaland the question of what can we actually learn from simulations remains. We will address this question together with future numerical prospects and tests in the following section.

Learning from Simulations and Future Tests
First thing to learn is how not to model WDM, or at least how not to trust strong claims coming from inaccurate, unrefined simulations. There are, however, some hints, even if just qualitative, as to what implications a WDM particle might have on the structure formation and the internal structure of galaxies. Several scenarios in which the qualitative differences between the models can be explored in simulations and can be tested by observations are summarized below: • Since at larger scales the structure forms differently than in the CDM scenario and baryonic physics may put a twist on the structure formation, much higher resolution simulations with added baryonic processes should be performed. • The early presence of super massive black holes can be tested in this scenario. • Quantum pressure of the WDM particle should be included in simulations, since it has a significant effect on the cores. • At smaller scales, even though PSD cores can be too small for the current simulations resolution to distinguish from CDM, they do occur naturally in WDM and thus when considering quantum pressure effects and baryonic physics they would surpass the core sizes found in CDM + baryons models. • The formation of disk galaxies can be studied in this scenario, where baryons can condense coherently in a halo with a smooth potential. • Since WDM halos have visible shells and caustics, the effects of baryons should be explored in this scenario and the enhancement in the annihilation signal should be studied. • Voids and halos formed in voids may provide interesting insights into the qualitative differences between the models, differences that can be tested by observations.
Ultimately, while we cannot place an absolute constraint on the WDM particle from simulations, we can rule out some ranges of velocities (masses) that are too warm to account for the structure observed across several models. Unfortunately for detection, these loose constraints on the velocity of the particle cannot be translated into strong mass constraints in the absence of a universal model for WDM particle production.
It is important to note that while here we focused on the theoretical particle physics models that are more often considered in the simulations studies, there are continuous efforts in the field of particle physics to provide candidates for the dark matter particles, both in the theoretical and the experimental research.

Concluding Remarks
Theoretical studies show that WDM particles in the keV range can explain both the large-scale structure and the small-scale structure observations, when considering the WDM inherent properties-free-streaming length, velocity dispersion (corresponding to a particle mass in the keV regime) and quantum pressure, e.g., [11,12]. Numerical simulation studies, however, while providing an important tool to test some of these effects, have a long way to go in describing accurately some of the warm dark matter properties and their influence on structure formation and evolution.
We reviewed several numerical simulation studies that have attempted to model particles in what is generically called the warm dark matter regime. Focusing on the approximations made in simulations, both in the context of the physics employed and the numerical methods, we conducted a critical analysis of the results and the claims. Following a brief overview of the varied theoretical models of particle production we discussed the various parameters that are important in the attempt to constrain a model and how these parameters vary between the models assumed. Furthermore, we discussed the non-trivial implementation of these parameters in simulations and emphasized the importance of a proper numerical treatment of the physics involved, pointing out the limitations of the current numerical studies. Starting with the early structure formation mechanisms and continuing with the evolution of structure and the 'final' halo properties of simulated WDM distributions, we discussed the main similarities as well as the differences reported in different numerical studies. We focus on regimes in which qualitative differences with respect to CDM simulations are more apparent, i.e., structure formation mechanisms, shells and caustics in halos, internal structure of galaxies, stressing out the importance of studying further these regimes in a manner suitable for properly describing WDM properties.
Finally, we addressed several confusing claims found in the literature regarding both the warm dark matter properties and their observable signatures, and provided several insights on where the attention should be focused on further developing more accurate simulation studies that properly describe warm dark matter models.
Funding: This research received no external funding.

Acknowledgments:
The author would like to thank Norma G. Sanchez for organizing this Special Issue and is grateful for the paper submission invitation and for many useful discussions and comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: