Protostellar disk formation by a non-rotating, non-axisymmetric collapsing cloud: model and comparison with observations

Planet-forming disks are fundamental objects thought to be inherited from large scale rotation, through the conservation of angular momentum during the collapse of a prestellar dense core. We investigate the possibility for a protostellar disk to be formed from a motionless dense core which contains non-axisymmetric density fluctuations. The rotation is thus generated locally by the asymmetry of the collapse. We study the evolution of the angular momentum in a non-axisymmetric collapse of a dense core from an analytical point of view. To test the theory, we perform three-dimensional simulations of a collapsing prestellar dense core using adaptative mesh refinement. We start from a non-axisymmetrical situation, considering a dense core with random density perturbations that follow a turbulence spectrum. We analyse the emerging disk comparing the angular momentum it contains with the one expected from our analytic development. We study the velocity gradients at different scales in the simulation as it is done with observations. We show that the angular momentum in the frame of a stellar object which is not located at the center of mass of the core is not conserved, due to inertial forces. Our simulations of such non-axisymmetrical collapse quickly produce accretion disks at the small scales in the core. The analysis of the kinematics at different scales in the simulated core reveals projected velocity gradients of amplitudes similar to the ones observed in protostellar cores, and which directions vary, sometimes even reversing when small and large scales are compared. These complex kinematics patterns appear in recent observations, and could be a discriminating feature with models where rotation is inherited from large scales. Our results from simulations without initial rotation are more consistent with these recent observations than when solid-body rotation is initially [abridged]


Introduction
Protoplanetary disks are rotationally supported structures that formed around young stars (Li et al. 2014;Dutrey et al. 2014;Testi et al. 2014). It is currently believed that the rotation of these disks is inherited from large scales -a few thousands of astronomical units, the scale of the parent prestellar dense core.
During the gravitational collapse of the core, if the angular momentum is conserved, the infalling material naturally forms a rotation dominated structure at the small scale of a hundred astronomical units. The rotation of the disk is thus inherited from the large scale angular momentum, and as a consequence, the velocity gradients at large and small scales are correlated. This scenario is extensively studied in the literature and in particular a majority of collapse calculations starts with a prescribed rotation profile (see for example Bate (1998); Matsumoto & Hanawa (2003); Machida et al. (2005); Hennebelle & Fromang (2008)). While reasonable, this scenario leads to the question at which scale is the angular momentum being inherited from and how exactly this happens. Another frequently configuration consists in a cloud with a turbulent velocity field imprinted initially (Bate et al. 2003;Goodwin et al. 2004a,b;Dib et al. 2010;Hennebelle et al. 2016;Gray et al. 2018;Kuznetsova et al. 2019). In this context as well, it has been found that disks quickly form. The usual interpretation is that the angular momentum is initially present because of the turbulence.
Observationally, the kinematics of the dense gas in both prestellar cores and protostellar envelopes has been studied thanks to the analysis of molecular line emission and shown to harbor velocity gradients at scales 0.01 to 0.1 pc, interpeted as rotation of cores (see the early works by Goodman et al. (1993); Ohashi et al. (1997); Chen et al. (2007)). The values of spe-Article number, page 1 of 15 arXiv:2002.02144v1 [astro-ph.GA] 6 Feb 2020 A&A proofs: manuscript no. ms cific angular momentum measured inside protostellar envelopes at scales a few thousand astronomical units are on average one order of magnitude lower than the ones observed at larger scales in starless structures (a few 10 −4 km.s −1 .pc, see for example Belloche (2013), Yen et al. (2015a) and the very recent work by Gaudel et al. (2020)). More puzzling, however, were the observations showing that some protostellar cores show an apparent disorganization or even reversal of their velocity pattern, sometimes interpreted as the contribution of infall motions to the projected velocity field (Tobin et al. 2011;Harsono et al. 2014) or as counter rotation (Tobin et al. 2018;Takakuwa et al. 2018). Recent high dynamic range observations of a sample of 12 lowmass Class 0 protostars (in the calypso sample) by Gaudel et al. (2020) exhibit systematic dispersion of velocity gradients between disk's and envelope's scales, questioning the presence of large scale rotation. Moreover, observations of the specific angular momentum contained in T-Tauri disks suggest values larger, by about one order of magnitude, than the specific angular momentum observed in the low-mass protostellar cores at scales a few thousands astronomical units (Simon et al. 2000;Kurtovic et al. 2018;Pérez et al. 2018). These observations are hence difficult to reconcile with a simple picture of rotating-infalling protostellar envelope which conservation of angular momentum naturally produces a rotationally-supported disk in its center, and new models should be developed to reproduce these observations as well.
To tackle these issues we investigate a scenario that also leads to the formation of a protostellar disk. Similar ideas than the ones exposed in this paper has already been developed in the context of spiral galaxies formation by Hoyle (1949); Sciama (1955); Peebles (1969). In the context of protoplanetary disk formation this paper is meant to be exploratory, so we used minimal physical ingredients. Our simulations are thus purely hydrodynamics, except in part 4.6. We start from an extreme scenario, considering a perfectly motionless dense core with nonaxisymmetric density perturbations. The gravitational collapse of this core is thus non-axisymmetric. We show analytically and numerically that this non-axisymmetry leads to the possibility to generate rotation locally. As we start from an extreme motionless scenario, without considering velocity fluctuations, our model is not fully physical. Despite this fact, we then analyse the velocity gradients in our simulations and we reproduce the observational results from Gaudel et al. (2020) about the dispersion of velocity gradients. We show that the specific angular momentum step coincides with Belloche (2013) and Gaudel et al. (2020) results. Our model thus exhibit good agreement with observational constraints on kinematics.
The plan of the paper is as follows: in the second part we present the theory that motivates our study, in the third part we present the numerical methods we used to investigate our problem, in the fourth and fifth parts we present the results obtained and their discussion, and the sixth part is the conclusion.

