Galactic disc heating by density granulation in fuzzy dark matter simulations

Fuzzy dark matter (FDM), an attractive dark matter candidate comprising ultralight bosons (axions) with a particle mass $m_a\sim10^{-22}$ eV, is motivated by the small-scale challenges of cold dark matter and features a kpc-size de Broglie wavelength. Quantum wave interference inside an FDM halo gives rise to stochastically fluctuating density granulation; the resulting gravitational perturbations could drive significant disc thickening, providing a natural explanation for galactic thick discs. Here we present the first self-consistent simulations of FDM haloes and stellar discs, exploring $m_a=0.2-1.2\times10^{-22}$ eV and halo masses $M_\text{h} = 0.7-2.8\times10^{11}$ M$_\odot$. Disc thickening is observed in all simulated systems. The disc heating rates are approximately constant in time and increase substantially with decreasing $m_a$, reaching $dh/dt \simeq 0.04$ ($0.4$) kpc Gyr$^{-1}$ and $d\sigma_z^2/dt \simeq4$ ($150$) km$^2$s$^{-2}$Gyr$^{-1}$ for $m_a=1.2$ ($0.2$) $\times10^{-22}$ eV and $M_\text{h} =7\times10^{10} \text{M}_\odot$, where $h$ is the disc scale height and $\sigma_z$ is the vertical velocity dispersion. These simulated heating rates agree within a factor of two with the theoretical estimates of Chiang et al., confirming that the rough estimate of Church et al. overpredicts the granulation-driven disc heating rate by two orders of magnitude. However, the simulation-inferred heating rates scale less steeply than the theoretically predicted relation $d\sigma^2_z/dt \propto m_a^{-3}$. Finally, we examine the applicability of the Fokker-Planck approximation in FDM granulation modelling and the robustness of the $m_a$ exclusion bound derived from the Galactic disc kinematics.

Persistent local perturbations in the background gravitational potential, agnostic to the sources, can lead to significant stellar heating over a Hubble time (e.g.Chandrasekhar 1942;Chavanis 2013).The Galactic thick disc populations, with stellar ages ≳ 8 Gyr, exhibit substantially higher velocity dispersion than the kinematically cold and young thin disc populations (e.g.Bland-Hawthorn & Gerhard 2016;Miglio et al. 2021).This morphological evolution of disc thickening is increasingly well constrained (e.g.Ting & Rix 2019;Mackereth et al. 2019;Sharma et al. 2021), albeit uncertain in the underlying physical mechanisms.Beyond the MW, the formation of old and kinematically hot thick disc components (Yoachim & Dalcanton 2005, 2008) are also observed in external galaxies across wide ranges of Hubble types, stellar masses, stellar mass-to-light ratios, and thick-to-thin disc mass ratios (Tsikoudi 1979;Burstein 1979;van der Kruit & Searle 1981;Dalcanton & Bernstein 2002;Yoachim & Dalcanton 2006;Comerón et al. 2014).
The exact thick disc formation pathway is still debated in the liter-ature, and various qualitatively distinct models have been proposed and can be largely grouped into two classes.In the first category, disc stars are kinematically cold at birth, and the morphologically thin disc is dynamically heated over time to form the thick disc.Some possible sources of stellar heating include the gravitational interactions of disc stars with infalling satellites (Toth & Ostriker 1992;Quinn et al. 1993), giant molecular clouds (Spitzer & Schwarzschild 1951;Lacey 1984), spiral arms (Barbanis & Woltjer 1967;Carlberg & Sellwood 1985), or the Galactic bar (Minchev & Famaey 2010).These secular processes could also drive radial migration of disc stars (Schönrich & Binney 2009;Loebman et al. 2011).The second class of models advocates that, via violent and rapid gravitationalinstability-driven heating, disc stars are already kinematically hot during the initial assembly process to produce a morphologically thick disc.This scenario could be achieved via gas-rich major mergers (Brook et al. 2004), early turbulent and clumpy phases of disc evolution (Kroupa 2002;Bournaud et al. 2009), clustered star formation (Kroupa 2002), and/or direct accretion of tidally stripped satellite galaxies (Abadi et al. 2003;Martin et al. 2004).Chiang, Ostriker & Schive (2023), henceforth C23, recently explored another possible Galactic thick disc formation scenario in a Fuzzy Dark Matter (FDM) universe.As a promising alternative to cold dark matter (CDM), the FDM paradigm reproduces the largescale successes of the ΛCDM cosmological model (Schive, Chiueh & Broadhurst 2014a).The ultralight mass scale   ∼ 10 −22 eV yields a kpc/sub-kpc-size de Broglie wavelength and could address several small-scale discrepancies of ΛCDM (e.g.Schive et al. 2014b;Hui et al. 2017;Marsh & Pop 2015;Calabrese & Spergel 2016;Chen et al. 2017;Leung et al. 2018;Wasserman et al. 2019;Amruth et al. 2023); on the flip side, this mass range is disfavoured in some other analyses (e.g.Bernal et al. 2018;Schutz 2020;Hayashi et al. 2021;Rogers & Peiris 2021;Nadler et al. 2021;Powell et al. 2023; see however C23).
Quantum wave interference inside an FDM halo gives rise to ubiquitous de-Broglie-scale density granulation that are stochastically distributed and undamped on cosmological timescales (e.g.Schive et al. 2014a;Veltmaat, Niemeyer & Schwabe 2018).In an MW-sized host halo, these locally fluctuating granules weight comparably to giant molecular clouds ∼ 10 6 M ⊙ .By incorporating the full baryon and dark matter distributions of the MW together with stellar disc kinematics inferred from Gaia, APOGEE, and LAM-OST surveys, C23 provided a detailed estimate of the Galactic disc heating rate based on the Fokker-Planck approximation formalism first developed in Bar-Or, Fouvry & Tremaine (2019).We found that granulation-and subhalo-induced stellar heating can quantitatively reproduce both the observed 'U-shaped' disc vertical velocity dispersion profile (Sanders & Das 2018) and the age-velocity dispersion relation in the solar neighbourhood (Sharma et al. 2021).This particular observable places a conservative exclusion bound   ≳ 0.4 × 10 −22 eV and favours   ≃ 0.5-0.7 × 10 −22 eV.
The main aim of this paper is to numerically verify the analytical estimates of FDM-granulation-driven stellar heating rate derived in C23.We present the first self-consistent simulations of FDM haloes with embedded -body stellar discs, covering a range of FDM particle masses   = 0.2-1.2× 10 −22 eV and host halo masses  h = 0.7-2.8× 10 11 M ⊙ .The granulation-driven stellar disc thickening is present in all simulated systems and could provide a natural mechanism for the thick disc formation observed in external galaxies with assembly histories likely dissimilar from that of the MW.
We begin in Section 2 with a brief review of the analytical framework for estimating FDM granulation-driven disc heating rates presented in C23.An overview of the simulation setup is presented in Section 3, where the suitable initial condition construction methods for FDM haloes and galactic discs are detailed.Section 4.1 presents the simulation results, with a particular focus on the galactic disc morphological evolution in FDM haloes with different   and  h .In Section 4.2, the halo-granulation-induced stellar heating rate measured in self-consistent FDM simulations is compared directly against the analytical predictions of C23.In Section 5, the relevant implications on the intrinsic uncertainty of hydrostaticequilibrium-based   in observed external disc galaxies are discussed.We conclude in Section 6. Technical details are organised as follows: FDM halo-disc co-relaxation (Appendix A), numerical convergence tests (Appendix B), and heating in CDM halo-disc simulations (Appendix C).
In this paper, cylindrical coordinates (, ) denote the axial distance along and vertical height from the mid-plane in a galactocentric frame;  ≡ √︁  2 +  2 gives the radial distance from the instantaneous galactic centre.For notational convenience, we introduce a dimensionless FDM particle mass  22 ≡   /10 −22 eV.

