$\sigma(E,E^{\prime})\subseteq\tau_{\Sigma^{\prime}}$, $S$ is a
$\tau_{\Sigma^{\prime}}$-closed subset of $K$ and hence $S$ is
$\tau_{\Sigma^{\prime}}$-compact. Therefore, by Lemma 4.4, $S$ is a
$\tau_{\Sigma^{\prime}}$-null sequence. Thus the space
$(E,\tau_{\Sigma^{\prime}})$ has the $p$-Schur property.

(ii) Taking into account (i) we have to show only that if $E$ has the strict
$DP_{p}$ property, then $E$ has the $DP$ property. Let $K$ be an absolutely
convex weakly compact subset of $E$. Fix an operator $T$ from $E$ into a
Banach space $L$ which transforms bounded sets into relatively weakly compact
sets. By (v) of Theorem 4.2, it suffices to show that $T(K)$ is a compact
subset of $L$. Since $L$ is a Banach space, it suffices to prove that every
sequence $\\{T(x_{n})\\}_{n\in\omega}$ in $T(K)$ has a subsequence converging
to a point of $T(K)$. As $E$ is weakly sequentially $p$-angelic, $K$ is a
weakly sequentially $p$-compact subset of $E$, and hence the sequence
$\\{x_{n}\\}_{n\in\omega}$ has a subsequence $\\{x_{n_{k}}\\}_{k\in\omega}$
which weakly $p$-converges to some point $x\in K$. Since $E$ has the strict
$DP_{p}$ property, $T$ is $p$-convergent and hence $T(x_{n_{k}})\to T(x)\in
T(K)$. Thus $T(K)$ is a compact subset of $L$, and hence $E$ has the $DP$
property.

(iii) By Theorem 4.9, to show that the space $E$ has the sequential
$DP_{(p,q)}$ property it suffices to prove that (1) the topology
$\tau_{S_{q}}$ of uniform convergence on weakly $q$-summable sequences of
$E^{\prime}_{\beta}$is compatible with the original topology $\tau$, and (2)
the space $(E,\tau_{S_{q}})$ has the $p$-Schur property. By (v) of Lemma 4.8,
we have
$\sigma(E,E^{\prime})\subseteq\tau_{S_{q}}\subseteq\tau_{\Sigma^{\prime}}$.
Now (i) of Lemma 4.1 implies that $\tau_{S_{q}}$ is compatible with $\tau$ and
hence (1) is satisfied. As we showed in (i), the space
$(E,\tau_{\Sigma^{\prime}})$ has the $p$-Schur property. Since
$\sigma(E,E^{\prime})\subseteq\tau_{S_{q}}\subseteq\tau_{\Sigma^{\prime}}$, it
follows that $(E,\tau_{S_{q}})$ has the $p$-Schur property, so (2) holds true.

(iv) Taking into account (ii) it suffices to show that $E$ has the strict
$DP_{p}$ property. To this end, by Theorem 4.5, it suffices to prove that the
space $(E,\tau_{\Sigma^{\prime}})$ has the $p$-Schur property.

Suppose for a contradiction that $(E,\tau_{\Sigma^{\prime}})$ does not have
the $p$-Schur property. Since $\tau_{\Sigma^{\prime}}$ is compatible with
$\tau$, there exists a weakly $p$-summable sequence
$S=\\{x_{n}\\}_{n\in\omega}$ in $E$ which does not converge to $0$ in
$\tau_{\Sigma^{\prime}}$. Therefore there are $\varepsilon>0$, a strictly
increasing sequence $\\{n_{k}\\}_{k\in\omega}$ in $\omega$, and a sequence
$\\{\chi_{k}\\}_{k\in\omega}$ in some $K\in\Sigma^{\prime}(E^{\prime})$ such
that

$\big{|}\langle\chi_{k},x_{n_{k}}\rangle\big{|}\geq\varepsilon\quad\mbox{ for
every $k\in\omega$}.$ (5.5)

By the definition of $\Sigma^{\prime}(E^{\prime})$, the set $K$ is weakly
compact in $E^{\prime}_{\beta}$. Therefore, by the weak sequential
$q$-angelicity of $E^{\prime}_{\beta}$, there is a subsequence
$\\{\chi_{k_{j}}\\}_{j\in\omega}$ of $\\{\chi_{k}\\}_{k\in\omega}$ which
weakly $q$-converges to some $\chi\in K$. Since $S$ is weakly $p$-summable so
is its subsequence $\\{x_{{n_{k_{j}}}}\\}_{j\in\omega}$. Then (5.5) and the
sequential $DP_{(p,q)}$ property of $E$ imply

$\varepsilon\leq\big{|}\langle\chi_{k_{j}},x_{{n_{k_{j}}}}\rangle\big{|}\leq\big{|}\langle\chi_{k_{j}}-\chi,x_{{n_{k_{j}}}}\rangle\big{|}+\big{|}\langle\chi,x_{{n_{k_{j}}}}\rangle\big{|}\to
0\quad\mbox{ as $j\to\infty$}.$

This contradiction shows that $(E,\tau_{\Sigma^{\prime}})$ has the $p$-Schur
property.

Proposition 3.3(i) of [2] states that if $E$ is a barrelled quasi-complete lcs
with the $DP$ property, then $E$ has the sequential $DP$ property. The next
corollary essentially generalizes this result.

###### Corollary 5.13

Let $E$ be a quasibarrelled locally complete space. If $E$ has the $DP$
property, then $E$ has the strict $DP_{p}$ property and the sequential
$DP_{(p,q)}$ property for all $p,q\in[1,\infty]$.

###### Proof 34

It is clear that $E$ is $q$-quasibarrelled for all $q\in[1,\infty]$. By
Theorem 12.1.4 of [25], the space $E^{\prime}_{\beta}$ is locally complete.
Now (i) and (iii) of Proposition 5.12 apply.

###### Corollary 5.14

Let $p,q\in[1,\infty]$, and let $E$ be a quasibarrelled, locally complete
space such that $E$ is weakly sequentially $p$-angelic and
$E^{\prime}_{\beta}$ is a weakly sequentially $q$-angelic space. Then the
following assertions are equivalent:

  1. (i)

$E$ has the $DP$ property;

  2. (ii)

$E$ has the strict $DP_{p}$ property;

  3. (iii)

$E$ has the sequential $DP_{(p,q)}$ property.

###### Proof 35

The equivalence (i)$\Leftrightarrow$(ii) follows from (ii) of Proposition
5.12. The implication (i)$\Rightarrow$(iii) follows from Corollary 5.13, and
(iii)$\Rightarrow$(i) follows from (iv) of Proposition 5.12.

If $p=q=\infty$, the Eberlein–Šmulian theorem and Corollary 5.14 immediately
imply the following classical Grothendieck result.

###### Corollary 5.15

A Banach space $E$ has the $DP$ property if and only if it has the sequential
$DP$ property.

## References

  * [1] F. Afkir, K. Bouras, A. Elbour, S. El Filali, Weak compactness of almost $L$-weakly and almost $M$-weakly compact operators, Quest. Math. 44 (2021), 1145–1154.
  * [2] A.A. Albanese, J. Bonet, W.J. Ricker, Grothendieck spaces with the Dunford–Pettis property, Positivity, 14 (2010), 145–164.
  * [3] F. Albiac, N.J. Kalton, Topics in Banach Space Theory, Grad. Texts in Math., vol. 233, Springer, New York, 2006.
  * [4] T. Banakh, S. Gabriyelyan, Locally convex properties of Baire type function spaces, Journal of Convex Analysis, 28 (2021), 803–818.
  * [5] T. Banakh, S. Gabriyelyan, On free locally convex spaces, Filomat 36:18 (2022), 6393–6401.
  * [6] J.M.F. Castillo, F. Sanchez, Dunford–Pettis-like properties of continuous vector function spaces, Revista Mat. Complut. 6:1 (1993), 43–59.
  * [7] J. Diestel, _Sequences and Series in Banach Spaces_ , GTM 92, Springer, 1984.
  * [8] M.B. Dehghani, S.M. Moshtaghioun, On the $p$-Schur property of Banach spaces, Ann. Funct. Anal. 9 (2018), 123–136.
  * [9] R.E. Edwards, _Functional Analysis_ , Reinhart and Winston, New York, 1965.
  * [10] K. El Fahri, J. H’michane, The relationship between almost Dunford–Pettis operators and almost limited operators, Quest. Math. 39 (2016), 487–496.
  * [11] M. Fabian, P. Habala, P. Hájek, V. Montesinos, J. Pelant, V. Zizler, _Banach space theory. The basis for linear and nonlinear analysis_ , Springer, New York, 2010.
  * [12] J.H. Fourie, E.D. Zeekoei, On $p$-convergent operators on Banach lattices, Acta Math. Sinica 34:5 (2018), 873–890.
  * [13] S. Gabriyelyan, The Mackey problem for free locally convex spaces, Forum Math. 30 (2018), 1339–1344.
  * [14] S. Gabriyelyan, Locally convex properties of free locally convex spaces, J. Math. Anal. Appl. 480 (2019) 123453.
  * [15] S. Gabriyelyan, Locally convex spaces and Schur type properties, Ann. Acad. Sci. Fenn. Math., 44 (2019), 363–378.
  * [16] S. Gabriyelyan, Maximally almost periodic groups and respecting properties, Descriptive Topology and Functional Analysis II, J.C. Ferrando (ed.), Springer Proceedings in Mathematics & Statistics, 286 (2019), 103–136.
  * [17] S. Gabriyelyan, Pełczyński’s type sets and Pełczyński’s geometrical properties of locally convex spaces, preprint.
  * [18] S. Gabriyelyan, J. , Dunford–Pettis type properties for function spaces, Revista Matemática Complutense 33 (2020), 871–884.
  * [19] I. Ghenciu, The $p$-Gelfand–Phillips property in spaces of operators and Dunford–Pettis like sets, Acta Math. Hungar. 155:2 (2018), 439–457.
  * [20] I. Ghenciu, P. Lewis, _The Dunford–Pettis property and the Gelfand–Phillips property, and $L$-sets_, Colloq. Math. 106 (2006), 311–324.
  * [21] L. Gillman, M. Jerison, _Rings of continuous functions_ , Van Nostrand, New York, 1960.
  * [22] W. Govaerts, A productive class of angelic spaces, J. London Math. Soc. 22 (1980), 355–364.
  * [23] A. Grothendieck, Sur les applications linéaires faiblement compactes d’espaces du type $C(K)$, Canad. J. Math. 5 (1953), 129–173.
  * [24] J. Horváth, _Topological Vector Spaces and Distributions, I_. Addison-Wesley, Reading, Mass, 1966.
  * [25] H. Jarchow, _Locally Convex Spaces_ , B.G. Teubner, Stuttgart, 1981.
  * [26] A.K. Karn, D. P. Sinha, An operator summability of sequences in Banach spaces, Glasgow Math. J. 56 (2014), 427–437.
  * [27] G. Köthe, _Topological vector spaces_ , Vol. I, Springer-Verlag, Berlin, 1969.
  * [28] L. Li, D. Chen, J.A. Chávez-Domínguez, Pełczyński’s property $(V^{\ast})$ of order $p$ and its quantification, Math. Nachrichten 291 (2018), 420–442.
  * [29] A.A. Markov, On free topological groups, Dokl. Akad. Nauk SSSR 31 (1941) 299–301.
  * [30] L. Narici, E. Beckenstein, _Topological vector spaces_ , Second Edition, CRC Press, New York, 2011.
  * [31] R. Nouira, D. Lhaimer, A. Elbour, Some results on almost $L$-weakly and almost $M$-weakly compact operators, Positivity 26 (2022), No 39.
  * [32] A. Pełczyński, Banach spaces on which every unconditionally converging operator is weakly compact, Bulletin de l’Académie Polonaise des Sciences, 10 (1962), 641–648.
  * [33] P. Pérez Carreras, J. Bonet, _Barrelled Locally Convex Spaces_ , North-Holland Mathematics Studies 131, North-Holland, Amsterdam, 1987.
  * [34] J.D. Pryce, A device of R.J. Whitley’s applied to pointwise compactness in spaces of continuous functions, Proc. London Math. Soc. 23 (1971), 532–546.
  * [35] H.H. Schaefer, _Topological Vector Spaces_ , Macmillan, New York, 1966.


