Protostellar disks subject to infall: a one-dimensional inviscid model
and comparison with ALMA observations

Figure 2 is a schematic of the model and is best studied along with Figure 1. Figure 2a is suggestive of the actual situation. There is a disk-surface shock (colored blue) at height $H_{\mathrm{shock}}(r)$ which reduces the velocity component normal to it. This shock is followed in the normal direction by a thin cooling layer whose properties have been studied in detail (Neufeld & Hollenbach, 1994, hereafter NH94). The "post cooling-layer" (pcl) height is denoted $H_{\mathrm{pcl}}(r)$ and will be assumed to equal $H_{\mathrm{shock}}(r)$ owing to the thinness of the cooling layer. The shock surface terminates where the normal component of the infall velocity goes to zero. Due to finite disk height, this point is slightly outward of the centrifugal radius where $u_{z}\to 0$ at the midplane. At this point, the shock abuts a streamsurface (green dotted line) along which the flow is tangent by definition. Strictly speaking, the shock ends where the normal velocity becomes subsonic rather than zero; this distinction is ignored in what follows. Our model, shown in panel (b), idealizes the above picture to a razor-thin disk. The UCM infall flow (Appendices A & B) evaluated at the disk midplane, provides boundary conditions to the vertically integrated disk equations given in the next subsection. The infall into the disk has two parts: (i) For $r<r_{\mathrm{c}}$, there is both a vertical and radial component, infall being more vertical closer to the star and more radial, less vertical, and more dense as $r\to r_{\mathrm{c}}$. (ii) At $r=r_{\mathrm{max}}$, the outer radial boundary of the simulation, radial infall is applied as a boundary condition (Appendix C).

Assuming symmetry about the midplane, the vertically integrated density and pressure are

where $Z(r,t)$ is the upper disk surface defined as

Here $H_{\mathrm{shock}}^{+}>0$ denotes a height just above the disk-surface shock where infall quantities are specified. For $r>r_{\mathrm{c}}$, the surface $Z=S(r,t)$ is the streamsurface that joins the shock.

Performing a vertical integration of the conservation equations for mass, radial and angular momentum, we obtain

where each equation has been multiplied by $r$ to allow a finite-volume numerical treatment after averaging over a computational cell $[r_{i-1/2},r_{i+1/2}]$ (see §2.4). The terms with a ‘1’ subscript are pre-shock quantities that arise from $z$-integration of $\partial/\partial z$ terms from $z=0$ to $Z(r)$. Appendix E discusses a subtlety in the vertical integration that can be ignored for the present assumption of a razor thin disk, but should be kept in mind when finite disk thickness is accounted for.

Note that in its usual form, the right-hand-side of the radial momentum equation (multiplied by $r$) has a $-r\partial\mathcal{P}/\partial r$ term. We have written this as

and moved the $-\partial\left(r\mathcal{P}\right)/\partial r$ to the left hand side, leaving a $\mathcal{P}$ on the right hand side of (9). This was done in order that radial fluxes have the same form as for Cartesian coordinates which permits direct use of shock-capturing methods developed for Cartesian coordinates.

A key assumption of the vertical integration procedure is that the incoming mass and momenta (and below, internal energy) are instantaneously mixed vertically into the disk. Also, note that turbulent viscosity has not been included in the above equations.

The equation of state for internal energy (per unit volume) is taken to be $e=p/(\gamma-1)$, where $\gamma\equiv c_{p}/c_{v}$. For typical temperatures in our simulation, which range from $T=10$ and 100 K, $\gamma$ varies between $5/3$ and $7/5$ as the rotational mode of H₂ is excited. For simplicity we take $\gamma=7/5$.

For axisymmetric flow, the internal energy satisfies

Since vertical compressional heating $-p\partial u_{z}/\partial z$ involves the product of a Heaviside and $\delta$ function at the disk-surface shock and there is a cooling layer following the shock, we take a slightly different approach to vertically integrating (12) than was done for the mass and momentum equations. The cooling layer is thin: using the NH94 code we find that its column density varies between $0.002$ and $0.005$ gm cm^-2 for the simulation parameters. This allows one to perform the integration up to $H_{\mathrm{pcl}}$, the post cooling-layer (subscript pcl) height with negligible loss of consistency with the shock height used for integrating the mass and momentum equations. This avoids having to deal with the pressure compression term at the shock. Making the assumption that the pressure compression term $-p\partial u_{z}/\partial z$ can be neglected within the disk compared to the radial pressure compression, vertical integration of (12) gives

where

For the second term on the left-hand-side of (13) we have used the fact that the mass flux $\rho u_{z}$ is preserved from the pre-shock state to the end of the post-shock cooling layer.

The shock is effectively a jump and immediately across it, the gas temperature rises to a large value. For instance, for our simulation parameters we have at $r/r_{\mathrm{c}}=0.2$ a pre-shock velocity of $(u_{z})_{1}=4.64$ km s^-1 and a density of $n_{\mathrm{H}}=6.6\times 10^{7}$ (see Figure 13). For a pre-shock temperature of $T=20$ K, the NH94 code then gives $T=1066$ K immediately behind the shock. The dominant emission in the cooling layer following the shock is from rotational transitions of H₂O and CO (see Figure 9, top left panel in NH94). However, these details are not important for specifying the temperature at the end of the cooling layer. Given that at the end of the cooling layer, gas and grain temperatures have equilibrated, we set $T_{\mathrm{pcl}}$ equal to an estimate for the minimum grain temperature provided by NH4 (their Equation 43) which we slightly modify by adding a term for the flux from the surrounding cloud:

Equation (15) assumes that all of the vertical kinetic energy flux is annulled across the shock and radiated. The above neglects pre-heating of the gas by stellar and shock irradiation (Chick & Cassen, 1997).

For radiative cooling $\mathcal{Q}_{\mathrm{cool}}$ and stellar heating $\mathcal{Q}_{\mathrm{heat}}$, the formulation due to Nakamoto & Nakagawa (1994) is used:

