© The Authors 2025

1 Introduction

Protostars harbor a rich and complex chemistry, both in their envelope and their protoplanetary disk (e.g., Jørgensen et al. 2016; McGuire 2022). Protostars are however chemically diverse (Lefloch et al. 2018), highlighting the importance of understanding the chemical and physical processes that drive this diversity.

Interferometers such as the Atacama Large Millimeter/ submillimeter Array (ALMA) and the NOrthern Extended Millimeter Array (NOEMA), with their high angular resolution (<1″), can probe the chemical content of protostellar sources at the scale of protoplanetary disks. Coupled with spectrometers, those instruments deliver observations in the form of hyperspectral cubes. Although rigorous analysis may distinguish the origin of molecular emission (the disk or the outflow for example, see e.g., Tychoniec et al. 2021), interferometric observations remain complex to analyze. Chemical models are a complementary tool to understand the origin of molecules in such regions. By carrying out radiative transfer calculations on simulated chemical abundance maps, we can produce synthetic (i.e., simulated) observations (as in Nazari et al. 2022), which allows for a more direct comparison between models and interferometric observations.

Such models require physical conditions, typically density and temperature, as input for their calculations. They are usually provided by analytical physical models (e.g., Visser et al. 2009; Gavino et al. 2021; Murillo et al. 2022) or by magnetohydrodynamics (MHD) simulations (e.g., Coutens et al. 2020; Navarro-Almaida et al. 2024). MHD simulations however require large amounts of resources. Simulating the collapse of a prestellar core, followed by the formation and evolution of a central object and a protoplanetary disk, is extremely challenging due to the wide range of spatial and temporal scales needed to describe those phenomena (Vaytet et al. 2018; Ahmad et al. 2024). Analytical and semi-analytical models require much less computational resources, allowing for broader accessibility and making parametric studies affordable. These models rely on simplifying assumptions and cannot fully capture the complexity of the underlying physics. On the other hand, they provide clear relationships between key physical parameters and their macroscopic results (e.g., the impact of the magnetic mass-to-flux ratio on the size of protoplanetary disks), which greatly facilitates the interpretation.

In this paper, we present the Analytical Protostellar Environment (APE) code¹. APE simulates density, temperature, extinction, and dust-size distribution maps in a protostellar environment, including the central object, the envelope, the disk, and the outflow. The analytical model can simulate the whole system evolution throughout the prestellar core, Class 0 and Class I phases, at a fraction of the cost of MHD simulations. In addition to the physical model, we provide scripts and interfaces to other publicly available codes to perform chemical simulations (Nautilus, Ruaud et al. 2016), radiative transfer calculations (RADMC-3D, Dullemond et al. 2012), and synthetic imaging (Imager²). The aim of APE is to provide a fast and accessible way to generate maps of chemical abundances and their simulated emission.

Section 2 is dedicated to the description of APE, with Sect. 2.1 explaining in more detail the general purpose of the code. In Sect. 3, we show an example of application. We conclude in Sect. 4.

Fig. 1 Schematic of an APE snapshot simulation.

2 The APE code

2.1 General description

APE can be used in two modes: a “snapshot” mode, in which it generates a map of density and temperature, and a “particle” mode, which calculates the trajectory of a virtual particle in the evolving system. The initial condition (t = 0) is a critical Bonnor-Ebert sphere (Ebert 1955; Bonnor 1956) for which we derived analytical evolution equations for its dynamics and density profile. After a free fall time, a central object, a disk and an outflow are created. The mass accreted in the center is distributed between the central object and the disk, from which the code computes the radius of the central object, luminosity and temperature, and the disk density and temperature profiles. Dust grains are included in the form of an evolving size-distribution. A schematic of the simulated region is displayed in Fig. 1. All APE input parameters are listed in Table 1 with their description.

2.1.1 Working modes

Snapshot mode. In snapshot mode, APE generates a map of a protostellar system with an envelope, a disk, a central object and an outflow. The code uses a 2D grid in spherical coordinates, assuming a cylindrical symmetry and a symmetry with respect to the mid-plane of the disk. The region simulated, in spherical coordinates, is then [0 ≤ r ≤ l_box, φ= 0, 0 ≤ θ ≤ π/2], with l_box the size of the simulated region (input parameter). The user can either choose a linear or a logarithmic grid spacing in radius. In each cell, the flow variables (density, temperature, velocities and the grain size-distribution) are calculated as a function of the model parameters input by the user: the mass of the initial Bonnor-Ebert sphere, the magnetic mass-to-flux ratio and the age of the central object. An example of a snapshot is displayed in Fig. 2. Generating the full map takes typically 1 second for 10 000 grid cells on a regular laptop.

Particle mode. The purpose of the particle mode is to emulate hydrodynamics simulations at a fraction of the cost. The user chooses an initial position of the particle at t = 0 (in the initial Bonnor–Ebert sphere). The trajectory of the particle is then calculated over time, until the particle reaches the center (the central object if it is created), the outflow or until a time defined by the user. The result is a complete density and temperature history of the particle. This computation can be done in reverse, with the user choosing the age of the system, like in snapshot mode, and the final position of the particle at this time. The trajectory is then calculated backward until the system reaches its initial condition at t = 0. Five to thirty seconds are needed to compute the full trajectory.

Grid of particles. The snapshot and particle modes can be combined to compute the full history of the gas displayed in a snapshot. In that case, a map is created like in snapshot mode, and a particle placed in each grid cell. APE then computes the trajectory in reverse of all those particles until t = 0. The density and temperature history of those particles can then be used as an input for a chemical modeling code, to ultimately produce maps of chemical abundances.

Fig. 2 Top panel: density map generated by APE. Bottom panel: zoom of the inner 200 au, with a reduced color sampling.

2.1.2 Radiative transfer

APE has been designed to work with the radiative transfer code RADMC-3D (Dullemond et al. 2012). In snapshot mode, the output files of APE can be directly used as input for RADMC-3D with opacities calculated directly from the dust grain model of APE. That way, the user can either generate a temperature map more accurate than the analytical descriptions provided in APE, or address any science topic for which radiative transfer calculations are needed. APE also offers the possibility to use temperatures from RADMC-3D simulations in particle mode. Since running radiative transfer calculations at each time-step would be costly, we provide scripts to generate snapshots at regular time intervals, on which RADMC-3D is run to produce temperature maps at those times. Subsequent APE simulations in particle mode then use those maps to interpolate (in time and space) the temperature at every time-step in the particle trajectory.

2.2 The envelope

The envelope initial condition is a critical Bonnor-Ebert sphere (Ebert 1955; Bonnor 1956) in free fall.

2.2.1 Base of the model

The Bonnor-Ebert sphere is a spherically symmetric isothermal cloud of gas in which the gravity and the thermal pressure are at equilibrium. The density profile, ρ, is normalized by the central density, ρ_c, $$D=\frac{\rho }{{\rho }_{\mathrm{c}}},$$(1)

while the spherical radius coordinate, r, is normalized as $$x_{\mathrm{b}\mathrm{e}}=r\frac{\sqrt{4\pi G{\rho }_{\mathrm{c}}}}{c_{\mathrm{s}}},$$(2)

with G the gravitational constant and c_s the sound speed. The equilibrium between thermal pressure, P, and gravity potential, Φ, reads $$\mathrm{\nabla }\mathrm{\Phi }=-\frac{1}{\rho }\mathrm{\nabla }P.$$(3)

We can combine Eq. (3) with the equation of state for an isothermal perfect gas $$P=\rho c_{\mathrm{s}}^2,$$(4)

to obtain the relation $$\rho ={\rho }_{\mathrm{c}}{e}^{-\frac{\mathrm{\Phi }}{c_{\mathrm{s}}^2}}.$$(5)

Finally, combining the Poisson equation $$\mathrm{\Delta }\mathrm{\Phi }=4\pi G\rho ,$$(6)

with the divergence of Eq. (3) in spherical coordinates results in the Lane–Emden equation $$\frac{{\mathrm{d}}^2\mathrm{\Phi }}{\mathrm{d}{x}_{\mathrm{b}\mathrm{e}}^2}+\frac{2}{x_{\mathrm{b}\mathrm{e}}}\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}{x}_{\mathrm{b}\mathrm{e}}}=e^{-\frac{\mathrm{\Phi }}{c_{\mathrm{s}}^2}}.$$(7)

The normalized radius of a critical Bonnor-Ebert sphere is x_be = x_crit ≈ 6.451, above which no hydrostatic equilibrium is possible. The gravitational potential and density profiles, Φ(x_be) and ρ(x_be), are integrated using the Lane-Emden Equation (7), with Φ(0) = 0 and dΦ/dx_be(0) = 0. The resulting Bonnor-Ebert density profile, D(x_be), is represented in Fig. 3.