* [29] N. E. Frangos, S. D. Vrontos & A. N. Yannacopoulos, Ruin probability at a given time for a model with liabilities of the fractional Brownian motion type: A partial differential equation approach., Scandinavian Actuarial Journal 2005, 4 (2005): 285–308.
  * [30] N. Ghoussoub & Y. Guo, On the partial differential equations of electrostatic MEMS devices II: dynamic case, Nonlinear Diff. Eqns. Appl. 15 (2008) 115–145.
  * [31] Y. Guo, Dynamical solutions of singular wave equations modeling electrostatic MEMS, SIAM J. Appl. Dyn. Syst., 9 (2010), 1135–1163.
  * [32] J.-S. Guo, B. Hu & C.-J. Wang, A nonlocal quenching problem arising in micro-electro mechanical systems, Quarterly Appl. Math., 67 (2009), 725–734.
  * [33] J.-S. Guo and N.I. Kavallaris, On a nonlocal parabolic problem arising in electrostatic MEMS control, Discrete Contin. Dyn. Syst., 32(5) (2012), 1723–1746.
  * [34] J.-S. Guo, N.I. Kavallaris, C.-J. Wang & C.-Y. Yu , The bifurcation diagram of a micro-electro mechanical system with Robin boundary condition, Hiroshima Mathematical Journal (to appear).
  * [35] H. Gimperlein, R.He & A. A. Lacey, Quenching for a semi-linear wave equation for micro-electro-mechanical systems, Proc. R. Soc. A. 478 2022049020220490 http://doi.org/10.1098/rspa.2022.0490
  * [36] S. Jarohs & T. Weth, On the strong maximum principle for nonlocal operators, Mathematische Zeitschrift 293 (2019), 81–111.
  * [37] I. Karatzas & S. Shreve, Brownian motion and stochastic calculus, Vol. 113 of Graduate Texts in Mathematics, Springer-Verlag, New York, 2nd edition, 1991.
  * [38] N. I. Kavallaris, Explosive solutions of a stochastic nonlocal reaction-diffusion equation arising in shear band formation, Math. Methods Appl. Sci. 38(16) (2015), 3564–3574.
  * [39] N. I. Kavallaris, Quenching solutions of a stochastic parabolic problem arising in electrostatic MEMS control, Math. Methods Appl. Sci. 41 (3) (2018), 1074–1082.
  * [40] N.I. Kavallaris, T. Miyasita and T. Suzuki, Touchdown and related problems in electrostatic MEMS device equation, Nonlinear Diff. Eqns. Appl., 15 (2008), 363–385.
  * [41] N. I. Kavallaris, A. A. Lacey, C. V. Nikolopoulos and D. E. Tzanetis, A hyperbolic non-local problem modelling MEMS technology, Rocky Mountain J. Math., 41 (2011), 505–534.
  * [42] N. I. Kavallaris, A. A. Lacey, C. V. Nikolopoulos and D. E. Tzanetis, On the quenching behaviour of a semilinear wave equation modelling MEMS technology, Discrete Contin. Dyn. Syst., 35 (2015), 1009–1037.
  * [43] N.I. Kavallaris, A.A. Lacey and C.V. Nikolopoulos, On the quenching of a nonlocal parabolic problem arising in electrostatic MEMS control, Nonlinear Analysis, 138 (2016), 189–206.
  * [44] N.I. Kavallaris & T. Suzuki, Non-Local Partial Differential Equations for Engineering and Biology: Mathematical Modeling and Analysis, Mathematics for Industry Vol. 31 Springer Nature 2018.
  * [45] N.I. Kavallaris & Y. Yan, Finite-time blow-up of a non-local stochastic parabolic problem, Stoch. Proc. Applications, 130(9), (2020), 5605–5635 doi.org/10.1016/j.spa.2020.04.002.
  * [46] F.C. Klebaner, Introduction to Stochastic Calculus with Applications, 2nd edition, Imperial College Press, London, 2005.
  * [47] L. Li, Z. Chen & T. Caraballo, Dynamics of a stochastic fractional nonlocal reaction–diffusion model driven by additive noise, Discrete and Continuous Dynamical Systems - S. doi: 10.3934/dcdss.2022179.
  * [48] J. López-Mimbela & Pérez, Global and nonglobal solutions of a system of nonautonomous semilinear equations with ultracontractive Lévy generators, J. Math.Anal.Appl. 423 (2015) 720–733.
  * [49] G.J. Lord , C.E. Powell & T. Shardlow An Introduction to Computational Stochastic PDEs, Cambridge University Press: Cambridge, UK, 2014.
  * [50] B. Maslowski & D. Nualart, Evolution equations driven by a fractional Brownian motion, J. Funct. Anal. 202 (2003), 277–305.
  * [51] Y. Mishura, Stochastic calculus for fractional Brownian motion and related processes, Springer Lecture Notes in Mathematics, 1929, Springer, Berlin (2008).
  * [52] T. Miyasita, On a nonlocal biharmonic MEMS equation with the Navier boundary condition, Sci. Math. Jpn. 80(2) (2017), 189–208.
  * [53] T. Miyasita, Convergence of solutions of a nonlocal biharmonic MEMS equation with the fringing field, J. Math. Anal. Appl. 454(1) (2017), 265–284.
  * [54] T. Miyasita, Global existence of radial solutions of a hyperbolic MEMS equation with nonlocal term, Differ. Equ. Appl. 7 (2) (2015), 169–186.
  * [55] D. Nualart,  The Malliavin Calculus and Related Topics, Springer Verlag, Berlin, 2006.
  * [56] D. Nualart & P. Vuillermot, Variational solutions for partial differential equations driven by a fractional noise, J. Funct. Anal. 232 (2006), 390–454.
  * [57] D. Nualart, Fractional Brownian motion: stochastic calculus and applications. In International Congress of Mathematicians (Vol. 3, pp. 1541-1562). Eur. Math. Soc., 2006.
  * [58] E. Nuxoll, BioMEMS in drug delivery, Adv. Drug Deliv. Rev. 65 (2013), 1611–1625.
  * [59] A. Pazy ,Semigroups of linear operators and applications to partial differential equations Applied Mathematical Sciences, 44 Springer-Verlag, New York, 1983.
  * [60] J.A. Pelesko & D.H. Bernstein, Modeling MEMS and NEMS, Chapman Hall and CRC Press, 2002.
  * [61] J.A. Pelesko & A.A. Triolo, Nonlocal problems in MEMS device control, J. Eng. Math. 41 (2001) 345–366.
  * [62] H. Schaefer, Banach lattices and positive operators, Springer 1974.
  * [63] M. Yor, Exponential Functionals of Brownian Motion and Related Processes,  Springer 2001.
  * [64] P. Yager, T. Edwards, E. Fu et al., Microfluidic diagnostic technologies for global public health, Nature 442 (2006), 412–418.
  * [65] M. Younis, MEMS Linear and Nonlinear Statics and Dynamics, Springer, New York, 2011.
  * [66] J. Xu & T. Caraballo, Dynamics of stochastic nonlocal partial differential equations, Eur. Phys. J. Plus 136, 849, (2021). https://doi.org/10.1140/epjp/s13360-021-01818-w.
  * [67] J. Xu & T. Caraballo, Dynamics of stochastic nonlocal reaction-diffusion equations driven by multiplicative noise, Analysis and Applications https://doi.org/10.1142/S0219530522500075.
  * [68] M. Zähle, Integration with respect to fractal functions and stochastic calculus II Math. Nachr. 225, (2001), 145–183.


# Nested solitons in two-field fuzzy dark matter