THEORY REVIEW
In this section, we briefly review the analytical formalism of C23 for estimating FDM-granulation-driven galactic disc heating rates.Consider an axisymmetric, locally isothermal galactic disc embedded in a spherically symmetric background halo potential Φ bg .The disc three-dimensional density profile  d (, ) and scale height ℎ() are determined by solving the appropriate Poisson equation and vertical hydrostatic equilibrium condition (e.g.Lacey & Ostriker 1985) where Φ tot ≡ Φ d + Φ bg , and the effective background density reads where  cir denotes the total circular velocity.Exact closed-form solutions to Eq. (1) exist only in either the self-gravity-dominated (SGD) limit  d ≫  eff bg or the background-dominated (BGD) limit with corresponding disc scale heights (C23) The transition between these two limits can be estimated by equating the two expressions of Eq. ( 4) where  ≪ 1 (≫ 1) corresponds to the SGD (BGD) limit.For individual disc stars, the vertical oscillation period reaching the maximum vertical displacement ± max away from the mid-plane  = 0 takes the form Quantum-wave-interference-induced stochastic density fluctuations inside an FDM halo, under a suitable set of assumptions listed below, can be treated as classical quasi-particles with a typical effective mass of (Bar-Or et al. 2019) where  h denotes the local halo density.The de Broglie wavelength  dB is defined as 1 where ℏ is the reduced Planck constant.
where the dimensionless factor () absorbs the radius dependence of all the non-  velocity components appearing in the Fokker-Planck formulation.For the MW stellar disc, () ≃ 3-7 varies by about a factor of two across  = 2-20 kpc (C23).The first-order and second-order diffusion coefficients describing the dynamical evolution of a test star  orbiting in a granular FDM halo are (Bar-Or et al. 2019) where and therefore  = ℏ/( √ 2   2 h ).In the limit  gra ≫ , the ensemble-averaged change rates in the disc vertical velocity dispersion squared in the SGD and BGD limits are (C23) 1 It is worth stressing that the genuine halo one-dimensional velocity dispersion  h (i.e.kinetic energy sourced by the turbulent motion) differs by definition from the total effective velocity dispersion  Jeans (i.e.kinetic energy + quantum internal energy) obtained from solving the Jeans equation.
One must take equipartition of energy  h =  Jeans / √ 2 into account for an FDM halo, unlike the CDM counterparts; see Dutta Chowdhury et al. (2021) and C23.This physical distinction is crucial in accurately estimating the FDM granulation-driven disc heating rate that scales as ∝  −6 ℎ ; see Eq. ( 14).
For a sufficient large Coulomb factor ⟨⟩  / ≫ 1, the stellar heating rate ⟨H ⟩  is comparatively insensitive to ⟨⟩  and roughly scales with M. The associated characteristic heating timescale driven by FDM halo density granulation is to leading order identical in both the SGD and BGD limits.The granulation-driven disc heating rate H ∝  −1 heat is sensitive to the halo attributes  h and  h .In contrast, the   -dependence enters H only through the vertical oscillation period in the SGD limit and is logarithmically suppressed.The disc vertical heating rate is thus independent of   () to zeroth order.For a given FDM halo (i.e.fixing   ,  h , and  h ), the approximately time-independent H implies that  2  and the disc scale height ℎ ∝  2  in the SGD limit Eq. ( 4) grow linearly with time, a prediction to be compared with simulation results in Section 4.2.

SIMULATION METHODS
We describe the initial condition construction of FDM haloes in Section 3.1 and galactic discs in Section 3.2.The adaptive mesh refinement (AMR) criteria suitable for simulating dynamical coevolution of FDM haloes and -body galactic discs are detailed in Section 3.3.We outline the analysis of disc properties in Section 3.4.

Initial conditions of FDM haloes
The FDM halo initial conditions are constructed using the algorithm presented in Lin et al. ( 2018) that, provided an input density profile, iteratively solves for self-consistent wave distribution function from the Schrödinger-Poisson equations.The distribution function for each wave function eigen-state follows the King model (King 1966) which dependents merely on the eigen-energy .The parameters , , and   are chosen properly for each halo, such that the final density profile of the outer halo can roughly follow the assigned NFW profile (Navarro et al. 1996) with halo concentration parameter  h and halo virial mass  h given as input, meanwhile the total FDM halo mass inside the virial radius  vir can be kept as close to the desired  h as possible.We also confirm that the selected parameters can render a convergent final density profiles after sufficient iterations.Alternatively, one can also use other distribution functions to construct the FDM halo, such as the distribution function obtained from the Eddington formula (see Dalal et al. 2021).In this work,  vir is defined as the radius within which the mean enclosed density is 347 times the critical density.We construct four FDM haloes of identical virial masses  h = 7 × 10 10 M ⊙ with different FDM particle masses  22 = 0.2, 0.4, 0.8, and 1.2, all yielding  vir ≃ 102.6 kpc.Fixing  22 = 0.4, we add two additional cases with heavier halo masses  h = 1.4 × 10 11 M ⊙ and 2.8 × 10 11 M ⊙ that have  vir 129.3 kpc and 162.9 kpc, respectively.The adopted halo properties are listed in Table 1.
The physical properties of the central soliton core follow the soliton core-halo relation of Schive et al. (2014b) and the soliton scaling relation where  0 is the time-averaged peak density of the soliton and  c is the half-density radius  soliton ( c ) =  0 /2.The soliton core extends out to ≃ 3.3 c (Chiang et al. 2021) 1 (solid), compared with an NFW profile with  h = 7.0 × 10 10 M ⊙ and  h = 8 (dashed).For a fixed halo mass  h , the central soliton core exhibits higher peak density and smaller core radius  c for heavier FDM particle mass  22 .Beyond the characteristic density profile transition at ≃ 3.3 c , the halo density distribution roughly follows the NFW profile.
are presented in Fig. 2. The ubiquitous granular structures in each halo density slice arise from the quantum wave interference.The granules have characteristic length scales comparable to the local de Broglie wavelength Eq. ( 8).The typical physical sizes of density granulation increase for smaller FDM particle masses  22 .In this work, we quantify the disc heating effect caused by these fluctuating granular structures.

Initial conditions of galactic discs
The stellar discs are constructed using GALIC (Yurin & Springel 2014), a code that generates -body halo-disc initial conditions by iteratively optimizing the velocities of populated particles to reach an approximate equilibrium solution to the collisionless Boltzmann equation.However, the predefined analytic density distribution functions implemented in GALIC only support CDM haloes parameterised by the Hernquist profile.We have thus modified the source code to take any shell-averaged FDM halo density profiles constructed in Section 3.1 as inputs, and output initial conditions of galactic discs and -body haloes that match the desired density profiles.After discarding the halo particles, we introduce the disc component to the dynamically relaxed FDM halo and perform selfconsistent halo-disc simulations as outlined in the flowchart Fig. 3.
The initial axisymmetric stellar discs follow the canonical exponential surface density profile where  d denotes the disc scale radius, and the disc total mass is The initial three-dimensional disc density profile is taken to be (cf.Eq. ( 3)) where ℎ() is the disc scale height profile.In all the six simulated FDM haloes, the discs have the same total mass  d = 3.16×10 9 M ⊙ , scale radius  d = 3.0 kpc, and initially radius-independent scale  1.The de Broglie wavelength  dB , Eq. ( 8), and consequently the FDM density granulation length scale decrease for higher FDM particle mass  22 and larger virial halo mass  h .
height ℎ = 0.15 kpc.Each -body galactic disc is constructed with a total number of  d identical particles; each disc star particle hence weights  d =  d / d .The adopted  d in each production run listed in Table 1 has been numerically verified to have sufficiently high mass resolution for numerical convergence (see Appendix B for details).
The initial disc velocity dispersion profile is obtained by solving the cylindrically symmetric Jeans equations where Φ tot ≡ Φ d + Φ bg denotes the total gravitational potential sourced by both the stellar disc and the background halo.  is the azimuthal velocity of disc stars and related to the total circular velocity defined in Eq. ( 2) via  2 cir =  2  −  2  ; ⟨⟩ is the ensemble average.We let   =   =   when the disc is generated, where   ,   , and   are the vertical, radial, and azimuthal stellar velocity dispersion.The velocity distribution is set to be Gaussian.
Although each constructed disc is initially in equilibrium with the spherical and smoothed background halo potential, the disc configuration can still significantly evolve during the initial disc-halo co-evolution due to following three sources of perturbations: (1) The additional of a disc component causes an oscillatory quadrupolar distortion in the initially spherical host halo potential.(2) At a fixed time slice, FDM haloes are filled with local granular fluctuations that deviate from the assumption of spherical symmetry used in GALIC.
(3) These local density fluctuations are time-varying.The resulting potential perturbations can temporarily destabilise the initial disc configuration.We therefore need to relax the system further before any disc properties can be reliably measured (see Appendix A for details).The disc stability can be quantified by the Toomre  parameter (Toomre 1964) where  denotes the disc epicyclic frequency.The disc is said to be Toomre stable if  > 1 (or unstable if  < 1).Dynamical evolution of FDM haloes and stellar discs are carried out using gamer-2 (Schive et al. 2018), a GPU-accelerated gridbased simulation code with adaptive mesh refinement (AMR).The flowchart of the entire simulation setup is presented in Fig. 3.We define the relaxation time  rel , where  rel = 0 Gyr corresponds to the start of each self-consistent simulation.During this phase of initial relaxation where ( ≤ 15 kpc) ≲ 1, disc spiral arm structures and local density clumps quickly develop (see Fig. A1), causing the disc rotation curve to fluctuate considerably with time.Due to the aforesaid sources of potential perturbations, the disc heating rate in this period is significantly higher than that in the post-relaxation phase, making the robust inference of the granulation-driven disc heating rate difficult.This rapid increase in   () in turn raises the values of () across the entire disc.
After about half a free-fall time 0.5 ff ≃ 1.05 Gyr, the spiral arms and local density clumps become less pronounced; furthermore, both the disc rotation curve and halo density profile stabilise.We thus measure the granulation-driven disc heating only after  rel ≥ 1.4 Gyr when ( ≤ 15 kpc) > 1 is always satisfied, and define In each of the six cases listed in Table 1, we first dynamically co-evolve the GALIC-generated stellar disc and the relaxed FDM halo over  rel = 0-1.4Gyr.At the end of this corelaxation phase that we define as  sim = 0 Gyr, the disc configuration reaches a quasi-equilibrium state; we continue each simulation for another 2.1 Gyr and quantify the granulation-driven disc heating rate.Hence  rel ≡  sim + 1.4 Gyr for all cases.
sim =  rel − 1.4 Gyr for all cases.Each system at  sim = 0 Gyr is regarded to have achieved dynamical quasi-equilibrium, where nongranulation-driven disc heating sources have been minimised.Fig. 4 shows the edge-on (top) and face-on (bottom) views of the relaxed initial condition of the disc (i.e. at  sim = 0 Gyr) hosted by the FDM halo with  22 = 0.4 and  h = 7 × 10 10 M ⊙ .The disc-relaxed initial conditions of all six cases now look very distinguishable, as shown in the leftmost column of Fig. 6.This post-relaxation dischalo configuration is adopted as the suitable initial condition in all our simulation runs, and further details on the relaxation process can be found in Appendix A. For comparison, in Appendix C, we have also simulated and examined a CDM halo-galactic disc system constructed by GALIC, following the same numerical setup (Fig. 3).