The axisymmetrical case
In the introduction, we evoked the conservation of angular momentum during the collapse of a dense core, leading to the formation of a disk. However, the angular momentum is correctly defined only in a given frame and with respect to a given point. Let's consider the angular momentum calculated in the simulation box frame R, with respect to the center O of this box 1 . It 1 In fact, in relation to any fixed point of the simulation.
is computed as follows, with the summation referring to the different cells i of the simulation, m i and M i being respectively the mass and position of each cell i.
This momentum is conserved in virtue of the fundamental law of evolution of the angular momentum in a Galilean frame: As no external force is applied on the system, the angular momentum σ O | R is conserved. In a simple axisymmetrical case, this momentum coincides with what we will call the momentum of the disk. Indeed, during the collapse, a disk forms in the center of the box, thus σ O | R represents the angular momentum computed in the frame of the disk, in relation to the center of the disk.

Non-axisymmetric configuration
In a non-axisymmetrical case, σ O | R is no longer a relevant quantity to study the disk formed in the simulation. In fact, to measure the angular momentum in protostellar disks, the reference point with respect to which the angular momentum is computed is the center of the disk, and the velocities considered are those in relation to the center of the disk, deducted from those projected on the line of sight (Belloche 2013). In the axisymmetrical case, the center of the formed disk remains motionless at the center of the simulation box. In the non-axisymmetrical case, the disk is not formed at the center of the simulation box, and have a proper motion. We thus have to consider the angular momentum computed in the frame R of the disk, in relation to the center C of the disk: We show in appendix A.1 that for an initial condition where all cells are at rest in R, σ C | R can be expressed as: where M = i m i is the total mass of the system, and G is the center of mass 2 . Furthermore, the time derivative of Eq. (4) gives: This equation can be interpreted as the variation of the angular momentum σ C | R due to the torques of inertial force that apply on each cell of the simulation in the non Galilean frame R (see appendix A.2 for detailed development). As the center of mass G and the center of the formed disk C do not coincide, the angular momentum σ C | R is not expected to be conserved. The point C being the accretion center of the system, the matter collapses toward it. As the angular momentum (in R with respect to C) of this matter does not vanish, a rotationally supported structure forms around C.
In the axisymmetrical case, the three points C, G and O coincide, as well as the two frames R and R . Thus σ C | R = σ O | R and the angular momentum is therefore conserved through the temporal evolution of the structure. In the non-axisymmetrical case, Eq. (4) shows that the angular momentum σ C | R is equal to the angular momentum of the point C to which the whole simulation mass has been allocated, computed with respect to G, in the frame R of the simulation box. Eq. (5) shows that this angular momentum is not conserved in the general case. We consider here the extreme case where every cells are initially at rest in the simulation frame R. σ C | R is thus initially null, C being the accretion center of the system. As matter is falling toward the center, the transversal velocity of this matter has to increase for the angular momentum to grow. As a result, a rotational structure naturally forms around the accretion center, without violating the conservation of the angular momentum in the simulation frame R. These ideas are similar than the ones of Peebles (1969) who showed that the angular momentum of spiral galaxies can be understood as a consequence of tidal torques acting during the gravitational collapse.
In our context of protostellar disk formation and before analysing in details our simulations in part 4.2, we can provide an order of magnitude for the expected momentum σ C | R . The simulations show that the distance between the points C and G is about a thousand astronomical units, which corresponds roughly to a tenth of the dense core radius -875 AU. Let's take for dGC dt the order of magnitude of the sound speed, 0.2 km.s −1 . The total mass inside the simulation being 2.5 M , we obtain: σ C | R ∼ 2.5 M . 10 3 AU . 0.2 km.s −1 ∼ 10 47 kg.m 2 .s −1 As will be seen later, this order of magnitude is consistent with the values given by our analyses (see Fig. 3) and stresses good agreement with Belloche (2013) and Gaudel et al. (2020).