The mass of the cloud of radius r_cld is given by $$\begin{array}{rl}{M}_{\mathrm{c}\mathrm{l}\mathrm{d}}& ={\int }_0^{r_{\mathrm{c}\mathrm{l}\mathrm{d}}}4\pi r^2\rho (r)dr,\\ & =4\pi {\rho }_{\mathrm{c}}{\left(\frac{c_{\mathrm{s}}}{\sqrt{4\pi G{\rho }_{\mathrm{c}}}}\right)}^3{\int }_0^{x_{\text{crit\hspace{0.17em}}}}D(x)x_{\mathrm{b}\mathrm{e}}^{2}d{x}_{\mathrm{b}\mathrm{e}},\\ & =I_m{c}_{\mathrm{s}}^3{G}^{-\frac{3}{2}}(4\pi {)}^{-\frac{1}{2}}{\rho }_{\mathrm{c}}^{-\frac{1}{2}},\end{array}$$(8)

where I_m is the value of the integral on x that approximates to 15.7 when integrating the whole profile. Finally the central density is given by $${\rho }_{\mathrm{c}}=I_m^2{c}_{\mathrm{s}}^6{G}^{-3}(4\pi {)}^{-1}{M}_{\mathrm{c}\mathrm{l}\mathrm{d}}^{-2}=7.78\times {10}^{-19}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}{\left(\frac{M_{\mathrm{c}\mathrm{l}\mathrm{d}}}{M_{\odot }}\right)}^{-2}.$$(9)

We assumed that the enclosed mass is invariant with time for a given fluid particle (equivalent to collapsing spherical shells that cannot cross each-other). The free fall velocity of a particle at radius r would then be $v_r(r)=-\sqrt{2G{M}_{\text{in}}(r)/r}$, where M_in(r) is the mass enclosed at radius r. However, several physical processes support the gas against gravity, effectively reducing the infall velocity: rotation, thermal pressure, magnetic fields and turbulence. The effective free fall time can then be increased by a factor ~1.3–2 with rotation and thermal pressure alone (see e.g., Gabbasov et al. 2017), and by a factor ≳2 between moderate and strong magnetic fields (Masson et al. 2016). We therefore assumed a factor 2 reduction in radial velocity and adopt $$v_r(r)=-\sqrt{\frac{G{M}_{\mathrm{i}\mathrm{n}}(r)}{2r}}.$$(10)

Solving the differential equation $$\frac{dr}{dt}=v_r(r),$$(11)

yields the following analytical expression $$r(t)={[r_0^{\frac{3}{2}}-\frac{3}{2}\sqrt{\frac{G{M}_{\mathrm{i}\mathrm{n}}\left(r_0\right)}{2}}t]}^{\frac{2}{3}},$$(12)

where r₀ is the radius of the fluid particle at t = 0.

Fig. 3 Normalized density profile of a Bonnor-Ebert sphere.

2.2.2 Computing the flow variables

The aim is to compute the density and velocity at any given time, t, and radius, r, in the envelope. Rewriting Eq. (12) at a time t and radius r, the equation $$\frac{1}{r_0}{[r_0^{\frac{3}{2}}-\frac{3}{2}\sqrt{\frac{G{M}_{\mathrm{i}\mathrm{n}}\left(r_0\right)}{2}}t]}^{\frac{2}{3}}-\frac{r}{r_0}=0,$$(13)

is solved for r₀(r, t) by dichotomy. With the initial position of the fluid particle, we obtain the enclosed mass by integrating Eq. (8), and the radial velocity with Eq. (10). To calculate the density, we considered the compression of a thin spherical shell contained between radii r and r + dr. The density is then computed using the initial density of the fluid particle ρ(r₀(r, t), t = 0), from the initial Bonnor-Ebert density profile, and the volume compression of the spherical shell between [r₀(r, t); r₀(r + dr, t)] and [r, r + dr]. We therefore have $$\rho (r,t)=\rho (r_0(r,t),t=0)\frac{r_0(r+\mathrm{d}r,t{)}^3-r_0(r,t{)}^3}{(r+\mathrm{d}r{)}^3-r^3}.$$(14)

In the code, we take dr = r/100. While the Bonnor–Ebert density profile is asymptotically proportional to r⁻², the slope transitions toward r^−1.5 as the collapse proceeds, similarly to other free fall analytical models (e.g., Cassen & Moosman 1981; Hartmann 2009). Figure 4 shows the density profile of the envelope for a 2 M_⊙ cloud at time t ≈ 192 kyr (free fall time of 152 kyr). To validate our envelope model, we also present a comparison between the envelope model of APE and a MHD simulation in Appendix A.

The envelope rotation velocity has no impact on the physical model, but can be relevant for synthetic observations. We assumed that the envelope is initially rotating at an angular velocity of Ω₀, which is chosen by the user. The angular velocity at a radius r and time t is then simply calculated assuming angular momentum conservation $$\mathrm{\Omega }(r,t)={\mathrm{\Omega }}_0{\left(\frac{r_0(r,t)}{r}\right)}^2.$$(15)

Fig. 4 Density profile of a collapsing Bonnor–Ebert sphere following our model, for a 2 M⊙ sphere at time t ≈ 192 kyr (in blue). The dashed black line represents a $r^{-\frac{3}{2}}$ profile.

2.2.3 Maximum time

Assuming the cloud has a normalized radius of x = x_crit = 6.451, we can compute the cloud radius r_cld with Eq. (2). The edge of the cloud reaches the center at time t_max, with $$r\left(t_{max}\right)={[r_{\mathrm{c}\mathrm{l}\mathrm{d}}^{\frac{3}{2}}-\frac{3}{2}\sqrt{\frac{G{M}_{\mathrm{c}\mathrm{l}\mathrm{d}}}{2}}{t}_{max}]}^{\frac{2}{3}}=0.$$(16)

Therefore $$t_{max}=\frac{2\sqrt{2}}{3\sqrt{G{M}_{\mathrm{c}\mathrm{l}\mathrm{d}}}}{r}_{\mathrm{c}\mathrm{l}\mathrm{d}}^{\frac{3}{2}}.$$(17)

Replacing r_cld with Equation (2), M_cld with Eq. (8), and ρ_c using the free fall time expression $$t_{\mathrm{f}\mathrm{f}}=\sqrt{\frac{3\pi }{32G{\rho }_{\mathrm{c}}}},$$(18)

we obtain $$t_{max}=\frac{8x_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}^{\frac{3}{2}}}{3\sqrt{3}\pi \sqrt{I_m}}{t}_{\mathrm{f}\mathrm{f}}\approx 2t_{\mathrm{f}\mathrm{f}}.$$(19)

Beyond 2 free fall time, we therefore considered by default that the envelope is depleted. Should the user wish to extend the envelope lifetime beyond 2 free fall times, APE provides the possibility to do so by changing an “effective” x_crit which only affects the extent of the initial Bonnor-Ebert sphere (not the central density or the free fall time).

2.2.4 Envelope temperature

In APE, the temperature may be computed using the radiative transfer code RADMC-3D (Dullemond et al. 2012), which provides a temperature at each location for each grain size. We also provide a size-independent approximation of the temperature for faster calculations, which we describe hereafter. The initial envelope temperature T_mc is set by the user with the temp_molecular_cloud parameter. Once the central object is ignited and starts radiating from its surface, the temperature rises. Assuming the envelope is optically thin, the total energy flux is conserved when integrated over a sphere. Therefore, writing ϕ(r) the radiative power received per surface unit at radius r, the quantity ϕ(r)4πr² is constant with respect to r. Assuming the central object is a black body, the radiative flux emitted per surface unit at the central object surface is $${\varphi }_{\ast }\left(R_{\ast }\right)={\sigma }_{\mathrm{B}}{T}_{\ast }^4,$$(20)

with R_* and T_* the radius and temperature of the central object, and σ_B the Stefan-Boltzmann constant. The total radiative flux crossing a sphere of radius r is then equal to the flux emitted by the central object $$4\pi r^2\varphi (r)=4\pi R_{\ast }^2{\varphi }_{\ast }\left(R_{\ast }\right).$$(21)

Assuming a background radiation from the surrounding molecular cloud at temperature T_mc, a dust grain at a distance r then receives the radiative power per surface unit $${\varphi }_{\text{grain\hspace{0.17em}}}(r)=\varphi (r)+{\sigma }_{\mathrm{B}}{T}_{\mathrm{m}\mathrm{c}}^4,$$(22)

At thermal equilibrium, assuming that the grains are black bodies, the received and emitted flux are equal. $${\varphi }_{\text{grain\hspace{0.17em}}}(r)={\sigma }_{\mathrm{B}}{T}_{\text{grain\hspace{0.17em}}}(r{)}^4.$$(23)

Combining Eqs. (20), (22) and (23), we obtain the grain equilibrium temperature independent of the grain size. $$T_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}(r)={[{\left(\frac{R_{\ast }}{r}\right)}^2{T}_{\ast }^4+T_{\mathrm{m}\mathrm{c}}^4]}^{\frac{1}{4}}.$$(24)

