Fuzzy Dark Matter Dynamics and the Quasiparticle Hypothesis

Dark matter may be composed of ultra-light bosons whose de Broglie wavelength in galaxies is of order 1 kpc. The standard model for this fuzzy dark matter (FDM) is a complex scalar field that obeys the Schr\"odinger-Poisson equations. The wavelike nature of FDM leads to fluctuations in the gravitational field that can pump energy into the stellar components of a galaxy. Heuristic arguments and theoretical analyses suggest that these fluctuations can be modelled by replacing FDM with a system of quasiparticles (QPs). We test this hypothesis by comparing self-consistent simulations of a Schr\"odinger field with those using a system of QPs in one spatial dimension. Simulations of pure FDM systems allow us to derive a phenomenological relation between the number of QPs that is required to model FDM with a given de Broglie wavelength. We also simulate systems of FDM and stars and find that the FDM pumps energy into the stars whether it is described by QPs or a Schr\"odinger field with the FDM adiabatically contracting and the stellar system adiabatically expanding. However, we find that QPs overestimate dynamical heating.


INTRODUCTION
In the fuzzy dark matter (FDM) scenario, dark matter is a scalar field ψ whose mass m is so small that its de Broglie wavelength λ = h/mv is on the order of a kiloparsec.It therefore exhibits wavelike behavior on galactic scales Hu et al. (2000), which could help distinguish it from ordinary cold dark matter (CDM).Indeed, it was proposed to solve purported problems with CDM such as the apparent deficit of satellites in Milky Way-sized galaxies.Though these problems may actually reflect gaps in our understanding of baryonic physics (see, for example, Bullock & Boylan-Kolchin (2017)) the level of interest in FDM remains high in part because FDM dynamics can lead to phenomena that allow one to distinguish it from other dark matter candidates (see Hui et al. (2017); Niemeyer (2020); Hui (2021) for recent reviews).In particular, the interference patterns in an FDM field result in rapid, large-scale fluctuations in the gravitational potential and hence the exchange of energy between the stellar and dark matter components of a galaxy (Schive et al. 2014;Dutta Chowdhury et al. 2021;Yavetz et al. 2022;Dutta Chowdhury et al. 2023).
Typically, FDM is treated as a non-relativistic, complex field ψ that obeys the Schrödinger equation with a potential given by the Newtonian gravitational potential Φ. Furthermore, FDM contributes ρ ψ ∝ |ψ| 2 to the density and thus acts as a source for the gravitational potential through Poisson's equation.We refer to this model of FDM as a Schrödinger field (hereafter SF).In the limit where the de Broglie wave-length is much smaller than the scales of interest, FDM behaves like ordinary collisionless dark matter.In fact, one can define an effective phase space distribution function that obeys the collisionless Boltzmann equation in the limit where the de Broglie wavelength goes to zero (see Widrow & Kaiser (1993) and Section 2.3 below).
Much of our intuition regarding the dynamics of stars in FDM-dominated galaxies comes from the hypothesis that FDM behaves like a system of quasiparticles (QPs).The transfer of energy between FDM and stars can therefore be thought of as arising from two-body interactions between the light stars and heavy QPs.Dimensional analysis suggests that the mass mQP of a QP is given by the de Broglie volume ∼ λ 3 times the local dark matter density (Hui et al. 2017;Church et al. 2019;Hui 2021).Alternately, one can write mQP ≃ (ℏ/m) 3 f (x, v) where f is the phase space density or distribution function (DF) for the dark matter and (ℏ/m) 3 is the phase space volume of a QP.In other words, the number density of QPs in phase space is constant while the mass of each QP is determined by the local DF.If rt and vt characterize the extent of the dark matter distribution in space and velocity, then the phase space volume occupied by the dark matter will be O r 3 t v 3 t and the number of QPs will be NQP ≃ (mrtvt/ℏ) 3 .Theoretical arguments based on the QP hypothesis have been used to estimate the time scale for FDM to disrupt star clusters and wide binary stars, thicken stellar discs, and draw energy out of the orbits of supermassive black holes (see Hui et al. (2017) and references therein).Amorisco & Loeb (2018) implemented a version of the QP hypothesis in their numerical simulations of the disruption of stellar streams by FDM.
Ultrafaint dwarf galaxies are a particularly attractive arena for studying the effects of FDM on stellar dynamics because they are dominated by dark matter and because the dynamical times are very short.Indeed, their very existence may rule out large regions of FDM parameter space.Dalal & Kravtsov (2022) reach this conclusion using both a heuristic argument and N-body simulations.For the former, they assume that the density fluctuations in a SF are of order unity on the scale of the de Broglie wavelength and show that the change in the variance ∆σ * of stars over time t in an FDM halo can be written where σ dm is the velocity dispersion of the dark matter and r 1/2 is the stellar half-mass radius.Dalal & Kravtsov (2022) use equation 1 to estimate the time it would take for σ 2 ⋆ to double in the UFDs Segue 1 (Belokurov et al. 2007) and Segue 2 (Belokurov et al. 2009) as a function of m.In doing so, they estimate the sensitivity of stars in these two UFDs to heating by FDM.Of course, both σ⋆ and r 1/2 are also changing and therefore equation 1 should be properly read as a contribution to the Fokker-Planck equation for the stellar DF.If changes in the stellar system are adiabatic, then the virial relation between σ⋆, r 1/2 and M will be approximately satisfied at all times.Whether σ⋆ increases, decreases, or stays the same will therefore depend on the radial density profile of the system and whether it changes with time.Dalal & Kravtsov (2022) obtained more quantitative results by simulating FDM-star galaxies.However, they treated the gravitational potential as fixed so that the FDM field could be written as a superposition of energy eigenstates, which evolved independently via a unitary, analytic transformation.This greatly reduced the computational complexity of their simulations though they forfeited self-consistency in computing the gravitational potential as well as energy conservation.
Recently, Dutta Chowdhury et al. (2023) conducted selfconsistent simulations of stellar systems that were embedded in FDM halos.The FDM halos were extracted from the cosmological simulation described in Schive et al. (2014) while the stars were introduced by treating them as a testparticle sample drawn from an equilibrium distribution function, which in turn was constructed using the time and azimuthally-averaged FDM density profile.Since the stars in their simulations accounted for less than 0.1% of the total mass, this procedure yielded a system in approximate dynamical equilibrium.Dutta Chowdhury et al. (2023) found that the stars were rapidly heated by the fluctuating FDM potential.In particular, over 7.5 Gyr the characteristic size of the stellar components increased by an order of magnitude while the velocity dispersion increased by a factor of 2 − 3. The results therefore confirmed the arguments in Dalal & Kravtsov (2022) though they cautioned that their simulations made a number of assumptions that might not apply to actual UFDs.El-Zant et al. (2019); Bar-Or et al. (2019) and Bar-Or et al. (2021) provided more formal treatments of stellar heating by FDM within the Fokker-Planck formalism.These studies confirmed the correspondence between heating of stars by QPs and by a SF.However, these analyses were carried out for homogeneous systems and were therefore not directly applicable to isolated systems of FDM and stars.
The goal of this paper is to test the QP hypothesis for iso-lated, self-gravitating systems by running side-by-side SF and QP simulations.We consider pure SF and pure QP systems as well as systems that have both dark matter and stars.All components are live and the gravitational field is determined self-consistently.All of our simulations are done in one spatial dimension and assume plane symmetry.This choice is made to keep the computational complexity at a manageable level since the complexity in simulating a self-gravitating SF in three dimensions is very high when the de Broglie wavelength is small.Furthermore, one can easily visualize the full position-velocity phase space in one dimension.
It is worth reflecting on the applicability of our results to three dimensions.Self-gravitating systems in one dimension share many properties with their three dimensional counterparts.In both cases, equilibrium systems can be set up via the Jeans theorem as described in Binney & Tremaine (2008) and below.We choose the lowered isothermal plane as initial conditions (Weinberg 1991), which is the one dimensional analog of the well-known King model (King 1966).We further note that both one and three dimensional systems exhibit vibrations about these equilibrium states (see Kalnajs (1973); Mathur (1990); Weinberg (1991); Widrow & Bonner (2015) for the one dimensional case and Binney & Tremaine (2008) and references therein for the three dimensional case).On the other hand, systems that start from cold initial conditions will undergo gravitational collapse and phase mixing regardless of their dimensionality.(See Fillmore & Goldreich (1984) for a nice example in the context of cosmological structure formation.)Of course, the force-law in a plane-symmetric one-dimensional system is very different from the three dimensional force-law.In particular, the force between two infinite planes (the one-dimensional analog of a three-dimensional particle) is independent of the separation between the particles.The scattering of two planes therefore occurs instantaneously when they pass one another and the force changes sign.Thus the details of two-body relaxation will be rather different in one and three dimensions though the essential physics of the process, which involves the transfer of energy between heavy particles and a sea of light particles, is the same.
At first glance, it is not obvious whether standard N-body methods are appropriate for QPs.In particular, the density fluctuations associated with FDM appear and disappear in a stochastic manner.Can these fluctuations really be described by particles that follow continuous Newtonian orbits?One might imagine a scheme in which QPs are randomly created and destroyed but adding such complications would seem to defeat the purpose of using QPs to model FDM.
A further examination of equation 1 suggests that a simple N-body scheme may actually capture the dynamics of QPs.We define the crossing time of the system as tcr ≡ r 1/2 /σ * .The fractional change in the stellar variance over a single crossing time is then (2) The quantity (mr 1/2 σ dm ) 3 is, up to a constant of order unity, equal to the phase space volume of the system while ℏ 3 is the volume of a quantum cell in phase space.Thus, the final factor in equation 2 is equal to 1/N cell , where N cell is equal to the number of cells in the system, again up to a constant of order unity.This expression is reminiscent of the usual two-body relaxation formula if we identify N cell with the number of "particles" in the system (Binney & Tremaine 2008) and motivates the idea that QPs have equal phase space volume.
The mass of each QP is then given by the local value of the DF multiplied by this volume.
As we'll see, the QP hypothesis works extremely well for systems that stay close to their equilibrium distribution.For example, the spectrum of vibrations about an equilibrium state is remarkably similar whether the system is modelled as a SF or system of QPs.Conversely, when the system evolves significantly away from its initial state, the behaviour of a SF and a system of QPs can be very different.In particular, we find that over long times, QPs are more efficient at heating stars than a SF, a result that may have implications for understanding the evolution of dwarf galaxies in an FDM cosmology.
In Section 2, we describe various numerical methods used in our simulations.In Section 3, we present results for a pure FDM system modelled either as a SF or QPs.In both cases, we find that the system vibrates due to fluctuations in the density and potential.The amplitude of the vibrations can be characterized by fluctuations in the virial ratio V = 2K/W where K and W are the kinetic and potential energies of the system.We find that the fluctuations in V are approximately proportional to λ 1/2 or, alternatively, N −1/2 QP , as one would expect from root-N statistics.These results allow us to derive a correspondence between NQP and m.We also examine the power spectra for the gravitational force in each of the simulations.Results from simulations that include stars and either FDM or QPs are given in Section 4. We show that in both cases, the stars gain energy at the expense of the FDM or QPs.Nevertheless, the stars maintain approximate virial equilibrium throughout with no long term drift in V. We conclude in Section 5 with a summary of our results and some thoughts on their applicability to systems in three dimensions.