Code and numerical parameters
To study whether the asymmetry of a gravitational collapse can be sufficient to form a prestellar disk, we carry out a set of hydrodynamics simulations with Ramses (Teyssier 2002). This numerical Eulerian code uses Adaptative Mesh Refinement (AMR) technique to enhance resolution locally, where it is needed, on a Cartesian mesh. Our refinement criterion is based on Jeans length such as each local Jeans length is described by at least 40 cells. We use ten levels of AMR, from 8 to 18, leading to a spatial resolution which goes from 270 AU to 0.26 AU.

Initial conditions
We consider a 3D cubic box with sides of 0.33 pc (about 70 000 AU). At the middle of the box we place a 2.5 M sphere of gas which diameter is the quarter of the box length and which radial density profile is flat. This sphere acts as a model of prestellar dense core. The rest of the box is filled with an envelop of gas whose density is constant and equal to a thousandth of the mean density of the gas sphere. Initially all velocities are set to zero so that the angular momentum in the box frame with respect to any motionless point is initially null. The alpha parameter -thermal over gravitational energy ratio -of the gas sphere is 0.35. We use a barotropic equation of state for the gas: T and ρ being the temperature and density of the gas, T 0 = 10 K, the critical density ρ c = 10 −13 g.cm −3 and the adiabatic index γ = 1.4. As we want to study the effect of asymmetrical gravitational collapse, we break the symmetry of the cloud introducing random density perturbations in the dense core. To mimic roughly the physics of the interstellar gas, we based the probability distribution of the density perturbations on the one of the turbulence (see for example chapter 3 of Hennebelle & Falgarone (2012) for a review on turbulence in interstellar clouds). In the Fourier space the perturbations spectrum matches with the power spectrum of the velocity for Kolmogorov turbulence 3 , and the phases are randomly chosen. We can thus write the density of each cell i of the prestellar core as : where d 0 represents the mean density of the prestellar core, δρ i is the value of the perturbation at the considered point, and A ∈ [0, 1] is an internal parameter allowing us to control the amplitude of the perturbation. To ensure that the density stay positive everywhere, it is necessary that δρ i ∈ ]−1, +∞[ . To satisfy this condition, we modified all the values of δρ i < −1 to bring them to −0.99. The rotation axis z d of the disk is not overlapping with any simulation box axis. We see that a disk forms and grows in spite of the fact that initially the angular momentum is null.
To assure that the mean density of the prestellar core stay constant when vaying A, we renormalised 4 the value of d i defined in Eq. (7) in each cell : where the operator < . > i represents the mean value over all the cells i of the prestellar dense core. As the mean value of δρ i is not null, this operation warrants that we can change the amplitude of the perturbations without modifying the mean value of the density.

Choice of a time reference
Since the free-fall time depends on the density, it can slightly vary from a level of perturbation to another. As the grid is initially coarse, the first time-steps of the simulation are much larger than the later ones after collapse occurred, when the level of refinement is higher. These two effects lead to a bad description in time at the beginning of the simulation. We are thus compelled to choose a time reference from which the ages are computed. We based our time reference on the maximum density inside the simulation box. We take for time reference the moment when the maximum density reaches 10 −13 g.cm −3 . It corresponds to the limit density beyond which the compression of the dense gas changes from isothermal to adiabatic (Larson 1969).

Formation of a disk
The main result of this study is, as we expected from our theoretical development, the formation of a disk within our simulations. We show a sequence of images of the disk growing with time in the simulation with 50% of density perturbations in Fig. 2. This disk is very similar to the ones found in many studies (Bate et al. 2003;Matsumoto & Hanawa 2003;Goodwin et al. 2004a;Hennebelle & Fromang 2008;Gray et al. 2018). As can be seen it presents prominent spiral arms that transport angular momentum. To verify that this structure is a disk, we made sure to check that the formed structure is rotationally supported. In order to keep our model as simple as possible, we did not use sink particles (Krumholz et al. 2004;Bleuler & Teyssier 2014) in our simulations presented in this part and analysed in parts 4.2 and 4.3. As a consequence, the dense gas accumulates in the disk, making it self-gravitating. Thus the azimuthal velocity profile shows substantial deviations from the keplerian profile, but the structure is still rotationally supported, as the azimuthal velocity is much larger than the radial velocity inside the disk, typically by a factor 50. We have however also showed that disks form in simulations with sink particles (see part 4.4) and we verified that the formed structures match very well the keplerian profile. We run a set of simulations which differ only by their perturbation level, from 10% to 60%. In all of these simulations we observe the formation of a protostellar disk. The comparison with theory, studied in details for the simulation with 50% of density perturbations in part 4.2 gives clues to assess the trust level of our simulations (for the other levels of perturbation, the results are visible in appendix B). For the simulations over 20% of perturbations, the comparison showed that the formation of the disk can be trusted. For the simulations with low perturbation level, typically less than 20%, the agreement with theory is less good and the numerical errors are higher.

Comparison with theory
We showed in part 2.2 that a rotational structure should emerge from an asymmetric gravitational collapse, even without any initial motion. To compare this theoretical prediction with our simulations, we exploit Eqs. (3) and (4). The first expression simply expresses the numerical way to compute the total angular momentum in the simulation (in the frame R , with respect to the point C, what we will not mention anymore). Let's name this quantity σ num . The second one is the analytic expression of the same quantity. It is equivalent to the first one if a set of hypothesis is verified (see appendix A.1), which is the case in our sim-ulations. Let's name this quantity σ an . The equality of the two quantities σ num and σ an ensures that we are observing the physical phenomenon we described. The differences could be due to numerical errors. In the first place, we verify numerically the equality: The comparison between σ num and σ an is represented on the top panel of Fig. 3 for the simulation with 50% of perturbations. The slight relative difference -globally less than 2% -shows that the equality of these two quantities is numerically consistent. Furthermore, the relative difference does not keep the same sign during the temporal evolution. It shows that this difference is not a systematic error on one of the two quantities.
Let's name |∆σ| = |σ num − σ an |. To be entirely sure that the disk formed in our simulation is not a numerical artifact, we verified that the angular momentum it contains is larger than |∆σ|.
In the most pessimistic scenario where |∆σ| would be entirely concentrated in the disk, this ensures that |∆σ| is not sufficient to explain the presence of the disk. The angular momentum contained in the disk is computed as σ num , but in the restricted area of the simulation corresponding to the disk 5 . The comparison between |∆σ| and the angular momentum in the disk is visible on the bottom panel of Fig. 3. As |∆σ| is smaller than the angular momentum in the disk, and as ∆σ is switching sign over the temporal evolution of the system, it confirms that the formed disk is the result of the physics described part 2.2. For the other levels of perturbations, the results are presented in appendix B.

Analysis of velocity gradients
In this section we analyse our simulations from an observational point of view to highlight whether or not our model succeed to describe some features of real observations. As the rotation emerges from the asymmetrical gravitational collapse in our model, it is not straightforwardly inherited from larger scales. This lack of connection between disk and envelope scales leads to important features. We compute three quantities at different scales that are accessible to observations: the direction and amplitude of velocity gradients and the specific angular momentum.
To compute direction and amplitude of the velocity gradients at a scale R we follow the method of Goodman et al. (1993). We consider a cube of side 2R around the disk, aligned with the three main axis of the simulation. We consider these three main axis as our lines of sight to compute maps of projected velocities weighted by density, on a depth of 2R. Then we fit these maps with a solid-body rotation profile: with v 0 the systemic velocity of our object, ∆x and ∆y the vertical and horizontal dimensions of our projected velocities map. The magnitude of the velocity gradient is thus Ω = √ a 2 + b 2 and its direction is given by θ = arctan a b . The specific angular momentum at a scale R is thus given by j = R 2 Ω.  (3)) in blue, and the relative difference between σ num and σ an (see Eq. (4)) in purple for each output. Bottom: the momentum in the disk (in violet) compared to the absolute value of the difference between σ num and σ an (in beige). We see that ∆σ σ is less than 2.5% and that this difference is smaller than the momentum contained in the disk. The disk is thus resulting in the physics we described part 2.2.
In an axisymmetrical model with initial rotation, velocity gradients at small and large scales are perfectly aligned. In our model, due to the lack of initial rotation, it is interesting to see how velocity gradients at different scales are organised. The top panel of Fig. 4 shows the angular variation of velocity gradients according to the probed scale, relatively to the velocity gradient at the disk scale. Velocity gradients in the disk and in the envelope are misaligned. For the y projection, the velocity gradient in the envelope is even reversed in comparison to the small scale gradient, whereas for the x projection the velocity gradients make a complete turn from the disk to the envelope scales.
The amplitude of the velocity gradient at different scales is showed on the mid panel of Fig. 4. In the three projections the amplitude profile from a hundred astronomical units roughly follow a power law j ∝ R −1.8 . It shows that these gradients are detectable in real observations, as the choice of appropriate molec- ular lines allows to detect gradient amplitude about 1 km.s −1 .pc −1 in a solar type star-forming core at a distance of 200 pc. Once we have computed the amplitude of these gradients, we can determine the specific angular momentum. The evolution of this quantity at different scales is presented on the bottom panel of Fig. 4. This figure shows that for the three main lines of sight of the simulation, the specific angular momentum do not vary so much through the different scales, excepted some peaks correlated to abrupt changes in angular direction of velocity gradients. Furthermore, this quantity is roughly constant in the envelope, at the scale of the hundred and thousands astronomical units, with a mean value of about 3 .10 −4 km.s −1 .pc. In the discussion section, we compare the projected kinematics properties of our modeled core to observations in protostellar envelopes.
To determine the origin of these velocity gradients, we compute the maps of projected velocities along the line of sight taking only the radial part of the velocity with respect to the disk center on the one hand, and the orthoradial part of the velocity on the other hand. We compute the velocity gradients directions over the different scales for the radial and orthoradial part of the velocity as we did for the total velocity gradients. We compared the angular deviation between the total velocity gradient and the radial velocity gradient -∆θ rad -, and between the total velocity gradient and the orthoradial velocity gradient -∆θ orthorad . The results for different levels of perturbation taken at similar times are visible Fig. 5. For the simulation with 10% of perturbations, the total and radial velocity gradient directions are very close over the probed scales 6 . As the level of perturbation increases, the radial and total gradients begins to misalign -|∆θ rad | increase -for scales larger than 400 AU, whereas the orthoradial and total gradients directions gets closer -|∆θ orthorad | decrease. These results depend on the line of sight choosen to compute the projected velocity, and to a lesser extent on the evolutionary stage of the simulations. What is a common feature is that for the small perturbation levels, the radial and total velocity gradients directions are very close over all the probed scales.

