Astronomy & Astrophysics manuscript no.

aa54938-25 ©ESO 2025

August 8, 2025

DAWN
I. Simulating the formation and early evolution of stellar clusters with Phantom
N-Body
Y. Bernard1 , E. Moraux1 , D. J. Price2 , F. Motte1 , F. Louvet1 , and I. Joncour1

1
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
2
School of Physics and Astronomy, Monash University, Vic. 3800, Australia
arXiv:2508.05296v1 [astro-ph.GA] 7 Aug 2025

Received April 1, 2025; accepted August 4, 2025

ABSTRACT

Context. The star formation process produces hierarchical, clustered stellar distributions through gravoturbulent fragmentation of
molecular clouds. Simulating stellar dynamics in such an environment is numerically challenging due to the strong coupling between
young stars and their surrounding gas, and the large range of length and time scales.
Aims. This paper is the first of a suite aimed at investigating the complex early stellar dynamics in star-forming regions, from the
initial collapse of the molecular cloud to the phases of complete gas removal. We present a new simulation framework. This advanced
framework is the key to generating a larger set of simulations, enabling statistical analysis, which is mandatory to address the stochastic
nature of dynamical interactions.
Methods. Methods originating from the stellar dynamics community, including regularisation and slowdown methods (SDAR), have
been added to the hydrodynamical code Phantom to produce simulations of embedded cluster early dynamics. This is completed by a
novel prescription of star formation to initialise stars with a low numerical cost, but in a way that is consistent with the gas distribution
during the cloud collapse. Finally, a prescription for H ii region expansion has been added to model the gas removal.
Results. We have run testcase simulations following the dynamical evolution of stellar clusters from the cloud collapse to a few Myr.
Our new numerical methods fulfil their function by speeding up the calculation. The N-body dynamics with our novel implementation
never appear as a bottleneck that stalls the simulation before its completion. Our new star formation prescription avoids the need
to sample individual star formations within the simulated molecular clouds with a large number of SPH particles. Overall, these
new developments allow accurate hybrid simulations in minimal calculation time. Our first simulations show that massive stars
largely impact the star formation process and shape the dynamics of the resulting cluster. Depending on the position of these massive
stars and the strength of their H ii regions, they can prematurely dismantle part of the cloud or trigger a second event of cloud
collapse, preferentially forming low-mass stars. This leads to different stellar distributions for numerical simulations with similar
initial conditions and confirms the need for statistical studies. Quantitatively, and despite the implementation of feedback effects, the
final star formation efficiencies are too high compared with those measured in molecular clouds of the Milky Way. This is probably
due to the lack of feedback mechanisms other than H ii regions, in particular jets, non-ionising radiation, or Galactic shear.
Conclusions. Our new Phantom N-Body framework, coupled with the novel prescription of star formation, enables efficient simulation
of the formation and evolution of star clusters. It enables the statistical analysis needed to establish a solid theoretical framework for
the dynamical evolution of embedded star clusters, continuing the work done in the stellar dynamics community.
Key words. Star formation — Simulations — Embedded clusters — Stellar dynamics

1. Introduction Observations have shown that this clustered state is transient

for 90% of stars (Lada & Lada 2003), and that only 10% of new-
born stars remain in open clusters after 10 Myr. This dynami-
Star formation occurs within turbulent, collapsing molecular cal evolution is fast and complex, and a better understanding of
clouds in the Galaxy. The gravitational collapse of these tur- this evolution is needed to constrain environmental properties re-
bulent clouds produces interconnected filaments and hubs that quired in other fields. For example, the statistics of fly-by events
amplify density anisotropies. Observations with the Spitzer tele- during the early evolution of star-forming regions could help to
scope have revealed that the star formation process takes place constrain their effect on planet formation within protoplanetary
preferentially inside these dense structures in star-forming re- disks (Cuello et al. 2023).
gions (Evans et al. 2009; Gutermuth et al. 2009; André et al. Pure N-body simulations (Aarseth & Hills 1972; Goodwin &
2014). Denser regions with similar levels of turbulent and mag- Whitworth 2004; Bastian & Goodwin 2006; Parker et al. 2014;
netic support are more sensitive to fragmentation, leading to a Wright et al. 2014) have tried to explain the rapid dissolution of
dominant ‘clustered’ star formation channel. Consequently, the natal star clusters with stellar dynamics. These simulations mod-
distribution of newborn stars is inherited from the distribution of elled young stellar clusters using an initially fractal spatial distri-
filaments and hubs/ridges (Schneider et al. 2010; Peretto et al. bution, consistent with observational constraints (Cartwright &
2013; Motte et al. 2018), and appears to be structured, with mul- Whitworth 2004). The results indicate that dynamical instabili-
tiple clusters embedded in their birth nebula/clump. ties, such as close encounters or unstable triples, can lead to the
Article number, page 1 of 21
 A&A proofs: manuscript no. aa54938-25

complete dissolution of an initial cluster or to the conservation of dynamics community into their hydrodynamic codes to perform
a strongly bound cluster of stars at the end of the simulation. Fur- what is called a hybrid hydro/N-body simulation.
thermore, the statistical robustness of these findings is ensured These simulations aimed to bridge the gap between previous
by the averaging of hundreds of simulations to draw conclusions. work by the stellar dynamics community and star cluster for-
However, pure N-body simulations ignore one key component mation simulations, including stellar feedback. They provided
of this dynamical evolution: the strong links and interactions be- a framework to explain the dynamical evolution of young star
tween newborn stars and their gaseous environment. Previous clusters influenced by their environment. However, they did not
N-body studies have attempted to account for the gas potential follow the statistical approach of previous pure N-body clus-
using simple radial analytical formulae with a decay in time to ter simulations, which is needed to identify typical evolutionary
model the gas expelled from embedded clusters (Kroupa et al. tracks followed during the stochastic dynamics of stars in a clus-
2001). More recent simulations have focused on the early phases ter.
of star formation, revealing the complex and time-variable na- In this paper, we show how methods developed for acceler-
ture of this coupling that influences both the dynamics and the ating N-body simulations and efficient star formation and stellar
formation of stars.(Verliat et al. 2022; Fujii et al. 2021; Grudić feedback prescriptions can be incorporated directly into a 3D
et al. 2021; González-Samaniego & Vazquez-Semadeni 2020). hydrodynamics code, adding missing physics to pure N-body
Hence, only direct methods can model star formation and stellar simulations. Such a framework allows us to initialise young
feedback simultaneously to follow the stellar dynamics during embedded clusters directly from star formation in a molecular
the early evolution of clusters in star-forming regions. cloud collapse, leading to self-consistent cluster initial condi-
On the other hand, numerical studies focusing on star for- tions linked with the natal gaseous environment. Their evolu-
mation have naturally characterised this process using a direct tion can therefore be followed through the direct influence of the
hydrodynamic model. Larson (1969) and Boss & Bodenheimer gas, itself influenced by gravity and stellar feedback. We aim to
(1979) focused on a single, isolated clump of gas gravitation- follow this evolution up to the emergence of clusters from their
ally collapsing to form single or binary stars. Such works tried to natal cloud. Our new framework has been designed to be fast
self-consistently follow the collapse from the initial dense clump enough to enable the statistical approach necessary to study the
up to the onset of nuclear fusion (Bate 1998; Bate et al. 2014; dynamical evolution of embedded clusters.
Wurster et al. 2018). Recent studies (e.g., Ahmad et al. 2024; The paper is organised as follows. First N-body algorithms
Mauxion et al. 2024) continue to refine models of protostar for- newly implemented in the code Phantom (Price et al. 2018) are
mation, incorporating more realistic physical processes to pre- presented in Sect. 2. Section 3 describes the new star forma-
dict phenomena such as the formation and evolution of proto- tion prescription (Sect. 3.1), the implementation of H ii region
planetary disks that are essential to planet formation. feedback (Sect. 3.2), and the fiducial model studied in this work
However, the star formation community rapidly understood (Sect. 3.3). Simulation results of this model are described and
that the environment is key to explaining stellar properties. discussed in Sect. 4. Section 5 discusses numerical performance
Hence, it became critical to perform simulations at the molecular and accuracy. Finally, Section 6 gives our conclusions and sum-
cloud scale while modelling the formation of individual proto- marises future prospects. Paper II will then be dedicated to a
stars. Bate et al. (1995) developed a method called sink particles deeper investigation of our new star formation prescription, fo-
that hides dense parts of a collapsing clump inside a massive cusing on the mass distribution and multiplicity statistics.
particle capable of accreting infalling material without computa-
tionally expensive hydrodynamic calculations. This method un-
locked the ability to follow the evolution of star-forming regions 2. Phantom
beyond the point of initial collapse. As a consequence, sink par-
ticles are now used in almost every star formation simulation. and N-body algorithms Simulations presented in this work were
Early studies (Bate et al. 2002; Bate & Bonnell 2005) simu- performed with the 3D Smoothed Particle Magnetohydrodynam-
lated small, isolated molecular clouds with a diameter of a frac- ics (SPMHD) code Phantom (Price et al. 2018). Smoothed par-
tion of a parsec and an initial mass of tens of solar masses. Later ticle hydrodynamics (SPH; Lucy 1977; Gingold & Monaghan
developments enabled the modelling of radiative feedback from 1977) solvers use a Lagrangian approach to solve fluid equa-
protostars on the surrounding gas, which was found to be es- tions by discretising the medium into particles of fixed mass that
sential to produce a more convincing stellar mass distribution follow the fluid motion. Fluid forces are derived by calculating
(e.g., Bate 2009; Hennebelle et al. 2020). Other studies tried to quantities on these particles with neighbour weighted sums using
identify links between molecular cloud initial conditions, such as smoothing kernels (Price 2012). We chose Phantom because of
the temperature, magnetic field (Myers et al. 2014; Lee & Hen- its versatility. It is capable of modelling a variety of astrophysical
nebelle 2019) or initial turbulent state (Joos et al. 2013; Liptai systems with its numerous physics modules, like an MHD solver
et al. 2017), and stellar properties like the mass or the spatial dis- (ideal and non-ideal), self-gravity, dust, chemical networks, sink
tribution. Most Recent studies have attempted to model physical particles and more. However, this code was not ready yet to per-
processes in star-forming regions more realistically. Specifically, form embedded cluster simulations at the beginning of this work.
stellar feedback is thought to play a major role in both star for- While a few studies on star formation have already been carried
mation and the global evolution of molecular clouds (Krumholz out (Liptai et al. 2017; Wurster & Bonnell 2023), these studies
et al. 2014). Verliat et al. (2022) found that protostellar jets can stopped the evolution of their star-forming region at around the
regulate accretion around dense objects by injecting turbulent initial free fall time of the cloud τff . This end time is limited by
support inside the surrounding gas. State-of-the-art simulations the large computational cost of such models. One of the growing
attempt to model the evolution of giant molecular clouds through constraints is the gravitational dynamics between newly formed
the entire star formation process, from the initial collapse to the stars when no softening is applied, which can then form in close
gas removal mainly driven by H ii region expansion produced by binaries or multiple systems. Stars in such systems are subject to
massive stars (Guszejnov et al. 2022). For instance,Wall et al. strong interactions that result in small time-step constraints that
(2019) and Fujii et al. (2022) reused methods from the stellar can dramatically slow down a simulation.
Article number, page 2 of 21
 Y. Bernard et al.: DAWN

