© The Authors 2023

1 Introduction

Observations and numerical magnetohydrodynamical (MHD) models of the formation of massive protostars along with their accretion disks and collimated jets have made undeniable progress over the last decade. Even though there are several observational examples of jets from massive protostars (see e.g., Purser et al. 2016), direct observations and estimations of the radii of their associated accretion disks are rarer. The source for the HH 80–81 radio jet, for example, has an accretion disk with an estimated radius of 950–1300 au, and an inner radius of ≲170 au, with the mass of the protostar in the range of 4–18 M_⊙ (Girart et al. 2017; Carrasco-González et al. 2010).

Recently, a new generation of multiscale observations of the star-forming region IRAS 21078+5211 (Moscadelli et al. 2021, 2022) have revealed the existence of a disk-jet system around a protostar of 5.6 ± 2 M_⊙. Emission peaks from molecular rotational transitions of CH₃CN and HC₃N trace the kinematical footprint of a Keplerian-like accretion disk of around 400 au in diameter surrounding the protostar. For the first time, the direct study of the launching mechanism of the jet has been possible, thanks to water maser observations performed with very long base interferometry (VLBI; Moscadelli et al. 2022), which traced individual streamlines in the jet and reached a resolution of 0.05 au.

In Moscadelli et al. (2022), we modeled the disk-jet system in IRAS 21078+5211 with a resistive-radiation-gravito-magnetohydrodynamical simulation of a forming massive star starting from the collapse of a cloud core. We followed the evolution of the protostar, the formation of the accretion disk, and the launching, acceleration, collimation, propagation, and termination of the magnetically driven outflows, which self-consistently form within the computational domain. Thanks to the use of an axisymmetrical grid, we were able to reach unprecedented resolutions in the regions surrounding the forming massive star (up to 0.03 au), allowing for the comparison with the VLBI data (which has a resolution of 0.05 au). We found good agreement between the simulated and observed disk size (given the protostellar mass) as well as the main features of the jet. The simulation was part of a larger study that examines the effects of the natal environment on the formation of a massive star. We present the full series of simulations in this article, with a special focus on the dynamics of the accretion disk, while the dynamics of the magnetically driven outflows are the focus of an upcoming article (Oliva & Kuiper 2023, hereafter Paper II).

Previously, the numerical simulations performed by Kölligan & Kuiper (2018) showed the self-consistent magneto-centrifugal launching of the jet and a tower flow driven by magnetic pressure in the context of massive star formation. The present study builds on their work by improving the treatment of the thermodynamics of the gas and dust with the inclusion of radiation transport with the gray flux-limited diffusion approximation (Kuiper et al. 2020), and increasing the resolution of the grid. Those changes enable us to examine the vertical structure of the accretion disk. Previous numerical studies of the formation of a massive star in an environment with magnetic fields did not consider magnetic diffusion (e.g., Myers et al. 2013; Seifried et al. 2011) or approximate the thermodynamics of the system with the use of a barotropic (e.g., Machida & Hosokawa 2020) or isothermal (Kölligan & Kuiper 2018) equation of state, instead of solving the radiation transport equations. Recent studies by Mignon-Risse et al. (2021) and Commerçon et al. (2022), which include ambipolar diffusion as the nonideal MHD effect and as well as radiation transport, find indications of a magneto-centrifugally launched jet (although the grid used in that study does not properly resolve the launching region), and smaller disks with stronger magnetic fields. The question of the size of the disk is explored in more detail in the present article, as we find an opposite trend in our results. For a further literature review, we refer the reader to Sect. 7.

Because of computational costs, previous efforts have focused on only a select number of natal environmental conditions for the massive (proto)star. Our setup, in turn, allows us to explore a wide range of natal environmental conditions and their impact on the outcome of the formed massive star and its disk-jet system. In Sect. 2 we describe the model, physical effects considered, and the parameter space of the simulations. We present an overview of the evolution of the system based on the fiducial case of the parameter space in Sect. 3. In Sect. 4, we focus on the dynamical processes observed in the accretion disk, that is to say the effect of the magnetic field, Ohmic dissipation, and radiation transport. In Sect. 5 we explore the dependence of the disk size and evolution with the initial conditions for the gravitational collapse, and finally in Sect. 7 we offer a comparison of our results with previous numerical studies.

2 Model and method

2.1 Description of the model

We model axisymmetrically the gravitational collapse of a cloud core of mass M_C and a radius of R_C = 0.1 pc (cf. Fig. 1), which has an initial density profile of the form $$\rho \left(r,t=0\right)={\rho }_0{\left(\frac{r}{r_0}\right)}^{{\beta }_{\rho }}.$$(1)

The constant r₀ is chosen to be 1 au and ρ₀ is determined by computing the mass of the cloud with the integral of the density profile over 0 < r < R_C. The cloud is in slow initial rotation, given by an initial angular velocity profile of the form $$\text{Ω}\left(R,t=0\right)={\text{Ω}}_0{\left(\frac{R}{R_0}\right)}^{{\beta }_{\text{Ω}}},$$(2)

where R is the cylindrical radius and R₀ is chosen to be 10 au. The constant Ω₀ is fixed by computing the ratio of rotational to gravitational energy ζ ≡ E_r/E_g from the initial density and rotation profiles; ζ is then a parameter of the simulation.

The cloud core is threaded by an initially uniform magnetic field directed parallel to the rotation axis. Its magnitude B₀ is determined by the normalized mass-to-flux ratio $\overline{\mu }$, $$\overline{\mu }\equiv \frac{M_C/{\Phi }_B}{{\left(M_C/{\Phi }_B\right)}_{\text{crit}}}.$$(3)

The denominator of the previous expression represents the critical value of the mass-to-flux ratio, that is the value for which the gravitational collapse is halted by the magnetic field under idealized conditions. We take the value from Mouschovias & Spitzer (1976). The magnetic flux Φ_B is simply computed as the flux across the midplane, $B_0\cdot \pi R_C^2$, which in magnitude is the same flux that enters or leaves the spherical cloud.

2.2 Physics

The simulations solve for the equations of magnetohydrodynam-ics, with the addition of Ohmic resistivity as a nonideal effect, self-gravity, the gravitational force from the forming massive star, and radiation transport for the thermal emission from the gas and dust, on an axisymmetric grid. Next, we describe the set of equations solved.

The dynamics of the gas and dust are followed by numerically solving the following system of conservation laws, using the Pluto code (Mignone et al. 2007): $${\partial }_t\rho +\nabla \cdot \left(\rho \upsilon \right)=0,$$(4) $${\partial }_t\left(\rho \upsilon \right)+\nabla \cdot \left[\rho \upsilon \otimes \upsilon -\frac{1}{4\pi }B\otimes B+P_t|\right]=\rho a_{\text{ext}},$$(5) $${\partial }_{t}B+\nabla \times \left(cℰ\right)=0,$$(6) $$\begin{array}{l}{\partial }_t\left(E^K+E^{\text{th}}+E^B\right)\hfill \\ +\nabla \cdot \left[\left(E^K+E^{\text{th}}+P\right)\upsilon +cℰ\times B\right]=\rho \upsilon \cdot a_{\text{ext}},\hfill \end{array}$$(7)

which correspond to the continuity equation, momentum equation, induction equation and energy conservation equation, respectively. The gas is weakly ionized and has density ρ, velocity υ and magnetic field B. The magnetic field must also satisfy the solenoidality condition ∇ • B = 0. The total pressure $P_t=P+\frac{1}{8\pi }{B}^2$ is composed of the thermal and magnetic pressures. The electric field ℰ can be directly substituted by the expression $$cℰ=-\upsilon \times {\rm B}+\eta \nabla \times {\rm B},$$(8)

where η(ρ, T) is the Ohmic resistivity. The energy density of Eq. (7) is separated into its components: kinetic (E^K), thermal or internal (E^th) and magnetic (E^B). The acceleration source term a_ext is a sum of the gravity of the forming star a_g⋆ = −GM_⋆+/r² e_r, the self-gravity of the gas a_sg and the viscosity source term a_v.

The self-gravity acceleration source term is given by a_sg = −∇Φ_sg. The self-gravity potential Φ_sg is found by solving Pois-son’s equation ∇²Φ_sg = 4πGρ with a diffusion Ansatz, according to Kuiper et al. (2010). Self-gravity is inherently a three-dimensional effect, and a cloud with the characteristics we consider here is expected to produce an accretion disk with spiral arms. Since we perform the simulations on a two-dimensional axisymmetric grid, we mimic the missing angular momentum transport by the gravitational torques with the addition of a shear viscosity source term a_v = ∇Π/ρ, where the viscosity tensor Π is fixed with the α-parametrization of Shakura & Sunyaev (1973) and no bulk viscosity is considered. This approach was extensively studied in Kuiper et al. (2011), and we refer the reader to that article for more details on the computation of the viscosity tensor, as well as comparisons of the efficacy of angular momentum transport using this approach in comparison to nonviscous 3D models.

We treat radiation transport with the gray flux-limited diffusion approximation, by using the module Makemake described in ample detail in Kuiper et al. (2020). Specifically, we point the reader to Sects. 2.3.2 and 2.3.4 in that reference. In a nutshell, two equations are solved simultaneously for treating radiation transport. The first one is the zeroth-moment of the radiation transport equation within a given volume, which essentially says that the net emitted energy density (emission minus absorption) that is not lost through the boundary as a radiative flux, has to be stored as a local radiative energy density. This equation is simplified by assuming that the radiative flux can be written as a diffusion term in terms of the radiation energy density (hence flux-limited diffusion). The second equation describes the changes to E^th: the net absorbed energy density (absorption minus emission) has to be stored as internal energy in the gas. Both energy fields are solved without previous assumption of equilibrium between them; this is the so-called two-temperature approach Commerçon et al. (2011). In all but one of the simulations (see Sect. 5.2 for details), a constant value of the opacity for the dust and gas of 1 cm² g⁻¹ was used. The initial dust-to-gas mass ratio is set to 1%. The irradiation from the forming massive star is not considered in this study.

We use the resistivity model by Machida et al. (2007; which is based on a numerical study by Nakano et al. 2002) $$\eta =\frac{740}{X_e}{\left(\frac{T}{10\text{\hspace{0.17em}K}}\right)}^{1/2}\left[1-\tanh \left(\frac{n_{\text{H}}}{{10}^{15}{\text{cm}}^{-3}}\right)\right]{\text{cm}}^2{\text{s}}^{-1}\text{,}$$(9)