Assuming the gas and grains reach thermal equilibrium quickly, we write T(r) ≈ T_grain(r). This model does not account for the shadow cast by the disk, which lowers the temperature of the envelope near the midplane. However, this temperature is simply an estimate that does not affect the density or dynamics of the envelope. This relation provides a good approximation to the envelope temperature calculated by the radiative transfer computations of RADMC-3D, for r ≳ 50 au. Closer to the central object, the envelope temperature is higher along the walls of the outflow cavity. For a precise temperature map, we encourage users to perform the radiative transfer calculations with the tools provided with APE.

2.3 The central object

The typical density at which the first hydrostatic core collapses to form the central object is typically 10⁻⁸ g cm⁻³ (Larson 1969; Vaytet et al. 2018). This density is reached at around t = t_ff after the start of the collapse. We therefore set the creation of the central object and the disk at t = t_ff.

2.3.1 Mass and radius

We defined the total accreted mass M_acc as the mass of the envelope that has collapsed within the innermost cells (typically on a scale of 1 au). This mass is distributed between the disk M_disk and the central object M_*. Adams & Shu (1986) derive $$M_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}=(1-\eta ){\mathcal{M}}_{\mathrm{d}}{M}_{\mathrm{a}\mathrm{c}\mathrm{c}},$$(25)

with η a free parameter, and $${\mathcal{M}}_{\mathrm{d}}=\frac{1}{3}{u}_{\ast }{\int }_{u_{\ast }}^1(1-u{)}^{\frac{1}{2}}{u}^{-\frac{4}{3}}du,$$(26)

where u_* = R_*/r_disk. Following Young & Evans (2005), we take η = 0.75 so that the ratio M_disk/M_* remains close to the observationally estimated value of ≈1/3 (Li 2002; Kratter et al. 2008). The star mass is therefore $$M_{\ast }=M-M_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}=(1-(1-\eta ){\mathcal{M}}_{\mathrm{d}})M_{\mathrm{a}\mathrm{c}\mathrm{c}}.$$(27)

The radius of the central object, R_*, is controlled by the increase in entropy at the accretion shock on its surface, as well as the extent of its convection zone and the deuterium-burning zone. This calculation was performed in detail by Hosokawa & Omukai (2009) for low- and high-mass central objects. We used tabulated results from their simulations that provide the evolution of the radius and luminosity of young stars as a function of their mass and accretion rate. The instantaneous accretion rate $\dot{M}$ is calculated at any time t as the variation in the accreted mass M_acc between times 0.99t and 1.01t. We assumed that the central object receives 75% of this mass to match the approximate mass ratio between the disk and the central object. We then paired the evolution tracks of Hosokawa & Omukai (2009) with Equation (27) and determine the corresponding M_* and R_* by dichotomy and interpolation of the tables. The upper panel of Fig. 5 shows the mass–radius relation for several accretion rates. The radius generally increases with mass, although central objects undergo a phase of swelling and contraction above 3 M_⊙. Those relations have been computed up to 10 M_⊙ for $\dot{M}={10}^{-6}{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$ and 90 M_⊙ for $\dot{M}={10}^{-3}{\mathrm{M}}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}$.

2.3.2 Luminosity and temperature

We interpolated the central object luminosity L_* from the tabulated results of Hosokawa & Omukai (2009), using the calculated accretion rate and mass of the central object. L_* includes both the photometric luminosity and the accretion luminosity. The bottom panel of Fig. 5 displays both the photometric and total luminosities as a function of M_* for different accretion rates. The effective temperature of central object is then calculated from L_* and R_* as $$T_{\ast }={\left(\frac{L_{\ast }}{4\pi {\sigma }_{\mathrm{B}}{R}_{\ast }^2}\right)}^{\frac{1}{4}}.$$(28)

2.4 The disk

2.4.1 The magnetic field and disk radius

Magnetic fields are important for the formation and evolution of protoplanetary disks. The field is coupled to the gas and gets twisted by its rotation. That creates a magnetic tension force that opposes the rotation, which is called magnetic braking. In an ideal magnetohydrodynamics (MHD) description, the magnetic field is perfectly coupled to the gas. The field piles-up in the central region and creates a strong braking, preventing the formation of a rotationally supported disk (RSD). The non-ideal MHD description, which accounts for dissipative and decoupling processes of the magnetic field, is needed to prevent such a strong braking and the unbounded pile-up of the magnetic field.

Lee et al. (2021) derived analytical formulae for the radius of RSDs under a full non-ideal magnetohydrodynamics (MHD) regime. They assume that at equilibrium, the accretion of angular momentum compensates exactly the loss by magnetic braking, and consider relevant non-ideal MHD effects in different regimes. In the present model, we do not explicitly calculate the magnetic field geometry, but we can estimate the strength of the magnetic field in the disk. As highlighted by Masson et al. (2016) and Vaytet et al. (2018), the magnetic field strength reaches a plateau in the disk region due to its dissipation by diffusive effects. In the envelope, the magnetic field intensity increases with the contraction of the gas. Following Nakano et al. (2002), we have $$B=\frac{1}{\lambda }4{\left(\frac{\pi k_{\mathrm{B}}T\rho }{\mu m_{\mathrm{p}}}\right)}^{\frac{1}{2}}=4.27\times {10}^{-2}\mathrm{G}{\left(\frac{\rho }{{10}^{-13}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{\frac{1}{2}}{\lambda }^{-1},$$(29)

where k_B is the Boltzmann constant, μ = 2.31, the mean atomic mass of molecules in the gas, and m_p the mass of the proton. λ is the mass-to-flux ratio, which is defined by Mouschovias & Spitzer (1976) as $$\lambda =\frac{M/{\mathrm{\Phi }}_{\mathrm{B}}}{{\left(M/{\mathrm{\Phi }}_{\mathrm{B}}\right)}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}},$$(30)

with M the cloud mass, ${\mathrm{\Phi }}_{\mathrm{B}}=\pi r_{\text{cld}}^{2}B$ the magnetic flux and $${\left(\frac{M}{{\mathrm{\Phi }}_{\mathrm{B}}}\right)}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\frac{0.53}{3\pi }\sqrt{\frac{5}{G}}.$$(31)

We therefore have λ ∝ B⁻¹. The mass-to-flux ratio is the main parameter influencing the disk size, and can take any value > 1 (Crutcher 1999). Numerical simulations generally consider λ values between 2 and 10 for magnetized cores. We calculated the magnetic field strength at ρ = 10⁻¹³ g cm⁻³, which is the formation density of the first hydrostatic core (Larson 1969), which coincides with the formation of the disk.

Lee et al. (2021) determine two possible outcomes that are the “para-disk” and the “ortho-disk” (see also Tsukamoto et al. 2015), depending on the relative directions of the disk rotation and the Hall drift of magnetic field lines. Ultimately, they show that the ortho-disk is short-lived and that a counter-rotating para-disk eventually forms, which is also observed in simulations (Tsukamoto et al. 2017; Zhao et al. 2020). Hennebelle et al. (2016) and Lee et al. (2021) derive the radius of the para-disk $$r_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}(t)=19.2\mathrm{a}\mathrm{u}{\left[\frac{{\eta }_{\mathrm{A}}}{{10}^{19}{\mathrm{c}\mathrm{m}}^2{\mathrm{s}}^{-1}}\right]}^{\frac{2}{9}}{\left[\frac{M_{\mathrm{a}\mathrm{c}\mathrm{c}}(t)}{0.1{\mathrm{M}}_{\odot }}\right]}^{\frac{1}{3}}{\left[\frac{B}{0.1\mathrm{G}}\right]}^{-\frac{4}{9}}.$$(32)

Using our magnetic field prescription, Eq. (32) becomes $$r_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}(t)=28.0\mathrm{a}\mathrm{u}{\left[\frac{{\eta }_{\mathrm{A}}}{{10}^{19}{\mathrm{c}\mathrm{m}}^2{\mathrm{s}}^{-1}}\right]}^{\frac{2}{9}}{\left[\frac{M_{\mathrm{a}\mathrm{c}\mathrm{c}}(t)}{0.1{\mathrm{M}}_{\odot }}\right]}^{\frac{1}{3}}{\lambda }^{\frac{4}{9}}.$$(33)

The ambipolar resistivity η_A depends on the physical and chemical environments. The ambipolar resistivity is of the order of 10²⁰ cm² s⁻¹ in the envelope, but decreases by a few orders of magnitude in higher density regions such as the disk (Marchand et al. 2016). However, the disk radius depends only weakly on the resistivity, with an exponent of 2/9. While it is possible to calculate it self-consistently, it would increase the number of parameters that r_disk depends on for a minor gain. Instead, we imposed a value of η_A = 10¹⁸ cm² s⁻¹ (consistent with MHD simulations, e.g., Wurster et al. 2016, 2017; Marchand et al. 2023; Lebreuilly et al. 2023), so that the user can act more directly on the disk size through the mass-to-flux ratio λ.