AMR criteria for halo-disc simulations
The dynamical evolution of disc-halo systems requires sufficient spatial and temporal resolution over the entire region of interest to prevent artificial disc thickening and outer mass infall of the FDM halo caused by numerical errors.In each simulation run, we adopt separate AMR schemes for the FDM halo and the galactic disc to ensure both components are always adequately resolved.If at a given location two schemes yield different refinement requirements, the more stringent (i.e. more finely resolved) condition is applied.
Since the local de Broglie wavelength  dB Eq. ( 8) and hence the granulation length scale decrease with increasing FDM particle masses  22 and halo velocity dispersion  h , FDM haloes with heavier  22 and/or at smaller radii require higher spatial resolution.The finest grid size is set to resolve the smallest density granulation size by at least 20 cells, and this finest refinement level covers the central halo region.The transition from the finest to the second finest refinement levels is set to occur roughly at 40 c , where  c is the soliton core radius defined in Eq. ( 16).The boundary between the second highest and the third highest refinement levels corresponds to the radius at which the granulation size doubles compared to the ).For all six cases simulated in this work, the relaxed disc configuration is taken as the suitable initial condition after which granulation-driven disc heating rate can be reliably measured.
smallest one across the entire simulation box, which occurs at  ≃ 60-80 kpc for  h = 7 × 10 10 M ⊙ haloes.Lastly, the transition from the third to the fourth refinement level marks the radius at which the granulation size increases fourfold compared to the smallest one.As an example, Fig. 5 shows the FDM halo density slice together with the grid refinement structures for the case  22 = 0.8.The innermost AMR boundary is located at  ≃ 40 c ≃ 20 kpc, and the transition from the second highest to the third highest refinement levels lies at  ≃ 60 kpc.
For the stellar discs, the finest refinement level is set to resolve the initial disc scale height vertically by a total of five to six cells, giving a cell size of 0.05-0.06kpc.In practice, the refinement criteria are determined by the particle number in each cell and dynamically updated to ensure that ≥ 97% of all disc particles are always resolved to this highest level throughout each simulation run.In general, the cell sizes at the highest refinement level differ under the halo and disc refinement criteria.As mentioned earlier, we always adopt the more stringent refinement criteria (smaller cell size).For the six cases listed in Table 1, the finest cell size is determined by the disc (halo) refinement scheme for the cases  22 = 0.2 and 0.4 ( 22 = 0.8 and 1.2).These AMR criteria are tailored towards the finite difference numerical scheme for FDM simulations (see Schive et al. 2014a for details).

Disc property analysis
The prime focus of this work is on the kinematic and morphological evolution of galactic discs with time.During the disc-halo coevolution, the disc dynamical centre traces the trajectory of the central soliton core (Fig. 1) and can differ from the centre of mass of the entire disc.Hence to reliably infer the stellar disc properties at a given time slice, we adopt the instantaneous location of the soliton peak density as the disc centre.With respect to this reference centre, the -axis is taken to be the instantaneous direction of total disc angular momentum to account for any possible halo-disc misalignment.
The disc properties are then computed in this cylindrical coordinate system, with  = 0 as the disc mid-plane.
To compute the ensemble-averaged disc scale height ℎ() and vertical velocity dispersion   , we group all disc particles into distinct annuli evenly spaced in radius .In each radial bin, the scale height ℎ() is defined as the vertical distance such that the ratio of mass enclosed within || ≤ ℎ() and the total mass equals tanh 1 ≃ 0.762, a definition consistent with Eq. ( 20).After projecting all disc particles onto the mid-plane, the vertical velocity dispersion   can be obtained by where   denotes the vertical velocity of individual star particles.

DISC HEATING: SIMULATIONS VS. THEORY
We quantify in Section 4.1 the rate of galactic disc thickening observed in the self-consistent FDM simulations.Section 4.2 examines the level of consistency between simulated and predicted granulationdriven disc heating rates.We discuss in Section 4.3 the applicability of the Fokker-Planck approximation in deriving the disc heating rate by C23.

Simulation results
We first examine the disc morphological evolution in self-consistent FDM simulations.Starting from the suitably relaxed disc-halo initial conditions at  sim = 0 Gyr (see Section 3.2 and Fig. 3), we dynamically co-evolve the disc-halo systems further to  sim = 2.1 Gyr.Fig. 6 shows the edge-on disc density projections of the six cases listed in Table 1 (top to bottom row) at  sim = 0.0, 0.7, 1.4, and 2.1 Gyr (first to fourth column).Disc thickening is observed in all simulation runs.For the same FDM halo mass  h , the level of disc thickening increases for smaller  22 , as larger de Broglie wavelengths Eq. ( 8) give rise to more massive granular structures and thus stronger disc heating (see Section 2 for the precise theoretical formulation).Similarly for a fixed  22 = 0.4, the disc heating rate increases for lighter halo masses  h that exhibit lower  h () and consequently larger de Broglie wavelengths.These overall qualitative trends are carefully assessed below.
The time evolution of azimuthally averaged disc scale heights ℎ(,  sim ) is compared in Fig. 7, showing a systematic increase in ℎ(,  sim ) with time.The wave-like fluctuating feature in all scale height profiles in the outer disc region  ≳ 10 kpc is caused by the persistent spiral structures and local density clumps.For haloes with more compact solitons ( 22 = 0.8, 1.2), the scale height profiles also form sharp peaks around the transition between solitons and haloes, which occurs at ≃ 3.3 c = 1.6 (1.1) kpc for  22 = 0.8 (1.2).To further quantify the disc thickening rate, Fig. 8 plots the ensembleaveraged scale height of four radial bins with Δ = 2 kpc in width centred on  = 4 kpc (blue), 6 (orange), 8 (green), and 10 kpc (red).In all six cases, the disc scale heights increase approximately linearly with time during  sim = 0-2.1 Gyr, which can be easily seen by comparing to the linear best-fit growth curves of the respective  = 6 kpc bin data (dashed lines).Furthermore, the level of disc thickening is sensitive to  22 .For FDM haloes with the same virial mass  h = 7.0 × 10 10 M ⊙ , the averaged disc scale height increases by ≃ 0.8 kpc within 2 Gyr for  22 = 0.2, and only by ≃ 0.08 kpc for  22 = 1.2.For the three haloes with  22 = 0.4 and varying  h , the disc thickening rate is slightly higher for smaller  h .
The granulation-driven heating can be most directly quantified by measuring the disc vertical velocity dispersion   (,  sim ).Fig. 9 shows the time evolution of   ( ≤ 15 kpc,  sim ) from  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue) for all six simulated cases listed in Table 1.As expected from the universal disc thickening seen in Figs. 6 and 7, we observe that disc particles all become kinematically hotter over time.For each individual profile, the disc vertical velocity dispersion increases more rapidly for smaller .Since the halo velocity dispersion  h between 3.3 c and 15 kpc is approximately constant for all the FDM haloes simulated in this work, the radial Although the relaxed initial conditions at  sim = 0 Gyr can non-trivially differ, all galactic discs experience granulation-driven stellar heating.For a fixed halo mass  h , the rate of disc thickening is higher for smaller  22 .For a fixed  22 = 0.4, thickening rates decrease slightly with increasing  h .
dependence of the disc heating rate is primarily sourced by  h ().
Indeed with decreasing , other things being equal, the granulationdriven potential perturbations increase in magnitude due to higher local halo density  h () (Fig. 1).
We next examine the time evolution of vertical velocity squared  2  of disc particles grouped in four 2-kpc-wide radial bins centred on  = 4 kpc (blue), 6 (orange), 8 (green), and 10 kpc (red), as shown in Fig. 10.The simulation data outputs are evenly spaced in time across  sim = 0-2.1 Gyr, and the slope between any two adjacent data points in  2  yields the instantaneous stellar heating rate  2  /.For a fixed value of  h , stellar discs in FDM haloes with smaller  22 experience stronger heating due to relatively enhanced granulation-driven gravitational perturbations.For the three cases with  22 = 0.4, we observe a gentle increase in the overall disc heating rates with decreasing halo masses  h = 0.7-2.8× 10 11 M ⊙ probed here.These trends are consistent with the disc scale height evolution shown in Figs.7 and 8. Across all six cases, the  2  profiles generally exhibit an approximately linear growth.The only exception is the  = 10 kpc bin in the  22 = 0.4 and  h = 2.8 × 10 11 M ⊙ host halo, which has a small bump around  sim = 0.6 Gyr, and its scale height also grows more rapidly compared to other radii as shown in Fig. 8.