where we use for simplicity n_H = ρ/µ_H as the number density of hydrogen nuclei (µ_H being the molecular weight of hydrogen) and the ionization fraction X_e is given by $$X_e=5.7\times {10}^{-4}{\left(\frac{n_{\text{H}}}{{\text{cm}}^{-3}}\right)}^{-1}.$$(10)

At all moments in time, the properties of the formed massive star are computed using the evolutionary tracks for high-mass stars calculated by Hosokawa & Omukai (2009) and the instantaneous values of the stellar mass and accretion rate.

2.3 Boundary conditions

The computational domain assumes symmetry with respect to the rotation axis and equatorial symmetry with respect to the midplane. This means that scalar quantities and the velocity are simply reflected across the z-axis and the midplane. In the case of the magnetic field, it is reflected across the z-axis, but across the midplane, the parallel component switches sign instead of the normal component. The same symmetries imply that the boundaries across the azimuthal direction are periodic. Both the inner and outer boundaries (i.e., the sink cell and the outermost part of the cloud) impose a zero gradient condition for the magnetic field and the azimuthal and polar components of the velocity; for the radial component of the velocity and the density, only outflow but no inflow is allowed. Therefore, all matter that goes inside of the sink cell is considered as accreted.

2.4 Initial conditions: Parameter space and numerical configuration

We present an extensive parameter study of the system, with the express aim of aiding in the understanding of the physical processes that intervene in the formation and evolution of the disk and the magnetically driven outflows. In total, we present results from 31 different runs. Instead of describing all the parameters in multiple tables, we show a summary of the parameter space as a diagram in the right panel of Fig. 1.

We used five different axisymmetrical grids with increasing resolution, that we name ×1, ×2, ×4, ×8 and ×16. The grid for the radial coordinate r increases logarithmically with distance, starting from the radius of the inner boundary (3 au in the fiducial case), to the radius of the cloud core (0.1 pc). We often refer to the inner boundary as the sink cell. Using the assumed symmetries of the system, the colatitude θ extends from 0 to π/2 radians, and it is divided linearly into cells such that they have approximately the same dimensions in both the radial and polar directions for each distance to the center of the cloud. Table 1 details the number of cells in each direction, the minimum cell size and a reference cell size at 1000 au for all the grids used in this study, and assuming the sink cell size of the fiducial case.

The parameters with direct relation to physical quantities we considered are: the mass of the cloud core M_C, the resistivity model used, the normalized mass-to-flux ratio $\overline{\mu }$, the initial density profile exponent ß_ρ, the initial rotation profile exponent ß_Ω, and the ratio of the rotational-to-gravitational energy ζ. The values for each parameter in the fiducial case are highlighted in blue in Fig. 1. The fiducial case was run on all the grids. Modifications of the fiducial case are then investigated while keeping the rest of the parameters constant. Because of this, we refer in the paper to each simulation by the modified parameter followed by the grid it was run on; for example, M_C = 150 M_⊙[×8] refers to a cloud of 150 M_⊙, with the rest of the parameters as in the fiducial case for grid ×8. The full set of values with direct relation to physical quantities was run on the base grid ×4, and some selected values, which are underlined in Fig. 1, were run on grid ×8 as well. During the preparation of this study, we also ran a part of this physical parameter space on the lower resolution grids, however, we do not include the results here because they are redundant. Details of the choice of the parameters, as well as the comparison of results obtained, are presented in Sect. 5.

The numerical convergence of the results was studied by changing the values of the sink cell radius, the strength of the α shear viscosity and the Alfvén limiter using the grid ×2 as a base. A more detailed explanation of these latter parameters, as well as the results of the convergence study can be found in Appendix A.

Fig. 1 Initial conditions, overview of the evolution of the system, and summary of the parameter space presented in this study.

3 Overview of the evolution of the system

Figure 1 shows a schematic overview of the features observed during the evolution of the system, based on the results of the fiducial case but generally present in most of the simulations of the study. We describe these characteristic features and processes in the following subsections case-by-case.

3.1 Gravitational infall epoch

As soon as the simulation begins, the cloud starts collapsing under its own gravity. The flow quickly becomes super-Alfvénic (the kinetic flow velocity is larger than the Alfvén speed), and the magnetic field lines are dragged by the infall accordingly. The result of this collapse is that the configuration of magnetic field lines adopts an “hourglass” shape, as described in Galli et al. (2006), and that has been observed in Beltrán et al. (2019).

3.2 Protostar

A massive (proto)star is formed at the center of the cloud, which reaches a mass of ~20 M_⊙ after 30 kyr. The luminosity of the (proto)star is dominated by the conversion of gravitational energy into radiation during the first 20 kyr, and only after that time the stellar luminosity becomes comparable to the accretion luminosity. We point out that both the stellar and accretion luminosities are not taken into account in this study, but we will report on these radiative feedback effects in future work.

3.3 Magneto-centrifugal epoch

Due to the conservation of angular momentum, gas coming from large scales rotates faster as it reaches the center of the cloud. At t ~ 5 kyr, enough angular momentum is transported onto the center of the cloud to form a Keplerian-like accretion disk as a result. The disk grows in time due to the continued infall of material from large scales.

The accretion disk can be morphologically divided into a thin and a thick layer. The thin layer of the disk is a region of the mid-plane consisting of material in rotation with speeds close to the Keplerian speed. The density of the thin layer increases toward the center of the cloud 10⁻¹⁶−10⁻¹¹ g cm⁻³. Enclosing the thin layer, there is a less dense thick layer (with typical densities of 10⁻¹⁴ g cm⁻³) that is also rotating.

Roughly at the same time as the disk forms, magnetically driven outflows are launched, which we classify according to their speeds and driving mechanisms. The jet is launched by the magneto-centrifugal mechanism (in a way similar to the model of Blandford & Payne 1982) and it reaches typical speeds higher than 100 km s⁻¹. Rotation drags the magnetic field lines, as it has been observed in Beuther et al. (2020) and Girart et al. (2013). This creates a magnetic pressure gradient that drives a slower and broader tower flow (cf. the model of Lynden-Bell 2003). Its typical speeds are on the order of 10 km s⁻¹, and it broadens with time. A discussion of the dynamical processes present in the launching of outflows is offered in Paper II.

On the interface between inflow along the thick layer of the disk and outflow along the cavity, an additional structure is identified, which we refer as the cavity wall. Densities at the cavity wall are higher than those in the cavity, but also higher than the surrounding infalling envelope and the thick layer of the disk. In terms of its dynamics, material from the cavity wall contributes episodically to both the inflow and the outflow.

Fig. 2 Keplerianity on the midplane, computed with the fiducial simulation on grid ×8.

3.4 Magnetic braking epoch

As magnetic field lines are dragged by rotation, the resulting magnetic tension exerts a torque on the gas that brakes it, infalling as a consequence of the angular momentum loss. The inclusion of Ohmic dissipation makes magnetic braking negligible at early times in the simulation, but at around t ~ 15 kyr, magnetic braking starts to affect the innermost region of the disk and the cavity wall (r ≲ 30 au). As a result, a magnetically braked region is formed, where the material is mostly infalling but where the disk densities are kept because of constant replenishment through material from the cavity wall.

The extraction of angular momentum from the inner region of the cloud by magnetic braking modifies or even interrupts the magneto-centrifugal mechanism that drives the jet, but not the tower flow, which is mostly present at large scales. The tower flow broadens over time because the continuous dragging of magnetic field lines by rotation increases the magnetic pressure gradient on top of the disk, which is itself growing in time as well. This outflow broadening mechanism offers an explanation for the earliest outflow broadening observed and discussed in Beuther & Shepherd (2005). An analysis and discussion of the effects of magnetic braking in outflows, as well as their propagation over time into the large scales of the cloud is offered in Paper II, which focuses on the outflow physics of the same simulation data analyzed here.

4 Dynamics of the disk

After giving an overview on the general evolution of the system, we aim at understanding the dynamical processes present in the disk in detail by taking the fiducial case and examining selected terms of the equations of motion. With the identification of these processes, we are able to perform a more effective and comprehensive comparison of results for the full parameter space, which we present in Sect. 5.

4.1 Radial dynamics in the disk

The forming massive star accretes mass through the accretion disk, which implies that the velocity in the disk must have a component in the radially inward direction. The torque necessary for the angular momentum transport in the disk is provided by gravitational torques of forming spiral arms, which is modeled here through the use of the α shear viscosity model. Nevertheless, υ_ϕ ≫ υ_R in the disk, and so, it is useful to examine the effects of the different forces that act in the disk under the assumption of radial equilibrium.

In cylindrical coordinates (R, ϕ, z), the radial component of the MHD momentum equation, setting θ = π/2 (midplane), assuming axisymmetry, is given by $$\frac{\partial \left(\right(\rho {\upsilon }_R\left)\right)}{\partial t}+\nabla \cdot \left(\right(\rho {\upsilon }_R\upsilon -\frac{B_{R}B}{4\pi }\left)\right)+\frac{\partial P_t}{\partial R}=\rho g_R+\frac{\rho {\upsilon }_{\varphi }^2-B_{\varphi }^2/\left(\right(4\pi \left)\right)}{R},$$(11)

where the divergence operator ∇· acts according to $$\nabla \cdot f\left(\right(R,z\left)\right)=\frac{1}{R}\frac{\partial \left(\right(R{f}_R\left)\right)}{\partial R}+\frac{\partial f_z}{{\partial }_z}$$(12)

onto the vector fields f(R, z) ∊ {ρυ_Rυ, B_RB}; ɡ is the total gravitational force from the star and self-gravity, and P_t is the sum of the thermal and magnetic pressures.

In order to examine the forces that define radial equilibrium in the disk, we consider a parcel of gas moving in a circular orbit, for which we set υ_R = 0 in Eq. (11) and obtain $$\frac{\rho {\upsilon }_{\varphi }^2-B_{\varphi }^2/\left(\right(4\pi \left)\right)}{R}=-\rho g_R+\frac{\partial P}{\partial R}+\frac{\partial }{\partial R}\left(\right(\frac{B^2}{8\pi }\left)\right)-\nabla \cdot \left(\right(B_{R}B\left)\right).$$(13)

For B = 0, P = 0, and considering the midplane, one recovers the well-known condition for gravito-centrifugal equilibrium (on the corotating frame), which defines the Keplerian velocity, ${\upsilon }_K=\sqrt{GM\left(r\right)/R};M\left(r\right)$ being the enclosed mass (dominated by the mass of the protostar for small r). The use of the enclosed mass is an approximation to the use of the full potential given by the Poisson equation.