Hoang Nhan Luu 1,2, Philip Mocz 3, Mark Vogelsberger 2, Simon May 4,5, Josh
Borrow 2, S.-H. Henry Tye 1,6, Tom Broadhurst 7,8,9
1Department of Physics and Jockey Club Institute for Advanced Study, The Hong
Kong University of Science and Technology, Hong Kong
2Department of Physics, Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA
3Lawrence Livermore National Laboratory, 7000 East Ave, Livermore, CA 94550,
USA
4Perimeter Institute for Theoretical Physics, 31 Caroline Street North,
Waterloo, ON, N2L 2Y5, Canada
5Department of Physics, North Carolina State University, Raleigh, NC,
27695-8202, USA
6Department of Physics, Cornell University, Ithaca, NY 14853, USA
7University of the Basque Country UPV/EHU, Department of Theoretical Physics,
Bilbao, E-48080, Spain
8Donostia International Physics Center, Basque Country UPV/EHU, San Sebastian,
E-48080, Spain
9Ikerbasque, Basque Foundation for Science, Bilbao, E-48011, Spain E-mail:
<EMAIL_ADDRESS>0000-0001-9483-1099 0000-0001-6631-2566
0000-0001-8593-7692 0000-0002-2781-6304 0000-0002-1327-1921
0000-0002-4386-0102 0000-0002-8785-8979

###### Abstract

Dark matter as scalar particles consisting of multiple species is well
motivated in string theory where axion fields are ubiquitous. A two-field
fuzzy dark matter (FDM) model features two species of ultralight axion
particles with different masses, $m_{1}\neq m_{2}$, which is extended from the
standard one-field model with $m_{a}\sim 10^{-22}\,{\rm eV}$. Here we perform
numerical simulations to explore the properties of two-field FDM haloes. We
find that the central soliton has a nested structure when $m_{2}\gg m_{1}$,
which is distinguishable from the generic flat-core soliton in one-field
haloes. However, the formation of this nested soliton is subject to many
factors, including the density fraction and mass ratio of the two fields.
Finally, we study non-linear structure formation in two-field cosmological
simulations with self-consistent initial conditions and find that the small-
scale structure in two-field cosmology is also distinct from the one-field
model in terms of DM halo counts and soliton formation time.

###### keywords:

fuzzy dark matter – numerical simulation – cosmology

††pubyear: 2023††pagerange: Nested solitons in two-field fuzzy dark
matter–Nested solitons in two-field fuzzy dark matter

##  1 Introduction

The nature of dark matter (DM) remains elusive since its first evidence was
found almost a century ago by Fritz Zwicky (Zwicky, 1937). The cold dark
matter (CDM) hypothesis, where DM is made of non-relativistic, collision-less
particles, has successfully explained the anisotropies of the Cosmic Microwave
Background and the large-scale structure of the universe; see Aghanim et al.
(2020) and Alam et al. (2021) for the most recent observational data from
Planck and eBOSS, see Vogelsberger et al. (2020) for a review on state-of-the-
art cosmological simulations of galaxy formation. From a microscopic
perspective, the models of weakly interacting massive particles (WIMPs), where
DM particles are presumably heavy and have interactions of electroweak scales
with Standard Model particles, have been the leading candidates for CDM in the
past decades. However, CDM is also facing several challenges due to
inconsistencies between $N$-body simulations and observations on small scales,
including the “core–cusp problem” (de Blok, 2010), the “too-big-to-fail
problem” (Garrison-Kimmel et al., 2014; Boylan-Kolchin et al., 2011), the
diversity problem of rotation curves (Oman et al., 2015) – together with the
null results of WIMP signals in many laboratory searches (Roszkowski et al.,
2018).

Even though baryonic feedback, e. g. from supernova explosions or stellar
winds (Governato et al., 2012; Di Cintio et al., 2014), may partially explain
these CDM small-scale crises, other dark matter models, e. g. warm dark matter
(Lovell et al., 2012), fuzzy dark matter (Hu et al., 2000), or self-
interacting dark matter (Spergel & Steinhardt, 2000), have emerged as
promising alternatives to CDM. Fuzzy dark matter (FDM), also known as wave
dark matter (Schive et al., 2014a; Hui, 2021) or ultralight axions (Marsh,
2016; Hui et al., 2017; Ferreira, 2021), is a particularly interesting DM
candidate to resolve the small-scale problems but still maintain CDM
predictions on large scales. FDM is typically comprised of extremely light
bosons with mass $m_{a}\sim 10^{-22}\,\mathrm{eV}$, which have de Broglie
wavelengths of astronomical scales, $\lambda_{\rm
dB}\propto(m_{a}v)^{-1}\sim\mathrm{kpc}$.

Therefore, in FDM cosmology, bound structures below kiloparsec scales are
washed out due to an effective “quantum” pressure acting against gravity. For
the same reason, FDM haloes form a central flat core called soliton (Schive et
al., 2014a; Schwabe et al., 2016; Mocz et al., 2017) rather than the cuspy
Navarro-Frenk-White (NFW) profile of CDM (Navarro et al., 1996). A typical
density profile of FDM haloes also features sophisticated interference
granules of density fluctuations (Schive et al., 2014a; Mocz et al., 2017)
around the central core. However, the simplest model of FDM, where $m_{a}$ is
the only parameter, has been tightly constrained by observational data from
the Lyman-$\alpha$ forest (Armengaud et al., 2017; Iršič et al., 2017),
ultrafaint dwarf galaxies (Hayashi et al., 2021; Dalal & Kravtsov, 2022), and
strong lensing (Powell et al., 2023). The current consensus is that a single
FDM species of mass $10^{-22}\,\mathrm{eV}$ is unlikely to account for the
total dark matter budget while still addressing the small-scale issues of CDM
at the same time.

Several extensions of the most simplified FDM model have been widely explored
in literature. A few notable studies are self-interacting FDM (Mocz et al.,
2023) where the quartic self-coupling term is included, vector fuzzy dark
matter (Amin et al., 2022) where FDM is a higher-spin field, and mixed cold
and fuzzy dark matter (Schwabe et al., 2020). From a theoretical perspective,
the ultralight mass of the FDM field can be generated in the context of string
compactifications. It is expected that there may exist more than one species
of ultralight axions with masses spanning many orders of magnitudes as in the
string axiverse scenario (Svrcek & Witten, 2006; Arvanitaki et al., 2010).
Thus, the multiple-field FDM model is naturally motivated by string theory.
Along this direction, the multiple-field FDM model with a specific focus on
the two-field case has been proposed in Luu et al. (2020).

We previously found that the ground state of a time-independent two-field FDM
system is a nested configuration of two solitons, each from one component
field, which is distinct from the usual soliton in the one-field FDM model.
Ever since our work, there has been increasing interest in the two-field FDM
model, starting with Eby et al. (2020) and Guo et al. (2021) where the authors
studied the semi-analytical profile of two-field solitons. More recently,
Huang et al. (2023) performed the first two-field cosmological simulations
with CDM initial conditions; Gosenca et al. (2023) investigated stellar
heating in two- and higher-field FDM haloes; Glennon et al. (2023) simulated
soliton collision of two-field FDM with self-interactions and inter-field
interactions; Jain et al. (2023) studied Bose-Einstein condensation in the
kinetic regime. These studies showed that the phenomenology of the two-field
FDM model can be complex, novel, and compelling. Importantly, it has the
potential to retrieve the initial promises of the one-field FDM model.

In this work, we carry out non-cosmological and cosmological simulations of
the two-field FDM model to study the final state of virialized haloes. Our
main goal is to verify whether the nested soliton, previously predicted in Luu
et al. (2020) can emerge and become stable in numerical simulations. The paper
is organized as follows. In Sec. 2, we review the theoretical background of
the two-field Schrödinger–Poisson equations. Sec. 3 describes the numerical
method for our simulations and ground-state solver in detail. Sec. 4 and Sec.
5 set up several simulations and then report the main results and new findings
of our paper. Finally, conclusions and potential directions for future work
are discussed in Sec. 6.

##  2 Theoretical Background

In this section, we introduce the time-dependent and time-independent
Schrödinger-Poisson equations of the two-field FDM model in Sec. 2.1 and Sec.
2.2, respectively.

###  2.1 Two-field Schrödinger–Poisson equations

In the non-relativistic limit, the minimal (non-self-interacting) system of
two-field FDM with different masses is governed by gravitationally coupled
Schrödinger–Poisson equations. In physical coordinates in the Newtonian gauge,
these equations are given by

$\displaystyle i\hbar\dfrac{\partial\psi_{1}}{\partial
t}+\dfrac{3}{2}H\psi_{1}=-\dfrac{\hbar^{2}}{2m_{1}}\nabla^{2}\psi_{1}+m_{1}\Phi\psi_{1},$
(1) $\displaystyle i\hbar\dfrac{\partial\psi_{2}}{\partial
t}+\dfrac{3}{2}H\psi_{2}=-\dfrac{\hbar^{2}}{2m_{2}}\nabla^{2}\psi_{2}+m_{2}\Phi\psi_{2},$
$\displaystyle\nabla^{2}\Phi=4\pi
G\left(|\psi_{1}|^{2}+|\psi_{2}|^{2}-\bar{\rho}\right),$

where $\Phi$ is the gravitational potential, $\psi_{1}$ and $\psi_{2}$ are the
wavefunctions of the FDM field with masses $m_{1}$ and $m_{2}$, respectively.
The first field $\psi_{1}$ is always chosen as the lighter FDM field in our
convention throughout this study. The density of each field is given by
$\rho_{i}=|\psi_{i}|^{2}$, where $i=1,2$. The total density is obtained by
summing individual densities, $\rho=\rho_{1}+\rho_{2}$, while $\bar{\rho}$ is
the corresponding mean of the total density. For convenience, we define a
parameter to quantify the ratio of the heavy-field mean density to the total
density as

$\displaystyle\beta_{2}\equiv\dfrac{\bar{\rho}_{2}}{\bar{\rho}_{1}+\bar{\rho}_{2}}=\dfrac{M_{1}}{M_{1}+M_{2}}\quad{\rm
where}\quad\bar{\rho}_{i}=\dfrac{M_{i}}{L_{\rm box}^{3}},$ (2)

and $M_{i}$ denotes the total mass of each field enclosed within a volume with
length $L_{\rm box}$ along each spatial dimension. As the two-field
Schrödinger–Poisson equations (1) are highly non-linear due to gravitational
couplings, numerical simulations are mandatory to study the system’s dynamics.
Simulation techniques for the Schrödinger-Poisson equations will be reviewed
later in Sec. 3.1.

###  2.2 Time-independent Schrödinger–Poisson equations