Fig. 5 Upper panel: mass–radius relation of central objects for different accretion rates, from the tabulated results of Hosokawa & Omukai (2009). Bottom panel: same as upper panel for the mass-luminosity relation. The solid lines represent the total central object luminosity L*, while dashed lines represent only the photometric luminosity (excluding the accretion luminosity).

2.4.2 The disk profile

We assumed an α disk framework (Shakura & Sunyaev 1973), where α is a free-parameter chosen by the user. At a given time, t, and cylindrical radius, r_c, the viscosity of the gas is given by $$\nu (r_{\mathrm{c}},t)=\alpha c_{\mathrm{s}}(r_{\mathrm{c}},t)H(r_{\mathrm{c}},t),$$(34)

where $$c_{\mathrm{s}}(r_{\mathrm{c}},t)=\sqrt{\frac{k_{\mathrm{B}}{T}_{\mathrm{m}}(r_{\mathrm{c}},t)}{\mu m_{\mathrm{p}}}},$$(35)

is the sound speed evaluated at the mid-plane, $$H(r_{\mathrm{c}},t)=\frac{c_{\mathrm{s}}}{{\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{c}},t)},$$(36)

is the scale height of the disk at a radius r_c, and $${\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{c}},t)=\sqrt{\frac{G{M}_{\ast }(t)}{r_{\mathrm{c}}^3}},$$(37)

is the Keplerian angular velocity.

We used a power-law with an exponential cut-off to describe the surface density profile (Lynden-Bell & Pringle 1974) $$\mathrm{\Sigma }(r_{\mathrm{c}},t)={\mathrm{\Sigma }}_0(t){\left(\frac{r_{\mathrm{c}}}{r_{\text{disk\hspace{0.17em}}}(t)}\right)}^{-1}{e}^{-r_{\mathrm{c}}/r_{\text{disk\hspace{0.17em}}}(t)},$$(38)

where Σ₀(t) is calculated as a function of the disk mass $$\begin{array}{rl}{M}_{\text{disk}}(t)& ={\int }_{r_{\mathrm{i}}}^{r_{\text{disk}}}{\int }_0^{2\pi }\mathrm{\Sigma }(r_{\mathrm{c}},t)r_{\mathrm{c}}\mathrm{d}{r}_{\mathrm{c}}\mathrm{d}\varphi \\ & =2\pi {\mathrm{\Sigma }}_0(t)r_{\text{disk}}(t{)}^2(e^{-r_{\mathrm{i}}/r_{\text{disk}}(t)}-e^{-1}).\end{array}$$(39)

In APE, the user can choose to either leave the exponential cutoff in Eq. (38), or to sharply cut the disk at a radius r_disk. The volume density profile at the mid-plane can then be calculated as $${\rho }_{\mathrm{m}}(r_{\mathrm{c}},t)=\frac{\mathrm{\Sigma }(r_{\mathrm{c}},t)}{H(r_{\mathrm{c}},t)\sqrt{2\pi }},$$(40)

and the vertical profile $$\rho (r_{\mathrm{c}},z,t)={\rho }_{\mathrm{m}}(r_{\mathrm{c}},t)\exp (-\frac{z^2}{2H^2(r_{\mathrm{c}},t)}).$$(41)

In addition to the Keplerian angular velocity (Eq. (37)), the gas in an α disk slowly drifts inward due to the viscosity, with a velocity $$v_r(r_{\mathrm{c}},t)=-\frac{3}{2}\alpha \frac{c_{\mathrm{s}}^2(r_{\mathrm{c}},t)}{r_{\mathrm{c}}{\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{c}},t)}.$$(42)

2.4.3 The disk temperature

Many variables depend on the temperature of the disk, through the value of the sound speed. To calculate the temperature of the mid-plane T_m(r_c, t) (dropping the (r_c,t) notation hereafter), we followed the method of Hueso & Guillot (2005). The temperature value is a balance between the stellar and cloud ambient irradiation, the viscous dissipation of kinetic energy, and the cooling $${\sigma }_{\mathrm{B}}{T}_{\mathrm{m}}^4=\frac{1}{8}(\frac{3}{8}{\tau }_{\mathrm{R}}+\frac{1}{2{\tau }_{\mathrm{P}}})\mathrm{\Sigma }(r_{\mathrm{c}},t)\nu (r_{\mathrm{c}},t){\mathrm{\Omega }}_{\mathrm{K}}^2(r_{\mathrm{c}},t)+{\sigma }_{\mathrm{B}}{T}_{\ell }^4.$$(43)

The first term on the right-hand-side is a sum of optically thin and optically thick contributions (Nakamoto & Nakagawa 1994), with τ_R and τ_P the Rosseland and Planck mean optical depths, respectively. Similarly to the envelope, T_ℓ represents the temperature equilibrium between the stellar irradiation T_ℓ,1 and the ambient cloud temperature T_mc = 10 K, so that $$T_{\ell }^4=T_{\ell ,1}^4+T_{\mathrm{m}\mathrm{c}}^4.$$(44)

Temperature T_ℓ,1 is a function of the temperature of the central object T_* and radius R_*. It is given by Adams et al. (1988); Ruden & Pollack (1991) and reads $$T_{\ell ,1}=T_{\ast }{[\frac{2}{3\pi }{\left(\frac{R_{\ast }}{r_{\mathrm{c}}}\right)}^3+\frac{1}{2}{\left(\frac{R_{\ast }}{r_{\mathrm{c}}}\right)}^2\left(\frac{H}{r_{\mathrm{c}}}\right)(\frac{d\ln H}{d\ln r_{\mathrm{c}}}-1)]}^{\frac{1}{4}},$$(45)

where d ln H/d ln r = 9/7 is assumed to avoid numerical instability (Hueso & Guillot 2005).

As Hueso & Guillot (2005), we assumed $${\tau }_{\mathrm{P}}=2.4{\tau }_{\mathrm{R}},$$(46)

and we took $${\tau }_{\mathrm{R}}=\frac{{\kappa }_{\mathrm{R}}\mathrm{\Sigma }}{2}.$$(47)

Rosseland opacities κ_R are calculated from the DSHARP opacities (Birnstiel et al. 2018) using the local dust grain size-distribution (see Sect. 2.6). The coefficients are averaged over the size-distribution $${\kappa }_{\mathrm{R}}=\frac{{\int }_{a_{min}}^{a_{max}}n(a)m(a){\kappa }_{\mathrm{R}}(a)\mathrm{d}a}{{\int }_{a_{min}}^{a_{max}}n(a)m(a)\mathrm{d}a},$$(48)

where a_min and a_max are the minimum and maximum grain sizes of the distribution, n(a) the number density of grains of size a and m(a) their mass. Here, the choice of opacities effectively influences the temperature of the midplane of the disk where the viscosity is high (typically r ≲ 30 au). In turn, this impacts the local scale-height, H, and thus the density in the same regions (with subsequent consequences on chemical abundances). For example, higher opacities would increase the temperature and the scale-height, but decrease the density of the gas. The DSHARP opacities are the default in APE, but the user can also generate their own opacity table (the process is explained in the APE user manual provided with the code).

Grains start to sublimate at temperatures ≳750 K, effectively reducing the opacity. Based on Lenzuni et al. (1995) and similarly to Marchand et al. (2016), we considered that the grain number decreases in several stages, representing the successive sublimation of carbon, silicate and aluminum material. We assumed that grains are made of a unique material, and we modeled this phenomenon by considering a linear decrease in a given material between two temperatures T₁ and T₂. Those temperatures and relative abundance of the materials are summarized in Table 2. The temperature at which all grains have sublimated, T_evap = 1700 K, is imposed as the maximum temperature in the disk. This impacts only the innermost regions of the disk (typically ≲5 au).

The temperature Equation (43) can then be rewritten in the form $${\sigma }_{\mathrm{B}}{T}_{\mathrm{m}}^4-(Y_1(r_{\mathrm{c}},t){\kappa }_{\mathrm{R}}\left(T_{\mathrm{m}}\right)+\frac{Y_2(r_{\mathrm{c}},t)}{{\kappa }_{\mathrm{R}}\left(T_{\mathrm{m}}\right)})T_{\mathrm{m}}-Y_3(r_{\mathrm{c}},t)T_{\mathrm{m}}^{\frac{1}{2}}-Y_4(r_{\mathrm{c}},t)=0,$$(49)