4.1.1 Keplerianity

The degree of gravito-centrifugal equilibrium in the whole disk can be measured by computing the ratio of the centrifugal force density to the gravitational force density: $$\frac{a_R^{\text{c}}}{\rho g_R}=\frac{\rho {\upsilon }_{\varphi }^2}{r\sin \theta }\cdot \frac{r^2}{\rho GM\left(\right(r\left)\right)\sin \theta }=\frac{{\upsilon }_{\varphi }^2}{{\upsilon }_K^2{\sin }^3\theta }.$$(14)

This motivates us to define the offset from Keplerianity as offset from Keplerianity $$\text{offset\hspace{0.17em}from\hspace{0.17em}Keplerianity\hspace{0.17em}=\hspace{0.17em}}\left|\right|\frac{{\upsilon }_{\varphi }}{{\upsilon }_K{\sin }^{3/2}\theta }\left|\right|-1,$$(15)

which in the midplane corresponds to the fractional difference of the azimuthal velocity with respect to the Keplerian value.

Figures 2 and 3 show the offset from Keplerianity presented as a percentage, in such a way that the negative values indicate sub-Keplerianity, and the positive values, super-Keplerianity. In the midplane, there is a Keplerian region within ±25% for all times, which corresponds to the thin accretion disk. The Keplerianity drops outside of the disk, which corresponds to the infalling envelope. At the exterior boundary of the disk, the mostly rotating material encounters the mostly-infalling material from the region of the envelope, creating a centrifugal barrier that is seen as a slight super-Keplerianity in Fig. 2. Even though we ignore the effect of viscosity in the analysis of this section, we note that, as it transports angular momentum outward, its effect is to increase positively the deviation from Keplerianity. The rest of the deviations from Keplerianity in the disk are discussed in the following sections.

Panels b and f of Fig. 3 contrast the values of the offset from Keplerianity in the midplane with values in planes parallel to it, which slice through the thick layer of the disk and reveal its dynamical structure. In general, the thick layer also shows Keplerianity within ±25%, but it decreases with altitude because the cylindrical-radial component of gravity also decreases with altitude. The regions close to the cavity become increasingly super-Keplerian, while the regions close to the infalling envelope become increasingly sub-Keplerian.

Fig. 3 Dynamics of the thick layer of the disk. The panels of the figure are grouped in columns such that they show: the density structure of the thick layer (a and e), the offset from Keplerianity at three different heights in the thick layer (b and f), the contributions to the specific energy (c and g), and the vertical balance of forces (d and h). Data from the fiducial simulation on grid ×8 was used for this figure.

Fig. 4 Contributions to the specific energy in the thin layer of the disk, calculated with the fiducial simulation on grid ×8.

4.1.2 Specific energies

We study the dynamics of the thin and thick layers of the disk by computing the terms of the Bernoulli equation, that is, the contributions to the specific energy of the material. This allows us to isolate the role of each term in the equations of magnetohy-drodynamics and compare its relative importance. The results of the energy calculations are available in Fig. 4 for the midplane, and Fig. 3 for the plane z = z₁. The contributions to the specific energy we considered are:

• the gravitational specific energy $$e^{\text{grav}}=\frac{GM\left(r\right)}{r},$$(16)

• the thermal specific energy $$e^{\text{th}}=\frac{P}{\rho \left(\text{Γ}-1\right)},$$(17)

• the contribution to the specific kinetic energy by the azimuthal component of velocity $$e_{\varphi }^K=\frac{{\upsilon }_{\varphi }^2}{2},$$(18)

• the contribution to the specific kinetic energy by the spherical-radial component of velocity $$e_r^K=\frac{{\upsilon }_r^2}{2},$$(19)

• the contribution to the specific magnetic energy by the toroidal component of the magnetic field $$e_{\text{tor}}^B=\frac{B_{\varphi }^2}{8\pi \rho },$$(20)

• the contribution to the specific magnetic energy by the poloidal component of the magnetic field $$e_{\text{pol}}^B=\frac{B_r^2+B_{\theta }^2}{8\pi \rho }.$$(21)

Initially, the specific gravitational energy increases as a function of r because it is given by the enclosed mass (see Fig. 4a). As the collapse progresses (Fig. 4b), the gravity of the forming massive star becomes more important, which means that the curve for e^grav has a point-source-like behavior close to the center of the cloud. At large scales, the enclosed mass of the envelope becomes important. As expected, the infalling envelope has $e_r^K>e_{\varphi }^K$, but the advection of angular momentum from large to small scales causes $e_{\varphi }^K$ to increase until $e_{\varphi }^K>e_r^K$, forming the disk. When the magnetic field lines are dragged by rotation in the disk, the toroidal contribution to the magnetic energy increases.

4.1.3 Thermal and magnetic pressure

The following two sections examine the contributions of the remaining terms of Eq. (13) to the equilibrium of cylindrical radial forces in the disk and their impact to the values of the deviation of Keplerianity shown in Fig. 2. First, we only consider the additional radial support that the thermal pressure gradient ∂P/∂R gives to the material in the disk, which be inferred by using Eq. (17) and Figs. 4 and 7. In agreement with the theory, we find that the thermal pressure gradient provides an additional small support against gravity and therefore makes material in radial equilibrium to appear to be slightly sub-Keplerian.

In Gaussian units, the magnetic pressure and magnetic energy density correspond to the same quantity, which means that we can use the specific energy density to study the magnetic pressure support in the disk. In the infalling envelope, the poloidal component of the magnetic field contributes the most to magnetic pressure (see, e.g., Fig. 4b), given that the initial magnetic field distribution is parallel to the rotation axis. When the disk is formed, the toroidal component of the magnetic field increases in general because the magnetic field lines are dragged by rotation. The toroidal magnetic field becomes dominant over the poloidal component, which additionally allows us to neglect the magnetic tension term ∇ • (B_RB) from Eq. (13) in most parts of the disk.

The magnetic pressure $P_B=\rho \left(e_{\text{tor}}^B+e_{\text{pol}}^B\right)$ decreases with distance in the infalling envelope (compare Figs. 4b and 7a), and overall in the disk as well, although the more complex structure of the magnetic field in the thin and thick layers of the disk means that there are regions of local increase with distance. Similarly to our argument with the thermal pressure gradient, υ_ϕ < υ_K for a negative magnetic pressure gradient. Therefore, the radius of the disk should be slightly larger compared to the nonmagnetic case, given that υ_ϕ decreases with distance and the angular momentum contained at larger R becomes sufficient to support a circular orbit. This argument is revisited in Sect. 5.3 after considering the effects of magnetic braking (Sect. 4.3) and magnetic diffusion (Sect. 4.2).

Additionally to magnetic pressure, the second term of the left-hand side of Eq. (13) reveals that the azimuthal component of the magnetic field reduces the effectiveness of the centripetal force required to keep a circular orbit. This can be measured by comparing the azimuthal component of the kinetic and magnetic energies: $$\frac{B_{\varphi }^2}{4\pi R}\cdot \frac{R}{\rho {\nu }_{\varphi }^2}=\frac{e_{\text{tor}}^B}{e_{\varphi }^K}.$$(22)

In the thin layer of the disk, the contribution of this term is negligible, but in the thick layer, they are comparable. In terms of Keplerianity, this term can cause the material to be slightly super-Keplerian; this can be observed when comparing panels f and g of Fig. 3.

4.2 Ohmic dissipation

According to the model of Ohmic resistivity by Machida et al. (2007), which we use, the resistivity increases progressively with density for the number densities we obtain here (n_H ≲ 10¹³cm⁻³). We computed the ratio of the magnitude of the induction and diffusion terms in the induction equation, which is analogous to the magnetic Reynolds number, as $${\mathrm{Re}}_m=\frac{\left|\right|\nabla \times \left(\right(\upsilon \times B\left)\right)\left|\right|}{\left|\right|\eta \nabla \times \left(\right(\nabla \times B\left)\right)\left|\right|},$$(23)

and show the results in Fig. 5, which are computed for two instants of time from the fiducial case. At t = 10 kyr, the plasma is fully diffusive in most regions of the thin layer of the disk; the thick layer is dominated by advection but it experiences magnetic diffusion in around one part per hundred; and the cavity wall is also partially diffusive. Later in time, at t = 22.75 kyr, only the highest density regions close to the massive protostar are dominated by diffusion, the thin layer experiences magnetic diffusion of around one part per hundred, and the thick layer is almost completely dominated by advection.

Ohmic resistivity acts mostly on the toroidal component of the magnetic field, as it can be seen by comparing the energy curves in Figs. 4 and 3. In the thick layer of the disk at late times (Fig. 3g) $e_{\text{tor}}^B$ roughly follows $e_{\varphi }^K$. This sets the picture for the case where magnetic diffusion is low. On the other hand, in the thin layer of the disk, at early times (Fig. 4b), $e_{\text{tor}}^B$ and $e_{\varphi }^K$ are clearly decoupled: the toroidal magnetic field component becomes low inside of the disk, corresponding to the region where diffusion dominates over advection. In the outer disk, where the density is lower, the resistivity is lower as well, and the toroidal component of the magnetic field is allowed to increase.

Fig. 5 Magnetic Reynolds number, computed with the fiducial simulation on grid ×8. For t = 10 kyr, the radius of the disk is around 80 au, while for t = 22.75 kyr, the exterior radius is ≈700 au and the interior radius is ≈20 au.

4.3 Magnetic braking

The momentum equation in the azimuthal direction, considering axisymmetry and circular motion, reduces to $$\frac{\partial \left(\right(\rho {\upsilon }_{\varphi }\left)\right)}{\partial t}=\frac{1}{R^2}\frac{\partial }{\partial R}\left(\right(R^2{B}_{\varphi }{B}_R\left)\right)+\frac{\partial }{\partial z}\left(\right(B_{\varphi }{B}_z\left)\right).$$(24)

The right-hand side of this equation is equal to zero if there are no magnetic fields, and it leads to the conservation of angular momentum. The terms on the right-hand side constitute the azimuthal component of the magnetic tension force, which exerts a torque that reduces angular momentum locally. This process is known as magnetic braking, and it can be understood as a consequence of the dragging of magnetic field lines by rotation (Galli et al. 2006).

