Steady Wind-Blown Cavities within Infalling Rotating Envelopes: Application to the Broad Velocity Component in Young Protostars

Wind-driven outflows are observed around a broad range of accreting objects throughout the Universe, ranging from forming low-mass stars to super-massive black holes. We study the interaction between a central isotropic wind and an infalling, rotating, envelope, determining the steady-state cavity shape formed at their interface under the assumption of weak mixing. The shape of the resulting wind-blown cavity is elongated and self-similar, with a physical size determined by the ratio between wind ram pressure and envelope thermal pressure. We compute the growth of a warm turbulent mixing-layer between the shocked wind and the deflected envelope, and calculate the resultant broad line profile, under the assumption of a linear (Couette-type) velocity profile across the layer. We then test our model against the warm broad velocity component observed in CO $J$=16--15 by Herschel/HIFI in the protostar Serpens-Main SMM1. Given independent observational constraints on the temperature and density of the dust envelope around SMM1, we find an excellent match to all its observed properties (line profile, momentum, temperature) and to the SMM1 outflow cavity width for a physically reasonable set of parameters: a ratio of wind to infall mass-flux $\simeq 4\%$, a wind speed $v_{\rm w} \simeq 30$ km/s, an interstellar abundance of CO and H$_2$, and a turbulent entrainment efficiency consistent with laboratory experiments. The inferred ratio of ejection to disk accretion rate, $\simeq 6-20\%$, is in agreement with current disk wind theories. Thus, the model provides a new framework to reconcile the modest outflow cavity widths in protostars with the large observed flow velocities. Being self-similar, it is applicable over a broader range of astrophysical contexts as well.