Some prominent features of FDM haloes, such as the stable ground-state
solution, can be found in an approximate limit. Once the field dynamics are
stabilized, the time-dependent part of the FDM fields can be factorized out,
$\psi_{i}(\mathbf{x},t)=\exp(-iE_{i}t/\hbar)\,\psi_{i}(\mathbf{x})$, where
$E_{i}$ is the “eigenenergy” of each field. Let us also assume that the system
is spherically symmetric,
$\psi_{i}(\mathbf{x})=\psi_{i}(r),\Phi(\mathbf{x})=\Phi(r)$, and omit the
Hubble flow as well as the mean density. The time-independent
Schrödinger–Poisson equations are then derived as follows

$\displaystyle\dfrac{\partial^{2}\psi_{1}}{\partial
r^{2}}=-\dfrac{2}{r}\dfrac{\partial\psi_{1}}{\partial
r}+\dfrac{2m_{1}}{\hbar^{2}}(m_{1}\Phi-E_{1})\psi_{1},$ (3)
$\displaystyle\dfrac{\partial^{2}\psi_{2}}{\partial
r^{2}}=-\dfrac{2}{r}\dfrac{\partial\psi_{2}}{\partial
r}+\dfrac{2m_{2}}{\hbar^{2}}(m_{2}\Phi-E_{2})\psi_{2},$
$\displaystyle\dfrac{\partial^{2}\Phi}{\partial
r^{2}}=-\dfrac{2}{r}\dfrac{\partial\Phi}{\partial r}+4\pi
G\left(|\psi_{1}|^{2}+|\psi_{2}|^{2}\right).$

The ground-state solution of this system is called a soliton, which is a cored
density profile that satisfies the following boundary conditions

$\displaystyle\psi_{1}(r=0)=\psi_{1}(0),\quad\psi_{2}(r=0)=\psi_{2}(0),$ (4)
$\displaystyle\psi_{1}(r\rightarrow\infty)\rightarrow
0,\quad\psi_{2}(r\rightarrow\infty)\rightarrow 0,$
$\displaystyle\dfrac{\partial\psi_{1}}{\partial
r}(r=0)=\dfrac{\partial\psi_{2}}{\partial r}(r=0)=0,$

and $\Phi(0)$ can be arbitrarily chosen. Integrating the density of each field
yields the corresponding total mass of the soliton

$\displaystyle
M_{s,i}\equiv\int\rho_{i}(\mathbf{x})d^{3}\mathbf{x}=4\pi\int\rho_{i}(r)r^{2}dr.$
(5)

In practice, the set of equations (3) can be further simplified thanks to the
scaling symmetry of the Schrödinger–Poisson equations. The Schrödinger–Poisson
equations are unchanged under a transformation of variables with parameter
$\lambda$ such that (Ji & Sin, 1994)

$t\rightarrow\lambda^{-2}t,\quad\mathbf{x}\rightarrow\lambda^{-1}\mathbf{x},\quad\psi_{i}\rightarrow\lambda^{2}\psi_{i},\quad\Phi\rightarrow\lambda^{2}\Phi.$
(6)

This scaling relation suggests that Eqs. (3) can be re-scaled into

$\displaystyle\dfrac{\partial^{2}\tilde{\psi}_{1}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}_{1}}{\partial\tilde{r}}+\left(\dfrac{m_{1}}{m}\right)^{2}(\tilde{\Phi}-\tilde{E}_{1})\tilde{\psi}_{1},$
(7)
$\displaystyle\dfrac{\partial^{2}\tilde{\psi}_{2}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}_{2}}{\partial\tilde{r}}+\left(\dfrac{m_{2}}{m}\right)^{2}(\tilde{\Phi}-\tilde{E}_{2})\tilde{\psi}_{2},$
$\displaystyle\dfrac{\partial^{2}\tilde{\Phi}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\Phi}}{\partial\tilde{r}}+|\tilde{\psi}_{1}|^{2}+|\tilde{\psi}_{2}|^{2},$

where a new set of dimensionless variables has been introduced with a tilde
sign. They are related to the original variables by

$\displaystyle r=\dfrac{\hbar^{2}}{2m^{2}GM}\tilde{r},\quad
E_{i}=\dfrac{2G^{2}M^{2}m^{2}m_{i}}{\hbar^{2}}\tilde{E}_{i},$ (8)
$\displaystyle\psi_{i}=\sqrt{\dfrac{2G^{3}M^{4}m^{6}}{\pi\hbar^{6}}}\tilde{\psi}_{i},\quad\Phi=\dfrac{2G^{2}M^{2}m^{2}}{\hbar^{2}}\tilde{\Phi}.$

Here, $M$ is a free scaling parameter, corresponding to $\lambda$ in (6). $M$
is the normalization factor of the soliton mass,
$\int|\tilde{\psi}_{i}|^{2}\tilde{r}^{2}d\tilde{r}=M_{s,i}/M$, and $m$
represents an effective FDM mass (which can, however, be freely chosen), which
is usually set to $m=\sqrt{m_{1}m_{2}}$ in our computations.

##  3 Numerical method

In this section, we review the pseudo-spectral method used for our simulations
in Sec. 3.1, then provide detailed guidelines to solve for the ground state of
two-field equations in Sec. 3.2.

###  3.1 Pseudo-spectral solver

The dynamical two-field FDM system in Eqs. (1) can be simulated using the
pseudo-spectral method extended for two fields with different masses. The
simulation volume is a cubic box discretized on a fixed Cartesian grid with
uniform resolution. Periodic boundary conditions are implicitly imposed in the
spectral solver. The wavefunction of each field is evolved unitarily for each
timestep $\Delta t$ as follows:

  * •

Step 1: Evolve $\psi_{1}$ and $\psi_{2}$ with the potential “kick” operator by
half a timestep

$\displaystyle\psi_{i}\longleftarrow\exp\left(-i\dfrac{\Delta
t}{2}\dfrac{m_{i}}{\hbar}\Phi\right)\psi_{i},$ (9)

  * •

Step 2: Evolve $\psi_{1}$ and $\psi_{2}$ with the kinetic “drift” operator by
a full timestep in Fourier space

$\displaystyle\psi_{i}\longleftarrow{\rm ifft}\left[\exp\left(-i\Delta
t\dfrac{\hbar}{2m_{i}}k^{2}\right){\rm fft}[\psi_{i}]\right],$ (10)

  * •

Step 3: Calculate the total density and update the potential (also in Fourier
space)

$\displaystyle\rho=|\psi_{1}|^{2}+|\psi_{2}|^{2},\quad\Phi\longleftarrow{\rm
ifft}\left[-\dfrac{1}{k^{2}}{\rm fft}\left[4\pi
G(\rho-\bar{\rho})\right]\right],$ (11)

  * •

Step 4: Repeat Step 1 with the updated potential.

Here “fft” and “ifft” represent the Fast Fourier Transform and its inverse. We
employ the FDM module of the magneto-hydrodynamics code AREPO (Springel, 2010;
Pakmor et al., 2016; Weinberger et al., 2020) originally implemented by the
authors of Mocz et al. (2017) and modify it accordingly for two-field FDM
simulations.

Since the pseudo-spectral algorithm converges at second order in time and
exponentially in space, the heavy-field (smaller) soliton is fully resolved
even when it overlaps only a few cells along each dimension. On the one hand,
the spacing requirement of the grid, $\Delta x$, needs to cover the finest
features of the heaviest field (Mocz et al., 2018; May & Springel, 2021), i.
e.,

$\displaystyle\Delta
x<\min\left[\dfrac{\pi\hbar}{v_{1,\max}m_{1}},\dfrac{\pi\hbar}{v_{2,\max}m_{2}}\right],$
(12)

where $v_{i,{\rm max}}$ denotes the maximum velocity of the corresponding
field. On the other hand, the timestep criterion needs to be satisfied for the
lightest field (Schwabe et al., 2016; May & Springel, 2021), i. e.,

$\displaystyle\Delta t<\min\left[\dfrac{4m_{1}}{3\pi\hbar}\Delta
x^{2},\dfrac{2\pi\hbar}{m_{2}|\Phi_{\rm max}|}\right].$ (13)

These two criteria often require simulations with very high temporal and
spatial resolution to resolve the inner core of both fields. Therefore,
multiple-field simulations are much more computationally challenging than one-
field simulations.

###  3.2 Two-field ground state solver

The dimensionless equations (7) are particularly simple to solve and useful to
scale up to different physical systems. For instance, in the one-field FDM
model, the scaling parameter $m$ can be chosen as the FDM mass $m_{a}$ and we
obtain a completely scale-free set of equations

$\displaystyle\dfrac{\partial^{2}\tilde{\psi}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\psi}}{\partial\tilde{r}}+(\tilde{\Phi}-\tilde{E})\tilde{\psi},\qquad\dfrac{\partial^{2}\tilde{\Phi}}{\partial\tilde{r}^{2}}=-\dfrac{2}{\tilde{r}}\dfrac{\partial\tilde{\Phi}}{\partial\tilde{r}}+|\tilde{\psi}|^{2},$
(14)

where the field indices have been omitted. As the above equations are
independent of FDM parameters, the ground state corresponds to a unique
solution of $\tilde{\psi}$. This solution can be freely scaled up (via the
parameter $M$), yielding the soliton profile approximated by (Schive et al.,
2014a)

$\displaystyle|\tilde{\psi}|^{2}(r)\propto\rho(r)\simeq 1.9\times
10^{9}\dfrac{{\rm
M}_{\odot}}{\mathrm{kpc}^{3}}\dfrac{(10^{-22}\,\mathrm{eV}/m_{a})^{2}(\mathrm{kpc}/r_{\mathrm{c}})^{4}}{\left[1+0.091(r/r_{\mathrm{c}})^{2}\right]^{8}},$
(15)

where $r_{\mathrm{c}}$ is the core radius defined as the radius where the core
density drops by half, $\rho(r_{\mathrm{c}})=\rho(0)/2$. The total mass of
each soliton is given approximately by $M_{s}\simeq
M_{s,0}\,(r_{\mathrm{c}}/\mathrm{kpc})\,(m_{a}/10^{-22}\,\mathrm{eV})^{2}$
where $M_{s,0}\equiv 2.2\times 10^{8}\,{\rm M}_{\odot}$. We note that this
profile only describes the ground-state solution of a one-field FDM system. In
contrast, the ground-state solution of a two-field FDM system described by
Eqs. (3) is not well approximated by this profile. As such, we will refer to
the one-field and two-field ground states as one-field soliton and two-field
soliton to distinguish them from now on.