Even though Ohmic resistivity reduces the toroidal component of the magnetic field (cf. $e_{t\text{or}}^B$. in Fig. 4b), it does not completely suppress it. Over time, the magnetic field lines are wound enough so that magnetic braking happens in the innermost parts of the disk. In Fig. 4c, which corresponds to t = 22.75 kyr, $e_{t\text{or}}^B$. is higher in the disk than when it is measured at t = 10 kyr, while $e_{\varphi }^K$ decreases over the same interval of time in the inner 20 au; those quantities show the reduction of angular momentum due to magnetic tension. In the face of the reduction of angular momentum, the material falls toward the protostar, as evidenced by the increase in $e_r^K$; this creates the magnetically braked region. The cavity wall and the thick layer of the disk are also affected by magnetic braking; material from several directions is then delivered to the magnetically braked region.

4.4 Vertical support

In the vertical direction (i.e., parallel to the rotation axis), the momentum equation reads $$\frac{\partial \left(\right(\rho {\upsilon }_z\left)\right)}{\partial t}+\nabla .\left(\right(\rho {\upsilon }_z\upsilon -\frac{B_{z}B}{4\pi }\left)\right)+\frac{\partial P_t}{\partial z}\rho {ɡ}_z$$(25)

with the same definition of the divergence as given in Eq. (12). Considering a parcel of gas moving in a circular orbit (υ_z = 0), and given that in the disk region B_z ≪ B_ϕ, Eq. (25) becomes $$\frac{\partial P}{\partial z}+\frac{\partial }{\partial z}\left(\right(\frac{B^2}{8\pi }\left)\right)=\rho {ɡ}_z.$$(26)

The ratio of the thermal to magnetic pressure is equivalent to the ratio of the specific thermal energy to the total magnetic energy. In the thin layer of the disk (Fig. 4), thermal pressure dominates over magnetic pressure, from which we conclude that the thin layer is supported mainly by thermal pressure (at late times, magnetic pressure becomes comparable but not much higher than thermal pressure). On the contrary, the thick layer of the disk is supported by magnetic pressure, as evidenced by Figs. 3d and h, where we present a direct comparison of both pressure gradients to gravity in a vertical slice that crosses the disk. As the vertical component of gravity vanishes at the mid-plane, both curves increase toward z = 0. In the thick layer of the disk, thermal pressure decreases, while the magnetic pressure gradient compensates gravity. At the cavity wall, thermal pressure increases because it has a higher density than the thick layer of the disk. Finally, in the cavity, magnetic pressure dominates over gravity in the vertical direction, which is expected in the case of a magnetically driven outflow.

Fig. 6 Distribution of the magnetic field strength as a function of density. This histogram was computed with the fiducial simulation run on grid ×16, for time t = 7 kyr. The color scale indicates the polar angle: the dark purple colors correspond to the midplane (with the disk corresponding to ρ ≳ 10−15 g cm−3) and the light blue colors correspond to the rotation axis (where the jet is located).

4.5 Effects of other forms of magnetic diffusion

In this study, the effects of ambipolar diffusion and the Hall effect are not considered, and this formally constitutes a caveat of our approach which we plan to address in a future study. In this section, we discuss the expected effects that other forms of magnetic diffusivity might have on the results reported in the previous sections, with a focus on the fiducial case. We examine every region of interest for the discussion: the accretion disk, the infalling envelope and the outflow cavity. The simulation catalog also includes an investigation of the effects of the resistivity model; the results can be used to further assess the missing magnetic diffusivity and are presented in Sect. 5.2.

First, we consider only the accretion disk. The curves for ambipolar diffusivity and Ohmic resistivity given in, for example, Marchand et al. (2022, 2016) and Tsukamoto et al. (2021) reveal that, at least for the range of number densities 10⁸ ≲ n_H ≲ 10¹³ cm⁻³, which are found in the thin layer of the accretion disk, the values of the diffusivities are similar within the uncertainties of the dust models considered for their derivation, differing by about a factor of ten. The distribution of the magnetic field strength with density is presented in Fig. 6, where we use a color scale that highlights with darker colors the midplane. For the densities of the accretion disk (ρ > 10⁻¹⁴ g cm⁻³), the magnetic field strength increases with density, until it reaches between 0.1 and 1G. When comparing the magnetic field distribution to the results of a recent study by Commerçon et al. (2022), who considered the effects of ambipolar diffusion but not Ohmic dissipation, we find consistent results for the thin layer of the disk. In the thick layer, however, the addition of higher, more realistic magnetic diffusivities might reduce the vertical magnetic pressure support. This means that thinner accretion disks than those reported here would be expected in a more realistic setting.

As for the effects of magnetic braking, we expect that ambipolar diffusion and the Hall effect reduce the magnetic Reynolds number, that is, the material would become more diffusive. Contrary to Ohmic dissipation, ambipolar and Hall diffusivities depend on magnetic field, and so the reduction could be of particular importance at later times, when considerable magnetic flux has been transferred to the central regions of the cloud core. This would cause a delay in the formation of the magnetically braked region in the inner disk (see also the discussion in Sect. 5.2). On the other hand, magnetic braking could play an important role in determining the total angular momentum transferred to the forming massive star, as well as the growth of HII regions by terminating the magnetically driven outflows (Martini et al., in prep.). An investigation of the role of all forms of magnetic diffusivity in the long-term evolution of the accretion disk and the jet is necessary. However, at least for early times (t ≲ 20 kyr) in the formation of the massive star, our results fit observational constraints (see Sect. 6).

As a second region of interest, we consider the infalling envelope. According to the curves by Marchand et al. (2022), ambipolar diffusivities tend to be higher at lower densities (n_H < 10⁸ cm⁻³). This suggests that we are missing key magnetic diffusion during the early phases of gravitational collapse, where densities are low. However, there are theoretical and observational results that suggest that the collapse phase is still reasonably modeled in our simulations. In this simulation series, we consider the case of weak magnetic fields (see the discussion on this assumption in Sect. 5.3). Because ambipolar diffusivities increase with magnetic field strength, we expect magnetic diffusion to play a stronger role in cases when the magnetic fields are initially strong. In the studies by Commerçon et al. (2022) and Mignon-Risse et al. (2021), who considered $\overline{\mu }=5$ and ambipolar diffusion but no Ohmic dissipation, the magnetic field distribution for low densities coincides well with what is depicted in Fig. 6. Based on this, we would not expect a strong qualitative divergence of our results for the gravitational collapse epoch (t ≲ 5 kyr) upon the addition of ambipolar diffusion to our models. From the observational point of view, the role of ambipolar diffusion as a necessary enabler for gravitational collapse in a magnetized cloud core was recently questioned in Ching et al. (2022). These latter authors found a supercritical low-mass prestellar core, which is not yet self-gravitating. This means that the magnetic field was diffused before the gravitational collapse takes place, that is, earlier than predicted by the classical picture of ambipolar diffusion turning a subcritical core supercritical and in a consistent way with our simplification.

Lastly, we take a look at the magnetic diffusivities in the outflow cavity, which is the focus of Paper II. Observations of protostellar jets have found evidence of ionized material (see e.g., Moscadelli et al. 2021; Rodríguez-Kamenetzky et al. 2017; Carrasco-González et al. 2021; Guzmán et al. 2010), possibly from shock ionization. As discussed in Marchand et al. (2016, 2022), ionization leads to a general decrease in the magnetic dif-fusivities. So far, the models of ambipolar diffusion and the Hall effect done in, for example, Marchand et al. (2016, 2022) and Tsukamoto et al. (2021) do not consider the shock ionization in the outflow cavity. Material in the outflow cavity and close to the massive protostar is very low density and is in presence of strong magnetic fields, which means that neglecting shock ion-ization would overestimate magnetic diffusivities coming from the ambipolar diffusion and Hall effect terms, when comparing to the curves by Marchand et al. (2022). By only considering Ohmic dissipation (for which the resistivity is low for low densities and independent of magnetic field strength), the material in the cavity behaves closer to ideal MHD theory, in agreement to what has been found from observational evidence (see Sect. 6 and Moscadelli et al. 2022). This can be seen in Fig. 6, where the magnetic field strengths in the cavity (light-colored dots) are able to increase while keeping the values constrained in the disk (dark-colored dots). A realistic treatment of the problem, however, would require the consideration of all forms of magnetic diffusivity (which would probably impact the disk dynamics) and the effect that thermal ionization and shock ionization have on the magnetic diffusivities corresponding to the outflow cavity.

Fig. 7 Density (panel a) and temperature (panel b) profiles on the mid-plane, for the fiducial simulation on grid ×8.

Fig. 8 Density (panel a) and temperature (panel b) for two different altitudes in the thick layers of the disk (blue lines), and the corresponding profile for the thin layer (black). The data corresponds to the fiducial simulation on grid ×8.

4.6 Density and temperature profiles

Figure 7 presents the density and temperature profiles for the midplane. The thin layer has densities over 10⁻¹⁶g cm⁻³, and they reach over 10⁻¹¹g cm⁻³ at the inner boundary. For t ≲ 15 kyr, the density profile of the thin layer follows an approximate power law. Later, the magnetically braked region is formed, corresponding to the high density feature observed forr ≲ 30au in Fig. 7a. Some of the material of the magnetically braked region is drawn from the thin layer of the disk: the density profile of the thin layer reveals the existence of a plateau close to its magnetic braking radius, which coincides as well with the decrease in the azimuthal contribution to the specific kinetic energy due to magnetic braking (Fig. 4c). However, another part of the material in the magnetically braked region corresponds to the inflows from the cavity wall and thick layer of the disk, both of which also experience magnetic braking and continuously replenish it.

The temperature profile of the thin layer follows roughly a power law of the form T ∝ r^−1/2, and it experiences a general increase over time. The transition between the disk and the envelope is accompanied by a change in the temperature gradient; the transition between the magnetically braked region and the thin layer of the disk is subtler. We remind the reader that the present simulations do not include the effect of irradiation from the star, for which we anticipate two additional features not present in Fig. 7b, namely, an increase in temperature due to the flux from the forming massive star, and the existence of the dust evaporation front (see Kuiper et al. 2010). Those features are only relevant for t ≳ 20 kyr.

Figure 8a presents the density profile at two different altitudes crossing the thick layer of the disk at the time t = 10 kyr (when M_⋆ = 3 M_⊙), accompanied by a comparison with the density in the midplane. Differently than the thin layer of the disk, the thick layer has a nearly uniform density of ~10⁻¹⁴g cm⁻³. The density drop for r ≲ 10 au corresponds to the outflow cavity. The mean density in the thick layer decreases slightly with time, however, it roughly remains within the same order of magnitude throughout the simulation.