where we have $$Y_1(r_{\mathrm{c}},t)=\frac{3}{64}\frac{k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}\alpha {\mathrm{\Sigma }}^2(r_{\mathrm{c}},t){\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{c}},t),$$(50) $$Y_2(r_{\mathrm{c}},t)=\frac{5}{96}\frac{k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}\alpha {\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{c}},t),$$(51) $$Y_3(r_{\mathrm{c}},t)={\sigma }_{\mathrm{B}}{T}_{\ast }^4\sqrt{\frac{k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}}\frac{1}{7{\mathrm{\Omega }}_{\mathrm{K}}(r_{\mathrm{c}},t)r_{\mathrm{c}}}{\left(\frac{R_{\ast }}{r_{\mathrm{c}}}\right)}^2,$$(52) $$Y_4(r_{\mathrm{c}},t)=\frac{2}{3\pi }{\sigma }_{\mathrm{B}}{T}_{\ast }^4{\left(\frac{R_{\ast }}{r_{\mathrm{c}}}\right)}^3+{\sigma }_{\mathrm{B}}{T}_{\mathrm{m}\mathrm{c}}^4.$$(53)

Solving Eq. (49) for T_m provides the mid-plane temperature profile, from which c_s(r_c, t), H(r_c, t), and ρ(r_c, z, t) can be deduced.

2.5 Outflow

The magneto-centrifugal forces can launch a collimated jet at the scale of the central object. This creates a conic cavity with an opening angle, γ, in which low-density gas escapes at velocities of several 10 s of km s⁻¹. We assumed that the opening angle varies with height, in accordance to observations (Velusamy & Langer 1998; Cantó et al. 2008). As in Visser et al. (2011), the cavity boundary is defined by $$r_{\mathrm{c}}(\mathrm{a}\mathrm{u})<r_{\mathrm{c}\mathrm{a}\mathrm{v}}={\left(\frac{z}{0.191\mathrm{a}\mathrm{u}}\right)}^{\frac{2}{3}}{\left(\frac{t}{t_{\mathrm{a}\mathrm{c}\mathrm{c}}}\right)}^2,$$(54)

for t > t_*, where t_acc = 2t_ff. Although the gas velocity in the jet does not affect the global structure, we assumed that this velocity is the escape velocity from the central object taken at the inner radius of the disk $$v_{\mathrm{e}\mathrm{s}\mathrm{c}}=\sqrt{\frac{2G{M}_{\ast }}{r_{\mathrm{i}\mathrm{n}}}},$$(55)

where $$r_{\mathrm{i}\mathrm{n}}=\sqrt{\frac{L_{\ast }}{4\pi {\sigma }_{\mathrm{B}}{T}_{\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{p}}^4}}.$$(56)

The density of the gas inside outflow cavity walls is much lower than that of the disk or the envelope. It is generally measured to be ~10³ to ~10⁶ cm⁻³ in both models and observations (e.g., Lefloch et al. 2015; Rivera-Ortiz et al. 2023), with densities of ~10⁴ cm⁻³ generally used in analytical models (e.g., Visser et al. 2009). Here, we assumed that the outflow density is 10⁴ cm⁻³ at z = 1000 au, and decreases with height as the cross-section of the cavity increases. $$n_{\mathrm{H}}(z)=n_{\mathrm{H},0}{\left(\frac{z}{z_0}\right)}^{-2},$$(57)

with n_H,0 = 10⁴ cm⁻³ and z₀ = 1000 au. The outflow is a low-opacity region through which the radiative flux of the central object can heat the envelope outside the cavity walls. The opening angle therefore impacts the location of the heated region, while the density affects the efficiency of this heating.

2.6 Dust grains

We included a treatment of the dust size-distribution to estimate the opacities for radiative transfer calculations. While the characteristic size of grains is sub-micron in the interstellar medium (Mathis et al. 1977), they grow by coagulation during star formation and can quickly reach radii larger than 100 μm in protoplanetary disks (Marchand et al. 2023). In APE, we have implemented the simulation results of Lebreuilly et al. (2023). They performed 1D non-ideal MHD simulations of protostellar collapses using the SHARK code (Lebreuilly et al. 2023) to track the evolution of the grain size-distribution by coagulation and fragmentation. This evolution is controlled by the collision rates between grains. Lebreuilly et al. (2023) take into account differential grain-grain velocities originating from the gas turbulence (Ormel & Cuzzi 2007), ambipolar diffusion (due to the electrical charges of grains, see Guillet et al. 2020) and Brownian motion (Bate 2022). Turbulence is efficient at growing large grains in high-density environments (as the disk), while ambipolar diffusion and Brownian motion efficiently remove small grains in the envelope. In the disk, the fragmentation limits the growth and repopulate the small-grains population.

We used the last output of the simulation PC1-FRAG of Lebreuilly et al. (2023) to associate gas densities with size-distributions. The simulation was originally performed using 100 bins in size between 7 nm and 1 cm, with a dust-to-gas mass ratio of 1%. This large number of bins may however result in expensive radiative transfer calculations, and we resampled the distribution into 20 bins between 16 nm and 3.9 mm. In APE, we provide the resampled distribution and the associated opacities for radiative transfer post-processing. Size-distributions at ρ = 10⁻¹⁷, 10⁻¹⁴ and 10⁻¹¹ g cm⁻³ are displayed in Fig. 6. In the envelope (represented by the solid line), the distribution is close to the initial Mathis, Rumpl, Nordsieck (MRN) distribution (Mathis et al. 1977), with the peak of the distribution shifting toward larger sizes as the density increases.

The results of Lebreuilly et al. (2023) are however only valid for a dust-to-gas mass ratio of 0.01. As users may desire different dust-to-gas mass ratios, we provide an alternative grain coagulation model. Marchand et al. (2021) show that, when considering only collisions from turbulence, the dust coagulation can be treated as a 1-dimension problem. Size-distributions depend then only on the initial distribution, and one variable that they name χ, which is a function of the density and temperature history of the grains. This variable also scales linearly with the dust-to-gas mass ratio. 3D simulations with non-ideal MHD by Marchand et al. (2023) present the evolution of χ and the dust size throughout the collapse. From their simulation, we derived a prescription for χ at a dust-to-gas mass ratio of 0.01 $$\chi (0.01)=\{\begin{array}{ll}{10}^{16}{\left(\frac{\rho }{3\times {10}^{-20}}\right)}^{\frac{1}{4}}& (\rho \le {10}^{-15}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3})\\ 1.34\times {10}^{18}{\left(\frac{\rho }{3.83\times {10}^{-12}}\right)}^{\frac{1}{15}}{\left(\frac{T}{10}\right)}^{\frac{1}{15}}& (\rho >{10}^{-15}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3})\end{array}.$$(58)

χ is then linearly scaled by the desired dust-to-gas mass ratio, dtg. $$\chi (dtg)=\chi (0.01)\frac{dtg}{0.01}.$$(59)

The prescription is discontinuous as the growth of grains is fast at densities ρ > 10⁻¹⁵ g cm⁻³ (Marchand et al. 2023). The outflow is composed of material lifted from regions in the vicinity of the disk, typically at densities close to ρ = 10⁻¹⁵ g cm⁻³ (Marchand et al. 2023); we therefore considered that the value of χ in the outflow matches the one at density ρ = 10⁻¹⁵ g cm⁻³. We have generated a tabulated size-distribution evolution path as a function of χ using the Ishinisan code (Marchand et al. 2021). We provide this evolution path with APE. The dust abundance is then read from the table using the local values of χ.

In APE, choosing a dust-to-gas mass ratio exactly equal to 0.01 is equivalent to choosing the coagulation model of Lebreuilly et al. (2023). Any other value of dtg switches the coagulation model to Marchand et al. (2023). The user can also choose to turn off the coagulation to only have a MRN distribution.

As grains grow, their total surface area decreases, which can impact the chemical evolution. However, to our knowledge, no chemical code currently supports an evolving grain size-distribution. Nautilus, in particular, supports only a fixed grain size with a fixed dust-to-gas mass ratio. A few attempts at emulating the evolution of grains in chemical networks have been published. For example, Iqbal & Wakelam (2018) implemented a fixed, full size-distribution in Nautilus. They found that differences in the chemical evolution between two size-distributions could be small as long as the total surface areas were close, but could become significant otherwise or in the presence of cosmicrays. More recently, Navarro-Almaida et al. (2024) modified Nautilus by computing an effective grain size at each step. They found that the main impact of the decrease in grain surface area is seen in dense regions, below 100 K, as many species are still in solid phase. APE has been designed to interface with any chemical code, preferentially with the public version of Nautilus. Tools and methods to account for an evolving grain size-distribution in chemical simulations are still being developed and, to this day, none are publicly available. Nonetheless, this is a crucial aspect of the chemistry of star formation, and APE is able to provide evolving size-distributions to future codes that include this physics. Until such implementations become available, we do not account for grain growth in our chemical runs.

Fig. 6 Grain size-distributions at various densities. The y-axis represents the mass density of the grains normalized by the mass density of the gas.

2.7 Extinction

APE has been designed to work with the three-phase kinetic chemistry code Nautilus (Ruaud et al. 2016) by generating necessary input files. Nautilus accounts for the impact of external UV emission on the chemical species, and requires an extinction value. We estimated the extinction, A_v, through the column density N_H2 of the gas between the point of interest and the outside of the cloud (Semenov et al. 2010) $$A_{\mathrm{v}}=\frac{N_{{\mathrm{H}}_2}}{1.59\times {10}^{21}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-2}}.$$(60)