Disc heating rates: Simulations vs. Theory
Having discussed the simulation results in Section 4.1 and reviewed the corresponding theoretical framework in Section 2, we present in this subsection the comparative analysis of simulated vs. predicted disc heating rates.To compare directly with the orbit-averaged theoretical predictions of disc vertical heating rate  2  / Eq. ( 12) in each radial bin, we first compute the linear slope Δ 2  /Δ sim of all adjacent data point pairs.Since the heating rate is approximately time-independent over Δ sim = 2.1 Gyr ≫ ( ≤ 10 kpc,  max ) ≃ 0.05-1.0Gyr, the time-averaged slope then gives the orbit-averaged disc heating rate in each radial bin.Fig. 11 compares the ensembleaveraged disc heating rates between numerical simulations (blue solid) and analytical estimates (red solid) 2 on the left -axis; the 2 Since our simulation setup does not include external mass infall (cf. the continuous disc accretion scenario considered in C23), the disc-mass-accretiondriven heating term  ln Σ  in the SGD limit, Eq. ( 12), is not considered when computing the theoretical heating rates here.The disc surface density can still fluctuate over  sim = 0-2.1 Gyr due to radial migration, as shown in Fig. A4.We have explicitly verified that this time-averaged change in Σ could have altered the theoretical ensemble-averaged disc heating rates only by at most MNRAS 000, 1-20 ( 2023

R (kpc)
Figure 7. Disc scale height profiles from  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue).Granulation-driven stellar disc thickening occurs at all radii of interest.The disc thickening rates are noticeably higher in FDM haloes with smaller  22 , fixing  h = 7 × 10 10 M ⊙ .For the cases with  22 = 0.4, the thickening rates are slightly higher in less massive haloes (see also Fig. 8).ratio of theoretical over simulated heating rates (green dashed) is labelled on the right -axis.
The simulation results generally agree well with the theoretical estimates.For  22 = 0.4, analytical estimates agree within 10-50% to the simulation measurements across all radial bins with different  h .However, the theory overpredicts (underestimates) the heating rate for  22 = 0.2 ( 22 ≥ 0, 8) by a factor of ≃ 1.5-2 at all radii (except for the soliton-occupied region for  22 = 0.2, at which the theory overpredicts the heating rate by more than a factor of five).We discuss in Section 4.3 some possible factors contributing to the heating rate discrepancies between theory and simulations.
We furthermore observe that the predicted  22 -and  hdependence of the disc heating rate H ∝  −1 heat ∝  −3 22  −6 h  2 h in Eq. ( 14) is qualitatively consistent with simulations.In essence, smaller  22 with fixed  h corresponds to a larger de Broglie wavelength Eq. ( 8), which gives rise to more massive granular structures Eq. ( 7) and hence higher disc heating rates.Similarly for a fixed  22 , since lighter halo masses  h correspond to smaller local halo density  h (Fig. 1) and velocity dispersion  h , the sharper  h dependence compared to that of  h in the disc heating rate H ∝  −6 h  2 h , implies stronger disc heating in less massive haloes.However, the quantitative heating rate scaling with respect to  22 found in simulations appears to be 'shallower' than the analytical prediction H ∝  −3 22 .Next concerning the -dependence of  2  /, Fig. 11 clearly shows that the disc heating rate increases with decreasing  in all six cases.Since the halo velocity dispersion  h between 3.3 c and 15 kpc is approximately constant for all the FDM haloes simulated in this work, the inner and outer disc regions in a given halo differ ≃ 10% in the outer radial bins for the case  22 = 1.2; this heating source becomes negligible for smaller  and in lighter  22 cases. in 2-kpc-wide radial bins centred on  = 4 kpc (blue), 6 (orange), 8 (green), and 10 kpc (red) over  sim = 0-2.1 Gyr.We observe that  2  increases roughly linear with time, and the granulation-driven disc heating rate is generally higher for smaller  22 and .In contrast, the disc evolved in a CDM halo shows negligible heating   2  / ≲ 1 km 2 s −2 Gyr −1 (Fig. C3).mainly in  h , varying by up to a factor of ≃ 4 between  = 4 kpc and 10 kpc (Fig. 1).Given the scaling relation H ∝  2 h , the larger  h () with decreasing radius is expected to yield higher local disc heating rates, consistent with the general trend in Fig. 11 (see also Figs. 9 and 10).
Having uncovered great simulation-theory consistency in the radial dependence, we now focus on the time dependence of  2  /.Across all six halo-disc systems, the approximate linear growth in  2  (Fig. 10) corresponds to a roughly time-invariant  2  /.As the relaxed halo profiles  h () and  h () remain stable throughout  sim = 0-2.1 Gyr with negligible secular evolution (see Fig. A5), the above observation implies a   -independent disc heating rate, consistent with the analytical prediction Eq. ( 14).To interpret this result in relation to the disc scale height growth via Eq.( 4), we first need to determine when either of the analytical expression ℎ() can be reliably applied.
The relative contributions of the disc self-potential Φ d and the non-disc background Φ bg to the total gravitational potential Φ tot can be effectively quantified by (  ,  eff bg , Σ) defined in Eq. ( 5), where  ≪ 1 (≫ 1) corresponds to the SGD (BGD) limit.Fig. 12 shows the time evolution of  from the end of the co-relaxation phase  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue) for all six cases listed in Table 1.We observe that in the three haloes with  22 = 0.4, the stellar discs are generally in the intermediate regime  ≃ 1.For the case  22 = 0.2 (0.8 and 1.2), the inner disc region remains in the BGD (SGD) limit for  sim = 0-2.1 Gyr.
Provided that the disc surface density Σ does not change significantly with time (Fig. A4), the analytical disc scale height Eq. (4) in the SGD limit ℎ ∝  2  naturally accounts for the approximate linear growth in disc scale height shown in Fig. 8. Contrastingly in the  / from self-consistent simulations (blue) and theoretical predictions (red) at  = 4, 6, 8, and 10 kpc.The ratio of theoretical to simulated heating rates (green) is shown on the right -axes.Outside of the soliton core  ≲ 3.3 c (gray shaded), theoretical estimates agree within about a factor of two to simulations at all radii for all six cases.In the three FDM haloes with  22 = 0.4, analytical estimates are highly consistent with the measured heating rates, agreeing within 10-50% at all radii.Interestingly, the theory overpredicts (underestimates) the heating rate for the case  22 = 0.2 (cases  22 = 0.8 and 1.2).The Fokker-Planck approximation breaks down and results in significant overestimation in the soliton-occupied region.We discuss the possible factors that contribute to the heating rate discrepancies in Section 4.3.
BGD limit (applicable for the case  22 = 0.2 at  = 4 kpc), Eq. ( 4) reduces to ℎ ∝   ∝  0.5 , suggesting a gentle decrease in the disc height growth rate ℎ/ with time.This mild flattening in ℎ( sim ) can indeed be identified in Fig. 8 (blue curve in the top-left panel) towards  sim ≃ 2.1 Gyr, but a unequivocal corroboration will require additional simulation data extending beyond  sim ≥ 2.1 Gyr.
Overall, the comparative analysis reveals a general agreement between the simulation measurements and the theoretical estimates.Nonetheless, a number of non-trivial quantitative differences have also been identified in the preceding discussion.In Section 4.3, we examine in detail the possible causes of this discrepancy and the validity of the Fokker-Planck approximation assumed in the analytical heating rate estimate of C23.