PRELIMINARIES
In this section we present the numerical methods used to set up and run simulations with a SF, a system of QPs, and stars.In all of our simulations we start from equilibrium initial conditions.To this end, we first write down a DF that solves the time-independent collisionless Boltzmann equation.For the SF, we set up initial conditions for the wave function ψ so that its phase space representation is approximately equal to the same equilibrium DF used for the QPs and stars.

Distribution Function
Consider a system in one dimension with DF f = f (x, v).In a static potential Φ(x) any function of the specific energy (i.e., energy per unit mass) E = 1 2 v 2 + Φ will be a solution to the time-independent collisionless Boltzmann equation.In this work, we take the initial DF to be that of the lowered isothermal plane (Weinberg 1991).This model is a truncated version of the isothermal plane first considered by Spitzer (1942) and Camm (1950) and has been used to describe the vertical structure of disc galaxies such as the Milky Way.In the original model, the phase space density went to zero only in the limit |x| → ∞ or |v| → ∞.In the lowered isothermal plane, the phase space density goes to zero at finite x and v and is therefore suitable as a starting point for numerical simulations.The DF for the lowered isothermal plane is given by where f0 is a normalization constant, σ is a velocity scale, and E0 is the energy cut-off.Note that f has physical dimensions of If we integrate the DF over the x − v volume, we obtain the total surface density Σ, which is the analog of the total mass in a three-dimensional system.The density ρ is found by integrating where vm(x) ≡ 2 (E0 − Φ(x)) is the maximum velocity at position x and um ≡ vm/ √ 2σ.Likewise, the variance of the velocity as a function of the potential is The density, force, and potential are found as functions of x by numerically integrating Poisson's equation.We define xt as the position at which the density goes to zero.The total surface density is then given by Σ = 2 x t 0 dxρ(x).We choose units in which G = σ = 1 and adjust f0 so that Σ = 1/π.We then define the unit of length to be x0 ≡ σ 2 /πGΣ = 1.The velocity dispersion is maximal at x = 0 with a value that depends on σ and E0 and decreases smoothly to zero as x approaches xt.
In the limit E0 → ∞, the model reduces to the isothermal plane (Spitzer 1942;Camm 1950).In that limit the potential and density are elementary functions of x and the velocity dispersion is a constant and equal to σ.For example, the potential for the isothermal plane is given by In what follows, we use E0 = 3.We then find that xt ≃ 2.00 and vt ≡ √ 2E0 ≃ 2.45.For x = 0, the density and velocity dispersion are ρ0 ≡ ρ(x = 0) = 0.204 and ⟨v 2 ⟩ 1/2 (x = 0) = 0.78.The time for a particle of energy E to complete one orbit in the x − v plane is which is a monotonically increasing function of E. We define the dynamical time as t dyn ≡ T (E → 0).Since Φ(x) ≃ 2πρ0x 2 for x ≪ 1, we have t dyn = π/ρ0 ≃ 3.92.Finally, the system occupies a phase space volume Ω = 4 which is just the usual action divided by m.Thus, if there are NQP QPs, each particle will occupy a phase space volume ΩQP = Ω/NQP .For E0 = 3, we find Ω ≃ 14.5.