Size of formed disks
We did not introduce until now any sink particle in our simulations. As a consequence, the gas cannot collapses at a smaller scale than our maximal resolution. This compels the gas to accumulate in the center of the disk, causing the disk to become autogravitating and leading to the apparition of the spiral arms visible in Fig. 2 to transport angular momentum inside the disk. As this accumulation of dense gas in the disk is not physicalthe gas should continue to collapse to the stellar scales -it is not correct to compute disk radius in these simulations. To handle this issue, we ran simulations with sink particles introduced when density reaches the threshold of 10 14 cm −3 . The disk grows but is less massive than previously, and match very well the keplerian profile. At some point the disk fragments. Figure 6 shows the disk at an advanced stage, shortly before fragmentation. The disk reaches radius of about 125 AU, which is large in comparison to observed disk radii (Maury et al. 2010;Tobin et al. 2015;Segura-Cox et al. 2016;Maury et al. 2019). To compute the disk radius in our simulations, the first step is to isolate what belongs to the disk. This selection is based on density and velocity criterion: a cell has to be dense enough and the orthoradial component of its velocity has to be at least two times larger than the Fig. 5: Angle between radial and total velocity gradients (∆θ rad , solid lines) and between orthoradial and total velocity gradients (∆θ orthorad , dashed lines), for different levels of perturbation ε. For the simulation with 10% and 20% of perturbations the total gradients directions follow the ones of the radial velocity gradients over all scales. For the simulations with more perturbations, ∆θ rad is higher from 400 AU and ∆θ orthorad gets lower. radial one. Once this selection done, we look at the maximum distance projected in the equatorial plane for a large number of angular sectors. The mean of these distances over the different angular sectors is taken as the disk radius. Figure 7 shows the evolution of disk radius over time for several simulations. The three green curves show the radius of the disk in the simulations with 10, 20 and 50% of perturbation level. It appears that in these three simulations the disks grow to reach radii around 200 AU in 20 to 30 kyr, before fragmenting. The disk size do not depends much on the perturbation level in the range of 10 to 50% of perturbations.