Because the density is nearly uniform, the temperature of the thick layer of the disk remains also relatively uniform, at around 200 K in comparison to the power-law temperature of the thin layer. This means that it should be possible in principle to observe the vertical stratification of the disk by distinguishing the line emission from each layer.

5 Dependence of the resulting disk properties on initial cloud properties

We perform a series of simulations to explore the large parameter space of initial conditions. We vary the fiducial model model described above in terms of the initial mass of the cloud core, its density and rotation profiles, its rotational-to-gravitational energy ratio, its initial magnetic field strength and probe the effect of different models of Ohmic resistivity. In this section, we discuss how those changes affect the processes discussed in Sect. 4. As a probe, we use the radial extent of the Keplerian-like disk (between the magnetic braking radius and the [outer] radius), defined as the region where the material is within ±25% Keplerianity, and the results of which are presented in Fig. 9 for the whole parameter space. Similarly, the mass of the forming massive star is presented in Fig. 10 for the whole parameter space. We also examine differences in the density structure of the disk for selected values of the parameters, and relate them to the dynamical information revealed by the specific energy; the results are shown in Fig. 11.

5.1 Dependence on cloud core mass

Cloud cores of masses of 50, 100, 150 and 200 M_⊙ were considered, which form stars ranging masses from 6 to 75 M_⊙ after 30 kyr of evolution (see Fig. 10a). Further evolution and additional physical effects which are not considered here would be necessary to determine the final mass of the star. When changing the initial core mass of the fiducial model, we also change the normalization of the rotation profile in order to keep the rotational-to-gravitational energy ratio as in the fiducial model. The mass and size of the cloud core determines its free-fall timescale $$t_{\text{ff}}={\left(\right(\frac{{\pi }^2{R}_c^3}{8G{M}_C}\left)\right)}^{1/2},$$(27)

which is 52.4 kyr for the fiducial case of M_C = 100 M_⊙ and RC = 0.1 pc. Panels A1 and A2 of Fig. 9 show that the disk and the magnetically braked region are smaller for the lower-mass cases when plotted against time, however, when investigated as a function of the stellar mass as a fraction of the initial mass of the cloud (panels A3 and A4), both the size of the disk and the magnetic braking radius coincide well. We find that resulting curves roughly follow the empirical linear relation $$R_{\text{disk}}\left(\right(m\left)\right)\approx R_{\text{onset}}+R_{\infty }\left(\right(m-m_{\text{onset}}\left)\right),$$(28)

where we define $$m\equiv \frac{M_{\star }}{M_C}.$$(29)

We note that R_onset = 3 au is the minimum radius of the disk we can detect in the simulation data (the size of the sink cell), m_onset = 0.016 is the corresponding normalized mass of the protostar for R_onset, and R_∞ = 6380 au can be interpreted as the approximate value of the radius of the disk at the hypothetical end of the collapse where m ~ 1. This latter value will never be reached in reality because of stellar feedback: Kuiper & Hosokawa (2018), for example, found that the disk lifetime is expected to be of only a few 10⁵ yr (see also Kee & Kuiper 2019 for the conditions for limiting accretion through the disk). A similar empirical disk size formula was previously offered in Kuiper et al. (2011).

We also verified that the mass of the protostar, expressed as a fraction of the mass of the cloud, scales with time presented as fraction of the free-fall timescale (cf. panels a and b of Fig. 10). This result indicates that despite the differences in the scale-dependent microphysics, namely, the radiative transfer, resistivity and self-gravity, the disk scales well with the free-fall timescale, which allows our results to be used with a more ample range of cloud masses than what is presented here. Defining τ ≡ t/t_ff, we also find an approximate fit for the mass of the protostar with a power law: $$m\left(\right(\tau \left)\right)\approx k_1{\tau }^{k_2},$$(30)

where k₁ ≈ 0.56 and k₂ ≈ 1.8.

Fig. 9 Radius of the disk (Rdisk) and its magnetic braking radius (Rmb) for different initial values of the mass of the cloud, magnetic field, density profile, rotational profile, rotational energy and resistivity models. The transparent lines in the panels that show Rmb indicate the full data, while the solid lines show the moving average. The colored dots in panel A4 indicate t = 30 kyr in each simulation. For all the panels in rows B and C, MC = 100 M⊙ was used.

5.2 Ohmic resistivity

The Ohmic resistivity model controls the diffusivity of the plasma, and therefore it has an impact on magnetic braking. We studied the effects of using the resistivity model of Machida et al. (2007) in two versions: the full expression dependent on both density and temperature, and the same expression but with a fixed temperature of T = 10 K, as it was used in the isothermal study of Kölligan & Kuiper (2018). We denote the first one as η(ρ, T) and the second one as η(ρ) given the near linear behavior of the resistivity in the range of densities for which we are interested. We estimate that the version that considers the temperature dependence differs at most by a factor of 100 in regions close to the protostar with respect to the version with the fixed temperature. We compare our results as well with the ideal MHD case (η = 0) and the nonmagnetic (hydrodynami-cal) case (η → ∞). In order to mimic the magnetic diffusivity expected in the disk region if ambipolar diffusion were taking into account, we also include a simulation with an artificially high Ohmic resistivity (ten times the temperature-independent value, denoted as 10η(ρ)). This factor is motivated by the difference in the curve for ambipolar diffusivity and Ohmic resistivity in Marchand et al. (2022), for disk densities. There is an additional difference between the simulations of this series: while most simulations of this study consider a constant dust and gas opacity of 1 cm² g⁻¹, the simulation with η(ρ, T) uses the dust opacity table from Laor & Draine (1993) and a constant gas opacity of 10⁻² cm² g⁻¹.

For the rest of Sect. 5, we present the size of the disk as a function of the stellar mass as a fraction of the cloud mass, instead of time, because the calculation of the offset from Keplerianity (which defines the disk) is dependent on the gravity of the (proto)star, and by fixing it we are able to compare and assess the effects of the other terms that influence the radial equilibrium of the disk, namely, the angular momentum and pressure support. Additionally, our results can be scaled as a result of the discussion in Sect. 5.1.

In the panel B1 of Fig. 9, we observe that the magnetically braked region develops earlier (when the mass of the star is lower, cf. Fig. 10c) in the ideal MHD case, as compared to the corresponding curves for the simulations that consider resistivity. The case with η(ρ) requires a more massive protostar (more time evolution) for the magnetic braking radius to appear, a trend that is confirmed with the simulations that consider higher resistivity: the temperature-dependent resistivity delays the appearance of magnetic braking, while the simulation with the artificial factor of ten in resistivity delays it even further. As expected, the disk in the case with no magnetic fields does not develop a magnetic braking radius. A disk is able to form at all in the ideal MHD case because the relatively weak magnetic field considered for this parameter scan.

From this evidence, we infer that the presence of resistivity delays the action of magnetic braking compared to the results only considering ideal MHD. Higher resistivities delay magnetic braking more, however, as discussed in Sect. 4.3, because magnetic diffusion does not completely suppress the toroidal magnetic field, the Lorentz force continues to grow until enough magnetic tension force is built to decelerate the material. As a result, the innermost parts of the disk lose their centrifugal support.

The disk radius (Fig. 9B2) is in general larger in the magnetic cases in comparison to the nonmagnetic case. Moreover, the disk radius of the run with no resistivity (ideal MHD) is marginally larger than the ones with resistivity. Magnetic braking removes angular momentum from the disk because of the torque produced by magnetic tension; therefore, a smaller disk would be expectable. However, magnetic tension is highest in the inner disk, where the magnetic field lines are wound the most, and not in the outer disk. As discussed in Sect. 4.1.3, magnetic pressure can provide an additional radial support for the disk. This means that in the ideal MHD case, where the toroidal magnetic field is high in the disk, the additional magnetic pressure supports a larger disk. Magnetic diffusion lowers the toroidal magnetic field, and therefore the disk does not have radial magnetic pressure support: this can be seen at early times (where M_⋆ < 5 M_⊙ or m = 0.05) where the simulations with η(ρ, T) and 10η(ρ) yield a disk radius that almost matches the one obtained in the nonmagnetic case. As the toroidal magnetic field increases over time due to the winding of magnetic field lines, the magnetic pressure support also increases and the disk becomes larger than the nonmagnetic case.

Given that magnetic braking delivers mass onto the protostar, we investigate the role of resistivity in the mass of the forming massive star. The corresponding curves for the mass of the protostar as a function of time are shown in Fig. 10c. The results for the resistive and nonresistive cases are very similar. The stellar mass, however, increases more in the nonmagnetic case in comparison to the magnetic cases. This happens because the protostar in the nonmagnetic case accretes material not only through the disk, but also through the infalling envelope along the bipolar directions. On the contrary, in the magnetic cases, there are magnetic outflows that impede accretion through the axis, resulting in a slightly lower mass.

Fig. 10 Mass of the protostar, corresponding to the mass in the sink cell, for several values of the (a–b) mass of the cloud, (c) resistivity model, (d) mass-to-flux ratio, (e) density profile, and (f) rotation profile and rotational energy. Panel b shows the same results as panel a, but normalized in such a way that scalability can be readily seen. For panels (c)–(f), MC = 100 M⊙ and tff = 53.73 kyr.

Fig. 11 Density structure (Col. 1) and specific energies (Col. 2) of the disk for selected values of the normalized values of the mass-to-flux ratio (rows B and C), and rotation profile (Col. D). All the snapshots were taken at t = 15 kyr of evolution. All cases consider a rotational-to-gravitational energy of 4%.

5.3 Initial magnetic field strength

We study values of the normalized mass-to-flux ratio ranging from $\overline{\mu }=2$ (strongest magnetic field) to $\overline{\mu }=40$ (weakest magnetic field). The results of the radial extent of the centrifugally supported disk are shown in the panels B1 and B2 of Fig. 9. The curve for $\overline{\mu }=2$ goes quickly to zero, which means that the structure formed is strongly sub-Keplerian and it is not classified as a Keplerian-like disk according to our criterion. A more detailed study of the specific energies revealed that the structure has a comparable radial and azimuthal velocities. Our result of no disk for $\overline{\mu }=2$ should be taken with caution, because we do not include ambipolar diffusion in our calculations. As ambipolar diffusivities increase with magnetic field (contrary to Ohmic dissipation), the cases with strong magnetic fields are missing magnetic diffusivity that may allow the disk to form (see the discussions in Sects. 4.5 and 7.3).