For the two-field equations in (7), we cannot eliminate the FDM masses with
scaling relations. Thus, the ground-state solution is dependent on at least
two parameters: the FDM mass ratio $m_{2}/m_{1}$ and the central wavefunction
ratio defined as

$\displaystyle\alpha_{2}\equiv\sqrt{\dfrac{\rho_{2}(0)}{\rho_{1}(0)}}=\dfrac{\psi_{2}(0)}{\psi_{1}(0)}=\dfrac{\tilde{\psi}_{2}(0)}{\tilde{\psi}_{1}(0)}$
(16)

since $\tilde{\psi}_{1}(0)$ can always be fixed to a constant, e. g.,
$\tilde{\psi}_{1}(0)=1$. We apply the shooting method to solve these
equations: (i) pick a not-too-large endpoint $\tilde{r}_{\rm end}$ and
arbitrarily choose some initial values for $\tilde{E}_{i}$; (ii) solve (7) as
an initial condition problem where
$\tilde{\psi}_{1}(0)=1,\tilde{\psi}_{2}(0)=\alpha_{2}$; (iii) adjust
$\tilde{E}_{1}$ and $\tilde{E}_{2}$ simultaneously such that
$\tilde{\psi}_{i}(\tilde{r}_{\rm end})\rightarrow 0$ with multi-dimensional
Newton-Raphson algorithm; (iv) repeat step (i) to (iii) until $\tilde{r}_{\rm
end}$ is large enough. There may exist several excited states of different
energies satisfying regular boundary conditions. Thus, we must impose the
condition $\tilde{\psi}_{i}(r)>0~{}\forall r$ to guarantee the ground-state
solution. We list in Tab. 1 some eigenenergies corresponding to the ground
state of the two-field systems for different values of $\alpha_{2}$ and
$\beta_{2}$ as an example.

Table 1: Dimensionless eigenenergies of two-field solitons for different density and central wavefunction ratios $\beta_{2}$ and $\alpha_{2}$ ($m_{1}=8\times 10^{-23}\,\mathrm{eV}$, $m_{2}=1.6\times 10^{-22}\,\mathrm{eV}$; same as Scenario B in Tab. 2). The first and second entries of each row correspond to $\tilde{E}_{1}$ and $\tilde{E}_{2}$, respectively. $\beta_{2}$ | $\alpha_{2}$ |  $\tilde{E}_{i}$ (exact) |  $\tilde{E}_{i}$ (approximate)
---|---|---|---
0.3 | 0.68 | 1.4761993842 | 1.4606887253
0.8739143554 | 0.9756256155
0.5 | 2 | 2.3898610754 | 2.3449988097
1.5042486006 | 1.5686446461
0.7 | 7.75 | 7.7300571770 | 7.6970072599
5.0941543699 | 5.1154890489

Although these steps are straightforward, we cannot apply them to solve the
two-field system with a high mass ratio, $m_{2}/m_{1}\gg 1$. In these limits,
we need to tune $\tilde{E}_{2}$ extremely precisely compared to
$\tilde{E}_{1}$ so that the solution does not blow up. As $\tilde{E}_{2}$
reaches floating-point precision, the shooting method stops converging without
arbitrary-precision arithmetic. The high-mass limit is, however, an ideal
condition to apply the flat density approximation, where the density
distribution of the light field is treated as constant with respect to the
heavy field assuming that their cores are very different in sizes,
$\tilde{r}_{c,2}\ll\tilde{r}_{c,1}$. This allows the use of what are
effectively the one-field equations (14), with an additional constant external
source. The approximation procedure also follows four steps:

  * •

Step 1: Solve for the heavy wavefunction with a constant density distribution
of the light field as an external source, setting $m=m_{1}$ and $M=\hat{M}$
(arbitrary), i. e.

$\begin{aligned}
&\dfrac{\partial^{2}\hat{\psi}_{2}}{\partial\hat{r}^{2}}=-\dfrac{2}{\hat{r}}\dfrac{\partial\hat{\psi}_{2}}{\partial\hat{r}}+(\hat{\Phi}-\hat{E}_{2})\hat{\psi}_{2}\\\
&\dfrac{\partial^{2}\hat{\Phi}}{\partial\hat{r}^{2}}=-\dfrac{2}{\hat{r}}\dfrac{\partial\hat{\Phi}}{\partial\hat{r}}+\left|\hat{\psi}_{2}\right|^{2}+\dfrac{1}{\alpha_{2}^{2}}\end{aligned}\;\;{\rm
with}\;\;\hat{\psi}_{2}(0)=1,$ (17)

where we have adopted an approximation
$\hat{\psi}_{1}(\hat{r})\simeq\alpha_{2}^{-1}$.

  * •

Step 2: Rescale the variables by setting $m=m_{2}$ and
$M=\alpha_{2}^{-1/2}(m_{2}/m_{1})^{3/2}\hat{M}$, resulting in

$\displaystyle\psi^{\prime}_{2}(r^{\prime})=\alpha_{2}\hat{\psi}_{2}(\hat{r}),\;\;r^{\prime}=\alpha_{2}^{-1/2}\left(m_{1}/m_{2}\right)^{1/2}\hat{r}.$
(18)

  * •

Step 3: Solve for the light wavefunction with an additional source of the
heavy field just found:

$\begin{aligned} &\dfrac{\partial^{2}\psi^{\prime}_{1}}{\partial r^{\prime
2}}=-\dfrac{2}{r^{\prime}}\dfrac{\partial\psi^{\prime}_{1}}{\partial
r^{\prime}}+(\Phi^{\prime}-E^{\prime}_{1})\psi^{\prime}_{1}\\\
&\dfrac{\partial^{2}\Phi^{\prime}}{\partial r^{\prime
2}}=-\dfrac{2}{r^{\prime}}\dfrac{\partial\Phi^{\prime}}{\partial
r^{\prime}}+|\psi^{\prime}_{1}|^{2}+\left|\psi^{\prime}_{2}\right|^{2}\end{aligned}\;\;{\rm
with}\;\;\psi^{\prime}_{1}(0)=1.$ (19)

  * •

Step 4: To return to our “usual choice” of $m=\sqrt{m_{1}m_{2}}$, rescale the
variables again by additionally setting $M=(m_{1}/m_{2})^{3/4}M^{\prime}$:

$\displaystyle\tilde{\psi}_{1}(\tilde{r})=\psi^{\prime}_{1}(r^{\prime}),\;\;\tilde{\psi}_{2}(\tilde{r})=\psi^{\prime}_{2}(r^{\prime}),\;\;\tilde{r}=\left(m_{2}/m_{1}\right)^{1/4}r^{\prime}.$
(20)

We then obtain an approximate solution of Eqs. (7).

Figure 1: Dimensionless wavefunctions of the ground state with
$m_{2}/m_{1}=2,\beta_{2}=0.5$ and $\alpha_{2}=2$ (see Tab. 1) solved by the
shooting (solid) and flat density approximation (dashed) methods.

The eigenenergies are not crucial in this method, but they can be
straightforwardly retrieved from the scaling relations:

$\displaystyle
E^{\prime}_{2}=\left(m_{1}/m_{2}\right)\alpha_{2}\hat{E}_{2},\;\;\tilde{E}_{i}=\left(m_{2}/m_{1}\right)^{1/2}E^{\prime}_{i}.$
(21)

These values provide a reasonable estimate of the exact energies (see Tab. 1)
and they can serve as an initial guess for the precise shooting method. Figure
1 compares the ground-state wavefunctions solved by both methods for Scenario
B with $\beta_{2}=0.5$. Even though the mass ratio is not high in this case,
the flat-density approximation (dashed curves) yields virtually
indistinguishable profiles from the precise ones (solid curves).

Regardless of the method used, we can restore the two-field soliton density by
computing the scaling factor $M$ based on the relation between dimensionless
and physical variables. More specifically, we require

$\displaystyle
M=\left[\dfrac{\pi\hbar^{6}(1+\alpha_{2}^{2})}{2G^{3}m_{1}^{3}m_{2}^{3}\rho(0)}\right]^{1/4}.$
(22)

to fit a two-field FDM soliton with central density
$\rho(0)=|\psi_{1}(0)|^{2}+|\psi_{2}(0)|^{2}$. Once $M$ is determined, the
physical radius and the physical wavefunction of both fields can be derived
according to (8).

Figure 2: The top panels show slices through a snapshot of the total, light, and heavy fields (from left to right) of a virialized two-field FDM halo at $t\sim 2.6\,\mathrm{Gyr}$. The bottom panels are the corresponding radial profiles of each field. The FDM mass of the heavy field is five times that of the light field, $m_{2}=5m_{1}=4\times 10^{-22}\,\mathrm{eV}$, while its abundance is four times less than the light field, $\beta_{2}=0.25$, as indicated in the figures. We distinguish the core radius of the best-fit one-field soliton (dashed), denoted as FDM, and two-field soliton (solid), denoted as 2FDM, as the former is overshooting the simulation profile of the light field. In contrast, the core radius of the heavy field is almost identical in both profiles (with dashed and solid green curves almost overlapping). The orange-dashed curve denotes the best fit of a NFW profile with the scale density $\rho_{s,0}\sim 9\times 10^{7}\,{\rm M}_{\odot}\mathrm{kpc}^{-3}$ and the scale radius $r_{s}\sim 1.85\,\mathrm{kpc}$ to the outer halo of the total field. Figure 3: Radial density profiles from several non-cosmological simulations calculated at $t\sim 2.4\,\mathrm{Gyr}$. (Upper row) Profiles of the light, $m_{1}=8\times 10^{-23}\,\mathrm{eV}$, and heavy field, $m_{2}=1.6\times 10^{-22}\,\mathrm{eV}$, in three simulations with varying mass fractions $\beta_{2}\equiv M_{2}/M_{\rm total}$ (Scenario B). The inset of the last panel shows the goodness-of-fit of the heavy-field profile (see definition in (24)) in other simulations with different $\beta_{2}$ and at the same moment. (Lower row) Profiles of the light and heavy fields in similar simulations with varying heavy-field mass $m_{2}$ and a fixed mass fraction $\beta_{2}=0.3$ (Scenario C). The heavy field cannot reach a stable ground state in some simulations where neither the two-field soliton nor the one-field soliton fits the simulation profile. Table 2: Setup of the soliton collision simulations in Sec. 4. Here the mass of the light field is $m_{1}=8\times 10^{-23}\,\mathrm{eV}$. The last column shows whether the two-field soliton is found at the end of each simulation. Scenario | $\beta_{2}$ | $M_{\rm total}$ | $m_{2}/m_{1}$ | $r_{c,1}$ | $r_{c,2}$ | Solion
---|---|---|---|---|---|---
$(10^{10}\,{\rm M}_{\odot})$ | (kpc) | (kpc) | formed
A | 0.25 | $1.03$ | 5 | 2 | 0.08 | Yes
B | 0.3 | $0.3$ | 2 | 1.64 | 0.95 | Yes
0.5 | 2.29 | 0.57 | Yes
0.7 | 3.82 | 0.41 | No
C | 0.3 | $0.3$ | 5 | 1.64 | 0.15 | No
6 | 0.11 | No
7 | 0.08 | Yes