Applicability of the Fokker-Planck approximation
Here we identify major factors that possibly contribute to the heating rate discrepancies between simulations and theory observed in Fig. 11: • The Effect of Soliton: The significant mismatch at inner radii for the case  22 = 0.2 is expected and attributed to the proximity to the central soliton core, given that the soliton-halo profile transition occurs at roughly 3.3 c ≃ 4.4 kpc as shown in Fig. 1.The Fokker-Planck formalism (Section 2) assumes that all background perturbers can be treated as statistically uncorrelated and spatially unconfined quasi-particles, which is appropriate only for modelling a large number of FDM density granules  gra ≫ 1. How- Figure 12.  profiles from  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue).As defined in Eq. ( 5),  provides an effective measure of the relative importance of the disc self-gravity and non-disc background potential such that  ≪ 1 corresponds to the self-gravity dominated (SGD) limit, while  ≫ 1 gives the background-dominated (BGD) limit.For the case  22 = 0.2, the disc lies almost entirely within the BGD limit.For the cases  22 = 0.4 ( 22 = 0.8 and 1.2), the systems are generally in the intermediate regime (SGD limit).Since analytical heating rate estimates exist only in the SGD and BGD limits (see Section 2), in the intermediate regime 0.5 <  < 1.5 we linearly interpolate theoretical heating rates from both limits; the accuracy of such an approach is discussed in Section 4.3.
ever, this assumption undoubtedly breaks down at sufficiently small radii  ≤ O (3.3 c ) where the gravitational potential perturbations are dominantly sourced by the soliton.As the ground-state solution to the governing Schrödinger-Poisson equations, this coherent soliton core executes confined random-walk like excursions around the host halo centre of mass (Schive et al. 2020;Dutta Chowdhury et al. 2021).The theory hence overpredicts the disc heating rate by 'mistakenly treating' the gravitational interactions with a single soliton as repeated encounters with a group of uncorrelated granules.A similar stellar heating rate overestimate near the central soliton core is also reported in Dutta Chowdhury et al. (2021), where the discrepancy was partially alleviated by introducing effective (suppressed) diffusion coefficients within  ≤ 2.3 c .
• Validity of Molecular Chaos in the Quasi-particle Treatment: The roughly factor-two theoretical heating rate overestimate in the outer three radial bins for the case  22 = 0.2 can partly result from the limited applicability of the molecular chaos assumption 3 .Here the granule effective radius is  gra ≃ 2.2 kpc ≳ ℎ(,  sim ) for 3 More precisely, in deriving the analytical expressions of the diffusion coefficient Eq. ( 10), Bar-Or et al. (2019) assumed that the characteristic correlation time of the background perturbers is negligibly small compared to all other dynamical timescales of interest.In the language of classical two-body relaxation, it is equivalent to the assumption that the velocities of colliding background particles are uncorrelated, and consequently the distribution functions of the subject particles and background perturbations are statistically independent (i.e. the assumption of molecular chaos). / 2 h of the disc star vertical motion to the FDM halo from  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue).Except for the soliton-enclosed region (grey shaded) for the case  22 = 0.2, the stellar discs all remain kinematically cold  2  / 2 h ≲ 1 compared to the respective host haloes.
sim = 0-2.1 Gyr (see Figs. 7 and 8), and the total number of FDM quasi-particles within 3.3 c ≤  ≤ 15 kpc and || ≤ ℎ should be less than  gra ≲ 35.The small number count  gra indicates that these density granules can factually be collisional, and consequently the velocities of colliding quasi-particles are no longer statistically uncorrelated.In addition, the cumulative energy transfer of these correlated potential perturbations can excite non-trivial bulk motion of disc stars, as opposed to all being directed to increasing the disc star random motion  2  .The analytical framework that treats FDM density granulation as uncorrelated quasi-particles thus expects to be less accurate for sufficiently small  gra and overestimate the genuine disc heating rate  2  /.Recently, Zupancic & Widrow (2023) explored the applicability of quasi-particle approximation in lieu of the fully self-consistent Schrödinger-Poisson dynamics.Their simulation results, albeit in one spatial dimension, suggest that these two treatments are consistent within the first ∼50 dynamical times of evolution, after which the quasi-particle formalism can non-trivially overpredict the stellar heating rates.However, whether this observation holds in threedimensional systems requires further tailored simulations and is beyond the scope of this work.
• Coulomb Logarithm: Another possible cause of the heating rate discrepancy at  22 = 0.2 lies in the ambiguous nature of the Coulomb logarithm definition ln Λ, which can be safely ignored only for a sufficiently large Coulomb factor Λ ≡  max / min ≫ 1 (e.g.Just & Peñarrubia 2005;Binney & Tremaine 2008).For  22 = 0.2, ln Λ ≃ 1.5-3.0within  ≤ 12 kpc is sufficiently small such that any physical corrections to ln Λ could sizeably impact the predicted heating rate.
• Heating Behaviour across the SGD and BGD Limits: Analytical expressions for the granulation-driven ensemble-averaged disc heating rates Eq. ( 12) exist only if the disc self-gravity either dominates (the SGD limit) or can be ignored relative to the background potential (the BGD limit).For a galactic disc with time-independent Σ(), the theoretical heating rate in the BGD limit is 1.5-2 times higher than that in the SGD limit.In the intermediate regime, where the heating rate cannot be computed directly, we have estimated the analytical heating rate by linearly interpolating the heating rates in both limits as discussed in Section 2. This comparatively uncertain regime applies to the three cases with  22 = 0.4 (Fig. 12), which could partly contribute to the ≲ 10-50% discrepancy between predicted heating rates and simulation measurements observed in Fig. 11.
Outside of this intermediate region 0.5 ≤  ≤ 1.5, however, comparatively greater heating rate overestimates (underestimates) are observed for the case  22 = 0.2 (0.8 and 1.2) in the BGD (SGD) limit.In all six disc-halo systems, there appears to a correlation between BGD/SGD limits and overestimate/underestimate of the disc heating rates.Furthermore in the SGD regime, the factor-two simulation-theory discrepancies in the cases  22 ≥ 0.8 can hardly be accounted for by the aforesaid uncertainties.It thus remains possible that the SGD limit itself has some unknown errors that leads to this non-trivial heating rate underestimation.Future simulations that explore the range  22 ≥ 1.2 would help clarify this uncertainty.
• Numerical (Artificial) Heating: As the granulation-driven disc heating decreases rapidly with increasing  22 , other sources of disc heating might become discernible in the simulation data.We assess the level of numerical (artificial) disc heating by conducting a disc-CDM-halo simulation run with  h = 7 × 10 10 M ⊙ (see Appendix C for more details).The heating rates in all four radial bins are below ≲ 1 km 2 s −2 Gyr −1 , which is negligible or sub-dominant compared to the FDM-driven heating rates as seen in Fig. 11.Hence under the sufficient particle and force resolutions adopted in this work (see Section 3.3 and Appendix B), we conclude that artificial heating is unimportant and can be safely excluded in our analysis.
• Evolution in the Disc Kinematic Temperature: In the Fokker-Planck formalism (Section 2), only perturbers travelling faster than the disc star particles contribute to diffusive heating (e.g.Binney & Tremaine 2008).Since the velocity distribution of FDM quasi-particles is assumed to be Maxwellian, the population of density granules capable of contributing to the disc heating shrinks as a stellar disc becomes kinematically hotter.Namely, the granulationdriven heating becomes ineffective if the stellar disc is kinematically hot relative to the host halo  2  / 2 h ≥ O (1).We examine the time evolution of  2  / 2 h from  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue) for all six disc-halo systems in Fig. 13.Except for the solitonoccupied regions (gray shaded), the stellar discs all remain kinematically cold relative to the respective host haloes.The disc heating rates hence should not be meaningfully impacted by the increase in  2  / 2 h , consistent with the nearly time-independent growth rate of   observed in Fig. 10 for  sim = 0-2.1 Gyr in all six cases.
• Other Physical Disc Heating Sources: The theoretical underestimate in disc heating rate for the two cases  22 = 0.8, 1.2 likely indicates the presence of other physical, non-granulation-driven heating mechanisms such as the outward radial migration (e.g.Sellwood & Binney 2002;Minchev & Famaey 2010), transient spiral modes (e.g.Barbanis & Woltjer 1967;Minchev & Quillen 2006), and/or large-scale bending instability (e.g.Rodionov & Sotnikova 2013).For these two galactic discs, orbital diffusion proceeds at a rate 10 −3 -10 −2  d per Gyr in all 2-kpc-wide radial bins, such that radial migration alone cannot account for the factor-two mismatch in predicted and observed disc heating rates (see also Figs. 9 and A4).Although quantifying these non-FDM disc heating mechanisms is beyond the scope of this work, it is worth noting that the level of theory-simulation discrepancy still remains within a factor of two even when the granulation-driven disc heating rate reaches as low as  2  / ≃ 3 km 2 s −2 Gyr −1 for the case  22 = 1.2 and  h = 7 × 10 M ⊙ .

IMPLICATIONS FOR EXTERNAL DISC GALAXIES
Section 5.1 compares the analytical disc scale height expressions with simulation measurements.Uncertainties in the hydrostaticequilibrium-based   inferences are discussed in Section 5.2.