Effects of initial rotation
To see the influence of initial rotation on our results, we ran four simulations in which a solid-body rotation velocity profile is imprinted in the initial conditions. We ran simulations with ε = 20% and 50%, and for each of these perturbation levels we choose two rotation levels 7 β = 0.25% and 1%. The evolution of the disk radius in these simulations is visible in Fig. 7, represented by the four curves from red to yellow. For these four simulations the evolution of the radius is similar. These disk grows rapidly before fragmenting between 2 to 7 kyr at a radius around 100 AU. This evolution is different from the ones of the disks in the simulation without initial rotation, in particular because the disk forms earlier and grows faster. While when there is no rotation, it takes about 15 kyr to get a disk bigger than 100 AU, this takes only 3 to 5 kyr when rotation is included.
We conducted in these simulations the same analysis of velocity gradients than in part 4.3. The results for the simulation with ε = 50% and β = 1% are shown in Fig. 8 and 9. In the top panel of Fig. 8, we see that the angular deviations of the velocity A&A proofs: manuscript no. ms Disk radius (AU) The three green curves are purely hydrodynamics simulations with different levels of perturbations ε. The four curves from red to yellow are also purely hydrodynamics, but include an initial solid-body rotation velocity profile, with different levels of perturbations and rotation β. The curves stop when the disk fragments. All the plotted curves have been smoothed by a sliding median.
gradients are less important than in the case without initial rotation. Only the face-on projection exhibit large deviations, but the amplitude Ω in this projection is much smaller than in the two edge-on projections, by a factor around 10. The resulting specific angular momentum is visible on the bottom panel of Fig. 8. For the two edge-on projections, it exhibit values that are larger than in the case without initial rotation, and that are larger than the ones deducted from observations. We also see that j ∝ R α with α 0.5 while observations revealed that in the inner few thousands AU, α 0 (see (Belloche 2013) for a review). The most recent observations by Gaudel et al. (2020) reveal α 0.3 ± 0.3 under 1600 AU.
Here also we computed the decomposition of the velocity in its radial and orthoradial components and conducted the same analysis of velocity gradients with these two components. The results are visible in Fig. 9 for the full velocity gradients, the radial velocity gradients and and orthoradial velocity gradients, in the simulation with ε = 50% and β = 1%, for the edge-on projection 1. The top panel shows the angle using the same arbitrary reference for the three gradients, whereas the middle and bottom panels show respectively the amplitude and the specific angular momentum deducted from these gradients. The joint analysis of the top and middle panel shows that for R < 10 2 AU and R > 3 .10 3 AU the orthoradial velocity gradient is dominant, as it is much larger in amplitude than the radial one and much closer to the full velocity gradient in angular position than the radial one is. This means that at large scale the full velocity gradient is tracing the rotation of the envelope. For 10 2 AU < R < 3 .10 3 AU, the contribution of the radial and orhoradial component are of the same order of magnitude in amplitude, which leads to the full velocity gradients to be misaligned with both the radial and orthoradial gradients at those scales. In the face-on projection the effect of the rotation is nearly invisible and the results are similar to the ones presented part 4.3.
The results are comparable for the four simulations introduced in this part and therefore the other cases are not shown for conciseness.