Equations of Motion
In the SF model, FDM obeys the Schrödinger-Poisson system of equations.In one dimension we have where we have introduced the dimensionless parameter R ≡ ℏ/mσx0 = ℏ/m.The QPs and stars are modelled as point particles whose orbits are determined by Newton's second law.
The density associated with ψ is given by m|ψ| 2 .For stars and QPs, we calculate the density ρp using a particle mesh scheme where the mass of each particle is assigned to the cell whose center is closest to the position of the particle.We then have Numerical evolution of FDM is handled by the kick-driftkick scheme that is now widely used in both cosmological simulations and simulations of isolated systems (Woo & Chiueh 2009;Edwards et al. 2018).The Poisson equation is solved using standard Greens function and FFT methods with zeropadding to handle the boundary conditions of an isolated system (Hockney & Eastwood 2021).

Initial Conditions
Our goal is to set up an initial state ψ(x) whose phase space representation is approximately equal to equation 3. We do so using the algorithm outlined in Widrow & Kaiser (1993).To make the connection between ψ and f we use the Husimi Q transformation (Husimi 1940) in which the effective DF is the absolute square of the windowed Fourier transform of ψ: where and N is a normalization constant.In the limit ℏ → 0, F obeys the familiar continuity and Jeans equations for a system of collisionless particles (Madelung 1926;Skodje et al. 1989;Widrow & Kaiser 1993).In general we choose η = (Rxt/vt), which gives equal resolution in space and velocity.
To evaluate equation 12, we set up a phase space grid with and where ∆x ≡ 2xmax/N and ∆v ≡ 2vmax/N .We set xmax and vmax to be 1.5 times larger than xt and vt so that the grid can also be used to map ψ in phase space even after the system has evolved away from its initial state.Equation 12then takes the form of a discrete Fourier transform provided we set N = 2xmaxvmax/πR.We next consider the ansatz where {R k } k is a sample of a random variable uniformly distributed on the complex unit-circle.A straightforward calculation shows that the phase space representation for ψ yields F ≃ f (Widrow & Kaiser 1993).
As discussed in the introduction and also El-Zant et al. (2019); Bar-Or et al. (2019) and Bar-Or et al. (2021), in three dimensions, we expect each QP to occupy a (positionvelocity) phase space volume of (ℏ/m) 3 .In one dimension, we therefore expect QPs to occupy an x − v phase space volume of ℏ/m.The surface density of each QP is then given by ℏ/m times the distribution function f .To initialize the QP distribution we uniformly sample the phase space volume Ω for points {(xj, pj)}j.The j th QP resides at the coordinates (xj, pj) and is given surface density N f (xj, pj) where the normalization factor N is adjusted so that the total QP surface density is equal to Σ.The QP distribution is then rather different from the normal distribution in N-body simulations where simulation particles, at least for a particular system such as a disc, bulge, or satellite, have the same mass.In that case, regions of high phase space density have a higher density of simulation particles.In the QP case, phase space is uniformly sampled but the particles have different masses, or rather different surface densities since we are in one dimension.In this regard, our QP simulations are similar to multi-mass N-body simulations, as discussed, for example, in Sigurdsson et al. (1995).For collisionless simulations, multimass methods afford one higher mass resolution for fixed computational cost.In our case, it's the collisional properties of the QPs that we wish to study while the scheme chosen is motivated by the physics of FDM.