Instantaneous disc scale height inference
Observationally, the fitted vertical structures of galactic disc density profiles are commonly assumed to follow either ∝ sech 2 (/ℎ) (e.g.Comerón et al. 2011) or ∝  − 2 /ℎ 2 (e.g.Dutta et al. 2009).However, even for locally isothermal discs in hydrostatic equilibrium, these two analytical scaling relations Eq. ( 3) are exact only if either the non-disc background potential or the disc self-gravity can be ignored relative to the other component (C23).In this subsection, we assess under what necessary conditions do the analytical formulae accurately predict the measured disc scale heights, as well as quantify the level of discrepancy when the conditions are not met (i.e. the total gravitational potential Φ tot is comparably sourced by the disc itself Φ d and the non-disc background Φ bg ).
As discussed in Section 4.2, the dimensionless number (  ,  eff bg , Σ) defined in Eq. ( 5) provides a quantifiable measure for the relative contributions of Φ d and Φ bg to Φ tot , where  ≪ 1 (≫ 1) corresponds to the SGD (BGD) limit.Prior to the disc-halo co-relaxation at  rel = 0 Gyr, the initial disc configurations (see Section 3.2) are all kinematically cold with  ranges from ≃ 0.25 at  = 4 kpc to 0.4-0.5 at  = 10 kpc.As each galactic disc thickens with increasing   over time, the effect of background potential can gradually dominate over the disc self-gravity in shaping individual disc stars' vertical oscillation motion Eq. ( 6), possibly leading to a transition from the SGD limit to the BGD limit.The post-corelaxation time evolution of each  profile is shown in Fig. 12 over  sim = 0 Gyr (light blue; equivalent to  rel = 1.4 Gyr) and 2.1 Gyr (dark blue).By  sim = 2.1 Gyr, the galactic discs are generally in the intermediate regime  ≃ 1, except for the inner disc region of the case  22 = 0.2 (0.8 and 1.2) lying in the BGD (SGD) limit.
Focusing on the following three cases with  22 = 0.2, 0.4, 1.2 (blue, orange, and green respectively) and  h = 7.0 × 10 10 M ⊙ , Fig. 14 compares the disc scale heights in the radial bins centred on  = 4 kpc (top panel) and 10 kpc (bottom panel), plotting values from either simulation measurements (solid) or analytical solutions (dashed, Eq. ( 4)).The latter is defined as the minimum of the two limits ℎ() = min(ℎ SGD , ℎ BGD ) 4 , since there exists no closed-4 The adoption of theoretical disc heights as the smaller analytical value in the two limits ℎ() = min(ℎ SGD , ℎ BGD ) is an (overly) optimistic choice that minimises the theory-simulation discrepancy.However even under such a 'best-case' scenario, the simulation(observation)-theory mismatch in the disc scale height inferences can still approach a factor of two as shown in Fig. 14, highlighting the need to explicitly verify the validity of Eq. ( 4) by carefully assessing the relative importance between Φ d and Φ bg when modelling disc vertical structures in external galaxies.If one instead assumes either of the analytical relation Eq. ( 4) across all  ≥ 0, the level of scale height overestimation is actually unbounded above and can exceed a factor of two.4)) at  = 4 kpc (upper panel) and 10 kpc (lower panel) for the cases  22 = 0.2 (blue), 0.4 (orange), and 1.2 (green) over  rel = 0-3.5 Gyr.Note that  sim ≡  rel − 1.4 Gyr.Analytical expressions generally overestimate the scale height, by up to a factor of two, and are sufficiently accurate only when the disc is properly in the SGD or BGD limit, as shown in Fig. 12; see Section 5.1 for details.
form expression for ℎ in the intermediate regime  ≃ 1. Across the entire simulation time span  rel = 0-3.5 Gyr, we observe that analytical estimates and simulation results agree well in the SGD limit at the onset of disc-halo co-relaxation  sim ≲ 0.25 Gyr.This general agreement remains for the case  22 = 1.2 at  = 4 kpc where the SGD limit  ≤ 0.5 remains applicable until  rel = 3.5 Gyr (see the top-right panel of Fig. 12 at  sim = 2.1 Gyr).At the other extreme where the BGD limit clearly applies, for the case  22 = 0.2 at  = 4 kpc, the theoretical and simulation-inferred scale heights differ by less than ≃ 15% for  rel ≳ 2 Gyr.In the intermediate regime  ≃ 1, however, adopting either the sech 2 (/ℎ) or  − 2 /ℎ 2 scale height analytical solutions 'undercounts' the total vertical gravity by about half.For the three cases with  22 = 0.4, Fig. 14 shows that the analytical predictions can overestimate the scale height by up to a factor of two.We discuss the relevant implications to observational inferences of disc properties in Section 5.2.

Observational uncertainties of disc scale height and 𝜎 𝑧
In external disc galaxies, radial profiles of stellar/gas disc vertical velocity dispersion   for edge-on (or scale height ℎ for face-on) systems cannot be directly measured.With observationally constrained baryonic surface density profile Σ() and/or effective background density  eff bg , it is customary to assume the galactic disc in question to be either self-gravitating (the SGD limit) or dictated by the non-disc background potential (the BGD limit) such that analytical hydrostatic-equilibrium solutions for isothermal discs exist and can be directly applied (e.g.van der Kruit & de Grĳs 1999; Leroy et al. 2008).The observationally inaccessible   () or ℎ() can then be directly inferred from either expression of Eq. ( 4), depending on which assumption is adopted (e.g.Kregel et al. 2002;Kasparova & Zasov 2008;Patra 2019;Das et al. 2020).However, in the intermediate regime 0.5 <  < 1.5 as commonly found in our simulations (Fig. 12), both the disc and non-disc background contribute comparably to the local total vertical gravity.By explicitly assuming the SGD (BGD) limit, the total vertical gravity is always undercounted by accounting for only the (non-)disc component.
The effect of this vertical gravity undercounting can be substantial in the inferred disc properties, as observed in our self-consistent disc-halo simulations.Analytical solutions Eq. ( 4) can overestimate the ensemble-averaged disc scale heights by up to a factor of two, as demonstrated in Fig. 14.In the intermediate regime 0.5 <  < 1.5 where both the SGD and BGD limits fail, 'blindingly' applying Eq. ( 4) can result in   () underestimated by up to a factor of √ 2 (or 2) in the SGD (BGD) limit.
This uncertainty applies to all external disc galaxies where the baryonic vertical velocity dispersion is indirectly inferred from ℎ() and Σ().That being said, the caveat is less of a concern when either the SGD limit (e.g. for ultra-thin discs studied by Matthews (2000) and Bizyaev et al. (2017)) or the BGD limit (e.g.stellar discs that are kinematically hot relative to the host haloes, as for the case  22 = 0.2 in Figs. 12 and 14) is guaranteed to hold.To partially address the uncertainty in this type of hydrostatic-equilibrium-based   inferences in observed external galaxies, estimating the value of  could serve as a self-consistency cross-check by requiring  ≲ 0.5 (≳ 1.5) if the SGD (BGD) limit is adopted.This source of inference error is expected to be important when  ≃ 1.