##  4 Soliton collision with two-field FDM

We perform non-cosmological simulations to investigate the features of two-
field haloes. Initially, the soliton cores with density profiles given by (15)
are randomly placed in a box of length $L_{\rm box}=50\,\mathrm{kpc}$ with
zero velocities. We then let them gravitationally interact, merge, and form a
virialized halo. The final state is captured at time
$t=2.4{-}2.6\,\mathrm{Gyr}$, when the two-field halo has become relatively
stationary. Table 2 summarizes the setup of the soliton collision simulations
in the following sections.

###  4.1 Nested soliton profile

Firstly, we collide 15 solitons of the light field, $m_{1}=8\times
10^{-23}\,\mathrm{eV}$, with 45 solitons of the heavy field with
$m_{2}=4\times 10^{-22}\,\mathrm{eV}$ in a simulation box with a resolution of
$1024^{3}$ cells. The core radius of the light-field solitons is manually set
to $r_{c,1}=2\,\mathrm{kpc}$. The core radius of the heavy-field solitons is
inferred from $r_{c,1}$ to match the mass of the light-field ones,
$r_{c,2}=r_{c,1}(m_{1}/m_{2})^{2}$, so that the mean density of the heavy
field accounts for 25% of the total density in the system, i. e.,
$\beta_{2}=0.25$. We refer to this simulation as Scenario A.

Figure 2 shows slices of the density distribution and corresponding radial
profiles that are averaged over angular directions in Scenario A. The origin
is chosen at the location with maximum total density. As expected, the halo of
the total field and the individual components show a central core surrounded
by interference granules in the simulation. Most notably, the core of the
total field exhibits a nested structure of two constituent solitons with a
transition around $0.16\,\mathrm{kpc}$. This nested profile is fitted well
with the ground-state soliton (red solid curve) of the time-independent two-
field Schrödinger–Poisson equations (3). Similarly, the predicted two-field
soliton profiles of the light field (blue solid curve) and heavy field (green
solid curve) are in excellent agreement with the inner cores of the associated
simulation fields, although the heavy-field soliton is only resolved by three
cells. The outer halo follows an NFW-like profile, similar to the one-field
FDM model, as expected.

The one-field profile of the light field (blue-dashed curve) with a larger
core radius does not fit the simulation data, unlike its two-field
counterpart. As this profile is determined by solving the Schrödinger–Poisson
equations for the light field only, the corresponding core radius is
noticeably wider due to the absence of gravitational attraction from the
heavy-field soliton. On the other hand, the one-field and two-field soliton
profiles of the heavy field are virtually identical since it only experiences
gravitational effects from an almost flat and low-density distribution of the
light field, i. e., $\rho_{2}(0)\gg\rho_{1}(0)$. The one-field profile of the
total field (red dashed curve), defined as the sum of two component one-field
solitons, is not compatible with the simulation data in a similar manner to
the light field.

In summary, the two-field halo profile contains a nested soliton with a
visible transition from one core to another and a second transition from the
soliton to the remaining NFW-like halo. Yet, the soliton-soliton transition is
not always visible due to time-dependent fluctuations of the heavy-field
granules which induce a slight density excess over the light-field soliton.
The soliton-halo transition is found at $r\sim 3-3.5r_{c,i}$ for the light and
heavy field, which agrees with previous studies (Schive et al., 2014b; Mocz et
al., 2017). As for the total field, although the soliton-halo transition can
be seen around $0.5\,\mathrm{kpc}$, it is generally ill-defined when the
light-field soliton blends with the interference granules of the heavy field
as mentioned above.

###  4.2 Soliton formation and instabilities

Even though we have seen the existence of a stable nested soliton as the
ground-state solution of two-field Schrödinger–Poisson equations, it may not
always be able to form with any value of the mass fraction. Recently, one
particular cosmological simulation of two-field FDM (Huang et al., 2023) with
$m_{2}=3m_{1}=3\times 10^{-22}\,\mathrm{eV}$ and $\beta=0.25$ has pointed out
that the heavy-field soliton cannot form under deep gravitational potential
fluctuations induced by the light field. We are, therefore, motivated to study
the formation of two-field solitons in Scenario B as follows.

We collide 10 solitons of each field with $m_{1}=8\times
10^{-23}\,\mathrm{eV}$ and $m_{2}=1.6\times 10^{-22}\,\mathrm{eV}$ in a
simulation box with a resolution of $512^{3}$ cells. The heavy-field solitons
are placed concentrically with the light-field ones to enhance collision
efficiency. The total mass is fixed, $M_{\rm total}=M_{1}+M_{2}=3\times
10^{9}\,{\rm M}_{\odot}$. The heavy-field density ratio varies from 30% to
70%, $\beta_{2}=\\{0.3,0.5,0.7\\}$, in three simulations. Thus, the core radii
of the initial solitons need to be adjusted in each simulation to satisfy the
density ratio required

$\displaystyle r_{c,1}$
$\displaystyle\simeq(1-\beta_{2})^{-1}(N_{s,1}M_{s,0}/M_{\rm
total})(10^{-22}\,\mathrm{eV}/m_{1})^{2},$ (23) $\displaystyle r_{c,2}$
$\displaystyle\simeq\beta_{2}^{-1}(N_{s,2}M_{s,0}/M_{\rm
total})(10^{-22}\,\mathrm{eV}/m_{2})^{2},$

where $N_{s,i}$ denotes the number of initial solitons from each FDM species
and $N_{s,1}=N_{s,2}=10$ in this case.

Figure 3 (upper row) shows radial profiles of the light and heavy field in
different simulations with $m_{2}=2m_{1}=1.6\times 10^{-22}\,\mathrm{eV}$ and
different mass fractions $\beta_{2}$. We observe distinctive results in each
simulation, which are discussed in order as follows. When the heavy field is
dominant, $\beta_{2}=0.7$, the two-field solitons fit the numerical simulation
well for both component fields. The one-field profile of the heavy field is
still aligned with its two-field profile, but a similar alignment is not seen
in the case of the light field. These findings completely agree with the
previous results in Scenario A. When the two fields are comparable in
abundance, $\beta_{2}=0.5$, the halo profile of each field still respects the
corresponding two-field ground state. However, both one-field profiles deviate
from the two-field ones. When the light field is dominant, $\beta_{2}=0.3$,
the inner core of this field is consistent with both one-field and two-field
solitons. In contrast, the inner profile of the heavy-field halo fails to fit
either of these solitons, which implies that the second soliton cannot form in
this case. Our results independently confirm that there exists a limit of a
mass fraction where the formation of the heavy-field soliton is inhibited, as
pointed out by Huang et al. (2023).

We estimate this threshold by running additional simulations which are set up
identically as above with $\beta_{2}$ varying continuously from $0.36$ to
$0.54$. We then calculate the goodness-of-fit of the heavy-field soliton which
describes how well it fits the corresponding simulation data, defined as

$\displaystyle\chi^{2}=\sum_{i=1}^{N_{\rm fit}}\dfrac{\log\rho_{{\rm
sim},2}(r_{i};\beta_{2})-\log\rho_{{\rm
sol},2}(r_{i};\beta_{2})}{\log\rho_{{\rm sol},2}(r_{i};\beta_{2})},$ (24)

where $\rho_{{\rm sol},2}$ and $\rho_{{\rm sim},2}$ are the two-field soliton
and the simulation profile of the heavy field, respectively; $N_{\rm fit}$ is
the number of simulation data points used counting from the halo center since
only these points are expected to fit the soliton. The inset of the upper
right panel in Fig. 5 shows $-\chi^{2}$ with respect to $\beta_{2}$ for
$N_{\rm fit}=3$ (orange) and $N_{\rm fit}=4$ (blue). A sudden drop from
$\chi^{2}\simeq 0$ implies no soliton found in the simulation data. It is,
hence, straightforward to infer the mass fraction threshold to be
approximately $\beta\simeq 0.44$ for Scenario B.

On the other hand, we must highlight two important properties of two-field FDM
haloes that are found in this set of simulations. Firstly, the one-field
profiles fail to describe two-field solitons, especially for the light field.
These discrepancies are particularly prominent when the heavy field becomes
more dominant ($\beta_{2}=0.5$ compared to $\beta_{2}=0.7$ in Scenario A) or
when the heavy-field mass is closer to the light-field mass ($m_{2}=2m_{1}$ in
Scenario B compared to $m_{2}=5m_{1}$ in Scenario A). As the heavy-field
soliton is more extended and massive in these two limits, its gravitational
pull on the light field is also larger and vice versa. The case
$\beta_{2}=0.7$ provides an example when the light soliton is so compressed
that its core is no longer flat while the case $\beta_{2}=0.5$ has the
solitons of both fields being squeezed under the gravitational potential of
each other. Thus, fitting solely one-field solitons to the numerical
simulation data could lead to the non-detection of solitons, which is a
misidentification of the real two-field ground state. Secondly, the threshold
of a mass fraction where the heavy-field core becomes disrupted is determined
by several characteristics of the two-field FDM system. The nested soliton can
form with a mass fraction of $\beta_{2}=0.25$ as seen in Scenario A, while the
same phenomenon does not happen in the case $\beta_{2}=0.3$ of Scenario B with
a higher mass fraction.