SF VERSUS QPS IN FDM SIMULATIONS
In this section we compare the evolution of SF and QP systems by running two sequences of simulations with different values of either R or NQP .In all of the simulations, the initial conditions are derived from equation 3 with E0 = 3, σ = 1, and f0 adjusted so that Σ = 1/π.The simulations are run up to T = 96 ≈ 24t dyn .We solve Poisson's equation on a mesh using a grid with N cells.We set N = 500 for QP simulations.As discussed in the introduction, the number of cells for the SF simulations is determined by the ratio of λ ∼ ℏ/mσ dm to the size of the system (∼ r 1/2 ).The timestep h is chosen to satisfy the Courant condition h = 0.5∆x/vt.Note that while this gives h ≈ 0.0025 for QPs, in the case of the SF the meshsize changes between simulations, and thus different values of h are used.This is a choice that allows us to lower the computational time of SF simulations at higher values of R.

Oscillations
Both SF and QP systems oscillate about their equilibrium configurations.To illustrate this, we choose two independent QP and SF systems with parameters NQP = 100, h ≈ 0.0025, N = 500 and R ≈ 0.067, h ≈ 0.0148, N = 84 respectively.In Fig. 1, we plot the time evolution of the kinetic and potential energies, K and W , as well as the virial ratio V ≡ 2K/W .The oscillations in K, W , and V are qualitatively very similar in the two systems.The dominant period is approximately P ≃ 3 ≃ 0.8t dyn while the beat patterns indicate  that there are multiple oscillations with different frequencies.These oscillations are characteristic of one-dimensional self-gravitating systems and have been studied by Kalnajs (1973); Mathur (1990); Weinberg (1991) and Widrow & Bonner (2015).
To further illustrate the similarities between SF and QP dynamics, we show the temporal power spectra PV (ω) of V in Fig. 2. The spectra are remarkably similar.In particular, there is a prominent peak at ω ≃ 2, which corresponds to a period of P ≃ 3, as expected from Fig. 1.In addition, there are secondary peaks at multiples of the first one as well as a continuous distribution with P ∝ ω −4 at large ω.We'll return to this last point below.
To derive a phenomenological relation between R and ΩQP we run an SF sequence with 0.003 < R < 0.15 and a QP sequence with 0.003 < ΩQP < 1.5.In the case of the latter, the Figure 3. Curve fits of the RMS amplitude's of oscillation in V, plotted against either Ω QP or R on log-log axes.The dotted black curve is the mean fit with the power (p or q) fixed to 1/2.The orange lines correspond to a sample of p and q values from the emcee routine.Each data point is the average of 15 trials for the given parameter.Similarly, the uncertainties of each data point are the standard deviation of each sample of trials.
number of QPs varies between 10 and 5000.In each simulation, we calculate δV ≡ ⟨(V − 1) 2 ⟩ 1/2 , that is, the root mean square deviation in V from unity over the course of the simulation.For each value of R or ΩQP , we run 15 simulations with different realizations of the initial conditions and compute the mean and standard deviation of δV.We then model δV as a power-law function of R or ΩQP : and The parameters a, p, b, and q are determined by fitting log V to a straight line using the Markov chain Monte Carlo sampler emcee Foreman-Mackey et al. (2013).The results are shown in Fig. 3 and the fit parameters are given in Table 1.
Since the fluctuations are seeded by Poisson noise from the initial conditions, we anticipate that δV will be proportional to R 1/2 or Ω 1/2 QP .This expectation is borne out in the SF simulation.In the case of the QPs, we find p ≃ 0.58 ± 0.015.Though the discrepancy from p = 0.5 is statistically significant, it amounts to only a ±15% change in δV over two orders of magnitude in ΩQP.Table 3.1 also includes values for a and b when p and q are fixed to 0.5 Setting δVFDM equal to δVQP leads to the power-law relation The joint and marginal probability distribution functions for α and β are shown in Fig. 4 while best-fit values for α and β are given in Table 1.