where

is the vertical optical depth. The Planck mean opacity, $\kappa_{\mathrm{P}}$, was obtained using the subroutines of Semenov et al. (2003). The opacities, given in cm² per gm of gas, are functions of gas density and temperature and for these we use the characteristic values

with $H=c_{\mathrm{s}}/\Omega$, $c_{\mathrm{s}}$ being the sound-speed. $F_{\mathrm{irr}}$ is the stellar flux and is given by

The factor of two dividing $F_{\mathrm{irr}}(r)$ in (16) accounts for the fact that half of the flux enters the disk while the other half is radiated to space.

The quantity $\gamma_{\mathrm{irr}}$ is the angle between a ray incident on a point $P=(r,z=H(r))$ on the scale height surface and the normal to the surface. Consulting Figure 3 we have that $\gamma_{\mathrm{irr}}=\theta_{\mathrm{i}}+\theta_{\mathrm{n}}$, where

are the inclination angles, relative to the midplane, of the incident ray and the surface normal, respectively. This treatment assumes that the radiation in each ray is absorbed at point $P$ where the ray meets the surface $(r,H(r))$. In reality, the absorption will take place over a distributed region. Stellar heating is applied to only those parts of the surface that face the star. It was found that when the numerical scale height surface was used, large oscillations resulted because star facing portions of the surface receive stellar irradiation, while non-star-facing portions receiving none. This increases the amplitude of oscillations. To avoid this, a quartic polynomial is fit to the scale height surface for the purposes of calculating $F_{\mathrm{irr}}$. The accuracy of the fit is assessed a posteriori; see the dashed line in Figure 5b. We use $L_{*}=2L_{\odot}$ which is the bolometric luminosity for L1527 IRS (Aso et al., 2017).

The transport equations (8)–(10) and (13) have the form

where

is a column vector of evolved quantities. The vector $\mathbf{F}(\mathbf{q})$ contains the radial advective fluxes of $\mathbf{q}$, and $\mathbf{S}(\mathbf{q})$ contains the rest of the terms.

The computational domain $[r_{\mathrm{min}},r_{\mathrm{max}}]$ is divided into $N$ equal cells of width $\Delta r$. Let a bar denote a cell average. Then averaging (21) over the $i$th cell gives:

where $\mathbf{H}_{i+\frac{1}{2}}$ is a numerical flux evaluated at a cell face. The numerical fluxes are evaluated using the second-order Kurganov & Tadmor (2000) scheme, specifically their equation (4.4) . The third-order TVD (total variation diminishing) Runge-Kutta scheme is used for time advancement. At $r_{\mathrm{min}}$, a one-way “diode” boundary condition is applied: if there is inflow into the domain ($u_{r}>0$) at the end of a sub-step, then we set $u_{r}=0$. At $r_{\mathrm{max}}$, radial infall is specified as described in Appendix C.

Simulation parameters (Table 2) were chosen to approximate L1527 IRS. According to Aso et al. (2017), $u_{\phi}=2.31$ km s^-1 at $r=74$ au (corrected for beam resolution) where the transition from a $j_{z}$-preserving $u_{\phi}$ to a Keplerian $u_{\phi}$ takes place. This implies that $M=0.45M_{\odot}$, which is fixed throughout the simulation. We choose the evolutionary epoch to be $t_{0}=10^{5}$ yr. Then, to obtain $M=0.45M_{\odot}$, (83) gives a cloud temperature of $T=20.1175$ K. From the discussion in the introduction, we know if Keplerian $u_{\phi}$ occurs at a certain radius in a $j_{z}$-preserving region, then that radius equals $r_{\mathrm{c}}$. Hence we want $r_{\mathrm{c}}=74$ au. This is achieved from (84) by setting the cloud rotation to be $\Omega_{0}=1.5015\times 10^{-13}$ s^-1.

The simulation is initialized with a very small surface density and zero velocities. Infall is then turned on. Time elapsed from the start of the simulation is denoted as $\tau$. To avoid very large imbalances between heating and cooling at the start, the pressure dilatation and radiation terms in the energy equation are ramped up for 3,000 yr to their full values.

How should one interpret the simulation time $\tau$? Since an infall field appropriate to a star of $M=0.45M_{\odot}$ is turned on at $\tau=0$ with a non-existent disk, $\tau$ cannot be interpreted as time since the start of collapse. If the simulation achieved a stationary state as $\tau\to\infty$, one might claim that the stationary state represented the true state of the disk at epoch $t_{0}$. Unfortunately, this is not the case here: the disk mass increases with $\tau$. The claim cannot be made even if the simulation arrived at a stationary state. This is because the true state of a disk at $t_{0}$ is very likely dependent on its past history, and our simulation does not trace the same history.

The best we can do is make a comparison at a $\tau$ value when the simulated disk mass is comparable to the observed disk mass. Nakatani et al. (2020, pg. 5) estimate that $M_{\mathrm{disk}}=0.26M_{\odot}$ for L1527 IRS. In our simulation this is achieved at $\tau\approx 40,000$ yr (Figure 5f). Fortunately, simulation velocities and temperature are insensitive to $\tau$ and disk mass.

Figure 4 shows the first set of radial profiles of various quantities at different times.

The radial velocity $u_{r}$ (Figure 4a) most clearly shows a radial shock at $r/r_{\mathrm{c}}\sim 1.5$, i.e, outward of the centrifugal radius, which propagates outward very slowly with decreasing speed. The shock reduces $|u_{r}|$ to a small value and $|u_{r}|$ further decreases as the centrifugal radius is approached. It was verified that the region $r>r_{\mathrm{shock}}$ (the IRE) is well described by zero mechanical energy ballistic flow. As stated in the introduction, the ballistic approximation breaks down inward of the radial shock since $u_{r}$ is subsonic there.

