Kinetic relaxation and nucleation of Bose stars in self-interacting wave dark matter

We revisit kinetic relaxation and soliton/Boson star nucleation in fuzzy scalar dark matter featuring short-ranged self-interactions $\mathcal{H}_{\rm int} = -\lambda|\psi|^4/2m^2$, alongside gravitational self-interactions. We map out the full curve of nucleation timescale for both repulsive ($\lambda<0$) and attractive ($\lambda>0$) short-ranged self-interaction strength, and in doing so reveal two new points. Firstly, besides the two usual terms, $\propto G^2$ and $\propto \lambda^2$, in the total relaxation rate $\Gamma_{\rm relax}$, there is an additional cross term $\propto G\lambda$ arising due to interference between gravitational and short-ranged self-interaction scattering amplitudes. This yields a critical repulsive interaction strength $\lambda_{\rm cr} \simeq - 2\pi Gm^2/v_{0}^2$, at which the relaxation rate is smallest and serves as the transition point between typical net attractive self-interaction ($\lambda \gtrsim \lambda_{\rm cr}$), and net repulsive self-interaction ($-\lambda \gtrsim -\lambda_{\rm cr}$). Secondly, while in the net attractive regime, nucleation time scale is similar to inverse relaxation time scale $\tau_{\rm nuc} \sim \Gamma^{-1}_{\rm relax}$, in the net repulsive regime nucleation occurs at a delayed time $\tau_{\rm nuc} \sim (\lambda/\lambda_{\rm cr})\Gamma^{-1}_{\rm relax}$. We confirm our analytical understanding by performing 3D field simulations with varying average mass density $\bar{\rho}$, box size $L$ and grid size $N$.