Figure 4: Snapshots of projected (comoving) density in four simulations at
$z=15.1,6.4,4.8,3.5$ to illustrate structure formation of different DM models
throughout cosmological history. These values of redshift are chosen to fully
capture the onset of soliton formation in each FDM simulation. The simulation
volume has a comoving side length of $L_{\rm box}=1.4\,h^{-1}\mathrm{Mpc}$.
The white box at $z=4.8$ in the bottom row marks the first structure in a two-
field FDM cosmology with $m_{2}=2m_{1}=2\times 10^{-22}\,\mathrm{eV}$ and
$\beta_{2}=0.5$. The close-up view of this region is shown in Fig. 5.

Figure 5: The upper panels show close-up slices through the center of the
first two-field halo formed at $z=4.8$ as seen in Fig. 4 (all quantities in
comoving coordinates). The lower panel shows the radial profile of this halo
in comoving coordinates, including a two-field soliton at the center (red) and
an outer NFW-like profile. The profiles of the light field (blue) and the
heavy field (green) are also displayed for comparison. Their solitonic cores
are just marginally resolved due to the limited resolution of our simulation.

However, comparing soliton formation in Scenario A and Scenario B is not
completely fair as there are other varying factors, in addition to $\beta_{2}$
and $m_{2}$, in these two scenarios such as the total mass $M$ (or
equivalently the mean density $\bar{\rho}$) and the initial conditions (core
radius, number of colliding solitons). Therefore, we set up Scenario C similar
to Scenario B but with the heavy-field mass varying,
$m_{2}=\\{4.0,4.8,5.6\\}\times 10^{-22}\,\mathrm{eV}$, and the density
fraction fixed to $\beta_{2}=0.3$ in three simulations. The core radius of
initial solitons also changes according to Eq. (23).

The final profiles are shown in the lower row of Fig. 3, where the heavy-field
soliton can only form with $m_{2}=7m_{1}=5.6\times 10^{-22}\,\mathrm{eV}$, but
not with $m_{2}=5m_{1}=4\times 10^{-22}\,\mathrm{eV}$ and
$m_{2}=6m_{1}=4.8\times 10^{-22}\,\mathrm{eV}$, for a fixed $\beta_{2}=0.3$.
The only difference in the last case with $m_{2}=7m_{1}$ is that its initial
state at $t=0$ consists of much denser solitons (see Eq. 15 and Tab. 2) which
enables the system to reach gravitational equilibrium faster. Once the stable
soliton of the heavy field has formed, it is no longer susceptible to tidal
disruption from density fluctuations of the light field. Thus, it is likely
that the formation of a nested soliton in the two-field FDM model is viable if
the heavy-field soliton formation has completed relatively earlier than the
light-field one. This condition can be realized even if the heavy-field
density is sub-dominant as long as the two FDM masses differ by orders of
magnitude, $m_{2}\gg m_{1}$.

##  5 Two-field FDM in cosmology

We perform DM-only cosmological simulations of the two-field FDM model with
self-consistent initial conditions.

###  5.1 Simulation setup

Four simulations are conducted in a comoving volume that has a length of
$L_{\rm box}=1.4\,h^{-1}\mathrm{Mpc}$ with a spatial resolution of $1024^{3}$
cells. The four simulations are set up as follows: (i) CDM via $N$-body
dynamics described by Vlasov-Poisson equations, (ii) one-field FDM with
$m_{a}=m_{1}=10^{-22}~{}\mathrm{eV}$ for the light field only (i. e.
$\beta_{2}=0$), (iii) one-field FDM with $m_{a}=m_{2}=2\times
10^{-22}~{}\mathrm{eV}$ for the heavy field only (i. e. $\beta_{2}=1$), (iv)
two-field FDM with $m_{1}$ and $m_{2}$ as given above and $\beta_{2}=0.5$. We
note that the mass fraction of two fields in this case is defined as the ratio
of cosmological density parameters, i. e.,

$\displaystyle\beta_{2}=\Omega_{\mathrm{a},2}/\Omega_{\mathrm{m}}=\Omega_{\mathrm{a},2}/(\Omega_{\mathrm{a},1}+\Omega_{\mathrm{a},2}),$
(25)

where $\Omega_{\mathrm{a},i}$ denotes the present-day density fraction of each
FDM field.

Initial conditions are generated from the same random seed for comparison at
$z=127$ using N-GenIC from Springel (2015) with a two-field matter power
spectrum calculated by the publicly available Boltzmann code aHCAMB from Luu
(2023)111The code aHCAMB, based on CAMB (Lewis et al., 2000) and axionCAMB
(Grin et al., 2022), was originally written for the axi-Higgs model. Thus, we
modify it to include additional perturbative equations of two axion fields
indirectly coupled via the metric tensor.. The heavy-field mass is assumed to
be $m_{2}=2\times 10^{-22}\,\mathrm{eV}$, two times heavier than the light-
field mass of $m_{1}=10^{-22}\,\mathrm{eV}$. Cosmological parameters are
chosen as constrained by Planck (Aghanim et al., 2020), i. e.,
$\Omega_{\mathrm{m}}=0.314,\Omega_{\Lambda}=0.686,h=0.6732,n_{s}=0.966$. One
exception is the primordial spectrum amplitude, which is boosted to
$A_{s}=7.5\times 10^{9}$ to enhance small-scale fluctuations. Thus, the
simulation box is more representative of an overdense region than an average
cosmological volume in the two-field FDM universe. The physical system in
these simulations evolved until redshift $z\simeq 3.5$, after which halo cores
become too small in the comoving grid to resolve.

###  5.2 Results and Discussion

Figure 4 provides a glimpse of the evolution of DM clustering in cosmological
simulations of four different DM models, as set up in Sec. 5.1. Each panel
shows the projected density distribution of DM structures which consists of
various DM haloes connected by cosmic filaments. As expected from previous FDM
simulations (Mocz et al., 2019; May & Springel, 2023), small-scale
fluctuations are smoothed out in all FDM models, compared to the scale-
independent clustering in the CDM model. The number of FDM haloes is, thus,
strongly suppressed due to a sharp cut-off below the Jeans scale in the
initial matter power spectra and the “quantum pressure” from FDM dynamics that
counteracts gravity. In terms of one-field simulations, the heavy-field
structures are more abundant and form earlier than the light-field ones
throughout cosmic history. The first structure of the heavy field emerges at
$z\simeq 6.4$, while the first light-field halo is not seen until $z\simeq
3.5$. We also identify four heavy-field haloes in total compared to only one
light-field halo at the final redshift.

On the other hand, the two-field FDM cosmology (with $\beta_{2}=0.5$) looks
like an intermediate state of the one-field models. Two haloes have formed by
the end of our simulation with the first one becoming visible at $z\simeq
4.8$. Zooming into this two-field halo in Fig. 5 (upper panel) reveals the
soliton core and wave interference patterns on scales of the de Broglie
wavelength of both fields. Since the total field is a superposition of two-
component fields, its core and granular structures are slightly denser and
fuzzier. We found that the radial profile of the central soliton exactly
follows the two-field ground state (lower panel), as in the non-cosmological
simulations, even though the solitons are just barely resolved. The nested
structure is not apparent in this case, owing to the low mass ratio of two
fields, i. e., $m_{2}=2m_{1}$.

More importantly, the haloes in our two-field simulation consist of solitons
of both light and heavy fields. These findings show that the light-field
soliton formation in the two-field simulation is greatly enhanced, as no
soliton appeared until the last redshift in the only-light-field simulation.
The light-field density likely started clustering faster thanks to deep
gravitational potentials induced by overdense regions of heavy-field haloes
which had formed beforehand in the two-field cosmology. Therefore, in the
limit of high mass ratios, $m_{2}\gg m_{1}$, we speculate that heavy-field
soliton formation may not be susceptible to the oscillating potential of the
light-field soliton in a realistic cosmological context, even with low-density
fractions $\beta_{2}$ between two fields, as its clustering commences
significantly earlier in the cosmic history. On the other hand, late-time DM
haloes may predominantly incorporate solitons of both fields in this context,
which poses a new challenge to explain the diverse core size of dwarf galaxies
by multiple FDM species (Luu et al., 2020; Pozo et al., 2023), where soliton
formation of each FDM species should be mostly independent.

##  6 Concluding remarks

As evidence has begun to build up against the one-field FDM model with mass
$m_{a}\sim 10^{-22}\,\mathrm{eV}$, the two-field FDM model with different
masses $m_{1}$ and $m_{2}$ offers an interesting and plausible alternative for
the FDM paradigm. In this study, we have conducted state-of-the-art
simulations to study structure formation in the two-field FDM model. With
simulations of soliton collision, we discovered that the virialized haloes
feature a nested soliton and an NFW-like outer profile. The nested soliton is
a stable configuration of two concentric cores of different sizes, which can
be solved as the ground state of two-field time-independent
Schrödinger–Poisson equations.

Unlike the one-field soliton profile, component solitons in the two-field
model are more compact due to mutual gravitational attraction. Depending on
the density ratio of the component fields, soliton formation may not occur in
case the heavy field is excessively dominated by the light field. Deep
gravitational potentials induced by soliton oscillations may be responsible
for this phenomenon (Huang et al., 2023). However, this threshold of soliton
instabilities is also sensitive to initial conditions and mass differences
between the two fields. A general pattern is that as the second field becomes
heavier, the nested soliton can form at a lower density fraction.

We also had a preliminary look into the large-scale structure of the two-field
FDM model. Despite the limited resolution, two-field soliton profiles are
found in our cosmological simulation. Since there are three viable profiles in
two-field FDM cosmology: two-field nested soliton, one-field light soliton,
and one-field heavy soliton, one might expect to find these different types of
haloes in a larger volume, which would explain the separate hierarchy of core
sizes in spheroidal dwarf galaxies and ultrafaint dwarf galaxies (Luu et al.,
2020). Unfortunately, every halo in our small-box simulation only contains a
two-field soliton, so we leave this interesting problem to a future study.