The mechanism of shock initiation ($\tau<1000$ yr) may be qualitatively described as follows. Consider a shock-free state where $u_{r}$ slows smoothly with decreasing $r$ due to angular velocity build-up, and eventually becomes subsonic. In the region where $u_{r}$ is subsonic, pressure waves can travel upstream up to the point where $u_{r}$ is sonic. A shock forms as waves pile-up at the sonic point. As the strength of the shock grows, it propagates into more supersonic flow until the Rankine-Hugoniot jump conditions (including radiative losses) for a steady shock are satisfied. To our knowledge, no simple analysis can predict where the shock will finally sit.

A better quantity to assess the radial flow in the post-shock, low $|u_{r}|$ region is the radial mass flux $\dot{M}(r)\equiv 2\pi r\Sigma u_{r}$. Figure 4b shows that $\dot{M}\approx-2.7\times 10^{-6}M_{\odot}\mathrm{yr}^{-1}$ in the IRE beyond the radial shock, and is conserved across the shock as it should be. From the shock inward it decreases to a small valued cusp at the centrifugal radius. In this region therefore, the surface density must build up with time, which it does and leads to a maximum in $\Sigma$ at $r/r_{\mathrm{c}}=1$. The fact that $\dot{M}(r)$ is small at $r=r_{\mathrm{c}}$ means that the region $r_{\mathrm{c}}>1$ feeds a relatively small amount of mass to the region $r_{\mathrm{c}}<1$, which accumulates mass via vertical infall.

For $0.3r_{\mathrm{c}}<r<0.9r_{\mathrm{c}}$, $\dot{M}(r)$ is relatively uniform at $\approx-0.8\times 10^{-5}M_{\odot}\mathrm{yr}^{-1}$. The dashed line shows the result of the inviscid formula (122) following Cassen & Moosman (1981, hereafter CM81) which accounts for the drag exerted by the sub-Keplerian infall on a Keplerian disk. Its agreement with the simulation is only fair. Much better agreement results (see Figure 4c) if we substitute the actual disk angular momentum, $\Gamma(r)\equiv j_{z}=u_{\phi}r$, into (118), which is valid for general but steady $\Gamma(r)$. The remaining differences with the simulation are attributed to unsteadiness of $\Gamma(r)$ and to the fact that the ratio in (118) becomes indeterminate as $r/r_{\mathrm{c}}\to 1$.

Figure 4d shows that $u_{\phi}\propto r^{-1}$ (angular-momentum preserving) for $r>r_{\mathrm{c}}$ and is nearly Keplerian for $r<r_{\mathrm{c}}$, however, as noted above there is sufficient deviation from Keplerian $u_{\phi}$ to account for a qualitative change in the disk $\dot{M}$.

Figure 4e shows the corresponding build-up of surface density. Its radial profile has a minimum at $r/r_{\mathrm{c}}\approx 0.6$, then a maximum at $r_{\mathrm{c}}\approx 1$ before decreasing exponentially for $r_{\mathrm{c}}<r<r_{\mathrm{shock}}$.

Motivated by its presentation in a model fit to observations (Ohashi et al., 2014) to be discussed later, Figure 4f shows the midplane number density $n_{\mathrm{H}_{2}}$ estimated using

which assumes a Gaussian (hydrostatically balanced) density profile for $z\in[-\infty,\infty]$. Here $m_{\mathrm{H}_{2}}$ is the molecular mass of H₂. A more accurate result for the region $r_{\mathrm{c}}<1$ which accounts for the fact that the Gaussian is truncated by the surface shock is given later and shows that the estimate represents an under-prediction, at least for $r_{\mathrm{c}}<1$.

At the inner boundary, where a one-way diode boundary condition is applied ($u_{r}$ is set to zero if it is incoming), $u_{r}$ alternates between periods of low outflow and zero flow. This would obviously not occur in reality. Since the outflow is subsonic, there is in reality one incoming characteristic wave (with speed $u_{r}+c_{\mathrm{s}}>0$) carrying information from the region $r<r_{\mathrm{min}}$ not in the computational domain. It is unclear how lack of this information in the simulation affects the disk within the computational domain. Since $\dot{M}(r)$ has a relatively long uniform region away from the computational boundary, one might conjecture that this should continue were it not for the computational boundary. If it did continue, the rate of loss ($\approx 0.09\times 10^{-5}M_{\odot}\mathrm{yr}^{-1}$) is small compared rate of growth of disk mass ($\approx 0.6\times 10^{-5}M_{\odot}\mathrm{yr}^{-1}$) and would diminish the rate of disk growth only slightly. Recall that the CM81 inviscid result that Keplerian $u_{\phi}$ produces a zero $\dot{M}(r)$ at $r=0$. If the conjecture is true that $\dot{M}(r)$ should remain uniform outside the computational domain, then this can only be due to a small deviation from Keplerian angular momentum. In reality, turbulent viscosity and outflows must also be taken into account.

Figure 5 presents the final set of diagnostics starting with the temperature in panel (a). The temperature has a sharp spike to $\sim 35$ K due to the radial shock, with a corresponding change in scale height $H(r)$ shown in panel (b). This temperature rise is insufficient to account for SO desorption in the ALMA observations and will be discussed in more detail in §4.4 by appealing to non-LTE (non local thermodynamic equilibrium) processes in the shock which are not captured in the simulation.

Figure 5c shows the Toomre parameter

where

is the epicyclic frequency for arbitrary rotation profiles. $Q<1$ implies gravitational instability (GI) to axisymmetric modes (assuming a purely rotational flow) and spiral features occur for $Q\lesssim 1.7$. The simulation has $Q<1$ at all $\tau$ in the region $1<r_{\mathrm{c}}<r_{\mathrm{shock}}$ where $u_{r}$ is small and $u_{\phi}\propto r^{-1}$. Axisymmetric GI is therefore possible in this region. For $\tau>10000$ yr a growing portion of the region $r<r_{\mathrm{c}}$ also develops $Q<1$. In this region $u_{r}$ is even smaller which renders the Toomre criterion valid. The sharp dip in $Q$ near the inner boundary is an artifact of the boundary condition; this was confirmed by running a simulation with $r_{\mathrm{min}}/r_{\mathrm{c}}=0.15$ instead of $r_{\mathrm{min}}/r_{\mathrm{c}}=0.20$.