SUMMARY AND CONCLUSIONS
In this work, we numerically quantify the stellar disc heating rates caused by FDM halo density granulation by performing the first self-consistent simulations of FDM halo and -body stellar disc, spanning   = 0.2-1.2× 10 −22 eV (equivalently  22 = 0.2-1.2) and halo virial masses  h = 0.7-2.8× 10 11 M ⊙ as listed in Table 1.The GALIC-constructed galactic discs weight  d = 3.16 × 10 9 M ⊙ and are resolved with  d = 0.8-1.6 × 10 8 equal-mass disc particles; prior to initial disc-halo co-relaxation, these discs have the same scale radius  d = 3.0 kpc and radius-independent scale height ℎ = 0.15 kpc.The main results are summarised as follows: • Disc thickening is observed in all disc-halo systems (Fig. 6), at rates ℎ/ ≃ 0.04-0.4kpc/Gyr (Figs. 7 and 8) increasing with smaller  22 and  h .The measured disc heating rates  2  / ≃ 4-150 km 2 /sec 2 (Figs. 9, 10, and 11) exhibit the same trend.Namely, for a fixed  h , discs hosted by haloes with smaller  22 have higher heating rates; lighter  h with fixed  22 exhibits higher heating rates.In each individual disc-halo system, the ensemble-averaged disc heating rates decrease monotonically with increasing radius (Fig. 11).Overall, the disc scale height ℎ and vertical velocity dispersion  2  increase approximately linearly with time after the halo-disc corelaxation (Figs. 8 and 10).
• The FDM granulation-driven disc heating rates quantified in simulations are compared directly with the Fokker-Planck-based analytical estimates of C23 (Section 4.2).The simulation measurements of the disc heating rates agree within a factor of two to the theoreti-cally predicted values for all six cases examined in this work (Fig. 11), except for the region within the soliton core  ≲ 3.3 c .For the three FDM haloes of  22 = 0.4, the theory-simulation discrepancy is less than 10-50% at all radii of interest.
• In individual cases, we first note that within the soliton core  ≲ 3.3 c , the Fokker-Planck formalism is invalid (Section 4.3) and expectedly yields a significant overestimation in the innermost radial bin  = 4 kpc of the case  22 = 0.2 (Fig. 11).At outer radii, the analytical prediction for  22 = 0.2 overestimates the disc heating rate within a factor of two.For the three cases of  22 = 0.4, we observe high level of theory-simulation consistency at all radii of interest.For the two cases  22 = 0.8 and 1.2, the theory generally underestimates the genuine disc heating rates (within a factor of two), suggesting that there could be additional heating sources not accounted for in the granulation-driven heating model.We have carefully verified that possible artificial heating due to inadequate numerical and/or improper disc-halo initial conditions is negligible in our FDM simulations by verifying numerical convergence (Appendix B) and comparing to a CDM simulation (Appendix C).
• The disc thickening and heating rates are positively correlated to the granule effective mass  gra ∝  −3 22  −3 h  h as defined in Eq. ( 7).This trend is confirmed across all three distinct regimes:  < 0.5 (the SGD limit, where the non-disc background is negligible), 0.5 <  < 1.5 (intermediate regime where the disc and non-disc background contribute comparably to the total vertical gravity), and  > 1.5 (the BGD limit, where the disc self-gravity is negligible).However, the heating rate scaling with  22 found in simulations appears to be 'shallower' than the analytical prediction H ∝  −3 22 , Eq. ( 14); see Section 4.2.
• Concerning the usual adoption of analytical disc scale heights Eq. ( 4) in modelling the observed external disc galaxies, these closedform solutions exist only in either the SGD or BGD limit and always undercount the total vertical gravity, as discussed in Section 5.1.As a result, the analytical predictions are higher than the true disc scale heights measured in self-consistent simulations (Fig. 14) by up to a factor of two.This implies that hydrostatic-equilibrium-based   inferences in observed external disc galaxies can be underestimated by up to a factor of √ 2 (or 2) in the SGD (BGD) limit; see Section 5.2.
In C23, a conservative exclusion bound  22 ≳ 0.4 was derived from the Galactic disc kinematics by requiring that predicted granulation-driven heating in the solar neighbourhood cannot exceed the observed ensemble-averaged disc velocity dispersion   ( ⊙ ) ≃ 22 km s −1 (Sharma et al. 2021).The fact that the theoretical disc heating rates are highly consistent with simulation data at  22 = 0.4 further supports the robustness of this FDM particle mass constraint.Relatedly, since the disc heating rates in simulations appear to decline less rapidly than the theoretical prediction H ∝  −3 22 with increasing  22 , we expect that the reported range  22 ≃ 0.5-0.7 in C23 favoured by the observed thick disc kinematics could be corrected upwards to be closer to  22 ≃ 1.0.We leave the careful cross-check of this prediction by performing self-consistent -body disc simulations in an MW-sized FDM halo to a future work.
On the other hand, the rough analytical calculations by Church et al. ( 2019) yield a disc heating rate larger by a factor of ∼ 360 than the estimates of C23 (see Sec. 4.2.4 therein), under the same disc and FDM halo parameters.To leading order, the discrepancy in H ∝  2 gra mainly stems from their overestimate of granulation mass  gra by O (10) and the simplified assumption that a disc is everywhere self-gravity dominated (SGD limit).Regarding the MW thick disc formation, Church et al. (2019) adopted an enormously large one-dimensional halo velocity dispersion  h = 200 km s −15 .Due to the strong dependence of granulation-driven disc heating rate H ∝  −3 22  −6 h  2 h on FDM parameters, their incorrectly quoted  h significantly reduced the predicted Galactic disc heating rate by a factor of ≃ 160, incidentally in the right direction to compensate for their theoretical heating rate overestimate (by a factor of ∼ 360 when adopting the more accurate halo velocity dispersion  h ≃ 86 km s −1 ).
Given the within-factor-two agreement between the heating rate predictions of C23 and the measurements in all six self-consistent disc-halo simulations performed in this work, we conclude that Church et al. (2019) overestimated the genuine granulation-driven disc heating rate by at least two orders of magnitude.To compare with observed Galactic disc kinematics, they arrived at a similar exclusion bound  22 ≳ 0.6 by requiring that the predicted disc heating over 12 Gyr cannot exceed a less up-to-date value of   ( ⊙ ) ≃ 32 km s −1 (Binney 2010), inferred indirectly from the best-fit analytical distribution functions to the survey data collected by Gilmore & Reid (1983).Although Church et al. (2019) still arrived at a  22 constraint reasonably close to that of C23, it should be stressed that their exclusion bound was derived on the bases of rather inaccurate assumptions and estimates.
Lastly, it is worth emphasising that the present work aims to quantify the FDM-granulation-driven disc heating rates in self-consistent FDM-baryon simulations, providing an independent cross-check against the analytical heating rate estimates.As discussed in Section 3.2, we adopt the kinematically cold, thin discs appropriate in the ΛCDM cosmology as the initial conditions to perform further FDM halo-stellar disc co-relaxation.It remains to be investigated whether such cold and thin stellar discs could naturally form in the first place in an FDM cosmology.Existing FDM cosmological simulations with baryonic feedback (Mocz et al. 2019;Kulkarni et al. 2022) still lack the mass resolution and large-sample statistics at sufficiently low redshifts to address this question.That being said, as the granulation-sourced disc heating rates depend primarily on the FDM halo attributes and are comparatively insensitive to the instantaneous disc properties (Figs. 10,11), we expect the general conclusions of this work to hold irrespective of the adopted disc initial conditions.
Deciphering the thick disc formation processes in the MW and nearby external disc galaxies is indispensable for gaining a more complete picture of galaxy evolution.Having examined analytically and corroborated numerically the granulation-driven disc heating rates in the MW (C23) and sub-MW mass haloes in this work, we argue that, in an FDM cosmology, this heating mechanism provides a robust and ubiquitous thick disc formation pathway, comparatively insensitive to the detailed disc morphology and assembly history as required in some other proposed thick disc formation models.Beyond the Local Universe, the recently discovered high-redshift disc galaxies by JWST (Nelson et al. 2023) also present exciting opportunities to further put to test various proposed thick disc formation models.Future works could perform a detailed analysis of FDM granulationdriven heating over a larger sample of observed disc-halo systems to place a more stringent exclusion bound on  22 , or examine the quantitative (dis)agreement between the Gaia phase spiral and the out-of-equilibrium features caused by FDM density granulation in MW-sized haloes. (,  rel ) in 2-kpc-wide radial bins centred on  = 4 kpc (blue), 6 (orange), 8 (green), and 10 kpc (red) over  rel = 0-3.5 Gyr, where the stellar discs are sampled by a total number of  d = 2 × 10 7 (dotted), 8 × 10 7 (dashed), and 1.6 × 10 8 (solid) equal-mass particles.The host FDM halo has  22 = 0.2 and  h = 7 × 10 10 M ⊙ , and the maximum spatial resolution is 0.062 kpc (see Table 1).Numerical convergence is achieved among all three runs with varying particle resolutions  d ≥ 2 × 10 7 .
d = 2.0 × 10 7 (dotted), 8.0 × 10 7 (dashed), and 1.6 × 10 8 (solid), evolved in the FDM host halo with  22 = 0.2 and  h = 7×10 10 M ⊙ as an example.During the initial 0.5 Gyr co-relaxation as discussed in Appendix A, instability-driven rapid growth in  2  is observed in all radial bins across the three resolution-varying runs.The postco-relaxation heating curves  rel ≥ 1.4 Gyr from the highest to the lowest particle resolution runs agree within ±3 km 2 s −2 Gyr −1 (≲ 6% at outer radii and ≲ 1% at inner radii).This high level of consistency implies that the numerical convergence is achieved for  d ≥ 2.0 × 10 7 .
How the distribution of disc particle vertical velocities evolves with time provides an additional cross-check for numerical convergence.For an initially Gaussian velocity distribution centred at   = 0 km/s (see Section 3.2), the net contribution of FDM-granulation-driven stellar heating is diffusive and purely second-order.Hence the local velocity distribution should still remain Gaussian while undergoing continuous broadening with time, provided that the adopted particle resolution is sufficient to resolve the disc dynamics.Fig. B2 shows the time evolution of disc velocity probability density distribution in a 2kpc-wide radial bin centred on  = 6 kpc (i.e.including disc particles within  = 5-7 kpc); the disc is resolved with  d = 1.6 × 10 8 equalmass particles and evolved in the FDM halo with  22 = 0.2 and  h = 7 × 10 10 M ⊙ .For the disc particle with velocities within 3  , the simulation data (blue) are consistent with the best-fit Gaussian distributions (red), as the full width at half maximum increases with time due to the continuous granulation-driven diffusive heating.We have verified that this consistency condition is satisfied at all radii of interest under this particle resolution.
To examine the impact of particle resolution on the evolved vertical velocity probability density distribution, we compare in Fig. B3 the probability density distributions in 2-kpc-wide radial bins centred on  = 4 kpc, 6, 8, and 10 kpc at  rel = 3.5 Gyr, for stellar discs resolved with  d = 2 × 10 7 (yellow), 8 × 10 7 (green), and 1.6 × 10 8 (blue; identical to Fig. B2) equal-mass particles.The FDM host halo has  22 = 0.2 and  h = 7 × 10 10 M ⊙ .The Gaussian probability density distributions for all three runs are well preserved until the end of each simulation.Contrastingly in simulation runs with inadequate particle resolution  d (not shown here), we have observed skewed distributions that noticeably deviate from the best-fit Gaussian distributions (red).The presence of such distorted non-Gaussian features in the disc vertical velocity distribution allows one to promptly identify insufficient particle resolution adopted to simulate disc dynamics.
Lastly, at a fixed particle resolution, we have also performed simulations with different maximum spatial resolutions.For the FDM halo with  22 = 0.4 and  h = 7 × 10 10 M ⊙ and the stellar disc resolved with  d = 8.0 × 10 7 equal-mass particles, the two simulation runs with respective maximum spatial resolutions of 120 pc and 60 pc (adopted in the production run) yield ensemble-averaged disc vertical heating rates consistent within ≤ 10% at all radii of interest.We thus conclude that a maximum spatial resolution of 60 pc is sufficient in this halo-disc setup.