The two-field model can potentially alleviate several tensions between FDM
predictions and observational data. Hayashi et al. (2021) argued that
ultrafaint dwarf galaxy cores as FDM solitons imply $m_{a}\gtrsim
10^{-21}\,\mathrm{eV}$, incompatible with the $10^{-22}\,\mathrm{eV}$ mass for
dwarf spheroidal galaxies, which can be precisely explained by the heavy
field. Dalal & Kravtsov (2022) suggested that orbiting stars are perturbed as
they interact with the gravitational potential of density granules in FDM
haloes, which is particularly prominent for lighter FDM fields, so the mass
range $m_{a}\lesssim 3\times 10^{-19}\,\mathrm{eV}$ is strongly excluded. As
indicated by Gosenca et al. (2023), the presence of one (or more) heavier
field(s) helps to smooth out the density fluctuations in FDM haloes, hence
relaxing the bound. Lastly, FDM suppresses the linear matter power spectrum at
small scales, which alters the structure of the intergalactic medium (IGM) on
the scales probed by the Lyman-$\alpha$ forest of early-time quasars. By
relating the one-dimensional flux spectrum with the FDM power spectrum,
Armengaud et al. (2017) and Iršič et al. (2017) constrained
$m_{a}\lesssim\times 10^{-21}\,\mathrm{eV}$, although there exist several
uncertainties of IGM gas physics that should be better understood (Hui et al.,
2017). We expect a sufficiently large density fraction of the heavy field will
resolve the Lyman-$\alpha$ constraint.

Nevertheless, more detailed studies on the two-field FDM model in theory,
phenomenology, and simulation are needed to settle these issues. On the other
hand, since the nested soliton is exclusive to two (or multiple) species of
FDM, potential observations of this profile in future galaxy surveys will be a
smoking gun signature of the multiple-axion scenario.

## Acknowledgements

The simulations in this paper were run on the Engaging cluster supported by
MIT. H.N. thanks Mihir Kulkarni for help with setting up the initial
conditions of two-field cosmological simulations, Xuejian Shen, and Leo W. H.
Fung for useful discussions to improve the paper. P.M. acknowledges this work
was in part performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344,
Lawrence Livermore National Security, LLC. S.M. acknowledges support from the
National Science Foundation under Grant No. 2108931.

## References

  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astron. Astrophys., 641, A6
  * Alam et al. (2021) Alam S., et al., 2021, Phys. Rev. D, 103, 083533
  * Amin et al. (2022) Amin M. A., Jain M., Karur R., Mocz P., 2022, JCAP, 08, 014
  * Armengaud et al. (2017) Armengaud E., Palanque-Delabrouille N., Yèche C., Marsh D. J. E., Baur J., 2017, Mon. Not. Roy. Astron. Soc., 471, 4606
  * Arvanitaki et al. (2010) Arvanitaki A., Dimopoulos S., Dubovsky S., Kaloper N., March-Russell J., 2010, Phys. Rev. D, 81, 123530
  * Boylan-Kolchin et al. (2011) Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40
  * Dalal & Kravtsov (2022) Dalal N., Kravtsov A., 2022, Phys. Rev. D, 106, 063517
  * Di Cintio et al. (2014) Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, Mon. Not. Roy. Astron. Soc., 437, 415
  * Eby et al. (2020) Eby J., Leembruggen M., Street L., Suranyi P., Wijewardhana L. C. R., 2020, JCAP, 10, 020
  * Ferreira (2021) Ferreira E. G. M., 2021, Astron. Astrophys. Rev., 29, 7
  * Garrison-Kimmel et al. (2014) Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Kirby E. N., 2014, Mon. Not. Roy. Astron. Soc., 444, 222
  * Glennon et al. (2023) Glennon N., Musoke N., Prescod-Weinstein C., 2023, Phys. Rev. D, 107, 063520
  * Gosenca et al. (2023) Gosenca M., Eberhardt A., Wang Y., Eggemeier B., Kendall E., Zagorac J. L., Easther R., 2023, Phys. Rev. D, 107, 083014
  * Governato et al. (2012) Governato F., et al., 2012, MNRAS, 422, 1231
  * Grin et al. (2022) Grin D., Marsh D. J. E., Hlozek R., 2022, axionCAMB: Modification of the CAMB Boltzmann code, Astrophysics Source Code Library, record ascl:2203.026 (ascl:2203.026)
  * Guo et al. (2021) Guo H.-K., Sinha K., Sun C., Swaim J., Vagie D., 2021, JCAP, 10, 028
  * Hayashi et al. (2021) Hayashi K., Ferreira E. G. M., Chan H. Y. J., 2021, Astrophys. J. Lett., 912, L3
  * Hu et al. (2000) Hu W., Barkana R., Gruzinov A., 2000, Phys. Rev. Lett., 85, 1158
  * Huang et al. (2023) Huang H., Schive H.-Y., Chiueh T., 2023, Mon. Not. Roy. Astron. Soc., 522, 515
  * Hui (2021) Hui L., 2021, Ann. Rev. Astron. Astrophys., 59, 247
  * Hui et al. (2017) Hui L., Ostriker J. P., Tremaine S., Witten E., 2017, Phys. Rev. D, 95, 043541
  * Iršič et al. (2017) Iršič V., Viel M., Haehnelt M. G., Bolton J. S., Becker G. D., 2017, Phys. Rev. Lett., 119, 031302
  * Jain et al. (2023) Jain M., Amin M. A., Thomas J., Wanichwecharungruang W., 2023, arXiv e-prints
  * Ji & Sin (1994) Ji S. U., Sin S. J., 1994, Phys. Rev. D, 50, 3655
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Lovell et al. (2012) Lovell M. R., et al., 2012, MNRAS, 420, 2318
  * Luu (2023) Luu H. N., 2023, Phys. Rev. D, 107, 023513
  * Luu et al. (2020) Luu H. N., Tye S. H. H., Broadhurst T., 2020, Phys. Dark Univ., 30, 100636
  * Marsh (2016) Marsh D. J. E., 2016, Phys. Rept., 643, 1
  * May & Springel (2021) May S., Springel V., 2021, Mon. Not. Roy. Astron. Soc., 506, 2603
  * May & Springel (2023) May S., Springel V., 2023, Mon. Not. Roy. Astron. Soc., 524, 4256
  * Mocz et al. (2017) Mocz P., Vogelsberger M., Robles V. H., Zavala J., Boylan-Kolchin M., Fialkov A., Hernquist L., 2017, Mon. Not. Roy. Astron. Soc., 471, 4559
  * Mocz et al. (2018) Mocz P., Lancaster L., Fialkov A., Becerra F., Chavanis P.-H., 2018, Phys. Rev. D, 97, 083519
  * Mocz et al. (2019) Mocz P., et al., 2019, Phys. Rev. Lett., 123, 141301
  * Mocz et al. (2023) Mocz P., et al., 2023, Mon. Not. Roy. Astron. Soc., 521, 2608
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, Astrophys. J., 462, 563
  * Oman et al. (2015) Oman K. A., et al., 2015, Mon. Not. Roy. Astron. Soc., 452, 3650
  * Pakmor et al. (2016) Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal K., Zhu C., 2016, Mon. Not. Roy. Astron. Soc., 455, 1134
  * Powell et al. (2023) Powell D. M., Vegetti S., McKean J. P., White S. D. M., Ferreira E. G. M., May S., Spingola C., 2023, doi:10.1093/mnrasl/slad074
  * Pozo et al. (2023) Pozo A., Broadhurst T., Smoot G. F., Chiueh T., 2023, arXiv e-prints
  * Roszkowski et al. (2018) Roszkowski L., Sessolo E. M., Trojanowski S., 2018, Rept. Prog. Phys., 81, 066201
  * Schive et al. (2014a) Schive H.-Y., Chiueh T., Broadhurst T., 2014a, Nature Phys., 10, 496
  * Schive et al. (2014b) Schive H.-Y., Liao M.-H., Woo T.-P., Wong S.-K., Chiueh T., Broadhurst T., Hwang W. Y. P., 2014b, Phys. Rev. Lett., 113, 261302
  * Schwabe et al. (2016) Schwabe B., Niemeyer J. C., Engels J. F., 2016, Phys. Rev. D, 94, 043513
  * Schwabe et al. (2020) Schwabe B., Gosenca M., Behrens C., Niemeyer J. C., Easther R., 2020, Phys. Rev. D, 102, 083518
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel (2015) Springel V., 2015, N-GenIC: Cosmological structure initial conditions, Astrophysics Source Code Library, record ascl:1502.003 (ascl:1502.003)
  * Svrcek & Witten (2006) Svrcek P., Witten E., 2006, JHEP, 06, 051
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Rev. Phys., 2, 42
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, Astrophys. J. Suppl., 248, 32
  * Zwicky (1937) Zwicky F., 1937, ApJ, 86, 217
  * de Blok (2010) de Blok W. J. G., 2010, Advances in Astronomy, 2010, 789293


chatGPT: |  You’re welcome. If you have any other questions related to machine learning or artificial intelligence, don’t hesitate to ask!
Example No. |  5
User: |  Predict total arrival delay for each destination airports whose flights will start next Sunday
chatGPT: |  Great! To help me understand your problem better, can we first narrow down the scope by fixing the aggregation operation, filter attribute and filter operation? It will be easier to formulate the problem.
User: |  can you fill the slots based on the prompt provided
chatGPT: |  Sure! So, based on the prompt, the ‘aggregation operation’ is ‘total’, the ‘filter attribute’ is ‘day of week’, and the ‘filter operation’ is ‘equal’. Since you want to predict total arrival delay for each destination airport whose flights will start next Sunday, then the ‘target attribute’ should be ‘arrival delay’ and the ‘filter value’ should be ‘Sunday’. Is that correct?
User: |  yes.
chatGPT: |  Great! Now that we have all the necessary slots filled, we can formulate the problem as follows: Problem Statement: Forecast the total arrival delay of flights for each destination airport, whose flights are scheduled for next Sunday. Aggregation Operation: Total Target Attribute: Arrival Delay Filter Attribute: Day of week Filter Value: Sunday Does this make sense?
User: |  i do not need filter value, instead, i need filter operation.
chatGPT: |  Sorry for the confusion. As per the prompt, we can set ‘filter operation’ to ‘equal’ since we want to predict arrival delay at each destination airport only for flights scheduled on Sunday. So, the complete problem statement is: Problem Statement: Forecast the total arrival delay of flights for each destination airport, whose flights are scheduled for next Sunday. Aggregation Operation: Total Target Attribute: Arrival Delay Filter Attribute: Day of Week Filter Operation: Equal
User: |  correct.