Figure 5d shows the disk-surface shock height $H_{\mathrm{shock}}(r)$ calculated as described in Appendix H. This calculation extends radially up to end of the shock where $u_{z}$ at the midplane becomes subsonic. The figure shows that the ratio $H_{\mathrm{shock}}(r)/H(r)$ is always $>1$ and increases with $\tau$. Since the pre-shock density and vertical speed are constant with $\tau$, the density ratio across the disk-surface shock is also constant. On the other hand, the disk density increases with $\tau$. Hence, the height of the disk-surface shock must increase with $\tau$.

Figure 5e shows that the corresponding midplane number density, obtained from (129), is larger than the estimate plotted in Figure 4b based on the assumption of a Gaussian density profile not truncated by a shock.

Finally, in order to assess the possibility of shear driven turbulence, Figure 6 plots the average vertical shear

normalized by the local orbital frequency $\Omega(r)$ for $u\to u_{r}$ and $u_{\phi}$. Since these components of the infall velocity are continuous across the surface shock, equation (27) represents the vertical shear within the disk interior. One observes that the shear for the $r$-component is everywhere larger than for the $\phi$-component and its negative sign means that the negative $u_{r}$ velocity at the top of the disk is more negative than at the midplane. Therefore $\partial u_{r}/\partial u_{z}$ shear driven turbulence remains possibility.

The negative sign of the $\phi$-component reflects the fact that $u_{\phi}$ is sub-Keplerian at the top of the disk and nearly Keplerian at the midplane.

Sa14a detect a region of high SO emission suggestive of a ring at $r=100$ au; see Figure 1d in Sa17. It is thought that a radial shock is present at this location that sublimates SO from icy grains. Figure 5a shows that the temperature rises to $\approx 35$ K after the radial shock and quickly cools back down 30 K; this is insufficient for SO sublimation whose required temperature is 50–60 K (Sa17). Using the radex code, Sa17 infer a post-shock temperature of $194(^{+146}_{-60})$ K and an envelope temperature of $29(^{+26}_{-11})$ K (their pg. L80). To account for SO sublimation, in §4.4 we will follow Aota et al. (2015) and consider non-LTE effects in a thin post-shock layer.

In the simulation, after the shock, $T$ drops to the pre-shock temperature which is consistent with the freeze-out of SO inferred from the observations. Further inward, the temperature rises as the star is approached which is consistent with the reappearance of SO emission in observations.

Finally, and crucially, Sa14a place the radial shock at the centrifugal barrier which is at half the centrifugal radius: $r_{\mathrm{shock}}=0.5r_{\mathrm{c}}$. On the other hand, the simulation gives $r_{\mathrm{shock}}\approx 1.5r_{\mathrm{c}}$. This discrepancy is addressed in §§4.5 and 4.6.

Ohashi et al. (2014, hereafter O14) obtained a transition from an angular momentum preserving $u_{\phi}$ to a Keplerian $u_{\phi}$ at $r=54$ au; following our previous reasoning, this is also the location of the centrifugal radius. Figure 7, reproduced from O14, shows the result of their Model 2 constructed to fit their observations. This model fits their line profiles much better than their Model 1. Figure 7 shows a radial shock located at $r_{\mathrm{shock}}=120$ au, i.e., outward of the centrifugal radius, in agreement with our simulations. Our pre-shock value of $u_{r}\approx-2$ km s^-1 is close to the pre-shock value in Figure 7. One difference is that in the O14 model, the midplane density jump does not coincide with the shock (as it should physically) but occurs at $r=r_{\mathrm{c}}=54$ au. In the O14 model, $\log_{10}n_{\mathrm{H}_{2}}$ jumps from 7.4 to 9.8 across the shock while in the simulations it jumps at the shock from $\approx 6.5$ to 8 at $\tau=40,000$ years. The O14 post-jump value is actually closer to our peak value which occurs at $r=r_{\mathrm{c}}$.

From an LVG (Large Velocity Gradient) analysis, O14 conclude that the gas temperature in the region of SO emission is $\approx 32$ K. This is comparable to the post-shock temperature the simulation. Since this is smaller than the SO sublimation temperature of 40–60 K, O14 suggested (following Aota et al. (2015)) that SO sublimation occurs in a thin layer behind the shock.

In ALMA cycle 1, (Aso et al., 2017, hereafter A17) (see their Figure 5) identified a region of $u_{\phi}\sim r^{-1.22}$, i.e., approximately $u_{\phi}\sim r^{-1}$ which is $j_{z}$ preserving. This transitions to a Keplerian $r^{-0.5}$ rotational velocity profile at $r\approx 74$ au (corrected for beam resolution). As mentioned earlier, the location of Keplerian angular velocity in a $j_{z}$ preserving region gives the centrifugal radius. Therefore, $r_{\mathrm{c}}=74$ au. This is in line with our simulations which indicate a transition from $u_{\phi}\propto r^{-1}$ to $u_{\phi}\propto r^{-1/2}$ at $r=r_{\mathrm{c}}$; see Figure 4d.

Like O14, A14 obtain a density jump in their best fit model at a location that is close to the edge of the Keplerian disk. Specifically, the edge of the Keplerian disk is at $r=74$ au while the density jump is at $r=84^{+16}_{-24}$ au.

The issue of SO desorption due to the putative radial shock of L1527 was taken up by Aota et al. (2015, henceforth Ao15). They performed 1D non-LTE calculations of shock structure accounting for gas cooling by collisions with dust, radiative emission from the dust, and SO desorption via thermal heating and sputtering. They concluded that for our value of 2 km s^-1 for the pre-shock velocity, the required gas density for SO desorption is $n_{\mathrm{H}}\gtrsim 10^{9}$ cm^-3; see their Figure 7. On the other hand, the (midplane) pre-shock density in our case is $n_{\mathrm{H}}=6.3\times 10^{7}$ cm^-3, a factor 16 smaller. In the following discussion we shall assume that SO desorption nevertheless occurs with the parameters given by the simulation, and then go on to calculate the SO column density $N_{\mathrm{SO}}$ for comparison with the ALMA observations of S14a.