In Phantom, gravitational forces between point masses (as kick operation for a particle i is then
stars) are computed using a direct method. This direct algorithm
δt
is coupled to a Leapfrog symplectic integrator to follow the dy- ai∗ = ai + gi . (4)
namics of stars. This integrator is only 2nd order, which needs 48
short timesteps to reach sufficient accuracy. This combination of This formula is not trivial to compute. First, all pairwise acceler-
methods is inefficient and can be inaccurate during the early evo- ations are required to compute each gradient term, which means
lution of embedded clusters. As the number of steps to integrate that they need to be evaluated beforehand. During the middle
a cluster up to a few million years is prohibitive, it also implies kick phase, there is therefore one acceleration and one gradient
an accumulation of truncation errors that can drastically change acceleration evaluation consecutively to compute the final ac-
the evolution of close binaries or multiples, which are key in the celeration to kick the velocity of each body in the system. This
overall cluster evolution. A hard binary disruption due to inac- method turned out to be expensive to implement in the current
curacy in the integration scheme can lead to the artificial fast ex- Phantom framework. Indeed, gradient accelerations in a simu-
pansion of the host cluster. These numerical effects motivate new lation where softening is applied between stars and gas need to
implementations written in Phantom that reuse methods devel- take this softening into account, which is possible but numer-
oped in collisional N-body codes. The SDAR method proposed ically complex. Instead, Omelyan (2006) proposed to compute
by Wang et al. (2020b) and the forward symplectic integrator this modified acceleration by extrapolating it from a slightly
(FSI) recently presented in Rantala et al. (2021) and Dehnen & shifted position of each body, which depends on the standard
Hernandez (2017). For the sake of completeness, we describe acceleration. Formally, it is written as
our implementation from these original papers in App. A.1.
δt2
ai∗ = a(ri∗ ); where ri∗ ≡ ri + a(ri ) (5)
24
2.1. Forward Symplectic integrator: FSI
and a is the function to compute the acceleration for each particle
Symplectic integrators preserve conservation properties. They in the system. One can see in this formula that there are still two
can conserve angular momentum to round-off error while con- acceleration evaluations. However, it reuses the same function
straining the energy error inside a small range, depending on the in both cases, with a shifted position in the latter one. Omelyan
truncation order and the timestep choice. They achieve that by (2006) and Chin & Chen (2005) showed that this new formula
rewriting the time evolution operator. This new formulation im- gives results close to the original. This method is far easier to im-
plies that the integrated Hamiltonian is slightly different from plement in an existing framework. Hence, we chose this method
the true one by what is called the error Hamiltonian. The art instead of the original one.
of building high-order symplectic integrators lies in minimising In terms of precision, FSI outperforms all other standard
this error term (Yoshida 1990). The Forward Symplectic Inte- fourth-order symplectic schemes like Forest-Ruth integrator
grator (or FSI) (Dehnen & Hernandez 2017; Rantala et al. 2021) (Forest & Ruth 1990) and even improved ones like Position
achieves this minimisation up to fourth order with the unique extended Forest Ruth like integrator (PEFRL)(Omelyan et al.
property of producing no backwards substeps inside its scheme, 2002) (See Appendix A.1). Chin & Chen (2005) argued that
which is essential if dissipative forces are in the scheme. More even with the same order of truncation error, the forward scheme
details on how to construct symplectic integrators are presented will always outperform standard schemes because others evalu-
in Appendix A.1. That aside, one can construct the new time ate forces far from the true trajectory of the object. By contrast,
evolution operator derived from this special method using forward schemes evaluate forces very close to the true trajec-
δt2 δt2
tory with the same time step resolution. In the end, standard
eδt H = eδt (T +U+ 48 G) ≈ e 6 δt U e 2 δt T e 3 δt (U+ 48 G) e 2 δt T e 6 δt U ,
1 1 2 1 1
(1) schemes need smaller time steps to achieve the same resolution.
FSI seems to be the best option for a new integration algorithm
where H is the Hamiltonian, T the kinetic energy, U the po- of point masses in Phantom. It has one of the best precisions for
tential and G an extra potential. All the coefficients in this time fourth-order schemes with limited supplementary force evalua-
evolution operator are strictly positive. The middle operation is tion in one global step. Its forward nature assures that nothing
modified with what is commonly called a gradient force eval- can go wrong if dissipation is present in the system, which is
uation produced by a slightly different potential U + δt48 G (cf.,
2
often the case in astrophysical hydrodynamical simulations.
Rantala et al. 2021). Knowing the definition of G and for a sys-
tem composed only of N point masses coupled by gravitational
interactions where aj is the acceleration sensed by each mass, 2.2. Slow down Algorithmic regularisation: SDAR
the gradient acceleration can be written as 2.2.1. Regularised integration
N
1 ∂G ∂ X The FSI presented above is an accurate integrator that already
gi = − = ∥a j ∥. (2) outperforms second-order leapfrog schemes. It is then possible
mi ∂ri ∂ri j to limit the number of steps per orbit needed to sample a good
trajectory of each body. Typically, this number of sample points
Expanding the partial derivative in this equation leads to the is set from 300 to 500 with the previous Leapfrog integrator.
common formula of the jerk computed in Hermite fourth-order With the new FSI, we can easily divide this number by a factor
integrator (Makino & Aarseth 1992) and gives of four while retaining better accuracy.
N X N However, updating the main integrator is not enough to re-
X Gm j h 2 i
solve the stellar dynamics within an embedded cluster in a rea-
gi = 2 ri j ai j − 3(ri j · ai j )ri j , (3)
i=1 i, j
ri j sonable computation time. Stellar densities are high in cluster
environments, and binaries and multiple systems represent a sig-
where ai j and ri j are the relative acceleration and position be- nificant fraction of the population (Offner et al. 2023). These sys-
tween particles i and j. The modified acceleration of the middle tems can produce strong computational constraints during close
Article number, page 3 of 21
 A&A proofs: manuscript no. aa54938-25

encounters or with short-period binaries. For example, a binary However, the idea is then to use a specialised algorithm to pro-
with 10 au separation and a total mass of 2 M⊙ has a period duce regular results. There are two main algorithms to achieve
around 20 yr. It could then produce a time step constraint around that. The LogH method presented in Rantala et al. (2023)
δt = 0.33 yr when the typical free fall time of a molecular cloud or Aarseth (2003) and the time-transformed leapfrog (TTL)
in this study is around 5.5 Myr, which gives a large dynamic (Mikkola & Aarseth 2002). We adopt the latter method in this
range. This simple example does not capture the stringent time work, which is described in detail in Appendix A.2. Briefly, it is
step constraints that could appear in a young cluster. Day orbit possible to construct a new symplectic integrator using the mod-
period binaries and high eccentric ones make direct computation ified system of equations. The change of variable produces a new
infeasible. integrated momentum-like variable W defined as
The stellar dynamics community has already tackled these
high dynamical scale problems to model the long-term evolu- N
dW X 1 dri ∂Ω
tion of stellar clusters containing a few thousand bodies to a few = · . (7)
tenths of a million. To meet this challenge, sophisticated meth- ds i=1
Ω dt ∂ri
ods have been developed to manage the two main issues at large
dynamic ranges, which are the gravitational singularity during This new variable will then have the same behaviour as ve-
pair interaction and the secular evolution of multiple systems. locities in a symplectic scheme. The Leapfrog method that
First, gravitational pair interaction is a function of the inverse lends its name to this method is not precise enough in such
square distance between two bodies. It produces a singularity at applications. Two solutions have been used to reach a satis-
the distance d = 0. This singularity can strongly constrain the factory accuracy. First, the original algorithm from Mikkola’s
simulation time step. During close encounters or periapsis ap- code and the new code BIFROST (Rantala et al. 2023) used
proaches in a highly eccentric binary, the velocity of the two the Gragg–Bulirsch–Stoer (GBS) extrapolation method (Gragg
bodies can dramatically increase due to the inverse square dis- 1965; Bulirsch & Stoer 1966) to ensure a user-chosen relative
tance proportionality of their interactions. Their true trajectory error. The second solution used in PeTar (Wang et al. 2020a),
will be missed if there are not enough sample points. Worse, another stellar dynamics code, is to integrate the system with a
the pair interaction can produce false runaway bodies. In colli- high-order symplectic scheme from Yoshida (1990). We chose
sionless models as cosmological ones like Gadget 4 (Springel the latter method for its simplicity. However, GBS extrapolation
et al. 2021), the interaction is softened in a certain radius around will be tested in further work.
each point mass to avoid these strong pair interactions. However, The essence of this method can be understood from the time
softening methods are prohibited in collisional N-body codes, as transformation (Eq. 6). The link between the new integration
the collisional dynamics govern the future evolution of a cluster. variable s and the real time t depends on the potential strength.
The immediate solution could be to increase the number of steps. This way, s can have constant large time steps as the integrated
However, this would also produce large truncation and round-off time t will scale depending on the potential strength. It means
errors due to the accumulation of steps in addition to being com- that when the gravitational forces are strong (close to the singu-
putationally expensive. larity), the time integration will be slow and algorithmically ad-
Various methods have been proposed in the literature for fifty justed to the potential strength. That way, singularity approaches
years now (Kustaanheimo et al. 1965; Mikkola & Aarseth 2002) are not an issue anymore even if it is still present in equations.
to manage these two issues. All of them are based on the same
idea: regularisation of gravitational interactions. The idea is to
hide the singularity. One manipulates the equations of motion by 2.2.2. Slowed down binaries
applying spatial and time transformations on the system of equa-
tions to make the 1/d2 term in the force evaluation vanish. The The last challenge to tackle is the secular evolution of perturbed
first regularisation method used in this field was proposed by binaries and multiple systems. Regularised integrators can solve
Kustaanheimo et al. (1965) (KS method). It transforms binary collisional trajectory issues. Nevertheless, they still need a few
gravitational motion into a simple harmonic oscillator, meaning sample points to integrate orbits. It means that very hard, eccen-
that the time step can be constant during a whole period of the bi- tric, and weakly perturbed binaries or stable triple systems can
nary, whatever its eccentricity. The details of this transformation bottleneck the computation, as the dynamic time range between
can be found in the original paper and in Aarseth (2003). It has a cluster and these systems can be extremely large. Because of
been successfully used in several collisional codes in the past, their hardness, these binaries can be considered unperturbed, and
such as the series of Aarseth’s N-Body codes (Aarseth 1999). thus their motion can be found analytically using Kepler solvers.
However, its complexity makes it hard to implement. In addi- It resolves only a part of the issue, as stable triples are com-
tion, KS regularisation of multiple systems requires specialised monly found in simulations and close encounters in dense stellar
methods, which add up to a rather complex method. More recent regions can greatly perturb binary systems. Mikkola & Aarseth
codes (Wang et al. 2020a; Rantala et al. 2023) preferred a dif- (1996) proposed an elegant solution with the cost of losing the
ferent approach to reach the same results, using what is called orbital phase of these hard binaries. The slowdown method pro-
an algorithmic regularisation (AR) (Mikkola & Tanikawa 1999; poses to slightly modify the original Hamiltonian of a weakly
Preto & Tremaine 1999). This work follows the same approach. perturbed binary according to
AR methods apply a time transformation to the system
1
ds = Ω(r)dt, (6) H= Hb + (H − Hb ), (8)
κ
where s is a new time integration variable that is linked with the
real time t by a time transformation function Ω(r). This function where κ is called the slow down factor and Hb is the unperturbed
is often taken as the gravitational potential of the system. binary Hamiltonian. This equation means that only the binary
By contrast to KS, it leaves the spatial coordinates un- dynamics are scaled down by a factor κ. Original velocities are
changed. The singularity is still present in a pair interaction. conserved by this scaling. Equations of motion derived from this
Article number, page 4 of 21
 Y. Bernard et al.: DAWN

Hamiltonian are given by the SDAR method is used to integrate the motion of stars but
not SPH gas particles. Thus, the new methods are applied to a
dvi Gmc (ri − rc ) ∂U(ri ) specific part of the system and not the entire one.
= − − , (9)
dt κ∥ri − rc ∥3 ∂ri The original Phantom integration scheme was derived from
dri vi the reference system propagator algorithms (RESPA) (Tucker-
= , (10) man et al. 1992) which split forces applied on particles in the
dt κ
system into two categories, the “slow” and “fast” forces with re-
where U is an arbitrary external potential and i and c index points spectively weak and strong time step constraints. This splitting
on one binary member and its companion, respectively. One can allows the construction of a two-level leapfrog scheme according
see that the forces due to gravitational interactions between the to
primary and the secondary are scaled down by κ. The original
orbital period of the binary is divided by κ. It is then straight- v = v + 12 ∆tsph asph (11)
forward to understand that the time constraint on the binary will 
 v = v + 12 ∆tfast afast
then be relaxed by the same factor.

 r = r + ∆tfast v



Equations 9 and 10 show that for a given binary system, its N substeps  (12)

 afast = afast (r)
orbital parameters are conserved by the scheme, despite pertur-

 v = v + 1 ∆t a