Effects of the magnetic field
For the sake of completeness, we finally add magnetic field 8 in our simulations to see its effects on the formed disk. The magnetic field is treated under the ideal MHD approximation. In these simulations a disk still forms even in the absence of initial large scale rotation. Figure 10 shows the appearance of the disk  at the same time than Fig. 6. Clearly the two disks are qualitatively very different. In the purely hydrodynamics case, the disk is big and massive, with sharp edges, whereas in the MHD case, the disk is smaller and less massive (smaller column density).
In the MHD case it is hard to define properly a disk radius. In fact, when looking at Fig. 10, the disk is buried in a very filamentary structure that has the same order of magnitude in column density as the outer part of the disk. At some places these filaments verify the velocity criteria set out in part 4.4 to define the belonging to the disk. Here we are confronted to a definition problem. In hydrodynamics simulations, the disk has sharp edges and it is clear to define by eye what belongs to the disk or not. In MHD simulations it is really hard to do so and therefore we do not attempt here to present a quantitative comparison.
In spite of this difficulty, the comparison between Figs. 6 and 10 reveals that in the MHD case, the disk is much smaller. For the MHD simulation without initial rotation, the disk stops growing after 20 kyr and does not fragment like in the hydrodynamics case. For the simulation with initial rotation, the disk grows more rapidly but also stabilises at the same size between 40 to 70 AU depending on the definition of what belongs to the disk. We believe that this is the signature of the magnetic braking occurring in magnetised disk.