B. Peculiar appearance of high frequency modes for the repulsive case 13 Understanding the nature of dark matter (DM) is one of the main quests of modern cosmology.It could be multi-faceted in the sense that there are many degrees of freedom in the whole dark sector, for instance the String theory Axiverse [1][2][3] or other confined sector(s) (e.g.see [4,5], and also [6]).Or it could be that there is a dominant degree of freedom, such as the QCD axion [7][8][9][10][11][12], that comprises all (or most) of the dark matter.Furthermore, while the DM appears to interact only gravitationally with the Standard Model degrees of freedom (or very weakly if it does otherwise), it can still have appreciable non-gravitational self-interactions (nGSI) besides the usual gravitational self-interactions (GSI).Such is the case even for the above mentioned examples.
Of particular interest to us in this paper, is the phenomenon of kinetic relaxation and associated nucleation of Bose stars within a bath of DM waves [34][35][36][37][38][39][40][41][42].The term "kinetic" implies two key aspects: (a) The self-interactions in the field are small, allowing wave modes to freely evolve (at leading order) with the non-relativistic dispersion relation ω k = k 2 /2m.This enables a kinetic treatment of the mode occupation number function; (b) the size of the 'box' (∼ the size of a DM halo for practical purposes), is much larger than the typical fluctuation scale ℓ dB ∼ π/ k ∼ π/(mv 0 ) in the bath of DM waves.The process of kinetic relaxation is attributed to these self-interactions of the DM field which although small, over large time scales τ relax ≫ ω −1 k drive the occupation number function to develop increasing support towards smaller wavenumbers k → 0. See [34] and [41] for a relevant discussion for the cases of point-like quartic self-interactions and gravitational self-interactions respectively.Once enough particles condense into lower momentum states, their collective net attractive self-interaction becomes strong enough to counter balance their wave pressure resulting in the nucleation of a Bose star.
In this paper, we focus on investigating kinetic relaxation and subsequent Bose star nucleation for a single scalar Schrödinger field with both GSI and point-like quartic nGSI.Employing wave-kinetic Boltzmann analysis and 3D simulations, we demonstrate the presence of a previously overlooked cross term ∝ Gλ in the rate of relaxation Γ relax .(Here G denotes Newton's constant and λ represents the point-like nGSI strength).It arises due to interference between the gravitational and point-like self-interaction scattering amplitudes.The presence of this cross term gives rise to a critical nGSI (repulsive) strength λ cr ≃ −(2πG)m 2 /v 2 0 , at which the rate of relaxation reaches its minimum value (corresponding to maximum nucleation time).This critical value also serves as the transition point from typical net (contributions from both gravity and short-ranged self-interactions) attractive to repulsive self-interactions.
Because of the presence of gravitational selfinteraction, kinetic relaxation is generally accompanied with nucleation of spatially localized clumps/Bose stars, with their nucleation times dependent on the nature of the short-ranged self-interactions -attractive or repulsive.For λ ≳ λ cr , the net typical self-interaction is attractive, and nucleation happens quickly after relaxation.On the other hand for λ ≲ λ cr , the net typical self-interaction is repulsive and nucleation gets delayed.We will study relaxation and nucleation of Bose stars, and also discuss their eventual fate. 1 However we will not dwell into a careful analysis of the growth rate of these nucleated stars.See [43][44][45] for the gravity only (λ = 0) case.
The rest of the paper is organized as follows: Starting with the basic model of fuzzy scalar DM carrying both GSI and point-like nGSI in sec.II, we describe the associated wave kinetic Boltzmann equation for the evolution of the occupation number function in sec.III.Highlighting the presence of the cross term (that gives rise to λ cr ), we estimate the total rate of kinetic relaxation/condensation.In sec.IV we discuss the two cases of λ ≳ λ cr and −λ ≳ −λ cr and write down the associated nucleation time scales of spatially localized bound objects.In sec.V we discuss our 3D simulations and compare our analytical estimates with them.We also discuss eventual behavior of Bose clumps observed in simulations.Finally in sec.VI, we summarize our work and also compare our results with the existing literature on this subject.In appendix A we discuss statistical convergence of our simulations, and in appendix B we discuss a peculiarity observed in the case of repulsive short-ranged self-interactions, over longer time scales as compared to nucleation.
Conventions: Unless stated otherwise, we will work in units where ℏ = c = 1.

II. MODEL
Ignoring Hubble flow (for we are interested in sufficiently sub-horizon dynamics), the evolution of the cold/non-relativistic fuzzy scalar dark matter with both GSI and short-ranged quartic nGSI, can be described using mean field theory.The dark matter field ψ obeys the following non-linear Schrödinger (Gross-Pitaevskii) equation: Here G is the Newton's constant, and λ is the point like self-interaction strength.In our convention, λ > 0 and λ < 0 dictate attractive and repulsive selfinteraction respectively.To obtain the form Eq. ( 1), we have plugged the self-gravitational potential, Φ = 4πG ∇ −2 (mψ * ψ − ρ) ≡ 4πGm∇ −2 /0 ψ * ψ, in the usual Schrödinger-Poisson system of equations.The ∇ −2 /0 denotes exclusion of the homogeneous part of the number density field ψ * ψ.In Fourier space with the decomposition ψ(x, t) = (2π) −3 dk e −ik•x Ψ k (t), the Schrödinger equation becomes where and it is understood that k ̸ = ℓ ̸ = 0 in the above.For later convenience, it is also useful to write down the Hamiltonian density (in physical space) for the mean field ψ: Here the different terms in the above can be attributed to the wave-pressure H wp = |∇ψ|2 /2m, gravitational self-interaction H gr = mΦ|ψ| 2 , and short-ranged selfinteraction The Gross-Pitaevskii (GP) equation, being non-linear, renders it difficult to analyze and study the evolution of the ψ field in generality.However for the purposes of kinetic relaxation leading to nucleation of localized Bose clumps however, wave kinetic Boltzmann analysis can be performed which we discuss next.To test and verify our analytical understanding, we perform 3D field simulations which we discuss in a later section.

III. WAVE KINETICS AND RELAXATION
For kinetic relaxation in wave dynamics, we can study the evolution of the mode occupation number function is the volume), which is nothing but the Fourier transform of the 2-point volume averaged field correlator ζ(x, t) = V −1 dy ψ * (y, t) ψ(y + x, t).Under random phase approximation with weak interactions, the relevant wave-kinetic Boltzmann equation can be derived.See for instance [46].For a derivation for the general case of arbitrary number of fields and 2 body interactions, see [41].For the scalar case at hand, characterizing the dependence of the occupation number functions on wavenumbers as f k/m , the wave-kinetic equation takes the familiar form Here v = k/m and ṽ = p/m are the incoming "velocities" in the 2-wave interaction, and the quantities in the 1-dimensional Dirac delta function are the free wave energies E k = k 2 /2m.The quantity dσ is the effective differential cross section.The cubic nature of the terms in the right hand side bracket (∼ f x f y f z ), usually understood as Boltzmann enhancement terms, arise due to the wave-mechanical nature of the system (1) and are crucial for the phenomenon of Bose condensation.Last but not the least, it is the form of the differential cross section that appears in the wave-kinetic equation, ∼ |T | 2 , that is of utmost importance for our discussion.The scatter-ing amplitudes due to the different kinds of 2-body interactions (here gravity and point-like self-interactions), are added first and then squared : What appears in the differential cross section is |T | 2 where T = T G + T λ (c.f.Eq. ( 3)), and 2 and T λ = −λ/m 2 are real).This can be attributed to the wave dynamical nature of the GP system.The above equation (5), after integration over the Dirac deltas, can be re-written in terms of the incoming and outgoing relative velocities u = u n and u ′ = u n′ respectively 2 , by redefining p/m = k/m − u n and q/m = ℓ/m − u n′ : Let us briefly discuss the different terms in the differential cross section explicitly.Broadly speaking, there are two types of interference terms that arise.One is the interference between the t and u channels (relevant mainly for the gravitational interaction), and the second is the interference between the two different types of interactions (gravitational and short-ranged).See fig. 1 for a pictorial representation.
For GSI only (λ = 0) case, the contribution from the t and u channels are the first two terms ∝ is due to their mutual interference (as also discussed in [41]).Note that the sole contributions from the t and u channels are identical: The full integral with the | n′ + n| −4 term is identical to that with the | n′ − n| −4 term.The sole contributions give rise to the Rutherford scattering cross section, carrying a logarithmic IR divergence (aka the Coulomb logarithm), while the interference term becomes sub-dominant in the large log limit and can be omitted.
For nGSI only (G = 0) case, contributions from t and u channels are identical to their mutual interference one, and goes as λ 2 .This is simply due to the interaction being a contact/point interaction.
Importantly when both of the interactions are present, their respective scattering amplitudes (for either of the two channels) are added first and then squared.All the terms ∝ Gλ, while giving identical contributions, characterize the interference between the two types of interactions.Splitting the contributions from GSI, nGSI, and their interference, we have the following wave-kinetic Boltzmann equation Here w = u( n′ − n)/2, and Λ = log(mv 0 L) is the aforementioned Coulomb logarithm with v 0 and L equal to typical velocity and box size (or halo size for physical considerations) respectively.While the cross term and nGSI term follow straightforwardly from Eq. ( 6), the Rutherford scattering collision term C GSI is obtained after an eikonal approximation and was derived explicitly in [41].Also see [46,47] for the same equation for a scalar field.Now in order to get a typical estimate for the total relaxation rate Γ relax ≡ 1 fv ∂fv ∂t , we can replace different quantities in the three collision terms with their appropriate scalings.Replacing angular volume dΩ → 4π, typical relative velocity 0 /n, and finally the occupation number function f v → (2π)3/2 ρ/(m 4 v 3 0 ), the total relaxation rate is parameterized as in an elastic collision, i.e. |u| = |u ′ | ≡ u.
Our scaling of the occupation number is dictated by Gaussian initial condition (see Eq. ( 12) ahead) which we shall use to perform simulations, described in the next section.In general, α 1 , α 12 , and α 2 are positive O(1) coefficients that would depend on the specific initial conditions.Eq. ( 8) is our master formula for the relaxation rate.The value of λ around which the relaxation rate becomes smallest, is easily estimated to be where β = α 12 /α 2 ∼ O(1), and the associated (minimum) rate is 3 Notice that this critical value of λ cr , can also be obtained from the GP equation ( 1) by balancing the gravitational term with the self-interaction term together < l a t e x i t s h a 1 _ b a s e 6 4 = " A g e 4 7 F z H m c e T h C I 7 h F D y 4 h g r c Q x X q w A D h G V 7 h z X l 0 X p x 3 5 2 P e m n O y m U P 4 A + f z B 3 T F j L k = < / l a t e x i t > + < l a t e x i t s h a 1 _ b a s e 6 4 = " A g e 4 7 F z m c e T h C I 7 h F D y 4 h g r c Q x X q w A D h G V 7 h z X l 0 X p x 3 5 2 P e m n O y m U P 4 A + f z B 3 T F j L k = < / l a t e x i t > + < l a t e x i t s h a 1 _ b a s e 6 4 = " W a L 0 t W j l t r k u x l q V + K u q J E p v 8 n Y d s x I c P 8 g J / n z r 9 w 4 9 4 O T I y C 3 q T J u e f e k 9 5 7 r h 0 w Pictorial / Feynman graph representation of all the terms appearing in the differential cross section in Eq. ( 6).
The total contribution to the interaction rate (left hand side of the equality), is the square of the sum of both gravitational amplitude (top two graphs) and point self-interaction amplitude (bottom graph ×2).For gravitational interaction, there are two distinct channels (t and u).Their mutual interference, as compared to their sole contributions, becomes subdominant in the large log limit (leading to Rutherford scattering result).More importantly, interference between scattering amplitudes of the two different interactions matters.See main text for details.
with replacing the exchange momenta by its typical value |k − ℓ| 2 ∼ 2(mv 0 ) 2 .This criticality marks the transition point from attractive to repulsive net typical self-interactions: For λ ≳ λ cr , typical interactions within the bath of DM waves are attractive since typical T is negative, whereas for −λ ≳ −λ cr they are repulsive since typical T is positive.