We assumed a power-law density profile in the envelope, ρ(r) ~ r^s (with s converging toward −1.5, see Sect. 2.2.2). After estimating s, an integration is performed toward the outer-edge of the envelope to obtain the column density $$N_{{\mathrm{H}}_2,\mathrm{e}\mathrm{n}\mathrm{v}}(r)={\int }_r^{r_{max}}\rho \left(r^{\mathrm{\prime }}\right)\mathrm{d}{r}^{\mathrm{\prime }},$$(61) $$=\frac{\rho (r)}{s+1}\frac{r_{max}^{s+1}-r^{s+1}}{r^s},$$(62) $$=\exp [\ln (\rho (r))+\ln (r)-\ln (-(s+1))+\ln (1-{\left(\frac{r_{max}}{r}\right)}^{s+1})].$$(63)

Equations (62) and (63) are equivalent, but we employed the latter in APE to avoid numerical precision issues. We have added a negative sign to the last two terms of Eq. (63) because s + 1 < 0.

In the disk, for the gas at cylindrical radius r_c and height z, the density profile is integrated vertically until one scale-height (dropping the dependency in t): $$N_{{\mathrm{H}}_2,\text{disk}}(r_{\mathrm{c}},z)={\int }_z^{H\left(r_{\mathrm{c}}\right)}{\rho }_{\mathrm{m}}\left(r_{\mathrm{c}}\right)\exp (-\frac{z^{\mathrm{\prime }2}}{2H^2\left(r_{\mathrm{c}}\right)})\mathrm{d}{z}^{\mathrm{\prime }},$$(64) $$=\sqrt{\frac{\pi }{2}}{\rho }_{\mathrm{m}}\left(r_{\mathrm{c}}\right)H\left(r_{\mathrm{c}}\right)[1-\mathrm{erf}\left(\frac{z}{\sqrt{2}H\left(r_{\mathrm{c}}\right)}\right)],$$(65)

with erf the error function. We added a flat value to account for the external cloud extinction. This value is a parameter, which can be chosen by the user. Finally, the total extinction is the sum of the three contributions $$A_{\mathrm{v}}(r,r_{\mathrm{c}},z)=A_{\mathrm{v},\text{env}}(r)+A_{\mathrm{v},\text{disk}}(r_{\mathrm{c}},z)+A_{\mathrm{v},\text{cloud}}.$$(66)

2.8 APE workflow

Unlike hydrodynamics simulations, determining the physical properties at a specific location in APE does not require the description of the whole history of the gas. For instance, the density at location (r,θ) at time t is calculated with no knowledge of the gas surrounding this location or the system at earlier time points. The calculation solely relies on key characteristics of the central object (mass, radius, luminosity, temperature) and the disk (mass, size, mid-plane profile). This feature significantly accelerates the computation: a snapshot at time t can be generated without computing the whole gas history since t = 0, and a particle trajectory can be traced without computing the whole map at each time-step. Hence, the grid resolution does not influence the physical result (although it impacts the post-processing calculations like radiative transfer and synthetic observations). The APE workflow is summarized in Fig. 7.

2.9 The trajectory of particles

In particle mode, particles start with an initial user-defined position and move according to the local velocity field. The default time-step is dt₀ = 200 years (this value can be changed by the user). However, the time-step is limited when the velocity of the particle is high or if the particle is close to the edges of the simulation box. The time-step is calculated as $$dt=min(\frac{1}{f}\frac{r}{v_r},d{t}_0),$$(67)

with r the spherical radial coordinate, v_r the radial velocity, and f an arbitrary safety factor. Dynamically, this factor is especially important in the inner envelope, where radial velocities can reach several km s⁻¹. Conversely, the radial velocity in the disk is low (~10⁻³ km s⁻¹), and the particle trajectory is little sensitive on the factor f. However, the value of the time-step may influence the chemical evolution of the particle. We performed a convergence study on the influence of f upon chemical abundances. We ran APE with a particle starting at position (x = r_c = 5 au, z = 1 au) (in the disk) at t = t_ff + 40 kyr, in reverse mode until t = 0, with dt₀ = 200 years. The system parameters are the same as the fiducial case presented in Sect. 3.1.2, with a cloud mass of 2 M_⊙ and a mass-to-flux ratio λ = 5. The particle spends ~150 kyr in the envelope and ~40 kyr in the disk. This final position has be chosen so that the particle travels across a wide range of physical conditions, from the envelope to the inner disk, and is thus more impacted by changes in the time-step. We investigated values of f of 10, 100, 1000 and 10 000. Dynamically, the trajectories are very close, with less than a 0.4% difference in their final positions, as listed in Table 3. We ran the Nautilus code on those trajectories using the publicly provided chemical network. In Fig. 8, we display box plots of the relative differences of abundance for all chemical species in the network, using f = 10 000 as the reference case. For f=10 and 100, half of the species have less than a 4% difference with respect to the reference case. Over 95% of species have less than a 15% difference. Nearly all species have less than a 50% difference, except 4 species in gas and solid phase for f=10 (FeOH, HFeOH, HOFeOH, and FeHOHOH), and 2 species for f=100 (NH₂PO and HOFeOH), while all within a factor 2. The case f=1000 is well converged with all species displaying less than 5% difference with f=10 000, except HCONO with a 8% difference and NH₂PO with a 90% difference. For chemical convergence, we therefore recommend to use 1000 ≥ f > 100, but choosing f=10 results in a ~50% maximum error on abundances while requiring much less computational resource. We reach the same conclusions for a particle starting at (x=1 au, z=5 au) and staying exclusively in the envelope.

Fig. 7 Workflow of the APE code.

Fig. 8 Distributions of the relative differences in abundances for several values of f, using f=10 000 as reference. The horizontal black lines represent the medians of the distributions, the colored boxes range from the first to the third quartiles, and the vertical black lines span from the minimum to the maximum values. For the f = 1000 case, the colored box is not visible as its upper limit stands at a value of 0.001.

2.10 Producing synthetic observations

The purpose of APE is to generate physical conditions to run chemical models and synthetic observations. Producing synthetic observations requires simulating the telescope view of an emission map, which is itself computed by radiative transfer calculations from abundance maps. While the users of APE are free to use any software to perform those steps, we only provide scripts to the publicly available codes Nautilus (Ruaud et al. 2016), RADMC-3D (Dullemond et al. 2012), and Imager³ to facilitate the process.

For molecular line observations, the first step is the production of chemical abundance maps. The “grid of particles” mode of APE places virtual particles in each cell of the map at a given time t_final (generated in the same way as in a snapshot), and computes their trajectory backward in time until the initial condition t = 0 is reached. The history of each particle can then be used as input for Nautilus, alongside initial abundances and a chemical reaction network, to compute its chemical evolution. The chemical abundance maps can then be reconstructed into a grid, similarly to a density snapshot, from the final abundances of all particles at time t_final.

The next step is the computation of emission maps. To that end, RADMC-3D needs information about the grid data (automatically generated by APE), spectroscopic data (user-provided) and observational parameters. Those parameters are user-defined and include the observation wavelength (alternatively which molecular transitions to observe), and the source inclination. RADMC-3D then reconstructs the 3D spatial distribution of the molecule from the axisymmetric APE data, and produces an emission cube. Finally, Imager is run on the emission map with the “simulate” function to produce cubes of observations.

With APE, we provide a script, available both in Python and Julia, in which the user can choose all observation parameters: the species, the frequency bandwidth, the spectral resolution, the inclination, the source distance and declination, the ALMA configuration, the observing time, and the noise. The script is an interface between the outputs of APE (once the chemical abundance maps have been generated), RADMC-3D, and Imager. The final output is a FITS file with the clean cube, optionally in which the continuum has been subtracted.

3 Application

3.1 Computing abundance maps

3.1.1 Method

In this section, we illustrate an application of APE, which is to compute the abundance map of several molecules, using the “Grid of particles” mode described in Sect. 2.10. The temperature is interpolated in time and space along the trajectory of particles using RADMC-3D as explained in Sect. 2.1.2.

3.1.2 Simulation setup

We considered a cloud with an initial mass of 2 M_⊙, and produce a snapshot 150 kyr after the formation of the central object, which occurs at t = t_ff ≈ 152 kyr. The disk size is determined by the mass-to-flux ratio, which we chose to be λ = 5 as it is a standard value in protostellar collapse simulations (e.g., Masson et al. 2016; Wurster et al. 2016; Marchand et al. 2023). The dust-to-gas mass ratio is 1% and we enabled dust coagulation. The initial angular velocity of the cloud is set at Ω₀ = 2 × 10⁻¹⁵ s⁻¹. The grid is composed of 5625 cells, with a logarithmic radial sampling of 75 points between r = 1 au and r = 1000 au, and a uniform poloidal sampling of 75 points between θ = 0 and θ = π/2. We used a limiting factor f = 300 on the time-step (see Sect. 2.9).