2 fast fast
bations that will be discussed shortly. It also means that even if
the apparent period is κ times longer than the original one, the asph = asph (r) (13)
true orbital parameters are conserved by the scheme. The only v=v+ 2 ∆tsph asph ,
1
(14)
change comes from the external perturbation.
Wang et al. (2020b) showed that applying this method on one where asph corresponds to SPH accelerations (pressure, viscos-
slowed-down period Psd = κP gives the same velocity modifica- ity, self-gravity...). ∆tsph is the time step associated with them and
tion as if the system was integrated during κP. It means that the used in the first leapfrog level. The drift phase of the first-level
slowdown method averages perturbations applied to the binary leapfrog is replaced by multiple substeps of another leapfrog in
on a slowed-down period Psd . Nonetheless, it is only applicable association with a smaller time step ∆tfast . This sub-timestep is
if the perturbation is weak enough. Thus, it means that we should controlled by faster accelerations contained in afast such as ex-
choose κ carefully to preserve the accuracy of the averaging. On ternal potential, sink-gas, sink-sink accelerations, etc...
intermediate perturbed binaries, this averaging property is rele- The methods we presented above relate to stars, which cor-
vant only for a few orbits. By contrast, external perturbations can respond to point masses in Phantom (sink particles without soft-
be almost negligible for a hard binary. In that case, the averaging ening). It means that only the second level needs to be updated,
property can be accurate enough on several orders of magnitude as it is where sink particles are updated during integration. The
of the original orbital period. first upgrade was to switch the coarse Leapfrog integrator for the
The generalisation to a few-body problem is straightforward more accurate FSI presented above. The substep integrator can
and detailed in Wang et al. (2020b). Only two changes are then be rewritten as
needed: i) the bulk motion of a whole binary in an N-body sys- 1
tem is not slowed down and ii) the perturbation on a binary of the v = v + ∆tfast afast (15)
6
system is now the summation of each direct gravitational force 1
produced by each perturber inside the system. As shown above, κ r = r + ∆tfast v (16)
quantifies the number of original periods where the perturbation 2
applied to the binary can be averaged. On a highly perturbed ∆tfast
2
ãfast = afast (r∗ ); r∗ = r + afast (r) (17)
binary, it is impossible to average the perturbation on multiple 24
periods, which means κ = 1. A hard binary can stay unperturbed 2
during a whole simulation. The averaging is not risky at all, and v = v + ∆tfast ãfast (18)
3
then κ can stay at a high value. However, one can see that this fac- 1
tor needs to scale with the strength of the perturbation. A special r = r + ∆tfast v (19)
2
method has been added to our implementation to modulate κ de- afast = afast (r) (20)
pending on external perturbations following Mikkola & Aarseth
1
(1996). A detailed overview is presented in Appendix A.3, and v = v + ∆tfast afast , (21)
a unit test performed to validate the implementation is presented 6
in Appendix C. where ãfast is the extrapolated gradient acceleration needed to
The Slow Down method can easily be mixed up with Algo- properly cancel 2nd order error terms in the FSI. We compute
rithmic regularisation, which then gives what is called a SDAR this gradient acceleration by evaluating forces twice in a row
integration method (Wang et al. 2020b). This method gives an with a slightly shifted position during the second evaluation.
efficient integrator for multiple systems, which are common in The SDAR method is not meant to integrate the motion of an
young stellar clusters. It can be used in combination with stan- entire cluster. It is specialised in few body evolution. Detection
dard integrators such as the FSI presented above to form a com- of subgroups of close stars is therefore mandatory to apply the
plete stellar dynamics framework capable of evolving an entire SDAR method only on small groups. Gas particles are excluded
cluster of stars for a long time. from such integration methods as their gravitational softening
and hydrodynamic forces can introduce unexpected behaviours
2.3. Phantom new N-body Implementation in the integration. The group finding algorithm is directly derived
from the one proposed in Rantala et al. (2023). Star neighbours
We implemented the integration methods described above into are identified by searching them within a fixed searching radius
the public version of Phantom. They could not fit directly into the rneigh around each star. The typical value used in this work was
main integration routine without a major refactor. Importantly, equal to 0.001 pc. A special criterion allows distant fly-by stars
Article number, page 5 of 21
 A&A proofs: manuscript no. aa54938-25

to enter a group if they enter the searching radius rneigh at the of stars. Instead, we adopt a low numerical resolution and use a
next step. It gives star formation scheme below this scale. This section provides a
r detailed description of this model
r ≥ rneigh and C < ∆tfast , (22)
v
where r and v are the relative position and velocity between the 3.1.1. Sink particles and resolution criteria
two stars tested, C is a safety factor usually set to 0.1. This neigh-
The star formation scheme is based on dynamic sink particle
bour search is then used to construct the adjacency matrix of
creation, already implemented in Phantom. It uses the criteria
the entire system of stars in the simulation. Subgroups are then
proposed in Bate et al. (1995). Collapsing clumps that reach a
identified by searching connected components associated with
user-defined threshold density ρc will be eligible to pass a series
this matrix using a depth-first search algorithm. On-the-fly bi-
of tests to verify that this clump will indeed quickly collapse. For
nary detection is applied to each newly formed subgroup, which
example, a clump needs to be at least in marginal virial equilib-
is essential for the slowdown method. This identification is done
rium. Formally, this gives
every substep to update subgroups as the system evolves.
Subgroups can then be integrated with SDAR methods sep- α j + βrot ≤ 1, (28)
arately from the entire system during the FSI drift phase. How-
ever, as explained in Rantala et al. (2023), each star in subgroups where α j and βrot are the ratios of thermal and rotational energy
still receives perturbations (from gas and point masses in Phan- over the gravitational energy of the clump, respectively. If the
tom) from their environment during the FSI kick phase. It then tests pass, the clump is then directly accreted onto a massive
gives the full new integration scheme in Phantom given by particle called a sink particle created at the centre of mass of the
clump. The sink particle has an associated volume defined by a
v = v + 21 ∆tsph asph (23) user-chosen accretion radius racc , within which gas particles are
v=v+ 6 ∆tfast afast
1 directly accreted.



r=r+ 2 ∆tfast v
An essential point that has long been debated in the literature


 1


rsub = SDAR(rsub , vsub ) (Bate & Burkert 1997; Bate et al. 2002; Hubber et al. 2006) is the





∆tfast choice of this accretion radius together with the critical density

 2
ãfast = afast (r∗ ); r∗ = r +

24 afast (r)



 value. The latter is typically chosen to resolve the initial Jeans
N substeps  v = v + 23 ∆tfast ãfast (24)

 mass of the clump to avoid numerical artefacts that could artifi-
rsub = SDAR(rsub , vsub )

cially fragment the collapsing clump. Thus, the maximal ρc and




r = r + 12 ∆tfast v


minimal racc are linked and depend on the mass resolution of the



afast = afast (r)


simulation. Following Hubber et al. (2006), this link is expressed



v = v + 16 ∆tfast afast


 as
asph = asph (r) (25) Nneigh m ≤ MJeans , (29)
v=v+ 2 ∆tsph asph ,
1
(26)
where Nneigh is the constant number of neighbours of a parti-
As the original Phantom integration scheme, the first integration cle in an adaptive smoothing length framework (around 60 in
level uses an individual time step scheme. However, the substep Phantom) and m is the mass of one SPH particle. Therefore, the
integration sticks with a global time step ∆tfast chosen to be the quantity Nneigh m gives the minimum mass contained in an SPH
most stringent constraint induced by fast forces. Subgroups only kernel surrounding one particle, thereby corresponding to the
evolve during the SDAR phase and are excluded from the stan- minimum resolved clump mass in the simulation. As mentioned
dard FSI drift routine. The SDAR integration routine uses the above, a clump that can turn into a sink needs to be resolved, i.e.,
TTL regularisation method in association with a 6th or 8th-order Mclump = Mjeans, min . Using equation (29) and the expression of
Leapfrog scheme derived from Yoshida (1990) following Wang the Jeans length
et al. (2020b). On top of that, a slowdown is applied to bina- s
ries inside subgroups. For a group composed of only one binary, πc2s
the slowdown factor κ is computed using outside perturbations λJeans = , (30)
Gρ
and is constant during the whole SDAR integration. For multi-
ple system groups, κ is computed with outside perturbations and One can find the maximum critical density:
perturbations induced by other group members. This κ needs to
!3 !2
be updated at each SDAR step as inside perturbations evolve fol- πc2s π
lowing the group dynamics. A correction of the new variable W ρc,max = . (31)
G 6Nneigh m
introduced in the TTL needs to be applied at each kappa update
to conserve all system properties. This correction has the form This maximum value can then be used to find the minimum ac-
of cretion radius. Eq. (29) can be rewritten using the Jeans length
as
W = W + (Ω(κ(t + dt)) − Ω(κ(t))), (27)
where Ω is the gravitational energy of the subgroup. λJeans ≥ 4h, (32)
where h is the smoothing length. Here, the number 4 is derived
3. Physical models from the size of the SPH kernel support which is equal to 2 for
the default cubic spline kernel in Phantom. The minimal accre-
3.1. Star formation prescription tion radius racc, min linked with ρcmax can then be expressed as
In order to save computational time and enable the statistical s
analysis of the evolution of embedded clusters, we do not fol- 1 πc2s
racc, min = . (33)
low the molecular cloud collapse all the way to the formation 2 Gρc,max
Article number, page 6 of 21
 Y. Bernard et al.: DAWN

Sink particles are often used to model one single star in for- radiation, influence it. These processes depend on the local en-
mation (Bate et al. 2002; Lee & Hennebelle 2018), which im- vironment and history of each core. Hence, it is very hard to
plies resolving the opacity limit for fragmentation around a few establish a canonical value. Dunham et al. (2015) estimated the
Jupiter masses (Low & Lynden-Bell 1976). This can only be accretion time for pre- and protostellar phases with differences
done for relatively small initial molecular clouds to keep the of a few orders of magnitude from 5 × 104 to a bit less than
computational time reasonable. For example, the initial mass of 1 × 106 years. These estimates are mainly for low mass objects;
the cloud in Bate & Bonnell (2005) was set to 50 M⊙ using 3.5 constraints on high mass cores are even looser. Here, we fix the
million SPH particles, allowing them to resolve the minimum accretion time tacc to an arbitrary median value of 0.5 Myr. This
Jeans mass in their simulations (0.0011 M⊙ ) with about 80 SPH is sufficiently long to accumulate material to form high-mass
particles. The sink particle accretion radius was correspondingly stars, but sufficiently short to ensure the dynamical stability of
set to 5 au, which is the typical size of the first Larson core star seeds inside the sink.
(Larson 1969). It ensures that only one star will form in a sink. Sinks can accrete during tacc , but stellar objects inside them
Following star formation down to this scale within an entire gi- do not form immediately after sink formation because the frag-
ant molecular cloud would require a tremendous number of SPH mentation process inside a gas clump is not instantaneous com-
particles and is currently beyond the reach of CPU-based codes pared to the accretion phase. Knowing if seeds are already cre-
like Phantom. ated is essential in a merger of two sinks. A qfirst guess for the
Other methods have been proposed to describe star forma- delay time could be the free-fall time τff = 32Gρ 3π
of the clump
tion within GMC following what was done in galaxy-scale sim-
ulations (Renaud et al. 2017). In these simulations, sink particles just before the sink formation, where ρ is the critical density for
are sufficiently large to represent small clusters of stars that form sink formation. In this work this value is around 1 × 10−18 g cm−3
together in the same volume. A similar approach recently devel- which gives a clump free fall time around 7 × 104 yr. The forma-
oped in Wall et al. (2019) considers sink particles as a stellar tion of gravitationally hard stellar cores may take longer. After
forge that produces stars from a list of masses sampled at the reaching a certain density where heating and feedback occur, the
birth of each sink. As gas accretes onto these sinks, stars are infall is not free fall anymore and a slightly higher value could
generated at sink positions following the list of masses. Liow then be more accurate. We set the formation time of stellar ob-
et al. (2022) proposed to group several sinks on the same list to jects in sinks tseed to 0.1 Myr. This time is analogous with tcollapse
produce more massive stars that cannot be generated by just one. in Bleuler & Teyssier (2014).
However, these methods introduce unwanted behaviour. For Newborn stars are then released in the cluster in formation
example, when stars are generated in a sink, the sink remains in after a sink reaches its maximum accretion time. To summarise,
the simulation with almost no mass in its reservoir. When ac- this new star formation procedure using sink particles in Phan-
cretion restarts, the sink position and velocity are mass-averaged tom is illustrated in Fig. 1 and proceeds as follows:
with new accreted gas. As the sink mass is near zero, its posi- 1. First, a sink is formed when a clump of gas reaches certain
tion jumps to the accreted gas position. Due to the number of conditions;
stars generated, sink particles in these simulations exhibited ar- 2. The sink accretes surrounding materials continuously ;
tificial Brownian-like motions. The erratic motion of these sinks 3. At the time tseed , the number of star seeds inside the sink is
can completely change the gas infall geometry and then perturb sampled randomly between 1 and 5. Accretion continues;
the spatial distribution of stars in the cluster. 4. When tacc is reached, accretion is over. The sink is dissolved
and young stars are released within its volume, keeping its
3.1.2. New star formation prescription from sink particles momentum. Star masses are set by sharing randomly the par-
ent sink mass reservoir. Positions and velocities are sampled
We propose a new prescription to model star formation using to respect the dynamical stability of the group as described
“large” sink particles where the fragmentation and accretion pro- below;
cesses are not yet over. Typically, accretion radii are set to a few 5. If sinks intersect each other by racc during their lifetime, they
thousand au in this prescription. In that sense, sinks can be com- merge to form one single sink. In that case, if star seeds have
pared to the dense cores observed in star-forming regions. already formed, then their number will be summed up into
However, these sinks are not necessarily cores from an ob- the merged one. If not, only the mass and the momentum of
servational point of view. They start their life as collapsing or both sinks are added to the merged sink.
marginally stable clumps. However, their future evolution is
based on the sub-sink scale physics that will be described in the As said above, stars released after accretion onto their
following. This evolution can diverge from what is called a core parental sink need to be distributed in their birth vicinity. Each
in observations. Sadavoy & Stahler (2017) studied YSOs within new star will have a mass, position, and velocity that need to
hot cores of a radius roughly equal to 10 000 au in the Perseus respect the main hypothesis of this prescription. There is conser-
molecular cloud. They concluded that only a few (up to 5) YSOs vation between sink mass and star mass. Therefore, sink mass is
formed in the surveyed cores where they saw multiples (24 em- shared between child stars. The simplest way to share sink mass
bedded cores). These results are used as the main hypothesis for between stars is a random sharing with a cut at 0.08 M⊙ to en-
this prescription. A sink particle is considered as a “core” of a sure the formation of stars and avoid brown dwarf and substellar
few thousand au where the fragmentation process will occur to mass objects. We implement this by sampling nseed − 1 random
form several (up to 5) YSOs. During a specified time the sink numbers between zero and one. That way, a unit segment is sub-
will continue to accrete infalling material that will feed the star divided into multiple sub-segments of random length, and their
“seeds” inside it. Nevertheless, physical processes inside a sink sum is always equal to one. Multiplying these segments by the
are not computed in the simulation (star seeds dynamics and ac- mass of the sinks then produces a set of masses randomly chosen
cretion). from the available mass in the reservoir. In addition to that, we
The accretion time is hard to constrain, especially as a wide resampled some points to ensure that no sub-segments are lower
range of processes, such as outflows, dynamical encounters and than 0.08 M⊙ .
Article number, page 7 of 21
 A&A proofs: manuscript no. aa54938-25