The other values of the mass-to-flux ratio confirm that magnetic braking is responsible for the formation of the magnetically braked region because higher initial magnetic fields build enough magnetic tension earlier. Additionally, the magnetically braked region is larger for stronger magnetic fields. The radius of the disk, on the other hand, is larger for stronger magnetic fields, which supports the hypothesis that higher magnetic pressure provides additional support against gravity and allows material from larger radii (with lower angular momentum) to be part of the disk. As for the accretion onto the protostar, Fig. 10d reveals that higher magnetic fields do tend to deliver more mass onto the protostar, due to stronger magnetic braking at late stages.

The rows A–C of Fig. 11 present a comparison of the morphological differences of the disk when it has evolved for t = 15 kyr under different values of the mass-to-flux ratio. For $\overline{\mu }=40$, the weaker magnetic field produces a disk with a thick layer that is smaller than the thin layer. Both regions of the disk are of the same size for $\overline{\mu }=20$, which confirms magnetic pressure as the nature of the vertical support in the thick disk. However, the part of the thin layer that is not enveloped by the thick layer of the disk becomes inflated by magnetic pressure, as the curve for the toroidal contribution to the specific magnetic energy reveals in Fig. 11B2: in the outer disk (r ~ 100 au), $e_{t\text{or}}^B$ becomes higher than e^th. Nevertheless, a similar analysis of the specific energies at later times reveals that the thick layer of the disk grows over time until it completely envelops the thin layer (which happens at around t ~ 19 kyr). At t ≈ 25.5 kyr, the thin and thick disk layers of the disk of the case with $\overline{\mu }=40$ become morphologically and dynamically similar to the results shown in panels A1 and A2 corresponding to the $\overline{\mu }=20$ at 15 kyr; the difference being the size of the disk and the mass of the protostar. For a stronger initial magnetic field ($\overline{\mu }=5$, row C of Fig. 11), the thick layer becomes thicker because of the increased magnetic pressure.

5.4 Initial density profile

A steeper initial density profile accelerates the accretion onto the protostar at early times because of the increased concentration of mass (and therefore gravity) in the center of the cloud. Apart from a density profile with ß_ρ = −1.5, we also tried a profile with ß_ρ = −2. The mass of the protostar (Fig. 10e) increases rapidly at early times with ß_ρ = −2, but the accretion rate decelerates over time, which is the opposite behavior observed for ß_ρ = −1.5. This behavior can be readily seen from the power-law fits to those curves: the normalization constant is similar in both cases, but while the fiducial case is described by m ∝ τ^1.80, the steeper initial density profile exhibits m ∝ τ^0.94. The radius of the disk (Fig. 9C2) is larger in the case of the shallower density profile when seen as a function of the stellar mass and its behavior is longer linear. However, given that the stellar mass increases more rapidly in the case of the steeper density profile, in this latter case, a disk of a given size is observed earlier in time. A similar comparison of stellar masses and time reveals that magnetic braking happens at roughly the same time for both simulations, because the magnetic field is wound by rotation in a similar way in both cases. While the magnetically braked region appears at M ≈ 3 M_⊙ when ß_ρ = −1.5, it appears at M ≈ 17 M_⊙ for ß_ρ = −2; however, the protostar reaches both masses at the same time, t ≈ 0.2t_ff ≈ 10.7 kyr. As a result, we conclude that a disk that does not exhibit a magnetically braked region at a given time can be obtained for a range of protostellar masses by adjusting the initial density profile.

5.5 Initial rotation profile

We investigated the effect of two scenarios: initial solid-body rotation and a rotation profile that scales with the cylindrical radius at Ω ∝ R^−ßΩ, with ß_Ω = ß_ρ/2 = −0.75. This choice of ß_Ω keeps the initial ratio of rotational to gravitational energy independent of the radius of the cloud. Additionally, the particular value of −0.75 produces a cloud for which the offset from Kep-lerianity is uniformly negative. We complemented this parameter scan with an investigation of the effects of using a lower initial rotational-to-gravitational energy ratio ζ of 1% instead of 4% as in the fiducial case.

An inspection between panels D1 and A1 of Fig. 11 reveals that the thick layer of the disk is flatter in the cloud with the steep rotation profile. The higher angular velocity at the center of the cloud allows for the formation of a larger accretion disk earlier in the simulation, in comparison with the cloud initially rotating as a solid-body. Panel C4 of Fig. 9 confirms this: the curves for the disk radius for ß_Ω = −0.75, ζ = {0.01,0.04} are in general higher than the corresponding curves for ß_Ω = 0, and resemble a power law rather than a linear function as it is the case for solid-body rotation. For a fixed value of ß_Ω, the disk becomes larger with more content of initial rotational energy. Additionally, both curves for ß_Ω = −0.75 show that the disk forms at very early times during the simulation compared to ß_Ω = 0: for the solid-body case, angular momentum conservation during the gravitational collapse provides the necessary increase of angular momentum in the center of the cloud in order to build up a disk, while the steep rotation profile already provides some of this angular momentum from the start. Interestingly, the cases ß_Ω = −0.75, ζ = 0.01 and ß_Ω = 0, ζ = 0.04 produce a similar disk radius. The mass of the formed protostar (Fig. 10f) is lower in the cases with a steep rotation profile and higher initial rotational energy. This is consequence of the earlier formation of a disk that is also larger (and therefore more massive), at it reduces the efficiency of accretion onto the protostar (in Oliva & Kuiper 2020, we utilized this fact to study disk fragmentation in the absence of magnetic fields). The magnetic braking radius of the disk is very similar for all cases when plotted as a function of the protostellar mass, however, given that the masses grow different in time, magnetic braking appears later in the case with ß_Ω = 0, ζ = 0.04.

From the dynamical point of view, a comparison between panels D2 and A2 of Fig. 11 reveals that in general, the different contributions to the specific energy have a more uniform radial gradient. Given that the offset from Keplerianity in the midplane can also be computed as $1-2e_{\varphi }^K/e^{\text{grav}}$, the disk formed with the steep rotation profile is more uniformly Keplerian. This is also manifested when comparing the curves for $e_r^K$ inside of the disk in both cases.

6 Comparison to observations

6.1 Natal environment of massive star formation

The parameter space in our simulation set was motivated by observational typical values found in regions of massive star formation. Surveys of cloud cores in the early stages of massive star formation have found power-law density profiles with exponents ranging in 1.5 ≲ −ß_ρ ≲ 2.6 (see, for example, Gieser et al. 2022; Beuther et al. 2018, 2002; Mueller et al. 2002; Hatchell & van der Tak 2003), with fragmenting cores having typical values of around ß_ρ ≈ −1.5. We take values in the extremes of this interval.

For the initial angular momentum, we start from the results of the survey by (Goodman et al. 1993), who found that the ratio of rotational to gravitational energy of cloud cores (ζ) is of a few per cent, with typical values of around 2%. Those authors fit linear gradients to the cloud cores that had signs of rotation, and found solid-body profiles. Our second choice for the initial rotation profile, ß_Ω = ß_ρ/2, is motivated by keeping ζ independent of the radius of the cloud core, as explained in Sect. 5.5.

The values of the mass-to-flux ratio are motivated by the supercriticality in star-forming cloud cores (e.g., Beuther et al. 2020; Girart et al. 2013) and the following analysis. In low-mass prestellar cores, the magnitudes of the magnetic field have been found to be on the order of a few to tens of micro Gauss (Ching et al. 2022; Troland & Crutcher 2008). As the precol-lapse conditions in the massive star formation case are currently not known, we consider two scenarios (see also the discussion in Machida & Hosokawa 2020) as follows. The first scenario assumes that the normalized mass-to-flux ratio is independent of mass, which means that magnetic field strengths of 100−1000 µG are required in the cloud prior to the gravitational collapse. On the other hand, the second scenario considers that the magnetic field strength of the high-mass star formation case is similar to the low-mass counterpart, which implies that the mass-to-flux ratio must be high.

The first case is investigated by using the low mass-to-flux ratios in our parameter space $\left(\overline{\mu }\sim 5\right)$, while the second case corresponds to high mass-to-flux ratios, including our fiducial case $\left(\overline{\mu }=20\right)$.

6.2 Size of the disk

The analysis in Sect. 5 has shown that the radius of the disk is more strongly determined by the initial density and rotation (angular momentum) profiles of the large scale initial condition, as compared to the magnetic field and resistivity models, within the ranges of those values expected from observations. Panel B4 of Fig. 9 shows that when M_⋆/M_C ~ 0.15, the average R_disk for the values of $\overline{\mu }$ we considered is around 850 au, with a variation on the order of ±200 au. However, the radius of the disk and its growth rate vary much more with the density profile and the rotation profile (panels C2 and C4 of Fig. 9). Smaller disks are produced with a lower and flatter initial distribution of angular momentum, and with shallower initial density profiles. This conclusion is crucial for the comparison of our results with previous studies, which is done in Sect. 7.

With regards to the available observational evidence of diskjet systems around forming massive stars, we see that in principle it is possible to obtain a variety of disk sizes that are not only dependent on the age of the system but also on the initial distribution of matter and angular momentum, and not so strongly determined by magnetic braking in the outer regions of the disk. For example, assuming that the molecular cloud has the typical values of M_C = 100 M_⊙ in a sphere of radius 0.1 pc of our fiducial case, we obtain a disk of radius ~ 1200 au when the mass of the protostar is M_⋆ ~ 10 M_⊙, if we assume a steep and fast initial rotation profile (ß_Ω = −0.75, ζ = 0.04), which would provide a possible initial configuration for the observations of HH 80–81 (Girart et al. 2017).

In Moscadelli et al. (2022), we utilized the fiducial simulation run on grid ×16 to model the accretion disk and jet of IRAS 21078+5211. For that particular system, observations of molecular rotational transitions had revealed the existence of a Keplerian-like accretion disk of around 200 au in size (Moscadelli et al. 2021). In order to find a model of our catalog that was consistent with the observations, in Moscadelli et al. (2022) we used the mass of the protostar, the radius of the disk, and the morphology of the ejected material to constrain the parameter space and the time elapsed. In that study, we found strong agreement between the velocity field of the observed water masers and the velocity streamlines predicted by the chosen model from our catalog. This good agreement in turn constitutes evidence that the results for the size of the disk as a function of the precollapse conditions – which we present in Sect. 5 – are in line with what is expected from observations. We note, however, that the wide streamline traced by the maser points (Fig. 2 of Moscadelli et al. 2022) is slightly closer to the assumed disk plane compared to the reference wide streamline from the simulations. This can be caused by a number of factors, for example, a slightly different assumed perspective between the simulations and the observations, but a possible explanation is that the disk in IRAS 21078+5211 is thinner than what is predicted in the simulations.