3.1.3 Chemical model

Each particle history is used as input for the Nautilus code (Ruaud et al. 2016), which performs three-phase gas-grain chemical calculations (gas phase – grain surface – grain mantle). The chemical network used in the Nautilus code contains about 800 species and ~9000 individual reactions (Manigand et al. 2021; Coutens et al. 2022; Agúndez et al. 2023). The grain surface and the mantle are both chemically active for these simulations. A sticking probability of 1 is assumed for all neutral species while desorption can occur by thermal and non-thermal (cosmic-ray, chemical desorption) processes including sputtering of ices by cosmic-ray collisions (Wakelam et al. 2021). Surface reactions formalism and more detailed description of the simulations can be found in Ruaud et al. (2015). The diffusion energies are set to be a fraction of the binding energy for the surface (0.4 times) and the mantle (0.8 times), with a diffusion barrier thickness of 0.25 nm.

Initially, we ran Nautilus for t_pre = 10⁶ years under constant physical parameters characteristic of prestellar core conditions: a density of n₀ = 10⁴ cm⁻³, a temperature of T₀ = 10 K, and an extinction of A_v = 5. We assumed a cosmic-ray ionization rate of ζ_CR = 1.3 × 10⁻¹⁷ s⁻¹. We started this run using the same abundances as Ruaud et al. (2018), which we detail in Table 4 (with the exception of Mg that is not included in the chemical network). In the rest of the paper, we call “initial abundances” the final abundances from this run, which serve as the initial conditions for the chemical simulations along each particle trajectory. Since Nautilus does not support yet an evolving grain size-distribution (although a first attempt has been recently published by Navarro-Almaida et al. 2024), we used a unique, constant grain size that we assumed to be 0.1 μm, with a dust-to-gas mass ratio of 0.01. The bulk density of the grains is 3 g cm⁻³, and their surface site density of 1.5 × 10¹⁵ cm⁻². In this work, we did not compute the chemical abundances in the outflow cavity, where they are set equal to the initial abundances. In this study, we considerd that the sole impact of the outflow is on temperature maps, which can in turn affect chemical abundances.

3.1.4 System evolution

At t = t_ff + 150 kyr, the envelope has nearly vanished with a remaining mass of 0.05 M_⊙. The central object has accreted 1.5 M_⊙, while the disk has the remaining mass of 0.45 M_⊙, and a radius of 92.3 au. Figure 9 displays the density and temperature map of the system zoomed on the inner 200 au. The protostellar radiation hardly penetrates into the outer regions of the disk, resulting in a local minimum of temperature.

The upper panel of Fig. 10 displays the density and temperature history of a particle located in the envelope at a distance of ~20 au from the central object at the final time of t = t_ff + 150 kyr. Figure 11 is the same for a particle ending in the disk. The envelope particle starts at a density of ~10⁴ cm⁻³, which increases only by one order of magnitude over 250 kyr. The density increase then proceeds faster as the particle draws close to the central object, to reach ~10⁸ cm⁻³ at the end of the simulation. The density of the disk particle starts at ~10⁵ cm⁻³ and slowly increases to ~10⁷ cm⁻³, until it reaches the inner envelope (≲100 au) at the time of formation of the central object (t − t_ff = 0). This particle then enters the disk and experiences a density jump, before slowly drifting inward until the final time. The temperature of both particles is initially constant at T = T_mc = 10 K. The particles are then close enough to the central object at its formation to experience a jump of temperature, at t − t_ff = 0 kyr. The disk particle immediately reaches T ≈ 20 K, then heats up to T ≈ 80 K at its entry in the disk. Its temperature then increases slowly while drifting inward, reaching T ≈ 120 K at the end of the simulation. The envelope particle is further away at the formation of the central object, and its temperature only jumps to T ≈ 12 K, and then heats up increasingly faster the following 150 kyr until it also reaches T ≈ 120 K.

Fig. 9 Zoomed map of the system at t = tff + 150 kyr. Colors represent the number density of the gas, while white lines indicate iso-contours of the temperature at T=50, 75, 100 and 300 K. The red contours delimitate the envelope, the disk and the outflow.

Fig. 10 Upper panel: density (blue) and temperature (red) history as a function of time for a particle ending in the inner envelope (x, y = 14.0, 14.0 au) at its last time-step at t = tff + 150 kyr. Bottom panel: time evolution of the gas phase abundance of several species for this particle.

Fig. 11 Same as Fig. 10 for a particle ending in the disk (x, y = 19.1, 2.0 au).

3.1.5 Chemical abundances

The density and temperature evolution are then used as input for the chemical calculations. As illustration, the bottom panels of Figs. 10 and 11 show the time evolution of the abundance (relative to H) of several species: carbon monoxyde (CO), carbon monosulfide (CS), cyanide (CN), formaldehyde (H₂CO), and methanol (CH₃OH) in the gas phase for the same particles as in the upper panels. The chemical evolution of the envelope particle (Fig. 10) is slow, with variations smaller than one order of magnitude, until a few 10 kyr before the end of the simulation. Close to the protostar, the abundances then vary rapidly: the abundance of CN decreases by 5 orders of magnitude, while the abundance of other molecules increase by 1 to 4 orders of magnitude. The disk particle (Fig. 11) starts its journey closer to the center of the cloud. At the formation of the central object, the jump in temperature causes a ~1 order of magnitude increase in the gas phase abundance of H₂CO and CH₃OH, a decrease by the same factor for CN, while CO and CS are little affected. CO sees a decrease in abundance until the disk entry, where it suddenly increases by 2 orders of magnitude and stays constant until the end. H₂CO and CS exhibit a similar behavior with jumps of 4 orders of magnitude. In the last kyrs of the simulation, the abundance of CS sharply decreases by 5 orders of magnitude. Conversely, the abundance of CN drops by 8 orders of magnitudes at the disk entry, and continues decreasing by 2 orders of magnitude over the following ~100 kyr. The density jump at the disk entry also causes a sudden decrease in CH₃OH gas phase abundance by 2 orders of magnitude. As the particle drifts inward in the disk, the increasing temperature provokes a slow desorption of CH₃OH in the gas phase over time. At the final time, all CH₃OH molecules present on grains have desorbed and their abundance is 5 orders of magnitude higher than before the disk entry.

We performed similar calculations for all particles in the simulation box, and use the final abundance of each particle to reconstruct abundance maps. We display resulting abundance maps at t = t_ff + 150 kyr in Fig. 12, for those five species in the gas phase. The abundance of CO varies negligibly across the envelope. The whole computed map has a temperature larger than 30 K, which is above its desorption temperature. Conversely, CH₃OH is mostly frozen on grains at temperatures below 100 K. There is a clear correlation between the regions of high methanol abundance and the temperature isocontours. This map only display the gas phase abundance of methanol; the total gas+ice abundance is roughly constant throughout the collapse. This map therefore effectively illustrates the thermal desorption of the species. H₂CO and CS exhibit a similar behavior, with an increased abundance along the outflow cavity. On the contrary, the abundance of CN decreases in the high-temperature regions, especially above 50 K.

3.2 Synthetic observations

3.2.1 Observation parameters

In this section, we present synthetic observations of the molecular emission of the species presented in Sect. 3.1: CO, CN, CS, H₂CO, and CH₃OH, using our simulation at t = t_ff + 150 kyr. For qualitative comparison, we mimicked the observation parameters used by Podio et al. (2020) to observe those five molecules toward the Class I protostar IRAS 04302+2247 with ALMA. We selected spectral windows in the band 6 of ALMA that were centered on the transition frequencies listed in Table 5. The first step is generating the emission maps using RADMC-3D. The main input parameters are the abundance maps, the inclination angle, and the spectral resolution, which we set at 141 kHz (~0.18 km s⁻¹). The spectral window spans from −10 to 10 km s⁻¹, which results in a sampling of ~110 channels. The emission map is then used as input for Imager to simulate ALMA observations. To match the source parameters, we assumed a 90 degrees inclination of the source (edge-on disk), a source distance of 160 pc, and a declination of +22 degrees. We simulated an observation with the C-5 configuration of ALMA with a 13 minutes integration time, resulting in a beam size of 0.32 × 0.25″. The noise in synthetic observations is purely numerical. The observing time therefore only affects the uv-coverage.

3.2.2 Results

We display the moment 0 map of the molecular emission in Fig. 13 for the inner 6″ × 6″ region. We see an X-shaped emission for CO, CS, and H₂CO. This structure corresponds to the hot regions along the cavity of the outflow, as displayed for those species in Fig 12. Those X structures are extremely similar to what is observed in IRAS 04302+2247 (see Fig. 1 of Podio et al. 2020). As in Podio et al. (2020), not much can be seen in the CN emission. However, the CH₃OH emission is faint in IRAS 04302+2247, while we clearly see a bright X-shape at the ~1″ scale in our simulation. The low bolometric luminosity of 0.43 L_⊙ measured toward IRAS 04302+2247 (Ohashi et al. 2023) may not be enough to allow a significant desorption of methanol in the envelope.