Force Power Spectra
In Fig. 5 we show power spectra for force fluctuations, PF (k), from SF simulations with R ≈ 0.0084, 0.017, 0.034, 0.067, 0.13 and from QP simulations with the corresponding values of NQP as given by the relation derived in the previous section.
In these simulations, we fix the number of grid points to N = 500, for a more direct comparison of the power spectra.To  3 and for the relation between Ω QP and R. The third column shows the values given p and q fixed at 1/2.calculate PF (k), we first determine the force fluctuations as a function of x and t by subtracting off the time-averaged force: We next compute the spatial Fourier transform, which yields In the case of QPs, PF (k) ∝ k −2 down to kt ≃ 2π/xt.The k −2 spectrum corresponds to Brownian noise.In one dimension, the gravitational force at position xp is proportional to the difference between the surface density for x > xp and the surface density for x < xp.Since the density can be modeled as white noise, its Fourier transform will be constant in k and therefore the Fourier transform for the force or surface density will be proportional to k −1 .Fig. 5 also provides a connection to the temporal Fourier transform of V found in Fig. 2. Since the collisionless Boltzmann equation is linear and first-order in both space and time derivatives, we expect that k ∝ ω for the two-dimensional (spatial-temporal) Fourier transforms of various quantities such as the potential and force.Furthermore, both the oneand two-dimensional Fourier transform of the force are −ik times the Fourier transform of the potential: and the ω −4 behaviour for the temporal Fourier transform of V is consistent with the k −2 behaviour for spatial Fourier transform of F .In our SF simulations the force power spectra exhibit a fall off for length-scales below λ.Our Fig. 5 is analogous to results from Dalal et al. (2021) (see the ∆t = 0 curves in their Figure 13).This damping of power below the de Broglie scale was one of the initial reasons FDM was introduced and meant that in the context of structure formation, FDM had some of the same properties as warm dark matter Hu et al. (2000).