7 Comparison with previous numerical studies

The present study is a continuation of the work started by Kölligan & Kuiper (2018). Their simulations correspond to our fiducial case but with a grid equivalent to our ×2 grid, and used an isothermal equation of state. An isothermal equation of state yields a much smaller pressure scale height of the disk or thinner disk, respectively. In combination with the rather coarse grid resolution, the thickness of the disk could not be resolved on the numerical grid, making an in-depth convergence study problematic. In that sense, the addition of radiation transport and higher resolutions has enabled us to resolve the thin layer of the disk, and clearly differentiate it from the surrounding region that we call the thick layer of the disk. Also, they do not report about the existence of the magnetically braked region; however, it is visible upon close examination of their Fig. 21 as a region where the vertical velocity is negative close to the center of the disk. This was originally thought to be a numerical effect caused by the additional mass created by the Alfvén limiter, however, we performed a parameter scan with several values of the Alfvén limiter, the highest values of which produce negligible artificial mass, and found a magnetically braked region in all of them.

7.1 Studies with ideal MHD

Banerjee & Pudritz (2007), Seifried et al. (2011), Myers et al. (2013) and Rosen & Krumholz (2020) conducted simulations of the formation of massive stars under the assumption of ideal MHD. Although Banerjee & Pudritz (2007) observed magnetically driven outflows, their disk is always sub-Keplerian. Myers et al. (2013) started from a cloud of 300 M_⊙, an initial velocity profile with supersonic turbulence and no rotation, and $\overline{\mu }=2$ with B_z ∝ R^−1/2. They included gray radiative transfer and radiative protostellar feedback, and used a 3D AMR grid with a maximum resolution of 10 au, as well as an additional isothermal high resolution run (Δx = 1.25 au). They only find a disk in their high resolution run; it has a radius of 40 au when t ~ 0.2t_ff, for which M_⋆ ~ 3.5M_⊙. It is not trivial to compare their results to ours because of the difference in initial velocity fields, and the use of an isothermal equation of state. Using the scalability of our results with the mass of the cloud, we find the radius of the disk to be 128 au for t = 0.2t_ff in the simulation for $\overline{\mu }=2$. However, if we take M_⋆/M_C and compute the corresponding fraction of the free-fall timescale, we get t/t_ff ~ 0.104, for which our disk has a radius of ~30 au, which is in line with the fact that supersonic turbulence delays the formation of an accretion disk as enough angular momentum has to assemble close to the forming star. Rosen & Krumholz (2020) also considered a cloud core with supersonic turbulence and radiation transport, however their coarse grid (minimum cell size of 20 au) only allows them to partially resolve a disk-like structure and therefore a direct comparison to our results is not possible.

Seifried et al. (2011) considered a setup similar to ours in terms of the mass and radius of the cloud, as well as its initial density and rotation profiles, but used a cooling function instead of radiation transport, and only considered the ideal MHD approximation. The code FLASH with an AMR grid of maximum resolution of 4.7 au was used. The authors explored several values of µt ranging from 2.6 to 26, but with a uniform plasma ß, and obtained a Keplerian-like accretion disk for $\overline{\mu }=26$. They observe a drop in Keplerianity for the inner parts of the disk in the case of $\overline{\mu }=10.4$, and credit it to magnetic braking, in agreement with our results; however, they do not observe a disk for $\overline{\mu }<5.2$ because the high density of the thin layer of the disk decouples the gas from the magnetic field im resistive MHD runs.

7.2 Studies including Ohmic resistivity

Matsushita et al. (2017) and Machida & Hosokawa (2020) use the same Ohmic resistivity model as we do (Machida et al. 2007), but both consider a barotropic equation of state instead of solving for radiation transport. In Matsushita et al. (2017), a disk was formed, however, their analyses are strongly focused on the magnetic outflows. Machida & Hosokawa (2020) consider cloud cores of radius 0.2 pc with masses ranging from 11 to 545 solar masses, an enhanced Bonnor-Ebert density profile, and slow initial solid-body rotation. If we expand our cloud to 0.2 pc, keeping the same density profile as in the fiducial case, the mass of the cloud becomes 283 M_⊙, which means that our setup is similar to their simulation series E. The authors also make a parameter scan with the normalized mass-to-flux ratio with values ranging from 2 to 20, and use a three-dimensional AMR grid with a maximum resolution of 0.62 au. After 10 kyr of evolution, they find disk radii of a few hundred astronomical units, which develop spiral arms and fragments for some configurations.

7.3 Studies including ambipolar diffusion

Mignon-Risse et al. (2021) and Commerçon et al. (2022) studied the formation of massive stars from 100 M_⊙ cloud cores of radius 0.2 pc, with a centrally condensed initial density profile; which roughly correspond to our configuration M_C = 50 M_⊙. The cloud cores are initially in solid body rotation and are threaded by a uniform magnetic field determined from the normalized mass-to-flux ratio $\overline{\mu }=5$ and 2. They treated radiation transport with the flux-limited diffusion approach and gray stellar irradiation; for magnetic diffusivity, they considered ambipolar diffusion using the model from Marchand et al. (2016) but no Ohmic resistivity, and ran their simulations on a three-dimensional AMR grid with the code RAMSES down to a maximum resolution of 5 au. In a nutshell, those studies constitute the ones which include the most similar physical ingredients but to our setup but with ambipolar diffusion instead of Ohmic dissipation as the resistive effect and utilizing a different grid method.

Both studies report on thin accretion disks which are vertically supported by thermal pressure, in agreement with our results; however, they do not observe a surrounding thick layer of the disk supported by magnetic pressure. Opposite to our findings, the authors find smaller disks with stronger magnetic fields, however, they do find a bigger disk for their ideal MHD run compared to the runs with ambipolar diffusion. In that case, the disk is not only bigger, but also inflated and therefore less dense, in agreement with our results. They attribute the finding of smaller disks to higher ionization in the outer disk, and therefore more magnetic braking. In turn, our results indicate that magnetic braking is highest in the inner disk, where the fluid and the magnetic field lines are dragged faster by rotation, although resistivity causes the dragging to be partial and delayed because of the dominance of magnetic diffusion. Moreover, in our results for the ideal MHD case we observe simultaneously the largest disk and the strongest magnetic braking, because it mostly affects the inner disk.

A comparison of the radius of the disk that we report here and the disk size reported in Commerçon et al. (2022) can only be made qualitatively and with several considerations in mind, apart from the difference in the nonideal MHD effect considered. First, the density profiles in both studies are different because of the presence of an initial density plateau in the inner ~4125 au in Commerçon et al. (2022) and Mignon-Risse et al. (2021), in contrast to our continuous power-law slope down to the smallest scales. This density plateau affects the gravitational collapse and therefore, the elapsed time and mass of the protostar cannot be directly compared. Second, the initial rotational to gravitational energy ratio of the cloud is not the same in the fiducial cases of both studies, but we can roughly compare our ζ = 0.04 case to the “fast” case studied in Commerçon et al. (2022), which has ζ = 0.05. This is of special importance because of what is shown in Fig. 9C4: the radius of the disk is strongly influenced by the initial rotation of the cloud, even more than the initial magnetic field (cf. Fig. 9B4) or resistivity (Fig. 9B2). Finally, the disk identification criterion is not the same in both studies: while Commerçon et al. (2022) use a set of criteria that involves the radial and vertical dynamics of the disk as well as a density threshold, we consider a purely dynamical criterion (radial equilibrium between gravity and centrifugal force) as it is more usual in observational studies. With this in mind, we note that in our simulation M_C = 50 M_⊙ [×4], the radius of the disk is around 470 au after 27 kyr, while their disk is around the same value after 51 kyr. However, the mass of the protostar at that time is lower in our case (4.6 M_⊙) than theirs (9.1 M_⊙); both the difference in time and mass can be reasonably attributed to the initial density profile. Additionally, while we find that the disk grows in size with time, most of the nonideal MHD simulations by Commerçon et al. (2022) find disks that stop growing at around 100 au.

Commerçon et al. (2022) and Mignon-Risse et al. (2021) do not report the existence of a magnetically braked region, for which we offer the following explanations, apart from the different criteria for disk identification we already discussed. First, ambipolar diffusion adds magnetic diffusivity, which may have the effect of delaying the formation of the magnetically braked region. In more detail: we report on the fact that higher Ohmic resistivities delay the onset of magnetic braking and with it, the formation of the magnetically braked region. Ambipolar diffusion also decouples the magnetic field from the flow (even though it may act in a different direction than Ohmic resistivity), which leads us to assume that a similar delaying effect on magnetic braking could be obtained upon the inclusion of ambipolar diffusion in our calculations. We note, however, that if this delay in the onset of magnetic braking is longer than the timescale of gravitational collapse of the cloud, the massive star and the disk may form without ever developing a magnetically braked region. Second, spatial resolution: their 3D AMR grid has a minimum cell size (maximum spatial resolution) of 5 au, while our grid has minimum cell sizes on the order of 10⁻¹ au close to the sink cell for the simulation series ×1, ×2 and ×4, and 10⁻² au for the simulation series ×8 and ×16. Third, their sink particle algorithm requires the definition of an accretion radius, set to four times the minimum cell size, that is 20 au, which is the size of the magnetically braked region for most of the time in our fiducial case. As a consequence, gas orbiting the massive protostar within the accretion radius or in the forming magnetically braked region cannot be resolved on their finest AMR level while it is resolved on our spherical grid.

Finally, we mention the study by Masson et al. (2016), who carried out simulations of low-mass star formation considering both ideal MHD and ambipolar diffusion. Although we cannot compare our results quantitatively with theirs, there are several qualitative similarities. The authors observe an accretion disk vertically supported by thermal pressure, with magnetic pressure dominating above the disk. Their figure 14 shows a structure with relatively high angular momentum enveloping the disk for both the diffusive and ideal cases, which is reminiscent of the thick layer of the disk we observe, however, due to the different definitions of the disk used in both studies, we cannot establish a direct correspondence. The same figure also shows that in the ideal MHD case, the inner parts of the disk are destroyed by magnetic braking at late times. Nevertheless, their corresponding run that considers ambipolar diffusion also exhibits for late times a region of low angular momentum (probably due to magnetic braking) in a conical region around the inner disk and that could deliver material to it, in a way that is reminiscent of the early stages of magnetic braking we observe in our simulations. Those results seem to support the idea that increasing the diffusivity leads to a delay but not complete suppression of magnetic braking in a magnetized disk embedded in a collapsing cloud.