Gas accretion gas reaches the condition to form a sink in the vicinity of a star,
the sink is still allowed to form.

Sink formation 3.2. Stellar feedback: H ii regions

tests
Newborn stars inside active star-forming regions are connected
with their natal environment. They inject both kinetic and ther-
mal energy into their surroundings (referred to as stellar feed-
back; Krumholz et al. 2014), which can influence future star
Sink particle formation. Examples of momentum-driven feedback are jets or
winds. Thermal feedback occurs from stellar irradiation.
Dense clump of gas
Irradiation consists of two main components: ionising, and
Accretion phase non-ionising radiation. Non-ionising radiation from direct stars
or dust-reprocessed light can heat the gas via absorption. It oc-
curs mostly at infrared wavelengths around dense cores hosting
protostars, where the gas starts to be optically thick (Bate 2009).
Ionising radiation is produced by massive stars (> 8 M⊙ ) inside
molecular clouds. Photons with energy larger than ELy = 13.2eV
Stars release
can directly ionise hydrogen atoms. Ionised gas around mas-
sive stars is heated to 10 000 K (the result of heating and cool-
Sink boundaries removed
ing equilibrium). It creates a strong temperature gradient with
cold neutral gas inside the region, which consequently produces
Fig. 1: The star formation process described by our new pre- a strong pressure gradient. Thus, ionising radiation from mas-
scription, from the clump transformation into sink particles to sive stars produces bubbles of hot gas around them that will ex-
stars release inside the simulation domain. This process is sub- pand inside the star-forming region, creating what is called an
divided into three parts. First, a collapsing clump of gas passes H ii region (Osterbrock & Ferland 2006). This is thought to drive
tests to be transformed into a sink particle. The sink particle can gas expulsion from star-forming regions, a key element of young
accrete during tacc . After reaching this time, the sink is dissolved stellar cluster evolution.
and its mass is shared to produce newborn stars at its last loca- N-body simulations (Kroupa et al. 2001) tried to model gas
tion. expulsion with simple methods like an exponential decay of a
gravitational potential inside the star cluster. Since we model
the gas directly, it is possible to model this feedback process by
heating the gas around massive stars inside the newborn cluster.
Finally, star velocities and positions should produce a bound Models trying to reproduce H ii regions already exist and have
system inside the sink volume. Many solutions exist to sample been used in numerical simulations from the Galactic scale to
these quantities and the sampling method will greatly influence the near-neighbourhood of massive stars. They can be modelled
the resulting cluster morphology. The choice thus plays a key with more or less precision depending on trade-offs that authors
role in the simulation. However, exploring such a large param- made between accuracy and performance in agreement with their
eter space would be computationally expensive and hard to im- physical constraints. For example, Hopkins et al. (2012) used a
plement on the first try. We opt for a simplistic approach using coarse model to implement H ii region feedback from entire star-
the same method that produces Plummer-like clusters in stellar forming regions to probe their effects on a global galaxy simula-
dynamics-oriented simulations (Aarseth et al. 1974). It can pro- tion without killing numerical performance. At the more sophis-
duce spherical, centrally peaked, virialised star distributions with ticated end, radiation transfer specialised codes are used in Dale
a very small development cost. After their formation, stars are et al. (2012) and Verliat et al. (2022) to precisely model ion-
not able to accrete material and are considered as point masses. ising photon propagation inside a specific star-forming region
to resolve the coupling between star formation and stellar feed-
Stars and sink particles can coexist during the simulation. back. However, such schemes come with a large numerical cost.
Their time integration is handled the same way in the code. We Since our aim is to be computationally efficient while conserv-
use the original method from Dehnen & Read (2011) as de- ing individual star resolution, we want to model each H ii region
scribed in Price et al. (2018) to find the minimum time step con- produced by massive stars inside one simulation.
straint imposed by sink-sink, star-star or sink-star interactions. To achieve this, we implemented a model similar to that
This one is used for the substep integration if it is the smallest of Fujii et al. (2021) and González-Samaniego & Vazquez-
of all forces computed during this phase. If these particles be- Semadeni (2020) to simulate individual massive star H ii regions
long to a subgroup, then the integration is decoupled and han- inside star-forming regions. This method is based on a distri-
dled by subgroup routines. To properly handle the collisional bution of the ionising photon rate produced by massive stars
dynamics of stars, no softening is allowed in sink-sink, star-star (> 8 M⊙ ) at each SPH particle around them, assuming ionisa-
or sink-star interactions. However, interactions between sink/star tion equilibrium inside the hot region and neglecting diffusion.
and SPH particles are softened as the latter cannot be considered If the ionising rate is sufficient to completely ionise a particle,
as point-mass particles. The softening length rsoft-gas is chosen the particle is tagged as ionised and its temperature and sound
to be the maximum of the gas smoothing length and the accre- speed are changed accordingly. If not, it means that the ionisa-
tion radius racc of sink particles (softening lengths are equal to tion front is nearby and all particles inside the H ii region have
rsoft = 4000 au for stars). Softening these forces avoids creating a been tagged as ionised.
large concentration of dense gas around stars like discs or infall, Each massive star (> 8 M⊙ ) ionising photon rate ∆N is com-
preventing a slowdown of the simulation. However, if a clump of puted using a fit function of the mass proposed in Fujii et al.
Article number, page 8 of 21
 Y. Bernard et al.: DAWN

(2021) and based on observational surveys. Then, at each simu- Eq. (33), which can then give a minimum mass resolution using
lation update, an iterative algorithm verifies if the nearest SPH Eq. (31) to resolve the collapse of clumps. The maximum critical
particles of a massive star are ionised or not. The algorithm com- density is equal to ρc = 1×10−18 g cm−3 using an accretion radius
putes the ionising photon rate to fully ionise the massive particle. of 4000 au. The coarsest mass resolution is around 0.009 M⊙ ,
Mc
Following Hopkins et al. (2012), it gives which gives a minimum SPH particle number of ∆M ≈ 1.1 mil-
mi lion. Preliminary tests showed that simulations close to this min-
δNi = αβ ni , (34) imum number gave unusual accretion behaviours. A low number
µm p
of particles cannot sufficiently retrieve small turbulence modes
where mi is the particle mass, µ the mean molecular weight, m p , causing faster artificial damping of turbulence, leading to early
the proton mass, αβ ≈ 2.3 × 10−13 cm3 s−1 the recombination rate, collapse. We therefore set the resolution to 3.5 million particles
and n p the number density. If δN p < ∆N the particle is ionised for this fiducial model. Finally, we set the sink accretion time to
and δN p is subtracted from ∆N. Then the algorithm repeats on 0.5 Myr, and the star seeds formation time to 0.1 Myr, as dis-
the next nearest particle until δN p > ∆N. Reaching this point, the cussed in Sect. 3.1.
H ii region is sampled. The algorithm is applied to each massive Sinks are created after the density exceeds the critical density
star in the simulation. To help the algorithm, nearest neighbour and the tests to ensure that the clump collapses pass. This means
lists of massive stars are constructed using the tree-walking rou- that regions of cloud exceeding the critical density are not nec-
tine from Phantom. essarily directly transformed into sinks. These tests avoid con-
In contrast to the method described by Bisbas et al. (2009) or fusing isothermal compression that can increase the gas density
Dale et al. (2012), we do not implement ray-tracing. Our method without gravitational collapse. However, such isothermal com-
generates H ii regions efficiently but also produces side effects pression can be expensive in numerical integrations. To acceler-
when the mass resolution is not high enough to open H ii regions ate the computation, one can force the sink creation after exceed-
in a very dense clump of gas or if a region starts to grow inside ing the critical density by a certain factor. In this model, the sink
void (due to free boundaries in SPH). The first side effect can creation forcing is set to foverride = 100ρc . Table 1 summarises
easily be avoided by increasing the resolution and/or attributing all parameters of this model.
a fraction of the ionising temperature to the nearest particle to
help the gas expand and consequently open the H ii region. In
the second side effect case, all ionising radiation will be focused Parameter Value
where the gas is, as the iterative algorithm does not allow ra- Mc 104 M⊙
diation escape. It could then potentially lead to an overestimate Rc 10 pc
of the H ii region. Our implementation follows what Fujii et al. αt 2.0
(2021) proposed to minimise this effect by limiting H ii regions T 10 K
to a size near the initial radius of the molecular cloud. Our tests to τff 5.3 Myr
ensure convergence of the model following Bisbas et al. (2015) ∆m 0.00286 M⊙
are presented in Appendix B. racc 4000 au
ρc 1 × 10−18 g cm−3
tacc 0.5 Myr
3.3. Fiducial model
tseed 0.1 Myr
We studied one specific fudicial model of a giant massive molec- foverride 100ρc
ular cloud collapse to characterise the new star formation pre-
scription and the resulting star cluster. This model starts from Table 1: Main parameters of the fiducial model. Mc is the total
an isolated isothermal massive turbulent molecular cloud. This initial mass of the cloud and Rc is its initial radius. αt is the
cloud is initially spherical with a radius Rc = 10 pc and a total initial thermal virial ratio and T is the isothermal temperature
mass Mc = 104 M⊙ . A turbulent velocity field is initially injected of the cloud. Its initial free-fall time is given by τff . Each SPH
into the cloud. We generated this using methods from Make- particle composing the total cloud has a mass of ∆m. racc is the
Cloud (Grudić 2021). It produces initial turbulence by applying sink accretion radius that is connected to the critical density ρc
a Gaussian random field filtered by a Larson-like power spec- via ∆m. Respectively to the new star formation scheme, sinks
trum (Larson 1981) on the particle velocity. Compressive and can accrete during a fixed time tacc and seeds inside it form after
solenoidal modes are introduced with 30 % and 70 % proportion, reaching tseed . Finally, sink creation tests can be overridden if the
respectively. This field is finally scaled to reach a turbulent virial density reaches a too high value. This override is controlled by
ratio αt = 2 (marginally stable condition). The cloud’s initial foverride .
temperature is set to 10 K, and the mean molecular weight is set
to µ = 2.35. This model ignores the magnetic field and external
driving (pressure and turbulence). As the cloud is supposed to
be isothermal, no cooling/heating processes are included, except
photoionisation.
For the star formation prescription, accretion radius racc and 4. Results
critical density ρcrit need to fulfil resolution criteria to ensure that
gravitational collapse is resolved up to the sink scale. These two In this study, we evolved three simulations of this fiducial model
values then need to be chosen in agreement with the mass reso- with different initial turbulence states to 6.5–7 Myr. We denote
lution of the simulation that gives the number of SPH particles. the three simulations as Embedded Cluster Formation 1, 2 and 3
Our sink particle accretion radii should represent the character- (hereafter ECF1, ECF2 and ECF3) in the subsequent sections.
istic size of a core in star-forming regions of a few thousand Each simulation was computed in around 14 000 CPU hours,
au (Sadavoy & Stahler 2017). Thus, we set the accretion ra- which corresponds to around 18 days of computation on one
dius to 4000 au. The critical density can then be computed using compute node (32 cores).
Article number, page 9 of 21
 A&A proofs: manuscript no. aa54938-25