DYNAMICAL HEATING OF STARS BY FDM
In this section we compare the dynamical effects of a SF verses QPs on a system of stars.For illustrative purposes, we assume that dark matter and stars each account for half the mass of an isolated system and that both components start with the same equilibrium DF used in the previous section.We run two independent simulations with 5 × 10 4 particles to represent the stars and either a SF or QPs for FDM.The particles representing the stars all have the same mass; their initial positions and velocities are initialized by sampling the DF using a simple accept-reject algorithm.We set R ≃ 0.067 for FDM.With this value our previous results suggest using NQP = 106 if p and q are taken as free parameters or NQP = 84 if p and q are fixed to 1/2.Here, we set NQP = 106.In both simulations, we fix the number of grid points to N = 500, and the time step to ∆t ≈ 0.0025.

Phase space evolution of stars and dark matter
In Fig. 6 we present a sequence of phase space snapshots for the SF and QP runs.In both cases, the evolution proceeds through three stages: an initial isothermal phase during which the DF is peaked at the origin of the x − v plane; an intermediate phase where the DF is disturbed by either the SF or QPs; a final phase when the DF is approximately constant for E less than the truncation energy.Interestingly enough, the edge of the distribution doesn't change by very much over the course of the simulation.Rather, there appears to be a redistribution of particles from small to large energies.
These points are further illustrated in Fig. 7.The energy distribution f (E) shows a transition from an isothermal or Maxwellian distribution at the start of the simulation to one that is approximately constant in energy for E ≲ 1.5 .The shift of particles from low-to-high energy is shown in greater detail by the scatter plots of Fig. 8.While the distribution of initial and final energies show an overall increase in energy, it is clear that the heating effect is more pronounced for the particles with E initial ≲ 1.
In addition to the phase space plots of the SF seen in Fig. 6, we plot the initial and final SF density in Fig. 9.The characteristic size of the SF shrinks in both position and velocity space, from root-mean-square values xrms ≃ 0.598 and vrms ≃ 0.823 to xrms ≃ 0.435 and vrms ≃ 0.569.This demonstrates the losses of both potential and kinetic energy in the SF over the course of the simulation.