8 Summary and conclusions

We have modeled the formation of a massive star with a series of 30 magnetohydrodynamical simulations including Ohmic resistivity, radiation transport from thermal emission of the dust and gas, and self gravity. The series of simulations covered a wide range of cloud masses, magnetic field strengths, density profiles, rotation profiles, and ratios of rotational to gravitational energy, which is in line with the currently estimated values from observations. We also performed a convergence study to test the robustness of our results.

After analyzing the fiducial case of our parameter study in depth, we found the following general features of the system:

• After the initial gravitational collapse, a Keplerian-like accretion disk was formed.

• The accretion disk is divided into two layers: a thin layer supported vertically by thermal pressure, and a surrounding thick layer supported by magnetic pressure. The thin layer of the disk appears only in simulations with sufficiently high resolution.

• At early times magnetic diffusion due to Ohmic resistivity is strong in the inner parts of the disk and it greatly reduces magnetic braking there during the magneto-centrifugal epoch (5 kyr ≲ t ≲ 15 kyr). As time progresses, and the magnetic field lines are continuously dragged by rotation, magnetic braking is observed in the innermost ~ 50 au of the disk for the fiducial case in our parameter space.

• Magnetic pressure can increase the size of the accretion disk.

When examining the full parameter space of initial conditions for the onset of gravitational collapse, we find that:

• Our results for the size of the disk and the mass gain of the protostar scale with the initial mass of the cloud, despite the nonscalability of self-gravity and the thermodynamics considered.

• The thickness of the thick layer of the disk is controlled by the initial magnetic field strength.

• For a cloud with an initial density profile ρ ∝ r^−1.5 and in solid-body rotation, the disk grows roughly linearly in size as R_disk ≈ [6380 M_⋆/M_C − 98] au. The stellar mass grows approximately as M_⋆ ∝ (t/t_ff)^1.5..1.9.

• The size of the disk is more strongly determined by the initial distribution of density and rotational energy in the cloud than by the strength of the magnetic field.

• Multiple initial configurations of the cloud can produce a given disk size and (proto)stellar mass. This means that observations of disk-jet systems constrain (as opposed to determine) the possible conditions for the onset of gravitational collapse, and more measurements (distribution and strength of magnetic fields, for example) are needed to break the degeneracies.

These results indicate that the initial conditions of massive star formation determine the features of the system composed of the disk, jet and protostar. In the second part of this study (Paper II), we perform a dynamical analysis of the magnetically driven outflows of the same dataset we present here.

Acknowledgements

We thank Richard Nies for his contributions to the analysis of part of the dataset at the early stages of the project. G.A.O. acknowledges financial support from the Deutscher Akademischer Austauschdienst (DAAD), under the program Research Grants – Doctoral Projects in Germany, and complementary financial support for the completion of the Doctoral degree by the University of Costa Rica, as part of their scholarship program for postgraduate studies in foreign institutions. R.K. acknowledges financial support via the Emmy Noether and Heisenberg Research Grants funded by the German Research Foundation (DFG) under grant no. KU 2849/3 and 2849/9.

Appendix A Numerical convergence

As part of the parameter scan performed for this study, we investigated the effects of the grid and other numerically relevant parameters on our findings. Here, we focus on checking the effect on the extent of the disk, that is its inner and outer radii, which we present in Fig. A.1.

Appendix A.1 Sink cell size

We tried sink cell radii of r_sc = 1, 3, and 10 au as our inner boundary. As shown in panels A1 and A2 of Fig. A.1, the results for the radius of the disk and the magnetic braking radius are very consistent for r_sc = 1 and 3 au, but larger inner and outer radii at late phases of the simulation for r_sc = 10 au. Of special interest is the fact that for the smaller sink sizes, both the size and the time of formation of the magnetically braked region region does not change, which supports the idea that the effect is physical and not due to the inner boundary conditions.

During the formation of the massive star starting from the 100 M_⊙ cloud, the stellar radius computed with the stellar evolutionary tracks increases up to a maximum of 150 solar radii, which correspond to approx. 0.7 astronomical units. The use of a smaller sink cell would not be accurate because it would require the solution of the full set of stellar structure equations, which we do not do. With this in mind, we conclude that a sink cell size of 1 au is close to the reasonable limit for the approximations considered in this study, and a sink cell size of 3 au is enough to replicate the results of the smaller sink.

Fig. A.1 Variation of the radius of the disk and the magnetic braking radius for different the inner boundary (sink cell size), resolution of the grid, the α-viscosity used to model gravitational torques and the Alfvén limiter.

Appendix A.2 Resolution

We conducted the study of the fiducial case in five different resolutions, each one doubling the amount of cells in each one of the directions (radial and polar). The simulation ×16 was only run until t = 13.45 kyr due to the high computational cost. The radius of the disk is very similar in all of the simulations. The magnetic braking radius appears somewhat later in the high resolution simulations and it shows more variability for the grid ×8. This could be due to the increased resolution of the effects of the outflows, which will be discussed in Paper II. However, during the comparison of the disk radii, it is important to consider that our definition of the size of the disk, based on a simple threshold to Keplerianity, is affected by the size of the cells, rounding to the nearest cell that satisfies the criterion.

Appendix A.3 Viscosity

The α-viscosity model to gravitational torques produced by spiral arms in the self-gravitating disk introduces an additional parameter. For this reason, we investigated what happens if the value of α is changed to 0.75 and 1.5 from its fiducial value of 1. Panels B1 and B2 in Fig. A.1 show relatively consistent results in all cases, even though with some local differences. The magnetic braking region appears when the star is of the same mass in all cases, and the radius of the disk averages approximately to the same value over time. Even though viscosity transports angular momentum, the qualitative behavior is different than the effect of resistivity and initial magnetic field, both of which significantly alter the formation time of the magnetically braked region and its size. This marks the difference between magnetic braking and angular momentum transport by viscosity.

Appendix A.4 Alfvén limiter

We employed an upper limit to the Alfvén speed $${\upsilon }_A=\frac{B}{\sqrt{4\pi \rho }},$$(A.1)

by introducing a locally variable density floor such that the Alfvén speed never exceeds the limit, which is defined a priori. This limiter was introduced to avoid having the Alfvén speed increase indefinitely. The time step remains within limits that allow us to cover the desired timescales in a reasonable computing time. This approach has the disadvantage that artificial mass is created in the process, which in theory could affect the dynamics of the gas. A more detailed discussion of the effects of the Alfvén limiter will be offered in Paper II on jet physics, however, we mention here that with the values chosen for the simulation series ×4 and ×8, 2000 and 1000 km s⁻¹ respectively, artificial mass below 0.002 M_⊙ is generated at the end of the simulation. This value is very small compared to the mass of the cloud and should not change the dynamics of the disk and the outflows. In turn, for simulation series ×2, an Alfvén limiter of 2000 km s⁻¹ produces an additional 0.2 M_⊙ of artificial mass, with lower values yielding over 1 M_⊙, i.e., more than one per cent of the initial mass of the cloud. With these values of the limiter and the artificial mass in mind, we examine the extent of the disk in panels B3 and B4 of Fig. A.1 as a probe for effects in the dynamics of the disk and find that only in the case of an exceptionally low Alfvén limiter of 200 km s⁻¹ there are appreciable changes to the inner and outer radii of the disk for late times.

References

All Tables

All Figures

	Fig. 1 Initial conditions, overview of the evolution of the system, and summary of the parameter space presented in this study.
In the text

	Fig. 2 Keplerianity on the midplane, computed with the fiducial simulation on grid ×8.
In the text

Fig. 3 Dynamics of the thick layer of the disk. The panels of the figure are grouped in columns such that they show: the density structure of the thick layer (a and e), the offset from Keplerianity at three different heights in the thick layer (b and f), the contributions to the specific energy (c and g), and the vertical balance of forces (d and h). Data from the fiducial simulation on grid ×8 was used for this figure.

In the text

	Fig. 4 Contributions to the specific energy in the thin layer of the disk, calculated with the fiducial simulation on grid ×8.
In the text

	Fig. 5 Magnetic Reynolds number, computed with the fiducial simulation on grid ×8. For t = 10 kyr, the radius of the disk is around 80 au, while for t = 22.75 kyr, the exterior radius is ≈700 au and the interior radius is ≈20 au.
In the text

	Fig. 6 Distribution of the magnetic field strength as a function of density. This histogram was computed with the fiducial simulation run on grid ×16, for time t = 7 kyr. The color scale indicates the polar angle: the dark purple colors correspond to the midplane (with the disk corresponding to ρ ≳ 10−15 g cm−3) and the light blue colors correspond to the rotation axis (where the jet is located).
In the text

	Fig. 7 Density (panel a) and temperature (panel b) profiles on the mid-plane, for the fiducial simulation on grid ×8.
In the text

	Fig. 8 Density (panel a) and temperature (panel b) for two different altitudes in the thick layers of the disk (blue lines), and the corresponding profile for the thin layer (black). The data corresponds to the fiducial simulation on grid ×8.
In the text

Fig. 9 Radius of the disk (Rdisk) and its magnetic braking radius (Rmb) for different initial values of the mass of the cloud, magnetic field, density profile, rotational profile, rotational energy and resistivity models. The transparent lines in the panels that show Rmb indicate the full data, while the solid lines show the moving average. The colored dots in panel A4 indicate t = 30 kyr in each simulation. For all the panels in rows B and C, MC = 100 M⊙ was used.

In the text

Fig. 10 Mass of the protostar, corresponding to the mass in the sink cell, for several values of the (a–b) mass of the cloud, (c) resistivity model, (d) mass-to-flux ratio, (e) density profile, and (f) rotation profile and rotational energy. Panel b shows the same results as panel a, but normalized in such a way that scalability can be readily seen. For panels (c)–(f), MC = 100 M⊙ and tff = 53.73 kyr.

In the text

	Fig. 11 Density structure (Col. 1) and specific energies (Col. 2) of the disk for selected values of the normalized values of the mass-to-flux ratio (rows B and C), and rotation profile (Col. D). All the snapshots were taken at t = 15 kyr of evolution. All cases consider a rotational-to-gravitational energy of 4%.
In the text

	Fig. A.1 Variation of the radius of the disk and the magnetic braking radius for different the inner boundary (sink cell size), resolution of the grid, the α-viscosity used to model gravitational torques and the Alfvén limiter.
In the text