4.1. Qualitative description The formation of massive stars, their spatial locations, and
their coherent contributions are key elements in the future of
Figure 2 shows all three simulations at different timestamps of a young stellar cluster that will be sculpted by H ii region ex-
their evolution, from the initial gas fall up to the star cluster pansion. Their locations and masses are function of the initial
emergence. Each row and column corresponds respectively to a turbulent velocity field and mass sharing inside sinks, which is
specific simulation (top to bottom) and timestamp (left to right). random by definition. It means that the formation of such clus-
Snapshots display the column density of the gas overlaid by ters is, like their evolution, highly stochastic. This confirms that
sinks and star positions. a statistical study is necessary to identify the main evolutionary
As expected, the early evolution (t = 3 Myr) of each cloud is paths of embedded clusters.
similar in each case. They start to diverge after the first star for-
mation events. The first massive star formation presents a turning
4.2. Star formation efficiency
point in their evolution. These massive stars (> 8 M⊙ ) open large
H ii regions around them. Their location and efficiency dramati- The gas mass conversion and its evolution in time can help char-
cally change the subsequent cloud collapse. acterise the star formation inside clouds. Dividing the converted
The first massive star is formed early in ECF1, around mass by the initial gas mass gives an efficiency factor referred
1.5 Myr before reaching the initial free fall time. The resulting to as the star formation efficiency (SFE). This coefficient is an
H ii region pushes the gas outward efficiently, inducing a signifi- indicator of the evolutionary stage of a star-forming region and
cant compression in the cold molecular gas, causing it to collapse is computed at various scales in observational works, from 0.01
and trigger secondary star formation. In the meantime, other re- to 10 pc (Alves et al. 2007; Evans et al. 2014; Louvet et al. 2014;
gions of the cloud continue their collapse and eventually pro- Könyves et al. 2015a). From an entire giant molecular cloud to a
duce other massive stars able to efficiently heat and ionise the smaller scale, like in the neighbourhood of dense cores inside a
cold gas with their strong radiation. Their accumulation accel- larger complex, the definition of the SFE is the same and is given
erates the gas expulsion, leaving behind a new cluster of stars by
exposed from their parental cloud. Triggered formation follow-
ing expanding bubbles is imprinted with their motion, resulting M⋆
SFE = , (35)
in a global expansion of the newly formed cluster. Mcloud + M⋆
ECF3 shows similar outcomes from ECF1, with a large ex- where M∗ and Mcloud are the total mass of stars and the total
panding H ii region hosting secondary formation at its outskirts. mass of gas inside the considered region, respectively. In our
However, gas expulsion is more efficient. The final cluster ap- study, star formation efficiency is computed at the cloud scale.
pears correspondingly looser and more spatially extended. It Hence, no external mass is added to the calculation, which may
did not produce massive stars early. Instead, the cloud collapse not be the case in observations due to external gas inflows and
had more time to concentrate its mass in a central region. This protostellar outflows.
stronger collapse produces multiple massive stars in a short time Figures 3, 4, and 5 show for each simulation the evolution
interval and in the same area. The coherent addition of H ii re- over time of key parameters linked to the gas mass conversion
gions leads to the most efficient gas expel and the highest cluster into stars, namely the star formation efficiency (top panel), the
expansion rate seen in this work. median values of sinks mass and accretion rate (bottom left and
ECF2 produces a highly off-centred intermediate massive centre), and the number of sinks and stars (bottom right). Each
star early in the simulation in comparison to the other two. Its simulation reaches a SFE exceeding 10%, with values between
H ii region expands slowly and pushes the gas gently, helping 14% and 18%. These efficiencies are high compared to those ob-
it to collapse naturally under its weight. Moreover, later-formed served in parsec-scale star-forming regions (see however Lou-
massive stars do not show coherent effects on their H ii regions vet et al. 2014). Such high efficiencies tend to be found at much
as in ECF3, even if these stars are the most massive ones in smaller scales (e.g., Könyves et al. 2015a, at 0.01 pc scales). This
all simulations. Instead, shocks are produced between H ii re- suggests that other feedback mechanisms, such as jets (Machida
gions, thereby facilitating the condensation of cold gas to an even & Matsumoto 2012), or other physics such as non-ionising radi-
greater extent. With the help of gravitational collapse, this leads ation and magnetic fields (Price & Bate 2009), may be needed to
to the formation of a very compact cluster that is not yet fully regulate the mass conversion. However, it is difficult to perform
revealed from its natal gas and where star formation is still on- a direct comparison as real molecular clouds are certainly not in
going at the end of the computation (∼ 6.6 Myr). complete isolation and can gather gas inflow, which will reduce
The triggered formation has already been seen and studied the SFE.
in numerical simulations. Dale et al. (2012) performed molecu- Most of the gas conversion occurs in ∼ 2 Myr, between 4
lar cloud collapse with H ii feedback relatively similar to ours. and 6 Myr (0.75 to 1.13 τff ). The first sinks appear around 3 Myr
They found in their simulation that ionising radiation from mas- (0.56 τff ) and are limited to a small number during one million
sive stars shapes the gas inside such molecular clouds by creat- years. The mean star formation rate (SFR), given by the mean
ing large bubbles and pillars, which are also observed here. They slope of the mass conversion evolution in time, is equal to 6.8 ×
also concluded that this radiation feedback induced a significant 10−4 M⊙ yr−1 for the first two runs. ECF3 shows a slightly lower
rise in star formation, which also corroborates what is found in value with 4.5 × 10−4 M⊙ yr−1 . All three mass conversion slopes
our simulations. However, their calculation stops rapidly (Model increase monotonously during the whole computation.
J) due to growing time step constraints on the stellar dynamics. It is difficult to see where star formation halts from the num-
Our new methods allow us to continue the calculation further. ber of sinks over time. In ECF2 and ECF3, the number of sinks
Thanks to this amelioration, we can now see that this triggered starts to decrease at the last few points. This, in addition to the
star formation can still happen far from the revealed cluster, cre- SFE slowdown in both simulations, suggests that star formation
ating an expanding distribution of low-mass stars around denser is about to stop. By contrast, the number of sinks is still increas-
clusters. ing in ECF1 and the star formation efficiency does not show any
Article number, page 10 of 21
 Y. Bernard et al.: DAWN

ECF1 ECF2 ECF3

15 15 15
time : 3.006 Myr time : 3.006 Myr time : 3.006 Myr

10 10 10

5 5 5

y [pc]
y [pc]

y [pc]
0 0 0

5 5 5

10 10 10

15 15 15
15 10 5 0 5 10 15 15 10 5 0 5 10 15 15 10 5 0 5 10 15
x [pc] x [pc] x [pc]

15 15 15
time : 4.714 Myr time : 4.714 Myr time : 4.714 Myr
10 10 10

5 5 5
y [pc]

y [pc]

y [pc]
0 0 0

5 5 5

10 10 10

15 15 15
15 10 5 0 5 10 15 15 10 5 0 5 10 15 15 10 5 0 5 10 15
x [pc] x [pc] x [pc]
15 15 15
time : 5.510 Myr time : 5.510 Myr time : 5.510 Myr
10 10 10

5 5 5
y [pc]

y [pc]

y [pc]
0 0 0

5 5 5

10 10 10

15 15 15
15 10 5 0 5 10 15 15 10 5 0 5 10 15 15 10 5 0 5 10 15
x [pc] x [pc] x [pc]
15 15 15
time : 6.327 Myr time : 6.610 Myr time : 6.914 Myr

10 10 10

5 5 5
y [pc]
y [pc]
y [pc]

0 0 0

5 5 5

10 10 10

15 15 15
15 10 5 0 5 10 15 15 10 5 0 5 10 15 15 10 5 0 5 10 15
x [pc] x [pc] x [pc]
0.5

1.0

1.5

2.0

2.5

3.0

3.5

0.5

1.0

1.5

2.0

2.5

3.0

3.5
0.5

1.0

1.5

2.0

2.5

3.0

3.5

Log(column density) [Msun/pc²] Log(column density) [Msun/pc²] Log(column density) [Msun/pc²]

Fig. 2: Snapshots showing the evolution of three fiducial clouds computed at multiple time steps. Each column corresponds to a
specific run (top: ECF1, middle: ECF2, bottom: ECF3) and rows give different epochs in each of them. The colour map gives the z
integrated column density of the gas in M⊙ pc−2 , while red circles correspond to sink positions and stars are coloured using realistic
black body temperature extrapolated from a luminosity-mass relation. ECF1 presents a protostellar cluster under the influence of
H ii regions distributed in various parts of the molecular cloud, like for instance Aquila. ECF2 presents a more centrally concentrated
protostellar cluster, which develops a cluster of H ii regions, like the region W43. ECF3 is more difficult to link with an existing
star-forming region. As all the gas has been removed very rapidly, the resulting group of stars looks like a more evolved association
or cluster like NGC 6611.
Article number, page 11 of 21
 A&A proofs: manuscript no. aa54938-25

slowdown. This late increase is only due to secondary star forma- who link lower star formation efficiency to efficient photoionisa-
tion in H ii regions outskirts. It could also be induced by numeri- tion feedback.
cal effects. It seems clear that simulations should be evolved fur-
ther in time to determine when star formation stops. (discussed
in section 5) 5. Discussion
A key event occurs between 4 and 4.5 Myr (0.75 to 0.84 τff ) 5.1. Numerical performance
for all calculations, when the star formation process produces
a sharp rise in the number of sinks produced in the cloud. Ex- In terms of computational efficiency, we found our new N-body
cept for ECF3, this moment coincides with the first formation of and H ii region implementations to be effective at speeding up the
massive stars. calculation when the stellar environment densifies. Our new sub-
The median values of sink mass and accretion rate show two group methods correctly handle every close encounter or binary.
different regimes at the same turning point of 4Myr for ECF1 and ECF2, which shows the densest cluster formed in this study, pro-
ECF2. In ECF1, before this time, the mean value of the median duces a maximum of 17 subgroups involving 38 stars at the same
for both plots is equal to 10 M⊙ yr−1 and 3 M⊙ . These values drop snapshot. Even if this number is small compared to the total
to 2 M⊙ yr−1 and 0.8 M⊙ just after the turning point. The number number of stars at the end of the simulation (N f inal = 1225), the
of sinks also drastically increases from less than ten to more than performance boost is significant, as only one of these subgroups
a hundred in less than a million years after this timestamp. This can impose a timestep that is orders of magnitude lower in the
sudden jump occurs due to the opening of an H ii region. substep loop. We indeed observe one bound group which has a
time step a hundred times smaller than the minimal timestep with
Expanding bubbles in these simulations greatly modify the no subgroup integration with a minimum separation of 4 au.
star formation regime making a sharp transition from a domi- To properly demonstrate this performance boost, we re-
nant massive formation channel to a low-mass one, as observed launched the simulation near the time of creation of this hard
in some evolved massive protoclusters (e.g. Armante et al. 2024). binary. This new calculation took 48 hours (walltime) to com-
This new channel produced stars in the outskirts of H ii bub- plete only one-fifth of the time integration until the next dump.
bles where sinks and gas rapidly decouple after their formation. In comparison, the whole time integration took 3.1 hours with
As sinks form, these only evolve through gravitational interac- our new implementation, which gives a gain of approximately
tion and lose pressure support that drives H ii bubble expansion. 80. The gain depends on the hardest binary in the system.
Thus, if this expansion accelerates, sinks do not have the time to With our methods, time step constraints set by individual
accrete material before entering an ionised region where accre- stars (not in subgroups) are always close to the CFL (Courant-
tion is no longer possible. Most of the converted mass after the Friedrich-Lewy) constraint (Courant et al. 1928) of the densest
turning point then goes mostly into low-mass sinks, explaining SPH particles in the cloud during the simulation, which sets the
the sharp transition. slow force minimum time step, i.e., the number of substeps is
In unperturbed locations of the cloud, massive stars still form always contained to a reasonable number (a few tens at maxi-
following the natural cloud gravitational collapse. A revival of mum). Hence, none of these parts of the calculation is a major
both accretion rate and mass can be seen to start at ∼ 0.95 τff in bottleneck.
ECF1 and ECF2. And it reaches a maximum at ∼ 1.05 τff that The main performance hit is due to the interactions between
coincides with the most massive sink formation in both simula- point masses (sink and star) and the gas. These are calculated
tions. Higher revival in ECF2 also agrees with previous obser- with a direct method in Phantom and always at the lowest time
vations (see section 4.1) which showed that H ii bubble expan- step of the scheme. In most applications of Phantom to date, this
sion was more gentle, leading to a forced collapse instead of a part of the algorithm was not the main bottleneck, as the sink
dismantling. ECF2 produces the most massive sink and then the number was not large enough for the O(Nsink Ngas ) interaction to
most massive star with a mass of 180 M⊙ and 45 M⊙ respectively. dominate the complexity. Such interactions essentially scale lin-
It seems logical that regions unaffected by H ii bubbles will col- early with particle numbers in our simulations, like in previous
lapse very rapidly and efficiently near the initial free-fall time. studies. As a result, this bottleneck makes the computation too
The behaviour of the ECF1 and ECF2 runs reminds the Aquila expensive to go further in time than 6 to 7 Myr for our applica-
and W43 protoclusters (Könyves et al. 2015b; Motte et al. 2003). tion.
In comparison, ECF3 did not produce early massive stars Since the number of point masses is close to three thou-
like the other two. This particular simulation shows a smoother sand in our simulations, these interactions start to dominate the
formation process up to τff , maybe resembling the star forma- computational expense. As said above, they are also computed
tion activity in nearby, low-mass star-forming regions. It fol- for every point mass at every time step, even when only a few
lows the gravitational collapse of the cloud that naturally acceler- gas particles are active in the individual time stepping scheme.
ates. Thus, the number of formed sinks increases smoothly. The It is a waste of computation time, as the time step constraints
lack of forcing by the H ii region prevents the sharp transition at of sink/star and the gas interactions are often two or three or-
0.75 τff as in the other runs. The acceleration of the collapse stops ders of magnitude higher than the smallest one. On average, this
sharply at ∼ 1.05 τff . As many massive sinks were able to cluster part of the calculation represented more than 80 % of the com-
in dense regions without any feedback, this led to the simultane- putation time of one global step when the number of sinks is
ous formation of massive stars in the same region. This burst of higher than ∼ 1000. This bottleneck is the main remaining lim-
massive star formation produces the sharpest regime transition itation that makes it difficult to carry out this type of simulation
observed in the simulations, permitting solely low-mass star for- on a statistical scale, which is the ultimate goal of the project.
mation until the end of the computation. However, this scenario Future work will focus on optimising this part of the calcula-
is probably not the sole possibility if early massive star forma- tions. For example, following what has already been done in
tion does not occur. Finally, the star formation efficiency of this collisional N-body codes (Mukherjee et al. 2021; Wang et al.
simulation is the lowest, which may be due to rapid gas removal. 2020a), these forces — as they do not strongly perturb the evo-
This may correlate with a conclusion made by Dale et al. (2012), lution of hard groups or close encounters — can be computed
Article number, page 12 of 21
 Y. Bernard et al.: DAWN