IV. NUCLEATION AND BEHAVIOR OF SOLITONS A. Nucleation
In general, the process of kinetic relaxation is characterized by an increasing support of the occupation number function f k at vanishing wavenumber.(For instance see [34] and [41] for discussions of nGSI and GSI cases respectively).This implies increasing field correlation over larger length scales with diminishing density fluctuations, i.e. field homogenization.A heuristic understanding of the subsequent nucleation of a spatially localized and bound clump, can perhaps be gained most easily from a particle physics perspective, together with recalling that λ cr also marks the transition from typical net attractive self-interaction to typical net repulsive self-interaction.As particles lose kinetic energy on account of self-interactions and move towards smaller momenta (condensate state), there comes a time when within some region, the collective net potential (due to both self-gravitational and short-ranged interactions) becomes comparable to wave pressure.The time scale of this process is nothing but the inverse relaxation rate Eq. ( 8), which in the case of net attractive self-potential λ ≳ λ cr , leads to 'immediate locking' of such a region into a bound clump (having negative energy).That is, τ nuc ≃ Γ −1 relax .Strictly speaking, this can be taken as a definition of τ nuc with Γ nuc = Γ relax , in which case the different α constant coefficients in the rate Eq. ( 8), are understood as such.
On the other hand for −λ ≳ −λ cr , relaxation cannot immediately lead to nucleation of a localized bound clump.This is because the net typical interaction is repulsive: the collective self-potential within density fluctuation regions is not binding yet.Over time though, more particles get driven towards the condensate phase, and eventually there arises a potential for a bound object to nucleate (within which net gravity can now compensate for both repulsive short-ranged interaction and wave-pressure).This gives τ nuc > Γ −1 relax .In general, we can therefore write the following where h(λ) is a threshold function (or the delay factor), that relates nucleation times to relaxation rates.As mentioned earlier, relaxation means field homogenization, and we expect the rate at which the system relaxes to be comparable to the rate at which density fluctuations decrease.The delay factor can then be estimated as the ratio of typical density fluctuation at relaxation, to that at nucleation, h ∼ δρ relax /δρ nuc .While we expect this to be order unity for net attractive case (first case of Eq. 10), for large repulsive strengths it should increase with increasing −λ.Below we estimate this scaling.Consider a region of typical size ∼ (mv ) −1 where the field would have 'locked' itself into a bound configuration upon relaxation/condensation, were the net potential was binding.However this is not the case yet, and we may balance the typical self-interaction energy density (mostly due to short-ranged interactions) relax /2m 4 , with the wave pressure within H wp ∼ v 2 0 δρ relax /2.This gives δρ relax ∼ m 4 v /λ.As relaxation continues (meaning more particles are driven towards low momenta state), the value of both density fluctuations δρ and typical size of fluctuation regions (mv) −1 change.The former decreases and the latter increases so as to maintain H self ∼ H wp .Then, nucleation is expected to occur when gravity can compensate for both the wave pressure and repulsive short-ranged self-interaction.That is, we may balance (the magnitudes of) all the three energy densities, H wp ∼ v Here we have inserted another constant coefficient α 3 that depends upon the initial conditions.Through simulations, we will confirm our estimate Eq. ( 10) (together with Eq. ( 8) and Eq. ( 11)), and also extract the different α coefficients for Gaussian initial conditions.