APPENDIX C: HEATING FROM CDM HALO
To better disentangle different perturbative sources present during the initial disc-halo co-relaxation (Section 3.2) and further assess the degree of possible artificial heating in our FDM halo-disc simulations (Section 4.3), we perform an additional simulation replacing the FDM halo with CDM using the same code gamer-2.The GALIC-generated spherical -body halo has a virial mass of  h = 7 × 10 10 M ⊙ and an NFW concentration parameter of  h = 8 (see Section 3.2 for more details).The stellar disc initial condition is identical to that in all FDM simulation runs; the density profile follows Eq. ( 20) with  d = 3.16 × 10 9 M ⊙ ,  d = 3.0 kpc, and ℎ = 0.15 kpc.The CDM halo and the stellar disc are respectively sampled with  h = 1.6 × 10 8 and  d = 8.0 × 10 7 equal-mass particles.The highest spatial resolution is 0.05 kpc. rel ≥ 1.4 Gyr, low relative differences indicate that both profiles remain approximately constant with time.Following the same setup as the FDM halo-disc simulations (Fig. 3), we define this time slice as the end of disc-halo co-relaxation  sim =  rel − 1.4 Gyr = 0 Gyr.Next we quantify the ensemble-averaged disc vertical heating rates over  sim = 0-1.5 Gyr, identical to the procedure outlined in Section 4.2.We observe in Fig. C3 that disc vertical heating rates in the CDM halo  2  / ≲ 1 km 2 s −2 Gyr −1 are negligibly small compared with the values observed in FDM simulations, as shown in Figs. 10 and 11.This result further confirms that the comparatively low heating rates observed in the two cases  22 = 0.8 and 1.2 (e.g.see Fig. 11) cannot not be meaningfully attributed to numerical noises, and the observed theory-simulation discrepancy is physical in origin.
As a concluding remark, the perturbations sourced by the CDM and FDM haloes are physically distinct.The former stems from the inconsistency in the initial condition, where the initial CDM halo is constructed with the disc and its density profile undergoes little changes during the simulation.The resulting disc heating saturates after  rel = 1.4 Gyr and is generally contained within the central  ≲ 10 kpc.In contrast, stellar discs in FDM haloes exhibit comparatively richer spiral arm structures (Fig. A1).The gravitational perturbations post-relaxation are dominantly sourced by the FDM granulation and show sustained disc heating at all radii (e.g.Fig. 6).
This paper has been typeset from a T E X/L A T E X file prepared by the author. in radial bins centred on  = 4 kpc (blue), 6 (orange), 8 (green), and 10 kpc (red) over  sim = 0-1.5 Gyr in a GALICconstructed CDM halo with  h = 7 × 10 10 M ⊙ and  h = 8.The measured heating rates   2  / ≲ 1 km 2 s −2 Gyr −1 at all radii are negligible compared to the granulation-driven disc heating in all six FDM simulations carried out in this work (cf. Figs. 10 and 11).We thus conclude that numerical heating is unimportant in our simulation setup.

Figure 1 .
Figure1.Initial shell-averaged density profiles of all six FDM haloes listed in Table1(solid), compared with an NFW profile with  h = 7.0 × 10 10 M ⊙ and  h = 8 (dashed).For a fixed halo mass  h , the central soliton core exhibits higher peak density and smaller core radius  c for heavier FDM particle mass  22 .Beyond the characteristic density profile transition at ≃ 3.3 c , the halo density distribution roughly follows the NFW profile.

Figure 2 .
Figure 2. Density slices of relaxed haloes at  rel = 0 Gyr for all six cases listed in Table1.The de Broglie wavelength  dB , Eq. (8), and consequently the FDM density granulation length scale decrease for higher FDM particle mass  22 and larger virial halo mass  h .

Figure 3 .
Figure 3. Flowchart of the disc-halo simulation setup.In each of the six cases listed in Table1, we first dynamically co-evolve the GALIC-generated stellar disc and the relaxed FDM halo over  rel = 0-1.4Gyr.At the end of this corelaxation phase that we define as  sim = 0 Gyr, the disc configuration reaches a quasi-equilibrium state; we continue each simulation for another 2.1 Gyr and quantify the granulation-driven disc heating rate.Hence  rel ≡  sim + 1.4 Gyr for all cases.

)Figure 4 .
Figure 4. Edge-on (top) and face-on (bottom) projections of the dynamically co-relaxed stellar disc density at  sim = 0 Gyr for the case  22 = 0.4 and  h = 7 × 10 10 M ⊙ .At this stage, the disc rotation curve and surface density profiles have largely stabilised against transient structures (see also Appendix A and Fig.A1).For all six cases simulated in this work, the relaxed disc configuration is taken as the suitable initial condition after which granulation-driven disc heating rate can be reliably measured.

)Figure 5 .
Figure 5.One quarter of a density slice for the halo with  22 = 0.8,  h = 7 × 10 10 M ⊙ (upper right panel of Fig. 2) with AMR grid overlaid.Under the disc-halo AMR scheme described in Section 3.3, the highest refinement level should resolve both the pre-relaxed stellar disc vertically by at least five cells (0.05-0.06 kpc in cell size) and the minimum granulation length scale by at least 20 cells.Here, the innermost AMR boundary occurs at the radius  ≃ 40 c ≃ 20 kpc.The second (third) AMR boundary transition, where the granulation length scale increases twofold (fourfold) compared to the global minimum, is located at  ≃ 60 kpc (≃ 100 kpc).

Figure 6 .
Figure6.Edge-on disc density projections of all six cases listed in Table1(top to bottom row) at  sim = 0.0, 0.7, 1.4, and 2.1 Gyr (first to fourth column).Although the relaxed initial conditions at  sim = 0 Gyr can non-trivially differ, all galactic discs experience granulation-driven stellar heating.For a fixed halo mass  h , the rate of disc thickening is higher for smaller  22 .For a fixed  22 = 0.4, thickening rates decrease slightly with increasing  h .

Figure 11 .
Figure 11.Ensemble-and time-averaged disc heating rates   2 / from self-consistent simulations (blue) and theoretical predictions (red) at  = 4, 6, 8, and 10 kpc.The ratio of theoretical to simulated heating rates (green) is shown on the right -axes.Outside of the soliton core  ≲ 3.3 c (gray shaded), theoretical estimates agree within about a factor of two to simulations at all radii for all six cases.In the three FDM haloes with  22 = 0.4, analytical estimates are highly consistent with the measured heating rates, agreeing within 10-50% at all radii.Interestingly, the theory overpredicts (underestimates) the heating rate for the case  22 = 0.2 (cases  22 = 0.8 and 1.2).The Fokker-Planck approximation breaks down and results in significant overestimation in the soliton-occupied region.We discuss the possible factors that contribute to the heating rate discrepancies in Section 4.3.

Figure 13 .
Figure13.Velocity dispersion squared ratios  2  / 2 h of the disc star vertical motion to the FDM halo from  sim = 0 Gyr (light blue) to 2.1 Gyr (dark blue).Except for the soliton-enclosed region (grey shaded) for the case  22 = 0.2, the stellar discs all remain kinematically cold  2  / 2 h ≲ 1 compared to the respective host haloes.

Figure 14 .
Figure14.Disc scale heights obtained from simulations (solid) or analytical formulae (dashed, Eq. (4)) at  = 4 kpc (upper panel) and 10 kpc (lower panel) for the cases  22 = 0.2 (blue), 0.4 (orange), and 1.2 (green) over  rel = 0-3.5 Gyr.Note that  sim ≡  rel − 1.4 Gyr.Analytical expressions generally overestimate the scale height, by up to a factor of two, and are sufficiently accurate only when the disc is properly in the SGD or BGD limit, as shown in Fig.12; see Section 5.1 for details.

Figure A1 .
Figure A1.Face-on density projections of the stellar disc hosted by the halo with  22 = 0.4 and  h = 7 × 10 10 M ⊙ during the initial co-relaxation process.At  rel = 0 (upper left panel), the initial disc configuration constructed by GALIC appears smooth and free of inhomogeneous substructures.At the onset of the co-evolution (upper right), the development of spiral arm structures and local density clumps rapidly destabilises the disc.As the dischalo system evolves further in time (lower left), these transient substructures gradually disappear from the inner halo region.At the end of the co-relaxation phase  rel = 1.4 Gyr (lower right), disc clumps have largely dissolved within  ≲ 10 kpc.The entire disc remains Toomre stable  ( ≤ 15 kpc) > 1 and the disc rotation curve becomes stable in time beyond this point (see Figs. A2 and A3).

Figure A6 .R
Figure A6.The density slices of the  22 = 0.4,  h = 7.0 × 10 10 M ⊙ halo during the co-relaxation process.Relative to the initial granular structures at  rel = 0 (upper left), the typical granulation sizes shrink during the corelaxation due to the introduction of the disc component.Beyond  rel ≳ 0.4 Gyr, the physical scales of FDM granulation do not evolve significantly with time, again indicating that the FDM haloes relax more rapidly than the corresponding disc components (cf.Fig.A1).

Figure C2 .
Figure C2.The relative difference of disc rotation curve (upper) and the surface density (lower) profiles in the CDM host halo during  rel = 0-3.5 Gyr between the  rel = 3.5 Gyr snapshot.Both profiles exhibit negligible evolution after  rel ≥ 1.4 Gyr.
,  max ), (4) cause net energy transfer dominantly sourced by weak encounters  max ≡  h /2 ≫  min ≡ ( dB /2)/2, and (5) exhibit a Maxwellian velocity dispersion  h , the granulationinduced stellar heating rate H can be analytically estimated via the Fokker-Planck approximation.The validity thereof and predicted H

Table 1 .
The FDM particle mass  22 and virial halo mass  h for all six cases simulated in this work.In each setup, the production run adopts the highest disc particle number  d and spatial resolution.Simulations with lower particle resolutions (denoted with lowercase letters) are carried out to verify numerical convergence, as discussed in Appendix B.