Discussion
From a very simple model, we showed that the angular momentum computed in the frame of the disk in relation to the center of the disk is not a conserved quantity. This is due to the non Galilean nature of the frame, provoked by the nonaxisymmetrical gravitational collapse. This collapse leads to the formation of a protostellar disk. In our system, the angular momentum that forms the disk is not a mere conservation of preexisting large scale angular momentum. The rotation can be generated locally by the asymmetry of the collapse. This result is in agreement with the early work by Peebles (1969). This mechanism responsible for the formation of spiral galaxies is thus also able to account for the formation of protostellar disks. Thereby, the conservation of a pre-existing angular momentum at large scales might not be the only mechanism responsible for protostellar disk formation, and the assumption that small scale rotation is being inherited from the large scales may not be correct, at least for some systems.
The disks formed in the simulations with sink particles reach radius larger than a hundred astronomical units until they fragment. As these simulations start with no angular momentum, but only require moderate density perturbations -which are expected to be present in real cores -, this study gives a sort of "minimal hydrodynamic radius" that a disk should reach in a hydrodynamics collapse. As most of the observed disks are much smaller (Maury et al. 2010;Tobin et al. 2015;Segura-Cox et al. 2016;Maury et al. 2019), it implies that angular momentum extraction processes are at work in real disks. In the traditional picture of axisymmetric collapse, the size of the disk directly depends on the initial rotation. The presence of small disks in observations could therefore be interpreted as due to low rotation levels. In the present study, we show that in the nonaxisymmetric case, even modest levels of density perturbations lead to the formation of large hundred of AU size disks. These latter appear to be natural outcome of the gravitational dynamics and local conservation of angular momentum. These are generic features. Therefore the most natural way to reconcile these observations with the generic nature of big hydrodynamical disks is to invoke the presence of magnetic braking that we know is operating in cores (Li et al. 2014;Hennebelle & Inutsuka 2019). Our MHD simulation starting from the same initial conditions but just adding a magnetic field shows that the outcoming disk is much smaller and do not fragments, as expected from these previous studies.
The magnitudes of the velocity gradients from our modeled core do not depend much on the chosen projection, although a dip is observed in the z projection at scales of 600 AU (see mid panel of Fig. 4) which results from a reversal in angular direction of velocity gradients. Their amplitude is roughly consistent with velocity gradients observed in protostellar cores using dense gas tracers (see for example Chen et al. (2007); Belloche & André (2004) who recover amplitudes from 0.1 to 10 km.s −1 .pc −1 in protostellar cores at scales of 5000 AU), but also with the recent results of Gaudel et al. (2020) probing inner regions of Class 0 envelopes, down to a hundred AU. The specific angular momentum j in our modeled core -computed from the velocity gradients following a similar methodology as the one extracting j from observations -has a roughly constant value of about 3 .10 −4 km.s −1 .pc between 10 2 AU and 10 3 AU (see bottom panel of Fig. 4). This order of magnitude is similar to typical specific angular momentum values recovered from observations in solar-type protostellar cores: see Yen et al. (2015a) for observations in 7 Class 0 protostars and the work of Gaudel et al. (2020) in 12 Class 0 protostars giving a specific angular momentum of 5 .10 −4 km.s −1 .pc below 1400 AU. The few observational constraints on the spatial distribution of specific angular momentum in prestellar cores and protostars suggest a scenario where local specific angular momentum is following a powerlaw R 1.6 in starless structures at scales larger than 6000 AU, while it is constant and conserved within collapsing star-forming cores (see Belloche (2013) for a review). However, Yen et al. (2015b) found that the decreasing trend of j(R) observed at large scales propagates down to radii smaller than 5000 AU. Gaudel et al. (2020) resolved the break radius around 1600 AU where it stabilizes with a weak dependence on radius j ∝ R 0.3±0.3 between 50 and 1600 AU. These observations are very well matched by the properties of specific angular momentum in our simulations with density perturbations and no rotation. Our model shows a dependency of Ω along a R −1.8 power-law, which translates ( j = R 2 Ω) as a j ∝ R 0.2 power-law, centered on the value of 3 .10 −4 km.s −1 .pc, which coincides with the recent observational results cited earlier.
However, the analysis of velocity gradients also reveals a large dispersion of the gradients directions between the disk and envelope scales. Depending on the viewing angle, the observed velocity gradients in the envelope can be very misaligned with respect to the small-scale gradient resulting from disk rotation. The projected kinematics can even produce a complete reversal of the gradient at some scales: for example in top panel of Fig. 4, the y projection (red curve) shows a gradient about 200 • from the disk scale gradient at scales 2000 AU. Since the rotation is generated locally in our model, it is not surprising that small and large scales gradients are unrelated and hence misaligned. In our model there is no large scale rotation, suggesting that observations of strongly misaligned gradients in protostellar cores could be due to non-axisymmetric gravitational collapse, rather than to global rotational motions of the core. Such scenario could explain, for example, the reversal of velocity gradient observed in the L1527 protostar (Tobin et al. 2011), or more recently our findings that the rotation of the small-scale disk in the low-luminosity protostar IRAM04191 (Maury et al. in prep) is opposite to the large-scale gradient at 2000 AU scales in the envelope, interpreted as core rotation (Belloche & André 2004). This observation of a "counter-rotating" disk is incompatible with simple axisymmetrical collapse models in which the disk would have simply formed because of the conservation of a large scale angular momentum in a rotating core. Tsukamoto et al. (2017) developed MHD models including Hall effect that form counter-rotating envelopes at the upper region of a pseudodisk under some conditions about the alignment between magnetic field and initial rotation. This thin counter-rotating layer is located close to the pseudo-disk, typically between 50 and 200 AU, and are thus unable to explain counter-rotation in the outer envelopes, at scales larger than 1000 AU, while it is a natural feature of our model.
The analysis of the radial and orthoradial part of the velocity and their contribution to the total velocity gradient directions shows that for small perturbation levels (typically less than 20%), the radial and total velocity gradients directions are very close over all the probed scales whereas the orthoradial component is misaligned. As the perturbation level increase, for scales larger than 400 AU the radial component is more and more misaligned whereas the orthoradial one tends to be closer to the total velocity gradient direction. These features depends on the choice of the line of sight. A common point which seems independent from the choice of a line of sight is that for small perturbation levels, the radial component is very close to the total velocity gradient direction. As the simulations start with no global rotation, this analysis shows that the notion of radial and orthoradial though as infall and rotation is less and less clear as the level of perturbations increases, which means as the deviation from axisymmetry grows. With a small perturbation level, the deviation from axisymmetry is low and the velocity gradients are mainly due to infall motions. With large perturbation levels, the deviation from axisymmetry is high and the velocity gradients are resulting from the complex dynamics of the self-gravitating gas.
When adding global rotation to the initial conditions, we found that the formed disks grow more rapidly and fragment earlier. Analysing the kinematics of the gas in those simulation with initial rotation, we show that the angular deviations are lower than in the case without initial rotation, which is less in agreement with Gaudel et al. (2020) observational results. From this analysis, the inferred specific angular momentum exhibits a slope and a mean value larger than in the case without initial rotation. As for the angular deviation, this result is less in agreement with observations than the case without initial rotation.