Evolution of the half-mass radius and velocity dispersion
As discussed in the introduction, Dalal & Kravtsov (2022) argue that large regions of FDM parameter space are ruled out by the existence of dark matter dominated UFDs such as Segue 1 and Segue 2. The essence of their argument is that if the dark halos of these systems were composed of FDM with mc 2 ≲ 3 × 10 −19 eV, then the stellar dispersion σ⋆ and projected half-light radius R 1/2 would grow over time to the point of being inconsistent with observations.This argument is illustrated in their Figure 2, which shows R 1/2 and σ⋆ as functions of time for various values of m and for different initial conditions for the stars.They find that both R 1/2 and σ 2 ⋆ roughly double over ∼ 10 Gyr.Since these systems have dynamical times of order 1 Myr, one can assume that approximate virial equilibrium is maintained.Thus, the virial mass one would infer from the stars, Mvir ≃ 4σ 2 ⋆ R 1/2 (Wolf et al. 2010), will have increased by a factor of four, which is consistent with having a stellar system that expands by a factor of two within a cuspy dark halo.
We now carry out a similar analysis in our star + SF/QP simulations.We first note that for our initial conditions, σ⋆ ≃ 0.842 and x 1/2 ≃ 0.423 where, by definition, half of the total surface density is contained in the region |x| < x 1/2 .Therefore, one can define a virial surface density, in analogy with the virial mass, as Σvir ∝ σ 2 * /x 1/2 .The constant of proportionality, which is obtained by setting Σvir = 1/π, the total surface density, and σ * and x 1/2 to their initial values, is 0.19.
In Fig. 10 we compare the evolution of x 1/2 and σ⋆ for stars embedded in a SF halo and stars embedded in a halo of QPs.For each case, we run four simulations for ≃ 375t dyn .We find that in the star+QP simulations, x 1/2 and σ⋆ increase steadily by 30% and 20%, respectively.The virial surface density Σvir and hence the effective mass are roughly constant during the course of the simulation.
For t ≲ 50t dyn the evolution of x 1/2 and σ⋆ in the FDM and QP simulations are very similar.However, at later times, x 1/2 and σ⋆ are nearly constant in the stars + SF simulations and the overall increase in x 1/2 and σ 2 ⋆ is about half what it is in the stars + QPs simulations.To further explore the differences between QP and FDM simulations, we rerun the simulations with twice as many QPs.As expected, the change in x 1/2 and σ⋆ is decreased relative to the original QP simulation though it is still higher than the change found in the FDM simulations.Moreover, the heating at earlier times is of a lower rate.

SUMMARY AND DISCUSSION
Our aim in this paper has been to test the hypothesis that FDM can be treated as a system of QPs using standard N-body techniques.The main difference between QPs and the N-body systems normally encountered in cosmology and galactic dynamics is in how the QP masses are assigned.In our implementation of the quasiparticle hypothesis, QPs are uniformly distributed in phase space with a mass proportional to the DF at their initial location.
Our results regarding the validity of the QP hypothesis are mixed.Isolated systems of dark matter exhibit small oscillations about an equilibrium state whether they are modeled by QPs or a SF.Furthermore, the amplitude of oscillations in the virial ratio and the frequency power spectra of these oscillations are very similar so long the relation between NQP and R given by equation 18 is satisfied.The main difference is in the force fluctuation power spectrum PF (k) in that the SF power spectrum shows cut-off for scales below the de Broglie wavelength, which isn't present in the power spectra from QP simulations.
We also showed that FDM can dynamically heat stars whether it is described as a SF or system of QPs.In fact, at early times, the rate at which the half-mass radius and velocity dispersion increase is very similar.However, at later times, the QPs appear to be more efficient at heating the stars than the SF.
We contend that QPs, at least as we've modelled them, do best as a proxy for a SF in situations where the structure of the FDM component does not change appreciably with time.For example, QPs should provide a useful substitute for a SF in problems such as the disruption of stellar streams Amorisco & Loeb (2018); Dalal et al. (2021) or the heating of a disc in an FDM halo Hui et al. (2017).QPs can be reliably used in systems where dark matter dominates the potential and the stars can be treated as test particles, such as the analysis of UFDs in Dutta Chowdhury et al. (2023) where the stars make a negligible contribution to the potential.On  the other hand, there are many situations where the FDM distribution function changes significantly.We've given one example where the contributions to the potential from FDM and stars are comparable and where the FDM halo adiabatically condenses as the stellar system adiabatically expands.In these situations, the transfer of energy from a SF to stars may be significantly different than that for a system of QPs and the QP hypothesis might lead to erroneous conclusions.Of course, these conclusions were reached through one dimensional simulations.Though the essential physics of dynamical heating in one dimension is the same as in three dimensions, the details are very different.Ultimately, it will take three dimensional numerical experiments similar to the ones performed here to test the applicability of the QP hypothesis.The dashed black curve corresponds to a simulation using double the number of QPs.Each curve is the average over 4 simulations.Similarly, the light bands represent the standard deviation across the 4 simulations for each regime.Note that, for the sake of visibility, all curves have been smoothed using a Gaussian filter.
analysis in Section 3 was performed using the Markov chain Monte Carlo sampler emcee Foreman- Mackey et al. (2013).Figure 4 was produced using the python package corner.py(Foreman-Mackey 2016).The data for the figures and the code will be shared on reasonable request to the authors.