INTRODUCTION
Massive outflows are observed everywhere in the Universe, ranging from individual forming stars through galacticscale events. When supersonic stellar or galactic winds interact with the surrounding medium, be it the molecular envelope around forming stars or the intergalactic gas, they are observed to impart momentum and energy and entrain a slower-moving massive outflow. The actual entrainment mechanism and efficiency, however, remain poorly understood and highly debated, both because of a lack of strong observational constraints as well as a relative paucity of theoretical predictions against which to test observations.
A specific example of entrainment takes place when an accreting protostar launches a highly collimated jet, possibly surrounded by a wider-angle disk wind, carving out a large and slow massive outflow cavity into the parent cloud (Frank et al. 2014). When the Herschel Space Observatory (Pilbratt et al. 2010) started observing protostars in H 2 O and high-J CO rotational transitions, it quickly became clear that the dominant source of the emission was from molecular outflows (e.g., van Dishoeck et al. 2011;Kristensen et al. 2012Kristensen et al. , 2017. It also became clear that these emission lines highlight a different outflow component from the low-J CO transitions observed from the ground, such as J=2-1 and 3-2 (e.g., Yildiz et al. 2013). This Herschel -bright outflow component has both a significantly higher temperature 200 − 500 K and a larger line width at half maximum (FWHM) ≥ 30 km s −1 compared to low-J CO line profiles, where it only appears as a faint pedestal in very deep integrations (Margulis & Snell 1989). Accordingly, it was labeled the "broad" outflow component (Kristensen et al. 2012Mottram et al. 2014Mottram et al. , 2017. The physical origin of the "broad" warm outflow component, and its relation to both the slower cold outflow, seen in low-J emission, and the faster protostellar jet or wind is not clear. Two hypotheses have been put forward: either this broad component arises within a warm and dusty disk wind (Panoglou et al. 2012;Yvart et al. 2016) or it arises where ambient material is currently being entrained into the outflow by the protostellar wind, for example through non-dissociative shock waves (Kristensen et al. 2012Mottram et al. 2014Mottram et al. , 2017. While detailed, dynamical, and thermo-chemical predictions exist for the disk wind models, which reproduce the observed H 2 O emission (Yvart et al. 2016), only limited model predictions exist in the literature for the entrainment scenario, and it thus remains a hypothesis. The underlying physical issue is not a problem reserved for protostellar outflows, but remains an uncertainty for outflows in general.
A first type of entrainment scenario proposes that outflows are entrained by large jet bowshocks. These models predict substantial warm molecular material at intermediate velocities (e.g. Raga & Cabrit 1993;Downes & Cabrit 2003), but the resulting outflow cavities have been deemed too elongated compared to observations (Ostriker et al. 2001). To avoid this potential issue, a second type of entrainment scenario proposes that entrainment is dominated instead by a wide-angle wind. A particularly popular version of this "wind-driven" scenario assumes instantaneous full mixing between the shocked isotropic wind and the shocked envelope material. Due to the complete mixing, the cavity retains a primarily radial outflow motion with a roughly constant expansion speed over time (Li & Shu 1996;Lee et al. 2000), except very close to the disk mid-plane where the cavity remains trapped near the outer disk radius (Wilkin & Stahler 2003;Mendoza, Cantó, & Raga 2004;López-Vázquez, Cantó, & Lizano 2019). Such a radially expanding wind-blown cavity, however, grows too quickly: with expansion speeds ≥ 10 km/s similar to those observed in the broad component, it exceeds the typical radius ≤ 3000 au of protostellar outflow cavities (see e.g. Lee et al. 2015;Gueth, Guilloteau, &Bachiller 1998, for HH212 andL1157, respectively) in only a few 1000 yrs (Shang et al. 2006;López-Vázquez, Cantó, & Lizano 2019). This timescale is much shorter than the typical age accepted for Class 0 outflow sources (10 4 − 10 5 yrs, e.g., Kristensen & Dunham 2018, and references therein).
In order to circumvent this "age" problem, in this paper we consider wind-driven cavities with partial mixing, instead of full mixing, and explore stationary solutions for the cavity shape, formed as the fast stellar wind is obliquely deflected by the envelope and forced to flow along the cavity wall, instead of radially outwards. Indeed, numerical simulations of a spherical wind propagating into a rotating and infalling slab show that when mixing is inefficient, the cavity "flanks" quickly converge to a quasi-steady shape (Delamarter, Frank, & Hartmann 2000). A modest outflow cavity width can then be maintained over the whole duration of the Class 0 phase. Such a configuration is also prone to the development of a turbulent mixing-layer at the contact discontinuity between the shocked wind and the shocked envelope gas, as they slide past one other (Raga, Cabrit, & Cantó 1995). In this paper we will therefore follow a deliberate path to testing whether such a mixing-layer might explain the broad spectral component observed around protostars by Herschel /HIFI, and at the same time the observed cavity sizes.
Steady wind cavity solutions were first computed by Barral & Cantó (1981) for an isotropic wind expanding into a thick isothermal self-gravitating toroid. Smith (1986) showed that similarly elongated "flame-like" cavities could also be obtained for isotropic envelopes with a purely radial pressure profile p(r) ∝ r −n , provided that n < 2 and the wind is obliquely deflected at its closest point of impact (e.g. by a small-scale thin disk). Both of these early calculations show that the addition of a dense, disk-like, component along the horizontal axis can provide the required equatorial pinch to create steady, elongated outflow shapes similar to those observed around protostars. Similar physical conditions are expected for galactic-scale outflows (see e.g. Aalto et al. 2016). In this paper, we proceed one step further than these previous investigations by computing stationary solutions in a more realistic density and velocity distribution for the envelope, namely the more sophisticated Ulrich (1976) infalling and rotating solution. Despite an identical ambient density distribution, our cavity morphologies will strongly differ from the calculations of Wilkin & Stahler (2003), Mendoza, Cantó, & Raga (2004), and López-Vázquez, Cantó, & Lizano (2019) in that we assume weak mixing, instead of full-mixing, between the shocked wind and the shocked envelope, and we include the effect of thermal pressure in the envelope. These two ingredients allow the existence of stable stationary solutions on large scales, with pointed shapes at the pole.
In Section 2 we determine the stationary cavity shape formed by a wide-angle wind deflected by an infalling and rotating protostellar envelope. In Section 3 we consider the deflected wind material flowing along the cavity boundary, and the turbulent entrainment of envelope material within a mixing-layer. We then, in Section 4, compute the angular momentum and synthetic line profiles associated with material within the mixing-layer. Next, we quantitatively compare the model results against observations from Herschel of the broad CO component in the protostar Serpens-Main SMM1 (Section 5. Finally, in Section 6 we conclude with a recap of the implications of these results.

DETERMINATION OF THE CAVITY SHAPE
In this section we produced by an isothermal Ulrich (1976) infalling and rotating envelope model interacting with an isotropic wind. We determine the fundamental non-dimensional parameters and characteristic values and perform a stationary solution analysis in order to determine the range of cavity shapes produced. The shape of the thin shell formed by this interaction is determined by the (ram plus thermal) pressure balance between the wind and the envelope along with a "centrifugal term" due to the upward curving motion of the shocked wind layer.
The analysis in this section, as well as Section 3 and Section 4, is presented in dimensionless form in order to focus on the underlying self-similarity of the shape and the scaling relations underpinning the solutions, which have a general applicability for all manner of outflows. In Section 5, we adopt appropriate physical values in order to quantitatively determine the agreement between the model and observations, in the specific case of the warm CO outflow of a nearby protostar.

Infalling and Rotating Envelope Model
The Ulrich (1976) infalling envelope model generalizes from the case of an isothermal cloud by combining the spherical infall of envelope material under the force of gravity due to the mass of the central object M (Bondi 1952) together with a treatment of the centrifugal deflection of the flow due to initial solid body rotation. Thus, in the outer envelope, before angular momentum becomes dominant, the density structure is almost spherical with ρ(r) ∝ r −3/2 . At smaller radii, where centrifugal forces dominate, the flow streamlines are deflected toward the mid-plane, creating a disk-like structure inside of the fiducial radius Here, Γ ∞ is the specific angular momentum in the equatorial plane. Assuming a ballistic solution for the infalling material, the entire envelope solution can be described by a handful of parameters: the mass of the protostar, M , the size of the disk, r d , and the mass infall rate,Ṁ inf . Following Ulrich (1976) we therefore find for the density in the envelope where r is the spherical radius, θ is measured from the disk plane, and the subscript "0" denotes the initial value at very large distance from the origin. At any location (r, θ) in the infalling envelope, the initial angular origin of the streamline, θ 0 , can be obtained by solving Finally, the three components of the velocity of the infalling material at position (r, θ) are given by and Note, we have corrected the typographical error in equation 8 of Ulrich (1976) as noted by Tobin et al. (2012).

Wide-angle Wind
In this paper, we assume a spherical isotropic wide-angle wind of constant speed v w and mass-loss rateṀ w , with the density profile This approach enables a useful comparison with previous work (Barral & Cantó 1981;Smith 1986;Wilkin & Stahler 2003;Mendoza, Cantó, & Raga 2004;López-Vázquez, Cantó, & Lizano 2019) and an applicability to a wide range of astrophysical contexts. For example, while observations of young stars and Class 0 protostars show a strong and fast jet-like component along the axis, surrounded by a (seemingly) mostly empty lower-velocity outflow cavity, several MHD models predict that the jet may only be an "optical illusion" and may represent only the central densest core of a wider-angle wind, launched either from the inner disk edge ("X-wind" model, Shang, Shu, & Glassgold 1998) or from a larger portion of the disk surface ("D-wind" model, Cabrit, Ferreira, & Raga 1999).
Thus, in the context of protostellar outflows, the isotropic wide-angle wind provides an acceptable approximation to the X-wind in equatorial regions, where the interaction is the most critical to define the overall cavity shape (see below). We note that the addition of a strongly directed jet-like wind enhances breakout along the outflow axis and will thus modify the cavity shape in the polar regions. Comparison with observations should thus focus on regions close to the flow base, at wide angles to the flow axis.

Determining Fundamental Non-Dimensional Parameters and Characteristic Values
The trapping, breakout, and early evolution of the cavity formed by an isotropic wind colliding against an Ulrich (1976) infalling envelope was first calculated under the full-mixing hypothesis, including stellar gravity and various degrees of wind collimation (Wilkin & Stahler 2003), and envelope rotation (López-Vázquez, Cantó, & Lizano 2019). Full-mixing requires that the shell expansion is almost radial, hence the effect of thermal pressure in the envelope was neglected compared against the infall ram pressure in the frame of the expanding shell. A simplified study of the asymptotic shell expansion based on simple ram pressure balance, was conducted by Mendoza, Cantó, & Raga (2004). The authors show that it depends only on a single free parameter, namely the ratio of wind ram pressure to the fiducial infall ram pressure at r d , where is the Keplerian velocity at r d . Trapped solutions with sizes less than r d were found for λ < ∼ 1/2 (see Figures 5 and 8 from Mendoza, Cantó, & Raga 2004). For values of λ > 1/2, the cavity solutions were found to break out and expand forever, remaining pinched only along the disk mid-plane near r d . Similar results were found by Wilkin & Stahler (2003); López-Vázquez, Cantó, & Lizano (2019), with an additional breakout criterion on v w /v d for the polar cap to escape stellar gravity. For the X-wind and D-wind models currently favored in protostars, the denser and faster jet-like components along the axis will greatly facilitate breakout along the pole compared with the requirements for an isotropic wind (cf. discussion in Wilkin & Stahler 2003). Therefore, here we will assume that initial breakout has occurred and not consider this velocity constraint in our models.
In the present investigation, we are interested in finding steady asymptotic solutions to these breakout scenarios. For this, we assume instead that at most weak mixing occurs between the shocked wind and the shocked envelope material. Under this assumption, these two shocked and deflected layers will flow past each other, with an intermediate thin mixing-layer developing along the contact discontinuity between them.
A further difference between our models and those of Wilkin & Stahler (2003); Mendoza, Cantó, & Raga (2004); López-Vázquez, Cantó, & Lizano (2019) is that we take into account the role of thermal pressure in the envelope in confining the shell. This is the crucial element allowing to reach a steady configuration on large scales, instead of infinite expansion. Given that the density distribution in the envelope retains a modified, r −3/2 power-law, we anticipate that the resulting steady cavities will appear similar to the Smith (1986) elongated outflow cavities.
We thus introduce a characteristic scale length, r s , as the location where the ram pressure in the wind, ρ w v 2 w , is balanced by the thermal pressure in the equivalent spherically symmetric infalling envelope, ρ inf c 2 s . Solving this equality yields: where we define the useful characteristic velocity v 0 , fixed by the source ejection versus accretion physics, as Another region where thermal pressure will dominate infall ram pressure is near the equator, where the cavity is strongly pinched at r r d (see López-Vázquez, Cantó, & Lizano 2019) such that infall motions become almost parallel to the cavity walls. This introduces a second fundamental non-dimensional parameter in our model, Λ, which is proportional to the ratio of wind ram pressure to envelope thermal pressure at r d : where is the Mach number of the infall velocity at r d . We show in the next section that Λ only affects the cavity shape very close to the disk mid-plane, where it determines the initial foot-point and opening angle. While the disk centrifugal radius r d provides an appropriate scaling for the geometry at the flow base, we expect r s to provide an appropriate scaling for the geometry at large distances from the disk, where the cavity is confined by the thermal pressure in the envelope, rather than by the infall ram pressure.
Finally, we note that while trapped solutions confined by gravity on small scales, ≤ r d , were shown to be unstable, our steady shells on large scales, r d , are confined by thermal and ram envelope pressure and thus expected to be stable (see discussion in Wilkin & Stahler 2003). Moreover, our assumption of weak (instead of full) mixing allows a non-radial escape route for the bulk of the material reaching the shell, by either moving upward (wind) or downward (envelope). Such situations are much more stable against wind or infall variations than the full mixing case. This is supported by the robustness of the shell shape with respect to changes in initial or global parameters.

Calculating the Cavity Shape
To determine the location of the static boundary where the wind interacts obliquely with the infalling envelope, we follow the formalism of Matsuyama, Johnstone, & Hollenbach (2009) (their equations [2-5], which are derived in the appendix to that paper) which keeps track of both the mass and momentum flux deposited along the boundary by the shocked wind on the inner side and by the infalling material on the outer side.
The starting condition that ∂r/∂θ = 0 at the pole, used to compute cavity shapes with full mixing (eg. López-Vázquez, Cantó, & Lizano 2019), is no longer a requirement in the case of weak mixing, where flame-like shapes are allowed. Instead, we integrate from the disk mid-plane up, following Smith (1986). A starting location in the disk, which also explicitly sets the angle of incidence, is thus required in order to solve these equations. As noted above, López-Vázquez, Cantó, & Lizano (2019) find that breakout solutions lead to a strong pinch on the disk-plane near r r d . We have performed a detailed analysis of the possible angles of incidence allowed at the mid-plane as a function of r/r d < 1 for all combinations of λ and Λ (see Appendix A). We find that when the cavity is forced to meet the mid-plane at r/r d 1, the required angle of incidence at the mid-plane is such that the infalling streamlines approach the cavity wall from inside the cavity -an unphysical solution. As the mid-plane crossing approaches r d , however, there always exists a location where the angle of incidence of the cavity with the mid-plane is parallel to the infalling envelope streamlines. We therefore use this location as our foot-point for the cavity wall. Solutions close to this starting position quickly converge above the disk to the same surface; therefore, the exactness of the starting position is not critical for these models.
The left panel of Figure 1 plots the shape of the cavity scaled by r d , for a variety of values of Λ, while fixing λ = 1/2 (therefore Λ = M 2 d ). With this scaling, the envelope appears broadest and tallest when Λ is large due to the larger ratio of ram pressure in the wind to thermal pressure in the envelope at r d . As shown by the right panel of Figure  1, however, all the solutions with different Λ values are self-similar away from the base and actually have identical physical sizes in units of r s (defined in Equation 10). The height of the cavity is found to be Z max ∼ 2.5 r s , while the maximum cylindrical radius of the cavity is R max ∼ 0.16 r s . It is important to note that the self-similarity of these solutions breaks down near the base, where the relevant scaling length remains r d for all Λ.
The location and incidence angle at which the cavity intersects the mid-plane for all these solutions is set by requiring that the cavity interface be tangent to the infalling envelope streamlines. Figure 2 shows in detail the cavity shapes  near the disk surface, as well as the orientation of the infalling streamlines from the envelope for the cases investigated in Figure 1. From the figure, it is clear that the smaller Λ solutions intersect the base somewhat interior to the larger Λ solutions. These solutions therefore result in less interaction with the infalling envelope, as can be seen in the figure by noting the trajectories of the envelope streamlines.
Until now we have fixed λ = 1/2; however, breakout solutions should exist for other values of λ. In Figure 3 we show the cavity shape for Λ = 25 and a variety of λ. At the base, the Λ solutions are independent of λ while the maximum cylindrical radius and height increase with increasing λ, quickly asymptoting to a fixed solution. Furthermore, in Figure 4 we plot, as a function of λ, both the maximum cylindrical radius R max of the cavity (in units of r s ) and the ratio of the height to the cylindrical radius at this widest point in the cavity, for each of the Λ cases used in the previous figures. It is clear from these plots that for λ > 1/2 all the solutions become remarkably self-similar, with only a slight hint that the shapes are slightly broader for larger λ and smaller Λ. We further note that steady breakout solutions are found even for values of λ < 1/2 (down to λ 0.2). These solutions differ from the Mendoza, Cantó, Cavity profiles for fixed Λ = 25 and varying λ showing breakout solutions when λ > 0.2. The dashed lines represent streamlines for material flowing within the infalling envelope onto the disk of radius r d = 2 r s /Λ 2 . As discussed in the text, in units of r s the cavity shape becomes fixed for any Λ at large λ.  Ulrich (1976) infalling envelope for a range of λ and Λ values. We find steady breakout solutions even when λ < 1/2, to the left of the dotted vertical lines. Right: Ratio of the height of the cavity to the cylindrical radius of the cavity, measured at the location where the cavity is widest.
& Raga (2004) trapped solutions which are required a priori to have ∂r/∂θ = 0 along the vertical axis, creating a roundish "cap" strongly confined by infall ram pressure when λ < 1/2. Instead, our 0.2 λ < 1/2 cavities maintain flame-like shapes, such that infall ram pressure at the tip is strongly reduced by the highly oblique incidence there, and allows for breakout.
The fact that our numerical solutions are not highly dependent on the initial location of the interface at the disk surface (see Appendix A), nor on the exact values of λ or Λ, confirms that they are stable equilibrium solutions, as expected when the confinement is dominated by envelope pressure.

FLOWS ALONG THE CAVITY WALL
The shape of the cavity wall as a function of the two defining input parameters (λ, Λ) was shown in Section 2 to be close to self-similar as long as λ > 1/2, especially at large distance from the intersection with the mid-plane. Thus it is reasonable to expect that the flow of deflected material from either the wind or envelope side of the cavity wall can also be described in terms of a single simplified parametrization, with small deviations as a function of Λ due primarily to the slightly varying physical situation near the disk surface.
In order to keep track of the various flows, we divide the cavity wall surface into three components. We present a schematic of the various regions in Figure 5. First, there is the deflected shocked wind (denoted in the text by a subscript 1) that travels upward, parallel to and on the inside of the cavity surface. Second, there is the deflected infalling envelope (denoted in the text by a subscript 2) which travels downward, parallel to and on the outside of the cavity surface. In the absence of mixing across the surface, these two flows remain independent and can be fully described at each location by a mass and momentum flux explicitly determined by integration along the surface (see equations [2][3][4][5] in Matsuyama, Johnstone, & Hollenbach 2009). Third, a turbulent mixing-layer (denoted in the text by a subscript L) in which slow moving deflected envelope material is entrained upwards by the fast deflected wind, can develop at the contact discontinuity between these two flows.

Solutions Without a Mixing-Layer
We begin with solutions in which there is no mixing between the upward and downward deflected flows. Figure 6 plots the mass flux as a function of (scaled) height for both the interior (Ṁ 1 ; upward) and exterior (Ṁ 2 ; downward) layers. The results are shown for a single hemisphere, and scaled accordingly, as there is no explicit requirement within the model for symmetry about the disk plane. The right panel fixes λ = 1/2 and shows solutions for four values of Λ whereas the left panel fixes Λ and shows solutions for four values of λ. In these plots it is important to recognize that the mass-flux axes are scaled separately, with the out-flowing materialṀ 1 scaled to the total mass flux from the wind . .
. Figure 5. Schematic showing the various layers along the cavity wall, in which deflected wind material moves upward and deflected envelope material moves downward. A central turbulent mixing-layer with a linear "Couette-type" velocity gradient may grow between these two layers (c.f. Section 3.2). Mathematical notations used in the text for the velocity, mass-flux, sound speed, and density in each part of the flow are also shown for easy reference. over a hemisphere and the infalling materialṀ 2 scaled to the mass infall rate from the envelope over a hemisphere. Thus, while all the solutions are self-similar well above the mid-plane the relative importance of the upward versus downward mass flux depends on both λ and v w /v d (see Eqn. 8). In all cases the upward flowing surface asymptotes to the entire mass flux in the wind, as required, whereas the the cavity intercepts only a fraction of the infalling envelope, missing that part which lands on the disk surface between the cavity foot-point and r d . The trend with Λ, seen in the right panel, can therefore be understood as a direct consequence of the fact that smaller Λ solutions intercept the disk closer to the central source (see Figure 2). As shown in the left hand panel, all solutions at fixed Λ are approximately self-similar modulo the wind and infall mass flux scaling. A similar set of solutions is found for the momentum flux as a function of height along the interior (Π 1 =Ṁ 1 v 1 ; outward) and exterior (Π 2 =Ṁ 2 v 2 ; inward) cavity layers, as shown in Figure 7. The quantities shown in the figure are scaled independently to the fiducial momentum flux in the wind and infalling material. The relative scaling between these quantities, however, is explicitly λ (see Eqn. 8) and thus it is apparent that the downward momentum flux in the outer layer is always much less than the upward momentum flux in the inner layer except extremely close to the base. The left hand panel, again, shows that the solutions are approximately self-similar over a wide range of λ for fixed Λ.
Finally, the mean velocities (v 1 , v 2 ) for the two deflected flows along the surface can be calculated directly from the ratio of the respective momenta and mass fluxes. These are plotted in Figure 8. The quantities shown in the figure are scaled independently to v w and v d . The wind typically intersects the surface at an acute angle and thus the majority of the momentum from the wind is deposited in the deflected flow rather than contributing ram pressure to support the surface against the infalling envelope. Therefore, the magnitude of the velocity of the deflected wind along the  the shocked layers ( Figure 9). We note that on large scales the pressure may be well approximated by Next, we compute the ratio of the thickness of the layers H 1 and H 2 against the radial extent of the cavity, as a function of location along the cavity. In detail, the thickness of each deflected layer is equal to the surface density divided by the mass density within the layer. In the thin shell approximation, the density within each layer is set by pressure equilibrium through P = ρ 1 c 2 1 = ρ 2 c 2 2 , where c 1,2 are respectively the sound speed in the deflected wind and the deflected infalling envelope layer. Since P is proportional to the ram pressure from the wind, the scalings for the relative thickness of the two deflected layers simplify to and Figure 10 presents the relative thickness H/R of the two layers, normalized by their respective relevant scaling in each case. Except near the top of the cavity, where the solution converges back toward the axis of rotation, the deflected wind and envelope layers remain thin provided v w /c 1 > 10 and Λ (c s /c 2 ) 2 > 10, respectively. Furthermore, Figure 8 shows that the deflected wind suffers oblique shocks with speed v s = v w,⊥ 0.2 v w . Using the general expression for the maximum temperature reached behind a hydrodynamical shock of speed v s to set an upper limit on c 1 , it may then be determined that the condition for the shocked wind layer to remain thin is equivalent to v w /c w > 10, with c w the isothermal sound speed in the wind.
Similarly, Figure 8 also shows that the deflected envelope undergoes only small velocity jumps v s = v inf,⊥ < 0.1v d . Since we will always have c 2 ≥ c s , a conservative condition ensuring that the deflected envelope layer will remain thin, regardless of the value of v d , is simply that Λ > 10. This conservative condition on Λ becomes unnecessary, however, if v d is large enough for the velocity jumps to remain supersonic everywhere along the cavity wall. Using the expression for the maximum temperature behind a shock at v s to set an upper limit on c 2 , and rewriting the condition for a thin envelope layer as λ > 5 (c 2 /v d ) 2 (using Eqn. 12), we find that the inequality is then automatically fulfilled under our wind breakout condition λ > 0.2, regardless of the value of Λ. The combined results presented in Figures 6 -10 reveal that, in the absence of mixing, the deflected wind mass flowing upward along the boundary surface will be similar to the total mass flowing in the wind,Ṁ w , and that the magnitude of the momentum in this deflected flow will also be close to the wind momentum fluxṀ w v w .
Due to the large aspect ratio of the cavity, the bulk of the deflected wind flows roughly perpendicular to the disk and at a high velocity, v 1 ∼ v w . Alternatively, if mixing takes place between the momentum-rich outward flowing layer and the mass-rich infalling layer, the internal velocity structure of this turbulent mixing-layer should be significantly differentiated, as for example through a linear velocity gradient such as occurs in a Couette flow (e.g. Raga, Cabrit, & Cantó 1995). Such an occurrence will naturally produce a wider spread of velocities between v 2 and v 1 .

Solutions with a Mixing-Layer
A general formalism for the growth of a turbulent mixing-layer between two axisymmetric flows was derived by Raga, Cabrit, & Cantó (1995). Their main formulae included some ambiguities and typographical errors, and are therefore reproduced in corrected form in Appendix B. These authors assume that within the mixing-layer there is both a fixed temperature, referred to by its sound speed c L , and a fixed pressure P across the layer, varying only as a function of position along the flow. Across the layer they further assume that the velocity varies linearly, as in a Couette flow, bounded by the velocities of the fast v 1 layer and slow v 2 layer (see Figure 5). The coupled equations for the change in mass flux and momentum flux within the mixing-layer due to entrainment across the inner and outer boundaries are then solved so as to determine the entrainment required across each bounding surface in order to maintain the imposed Couette conditions.
As detailed by Raga, Cabrit, & Cantó (1995), entrainment occurs in two ways (see their equations 1 and 2): through the geometrical growth of the mixing-layer, intercepting a fraction of the flows on either side, and through "turbulent entrainment" on the slow-moving side (here the deflected envelope) as it is dragged into the mixing-layer by the fastmoving side. In this paper, we follow the Raga, Cabrit, & Cantó (1995) prescription for the turbulent entrainment velocity v ent = αc 2 2 /c L (see their equation 3), with a constant turbulent mixing parameter α. For simplicity, we will further assume constant values of c L , c 1 and c 2 at all positions.
At any location x along the surface, the linear velocity gradient across the turbulent layer assures that the massweighted mean velocity within the mixing-layer is v L (x) = [v 1 (x) + v 2 (x)] /2. At the same time, the ratio of the momentum fluxΠ L versus the mass fluxṀ L in the mixing-layer is skewed toward the higher velocities within the layer, such that To maintain this ratio as slow envelope materialṀ L2 is turbulently entrained across the outer boundary surface, a significant amount of fast-flowing shocked wind materialṀ L1 must also be entrained across the inner boundary, with the exact proportion set by the (changing) physical conditions along the surface. For the trivial case of a constant velocity, v 1 , in the fast moving layer and no motion on the slow-moving side (v 2 = 0), Eqn. 17 shows that the Couette flow requiresΠ L (x) = (2v 1 /3)Ṁ L (x). Since zero momentum can be provided from the stationary side, this condition is met by the mass entrainment rate from the fast moving wind side being exactly twice the mass entrainment rate from the stationary side (Cantó & Raga 1991).
A careful consideration of the interface solutions presented in Section 2 shows that for the self-similar cavities in this paper v 1 ∼ v w and |v 2 | v w , in all examined cases. This simplifies the general formulae described by Raga, Cabrit, & Cantó (1995); however, the changing radius of curvature and the steadily dropping pressure across the calculated wind-envelope surface conspire such that the detailed solution for entrainment must be calculated numerically and separately for all parameter pairs (λ, Λ). Fortunately, despite this somewhat more complicated geometry, Figure 11 shows that for all pairs (λ, Λ) the mass-flux entering the mixing-layer from the fast-flowing wind side,Ṁ L1 , remains close to twice the turbulent entrainment from the envelope side,Ṁ L2 .
The self-similarity of the boundary location also provides, for each (λ, Λ) pair, a scaling relation for the efficiency of the turbulent mixing solutions in terms of the physical parametersṀ w , v w , c L , as well as α. The turbulent entrainment into the mixing-layer from the slow-moving, envelope, side of the boundary,Ṁ L2 , can be found by integrating the turbulent entrainment along the boundary surface, x. That is: where ρ 2 (x) is the density of the deflected shocked envelope along the outer boundary and v ent = α(c 2 2 /c L ) is the parametrized entrainment velocity. This equation is exact for situations where the shocked ambient medium is static and remains an excellent estimate when |v 2 | v w . The result may be rewritten in terms of the pressure across the boundary surfaceṀ which further reduces by recognizing that the pressure at location x along the surface is set explicitly by the ram pressure: P (x) = a(x)Ṁ w v w where a(x) ∝ sin 2 γ/r 2 takes into account the varying angle of incidence between the isotropic wind and the boundary surface. Thus, Furthermore, given the almost fixed ratio between mass entrained into the mixing-layer from the deflected wind versus the deflected envelope (see Figure 11), the scaling forṀ L1 , andṀ L =Ṁ L1 +Ṁ L2 , should be the same as that foṙ M L2 . Utilizing this scaling as normalization, Figures 12 and 13 show the fraction of deflected wind and deflected envelope that is entrained into the mixing-layer as a function of height above the mid-plane, for a variety of (λ, Λ) pairs. Selfconsistent solutions require that these fractions remains less than unity, otherwise the reservoir of shocked material flowing along the cavity walls is not large enough to feed material into the mixing-layer at our assumed rates. From Figure 12, it is clear that for the wind side this constraint requires Similarly, for the envelope side ( Figure 13) the constraint is trivially met for the same physical parameters assuminġ M w < 0.5Ṁ inf , except at extreme heights, Z > r s , where the cavity shape converges to the axis of rotation. Furthermore, combining the information in Figures 12 and 13, and using the results of Figure 6, Figure 14 shows that the total mass flux in the mixing-layer for one outflow cavity lobe is almost independent of the (λ, Λ) pair. Within a factor of a few, the asymptotic value at high altitudes iṡ where α max is defined in Eqn. 21. We note that if α > α max , our physical model will not entirely break down. The mixing-layer will simply grow until it eventually engulfs all of the deflected wind layer andṀ L saturates at its maximum possible value ofṀ w /2. Without a fast laminar wind layer to enforce a Couette flow, however, the velocity field in the mixing-layer would become a gaussian velocity distribution peaked around a mean value v w /2, leading to line profiles much narrower than in a linear gradient Couette flow. and where Ω 2,φ = v 2,φ /(2πR) and Ω L,φ = v L,φ /(2πR). By integrating Eqn. 23 over R along the interface from the top of the cavity to the disk plane, we obtain Ω 2,φ , and hence v 2,φ . Subsequently by inputting the derived Ω 2,φ into Eqn. 24 and integrating it over R from the disk plane upward along the interface, we obtain Ω L,φ , and hence v L,φ . In the top panels of Figure 15, we show v L,φ (solid lines), v 2,φ , and v inf,φ as a function of R for fixed Λ and varying λ (left panels), as well as for fixed λ and varying Λ (right panels). We can see that for Λ > 20, v L,φ is no larger than a few tenths of v d and therefore can be considered as negligible compared to the bulk velocity of the gas in the mixing-layer (v L ≈ 0.5 v w ) when determining the shell shape and computing the observed line profiles.
The bottom panels of Figure 15 plot the specific angular momentum Rv φ for the same three velocity components as in the upper panels, as a function of Z. We can see that the specific angular momentum in the mixing-layer is virtually independent of λ and Λ. Since twice as much material is entrained from the (non-rotating) wind side than from the envelope side, the initial value of the specific angular momentum at the base is one third of that in the deflected envelope. As gas is advected upwards in the mixing-layer, this rotating material gets mixed with deflected ambient material of smaller specific angular momentum. However, since most of the mass entrainment occurs at Z ≤ 0.1 r s , the specific angular momentum in the mixing-layer remains close to its initial value, 0.15 v d r d .

Mixing-Layer Line Profile
Having developed a model for how material is entrained, the next step is to calculate the resulting line profile for direct comparison to observations. The model calculations in the preceding sections were all dimensionless and are thus applicable to essentially any type of entrainment irrespective of physical scale. Thus, with the appropriate scaling this model could be compared with outflows from protostars (e.g., Mottram et al. 2014;Kristensen et al. 2017) to extra-galactic outflows (e.g. Aalto et al. 2016Aalto et al. , 2017. The line profiles for the shocked wind layer and for the turbulent mixing-layer are generated independently under the optically thin assumption by computing the flux dF e (v LOS ) = e dV emitted by each elementary volume element dV , where e is the emissivity per unit volume of the molecular line of interest. At every height, z, along the surface, the  flux from each azimuthal interval, dφ, is added to the velocity bin corresponding to the line-of-sight velocity v LOS of the volume element. For the shocked wind layer, the velocity has a unique modulus v 1 , dependent only on z, whereas for the mixing-layer, the flux is evenly distributed in velocity between 0 and v 1 prior to projection. In Figure 16 we show example line profiles for our reference model presented in the previous section (λ = 1/2, Λ = 25). We scale the projected velocities by v w and integrate the emission up to 200 r d (0.6 r s ) from the base of the outflow. Four viewing angles are provided. The black curves show the line profiles for the shocked wind layer only, assuming emissivity proportional to density -mimicking the high density LTE regime. Except for the edge-on case, the emission always peaks near the projected wind velocity v w sin θ obs . This occurs because the bulk of the deflected wind flows at v 1 ∼ v w and is roughly parallel to the disk axis, due to the elongated shape of the cavity. The blue curves, on the other hand, show the computed line profiles for the mixing-layer only, assuming again that the emissivity is proportional to the density. As expected, the line profiles are much broader and flatter, peaking at zero and extending to a fraction of v w . Finally, the solid red curves show predicted line profiles from the mixing-layer when emissivity is proportional to the density squared, mimicking the low density limit. This emissivity condition increases the contribution of dense regions near the base, where the cavity opening angle is still large. Due to projection effects, it produces more extended line wings close to edge-on, and enhanced low-velocity emission when close to pole-on.
Our described model of partial entrainment of the wind, along with some of the exterior envelope, through a turbulent mixing-layer thus ensures that a fraction of the outflowing material is moving slowly due to the linear velocity gradient, Couette-type flow within the mixing-layer. Turbulent dissipation within the mixing-layer also provides a heating mechanism to make this material warmer than the shocked wind layer, hence brighter in high-J CO lines.
In the following section, we investigate whether such a model could explain at the same time the line profile, momentum, and temperature of the broad component observed in CO by Herschel towards the Serpens-Main SMM1 protostar, as well as the observed outflow cavity size, for reasonable envelope and wind parameters. Line profiles for the model with Λ = 25 and λ = 1/2. Each panel shows the result for a different direction of viewing as labeled, where θ obs is measured from the plane of the disk. Blue and red curves in each panel indicate the line profiles of material in the mixing-layer, under the assumption that emissivity of CO J = 16 -15 ( CO ) is proportional to the density and square of the density, respectively. For reference, we also show in each panel the emission profile of material in the shocked wind layer alone if no entrainment occurs, under the assumption that emissivity scales linearly to gas density (black curve).

APPLICATION TO THE BROAD COMPONENT OF CLASS 0 PROTOSTAR SERPENS-MAIN SMM1
As an illustration of the applicability of our formalism to a specific outflow case, we compare our model predictions with observations of the protostar Serpens-Main SMM1. This protostar is located in the Serpens Main cloud at a distance of 438 pc (Herczeg et al. 2019). The protostar has a bolometric luminosity of ∼ 100 L (Goicoechea et al. 2012), and is therefore on the border between low-and intermediate protostars. The envelope is correspondingly massive, with an estimated mass of ∼ 50 M (Enoch et al. 2008;Kristensen et al. 2012). When observed at high angular resolution, the protostar breaks up into multiple sources; however, the central most massive protostar is responsible for the primary outflow (Hull et al. 2016;Hull et al,. 2017;Le Gouellec et al. 2019). When observed in H 2 O and high-J CO emission with the HIFI instrument on Herschel, this source shows the brightest line intensity in H 2 O and CO J=16-15 in the sample of Kristensen et al. (2017). For this reason, the broad line component of SMM1 also has the highest signal-to-noise ratio. Hence, this protostar is a natural choice for a first comparison between the model presented in this paper and observational data.

Broad component line profile and wind velocity
When Herschel -HIFI started observing H 2 O emission toward protostars, one of the biggest surprises was that the velocity-resolved line profiles typically were dominated by a broad outflow component with a FWHM of 30 km s −1 (e.g., Kristensen et al. 2012Kristensen et al. , 2017Mottram et al. 2014). This line width was significantly larger than seen in low-J CO from the ground, e.g., J=3-2, where the FWHM is 15 km s −1 (Kristensen et al. 2012). It also became clear that when observing higher-J CO transitions with HIFI, the line profiles started resembling the H 2 O profiles so much so that the CO J=16-15 profiles are indistinguishable from the H 2 O profiles . That the CO profiles vary with excitation suggests that the change in shape is indeed due to excitation as opposed to chemistry. Furthermore, the change in profile shape is likely related to an increase in temperature, because when calculating the rotational temperature from the ratio between CO lines, the temperature increases from 100 K to ∼ 300 K (Yildiz et al. 2013;Kristensen et al. 2017). Thus, the higher-J CO lines, and by implication the similar H 2 O lines, trace a warmer, faster-moving, component of the protostellar outflow as compared to what is seen in low-J transitions , and this component is primarily seen as a broad outflow component in the velocity-resolved line profiles. Kristensen et al. (2012) and Mottram et al. (2014) speculated that, because of the higher temperature and velocity, this broad component is tracing gas closer to a shock front, possibly located where the protostellar wind shocks against the infalling envelope in an irradiated C-type shock. The colder gas, traced by lower-J CO lines, then would be the subsequently entrained swept-up ambient gas. Alternatively, the heating and entrainment process of the broad component could take place within a turbulent mixing-layer at the interface between the shocked wind and infalling envelope. This alternate scenario is investigated below, using the model results presented in the previous sections.
For the source SMM1, Kristensen et al. (2017) found that the CO J=16-15 line profile could be decomposed into three Gaussian components, one broad (FWHM ∼ 20 km s −1 ) and two narrower components (FWHM ∼ 8-10 km s −1 ). The narrower components are only seen in this high-J CO line and likely originate in shocks very close to the protostar , and they are not considered further here. Figure 17 compares the broad component, extracted from the CO J=16-15 line profile in SMM1 by Kristensen et al. (2017) after removal of the two narrower component Gaussian fits, with our model predictions. Excellent agreement is found for a mixing-layer with v w = 25 − 30 km/s, an inclination of θ obs = 30 • measured from the plane of the disk, and an emissivity proportional to density.
A few checks are in order to ensure that the model remains self-consistent. First, for the "thin shell" approximation to be satisfied, the wind velocity must remain large enough to produce a significant ram pressure which in turn confines the flow of shocked wind along the cavity surface. In Section 3.1, we showed that confinement requires v w /c w > 10 where c w is the isothermal sound speed in the wind. With our estimate for v w 30 km s −1 , we therefore require c w < 3 km s −1 , or a maximum wind temperature of 1000 µ K, with µ the mean molecular weight per particle (in a.m.u). This condition holds both in D-winds heated by ambipolar diffusion in Class 0 sources (Panoglou et al. 2012;Yvart et al. 2016) and in X-winds in the absence of mechanical heating (Shang et al. 2002). In Section 3.1, we also found that the deflected envelope layer will always remain thin when Λ > 10, regardless of the value of v d . We will verify that this condition holds in SMM1 in Section 5.5. Second, the model line profiles in blue that reproduce the observed line profile shape for SMM1 in Figure 17 are obtained only if a Couette linear velocity gradient exists across the mixing-layer. For this gradient to be maintained, the mass flow entrained in the mixing-layer from the wind side,Ṁ L1 , should not exhaust the available flux of shocked wind materialṀ 1 flowing along the shell. As discussed in Section 3.2, this constraint sets an upper limit on the turbulent entrainment coefficient in SMM1 α ≤ α max = (c L /v w ) 0.035 (see Eqn. 21 and Figure 12) where we have used our estimate of v w = 25-30 km s −1 from line profile modeling, and c L = 1 km −1 from the temperature 250 K inferred by multi-line CO analysis of the broad component in SMM1 (see Kristensen et al. 2017, and next section). We will verify below that this upper limit on α is still compatible with the observed momentum in the broad component of SMM1.

Outflow Cavity Size
In Figure 18 we compare a published CO outflow map of SMM1 (Hull et al. 2016) with our predicted self-similar cavity shape from Figure 1, for various values of the scaling parameter r s . Although only the inner region of the SMM1 outflow has been mapped at high angular resolution, the joint constraints on small and large scales indicate that r s must lie in the range 10, 000 − 40, 000 au. Hence, we adopt r s = 20,000 au as our fiducial value in the following.
The cavity physical scale r s requires a specific ratio of mass loss rate in the wind to infall rate in the envelope (see Since SMM1 is quite bright (100 L ) we adopt a fiducial sound speed in the envelope of c s 0.4 km s −1 corresponding to a temperature of 40 K. This value is consistent with radiative-transfer modeling of the dust emission from the envelope surrounding SMM1 (Kristensen et al. 2012), which recover a sound speed 0.3 -0.5 km s −1 in the envelope at the relevant physical scales, 350 -3000 au, shown in Fig. 18, which are also those of the Herschel/HIFI beam. With this value of c s , we find that the observed size of the outflow cavity in SMM1 can be reproduced with a wide-angle wind mass-flux on the order of 4% of the envelope infall rate, which is quite modest. In the following sections, we use our model and the observed momentum in the broad component to constrain the absolute value ofṀ w and then that ofṀ inf , through Eqn. 25.

Momentum in the mixing-layer
The two-sided momentum in the broad component of SMM1 was estimated from observations of CO J=3-2, 6-5, 10-9, and 16-15 taken with the JCMT, APEX, and Herschel -HIFI (Yildiz et al. 2013). The respective line profiles were rebinned to the same velocity scale and to channels of 3 km s −1 width. For each channel, a rotational diagram was constructed and, assuming LTE and optically thin emission, the rotational temperature and CO column density N CO were calculated. The rotational temperature was ∼ 250 K, irrespective of velocity. With this mass spectrum in place, the mass-weighted momentum Π BC of the broad component summed over both lobes is given by High-velocity CO jet and ionized outflow cavity . Low-velocity red-and blueshifted CO(J = 2 ! 1) from the ALMA data (red and blue contours, respectively), overlaid on the VLA 4 cm continuum image (grayscale). As in Figure  1, the image is centered on SMM1-a. The CO velocity ranges are 2 to 15 km s 1 (redshifted) and -20 to -5 km s 1 (blueshifted) relative to the v LSR of SMM1 of approximately 9 km s 1 . The contours are plotted at levels of 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95 ⇥ the peaks of the redshifted (3.76 Jy bm 1 km s 1 ) and blueshifted (4.16 Jy bm 1 km s 1 ) moment 0 maps. toward the extreme SE of the image, just beyond the SE radio knot (about 12 00 SE of SMM1-a; see Figure  1, left panel). These two CO features are not perfectly collinear; furthermore, the four main radio knots are not all collinear. Rodríguez-Kamenetzky et al. (2016) suggest that the precession of the jet could cause this noncollinearity; they attribute the precession to an unresolved 3 AU binary within SMM1-a and calculate a precession period of 20-30 yr and a precession angle of 10 (i.e., the angle between the central jet axis and the line of maximal could the c address this the emission 1.3 mm data Choi (2009 the 1.3 mm dust emissi 3.2-signifi (⇠ 3.5-4.0) Even in th dominated between th 3).   Figure 1). They have been rescaled to r s = 10, 000 au (magenta) and r s = 40, 000 (blue) in order to match the observed CO map profile.
where X CO is the (unknown) fractional abundance of CO molecules by number with respect to H nuclei in the broad component, and Π obs is the fiducial "observed" momentum assuming a standard interstellar CO abundance of 5 × 10 −5 . Equation 26 shows that the value of Π obs only depends on the observed CO intensity and excitation temperature, irrespective of the true X CO . It is therefore the quantity usually reported in observational papers. Using an updated distance, d = 438 pc, to the Serpens Main cloud (Herczeg et al. 2019), we recalculate the Yildiz et al. (2013) derived fiducial momentum inside a beam radius of R b = 5.5 = 2400 au to be Π obs 8 ± 2 × 10 −2 M km s −1 .
In our entrainment model, the two-sided momentum contained in the mixing-layer up to a distance z = ±R b is given by where η(z) is the normalized ratio plotted in Figure 14 and α max is defined in Eqn. 21. Since η(z) increases very slowly with height, the integral on z may be approximated as If the momentum in the broad component of SMM1 inside R b = 2400 au is provided by mixing-layer entrainment from a wide-angle wind, then 2Π L = Π BC . Using Equation 26 with Π obs 8 ± 2 × 10 −2 M km s −1 , and taking η b 1 (which we will verify in Section 5.5) we infer that the wide-angle wind must have a mass fluẋ

Infall rate
We can now compute the required envelope infall rate in our model, and compare with the infall rate independently suggested by dust envelope models. Combining Equations 25 and 28, we infer the required infall mass flux to reproduce both the outflow cavity size and the momentum in the broad component aṡ The required value is larger than typical infall rates for low-mass Class 0 protostars. It is in line, however, with that expected for sources with particularly massive envelopes such as SMM1, whose luminosity ∼ 100 L places it on the border between low-and intermediate-mass protostars. Using the estimated H 2 density at 1000 au, n 1000 , in the SMM1 dust envelope model of Kristensen et al. (2012) and rescaling by d 2 from d = 230 pc to 438 pc, we infer an "observed" envelope infall rate at R 1000 = 1000 au oḟ where the factor 1.4 accounts for the mass in the form of Helium. We note that M appears at the same power iṅ M inf andṀ env , hence its exact value, currently unknown in SMM1, does not matter for the comparison. There is therefore good agreement with our mixing-layer model as long as X CO , the CO abundance in the mixing-layer with respect to H nuclei, is close to the standard interstellar value of 5 × 10 −5 , and the turbulent entrainment parameter α is close to the maximum value to maintain a Couette flow, α max = c L /v w 0.03. We note that such a value of α matches very well with a model fit to supersonic mixing-layer experiments 1 , thus it appears physically plausible. With these values for X CO and α, the wind mass-flux that is required to provide the broad component momentum iṡ M w 6 × 10 −6 M yr −1 (see Eqn. 28).

Constraints on ejection / accretion ratio and disk radius
We have shown in Section 5.2 that the outflow cavity size in SMM1 can be reproduced with a modest ratio of wind mass flux to envelope infall rate of 4%. The disk accretion rate onto the central star, however, may be smaller than the envelope infall rate onto the disk. Assuming that the bolometric luminosity 100L of the SMM1 source (Goicoechea et al. 2012) is dominated by the accretion luminosity L acc GM Ṁ acc /R , and adopting stellar radii R on the birthline computed by Hosokawa & Homukai (2009), we infer a disk accretion rate ofṀ acc 10 −4 M yr −1 if M = 0.2 M , andṀ acc 3 × 10 −5 M yr −1 if M = 0.5 M . With the wide-angle wind mass-fluẋ M w 6 × 10 −6 M yr −1 derived in the previous section, the ratio of wind ejection rate to disk accretion rate is thus 0.06 − 0.2 for M = 0.2 − 0.5M . Such values are in the typical range predicted by D-wind and X-wind ejection models from accretion disks around young stars. Therefore, the wind mass-flux requirements in our model for SMM1 appear physically reasonable.
We next estimate the expected range of the parameter λ, the ratio of wind to infall ram pressure, for SMM1. From Eqn. 8 and Eqn. 10 we have 23 This is consistent with the condition λ 0.2 for which our cavity solutions break out and reach their full extent (see Figure 4), for the typical disk sizes in Class 0 sources . Interferometric continuum observations suggest that SMM1 possesses a particularly large and massive disk of 300 au (Enoch et al. 2009 This low value of Λ is consistent with our assumption of η b 1 on the scale z = R b 0.1r s of the Herschel beam (cf. the curves for η(z) in the right panel of Figure 14). Even with a large disk, r d 300 au, the condition Λ > 10 for the weakly shocked deflected envelope layer to remain thin is also fulfilled.

Temperature and Density in the mixing-layer
As an additional test of our model, we investigate whether the observed temperature, T L 250 K of the broad component in SMM1 suggested by multi-line CO analysis ) would be consistent with heating of the mixing-layer in our model by turbulent viscosity.
In principle, a full non-equilibrium thermo-chemical calculation should be performed as a function of position along the mixing-layer. Such a complex problem is, however, outside the scope of the present paper and is deferred to future work. For simplicity, we assume here that the temperature and chemistry in the mixing-layer have reached a steady-state on the scales observed by the Herschel HIFI beam, and check whether thermal equilibrium at T L 250 K could indeed be sustained.
Following Binette et al. (1999), we take a turbulent viscosity µ = (α/4)ρ L c L h with h the total thickness of the mixing-layer. A derivation of this expression for supersonic isothermal mixing-layers with a linear velocity profile (Couette flow) is given in Appendix C. We can then express the turbulent heating rate per unit volume in the mixinglayer as where we make use of (dv/dh) = v w /h with h =Ṁ L /(2πRρ L v L ), R is the local shell radius, P = ρ L c 2 L is the local pressure in the layer, and α max = c L /v w is defined in Eqn. 21.
Thermal equilibrium at constant T L will be maintained as long as where Λ exp is the rate of "expansion cooling" in the mixing-layer as the pressure P drops with altitude and Λ rad is the radiative cooling rate (both per unit volume). The contribution of H 2 formation is neglected in this analysis, as well as advection of thermal energy into the layer, since c 1 v w and c 2 v w . Under our isothermal hypothesis for the mixing-layer, the expansion cooling rate may be simply expressed as where x denotes position along the layer. We thus obtain where η is the normalized mass flux within the mixing-layer plotted in Figure 14. Therefore, in our model, the ratio of expansion cooling to viscous heating in the mixing-layer is independent of the turbulent entrainment efficiency α, and only scales with (c L /v w ) 2 . Furthermore, the remaining terms in this ratio are only weakly dependent on the values of λ and Λ (see Figures 9 and 14). On the typical scale z ≤ 0.1 r s encompassed by the Herschel /HIFI beam, we find that (Λ exp /Γ visc ) ≤ 200 (c L /v w ) 2 . With our fitted values of c L 1 km/s and v w 30 km/s for the broad component of SMM1, we infer that expansion cooling should be negligible with respect to viscous heating. Thus, we only need to compare the viscous heating rate with the radiative cooling rate. As noted by Kristensen et al. (2017, see their Figure 11), cooling by CO largely dominates over cooling by H 2 at temperatures of 250 K (for a standard CO/H 2 abundance ratio). We further assume that CO cooling is excited mainly by collisions with H 2 in the low-density limit (which we will verify a posteriori for SMM1). Denoting L 0 (T ) as the CO cooling rate coefficient (in erg s −1 cm 3 ) at temperature T , and X CO and X H2 as the CO and H 2 abundances relative to the total number density of H nuclei, n H = ρ L /(1.4 m H ), we have The ratio of turbulent heating to CO cooling is then independent of the mixing-layer density ρ L . With a mean layer velocity v L v w /2 (Couette flow), a typical value of 2Ṁ L within the Herschel beam of 2Π L /R b (see Eqn. 27), and 2Π L = Π BC where Π BC is the momentum in the warm CO broad component, this ratio can be expressed as It is remarkable that apart from X H2 there are no free parameters in this ratio, as all of the other factors are well constrained by observations of SMM1: The value of c L is fixed by the relative intensities of the high-J CO lines, indicating T L 250 K. The corresponding value of L 0 ∼ 3 × 10 −24 erg s −1 cm 3 at 250 K is set by molecular collision rate calculations (Neufeld & Kaufman 1993). The value of v w derives from our model fit to the CO(16-15) line profile in Figure 17. The cavity radius R 1000 au at z = R b = 2400 au derives from our model fitting of the outflow shape in Figure 18. The value of α α max is required for our model to be consistent with the dust envelope infall rate in SMM1 (see Section 5.4). Finally, the product Π BC X CO is equal to Π obs × (5 × 10 −5 ), where the value of Π obs = 8 × 10 −2 M km s −1 is fixed by the observed CO line profile intensity and excitation temperature in SMM1 (see Equation 26).
We conclude that if hydrogen is mostly in molecular form (X H2 0.5), and CO cooling is not far from the lowdensity regime, the ratio in Eqn. 38 is close to 1 for our mixing-layer model of SMM1, and thermal equilibrium can be maintained at the observed temperature 250 K of the broad CO J=16-15 component.
The low-density CO cooling expression applies only until L 0 n(H 2 ) 0.5 L LT E , where L LT E is the cooling rate per CO molecule in the high-density LTE regime. At 250 K, L LT E 10 −18 erg/s (Neufeld & Kaufman 1993) hence the validity extends to n(H 2 ) ≤ 1.7 × 10 5 cm −3 . To estimate the density in the SMM1 mixing-layer on the scale of the HIFI beam, we note that the shell pressure distribution on large scales is approximated by Eqn. 14. We infer the H nucleus density predicted in the mixing-layer at Z R b for the SMM1 model parameters If all hydrogen is in molecular form, we have n(H 2 ) = 0.5n H 0.7 × 10 5 cm −3 and the low-density regime of CO cooling assumed in Eqn. 38 is indeed justified for SMM1 on HIFI beam scales.

Summary and Discussion of the Model Fit to Protostar Serpens-Main SMM1
In summary, we have shown that our simple model of a turbulent mixing-layer across a static wind/envelope interface is able to reproduce successfully all of the observed properties of the broad CO outflow component discovered by Herschel /HIFI in the Serpens-Main SMM1 protostar, for a self-consistent and physically realistic set of parameters. The CO J=16-15 line profile shape and velocity extent are reproduced for a typical wind speed v w 30 km s −1 and a view angle of θ obs = 30 • to the disk plane (i.e. the median value expected for random inclinations). This wind speed is smaller than predicted for an X-wind from the innermost disk radius at 0.1 au (v w 150 km/s; see e.g. Shang, Shu, & Glassgold 1998) but remains compatible with a slow MHD disk wind launched from a few au in the disk (see e.g. Tabone et al. 2020). Next, the observed outflow cavity size on 300 − 3000 au scales, when combined with the estimated dust temperature in the envelope, requires a ratio of wind mass-flux to infall rate of 4%. With this imposed ratio, the observed CO-emitting momentum in the broad component (provided by wind entrainment) is consistent with the observed infall rate in the dusty envelope for a standard interstellar CO abundance and a turbulent entrainment coefficient α 0.03 (consistent both with our assumption of a Couette flow in the mixing-layer, and with mixing-layer laboratory experiments). The corresponding wind mass-flux then represents a fraction 0.06 − 0.2 of the disk accretion rate onto SMM1 (as determined from its bolometric luminosity), consistent with current disk wind models. Finally, the observed temperature in the broad CO outflow component of SMM1 is consistent with a balance between turbulent heating and CO cooling in the mixing-layer if H 2 is mostly in molecular form, which is very likely at such low temperatures. We also verify that the values of λ and Λ in our SMM1 model are consistent with the conditions for cavity breakout and the requirement of thin shells across the full range of disk radii expected in such a source, r d = 10-300 au.
A obvious next step for this modeling work would be to compute self-consistently the time-dependent evolution of temperature and chemistry through the wind shock and along the mixing-layer, using for example the molecular MHD disk wind models of Panoglou et al. (2012) and Yvart et al. (2016) as initial conditions. Such a calculation would provide an important check on our model requirement of an interstellar CO abundance in the mixing-layer of SMM1, to match independent constraints on the infall rate obtained from dust emission observations. It would also aid in the identification of the best tracer for the predicted narrow emission from the shocked wind layer (cf. black double-horned profile in Fig. 17).
Furthermore, since our model assumes a steady wind-blown cavity, it provides a natural explanation not only for the broad Herschel -bright CO component, but also for the narrow outflow cavity radii ≤ 3000 au observed at the Class 0 stage of 10 4 − 10 5 yrs, despite observed CO velocities on the order of 10 km/s. In contrast, for wind-driven shell models with full mixing, quasi radial outflow motions of the same amplitude lead to excessive cavity radii in only a few thousand years (see eg. Shang et al. 2006;López-Vázquez, Cantó, & Lizano 2019). Comparison over a larger sample of protostars with well characterized broad components and outflow cavities will be necessary to verify that self-consistent models can be found, as in SMM1, and to investigate how the required wide-angle wind properties would need to vary with source properties.

CONCLUSIONS
In this paper we have reconsidered the interaction of a wind expanding into a surrounding medium under the assumption of partial mixing across the boundary layer separating the shocked wind and envelope. Our solutions differ from conventional wind/envelope interaction models where instantaneous full mixing is assumed (eg. Li & Shu 1996;Lee et al. 2000;López-Vázquez, Cantó, & Lizano 2019) in that we produce static, rather than expanding, shells. To maintain the stationary shape, we allow the shocked and deflected wind to flow upward at close to v w along the interior of the cavity wall while the shocked and deflected envelope moves slowly downwards along the exterior of the cavity wall. A turbulent entrainment layer is thus able to form between these two deflected flows.
Specifically, we determine the shape of the stationary cavity formed when an isotropic wind interacts with an infalling and rotating (Ulrich 1976) envelope. The resulting model is then quantitatively compared with observations of the protostellar outflow from SMM1 in the Serpens Molecular Cloud.
The main results of our analysis are as follows: 1. The shape of the steady-state cavity (Section 2.3) is determined by two non-dimensional parameters, λ, the ratio of the wind ram pressure to the fiducial infall ram pressure (Eqn. 8), and Λ, the ratio of the wind ram pressure to the envelope thermal pressure at the edge of the disk (Eqn. 12) . We show that Λ sets the foot-point of the cavity at the disk plane ( Figure 2) and that breakout solutions require λ > 0.2, with the cavity shapes becoming self-similar for λ > 0.5 (Figure 3). In the self-similar regime, the size scaling of the cavity is determined by r s (Eqn. 10).
2. Under the assumption of no mixing (Section 3.1), the shocked and deflected wind moves upward along the cavity at close to the velocity v w , while the shocked and deflected envelope moves downward only slowly, except very near the base (Figure 7). Furthermore, away from the base the associated downward momentum flux is much less than the upward momentum flux (Figure 8).
3. Under the assumption of partial mixing within a turbulent layer between the upward and downward shocked deflected layers (Section 3.2), the overall amount of material brought into mixing-layer, from both sides, is directly proportional to the mass-loss rate in the wind multiplied by the entrainment efficiency α and v w /c L ( Figure 14). Furthermore, as previously shown by Cantó & Raga (1991), the mass entrainment from the upward, wind, side is roughly twice that of the downward, envelope, side (Figure 11), where the approximate proportionality is set by the assumption that across the mixing-layer the flow velocity profile is linear (i.e. a Couette flow).
4. The shape of the line profile produced by material flowing along the cavity wall strongly depends on which layer is responsible for the emission (Section 4.2). The upward, shocked wind layer moves fast, v ∼ v w , and has little curvature, resulting in a narrow profile peaked at the projected wind velocity. Alternatively, due to the Couette-type flow, emission from the mixing-layer is broad and peaks at rest velocity (Figures 16 and 17).
5. We find an excellent correspondence between the broad component of the CO J=16-15 line profile observed by Herschel towards the protostar Serpens-Main SMM1, and a mixing-layer model with v w = 25-30 km/s, a viewing angle 30 degrees from the disk plane, and an emissivity proportional to density (Section 5). Furthermore, taking α 0.03, a value which matches very well with experimental measurements of supersonic mixing, and assuming a standard CO abundance and a reasonable ratio of wind to infall rate of 4%, we find excellent quantitative agreement between the observed momentum in the CO broad component, the observed infall rate of SMM1, and the observed outflow cavity size (Section 5.4).
6. We compute the turbulent heating, expansion cooling, and radiative CO cooling within the mixing-layer, and show that their ratio is appropriate to keep the gas warm at the observed temperature T L 250 K in SMM1 (Section 5.6).
7. Finally, our model provides a natural explanation for the narrow outflow cavity radii observed at the Class 0 stage of 10 4 -10 5 yrs (Section 5.7). Unlike wind-driven shell models with full-mixing, in which radial motions quickly lead to large cavity sizes, our partial mixing solutions with a mixing layer separating the shocked wind and envelope produce a time-independent, steady cavity where observed velocities are parallel to the cavity walls, and thus do not lead to excessive expansion.
To summarize, we provide a model for the interaction between a wind and a surrounding envelope which potentially can be applied widely, from protostellar outflows to galactic-scale. The model produces steady-state cavities, broad line profiles peaked at the rest velocity, and constrains the turbulent entrainment efficiency. It therefore provides a new framework in which to interpret the observations of warm wind-driven outflows, and in particular to reconcile modest outflow cavity widths with the large observed flow velocities. While the model successfully reproduces a number of observational constraints for a single protostellar outflow, Serpens-Main SMM1, an obvious next step is to apply this analysis to a larger sample of protostellar sources in order to test its success; this will be done in a forthcoming publication. Envelope gas Figure A19:. The allowed angle of incidence of the wind-envelope interface at the disk plane as a function of R 0 , for models with different Λ values. The dotted black line indicates the angle of incidence of the infalling material at R 0 . Figure A20:. Cavity shapes near the base for an isotropic wind colliding with an Ulrich (1976) infalling envelope for different Λ values. The magenta and black lines correspond to the model of (Λ = 25, λ = 1/2) and (Λ = 200, λ = 10), respectively. For each model, the solid curve represents the 'fiducial' solution where the wind/envelope interface at the foot-point lies parallel to the local streamline of the material of the infalling envelope at the disk plane, whereas the dashed line represents a different solution where the base of the interface locates at a larger R 0 . For a given Λ, the solutions of different base position converge at a small distance from the disk plane (a few 0.01r d ). The thin dotted lines in the background indicate the streamlines of the material within the infalling envelope.
between the dashed line and each colored solid line in Figure 19). For a given Λ, when the interface foot-point is placed at somewhat larger R, the cavity quickly converges to the fiducial case within a small distance from the base, as is shown in Figure 20. The foot-point, however, cannot become arbitrarily close to r d without the infalling material crushing the wind and preventing a breakout solution, dependent on the value of λ. Thus the foot-point location is highly constrained, with larger λ values allowing a broader range of solutions at the base, all converging to self-similar solutions at altitude.

B. EQUATIONS FOR THE GROWTH OF THE MIXING-LAYER
Several typographical errors were present in the general equations from Raga, Cabrit, & Cantó (1995) describing the growth of the mixing-layer between two axisymmetric moving fluids of speed v 1 and v 2 < v 1 : in their equation (13), the term hP should have been −hP , while in their equation (16), the factor r c next to (h 1 + h 2 )P should not be present. Below, we reproduce their equations (15) and (16) where the latter typo has been corrected, and we use the subscript "L" to denote quantities in the mixing-layer, instead of the lower case letter "l", which was difficult to differentiate from the digit"1" in Raga, Cabrit, & Cantó (1995). Furthermore, for consistency with the notation in the main paper, here we refer to the velocity within the mixing-layer as v L whereas in Raga, Cabrit, & Cantó (1995) it is just v. All other notations are kept the same. Since we assume an isothermal mixing-layer with uniform sound-speed c L , we do not have to integrate their energy equation. Thus the system reduces to solving the following set of coupled equations for h 1 (x) and h 2 (x), which are the respective widths by which the mixing-layer encroaches into each fluid: In these equations, the subscript "1" denotes quantities pertaining to the fast fluid (in our case, the deflected wind), the subscript "2" pertains to the slow fluid (in our case, the deflected envelope), x is the distance along the flow, r c and P are the cylindrical radius and the pressure at the current point, primes denote derivatives versus x, and the mean velocities in the mixing-layer for a linear Couette flow are given by and The mass densities in the three layers are determined through transverse pressure equilibrium as ρ 1 c 2 1 = ρ 2 c 2 2 = ρ L c 2 L = P, where we assume here for simplicity that the sound speeds c 1 , c 2 and c L do not vary with position. An added complication for our paper is that our two fluids do not flow in the same direction. Fortunately, we always have v 2 < v d v 1 v w . We thus assume v 2 = 0 when integrating these equations upward along x, ie. that mass entrainment into the mixing-layer from the slow envelope side is largely dominated by the turbulent entrainment term.

C. VISCOUS DISSIPATION IN THE MIXING-LAYER
We consider the simplified case, relevant to the present paper, of an isothermal, supersonic mixing-layer of width h between two fluids with v 2 0 and v 1 v 2 , and with a linear velocity gradient across the flow direction (Couette profile) dv(y)/dy (v 1 /h) where v 1 changes weakly with position x along the flow.
We calculate Γ visc , the excess kinetic energy that needs to be locally dissipated by viscous turbulence per unit time and volume within the layer to maintain its internal linear Couette profile, as follows. The flux of kinetic energy flowing along the mixing-layer is given byĖ