Conclusion
We have shown that protostellar disks can emerge from a nonaxisymmetrical gravitational collapse in which there is no rotation initially. We ran purely hydrodynamical simulations of a collapsing dense core starting from an initial condition where all the cells are at rest, and we broke the symmetry of the problem by adding density fluctuations over a flat profile. We showed analytically that the angular momentum in the frame of the accretion center and with respect to the accretion center is not a conserved quantity due to the non Galilean character of the frame. This leads to the possibility of the formation of a disk and we demonstrated numerically that a disk indeed forms in these simulations.
We then analysed our simulations from an observational point of view, computing the velocity gradients at different scales as it is done with real observations, and deducing the amplitude of these gradients and the specific angular momentum over the different scales. The results we obtained for the value of the specific angular momentum matches the ones from Belloche (2013) and Gaudel et al. (2020).
This study then suggests a new paradigm for the formation of protostellar disks but does not replace the current one -the conservation of pre-existing angular momentum at large scale. It shows that even if this conservation of angular momentum is absolutely correct in axisymmetrical model, it is more complicated when the axisymmetry is broken, as the accretion center frame becomes non Galilean. In the extreme case we studied in which every cells are initially at rest, it shows that even without initial rotation the formation of a protostellar disk is possible. The formation of these disks in our model does no longer depends on specificities of large scales, but results in a more generic process only due to the density fluctuations of the gas. We showed that the different features of the model based on the analysis of veloc-ity gradients are realistic. This model could thus help understand the features observed in some objects, as the angular deviation of velocity gradients at different scales. It also helps to understand the formation of a disk within progenitors that does not seems to contain enough angular momentum at large scales to form a disk by conservation of the angular momentum during the collapse. We found that these models without initial rotation are in better agreement with observational results than the models with initial solid-body rotation.
Our model uses minimal physical ingredients: hydrodynamics, thermodynamics, gravity, and density fluctuations. It is thus robust in the sense that all these ingredients are always present in dense cores. In particular, we found that even for 10% of density perturbations initially, a disk forms. Protostellar disks appear to be natural features of a gravitational collapse as soon as it is not axisymmetric. As large disks are formed in this new scenario, we stress the necessity to evoke processes extracting angular momentum to lead to the formation of smaller disks, as found with observations. We showed that the presence of magnetic field reduces the size of the formed disk, in agreement with past studies (Hennebelle et al. 2016;Maury et al. 2018). .1: Analysis of numerical errors in the simulations with different levels of perturbations ε, without sink particle. The legend is the same for both panels. To improve visibility, the results have been smoothed. Left: the relative difference between σ num and σ an . Right: the relative difference between the momentum in the disk |σ disk | and |∆σ| = |σ num − σ an |. We calculate the momentum in the disk only when the disk is not fragmented. The simulation with 10% of perturbations shows non negligible amount of numerical errors, whereas all the other simulations can be trusted.