The NH94 shock code does not predict desorption but otherwise can be used to calculate 1D shock structure. We therefore ran the NH94 code to obtain the gas and dust temperature in a 1D shock with the pre-shock velocity and density given above, together with a pre-shock temperature of 25 K. The result is shown in Figure 8. The abscissa is the $N_{\mathrm{H}}$ column density measured from the shock front ($r_{\mathrm{shock}}$) and is defined as

The column density $N_{\mathrm{H}}(r)$ is a proxy for distance behind the shock. One observes from Figure 8 that $T_{\mathrm{gas}}\geq 100$ K for a column density of $N_{\mathrm{H}}=2.3\times 10^{20}$ cm^-2 behind the shock and $T_{\mathrm{dust}}\geq 50$ K for a column density of $N_{\mathrm{H}}=3.1\times 10^{20}$ cm^-2. This dust temperature is in the middle of the range 40–60 K for SO sublimation (e.g. Ohashi et al., 2014, and the references therein). The above values are comparable to those obtained in the non-LTE calculations of Ao15; see their Figures 3 and 4 for the gas and dust temperatures, respectively, and Figure 8 for the column density of warm gas. Our post-shock density $n_{\mathrm{H}}=8\times 10^{8}$ cm^-3 together with a warm dust column density of $N_{\mathrm{H}}=3.1\times 10^{20}$ cm^-2 implies that the length $\ell_{\mathrm{warm\ dust}}=0.026$ au.

Next, we use the time for SO to re-adsorb on to grains as given by Ao15 (their pg. 8, right column):

From this we can obtain the distance over which re-adsorption occurs as

where $|u_{r}|_{\mathrm{post}}=0.05$ km s^-1 is the post-shock velocity and $n_{\mathrm{H}}=8\times 10^{8}$ cm^-3 which gives $\ell_{\mathrm{adsorb}}=0.13$ au. Hence the total radial length of the SO region is

Assuming an SO abundance of $10^{-7}$ following Ao15, (31) implies an SO column density of

On the other hand, S14a state (last paragraph of their methods section) that the “column density of SO is well constrained to be $4\times 10^{14}$ cm^-2 for the central $1^{\prime\prime}\times 2^{\prime\prime}$ region”, which is comparable to our value. Note that our column density is along the radial direction whereas the S14a value is along the line of sight which is purely radial only along the central line of sight. Also note that if the observation measures both the region in front and behind the star, then our value should be multiplied by two for comparison with the observation. This brings the two values into excellent agreement, which is fortuitous given the many uncertainties.

Solving the midplane (ballistic) parabolic orbit equation (2) for $u_{r}$ and determining its maximum value, Sa14a (their Equation 6) obtain that

We note for future use that $|u_{r}^{\mathrm{max}}|$ occurs at the centrifugal radius, $r_{\mathrm{c}}$, and drops to zero at $r=r_{\mathrm{CB}}=0.5r_{\mathrm{c}}$.

Equation (34), together with the fact that $u_{\phi}^{\mathrm{max}}$ occurs at the centrifugal barrier, allowed Sa14a to determine the location of the centrifugal barrier from a position velocity (PV) plot of the envelope tracer cyclic-C₃H₂ (their Figure 2c). Specifically, along the stellar position, the line of sight velocity is purely radial for the nearly edge-on L1527 disk. From the PV plot, we have $u_{r}^{\mathrm{max}}=1$ km s^-1 (line B on the red-shifted side). On the other hand $u_{\phi}^{\mathrm{max}}$ in the envelope is seen to be $\approx 2$ km s^-1 at $r=100$ AU. Hence, relation (34) holds and one concludes that $r_{\mathrm{CB}}=100$ au. This coincides with the location of the radial shock (marked by SO emission) and justifies conflating the radial shock with a centrifugal barrier. Furthermore, given the relation $r_{\mathrm{c}}=2r_{\mathrm{CB}}$, Sa14a conclude that $r_{\mathrm{c}}=200$ au. Finally, from this value they conclude that the stellar mass is $M=(0.18\pm 0.05)M_{\odot}$, which is lower than the value $M=0.45M_{\odot}$ obtained by Aso et al. (2017) from the Keplerian velocity.

The Sa14a argument is indeed reasonable and cogent. However, it should be noted that since the radial Mach number, $M_{r}$, upstream of the radial shock must be supersonic, the radial shock cannot exactly coincide with the centrifugal barrier where $u_{r}=0$ by definition. The values of $r_{\mathrm{CB}}$ and $r_{\mathrm{shock}}$ can be close, however, provided that $M_{r}$ is large enough to produce the observed post-shock temperature.

It should be pointed out that if we restrict ourselves to our envelope, i.e., the region outward of our radial shock (given that c-C₃H₂ is an envelope tracer), then our analog of (34) is

which is quite different from (34). This is because our envelope is located so far outward that it does not access the higher $u_{\phi}$ region in the inner part of the disk. The ratio $1.55$ was obtained from Figures 4a and (d) at $\tau=50,000$ yr. From panel (a), we have $|u_{r}^{\mathrm{max}}|=2.12$ km s^-1at $r/r_{\mathrm{c}}=1.67$, just outward of the shock which is the end of the IRE. From panel (d), from the same location we read-off that $u_{\phi}^{\mathrm{max}}=1.37$ km s^-1. Their ratio is 1.55.

The next sub-section shows that even when we perform a simulation with initial conditions set to the Sa14a ballistic flow and allow only radial infall, the radial shock still positions itself at $r/r_{\mathrm{c}}\sim 1.5$.

Only radial infall is imposed and the initial velocity field is the Sa14a zero energy flow outside the centrifugal barrier ($r>r_{\mathrm{c}}/2):$