0.6 0.7 0.8 0.9

time [ ff] 1.0 1.1 1.2 1.3
1750 40
1500 35
1250 30

Mass [Msun]
M [Msun]

1000
25
750
20
500

250
15
0 10
3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0
Time [Myr]
0.5 0.6 0.7
time
0.8
[0.9ff] 1.0 1.1 1.2 0.5 0.6 0.7
time
0.8
[0.9ff] 1.0 1.1 1.2 0.5 0.6 0.7
time
0.8
[0.9ff] 1.0 1.1 1.2
median median stars
max max 103 sinks
103
Msink [Msun/Myr]

Msink [Msun]

102 101 102

count
101 101
100 100
100
3 4 5 6 3 4 5 6 3 4 5 6
time [Myr] time [Myr] time [Myr]
Fig. 3: Evolution over time of multiple key parameters linked to the gas-to-star mass conversion. The top panel shows the time
evolution of the mass of gas converted in sinks and stars in ECF1. The formation of massive stars is represented by vertical lines,
with their colour indicating the mass of the massive object. The bottom panels show the temporal evolution of the accretion rate
(left), mass median (centre) and the number of stars and sinks formed in the simulation (right). For the first two panels, the quartile
range (grey zone) and the maximum value (dashed line) are also shown. It allows us to see that the sharp transition at 4Myr occurs
on the majority of sinks. However, even after the change of regime, the cloud is still capable of producing some massive sinks with
a high accretion rate.

with the help of the Phantom kd-tree reducing the complex- with foverride equal to several orders of magnitude to find a good
ity to O(Nsink log(Ngas ) + Ngas log(Nsink )) and/or by only comput- compromise between efficiency and accuracy.
ing these forces at the first integration level. The latter solution For the sake of efficiency, this parameter was set to 100 in
would make it possible to avoid computing these interactions the simulations of this paper, but this is a relatively low value.
multiple times during the substep scheme, and computing them With this value, HII region expansions produced in the simu-
only once per step should save most of the unnecessary calcula- lations showed efficient sink formation at the outskirts of these
tions. regions. The compression induced by hot bubbles in the cold
collapsing gas can force this collapse and lead to a burst of sink
5.2. Effect of foverride on sink formation output production. However, it can also produce compression in non-
collapsing regions that can exceed the critical density. Therefore,
We investigate the effect of the foverride parameter used to force we find that a too low foverride parameter could overestimate this
the creation of a sink if an SPH particle reaches the critical den- burst of sink production, leading to an overestimated number of
sity without passing other physical tests detailed in Price et al. low-mass stars at the end.
(2018) or Bate et al. (1995). These tests check whether the gas Figure 6 shows histograms of the formation of sinks over
clump is indeed starting to collapse gravitationally, or whether time, and up to one freefall time, of two simulations with dif-
it is just experiencing hydrodynamical compression. If the gas ferent foverride values (ECF1 with foverride = 100 and ECF1 with
is isothermal, such compression can rise several orders of mag- foverride = 1000). ECF1 produced all its sinks, except one, via
nitude higher than the critical density. In our simulations, high forced creation. Looking at ECF1B, we see that the first sink
density implies high time constraints which will drastically slow creations are similar. Thus, the forced creation does not over-
down the computation. Moreover, all this extra cost is unneces- estimate the production of sinks during the initial collapse. The
sary, as the critical density is supposed to be close to the maxi- sink productions of ECF1 and ECF1B start to diverge after the
mum resolution of the simulation. The extra computation above first massive star formation. The total number of sinks produced
the critical density is not guaranteed to be physically accurate. after this timestamp is 263 (only forced) and 141 (79 standard
foverride is therefore here to force the creation of a sink after and 62 forced) for ECF1 and ECF1B, respectively. Compared
reaching to ECF1B, ECF1 overestimates by at least 1.9 times its sink
production at the outskirts of H ii regions with foverride set to
ρ = foverride × ρcrit , (36) 100. While some of the forced sink creations are physical, and
Article number, page 13 of 21
 A&A proofs: manuscript no. aa54938-25

0.6 0.7 0.8 0.9

time [ ff] 1.0 1.1 1.2 1.3
40
1750

1500
35
30

Mass [Msun]
1250
M [Msun]

1000 25
750
20
500

250
15
0 10
3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0
Time [Myr]
0.5 0.6 0.7
time0.9[ ff] 1.0
0.8 1.1 1.2 1.3 0.5 0.6 0.7
time0.9[ ff] 1.0
0.8 1.1 1.2 1.3 0.5 0.6 0.7
time0.9[ ff] 1.0
0.8 1.1 1.2 1.3
104 median median 103 stars
max max sinks
102
103
Msink [Msun/Myr]

102
Msink [Msun]

102

count
101
101 101
100
100
100
3 4 5 6 3 4 5 6 3 4 5 6
time [Myr] time [Myr] time [Myr]
Fig. 4: Same as previous figure applied on ECF2 simulation

not numerical artefacts, a higher value of foverride helps to min- scription hypothesis that permits only a few star seeds (less than
imise this numerical side-effect. Therefore, foverride should be set five in these simulations) in a sink. This means that the most
higher than 100 in future simulations while keeping numerical massive sinks need to share their mass with a large number of
efficiency in mind. stars during the release phase. Randomness during sharing is
smoothed in such cases as it needs to cut the parent mass in more
smaller parts. As a result, most massive sinks create a large num-
5.3. Sink merging ber of intermediate-mass stars. This is at odds with our model,
As mentioned in Sect. 3.1, sinks can collide with others during which predicts that the most massive stars are produced inside
their accretion phase. If two sinks collide and both encompass very massive sinks. Therefore, our random sharing of mass be-
their centre of mass, they fulfil the conditions for a sink merger. tween star seeds can bias the final star mass spectrum by under-
In such a case, the two sinks are merged to form a new sink that sampling high-mass stars and oversampling intermediate-mass
inherits the physical properties of its parents. The new position stars. Since less massive stars also imply less feedback, it can
and velocity are the barycentres between the two parents. This also have an impact on the star formation efficiency and dy-
ensures the conservation of linear momentum during the sim- namics inside the star-forming region. Ultimately, we need to
ulation. The angular momentum is also conserved as Phantom improve the merging process in relation to our new star forma-
tracks the spin of all sinks in the simulation. The sink properties tion scheme to make sure that the number of seeds inside a sink
provided by our new star formation prescription are also updated is always in agreement with the maximum value allowed, e.g.
during a merger (see section 3.1). In particular, the number of by adding star escapers after a merging and/or modelling seeds
seeds in the initial sinks is summed into the merged sink. merging inside our subgrid model.
Mergers are quite common in these simulations. The gravi-
tational pull of the cloud tends to bring sinks closer over time,
6. Conclusions
which inevitably creates collision trajectories between them.
Mergers are key to forming very massive sinks that can then cre- This paper presents a new implementation of direct N-body algo-
ate massive stars. This merging process can take place consecu- rithms as well as prescriptions for star formation and stellar feed-
tively on the same sink, leading to very massive sinks (the most back in the SPH code Phantom. Together, they unlock the pos-
massive one reaches 180 M⊙ in ECF2). Accretion rates measured sibility of studying embedded star cluster evolution from their
in these simulations are too low without merger events to ex- birth inside a molecular cloud to their unveiling from their natal
plain the formation of stars of 20 to 40 M⊙ , which are essential gaseous environment. The new algorithms come in two parts:
to rapidly remove the gas from the star-forming region with their First, our new star formation prescription was implemented,
H ii feedback. which provides distributions of stars self-consistently formed
In these cases, summing the parents’ star seeds can bring from molecular cloud collapse at low computational cost. It uses
their number to a few dozen in a single sink (50 at maximum in thousand-au-diameter sink particles that host star seeds within
ECF2). This can be seen as a contradiction with our main pre- them and are transformed into stars after passing some criteria.
Article number, page 14 of 21
 Y. Bernard et al.: DAWN

0.6 0.7 0.8 0.9

time [ ff] 1.0 1.1 1.2 1.3
40
1400

1200
35
1000 30

Mass [Msun]
M [Msun]

800 25
600
20
400

200
15
0 10
3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0
Time [Myr]
0.6 0.7 0.8
time
0.9
[ 1.0ff] 1.1 1.2 1.3 0.6 0.7 0.8
time
0.9
[ 1.0ff] 1.1 1.2 1.3 0.6 0.7 0.8
time
0.9
[ 1.0ff] 1.1 1.2 1.3
103 median median stars
max max 103 sinks

102
Msink [Msun/Myr]

Msink [Msun]

101 102

count
101
101
100
100
10 1 100
3 4 5 6 7 3 4 5 6 7 3 4 5 6 7
time [Myr] time [Myr] time [Myr]
Fig. 5: Same as previous figure applied on ECF3 simulation

70 70
standard standard
forced forced
60 60

50 50
formed sinks [N]

formed sinks [N]

40 40

30 30

20 20

10 10

0 0
0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.4 0.5 0.6 0.7 0.8 0.9 1.0
time [ ff] time [ ff]

Fig. 6: Histograms of the number of sinks over time. The red histogram gives the formation of sinks passing every test. The green
histogram gives the forced sinks. Vertical dashed lines represent massive star formation. ECF1 is on the left panel and the right
panel is the same initial conditions with foverride set to 1000 referred to as ECF1B.