Figure 1 .
Figure1.Kinetic and potential energies and the virial ratio as a function of time for the SF (left column) and QPs (right column).The upper panels show K (red) and W (blue) while the lower panels show V ≡ 2K/W .

Figure 2 .
Figure 2. Temporal power spectra of V for the SF (purple) and QPs (orange).

Figure 4 .
Figure 4. Corner plot showing the joint and marginal distributions of α and β.The blue lines represent the 50th percentile, while the dashed black lines are the 0.16 and 0.84 quantiles, representing ±1σ.

F
′ (k, t).The desired force fluctuation power spectrum is found by taking the time average of | F ′ | 2 .Note that PF is just the Fourier transform of the force auto-correlation function, which plays a central role in the Fokker-Planck analyses of Bar-Or et al. (2019); El-Zant et al. (2019), and Bar-Or et al. (2021).

Figure 5 .
Figure 5. Spatial power spectra of oscillations from the timeindependent acceleration field −∇Φ, for SFs with differing values of fuzziness R, and the corresponding number of quasiparticles N QP (R).The colours (dark blue, light blue, green, orange, red) correspond to the cited values of R, in increasing order.The dotted black line marks the truncation length xt, while the dashed vertical lines mark the respective de Broglie wavelengths λ.Note that the horizontal axes represent the wavelength 2π/k, where k is the spatial frequency.

Figure 6 .
Figure 6.Snapshots of the phasespace distributions of the QPs+stars (top) and SF+stars (bottom), taken at approximately 0, 7.5, 37.5, 75, 375 dynamical times.Particles representing the stars are shown as black points while the QP/SF portions are plotted in the upperright corners of the first and last snapshots.In the case of QPs, the colour signifies the mass.In the SF case, the colour represents the phase-space density given by the Husimi Transform, plotted as countours at 15 levels on a log-scale.

Figure 7 .
Figure 7. Initial (dotted) and final (solid) distributions of specific energies among particles in either mixed regime; QPs+stars (orange) or SF+stars (purple).These curves exclude the energies of the QP/SF component.

Figure 8 .
Figure 8. Scatter plot of the initial vs final specific energies of every particle heated by QPs (left) or SF (right).The color scale is linear in density in the E i − E f space.(The overall scale of the density is irrelevant.

Figure 9 .
Figure 9. Density of the SF at the start (dotted) and end (solid) of the simulation.Averaged between the left and right side of the box.

Figure 10 .
Figure10.Half-mass radius (top) and velocity dispersion (bottom) of the particles over time in either regime.Orange and purple correspond to 5×10 4 stars mixed with either QP or a SF, respectively.The dashed black curve corresponds to a simulation using double the number of QPs.Each curve is the average over 4 simulations.Similarly, the light bands represent the standard deviation across the 4 simulations for each regime.Note that, for the sake of visibility, all curves have been smoothed using a Gaussian filter.

Table 1 .
Fit parameters for Figure