For $r\leq r_{\mathrm{c}}/2$ we would like to impose Keplerian angular velocity but it does not match (37) at the interface. To eliminate the discontinuity, a flat bridge in $u_{\phi}(r)$ is introduced to connect the two profiles:

Instead of this ungainly equation, it is better to consult the black curves in Figures 9b and 9d. Initially $\Sigma=0.1$ gm cm^-2, and $T=T_{0}=20.1175$ K.

We note that from their observations of the protostar Barnard 5 IRS1, Velusamy & Langer (1998) suggest that an outflow can block infall from a biconical (or biparaboloidal) region such that most of the infall is radial. This can also be observed in magnetic collapse simulations (e.g., Zhao et al., 2016). The present simulation is also intended to mimic this situation, however, angular momentum transport by the magnetic field is not considered.

Figure 9a shows the Mach number, $(M_{r})_{\mathrm{ballistic}}$, in the S14a ballistic flow which was used as the initial condition. The interesting feature is that $(M_{r})_{\mathrm{ballistic}}$ is quite high close to the centrifugal barrier at $r_{\mathrm{CB}}=0.5r_{\mathrm{c}}$, and therefore a strong shock could plausibly sit close to $r_{\mathrm{CB}}$. For instance, $(M_{r})_{\mathrm{ballistic}}=-6$ at $r/r_{\mathrm{c}}=0.64$.

Figure 9b shows the radial velocity. Because in the initial condition, the inner disk has $u_{r}=0$, pressure waves propagate upstream in the subsonic region, and communicate to the oncoming flow the fact that $u_{r}$ has stagnated. As a result, a shock is produced which greatly reduces $|u_{r}|$ behind it. The shock propagates rapidly upstream up to $r/r_{\mathrm{c}}\sim 1.5$, after which it propagates slowly. This is similar to what we observed in the simulation with both vertical and radial infall.

A better view of the radial motion is obtained from the disk mass flow rate, $\dot{M}$; see panel (c) and focus on times after the flow has settled down ($\tau\geq 30,000$ yr). $\dot{M}(r)$ is continuous across the shock as it should be and gradually reduces to zero (due to absence of Cassen-Moosman drag) rather than to a finite value as in the previous simulation. As a result, the stellar mass accretion rate is zero.

Panel (d) shows that in the region of the initial $u_{\phi}$ bridge, the angular velocity changes to being nearly Keplerian The transition to $u_{\phi}\sim 1/r$ takes place at $r/r_{\mathrm{c}}\approx 0.5$ rather than at $r/r_{\mathrm{c}}\approx 1$ in the previous simulation. We conclude that the larger region of Keplerian $u_{\phi}$ in the previous simulation is aided by the sub-Keplerian angular momentum of the vertical infall.

Panel (e) shows that the surface density peaks at $r/r_{\mathrm{c}}=1$. This coincides with the location of inflection point in $\dot{M}(r)$ as expected from the continuity equation

A local maximum with respect to $r$ in the first term implies that $\partial_{r}^{2}\dot{M}(r)=0$. Due to the absence of both vertical infall and Cassen-Moosman drag, the inner disk is devoid of mass. This would be quite different if magnetic and/or viscous torques were present.

We will skip showing the temperature profile since it is very similar to the previous simulation.

Finally, Toomre $Q$ in panel (f) shows that the region $r_{\mathrm{c}}\gtrsim 0.75$ is susceptible to GI. The Toomre criterion will be valid inward of the radial shock since there is very little radial infall in this region.

The main conclusion is that our attempt to mimic the Sa14a ballistic flow failed to place the radial shock anywhere near the centrifugal barrier, $r_{\mathrm{CB}}=0.5r_{\mathrm{c}}$.

Finally, we would like to raise the possibility that with the addition of an appropriate physical effect that we have not included, the radial shock can position itself on the left side of the maximum ballistic $|M_{r}|$ in Figure 9, i.e, at $r_{\mathrm{shock}}<r_{\mathrm{c}}$ as in the S14a observations. We encourage investigators to consider various possibilities.

Cassen & Moosman (1981, CM81) pioneered the subject. They adopt Ulrich’s (1976) infall flow field (described in Appendix A and B) and use Shu’s (1977) inside-out collapse solution to set the initial radius (in the parent cloud core) of infalling particles as well as the mass infall rate, $\dot{M}_{0}$. CM81 apply the infall flow field evaluated at the midplane as a boundary condition to the vertically integrated thin disk equations. The angular velocity is assumed to be Keplerian and the radial momentum equation then provides an explicit solution for $u_{r}$. This is then substituted into the transport equation for surface density, leading to only one transport equation that needs to be solved. With the exception of Stahler et al. (1994), a similar procedure is applied in all the other works we discuss below. Recall that in the present model, we do not assume that $u_{\phi}$ is Keplerian and retain all relevant transport equations. We do not include a turbulent viscosity because, as stated earlier, our goal is to first understand the basic state and the turbulence generating mechanisms it supports.

In CM81, solutions are obtained for both the inviscid and viscous cases with uniform turbulent viscosity. The Keplerian assumption allows a local and analytical determination of the mass flow rate $\dot{M}$ through the disk. By local we mean that one does not need to solve a boundary value problem. CM81 were the first to point out the existence of a star-ward accretion flow even in the absence of viscosity; this arises because the sub-Keplerian angular momentum of the infall exerts a drag on the Keplerian flow in the disk.

Lin & Pringle (1990, hereafter LP90) have two transport equations. The first is a diffusion equation for surface density as in accretion disk theory without infall (Lynden-Bell & Pringle, 1974). This diffusion equation is obtained by inserting the following expression for radial velocity:

into the mass conservation equation. Here $\nu_{\mathrm{t}}$ is the turbulent viscosity, $\Gamma\equiv u_{\phi}r$ is the specific angular momentum, and a prime denotes differentiation with respect to $r$. Equation (42) is obtained from the angular momentum equation by (i) assuming that $\Gamma$ is time invariant or at least slowly varying compared to the time scale for viscous/turbulent diffusion; and (ii) neglecting Cassen-Moosman drag. To specify $\Gamma(r)$, LP90 assume a modified Keplerian balance that accounts for self-gravity of the disk. To specify $\nu_{\mathrm{t}}$, LP90 include the treatment of Shakura & Sunyaev (1973) as well as self-gravity contributions. The second transport equation in the LP90 model is the energy equation for temperature.

Hueso & Guillot (2005, hereafter HG05) adopt a model with a single transport equation, that for the surface density $\Sigma$. For the radial velocity that this equation requires, they specialize (42) to Keplerian $u_{\phi}$:

For the mass source term due to infall, HG05 do not use UCM infall flow field. Instead, they develop a source term based on the assumption that cloud material with a given angular momentum ends up in the disk at the radius where the Keplerian disk has the same angular momentum. Yang & Ciesla (2012), whose interest is in the distribution of refractory materials, adopt a similar model for disk dynamics.

The models of Visser et al. (2009, hereafter V09) and Visser & Dullemond (2010) are similar to HG05 but have some novel features. (i) The mass source term from infall is evaluated at the disk surface rather than at the midplane as in all other works, including ours. V09 define the disk surface to be where the disk density (assumed to be in hydrostatic balance) matches the density of the infall field. (ii) If the infall trajectory to a point on the boundary is obstructed by an outer part of the disk, V09 raise that point until it is no longer obstructed.

The set-up of the Stahler et al. (1994, hereafter St94) analysis is the most similar to ours: it is also inviscid, vertically integrated, and prescribes the UCM infall boundary condition. However, it does not include radial infall and assumes pressure-free steady flow. For their outer disk, St94 obtain a singularity in surface density at $r/r_{\mathrm{c}}=0.345$ where $u_{r}\to 0$, which they interpret as a dense ring. They explain (their pg. 348, second column) the mechanism of its formation as being similar to that of a centrifugal barrier, namely, increase of centripetal acceleration leads to a stagnation of $u_{r}$.

Note that the St94 centrifugal barrier is further inward than the Sa14a barrier. This is because vertical infall brings in inward radial momentum, and because Cassen-Moosman drag reduces $u_{\phi}$ which reduces the ability of a fluid element to oppose the pull of gravity.

This sub-section demonstrates that the singularity/dense ring is an artifact of their pressure-free assumption, which is equivalent to assuming that $u_{r}$ has infinite Mach number; this prevents upstream propagation of information. When finite pressure is allowed, $|u_{r}|$ will necessarily transition from being supersonic to subsonic before the stagnation point. This permits information to travel upstream and invalidates the one-way pressure-free integration procedure in St94. Intuitively, when a subsonic flow region (in an otherwise supersonic flow) senses a stagnation of $u_{r}$, it sends pressure waves upstream to halt the flow via a shock rather than allow all the mass to pile-up at the stagnation point. We show below that this is what happens when the pressure-free assumption is relaxed. A shock is created that rapidly propagates radially outward, greatly reducing $|u_{r}|$ behind it. The azimuthal velocity becomes nearly Keplerian behind the shock.

Let us begin by showing that we reproduce St94 when their set-up is mimicked. This is achieved by running isothermally with a uniform temperature of $T=1$ K to approximate the pressure-free assumption, and placing the inner boundary of the computational domain at $r_{\mathrm{min}}=0.355$, i.e., a little outward of the singularity where $|u_{r}|$ is still supersonic (the exit Mach number was found to be 17). The rest of the parameters are the same as for the L1527 simulation.

A steady state is reached quickly and Figure 10 shows that the agreement between the exact and computed solution is excellent except for small numerical oscillations in $\Sigma$ where the gradient becomes large. The ordinates are normalized in the same way as St94. The normalized surface density is

where $u_{\mathrm{K}}(r_{\mathrm{c}})$ is the Keplerian speed at $r=r_{\mathrm{c}}$. The specific angular momentum $\Gamma\equiv ru_{\phi}$ is normalized by the local Keplerian value $\Gamma_{\mathrm{K}}(r)$, and the radial velocity is normalized by $u_{\mathrm{K}}(r_{\mathrm{c}})$. The angular momentum $\Gamma$ is super-Keplerian. Note that the surface density $\Sigma\to\infty$ as $r/r_{\mathrm{c}}\to 0.345$ where $u_{r}\to 0$.

Next, the pressure-free assumption is relaxed by setting $T=20$ K in which case $u_{r}$ becomes subsonic in the computational domain before the singularity in reached. Figure 11 shows that a shock is created which propagates rapidly outward causing the solution to peel away from the St94 result. Behind the shock, the radial velocity is greatly reduced and $u_{\phi}$ becomes nearly Keplerian. If there were radial infall, the shock would come to rest where its propagation speed equals the radial infall speed.

Finally, for their inner disk ($r/r_{\mathrm{c}}<0.355$), St94 assume that the azimuthal velocity is Keplerian. Therefore, in this region the results are very similar to the inviscid results of CM81. In particular, there is an inward $u_{r}$ velocity induced by Cassen-Moosman drag as well as a contribution due to slow growth of stellar mass. When this $u_{r}$ is inserted into the mass conservation equation, one obtains a surface density in the inner disk that is similar to the CM81 solution except that CM81 set the homogeneous solution to the differential equation equal to zero, while St94 do not in order to match fluxes across the dense ring.

Simulations using a one-dimensional vertically integrated model of a protostellar disk with analytically prescribed infall (Ulrich, 1976; Cassen & Moosman, 1981) were performed targeting L1527. An implicit assumption of such models is that incoming flow quantities fully mix vertically with flow quantities in the disk at the radius of entry. Viscous (turbulent), magnetic, and gravitational torques were not included. One innovation compared to previous vertically integrated models is the inclusion of radial infall which is necessary for capturing the infalling rotating envelope (IRE) and a radial shock that separates the IRE from the main disk; radial infall is also a source of mass and momentum. Another difference is our use of all the unsteady transport equations without neglect of pressure. Most previous works assume Keplerian $u_{\phi}$ which eliminates many important effects such as pressure wave propagation, and reduces the problem to a single equation for surface density.