Second, we completely refactored the point mass time inte- slow-down method) to follow the true dynamics of binaries and
gration methods in Phantom. We replaced the existing second- multiples, which is an essential element of cluster dynamics.
order Leapfrog integrator with an accurate fourth-order Forward In this paper, we performed a first study of massive star clus-
Symplectic Integrator (FSI), which is now the default in the pub- ter formation to benchmark the method. We focused on a single
lic code. Its symplectic properties guarantee perfect conservation fiducial model starting from a massive turbulent molecular cloud
of angular momentum and good conservation of energy. Its for- of 104 M⊙ mass and 10 pc radius. Three simulations with differ-
ward properties also guarantee that there is no issue with dis- ent random initial turbulent conditions were computed with the
sipative forces in the code. In addition to this new integrator, following results:
we implemented specific methods (regularisation integration and
1. We found H ii regions — more precisely those formed first
in the simulation — to be the main game changers for the
Article number, page 15 of 21
 A&A proofs: manuscript no. aa54938-25

time evolution and spatial distribution of star formation in- Aarseth, S. J. 2003, in IAU Symposium, Vol. 208, Astrophysical Supercomput-
side these clouds. Their strength and location can completely ing using Particle Simulations, ed. J. Makino & P. Hut, 295
Aarseth, S. J., Henon, M., & Wielen, R. 1974, A&A, 37, 183
change the future of a star-forming region. We found that the Aarseth, S. J. & Hills, J. G. 1972, A&A, 21, 255
formation of an H ii region can lead either to a very loose, Ahmad, A., González, M., Hennebelle, P., & Commerçon, B. 2024, A&A, 687,
globally expanding cluster in the case of efficient ionisation, A90
or a very compact cluster when H ii regions are off-centred Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17
André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, in Protostars and
and help the gravitational focusing of the cloud. Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning,
2. As a corollary to our first conclusion, since the early for- 27–51
mation of H ii regions is a stochastic process, we found the Armante, M., Gusdorf, A., Louvet, F., et al. 2024, A&A, 686, A122
global evolution of such a cluster to be more stochastic than Baker, H. F. 1905, Proceedings of the London Mathematical Society, s2-3, 24
Bastian, N. & Goodwin, S. P. 2006, MNRAS, 369, L9
when just accounting for stellar dynamics. A suite of simu- Bate, M. R. 1998, ApJ, 508, L95
lations is thus needed to obtain statistically robust outcomes. Bate, M. R. 2009, MNRAS, 392, 1363
3. We found star-formation efficiency (SFE) in our simulations Bate, M. R. & Bonnell, I. A. 2005, MNRAS, 356, 1201
to be slightly high compared to observations, with the mass Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 336, 705
Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362
conversion of gas into stars also not finished at the end of Bate, M. R. & Burkert, A. 1997, MNRAS, 288, 1060
our calculations. However, the effect of our HII feedback on Bate, M. R., Tricco, T. S., & Price, D. J. 2014, MNRAS, 437, 77
the SFE corroborates previous studies. Adding new physical Bisbas, T. G., Haworth, T. J., Williams, R. J. R., et al. 2015, MNRAS, 453, 1324
processes, such as magnetic fields, and feedback processes Bisbas, T. G., Wünsch, R., Whitworth, A. P., & Hubber, D. A. 2009, A&A, 497,
649
including non-ionising radiation, jets, and cloud motions as- Bleuler, A. & Teyssier, R. 2014, MNRAS, 445, 4015
sociated with Galactic shear should help to reduce the SFE. Boss, A. P. & Bodenheimer, P. 1979, ApJ, 234, 289
Bulirsch, R. & Stoer, J. 1966, Numer. Math., 8, 1
These self-consistent simulations of the formation and early Campbell, J. E. 1896, Proceedings of the London Mathematical Society, s1-28,
evolution of embedded clusters improve previous pure N-body 381
Cartwright, A. & Whitworth, A. P. 2004, MNRAS, 348, 589
simulations. The aim of these is to scale up the number of sim- Chin, S. A. & Chen, C. R. 2005, Celestial Mechanics and Dynamical Astronomy,
ulations performed to produce statistically significant results on 91, 301
the dynamical evolution of clusters under the direct influence of Courant, R., Friedrichs, K., & Lewy, H. 1928, Mathematische Annalen, 100, 32
their natal environment. Paper II will aim to characterise the star Cuello, N., Ménard, F., & Price, D. J. 2023, European Physical Journal Plus, 138,
11
formation prescription outcomes, especially stellar mass distri- Dale, J. E., Ercolano, B., & Bonnell, I. A. 2012, MNRAS, 427, 2852
butions and multiplicity properties by exploring its parameter Dehnen, W. & Hernandez, D. M. 2017, MNRAS, 465, 1201
space with a larger set of simulations. This will potentially re- Dehnen, W. & Read, J. I. 2011, European Physical Journal Plus, 126, 55
quire refinements to generate a distribution of primordial bina- Dunham, M. M., Allen, L. E., Evans, Neal J., I., et al. 2015, ApJS, 220, 11
Evans, Neal J., I., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
ries in the actual prescription. On the other hand, another study Evans, II, N. J., Heiderman, A., & Vutisalchavakul, N. 2014, ApJ, 782, 114
will be performed on high-resolution star formation simulations Forest, E. & Ruth, R. D. 1990, Physica D Nonlinear Phenomena, 43, 105
to connect dynamical properties seen below our sink scale to Fujii, M. S., Hattori, K., Wang, L., et al. 2022, MNRAS, 514, 43
our prescription parameters. The new methods presented here Fujii, M. S., Saitoh, T. R., Hirai, Y., & Wang, L. 2021, PASJ, 73, 1074
Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 375
solve performance bottlenecks that have previously limited such González-Samaniego, A. & Vazquez-Semadeni, E. 2020, MNRAS, 499, 668
studies. New optimisations on gas-sink interactions may further Goodwin, S. P. & Whitworth, A. P. 2004, A&A, 413, 929
accelerate the computation. The aim is to carry out around a hun- Gragg, W. B. 1965, SIAM Journal on Numerical Analysis, 2, 384
dred simulations over an integration period of between 10 Myr Grudić, M. 2021, MakeCloud: Turbulent GMC initial conditions for GIZMO,
Astrophysics Source Code Library, record ascl:2106.011
and 20 Myr to have clusters with fully dispersed gas at final Grudić, M. Y., Guszejnov, D., Hopkins, P. F., Offner, S. S. R., & Faucher-
outputs. Ultimately, we planned to study the dynamical evolu- Giguère, C.-A. 2021, MNRAS, 506, 2199
tion of young clusters through statistical analysis by producing a Guszejnov, D., Markey, C., Offner, S. S. R., et al. 2022, MNRAS, 515, 167
large grid of simulations to explore a vast initial parameter space. Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2009, ApJS, 184, 18
Harris, A. & Tricco, T. 2023, The Journal of Open Source Software, 8, 5263
These simulations could then be compared with observations of Hennebelle, P., Commerçon, B., Lee, Y.-N., & Chabrier, G. 2020, ApJ, 904, 194
young clusters, including the latest data available from Gaia. Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3488
Hubber, D. A., Goodwin, S. P., & Whitworth, A. P. 2006, A&A, 450, 881
Acknowledgements. We acknowledge M. Vaileille-Manet for his participation
Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S. 2013, A&A, 554, A17
during the redaction of this paper. We thank the Australian-French Associ-
Könyves, V., André, P., Men’shchikov, A., et al. 2015a, A&A, 584, A91
ation for Research and Innovation (AFRAN) for financial support towards
Könyves, V., André, P., Men’shchikov, A., et al. 2015b, A&A, 584, A91
the 5th Franco-Australian Phantom users workshop held in Melbourne. DJP
Kozai, Y. 1962, AJ, 67, 591
acknowledges Australian Research Council funding via DP220103767 and Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699
DP240103290 and thanks IPAG for hosting and financial support from both Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, in Protostars and Planets
Monash and IPAG during his 2025 sabbatical in Grenoble where this paper was VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 243–266
completed. YB thanks LabEx OSUG@2020 (Investissements d’Avenir – ANR10 Kustaanheimo, P., Schinzel, A., Davenport, H., & Stiefel, E. 1965, 1965, 204,
LABX56) and IDEX for funding a 2-month mobility grant at Monash University. publisher: De Gruyter Section: Journal für die reine und angewandte Mathe-
All the computations presented in this paper were performed using the GRI- matik
CAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57
by Grenoble research communities. FM, FL and IJ acknowledge funding from Larson, R. B. 1969, MNRAS, 145, 271
the European Research Council (ERC) via the ERC Synergy Grant ECOGAL Larson, R. B. 1981, MNRAS, 194, 809
(grant 855130) and from the French Agence Nationale de la Recherche (ANR) Lee, Y.-N. & Hennebelle, P. 2018, A&A, 611, A88
through the project COSMHIC (ANR-20-CE31-0009). We used SPLASH and Lee, Y.-N. & Hennebelle, P. 2019, A&A, 622, A125
Sarracen for many of the figures presented above.(Price 2007; Harris & Tricco Lidov, M. L. 1962, Planet. Space Sci., 9, 719
2023) Liow, K. Y., Rieder, S., Dobbs, C. L., & Jaffa, S. E. 2022, MNRAS, 510, 2657
Liptai, D., Price, D. J., Wurster, J., & Bate, M. R. 2017, MNRAS, 465, 105
Louvet, F., Motte, F., Hennebelle, P., et al. 2014, A&A, 570, A15
Low, C. & Lynden-Bell, D. 1976, MNRAS, 176, 367
References Lucy, L. B. 1977, AJ, 82, 1013
Machida, M. N. & Matsumoto, T. 2012, MNRAS, 421, 588
Aarseth, S. J. 1999, PASP, 111, 1333 Makino, J. & Aarseth, S. J. 1992, PASJ, 44, 141

Article number, page 16 of 21

 Y. Bernard et al.: DAWN

Mauxion, J., Lesur, G., & Maret, S. 2024, A&A, 686, A253
Mikkola, S. & Aarseth, S. 2002, Celestial Mechanics and Dynamical Astronomy,
84, 343
Mikkola, S. & Aarseth, S. J. 1996, Celestial Mechanics and Dynamical Astron-
omy, 64, 197
Mikkola, S. & Tanikawa, K. 1999, MNRAS, 310, 745
Motte, F., Nony, T., Louvet, F., et al. 2018, Nature Astronomy, 2, 478
Motte, F., Schilke, P., & Lis, D. C. 2003, ApJ, 582, 277
Mukherjee, D., Zhu, Q., Trac, H., & Rodriguez, C. L. 2021, ApJ, 916, 9
Myers, A. T., Klein, R. I., Krumholz, M. R., & McKee, C. F. 2014, MNRAS,
439, 3420
Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in Astronomical Society
of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed.
S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275
Omelyan, I. P. 2006, Phys. Rev. E, 74, 036703
Omelyan, I. P., Mryglod, I. M., & Folk, R. 2002, Computer Physics Communi-
cations, 146, 188
Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and
active galactic nuclei
Parker, R. J., Wright, N. J., Goodwin, S. P., & Meyer, M. R. 2014, MNRAS, 438,
620
Peretto, N., Fuller, G. A., Duarte-Cabral, A., et al. 2013, A&A, 555, A112
Preto, M. & Tremaine, S. 1999, AJ, 118, 2532
Price, D. J. 2007, PASA, 24, 159
Price, D. J. 2012, Journal of Computational Physics, 231, 759
Price, D. J. & Bate, M. R. 2009, MNRAS, 398, 33
Price, D. J., Wurster, J., Tricco, T. S., et al. 2018, PASA, 35, e031
Rantala, A., Naab, T., Rizzuto, F. P., et al. 2023, MNRAS, 522, 5180
Rantala, A., Naab, T., & Springel, V. 2021, MNRAS, 502, 5546
Renaud, F., Agertz, O., & Gieles, M. 2017, MNRAS, 465, 3622
Sadavoy, S. I. & Stahler, S. W. 2017, MNRAS, 469, 3881
Schneider, N., Motte, F., Bontemps, S., et al. 2010, A&A, 518, L83
Sheng, Q. 1989, IMA Journal of Numerical Analysis, 9, 199
Springel, V., Pakmor, R., Zier, O., & Reinecke, M. 2021, MNRAS, 506, 2871
Suzuki, M. 1991, Journal of Mathematical Physics, 32, 400
Tuckerman, M., Berne, B. J., & Martyna, G. J. 1992, J. Chem. Phys., 97, 1990
Verliat, A., Hennebelle, P., González, M., Lee, Y.-N., & Geen, S. 2022, A&A,
663, A6
Wall, J. E., McMillan, S. L. W., Mac Low, M.-M., Klessen, R. S., & Portegies
Zwart, S. 2019, ApJ, 887, 62
Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020a, MNRAS, 497, 536
Wang, L., Nitadori, K., & Makino, J. 2020b, MNRAS, 493, 3398
Wright, N. J., Parker, R. J., Goodwin, S. P., & Drake, J. J. 2014, MNRAS, 438,
639
Wurster, J., Bate, M. R., & Price, D. J. 2018, MNRAS, 475, 1859
Wurster, J. & Bonnell, I. A. 2023, MNRAS, 522, 891
Yoshida, H. 1990, Physics Letters A, 150, 262

Article number, page 17 of 21

 A&A proofs: manuscript no. aa54938-25

Appendix A: Details of N-body methods where all e prefactors are called error coefficients and are func-
tions of (cki , dki ). The first two coefficients correspond to the sum
Appendix A.1: FSI derivation and accuracy of cki and dki , which by definition needs to equal unity. Thus, ec
The original Leapfrog algorithm used in Phantom can be easily and ed are equal to unity. One can refactor equation (A.7) using
derived using the Hamiltonian formalism to write the equation these equalities according to
of motion of the studied system. This equation can be written as N
Ysub

δt δt

dw log   e cki T
e dki U  = δt [H + Herr (δt)] = δtHneigh . (A.8)
= {w, H} (A.1) 
dt i=1
where w is the phase space vector of the system, H is the system This new form hints at what is needed for a symplectic inte-
Hamiltonian and {, } is the Poisson brackets operator. The formal gration. What is evolved in such a scheme is not the true Hamil-
solution of this equation over a certain time interval δt is then tonian H but a neighbour of this one Hneigh which differs from an
straightforwardly given by error called Herr . The goal of constructing high-order symplec-
w(t + δt) = eδt H w(t), (A.2) tic integrators is then to minimise this error quantity by setting
δtH the next error coefficients to zero. However, it has been proven
where e is an operator that evolves the system in time. All (Sheng 1989; Suzuki 1991) that after second order (cki , dki ) co-
symplectic integrators come from a particular way to approxi- efficients must contain some negative values to set the next order
mate this time evolution operator. Indeed, it can be reformulated coefficient to zero. However, negative time steps are problem-
into a product decomposition with a specific Hamiltonian that atic when integrating dissipative systems that are common in as-
can be split into two parts, for example, a kinetic part T and a trophysics. The solution proposed by Suzuki (1991) was to add
potential part U. It gives one of the first terms of the error Hamiltonian produced by a
Nsub
Y Leapfrog scheme into the true one. One can then erase all er-
eδt H = eδt (T +U) ≈ eδt cki T eδt dki U , (A.3) ror terms of the fourth order. Nonetheless, the part added to the
i=1 classical Hamiltonian still needs to be computed. It is possible to
where Nsub corresponds to the number of subdivisions of the time write the new Hamiltonian as
evolution operator. cki and dki are coefficients (also called sym- 1
plectic coefficients) that weigh the decomposition into multiple H = T + U + δt2 [U, [T, U]]. (A.9)
48
parts. The sums of these two series of coefficients need to be
equal to 1 to respect the decomposition. The resulting two op- In this equation, [U, [T, U]] is a composition of two commuta-
erator types can be associated with different steps in the time tors between U and T . Dehnen & Hernandez (2017) showed that
evolution. U potential is always associated with a force applied this commutator corresponds to an extra potential quantity in the
on the system using eδtU . eδtU is only applied on momentum co- Hamiltonian defined as
ordinates X N
1 ∂U
!2 XN

eδt U p(t) = p(t + δt). [U, [T, U]] = − =− mi ∥ai ∥. (A.10)

(A.4)
i
mi ∂ri i
This operation is called a kick. eδtT represent the inertial part of
the evolution. It will then only modify positions depending on Figure A.1a and A.1b illustrate the superior accuracy of the
velocities FSI method by showing the trajectory of one point mass around a
fixed gravitational potential during 10 orbits and its energy con-
eδt T x(t) = x(t + δt). (A.5) servation (logarithm of the relative energy) for one orbit, respec-
This operation is called a drift. A system can be integrated tively. FSI1 and FSI2 correspond, respectively, to the original
into a given time following a succession of alternating kick gradient force calculation and the extrapolation method. The ini-
and drift operations using this decomposition method. The num- tial eccentricity of the body is large (e = 0.7) to show the orbital
ber of kick and drift operations Nsub and a couple of coeffi- precession that can be introduced by numerical error. A large
cients (cki , dki ) give the integrator order of truncation. for the timestep is also deliberately chosen (around 60 steps per orbit)
Leapfrog, Nsub = 2 and (ck1 , dk1 ) = (1, 1/2) or (1/2, 1) and to amplify numerical imprecision. As expected, energy conser-
(ck2 , dk2 ) = (0, 1/2) or (1/2, 0) depending on the order of the vation and orbit precession are much improved with a fourth-
kick and drift operations. The number of operations and the sym- order scheme than with a simple leapfrog scheme. It also con-
plectic coefficients can be efficiently found (Yoshida 1990) with firms the Chin & Chen (2005) result that forward schemes have
the Campbell-Baker-Hausdorff (CBH) formula (Campbell 1896; better conservation properties.
Baker 1905), namely
1 Appendix A.2: Few body problem regularised by TTL
log(eX eY ) = X + Y + [X, Y]
2 (A.6) Let’s introduce a general N-body (only Newtonian) system that
1
+ ([X, [X, Y]] + [Y, [Y, X]]) + O(X, Y, [X, Y]). we want to regularise with the Time Transformed Leapfrog
12 (TTL) scheme. The equations of motion are given by
Substituting the left hand of equation (A.3) into equation (A.6),
N X
N
one finds X rj − rj
N sub r̈i = G mj = Fi (A.11)
Y δt ri3j
log( eδt cki T eδt dki U ) = δt(ec T + ed U + ecd [T, U] i=1 i, j
2
i=1 where N is the number of bodies and G the gravitational constant
δt and ri j is the distance between body i and j. For convenience
+ (eccd [T, [T, U]] + edcd [U, [T, U]]) G = 1. Now, let’s introduce the time transformation
12
+ O(T, U, [T, U]), (A.7) ds = Ω(r)dt, (A.12)
Article number, page 18 of 21
 Y. Bernard et al.: DAWN

(a) 10 orbits of a point mass around a fixed potential. Initial eccentric-

ity is set to 0.7 and the semi-major axis is 1. The amount of (spurious)
orbital precession can be seen to inversely correlate with integrator pre-
cision. (b) Logarithm of the relative energy of the particle on 1 orbit.

Fig. A.1: Evolution of a point mass around a fixed potential

where Ω(r) is an arbitrary time transformation function. To fol- where ck is the drift coefficient and δt is the time step. Then
low what is done in Mikkola & Aarseth (2002) and Wang et al.
dk δt Fi (ri∗ )
(2020b), this function can be defined as vi∗/2 = (A.21)
2 Ω(ri∗ )
N X
N N
X mi m j vi∗/2 ∂Ω(ri∗ )
Ω(r) = .
X
(A.13) W∗ = dk δt · (A.22)
i=1 i, j
ri j
i=1
Ω(r i∗ ) ∂ri∗

One can now set a new auxiliary variable, often called W = Ω, dk δt Fi (ri∗ )
vi∗ = (A.23)
to rewrite the system according to 2 Ω(ri∗ )
dri vi where dk is the kick coefficient and ∗ represents updated vari-
= (A.14) ables. As W depends on the velocities of each body, it is needed
ds W to update both variables simultaneously. To achieve that, the ve-
dvi Fi
= (A.15) locity kick is split in half with a W kick between the two parts.
ds Ω These two operations can be successively applied to the system
dt 1 to integrate its motion using any kind of symplectic integrator.
= ; (A.16)
ds W
where index i indicates one body of the system. The addition of Appendix A.3: Discussion on slow-down factor
a new variable gives a supplementary variable in the system that
can be extracted from the differential equations. As has been shown above, κ quantifies the number of original
periods where the perturbation applied to the binary can be av-
N
dW X ∂ri ∂Ω eraged. On a highly perturbed binary, it is impossible to average
= · , (A.17) the perturbation on multiple periods, which means κ = 1. On the
dt i=1
∂t ∂ri other side, a hard binary can stay unperturbed during a whole
simulation. The averaging is not risky at all, and then κ can stay
which can be expressed using our new “time” coordinate at a high value. That being said, one can see that this factor needs
N to scale with the strength of the perturbation. The solution pro-
dW X 1 ∂ri ∂Ω
= · . (A.18) posed originally by Mikkola & Aarseth (1996) is to compute the
ds i=1
Ω ∂t ∂ri relative tidal perturbation at the maximum apocenter ( 2a with
e = 1), giving
This new variable W acts like an associated momentum for the
time, which means that it will have the same behaviour as a ve- 8(m p + m s )a3 N
X mj
locity in a symplectic scheme. With that in mind, it is then pos- γ= , (A.24)
m p ms r3
j=1, j,p,s b j
sible to form one drift and one kick operation from this system,
given respectively by where a is the semi major axis of the binary,p and s point on the
vi binary components. rb j is the distance between any perturbers
ri∗ = ri + ck δt (A.19) and the binary centre of mass. κ value is then equal to
W
ck δt γ0
t∗ = t+ , (A.20) κ= , (A.25)
w γ
Article number, page 19 of 21
 A&A proofs: manuscript no. aa54938-25

where γ0 is a safety factor to adjust the precision wanted by the

STARBENCH solution
user, often chosen between 105 and 106 . Wang et al. (2020b) 2.25 Phantom derived
chose something slightly different for this quantity. In their work,
they proposed 2.00

m p ms N
X rb3 j 1.75
κ= . (A.26)
m p + m s (a(1 + e)) mj

r[pc]
j=1, j,p,s 1.50

Instead of using the maximum possible apocenter of the binary, 1.25

they chose the actual apocenter. The most different part of their
proposition is the summation. They directly sum the cubed dis- 1.00
tances scaled by the mass of each perturber. Even if these two ex- 0.75
pressions look similar, they have totally different behaviour. For
example, let’s consider a system composed of one hard binary 0.50
that can be slowed down and two single perturbers. To simplify
the system, each bodies have the same mass. One perturber is far 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
from the binary, producing weak perturbations. The other one is time[Myr]
on a hyperbolic trajectory with the binary close to its periapsis.
It then produces a strong perturbation. In theory, this scenario Fig. B.1: Time evolution of the HII region front radius. The green
should produce a slowdown factor near unity due to the close line is the analytical solution found in Bisbas et al. (2015). The
encounter with perturber two. The first expression will first sum red line is the Phantom-derived solution. The sudden drop after
the inverse cube distance, which will give a large number (since 0.6 Myr is due to the front shock reflection onto the boundary
the strong perturber dominates the summation). Taking the in- particles.
verse of this summation gives a number less than unity (which
means no slowdown in this framework). The first expression has

log density [g/cm3]

good behaviour, but with the second expression, the dominant
t=0.3 Myr
perturber is the weaker one, which will give a high κ value. This -20
expression is then potentially wrong. Therefore, in this work, we 2
stick with the original expression, which leads to more robust
behaviour by changing only the apocenter expression.
Finally, the κ value of each binary in the system needs to be -21
updated regarding the perturbation modification in the system.
y [pc]

However, this evolution needs to be slow to conserve the secular 0

evolution properties. Wang et al. (2020b) proposed a limit value
for κ that is controlled by the changing rate of the perturbation. -22
Formally, it is written as
c∥⟨rb j ⟩m ∥ -2
κmax = (A.27) -23
P∥⟨vb j ⟩m ∥
where c = 0.1 is a safety factor, P is the period of the binary, vb j
is the relative velocity vector between perturbers and the binary -2 0 2
x [pc]
CoM. ⟨⟩m refers to a mass-weighted average.
Fig. B.2: Density slice of a uniform density sphere where a HII
Appendix B: HII region unit test region expands from its centre. A central stellar object produces
many ionising photons that heat the gas around it. Two media
We performed a simple test in order to validate the HII feed- (hot and cold) are separated by a front shock that propagates
back implemented into this work, following what was proposed inside the cold sphere, sweeping up the gas. A denser shell of
in Bisbas et al. (2015). They described two main setups to check material forms just next to the ionisation front.
the validity of HII region expansion at two different time scales
(early and late) compared to analytical solutions. Such solutions
exist only in very specific and simple systems. The selected test front reflecting from boundary particles which compresses the
(late) tries to expand an HII region produced by a 40M⊙ star in HII bubble, forcing it to shrink suddenly. This phenomenon is
the centre of a homogeneous sphere of hydrogen atomic gas of not present in Bisbas et al. (2015) work. Nothing is detailed on
2.9pc in radius. Here, we set the cloud mass to 8 × 10M⊙ , which boundary conditions for SPH runs. Nonetheless, the satisfactory
was divided into 8 million SPH particles. The neutral tempera- agreement on the first part is sufficient to validate the implemen-
ture was set to 1000 K. Boundary particles are also mandatory at tation, as boundary conditions in large clouds have already been
the edge of the cloud to avoid relaxing into empty space, which tackled in 3.2.
could lead to an artificial acceleration of the HII bubble. Fig-
ure B.1 plots the HII region radius derived by Phantom, while
Figure B.2 shows a density slice of the cloud. Our implemen- Appendix C: SDAR unit test
tation agrees with the Bisbas et al. (2015) analytical solution
during 0.6Myr and comes close to the stagnation radius equal We performed a unit test taken from Wang et al. (2020b) to val-
to 2.31 pc. The sudden drop after 0.6Myr is due to the shock idate our SDAR implementation. This test consists of a hierar-
Article number, page 20 of 21
 Y. Bernard et al.: DAWN

40

kappa
0.95
e

20

0.9
0
0 50 100 150 0 50 100 150
t t

Fig. C.1: Time evolution of the inner binary eccentricity. This Fig. C.2: Time evolution of the slowdown factor κ. This one os-
eccentricity oscillates under the perturbation of the single star cillates on a shorter timescale that corresponds to the outer bi-
orbiting the binary in a very eccentric orbit. nary period. When the single star reaches the apoapsis of its or-
bit. κ reaches its maximum at this moment. Oppositely and as
expected, κ decreases during its descent to the periapsis. At this
chical triple system (binary-single) with a hard eccentric inner moment, slowing down the inner binary motion is not possible,
binary orbiting another single star. The outer orbit is eccentric as therefore κ = 1
well, with a semi-major axis ratio of 1000 between the inner and
outer binary. All orbital parameters are summarised in Table 1 of
Wang et al. (2020b).
Kozai-Lidov oscillations (Kozai 1962; Lidov 1962)) are ex-
pected on the inner eccentricity and inclination in such systems
as the outer star will perturb the inner binary at each periapsis.
Such a system is a good test to verify if our SDAR module can
correctly model this oscillation with and without the slowdown
method. Figure C.1 shows the eccentricity of the inner binary
over time. As expected, A clear oscillation between 0.9 and 1.
is visible, and both methods converge on the same solution. Fig-
ure C.2 shows the evolution of κ slow-down factor over time.
This one oscillates between around 50 and 1 depending on the
outer binary orbital phase. If the single outer star is at apoap-
sis, its perturbation on the inner binary is minimal; therefore, κ
is high. Its value decreases accordingly down to 1 as the star
reaches the periapsis of the outer orbit. In this configuration, the
perturbation is maximal. Hence, slowing down the inner binary
is not possible. These results also converge with what was found
in Wang et al. (2020b).

Article number, page 21 of 21