Fig. 12 Reconstructed maps of relative gas phase abundances of CO (top left), CN (top middle), CS (top right), H2CO (bottom left), and CH3OH (bottom right), 150 kyr after the formation of the central object. The black contours represent temperatures of 40, 50, 75,100 and 300 K. The red contours delimitate the envelope, the disk and the outflow.

Fig. 13 Moment 0 map of the molecular emission for an edge-on configuration, zoomed on the inner 6″. The white ellipse on the bottom left represents the beam size. From left to right, top to bottom: CO 2–1, CN 2–1, CS 5–4, H2CO 31,2−21,1, CH3OH 50,5−40,4.

4 Conclusion

In this work, we presented the Analytical Protostellar Environment (APE) code, which we made publicly available. APE is a tool that facilitates chemical simulations and synthetic observations of protostellar systems. The code contains a physical model able to simulate the gravitational collapse of a dense core, and the formation and evolution of a central object and a protoplanetary disk throughout the Class 0 and I phases. We also provide scripts and interfaces to the publicly available codes Nautilus (for chemical simulations), RADMC-3D (for radiative transfer calculations), and Imager (for synthetic imaging).

We illustrated the capabilities of APE by computing the abundance maps of CO, CN, CS, H₂CO, and CH₃OH in a Class I protostellar system. We also performed synthetic ALMA observations of those species assuming an edge-on source. The moment 0 maps of CO, CS and H₂CO show an X-shaped structure extremely similar to the emission seen toward the Class I protostar IRAS 04302+2247 (Podio et al. 2020).

Data availability

APE is publicly available on Bitbucket⁴. A user manual is provided with the code to detail its usage.

Acknowledgements

We thank the anonymous referee for their comments that helped us improve the clarity of the manuscript. We thank Takashi Hosokawa for kindly allowing us to publicly use and provide his tabulated simulation results. We also thank Ugo Lebreuilly for providing his simulation results for the dust evolution. This study is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 949278, Chemtrip).

Appendix A Benchmark of the envelope model

Since APE serves as an alternative to MHD simulations, we compare the results obtained by both methods in this section. In classical MHD simulations, the envelope collapses, then forms an opaque first hydrostatic core with a ≲20 au radius (Larson 1969). This core contracts adiabatically for at most a few thousand years, before collapsing into a central object. A circumstellar disk forms in and around the first core and grows in size. APE does not model the first core, nor the initial growth of the disk, and the computational cost of MHD simulations increases significantly at the formation of the central object (which is the main motivation for the creation of APE). For those reasons, we restricted our comparison to the envelope before the ignition of the central object.

We used the RAMSES code (Teyssier 2002) with its ideal MHD solver (Fromang et al. 2006). The initial condition is a 2 M_⊙ critical Bonnor-Ebert sphere, with an initial radial velocity of the gas matching that of APE $$v_r(r)=-\sqrt{\frac{G{M}_{\mathrm{i}\mathrm{n}}(r)}{2r}}.$$(A.1)

We performed a run with a moderate initial magnetic field (mass-to-flux ratio λ = 5). Once the initial condition is set up, the gas is subjected to the effects of gravity, thermal pressure and the Lorentz force. The temperature is imposed at T = 10 K and the extinction at A_v = 3.5.

We compared the evolution of particles ending at 50 au from the center in the midplane right before the ignition of the central object. The density history of the two particles is displayed in Fig. A.1. The initial ~100 kyr are nearly identical. The formation of a central object is however ~15% longer in the MHD simulation, which results in a slightly lower density after 100 kyr. Both particles then reach a density of ~4 × 10⁹ cm⁻³, at a time of ~152 kyr for the APE simulation, and ~175 kyr for the RAMSES simulation.

Fig. A.1 Time evolution of the density for a particle ending at 50 au just before the ignition of the central object, for the APE (solid line) and the RAMSES (dashed line) simulations.

We then ran Nautilus on the two particles, using the same chemical network and initial abundances as described in Sect. 3.1.3. As example, Fig. A.2 displays the abundance evolution of the same five species considered in Sect. 3.2: CO, CN, CS, H₂CO, and CH₃OH. The time axis is normalized by the formation time of the central object in each simulation. Similarly to the density, the abundances are nearly identical in both simulations for the first ~100 kyr (t/t_ff ≈ 0.6). As the abundances decrease, differences between the two cases increase up to an order of magnitude for the simplest species (CO, CN, and CS), but are less pronounced for the two more complex species (H₂CO and CH₃OH) with differences less than a factor of 2.

Fig. A.2 Time evolution of the gas-phase abundance of CO, CN, CS, H2CO, and CH3OH for APE (solid lines) and RAMSES (dashed lines). The time has been normalized with the formation time of the central object in each simulation.

In the following, we considered the full chemical network. At the final time-step, all species with non-negligible abundances (> 10⁻¹⁰) show differences lower than 50% between the two simulations (with 230 out of 240 being below a 10% difference). The only exception is gas-phase H, which is 10 times more abundant with APE (~10⁻⁹) than with RAMSES (~10⁻¹⁰). Including species with abundances above 10⁻¹⁵, 90% display differences smaller than 25%, and all remain within an order of magnitude.

Overall, APE successfully reproduces the chemical trends obtained from a MHD simulation, particularly for the most abundant species, whose abundances show strong convergence. This demonstrates that APE provides a reliable alternative for studying the chemical evolution of star-forming systems, at a fraction of the cost of full MHD simulations.

References

All Tables

All Figures

	Fig. 1 Schematic of an APE snapshot simulation.
In the text

	Fig. 2 Top panel: density map generated by APE. Bottom panel: zoom of the inner 200 au, with a reduced color sampling.
In the text

	Fig. 3 Normalized density profile of a Bonnor-Ebert sphere.
In the text

	Fig. 4 Density profile of a collapsing Bonnor–Ebert sphere following our model, for a 2 M⊙ sphere at time t ≈ 192 kyr (in blue). The dashed black line represents a $r^{-\frac{3}{2}}$ profile.
In the text

	Fig. 5 Upper panel: mass–radius relation of central objects for different accretion rates, from the tabulated results of Hosokawa & Omukai (2009). Bottom panel: same as upper panel for the mass-luminosity relation. The solid lines represent the total central object luminosity L*, while dashed lines represent only the photometric luminosity (excluding the accretion luminosity).
In the text

	Fig. 6 Grain size-distributions at various densities. The y-axis represents the mass density of the grains normalized by the mass density of the gas.
In the text

	Fig. 7 Workflow of the APE code.
In the text

Fig. 8 Distributions of the relative differences in abundances for several values of f, using f=10 000 as reference. The horizontal black lines represent the medians of the distributions, the colored boxes range from the first to the third quartiles, and the vertical black lines span from the minimum to the maximum values. For the f = 1000 case, the colored box is not visible as its upper limit stands at a value of 0.001.

In the text

	Fig. 9 Zoomed map of the system at t = tff + 150 kyr. Colors represent the number density of the gas, while white lines indicate iso-contours of the temperature at T=50, 75, 100 and 300 K. The red contours delimitate the envelope, the disk and the outflow.
In the text

	Fig. 10 Upper panel: density (blue) and temperature (red) history as a function of time for a particle ending in the inner envelope (x, y = 14.0, 14.0 au) at its last time-step at t = tff + 150 kyr. Bottom panel: time evolution of the gas phase abundance of several species for this particle.
In the text

	Fig. 11 Same as Fig. 10 for a particle ending in the disk (x, y = 19.1, 2.0 au).
In the text

	Fig. 12 Reconstructed maps of relative gas phase abundances of CO (top left), CN (top middle), CS (top right), H2CO (bottom left), and CH3OH (bottom right), 150 kyr after the formation of the central object. The black contours represent temperatures of 40, 50, 75,100 and 300 K. The red contours delimitate the envelope, the disk and the outflow.
In the text

	Fig. 13 Moment 0 map of the molecular emission for an edge-on configuration, zoomed on the inner 6″. The white ellipse on the bottom left represents the beam size. From left to right, top to bottom: CO 2–1, CN 2–1, CS 5–4, H2CO 31,2−21,1, CH3OH 50,5−40,4.
In the text

	Fig. A.1 Time evolution of the density for a particle ending at 50 au just before the ignition of the central object, for the APE (solid line) and the RAMSES (dashed line) simulations.
In the text

	Fig. A.2 Time evolution of the gas-phase abundance of CO, CN, CS, H2CO, and CH3OH for APE (solid lines) and RAMSES (dashed lines). The time has been normalized with the formation time of the central object in each simulation.
In the text