B. Eventual behavior
Once a spatially-localized Bose clump/soliton emerges, its subsequent evolution and long term dynamics depends on whether the short-ranged self-interactions are attractive or repulsive.The full spectrum of such solitons is extensively discussed in the literature.See e.g.[20,[48][49][50].To recapitulate some of the basic points that may suffice for our purposes, consider the energy landscape for objects of radius r s and mass M s in the theory.Using Eq. ( 4), the wave pressure energy, self-gravitational potential energy, and short-ranged selfinteraction potential energy are s /(r s ), and H self = −cλM 2 s /(m4 r 3 s ) respectively, with a, b, and c some positive coefficients that depend upon the exact profile of the object.The total energy is the sum of all three.
1. Attractive short-ranged interactions λ > 0 For this case, the energy vs radius curve (for a given mass M s ) has a local minima that corresponds to quasi stable negative energy (bound) states / solitons.It is separated from the runaway behavior towards small radii, ∼ −λ/r 3 s , by a barrier whose height decreases with increasing M s .The barrier disappears at a critical mass M s,crit ∝ (λG) −1/2 , beyond which the theory does not admit any quasi-stable bound states anymore.Starting in the kinetic regime and upon relaxation, a quasi-stable Bose clump nucleates and starts to accrete mass from its surroundings.Ultimately once it accumulates enough mass such that the energy barrier gets sufficiently low, and/or it 'breathes' rapidly enough so as to be able to probe beyond the energy barrier, it 'collapses'.This is because the region now prefers to lower its energy by transitioning on the runaway ∼ r −3 s curve.This is sometimes referred to as "Bosenova" (owing to its analogy with a type II supernova).While subsequent evolution of the object beyond this criticality requires a fully relativistic analysis and has been pursued in the literature [51] (also see [52] for an associated astrophysical phenomenology), the evolution leading upto this criticality (and even beyond until the wave pressure starts to become comparable to rest mass energy) is well captured by the non-relativistic treatment.For large attractive strengths λ ≳ −λ cr , the barrier is less significant and the objects quickly collapses upon nucleation.In our simulations we indeed observe this phenomenon (see fig. 3 in sec.V ahead).

Repulsive short-ranged interactions λ < 0
In the repulsive scenario there is no runaway domain since the energy for low radii is now −λ/r 3 s > 0. This renders the previous local minima stable (hence now global minima), corresponding to bound soliton states.The critical mass M s,crit ∝ (−λG) −1/2 serves as the transition point into the Thomas-Fermi regime [20,53,54].This is where the mass of solitons gets large enough to admit comparable amounts of self-gravitational and short-ranged interaction energy densities, with gradient pressure becoming sub-dominant.As a result, the radius starts to approach a constant r s ∼ √ −λ/(m 2 √ 4πG) (with the mass dependent correction term dying out as ∼ M −1 s ).Up until the mass becomes sufficiently large where GM s,relv ∼ r s and relativistic effects start to become important (see [55,56]), the theory then admits a set of "Chandrasekhar solitons" with masses ranging anywhere between M s,crit and M s,relv , and radii approximately around r s ∼ √ −λ/(m 2 √ 4πG). 4 Though the theory admits these stable Chandrasekhar solitons, understanding their evolution and long term behavior within the bath of DM waves is crucial, and has been extensively studied in the literature.See [59][60][61][62][63] for simulation setups using the fluid/Madelung equations instead of the Schrödinger field equation for scalar wave dark matter (with repulsive short-ranged selfinteraction).For our Fourier split simulation technique (which we discuss in the next section), we find that over longer time scales (after the nucleation of Bose clumps), the system reaches some sort of criticality when high frequency modes (near cutoff) start to appear in the simulation box.This leads to a breakdown of the simulation (along with the disruption of the clump), visible in the form of a checkerboard-like pattern.We present this peculiar artifact from our simulations in appendix B, although a detailed investigation of it is left for future work.

V. FIELD SIMULATIONS
To verify our analytical understanding of kinetic relaxation and associated nucleation of bound Bose stars, we have carried out a large suite (∼ 500) of 3D simulations of the GP system (1) with varying values of the nGSI strength λ.We evolve the GP system (1) with the following initial Gaussian function for the k-space Schrödinger field where θ k/m are random phases, uniformly distributed in (0, 2π), for every wavenumber k.Our numerical algorithm is based on the well known split Fourier technique/pseudo-spectral method [15,[64][65][66][67][68], and we have used both Python based and Matlab based codes to generate our simulation data.
As mentioned earlier, in order to be in the kinetic regime we require (i) interactions to be tiny as compared with the typical free wave evolution (occurring over time scales ∼ 2/mv 2 0 ), and (ii) the box size to be larger than the typical field fluctuation scale π(mv 0 ) −1 .Furthermore, we also impose the box size to be smaller than the gravitational Jeans scale associated with a incoherent bound halo ℓ J ∼ v 0 (π/Gρ) 1/2 , in order to avoid its formation within our simulation box.In this sense, our simulation box of a collection of DM waves with typical fluctuation scale ∼ π(mv 0 ) −1 , may be regarded as a region within a DM halo.In summary we require the following to hold true Kinetic regime : In our simulations, we work with dimensionless quantities, for which purpose we set G = 1/(8π) and m = 1.More explicitly, one can rescale different quantities in the fashion t → t/E, x → x/ √ m E, ψ → ψ E/ √ 8πGm and λ → λE/(8πGm 3 ) to get Eq.( 1) with both 8πG and m replaced by unity.Here E is any reference energy scale in the system (for instance E = mv 2 0 /2).The discretization in space is simply ∆x = L/(N x − 1) where L and N 3 x are the box size and number of grid points respectively, and the time discretization is ∆t = 2π(∆x) 2 m/(3η) with η ≥ 1. 5 In the split Fourier technique, the field evolution is split into a drift part where it is evolved solely due to the gradient term (free field evolution), and a kick part where it is evolved solely due to interactions.The Courant-Friedrichs-Lewy (CFL) condition ensures that the fastest process in the dynamics is captured appropriately.Hence, the fastest amongst the kick and drift processes, at any time iteration, sets the time discretization ∆t (e.g.see [67,68] for details).In the kinetic regime, the time discretization is always set by the free field evolution term ∼ (∆x) −2 /2m, and hence by space discretization as given above.
In all of our simulations we set v 0 = 1/ √ 2, and choose the box size and average mass density such that we are deep in the kinetic regime.Most of the simulations were performed with L = 40, L = 45, and L = 50 box sizes, and the average mass densities were chosen to be small enough such that the factor Γ −1 relax mv 2 0 /2 was at-least as large as ∼ 250, going all the way up-to even ∼ 4500.For robustness, we have performed simulations with different grid sizes N x = {192, 216, 256}, scanning over different λ/λ cr values.(We also performed simulations with N x = 150 and N x = 300 to test convergence of our results (see Appendix A).
To capture the formation of localized Bose clumps, we keep track of the mass density in the box, radially averaged (in k space) occupation number function f k , the associated volume averaged correlation function ζ(r), and the maximum mass density in the box ρ max . 6ucleation of a localized clump can be characterized by a change of trend of ρ max , wherein it starts to monotonically increase beyond just the statistical fluctuations that happen over short time scales.We record the corresponding times in all of our simulations, both by direct inspection and statistical methods such as moving average. 7n order to gauge the validity of our analytical estimate of nucleation times (c.f.Eq (10) with Eq. ( 8)), and to extract the different α coefficients, we split the data set into two, with λ = −2πGm 2 /v 2 0 being the splitting point.Below we elaborate on the statistical analysis we performed in the two regimes.

Net attractive interactions (λ ≳ λcr).
In order to test the λ dependence of our estimate, we construct the quantity r(λ) = log(mv 0 L)(τ nuc (0) − τ nuc (λ))/2τ nuc (λ) using Eq. ( 8) and Eq. ( 10).This gets rid of the ρ and L dependence, giving ).Solid gray curves are from the theory estimate Eq. ( 16) (c.f.Eq. (10) with Eq. ( 8) and Eq. ( 11)), where the different α coefficients are obtained from least square fitting as described in the main text.The different colored 1σ bars are from simulations (performed with Gaussian initial conditions (12)), with box sizes L = 36 (brown), L = 40 (pink), L = 45 (magenta), and L = 50 (blue), and varying average densities ρ.Here we have only plotted one theory curve for L = 50 (all curves for the four different box sizes lie practically on top of each other since the L dependence is quite mild).To show the effect of the interference term in the relaxation rate and the delay factor in the nucleation time, we have also plotted dotted and dashed gray curves.Respectively, these correspond to when the interference term from the relaxation rate is set to zero, and the delay factor in the nucleation time scale is set to unity.This delay factor is only relevant in the λ ≲ 2πGm 2 /v 2 0 case, and the dashed gray curve in the left panel is simply the extension of the main solid gray curve in the right panel.Bottom panel : Normalized ρmax (by their respective initial values) vs time curves for six different λ values, highlighted by colored points in the upper panel.The points of 'sudden' rise correspond to nucleation of respective localized Bose clumps.
Not only is the curve simple enough to do statistics with, this way we can also combine all of our simulation data (with different ρ and L).The analogous quantity for simulations is where hats denote simulation data and angle brackets denote averaging over all of the data (for a given λ value).
To extract the ratios α 12 /α 1 and α 2 /α 1 for the theory curve, we construct the cost function for least square fitting.Here N λ is the number of different simulations performed for a given λ value.Minimizing this cost function then fetches the optimal values for α 12 /α 1 , and α 2 /α 1 .For α 1 , we simply find the average of τnuc (0)/τ nuc (0), which we then use to get α 12 and α 2 from the previous two ratios.For our Gaussian initial condition (12), we found α 1 ≃ 0.8, α 12 ≃ 1.2, and α 2 ≃ 1.2.

Net repulsive interactions (−λ ≳ −λcr)
In this case, we expect nucleation to happen later than relaxation, given by Eq. ( 10) with the delay factor h(λ) in Eq. (11).We can use the previous case relationship τ nuc (λ ≳ λ cr ) ≃ Γ −1 relax , to test the scaling of h(λ) for the −λ ≳ −λ cr case.From simulations, we construct τnuc Γ relax with the three αs in the relaxation rate set to the ones obtained above.We then perform least square fitting by constructing the cost function similar to the previous case and minimizing it.For Gaussian initial conditions, we found α 3 ≃ 1.
With the above analysis and all the four α values obtained, upper panel of fig. 2 shows our main plot.We plot τ nuc (λ) (normalized by τ nuc (0)) as a function of λ(v 2 0 /2πGm 2 ): The error bars correspond to 1-σ fluctuations (owing to random and different initial condition for every simulation seed), with different colors corresponding to three different box sizes considered.In general, the agreement between analytical estimates and simulations is evident.Let us highlight our two main results: (a) The rising feature as λ goes from positive to negative, with a peak occurring around λ ≃ −2πGm 2 /v 2 0 ≃ λ cr , is a clear evidence of the interference term in the relaxation rate.To represent the effect of the interference term visually, we have also plotted a dotted gray curve (in the upper right panel of fig.2), which is equal to inverse of the relaxation rate with the interference term dropped.That is, inverse of Eq. ( 8) with the term ∝ Gλ set to zero; (b) To the left of the peak and increasing −λ, nucleation happens later than just the inverse relaxation rate.The delay factor h and the relaxation rate Γ relax scale as ∼ −λ and λ 2 (to leading order) respectively, resulting in the scaling of the nucleation time as (−λ) −1 (and not λ −2 ) to leading order.To highlight this, we have augmented the upper left panel of fig. 2 with just the relaxation time curve, i.e.Γ −1 rel , shown in dashed gray.(This is nothing but the solid gray curve on the right upper panel, extended towards the left upper panel).
The lower panel of fig. 2 shows moving averaged ρ max vs time curves for six simulations with different λ values.The unambiguous "sudden" rise in ρ max marks the nucleation of a localized object within which density grows over time. 8s visual examples, in fig. 3 we also present density projection snapshots for six different λ values at later times, showing the presence of nucleated Bose stars.For attractive short-ranged self-interaction λ > 0, nucleated Bose clumps eventually collapse into a Bosenova.This happens when it reaches the critical mass where it can no longer remain stable (see section IV B).

VI. SUMMARY AND DISCUSSION
In this paper we have investigated kinetic relaxation and associated nucleation times of Bose stars, in scalar fuzzy dark matter with short-ranged 2-body self-interactions.Starting with the wave-kinetic Boltzmann equation for the mode occupation number function (which we derived in an earlier work), we first highlighted the presence of a cross/interference term ∝ Gλ in the rate of relaxation Γ relax , alongside the usual two terms ∝ G 2 and λ 2 due to both gravitational and shortranged self-interaction individually.This is because of the wave-mechanical nature of the system: The rate depends on the total cross section, which is not just the sum of individual cross sections due to the different processes.Rather the scattering amplitudes due to all the processes must be added first, and then use its absolute square to get the cross section and associated rate of relaxation/condensation.
The presence of this cross term gives rise to a critical repulsive self-interaction strength λ cr ≈ −2πGm 2 /v 2 0 , around which the typical net self-interaction (due to both gravitational and short-ranged self-interaction), transitions between being attractive and repulsive, and the relaxation rate becomes smallest.Here k 0 = mv 0 is the typical wave-mode present in the system initially.
For nucleation times as a function of λ, we found that for the case of net attractive self-interaction λ ≳ λ cr , nucleation happens quickly upon relaxation, giving rise to the relationship τ nuc ≃ Γ −1 relax .One the other hand for net repulsive self-interaction −λ ≳ −λ cr , nucleation is delayed.This is because upon relaxation, short-ranged self-interaction dominate over gravitational self-interaction, preventing nucleation of a bound object.Over time as more particles are driven towards the condensate phase (equivalently, as field correlation length scale increases along with diminishing density fluctuations), a potential arises for the formation of a bound region where gravity can now overcome both       the wave pressure and short-ranged self-interaction.
The associated delay factor rises linearly with −λ, giving the nucleation time scale as τ nuc ≃ (λ/λ cr )Γ −1 relax .In summary, Eq. ( 10) along with Eq. ( 8) (with the delay factor give in Eq. ( 11)) is our main analytical estimate for the nucleation timescale of Bose stars, as a function of the short-ranged self-interaction strength λ.
To analyze this, we performed a large suit of 3+1 dimensional simulations of the Schrödinger -Poisson / Gross-Pitaevskii system (Eq.( 1)), for many different values of λ and different parameters such as the box size L and average mass density ρ.All of our simulations were carried out with Maxwell-Boltzmann distribution, with random phases for each value of the wavemode k (Eq.( 12)).Throughout most of our simulations, we kept track of the max density in the box, occupation number function f k (radially averaged f k in k space), the associated correlation function ζ(r), and projected mass density along some direction.By reading the times at which ρ max starts to monotonically rise beyond just the statistical fluctuations (together with making sure that a localized over dense region does appear in the simulation box around this time), we record the times of nucleation.Upper panel of Fig. 2 presents the comparison between simulations and analytical estimate.As examples, the figure is also appended (lower panel) with ρ max vs time curves for six different λ values.
While in this paper we have not analyzed our simulation data for the rate at which Bose stars accrete mass from their surroundings, we kept track of the eventual behavior of these objects (post nucleation), for many of our simulations.For the attractive case (λ > 0), we confirmed that the nucleated Bose stars eventually decay away.This is expected since there exists a maximum critical mass beyond which the star becomes unstable and collapses into a Bosenova.For instance see upper right snapshot in Fig. 3, when the nucleated star 'immediately' collapses.While the study of eventual dynamics and fate of such regions requires a full relativistic treatment, field dynamics up-to this point is well described by the nonrelativistic GP equation (e.g.see [51]).
For the repulsive case λ < 0, we found a peculiar decay behavior.We find that the nucleated clump eventually (over time scales longer than the nucleation time) reaches a type of criticality at which point very high frequency modes, passing through the clump and travelling along the three directions of the simulation box, appear in the system.See appendix B for some discussion.This could be an artifact of the periodic boundary conditions of the split-Fourier simulation setup, and if so, bringing into question its use to study long term dynamics of fuzzy dark matter with repulsive short-ranged self-interactions via such simulation setups.We leave a detailed investigation of this behavior for a separate work.

A. Comparison with earlier work
Let us now compare our results with some of the earlier work on the subject of kinetic nucleation of Bose stars.First, our results encompass the result of [35] for the gravity only (λ = 0) case, and is even in very good agreement with the order unity coefficient α 1 in the rate expression (besides the overall scaling with ρ, m, v 0 , L and G), obtained for Gaussian initial condition.Upon inclusion of short-ranged self-interaction (λ ̸ = 0 case), our results differ significantly from the existing literature [37][38][39].First, we find that there exists an interference term ∝ Gλ in the relaxation rate, which in fact is the leading order λ dependent term when short-ranged self-interaction is not dominating over gravitational self-interaction.Only in the scenario when the former is dominant, does the relaxation rate goes as λ 2 to leading order.Secondly, the nucleation time scale is not always equal to the inverse relaxation rate.While for λ ≳ λ cr , nucleation time scale is just the inverse relaxation rate, for the strong repulsive self-interaction −λ ≳ −λ cr , nucleation time is delayed by an extra factor of (λ/λ cr ).Therefore for the purposes of nucleation of Bose stars, only in the case of strong attractive short-ranged self-interaction, λ ≳ −λ cr , is it true that the nucleation time goes as λ −2 to leading order.For in the opposite case of strong repulsive shortranged self-interaction, −λ ≳ −λ cr , the nucleation time scale goes as λ −1 to leading order instead.

B. Implications
Our results could have important implications in the context of self-interacting fuzzy dark matter and various interesting phenomenon that it entails.The appearance of the interference term in the relaxation rate, and hence in the nucleation time scale of Boson stars, may modify results for some of the phenomenon such as recurrent axinovae [52], de-stabilization of gravitational atoms [69], among others.
In general, irrespective of the nature (attractive or repulsive) of point-like self-interaction, the interference term becomes the leading order λ dependent term (and hence extremely important), when |λ| is at best comparable to the critical value |λ cr |.As an example, even for the QCD axion we have λ qcd /|λ cr | ≃ 1.3(v 2 0 m 2 pl /f 2 a ), hence becoming comparable to, or less than |λ cr |, in cosmological environments with v 0 ≲ (f a /m pl ) ∼ 10 −5 .For instance this could be important in the study of axion miniclusters [70].
In this paper we have focused on kinetic nucleation via both gravitational and short-ranged self-interactions for a single scalar field.A natural generalization is to include multiple scalar fields with naturally different masses and 4-point interactions, or a single spin-1 field including density-density and spin-spin interactions [57,71], or even multiple spin-1 fields with extra Yang-Mills interactions [57], or a combination thereof.While there would necessarily be interference terms ∝ Gλ, and we expect similar scaling of nucleation time scales as presented in this work (as a function of λ), a detailed analysis of such cases is left for future work.The solid gray curves correspond to the analytical estimate Eq. ( 16) (c.f.Eq. ( 10) with Eq. ( 8) and Eq. ( 11)), with the O(1) α coefficients obtained using simulations performed with grid sizes 192 3 , 216 3 , and 256 3 as in the main text.
these waves circulate within the periodic simulation box.An artifact of this is wiping out of density fluctuations and eventual homogenization of the simulation box.
In fig. 5 we provide some simulation snapshots of this peculiarity, for three different values of λ.Notice the appearance of checkerboard like pattern in the right hand side panel.
In order to check if the phenomenon is an artifact of finite discretization, we simulated lower and higher resolution grids with the same initial conditions for long times.We found no clues if this is the case or not.As an example, for λ ≈ −4.2λ cr , with L = 40 and ρ = 0.01, in the lower resolution (128 3 ) and higher resolution grids (256 3 ) the nucleation times were ∼ 3600 and ∼ 3800 respectively.This shows decent convergence of the nucleation time.At the same time however, the onset of < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 p R L 0 u 3 K   these high frequency waves for both the grids were also similar (∼ 5300 and ∼ 5700 respectively).
We also performed two other tests to see if something can be learned about this phenomenon.In one test we put a pre-computed soliton in a bath of DM waves.We observed the same appearance of high frequency modes and checkerboard pattern developing, leading to disruption of the soliton.The times at which this happens depended upon the average mass density and total mass of the soliton.In yet another test, we performed a few runs with larger average mass densities, such that the associated gravitational Jeans scale ℓ J ∼ v (π/Gρ) 1/2 became smaller than the box size.Even in this case, upon formation of a halo, the same peculiar phenomenon appears.
While we have confirmed the presence of this feature in our (split-Fourier technique based-)simulations through multiple tests, we leave a detailed analysis of this peculiarity for future work.

T = m 2 < l a t e x i t s h a 1 _
r W e 8 W 7 3 A x x 5 x J e Y I S E 6 l h n K b o y 4 p J g R 1 T s S J E R 4 h A a k o 6 C P P C K 6 8 c y 0 K c w r p g / d gK v n S z h j F x U x 8 o Q Y e 4 6 q T F Y Q y 7 m E / C v X i a R 7 2 4 2 p H 0 a S + H j e y I 0 Y l A F M L g D 7 l B M s 2 V g B h D l V s 0 I 8 R M o O q e 6 U m G A t r 7 w K m q W i d V 2 8 e i r n K v e p H R l w C i 7 A J b D A D a i A R 1 A D D Y D Bp 3 a i Q e 1 c + 9 L P 9 J y e n 5 f q W q o 5 B r 9 C L 3 4 D 9 A u 7 Z w = = < / l a t e x i t > b a s e 6 4 = " K T 7 2 a 5 j w S h g k n d K B h s 7 D t m N C B v 7 C L 3 P n r 7 j R D k y i g t 6 k y e m 5 9 / T e e + q E j A p p m u + a v r a + s bm V 2 c 7 u 7 O 7 t H x i H R y 0 R R B y T J g 5 Y w D s O E o R R n z Q l l Y x 0 Q k 6 Q 5 z D Sd s b 3 S b 7 9 Q r i g g d + Q k 5 D 0 P D T 0 q U s x k o r q G 6 9 5 2 0 N y h B G L G 7 N + D d 7 C g u 1 y h O O y H V J Y g 9 5 T a R Z P x 7 A A b c L Y V N 2 y K 4 p / B M + r 1 T Z T k w 3 Q d 5 e U m M W J L G t L 6 h E B S 3 0 j Z x b N e c B V Y K U g B 9 K o 9 4 0 3 e x D g y

5 FIG. 2 .
FIG. 2.Upper Panel : Our main figure showing the nucleation time τnuc (normalized by the gravity only case) as a function of short-ranged self-interaction strength λ (normalized by the critical factor 2πGm 2 /v 2 0 ).Solid gray curves are from the theory estimate Eq. (16) (c.f.Eq.(10) with Eq. (8) and Eq.(11)), where the different α coefficients are obtained from least square fitting as described in the main text.The different colored 1σ bars are from simulations (performed with Gaussian initial conditions (12)), with box sizes L = 36 (brown), L = 40 (pink), L = 45 (magenta), and L = 50 (blue), and varying average densities ρ.Here we have only plotted one theory curve for L = 50 (all curves for the four different box sizes lie practically on top of each other since the L dependence is quite mild).To show the effect of the interference term in the relaxation rate and the delay factor in the nucleation time, we have also plotted dotted and dashed gray curves.Respectively, these correspond to when the interference term from the relaxation rate is set to zero, and the delay factor in the nucleation time scale is set to unity.This delay factor is only relevant in the λ ≲ 2πGm 2 /v 2 0 case, and the dashed gray curve in the left panel is simply the extension of the main solid gray curve in the right panel.Bottom panel : Normalized ρmax (by their respective initial values) vs time curves for six different λ values, highlighted by colored points in the upper panel.The points of 'sudden' rise correspond to nucleation of respective localized Bose clumps.
along with our simulation data.Here Λ = log(mv 0 L) is the Coulomb logarithm, and h(x) is unity for x ≳ −1 (right upper panel of fig 2) while linearly increasing as −x ≳ −1 (left upper panel of fig 2).Note that we have only plotted one curve for L = 50 (solid gray), since the dependence on L is very mild and renders different curves for different values of L practically on top of each other.
< l a t e x i t s h a 1 _ b a s e 6 4 = " O 48 x F 7 G O p Z v B s J z T R t w i U U a + w v Q = " > A A A C D H i c b V D L S g M x F M 3 4 r P V V d e k m W A R X d W Z 8 b o S i C 1 1 W s A / o T E s m k 2 l D k 8 y Q Z A q l 9 A P c + C t u X C j i 1 g 9 w 5 9 + Y t r P Q 1 g O B w z n n c n N P k D C q t G 1 / W w u L S 8 s r q 7 m 1 / P r G 5 t Z 2 Y W e 3 p u J U Y l L F M Y t l I 0 C K M C p I V V P N S C O R B P G A k X r Q u x n 7 9 T 6 R i s b i Q Q 8 S 4 n P U E T S i G G k j t Q t Fj 5 l w i O A V P I G e p p w o 6 H o J h b e Q t 9 z j f t t u u S Z l l + w J 4 D x x M l I E G S r t w p c X x j j l R G j M k F J N x 0 6 0 P 0 R S U 8 z I K O + l T 1 5 n t T c k n N e O r s / L Z a v s z p y Y B 8 c g C P g g A t Q B n e g A q o A g 0 f w D F 7 B m / V k v V j v 1 s c 0 u m B l M 3 v g D 6 z P H 8 t j m O 0 = < / l a t e x i t > = 3 ⇥ 2⇡Gm 2 /v 2 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " A p 0 c s D / i k j Z q k T S s N 5 T 1 J F T e y L Q = " > A A A B 8 n i c b V D L S g M x F M 3 U V 6 2 v q k s 3 w S K 4 K j P i a y M U 3 b i s Y B 8 w H U o m k 2 l D M 8 m Q 3 B H K 0 M 9 w 4 0 I R t 3 6 N O / / G t J 2 F t h 4 I H M 4 5 l 9 x 7 w l R w A 6 7 7 7 Z R W V t f W N 8 q b l a 3 t n d 2 9 6 v 5 B 2 6 h M U 9 a i S i j d D Y l h g k v W A g 6 C d V P N S B I K 1 g l H d 1 O / 8 8 S 0 4 U o + w j h l Q U I G k s e c E r C S 3 x M 2 G h F 8 g 9 1 + t e b W 3 R n w M v E K U k M F m v 3 q V y 9 S N E u Y B C q I M b 7 n p h D k R A O n g k 0 q v c y w l N A R G T D f U k k S Z o J 8 t v I E n 1 g l w r H S 9 k n A M / X 3 R E 4 S Y 8 Z J a J M J g a F Z 9 K b i f 5 6 f Q X w d 5 F y m G T B J 5 x / F m c C g 8 P R + H H H N K I i x J Y R q b n f F d E g 0 o W B b q t g S v M W T l 0 n 7 r O 5 d 1 i 8 e z m u N 2 6 K O M j p C x + g U e e g K N d A 9 a q I W o k i h Z / S K 3 h x w X p x 3 5 2 M e L T n F z C H 6 A + f z B + Y C k F w = < / l a t e x i t > = 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " b f F c G w 4 k W + d Z 2 E c h r o o 0 L L o J p c I = " > A A A C E n i c b V B N S 8 N A E N 3 4 W e t X 1 a O X x S I o S E l K / T g W v X i s Y F V o S p l s t 3 Z x s w m 7 E 7 G E / A Y v / h U v H h T x 6 s m b / 8 Z t m o N f D w Y e 7 8 0 w M y + I p T D o u p / O 1 P T M 7 N x 8 a a G 8 u L S 8 s l p Z W 7 8 w U a I Z b 7 N g 7 N b 2 8 s / u d 1 E h w c d V O h 4 g S 5 Y p N F g 0 R S j O g 4 H 9 o X m j O U I 0 u A a W F v p W w I G h j a F M s 2 B O / 3 y 3 / J R b 3 m H d T 2 z x r V 5 n E R R 4 l s k i 2 y Q z x y S J r k l L R I m z B y T x 7 J M 3 l x H p w n 5 9 V 5 m 7 R O O c X M B v k B 5 / 0 L M y 2 d 1 Q = = < / l a t e x i t > t ⇡ 2.4 ⌧nuc < l a t e x i t s h a 1 _ b a s e 6 4 = " 0 w u 1 L W B w 0 M b Y B 5 G 4 I 3 e / I 8 q Z 2 U v P P S 2 e 1 p s X w 1 j S N H 9 s k B O S Y e u S B l c k M q p E o Y e S T P 5 J W 8 O U / O i / P u f E x a F 5 z p z B 7 5 A + f z B 5 K i n I c = < / l a t e x i t > t ⇡ ⌧nuc < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 I I 6 s o s e 5 e L C G / e P J l s l + k + f N 0 + 9 2 z w e 6 r V R 0 b 7 A F 7 y B 6 x n L 1 g u + w t 2 2 N D J t g n 9 o 2 d s R / J 5 + R 7 8 j P 5 d R 5 d S 1 Y z 9 9 k / S H 7 / A e O d q z Y = < / l a t e x i t > t ⇡ 1.2 ⌧nuc < l a t e x i t s h a 1 _ b a s e 6 4 = " N + M 5 y k 8 8 m C 7 0 6 u l 2 5 H l S t w 3 r 3 D i 7 P y 1 X r v M 6 i m A P 7 I M j Y I E L U A F 3 o A p q A I N H 8 A x e w Z v 2 p L 1 o 7 9 r H N F r Q 8 p l d 8 A f a 5 w 8 l 1 Z m W < / l a t e x i t > = 1.2 ⇥ 2⇡Gm 2 /v 2 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " U 1 Z + W Q s O G M + t I i Z + q 5 s b 8 s + S 9 S Q = " > A A A C D 3 i c b V C 7 T s M w F H X K q 5 R X g J H F o g K x E J I I C g t S B Q O M R a I P q U k j x 3 F b q 8 5 D t l O p i v o H L P w K C w M I s b K y 8 T e 4 b Q Y o H M n S 0 T n n 6 v o e P 2 F U S N P 8 0 g o L i 0 v L K 8 X V 0 t r 6 x u a W v r 3 T E H H K M a n j m M W 8 5 S N B G I 1 I X V L J S C v h B I U + I 0 1 / c D 3 x m 0 P C B Y 2 j e z l K i B u i X k S 7 F C O p J E 8 / d J g K B w h e w m P T q E B H 0 p A I a D s J h T c w 7 N g n Q 8 / s 2 J 5 e N g 1 z C v i X W D k p g x w 1 T / 9 0 g h i n I Y k k Z k i I t m U m 0 s 0 Q l x Q z M i 4 5 q S A J w g P U I 2 1 F I 6 S 2 u t n 0 n j E 8 U E o A u z F X L 5 J w q v 6 c y F A o x C j 0 V T J E s i / m v Y n 4 n 9 d O Z f f C z W i U p J J E e L a o m z I o Y z g p B w a U E y z Z S B G E O V V / h b i P O M J S V V h S J V j z J / 8 l D d u w K s b Z 3 W m 5 e p X X U Q R 7 Y B 8 c A Q u c g y q 4 B T V Q B x g 8 g C f w A l 6 1 R + 1 Z e 9 P e Z 9 G C l s / s g l / Q P r 4 B K q e Z m Q = = < / l a t e x i t > = 0.6 ⇥ 2⇡Gm 2 /v 2 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " R b B x c b h n x f + P Z f b 4 z G I b P l p t 3 R 8 = " > A A A C E n i c b V B N S 8 N A E N 3 4 W e t X 1 a O X x S I o S E m K H z 0 W v X h U s C o 0 I U y 2 2 3 Z x s w m 7 E 7 G E / A Y v / h U v H h T x 6 s m b / 8 Z t 7 U G r D w Y e 7 8 0 w M y 9 K p T D o u p / O 1 P T M 7 N x 8 a a G 8 u L S 8 s l p Z W 7 8 0 S a Y Z b 7 F E J v o 6 A s O l U L y F FIG. 3. Density projection snapshots for 6 different λ values, at different times t in the respective simulations.Upper Panel : Snapshots for three λ values in the typical net attractive regime λ ≳ λcr ≃ −2πGm 2 /v 2 0 .In the rightmost snapshot, for λ ≈ 3 λcr, the nucleated Bose clump quickly collapses (within 5dt), shown in the smaller right corner image.Bottom panel : Snapshots for three different λ values for the other case of typical net repulsive self-interactions −λ ≳ −λcr ≃ 2πGm 2 /v 2 0 .

2 FIG. 4 .
FIG.4.Similar to the upper panel of figure2, for simulations performed with lower and higher resolution grids compared to the ones used in the main text.With box size L = 40, the 1σ bars represent simulations performed using 150 3 grid, while the solid points are from simulations using 300 3 grid.The solid gray curves correspond to the analytical estimate Eq. (16) (c.f.Eq. (10) with Eq. (8) and Eq.(11)), with the O(1) α coefficients obtained using simulations performed with grid sizes 192 3 , 216 3 , and 256 3 as in the main text.
S H 7 J J 9 4 p N T U i G X 5 I p U C S O P 5 J m 8 k j f n y X l x 3 p 2 P U e u M M 5 7 Z J n / g f P 4 A M Z G d 1 A = = < / l a t e x i t > t ⇡ 1.4 ⌧nuc < l a t e x i t s h a 1 _ b a s e 6 4 = " G X Z 8 o S D i 1 n R e W V I K G I D e g w f 6 J X 0 = " > A A A C D n i c b V D L S s N A F J 3 4 r P U V d e l m s B T c G J N Q H x u h 6 E K X F e w D m j R M J p N 2 6 O T B z K R Q S r / A j b / i x o U i b l 2 7 8 2 + c t l l o 6 4 G B w z n n c u c e P 2 V U S N P 8 1 p a W V 1 b X 1 g s b x c 2 t 7 Z 1 d f W + / I Z K M Y 1 L H C U t 4 y 0 e C M B q T u q S S k V b K C Y p 8 R p p + / 2 b i N w e E C 5 r E D 3 K Y E j d C 3 Z i G F C O p J E 8 v O 0 y F A w S v 4 E n F s B 1 J I y K g 7 a Q U 3 s K o Y 5 8 O P L N j e 3 r J N M w p 4 C K x c l I C O W q e / u U E C c 4 i E k v M k B B t y 0 y l O 0 J c U s z I u O h k g q 7 D O j b P 7 S q l 6 n d d R A I f g C B w D C 1 y A K r g D N V A H G D y C Z / A K 3 r Q n 7 U V 7 1 z 5 m 0 S U t n z k A f 6 B 9 / g D N + p l v < / l a t e x i t > = 4.2 ⇥ 2⇡Gm 2 /v < l a t e x i t s h a 1 _ b a s e 6 4 = " I p s g H J V P M 4 I d b A D R + 7 v e T 7 4 f z FIG. 5. Density projection snapshots for 3 values or repulsive short-ranged self-interaction strength λ, at two different times in the respective simulations.The grid size is 216 3 , and box length is 50.The left hand side snapshots are the same as presented in the main text in fig. 3. The right hand side snapshots are at later times, when the simulation has now been rendered 'unfaithful'.Notice the checkerboard like pattern in all of the snapshots in the right panel.