Except for a crucial disagreement, the simulation velocity and temperature profiles agree qualitatively with ALMA observations of L1527 (Ohashi et al., 2014; Sakai et al., 2014a, b, 2017; Aso et al., 2017). Specifically, the disk can be divided into three parts. (i) For $r>r_{\mathrm{shock}}$ (the IRE) there is radial infall with $u_{\phi}\propto r^{-1}$ which is angular momentum preserving. Since the mechanical energy very closely satisfies $E=0$ in this region, the flow is ballistic with parabolic particle trajectories. (ii) For $r_{\mathrm{c}}<r<r_{\mathrm{shock}}$ ($r_{\mathrm{c}}$ being the centrifugal radius), the radial velocity is greatly reduced by the radial shock but the angular velocity remains $\propto r^{-1}$. A separate non-LTE 1D shock calculation showed that the grain temperature rises to 60 K across the shock and remains near this value in a thin region. However, the pre-shock density was a factor of 16 smaller compared to the value necessary for SO desorption according to the study of Aota et al. (2015). Assuming that desorption nevertheless does take place gave a value of the SO column density close to that observed by Sakai et al. (2014a). Given the assumptions, this agreement is very likely, fortuitous and further studies should be performed. (iv) Inward of $r_{\mathrm{c}}$, the azimuthal velocity is close to Keplerian with an accretion mass flow, $\dot{M}(r)$, due to drag induced by sub-Keplerian vertical infall. $\dot{M}(r)$ is nearly uniform over a significant range of $r$ and this implies subtle deviations from Keplerian $u_{\phi}$.

The crucial difference is that the ALMA observations of Sakai et al. (2014a) are consistent with a ballistic maximum velocity relation (34) and imply that the radial shock is located at $r_{\mathrm{shock}}=r_{\mathrm{c}}/2$ and coincides with a ballistic centrifugal barrier, defined as the point of closest approach to the star of infalling parabolic orbits in the midplane. On the other hand, simulations give $r_{\mathrm{shock}}\approx 1.5r_{\mathrm{c}}$. This was the case even when we attempted to mimic the set-up of the Sakai et al. (2014a) ballistic argument by including only radial infall and using the ballistic flow in the initial condition.

Since the flow just upstream of the radial shock must be sufficiently supersonic to account for the post-shock temperature inferred from observations, the radial shock cannot coincide with the centrifugal barrier where $u_{r}=0$.

Nevertheless, since the Mach number of the ballistic flow can be large close to the centrifugal barrier (Figure 9a), it is possible that the radial shock is located close to the ballistic centrifugal barrier in reality. There might be a physical effect lacking in our simulation that would displace the shock from the region $r>r_{\mathrm{c}}$ to the region $r<r_{\mathrm{c}}$ where the ballistic Mach number is comparably large.

It was pointed out that that Model 2 of Ohashi et al. (2014) fit to their observations gives a jump in radial velocity at $r=120$ au, which is outward of their centrifugal radius at $r=54$ au. This is qualitatively in line with our simulations.

We noted (§5) that the radial shock in the magnetic core collapse simulations of Zhao et al. (2016) is also not in the vicinity of a centrifugal barrier and is in fact outward of $r_{\mathrm{c}}$. Therefore, if we are missing a shock-displacing physical effect, the same must be true of these simulations.

Some observers identify the shock location with a transition to Keplerian $u_{\phi}$. In the simulation, however, angular momentum is preserved across the radial shock at $r_{\mathrm{shock}}\approx 1.5r_{\mathrm{c}}$ and $u_{\phi}$ remains $\propto 1/r$ until $r=r_{\mathrm{c}}$ where the transition to Keplerian flow occurs. If the observers are correct, then there must be an angular momentum reduction process at the shock that the simulation does not include. In favor of the simulations, it may be stated that according to Aso et al. (2017), the transition to Keplerian flow occurs at 74 au while the radial shock is located at 100 au according to Sa14a.

Finally, the Toomre $Q$ parameter suggests gravitational instability in the outer disk where the radial velocity is small and therefore the Toomre analysis is valid.

1. 1.

   Is there a physical effect which neither we nor magnetic rotating-collapse simulations have included, that could move the radial shock to the left of the maximum ballistic Mach number in Figure 9a?

2. 2.

   The vertical structure of the disk should be elucidated, first via $(r,z)$ axisymmetric simulations. The vertically integrated model assumes instantaneous mixing of infall momentum with disk momentum. In really, this may occur by turbulence which may also lead to angular momentum transport and disk growth. As pointed out by Stahler et al. (1994) and shown in Figure 6, there is vertical shear of the planar velocity components $u_{r}$ and $u_{\phi}$. Axisymmetric and 3D simulations should be performed to study whether this shear leads to turbulence. Such shear is absent when there is purely radial infall. What is the structure of the flow and turbulence in this case?

3. 3.

   Magnetic core collapse simulations have shown that outflows impinge upon and shut-off infall at some locations on the disk surface. The physics of this process should be studied: for instance, what is the ram pressure of the outflow relative to the infall? Is the impingement turbulent and does it affect grain growth and transport? How does the lack of vertical infall together with an outflow affect disk structure and evolution?

4. 4.

   How valid is the UCM infall model when compared with magnetic collapse simulations. For instance, where vertical infall is present (i.e., not blocked by an outflow) is this infall sub-Keplerian as in the UCM model? Is the flow in the IRE described by zero energy parabolic orbits? A step toward such a comparison has been taken by (Lee et al., 2021, their §5.5.3) who compare the UCM mass source function with their simulations.

5. 5.

   The simulations suggested gravitational instability in the outer disk and where the radial velocity is small, which renders Toomre $Q$ criterion valid. This should be investigated even in the context of a planar $(r,\phi)$ treatment with self-gravity.