MHD simulations of the formation and propagation of protostellar jets to observational length scales

TL;DR

This paper presents detailed 2.5-D MHD simulations of protostellar jets, revealing their structure, dynamics, and scaling laws, and demonstrating how magnetic fields influence jet formation and propagation to observable scales.

Contribution

It introduces a comprehensive simulation framework that models the formation and evolution of protostellar jets, highlighting the roles of magnetic fields and launching mechanisms in reproducing observed jet properties.

Findings

Jet structure includes a core, sheath, and ambient medium.

Power-law relationships between magnetic field strength and jet speed.

Jets are launched by magnetic tower and bead-on-a-wire mechanisms.

Abstract

We present 2.5-D global, ideal MHD simulations of magnetically and rotationally driven protostellar jets from Keplerian accretion discs, wherein only the initial magnetic field strength at the inner radius of the disc, $B_{\rm i}$, is varied. Using the AMR-MHD code AZEUS, we self-consistently follow the jet evolution into the observational regime ($>10^3\,\mathrm{AU}$) with a spatial dynamic range of $\sim6.5\times10^5$. The simulations reveal a three-component outflow: 1) A hot, dense, super-fast and highly magnetised 'jet core'; 2) a cold, rarefied, trans-fast and highly magnetised 'sheath' surrounding the jet core and extending to a tangential discontinuity; and 3) a warm, dense, trans-slow and weakly magnetised shocked ambient medium entrained by the advancing bow shock. The simulations reveal power-law relationships between $B_{\rm i}$ and the jet advance speed, $v_{\rm jet}$, the average jet rotation speed, $\langle v_\varphi\rangle$, as well as fluxes of mass, momentum, and kinetic energy. Quantities that do not depend on $B_{\rm i}$ include the plasma-$\beta$ of the transported material which, in all cases, seems to asymptote to order unity. Jets are launched by a combination of the 'magnetic tower' and 'bead-on-a-wire' mechanisms, with the former accounting for most of the jet acceleration---even for strong fields---and continuing well beyond the fast magnetosonic point. At no time does the leading bow shock leave the domain and, as such, these simulations generate large-scale jets that reproduce many of the observed properties of protostellar jets including their characteristic speeds and transported fluxes.

Tables

Table 1. Table 1: Physical quantities as measured/inferred from observations of protostellar outflows taken from references listed in the text. The quantities include from top to bottom: v jet subscript 𝑣 jet v_{\rm jet} , advance speed of the jet; ⟨ v rot ⟩ delimited-⟨⟩ subscript 𝑣 rot \langle v_{\rm rot}\rangle , jet rotational speed; M ˙ ˙ 𝑀 \dot{M} , mass flux; p ˙ ˙ 𝑝 \dot{p} , linear momentum flux; L ˙ ˙ 𝐿 \dot{L} , angular momentum flux; and v entr subscript 𝑣 entr v_{\rm entr} , advance speed of entrained outflow.

$v_{\mathrm{jet}}$	100 – 1000 km s-1
$⟨v_{\mathrm{rot}}⟩$	5 – 25 km s-1
$\dot{M}$	${10}^{-9}$ – ${10}^{-5}{M}_{\odot }$ yr-1
$\dot{p}$	${10}^{-6}$ – $1.4\times {10}^{-4}{M}_{\odot }$ yr-1 km s-1
$\dot{L}$	${10}^{-6}$ – ${10}^{-5}$ $M_{\odot }$ yr-1 AU km s-1
$v_{\mathrm{entr}}$	1 – 30 km s-1

Table 2. Table 2: Estimates made from local simulations of the observational parameters listed in Table 1 , taken from Ouyed & Pudritz ( 1997a ); Ustyugova et al. ( 1999 ); Anderson et al. ( 2005 ); Fendt ( 2009 ); Staff et al. ( 2010 ); Sheikhnezami et al. ( 2012 ); Stepanovs & Fendt ( 2014 ) .

$v_{\mathrm{jet}}$	35 – 1,300 km s-1
$⟨v_{\mathrm{rot}}⟩$	km s-1
$\dot{M}$	${10}^{-8}$ – ${10}^{-5}{M}_{\odot }$ yr-1
$\dot{p}$	—
$\dot{L}$	$3\times {10}^{-8}$ – $9\times {10}^{-6}{M}_{\odot }$ yr-1 AU km s-1
$v_{\mathrm{entr}}$	—

Table 3. Table 3: Initial static grids used in all simulations (refinement ratio ν = 2 𝜈 2 \nu=2 ). All zones are square, and thus Δ ​ z = Δ ​ r = Δ Δ 𝑧 Δ 𝑟 Δ \Delta z=\Delta r=\Delta . The near powers-of-two for the grid dimensions is a consequence of requiring the number of zones—including boundary zones—in each dimension of grids 1–8 to be a multiple of the number of OpenMP threads used (typically 16).

Level	$z_{\max }$ (AU)	$r_{\max }$ (AU)	$\Delta$ (AU)
1	4070.4	249.6	1.6
2	508.8	124.8	0.8
3	254.4	62.4	0.4
4	127.2	31.2	0.2
5	63.6	15.6	0.1
6	31.8	7.8	0.05
7	15.9	3.9	0.025
8	7.95	1.95	0.0125
9	3.975	0.975	0.00625

Table 4. Table 4: Summary of simulations. t end subscript 𝑡 end t_{\rm end} and z end subscript 𝑧 end z_{\rm end} are the time and jet length at simulation end, respectively.

Simulation	A	B	C	D	E	F	G	H
${\beta }_i$	0.1	0.4	1.0	2.5	10	40	160	640
$B_i$ (G)	200	100	63.2	40	20	10	5	2.5
$t_{\mathrm{end}}$ (yr)	47	64	77	88	121	153	153	153
$z_{\mathrm{end}}$ (AU)	4070	4070	4070	4070	4070	3800	2770	2380

Table 5. Table 5: Summary of measured quantities from the simulations. Mass-weighted average quantities include: ⟨ v φ ⟩ delimited-⟨⟩ subscript 𝑣 𝜑 \langle v_{\varphi}\rangle , the jet rotation speed (all speeds in km s -1 ); ⟨ v z ⟩ delimited-⟨⟩ subscript 𝑣 𝑧 \langle v_{z}\rangle , the mass-weighted axial speed along the last 250 AU of the jet; v jet subscript 𝑣 jet v_{\rm jet} , the advance speed of the jet tip into the ambient medium; v entr subscript 𝑣 entr v_{\rm entr} , the mass-weighted average poloidal velocity of the ambient medium entrained by the bow shock; M ext = v jet / c s , ext subscript 𝑀 ext subscript 𝑣 jet subscript 𝑐 s ext M_{\rm ext}=v_{\rm jet}/c_{\rm s,ext} , the sonic Mach number of the advance speed relative to asymptotic external ambient medium (where c s = 0.77 subscript 𝑐 s 0.77 c_{\rm s}=0.77  km s -1 for T asym = 22 subscript 𝑇 asym 22 T_{\rm asym}=22\, K); M int = ⟨ v z ⟩ / ⟨ c s , int ⟩ subscript 𝑀 int delimited-⟨⟩ subscript 𝑣 𝑧 delimited-⟨⟩ subscript 𝑐 s int M_{\rm int}=\langle v_{z}\rangle/\langle c_{\rm s,int}\rangle , the internal sonic Mach number, where ⟨ c s , int ⟩ delimited-⟨⟩ subscript 𝑐 s int \langle c_{\rm s,int}\rangle is the mass-weighted average sound speed over the last 500 AU of the jet, and M f , int = ⟨ v z ⟩ / ⟨ a f , int ⟩ subscript 𝑀 f int delimited-⟨⟩ subscript 𝑣 𝑧 delimited-⟨⟩ subscript 𝑎 f int M_{\rm f,int}=\langle v_{z}\rangle/\langle a_{\rm f,int}\rangle , the internal fast magnetosonic Mach number, where ⟨ a f , int ⟩ delimited-⟨⟩ subscript 𝑎 f int \langle a_{\rm f,int}\rangle is the mass-weighted average fast speed over the last 500 AU of the jet. Data from along the magnetic field line anchored at r 0 = 1 subscript 𝑟 0 1 r_{0}=1\, AU in the disc include: v p , A subscript 𝑣 p A v_{\rm p,A} , v p , f subscript 𝑣 p f v_{\rm p,f} , v p , max subscript 𝑣 p max v_{\rm p,max} , the poloidal speeds at, respectively, the Alfvén, fast, and asymptotic points; ( z A , r A ) subscript 𝑧 A subscript 𝑟 A (z_{\rm A},r_{\rm A}) , ( z f , r f ) subscript 𝑧 f subscript 𝑟 f (z_{\rm f},r_{\rm f}) , the coordinates (all distances are in AU) of the Alfvén and fast points respectively; and s × subscript 𝑠 s_{\times} , the location along the field line where 2 ​ v φ = r ​ ω 0 2 subscript 𝑣 𝜑 𝑟 subscript 𝜔 0 2v_{\varphi}=r\omega_{0} (Sect.  5.5 ). Fluxes measured inside the TD at z = 1 , 000 𝑧 1 000 z=1,\!000 AU include: M ˙ ˙ 𝑀 \dot{M} , mass flux ( M ⊙ subscript 𝑀 direct-product M_{\odot}  yr -1 ); P ˙ ˙ 𝑃 \dot{P} , linear momentum flux ( M ⊙ subscript 𝑀 direct-product M_{\odot}  yr -1 km s -1 ); L ˙ ˙ 𝐿 \dot{L} , angular momentum flux ( M ⊙ subscript 𝑀 direct-product M_{\odot}  yr -1 AU km s -1 ); and K ˙ ˙ 𝐾 \dot{K} , kinetic energy flux (erg s -1 ). Quantities which follow a power law in B i subscript 𝐵 i B_{\rm i} include an estimate of the power-law index, α 𝛼 \alpha , in the last column. Uncertainties for all quantities are given by the standard-deviation of the data time-averaged over a period of ∼ similar-to \sim 1 yr (or > 400 absent 400 >400 time steps on the coarsest level), or 1 in the last of three significant digits, whichever is the greater.

	A	B	C	D	E	F	G	H		$\alpha$
${\beta }_i$	0.1	0.4	1.0	2.5	10	40	160	640
$B_i$ (G)	200	100	63.2	40	20	10	5	2.5
mass-weighted averages:
$⟨v_{\phi }⟩$	$21.0\pm 0.3$	$12.6\pm 0.4$	$9.1\pm 0.5$	$6.6\pm 0.6$	$4.3\pm 0.5$	$2.8\pm 0.6$	$2.2\pm 0.5$	$1.8\pm 0.4$		$0.67\pm 0.03$
$⟨v_z⟩$	$419.\pm 1.$	$316.\pm 1.$	$267.\pm 1.$	$234.\pm 1.$	$143.\pm 1.$	$113.\pm 1.$	$82.0\pm 0.4$	$64.0\pm 0.1$		$0.43\pm 0.01$
$v_{\mathrm{jet}}$	$417.\pm 1.$	$309.\pm 1.$	$252.\pm 1.$	$230.\pm 1.$	$166.\pm 1.$	$116.\pm 1.$	$85.6\pm 0.1$	$74.2\pm 0.1$		$0.42\pm 0.01$
$v_{\mathrm{entr}}$	$11.5\pm 0.1$	$7.84\pm 0.1$	$5.93\pm 0.02$	$3.9\pm 0.2$	$2.5\pm 0.4$	$2.01\pm 0.03$	$1.46\pm 0.01$	$1.10\pm 0.06$		$0.56\pm 0.01$
$M_{\mathrm{ext}}$	$536.\pm 1.$	$397.\pm 1.$	$324.\pm 1.$	$296.\pm 1.$	$214.\pm 1.$	$149.\pm 1.$	$110.\pm 1.$	$99.\pm 1.$		$0.40\pm 0.02$
$M_{\mathrm{int}}$	$22.2\pm 0.2$	$26.4\pm 0.4$	$27.7\pm 0.4$	$20.9\pm 0.3$	$23.\pm 4.$	$9.5\pm 0.1$	$11.9\pm 0.2$	$9.2\pm 0.2$		—
$M_{f,\mathrm{int}}$	$2.00\pm 0.01$	$2.52\pm 0.01$	$2.93\pm 0.01$	$3.99\pm 0.02$	$4.3\pm 0.1$	$6.02\pm 0.01$	$7.46\pm 0.01$	$6.29\pm 0.02$		—
$r_0=1$AU field line:
$v_{p,A}$	$161.\pm 1.$	$101.\pm 1.$	$73.0\pm 0.1$	$51.2\pm 0.2$	$29.6\pm 0.1$	$14.9\pm 0.1$	$7.5\pm 0.2$$†$	n/a$‡$		$0.76\pm 0.02$
$v_{p,f}$	$204.\pm 1.$	$129.\pm 1.$	$94.6\pm 0.1$	$69.5\pm 0.5$	$43.9\pm 0.1$	$27.9\pm 0.2$	$18.5\pm 0.2$	n/a		$0.667\pm 0.001$
$v_{p,\max }$	$297.\pm 1.$	$219.\pm 1.$	$166.\pm 1.$	$129.\pm 1.$	$87.1\pm 0.1$	$60.2\pm 0.1$	$33.2\pm 0.4$	n/a		$0.55\pm 0.01$
$r_A$	$11.6\pm 0.1$	$7.24\pm 0.01$	$5.46\pm 0.01$	$4.15\pm 0.01$	$2.84\pm 0.01$	$2.04\pm 0.01$	$1.66\pm 0.04$	n/a		$0.62\pm 0.02$
$r_f$	$18.9\pm 0.1$	$12.0\pm 0.1$	$9.16\pm 0.09$	$7.0\pm 0.1$	$5.06\pm 0.01$	$4.18\pm 0.01$	$4.4\pm 0.3$	n/a		$0.54\pm 0.02$
$z_A$	$106.\pm 1.$	$33.1\pm 0.1$	$17.0\pm 0.1$	$8.78\pm 0.01$	$3.26\pm 0.01$	$1.22\pm 0.01$	$0.51\pm 0.04$	n/a		$1.47\pm 0.02$
$z_f$	$315.\pm 2.$	$141.\pm 1.$	$68.\pm 2.$	$34.\pm 2.$	$13.3\pm 0.1$	$6.79\pm 0.04$	$4.9\pm 0.5$	n/a		$1.41\pm 0.04$
$s_{\times }$	$19.3\pm 0.1$	$8.56\pm 0.02$	$5.10\pm 0.01$	$3.12\pm 0.01$	$1.68\pm 0.01$	$0.982\pm 0.009$	$0.64\pm 0.01$	n/a		$1.10\pm 0.04$
jet fluxes at $z=1,000$ AU:
$\dot{M}$ $(\times {10}^{-6})$	$2.0\pm 0.1$	$1.5\pm 0.1$	$1.36\pm 0.07$	$1.16\pm 0.04$	$1.01\pm 0.06$	$0.82\pm 0.05$	$0.51\pm 0.04$	$0.039\pm 0.006$		$0.32\pm 0.03$
$\dot{p}$ $(\times {10}^{-4})$	$4.9\pm 0.6$	$2.8\pm 0.3$	$2.0\pm 0.1$	$1.37\pm 0.06$	$0.83\pm 0.06$	$0.46\pm 0.03$	$0.22\pm 0.02$	$0.017\pm 0.003$		$0.82\pm 0.03$
$\dot{L}$ $(\times {10}^{-4})$	$14.\pm 2.$	$5.\pm 1.$	$3.6\pm 0.4$	$1.8\pm 0.1$	$0.97\pm 0.04$	$0.36\pm 0.03$	$0.14\pm 0.01$	$0.0024\pm 0.0005$		$1.24\pm 0.05$
$\dot{K}$ $(\times {10}^{33})$	$96.\pm 18.$	$37.\pm 7.$	$20.\pm 2.$	$11.2\pm 0.7$	$4.6\pm 0.4$	$1.7\pm 0.1$	$0.60\pm 0.06$	$0.053\pm 0.009$		$1.37\pm 0.02$

Table 6. Table 6: Comparison of the speed at the fast point as predicted by Eq. ( 8 ) and as measured directly along the 1 AU field line.

$v_{p,f}$	A	B	C	D	E	F	G
Eq. (8)	205.8	127.3	94.10	69.43	43.94	27.6	19.4
1 AU field line	203.5	128.8	94.64	69.45	43.86	27.9	18.5
% difference	1.1	1.2	0.57	0.03	0.2	1.1	4.5

Table 7. Table 7: Estimates made from our global simulations for the observational parameters listed in Table 1 . Ranges for M ˙ ˙ 𝑀 \dot{M} , p ˙ ˙ 𝑝 \dot{p} , L ˙ ˙ 𝐿 \dot{L} , and K ˙ ˙ 𝐾 \dot{K} measured from the last 10 yr of data for each simulation at a height of z = 1 , 000 𝑧 1 000 z=1,\!000 AU above the disc. ‘Jet’ values are defined as regions interior to the TD (masked using M ≥ 5 𝑀 5 M\geq 5 ), while ‘entrained’ material is defined as the shocked ambient just beyond the TD (masked using M < 5 𝑀 5 M<5 ).

	Jet
$v_{\mathrm{jet}}$	74 – 420 km s-1
$⟨v_{\phi }⟩$	1.8 – 21 km s-1
$\dot{M}$	$2.9\times {10}^{-8}$ – $2.3\times {10}^{-6}$ $M_{\odot }$ yr-1
$\dot{p}$	$1.3\times {10}^{-6}$ – $6.8\times {10}^{-4}$ $M_{\odot }$ yr-1 km s-1
$\dot{L}$	$1.5\times {10}^{-7}$ – $1.7\times {10}^{-3}$ $M_{\odot }$ yr-1 AU km s-1
$\dot{K}$	$3.9\times {10}^{31}$ – $1.6\times {10}^{35}$ erg s-1
	Entrained
$v_{\mathrm{entr}}$	1.1 – 11.5 km s-1
$⟨v_{\phi }⟩$	—
$\dot{M}$	$4.8\times {10}^{-9}$ – $2.3\times {10}^{-7}$ $M_{\odot }$ yr-1
$\dot{p}$	$1.6\times {10}^{-8}$ – $1.9\times {10}^{-5}$ $M_{\odot }$ yr-1 km s-1
$\dot{L}$	$7.7\times {10}^{-10}$ – $1.4\times {10}^{-4}$ $M_{\odot }$ yr-1 AU km s-1
$\dot{K}$	$4.2\times {10}^{28}$ – $1.0\times {10}^{33}$ erg s-1

Equations

$$\eta (\psi )$$

$$l(\psi )$$

$$\Omega (\psi )$$

$$\epsilon (\psi )$$

$$\Omega (\psi )=\frac{v_{\mathrm{K},0}}{r_0}\mathrm{and}\epsilon (\psi )=-\frac{3v_{\mathrm{K},0}^2}{2},$$

$$\eta (\psi )=\sqrt{\frac{{\rho }_{\mathrm{A}}}{4\pi }}\mathrm{and}l(\psi )=r_{\mathrm{A}}^2\Omega (\psi )=v_{\mathrm{K},0}\frac{r_{\mathrm{A}}^2}{r_0}.$$

$$2\epsilon -2\varphi +r^2{\Omega }^2=v_{\mathrm{K},0}^2\left(-3+2\frac{r_0}{R}+\frac{r^2}{r_0^2}\right)=v_p^2+M_{\mathrm{A}}^2{a}_{\phi }^2=\frac{v_{\mathrm{p}}^2{a}^2}{a_{\mathrm{p}}^2},$$

$$v_{\mathrm{p},\mathrm{f}}=\sqrt{v_{\mathrm{K},0}{B}_{\mathrm{p},\mathrm{f}}}{\left[\frac{1}{4\pi {\rho }_{\mathrm{f}}}\left(\frac{r_{\mathrm{f}}^2}{r_0^2}+\frac{2r_0}{R_{\mathrm{f}}}-3\right)\right]}^{1/4},$$

$$v_{\mathrm{p},\mathrm{f}}\sim \sqrt{B_{\mathrm{p},\mathrm{f}}}\sim B_{\mathrm{i}}^{1/2},$$

$$v_{\mathrm{p},\max }={\left(\frac{{\Omega }^2(\psi )r_0^2{B}_{\mathrm{p}}(r_0)}{4\pi \eta (\psi )}\right)}^{1/3}\sim B_{\mathrm{i}}^{2/3},$$

$${\partial }_r{B}_{\phi }\frac{\Delta }{B}>\frac{1}{N},$$

$$\nabla p+\rho \nabla \varphi =0,$$

$$\text{g}=-\frac{1}{{\rho }_{\mathrm{h}}}\nabla p_{\mathrm{h}},$$

$${\rho }_{\mathrm{h}}={\rho }_{\mathrm{i}}{\left(\frac{r_{\mathrm{i}}}{\sqrt{r^2+z^2}}\right)}^{\frac{1}{\gamma -1}}\text{and}{p}_{\mathrm{h}}=p_{\mathrm{i}}{\left(\frac{{\rho }_{\mathrm{h}}}{{\rho }_{\mathrm{i}}}\right)}^{\gamma }.$$

$${\rho }_{\mathrm{h}}={\rho }_{\mathrm{i}}\{\begin{array}{ll}{\left(\frac{r_{\mathrm{i}}}{\sqrt{r^2+z^2}}\right)}^{\frac{1}{\gamma -1}},& r^2+z^2\ge {R_{\mathrm{sph}}}^2;\\ {\left(\frac{r_{\mathrm{i}}}{R_{\mathrm{sph}}}\frac{3{R_{\mathrm{sph}}}^2-r^2-z^2}{2{R_{\mathrm{sph}}}^2}\right)}^{\frac{1}{\gamma -1}},& r^2+z^2<{R_{\mathrm{sph}}}^2.\end{array}$$

$$A_{\phi }=A_0\frac{\sqrt{r^2+(z+z_{\mathrm{d}}{)}^2}-(z+z_{\mathrm{d}})}{r},$$

$$B_z=\frac{1}{r}{\partial }_r(r{A}_{\phi })\text{and}{B}_r=-{\partial }_z{A}_{\phi },$$

$$A_0=\frac{B_{\mathrm{i}}{r}_{\mathrm{i}}}{\sqrt{2-\sqrt{2}}},$$

$$B_{\mathrm{i}}=\sqrt{\frac{8\pi p_{\mathrm{i}}}{{\beta }_{\mathrm{i}}}}.$$

$$v_r=v_z{B}_r/B_z;v_{\phi }=v_{\mathrm{K}}+v_z{B}_{\phi }/B_z.$$

$$\begin{array}{rlrl}& & E_z(-z)& =E_z(z);\\ E_r(0)& =-v_{\mathrm{K}}{B}_z(0);& E_r(-z)& =2E_r(0)-E_r(z);\\ E_{\phi }(0)& =0;& E_{\phi }(-z)& =-E_{\phi }(z),\end{array}\}$$

$$c_{\mathrm{s}}^2=\gamma \frac{p}{\rho }=\left(\gamma -1\right)\frac{G{M}_{\ast }}{\sqrt{r^2+z^2}}=\left(\gamma -1\right)v_{\mathrm{K}}^2.$$

$$p_{\mathrm{i}}$$

$$\frac{{\rho }_{\mathrm{i}}}{⟨m⟩}$$

$$T_{\mathrm{i}}$$

$$c_{\mathrm{s},\mathrm{i}}$$

$${\tau }_{\mathrm{i}}$$

$$F_r=J_{\phi }{B}_z-J_z{B}_{\phi }\sim -\frac{B_{\phi }}{r}{\partial }_r(r{B}_{\phi }),$$

$$c_{\mathrm{s},\infty }=c_{\mathrm{s},\mathrm{i}}\sqrt{\frac{T_{\infty }}{T_{\mathrm{i}}}}=c_{\mathrm{s},\mathrm{i}}{\left(\frac{{\rho }_{\infty }}{{\rho }_{\mathrm{i}}}\right)}^{\frac{\gamma -1}{2}}=c_{\mathrm{s},\mathrm{i}}(10^{-6}{)}^{1/3}\sim 0.77{\mathrm{km}\mathrm{s}}^{-1},$$

$$⟨{\beta }_{\mathrm{tot}}⟩=\frac{8\pi ⟨p⟩}{⟨B_{\phi }^2+B_{\mathrm{p}}^2⟩},$$

$$f(z_{\mathrm{f}},r_{\mathrm{f}})\equiv {\left[\frac{1}{4\pi {\rho }_{\mathrm{f}}}\left(\frac{r_{\mathrm{f}}^2}{r_0^2}+\frac{2r_0}{R_{\mathrm{f}}}-3\right)\right]}^{1/4},$$

$$⟨v_{\phi }⟩=\frac{{\int }_V\rho v_{\phi }dV}{{\int }_V\rho dV},$$

$$⟨v_z⟩=\frac{{\int }_V\rho v_{z}dV}{{\int }_V\rho dV},$$

$$v_{\mathrm{jet}}\propto ⟨v_{\phi }{⟩}^{2/3}.$$

Full text

††pubyear: ?††pagerange: MHD simulations of the formation and propagation of protostellar jets to observational length scales–MHD simulations of the formation and propagation of protostellar jets to observational length scales

1 Introduction

One of the most important epochs in the early evolution of most stars is the short period during which it throws a small fraction of the accreted gas back into the interstellar medium (ISM) as a pair of collimated, bipolar, supersonic jets. This period, lasting typically $10^{4}$–$10^{5}$ years, is to the star’s entire lifetime what a few hours is to a human’s. Yet, in this single ‘afternoon’, the protostellar system manages to shed itself of sufficient angular momentum to enable significant accretion from the protoplanetary accretion disc onto the protostar. Without protostellar jets, stars as we know them would not exist.

Thought to be stars in their own right when first observed by Burnham (1890), protostellar jets only started to be appreciated for what they are by Snell et al. (1980). In this seminal work, a detailed bipolar outflow model for L1551 is described that is still basic to the modern view, and done without ever using the word jet111While the term ‘jet’ was first used in an astrophysical context by Baade & Minkowski (1954) in describing the optical ‘protrusion’ on M87, it did not enter the protostellar vernacular until Mundt & Fried (1983).. It is now known that protostellar jets reach lengths of 0.1–5 pc (Bally et al., 2007) and can transport $\sim\!10$% of the accreted mass and $\sim\!70$% of the angular momentum (Woitas et al., 2005) out of the protostellar system and back into the ISM. They can appear as straight, ballistic, highly supersonic flows complete with leading bow shocks (e.g. HH 34; Devine et al. 1997), or more like effluent from a smokestack; wide, twisted and with no particular evidence of a supersonic nature (e.g. HH 47; Hartigan et al. 2005).

Reipurth (1999) and Wu et al. (2004) catalogue some 1,000 Herbig-Haro (HH) objects and molecular outflows from protostellar objects, imaged with atomic (e.g., H$\alpha$, Oiii, Sii) and/or molecular (e.g., CO) line emission and from which important kinematical and dynamical quantities are measured. While observational properties of individual outflows vary widely, ranges of values for parameters most useful for constraining numerical magnetohydrodynamical (MHD) simulations can be established. These include the advance speed of the jet into the ISM ($v_{\rm jet}$), the speed of entrained material swept up by the bow shock leading the jet ($v_{\rm entr}$), the average rotational speed of jet material about its propagation axis ($\langle v_{\rm rot}\rangle$), fluxes in mass ($\dot{M}$), linear momentum ($\dot{p}$), and angular momentum ($\dot{L}$). Table 1 displays a summary of what is currently known of these parameters as reported by Hartigan et al. (1994), Reipurth & Bally (2001), Podio et al. (2006), McKee & Ostriker (2007), Ray et al. (2007), Coffey et al. (2008); Coffey et al. (2011) and Frank et al. (2014). While these observations provide an extensive and highly detailed picture of protostellar jets, the quantities listed in Table 1 are largely inferred and indirect, and should not be taken as hard limits.

For example, direct evidence of jet rotation remains somewhat controversial since only recently has observational resolution been sufficient (on the order of 10 AU; e.g. Coffey et al., 2008; Bjerkeli et al., 2016; Lee et al., 2017) to yield reliable radial profiles across the jet. Woitas et al. (2005) report line-of-sight velocity gradients which they interpret as rotation, though Soker (2005) suggest this might indicate an interaction of the jet with a warped disc, while Fendt (2011) suggests MHD shocks in a helical field. Still, it is widely believed that protostellar jets must rotate if they are to succeed in their presumed task of ridding the protostar of its angular momentum.

An important physical quantity missing from Table 1 is the magnetic field strength. Direct measurements of B in a protostellar jet remain elusive, with just two indirect measures reported to date (Ray et al., 1997; Carrasco-González et al., 2010) which, almost by definition, represent extreme cases. Still, the theoretical evidence for magnetic fields pervading protostellar jets is overwhelming, and it is nearly universally accepted as being a critical ingredient to jet dynamics (Hartigan et al., 2007). Certainly, strong fields are known to exist within the inner regions of protostellar discs ($\sim\!1$ kG; Donati et al., 2005), and it is difficult to imagine how this is not transported outward by the jet.

To a large extent, the base of a jet can be characterised by the presence of strong gravitational and magnetic fields, and rapid rotation. In such an environment, Blandford & Payne (1982), based on the ‘bead-on-a-wire’ model first suggested by Henriksen & Rayburn (1971), showed that the formation of a super-fast jet is virtually inevitable. In their model, a Keplerian disc is threaded with ‘frozen in’ vertical magnetic flux. As the disc rotates, magnetic field lines are twisted and, once the angle at which they emerge from the disc falls below the critical value of $60^{\circ}$, the centrifugal force overwhelms gravity and drives material outward as ‘beads sliding along a rotating wire’.

This model, also known as ‘magneto-centrifugal driving’, has precipitated a plethora of numerical simulations to investigate its consequences. Since a 2-D axisymmetric, ideal Keplerian disc is unstable to the magneto-rotational instability (MRI, Balbus & Hawley, 1992 and references therein; only in 3-D is the instability saturated, Stone et al., 1996), most simulations of magneto-centrifugally launched jets (e.g. Uchida & Shibata, 1985; Ustyugova et al., 1995, 1999; Meier et al., 1997; Ouyed & Pudritz, 1997a, b, 1999; Krasnopolsky et al., 1999, 2003; Fendt & Čemeljić, 2002; Vitorino et al., 2002; von Rekowski et al., 2003; Ouyed et al., 2003; Anderson et al., 2005, 2006; Porth & Fendt, 2010; Stute et al., 2014; Teşileanu et al., 2014; Staff et al., 2010, 2015) treat the accretion disc as a boundary condition, allowing the jet dynamics to be studied independently of the disc.

There are a number of 2-D studies which do include the disc as part of the simulations, even if in a somewhat idealised fashion (e.g. Casse & Keppens, 2002, 2004; Zanni et al., 2007; Tzeferacos et al., 2009; Murphy et al., 2010; Sheikhnezami et al., 2012; Fendt & Sheikhnezami, 2013; Stepanovs & Fendt, 2014, 2016; Suriano et al., 2017; Zhu & Stone, 2017; Bai, 2017; Suriano et al., 2018). These simulations typically use a magnetic resistivity to prevent excessive disc turbulence, and are more realistic by including the disc evolution self-consistently. However, they are much more expensive computationally because of the significantly shorter physical time scales in the disc and it is because of this we have chosen here to treat the disc as a boundary condition.

Because most of the magneto-centrifugal ‘action’ occurs near the inner radius of the disc, simulations must be performed at a resolution of 0.01 AU or less in order to resolve the important physics there. Thus, even the most ambitious of the works listed above have followed the jet to just 100 AU (Anderson et al. 2005), and more recently to 150 AU (Stepanovs & Fendt, 2014, 2016), above the disc. Accordingly, we refer to these simulations collectively as ‘local’ simulations.

There also exists a class of ‘global’ simulations which follow the gravitational collapse of isolated, magnetised molecular cloud cores including the formation of the protostar and accretion disc (e.g. Seifried et al., 2011, 2012; Tomida et al., 2013; Tomida et al., 2015; Masson et al., 2016; Kölligan & Kuiper, 2018). While these simulations include, by design, observational length scales, due to the extreme computational costs involved, they cannot include the sub-0.01 AU scales necessary to suitably resolve the physics of the jet launching mechanism for any substantial length of time.

Notably, the large scale difference between local simulations and observed jets ($10^{3}$–$10^{6}$ AU; e.g. Devine et al. 1997; Aso et al. 2015) makes direct comparisons impossible, and to make any comparison at all one must make severe assumptions on how local variables relate to global properties of the jet. As an example, most local simulations continue their calculations long after the leading bow shock or Alfvén wave has left the grid. To say nothing of the change to the (thermo)dynamics of the jet that the sudden loss of a confining bow shock must cause, a proxy for $v_{\rm jet}$ must be used. Typically, this is the speed at the Alfvén point ($v_{\rm A}$) or, if it is still in the domain, the fast magnetosonic point ($v_{\rm f}$). In the only simulation performed to date where the jet launching conditions are controlled and the leading bow shock remains within the computational domain (Ramsey & Clarke, 2011), we find that the jet continues to accelerate well beyond the fast point, and thus $v_{\rm A}$ and $v_{\rm f}$ are poor proxies for the final $v_{\rm jet}$.

Table 2 summarises estimates of the parameters in Table 1 made from the local simulations cited in the caption. Because of the assumptions and extrapolations inherent in these estimates, we offer them only as an ‘order-of-magnitude’ check with the observations. Notably, we are not aware of any estimates of $\dot{p}$ that can be gleaned from local simulations.

This work is a continuation of Ramsey & Clarke (2011). Here, we present eight 2.5-D axisymmetric global simulations in which the jet is followed from its launching point with $0.00625$ AU resolution to a length of up to 4,000 AU, well into the observational regime. Even still, this represents only about 1% of the age and length of the largest jets from class 0/I young stellar objects which, as we will see, puts some limitations on what can be inferred.

As ‘immature’ as our simulations may be, they are still global in nature (both resolving the region where the jet is launched, and following the jet to observational scales), and imply a dynamic range in length scale of $\sim 6.5\times 10^{5}$ for our most highly resolved simulation. A single-grid (4,$096\times 256$ AU) 2-D MHD simulation with a resolution of $0.00625$ AU would require $>100$ billion zones over 100 million time steps to complete. Thus, to perform these simulations, we have used the adaptive mesh refinement (AMR) MHD code, AZEuS (Ramsey et al., 2012). Each simulation differs from the others only in the strength of the magnetic field at the inner radius of the disc, $B_{\rm i}$, which is used to scale an initially force-free, global ‘hour-glass’ magnetic field distribution. The central gravitating mass ($0.5\,{\rm M}_{\odot}$), as well as the parameters governing the Keplerian disc (initially in gravito-centrifugal balance) and the coronal atmosphere (initially in hydrostatic balance) are the same for all simulations. The purpose of this study is to determine if this is sufficient to produce a jet with the right observational characteristics, and further, what role if any $B_{\rm i}$ has on determining observational and physical properties of the jet. Neither of these fundamental questions can be answered by local simulations.

Finally, a comment on the choice of axisymmetry is in order. Even with AMR and distributing the calculations over 16–24 CPU cores, some of the simulations discussed herein required more than six months to complete, and a fully 3-D treatment was simply impractical. Buoyed by the knowledge that many stellar jets appear axisymmetric (e.g. HH 34; Devine et al., 1997), a 2-D axisymmetric approach was adopted at the outset of this project.

Still, the cost in realism is undeniable. Even in systems with a high degree of apparent axisymmetry, the fluid is subject to all modes of Kelvin-Helmholtz (K-H) instabilities (e.g. Hardee & Clarke, 1995) on both the large- and small-scale, which take their toll as the jet propagates. On the large scale, Clarke (1993) showed that the otherwise perfectly stable nose-cone found in a 2-D axisymmetric magnetically-confined jet was periodically sloughed off to the side in 3-D, resulting in a blunter, more slowly propagating jet. On the smaller scale, (Clarke, 1996a) found that, in simulations of a propagating jet with a weak magnetic field, the 3-D jet once again propagated more slowly and formed a blunter bow shock. In this case, the numerous modes of K-H instabilities in 3-D—all but one unavailable in 2-D axisymmetry—effectively converts directed kinetic energy of outflow to turbulent and ultimately thermal energy in the expanding cocoon. Such considerations should therefore be borne in mind as discussion of the present simulations unfolds.

In Sect. 2, we review some of the relevant steady state theory which applies to portions of our numerical solutions. In Sect. 3, we describe briefly the numerical methodology and how the simulations are initialised. Sections 4 and 5 comprise the bulk of the paper in which the simulations are described and analysed in detail. Finally, conclusions are drawn in Sect. 6.

2 Steady state analysis

Analogous to Bernoulli’s constant for hydrodynamics, in a steady state ($\partial_{t}=0$)222Throughout this paper, we use the ‘abbreviated Leibniz notation’ for derivatives. Thus, $\partial_{t}\equiv\partial/\partial t$., axisymmetric ($\partial_{\varphi}=0$), ideal MHD fluid, there are four conserved quantities along a given magnetic field line, which themselves are contours of the flux function333$\psi=rA_{\varphi}$, where $A_{\varphi}$ is the toroidal component of the vector potential., $\psi$. In Gaussian cgs units, these are (Weber & Davis, 1967; Mestel, 1968; Pudritz & Norman, 1983; Pelletier & Pudritz, 1992; Spruit, 1996):

[TABLE]

where $\eta$ is the mass load, $l$ and $\Omega$ are, respectively, the specific angular momentum and angular speed of the field line (including terms describing the magnetic torque)444$\Omega(\psi)$ is sometimes referred to as the iso-rotation parameter (e.g. Fendt & Memola, 2001; Porth et al., 2011)., and $\varepsilon$ is the specific energy of the fluid. In axisymmetry, $\psi$ describes surfaces of constant magnetic flux, and eqs. (1)–(4) are therefore also constant on these flux surfaces. Henceforth, we refer to these as the ‘Weber-Davis (WD) constants’. Note that Eq. (4) is essentially Bernoulli’s constant generalised for MHD. Other variables include the density, $\rho$, the poloidal velocity, $v_{\rm p}$, the toroidal velocity, $v_{\varphi}$, the poloidal magnetic field, $B_{\rm p}$, the Alfvén Mach number, $M_{\rm A}=v_{\rm p}/a_{\rm p}$, the poloidal Alfvén speed, $a_{\rm p}=B_{\rm p}/\!\sqrt{4\pi\rho}$, the toroidal Alfvén speed, $a_{\varphi}=B_{\varphi}/\!\sqrt{4\pi\rho}$, the toroidal magnetic field, $B_{\varphi}$, the adiabatic sound speed, $c_{\rm s}=\!\sqrt{\gamma p/\rho}$, the ratio of specific heats, $\gamma$, and the gravitational potential of the protostar, $\phi=-GM_{*}/R$. Here, $R=\!\sqrt{z^{2}+r^{2}}$, is the spherical polar radial coordinate, $z$, $r$ and $\varphi$ are the cylindrical coordinates555Note the difference here between $\phi$, the gravitational potential, and $\varphi$, the cylindrical coordinate.. Equations (1)–(4) are a straight-forward extension from Eqs. (12, 13, 19, and 31) in Spruit (1996), where we set $f^{\prime}=\Omega$.

As an example, by definition, the mass flux (${\cal F}_{\rho}=\rho v_{\rm p}\delta A$) is conserved along a streamtube of cross section $\delta A$ and the magnetic flux (${\cal F}_{B}=B_{\rm p}\delta A$) is conserved along a magnetic flux tube. In the steady state when $v_{\rm p}\parallel B_{\rm p}$, streamlines are everywhere parallel to magnetic flux lines, and $\eta(\psi)={\cal F}_{\rho}/{\cal F}_{B}$ is constant along a field line. Arguments establishing the constancy of $l$, $\Omega$, and $\varepsilon$ along field lines in the steady state follow similar lines.

Following Spruit (1996), if $r_{0}(\psi)$ is the radial coordinate of the disc where a particular field line, $\psi$, is anchored (Fig. 1), one can evaluate $\Omega(\psi)$ and ${\varepsilon}(\psi)$ at the anchor point. At $(z,r)=(0,r_{0})$, $v_{\rm p}\sim 0\,\Rightarrow\,M_{\rm A}\sim 0$, $v_{\varphi}=v_{\rm K,0}=\sqrt{GM_{*}/r_{0}}$ (assuming the disc is Keplerian), $\phi(r_{0})=-v_{\rm K,0}^{2}$, and Eqs. (3) and (4) reduce to:

[TABLE]

assuming a cold fluid ($\beta\ll 1;c_{\rm s}\sim 0$). Conversely, $\eta$ and $l$ are most conveniently evaluated at the Alfvén point ($M_{\rm A}=1$) where:

[TABLE]

Local simulations treat the leading Alfvén torsional wave and bow shock as transients and, irrespective of the dynamical consequences, allow them to leave the computational domain. Thereafter, most local simulations reach some sort of steady state from which various comparisons with analytical theory are made. For example, being in a near-steady state, $\eta(\psi)$ from Eq. (1) is expected to be constant, and thus many investigations use $\eta(\psi)$ as a parameter to specify the nature of the outflow (e.g. Ouyed & Pudritz 1999; Anderson et al. 2005); if $\rho v_{\rm p}$ and $B_{\rm p}$ are changed in proportion to each other in a steady state jet, then the character of the outflow should remain unaltered.

Local simulations reaching a steady state typically show that the jet speed saturates at or just beyond the fast point, assuming that this point remains inside the grid. Those that report on asymptotic jet speeds find $v_{\rm p,max}\sim\eta^{\alpha}$, with the majority finding $\alpha<0$ (i.e. flow speed increases as the field strength increases or the poloidal momentum at the disc decreases; e.g. Anderson et al. 2005; Zanni et al. 2007; Porth & Fendt 2010), as is expected from steady state theory (Spruit, 1996)666A notable exception is Ouyed & Pudritz (1999), who find the opposite trend.. Indeed, it should come as no surprise that local simulations confirm various aspects of steady state theory, since the assumptions of no transients is common to both. As soon as $v_{\rm p}\parallel B_{\rm p}$ is realised over much of the computational domain, the conclusions of steady state theory become inescapable.

Since the leading bow shock and Alfvén torsion wave never leave the grid in our simulations, none reach a global steady state. However, regions near the disc of some simulations (more so for stronger $B_{\rm i}$) do reach (locally) a quasi-steady state (as confirmed by the constancy of $\eta$, $l$, $\Omega$, and $\varepsilon$ along field lines), and we exploit this observation in some of the analysis. Under the assumption that portions of the jet are in steady state, we can determine how the flow speed at the fast point, for example, varies with $B_{\rm i}$, and then attempt to relate this to $v_{\rm jet}$.

To this end, from Eqs. (3), (4), and (5), we can write:

[TABLE]

where $a^{2}=a_{\varphi}^{2}+a_{\rm p}^{2}=B^{2}/(4\pi\rho)$ is the fast speed squared when $c_{\rm s}^{2}=0$ (cold flow). Thus, at the fast point where $(z,r)=(z_{\rm f},r_{\rm f})$, $a_{\rm p}=a_{\rm p,f}=B_{\rm p,f}/\!\sqrt{4\pi\rho_{\rm f}}$, and $v_{\rm p}=a=v_{\rm p,f}$, Eq. (7) becomes:

[TABLE]

where $r_{\rm f}$ and $R_{\rm f}$ are, respectively, the cylindrical and spherical polar radial coordinates to the fast point, as shown in Fig. 1. We have verified this formula directly from our simulations, and find agreement to better than 3% along the field line anchored at $r_{0}=1$ AU with measures from other field lines in steady state giving similar results (e.g., Table 6 on page 6).

If we choose the same field line foot print, $r_{0}$, for each simulation, $v_{\rm K,0}$ becomes a constant and, in as much as the quantity in square brackets to the 1/$4^{\rm th}$ power in Eq. (8) depends weakly on $B_{\rm i}$, we might expect:

[TABLE]

since the magnetic field profile scales with $B_{\rm i}$. This can be contrasted with the asymptotic flow speed predicted for steady state flow and a purely radial magnetic field (e.g. Spruit 1996; Eq. 74):

[TABLE]

since $\eta(\psi)\sim B_{\rm p}^{-1}$. We return to these predicted dependencies on $B_{\rm i}$ in Sect. 5.3.

3 Numerical considerations

3.1 AZEuS

The simulations presented herein are performed with the adaptive mesh refinement (AMR) MHD code, AZEuS (Adaptive Zone Eulerian Scheme; Ramsey et al. 2012; http://people.virginia.edu/~jpr8yu/azeus), based on Version 3.6 of ZEUS-3D (Clarke, 1996b, 2010; http://www.ica.smu.ca/zeus3d). The ZEUS family of codes is among the best tested, documented, and most widely used astrophysical MHD codes available. Our version allows one to choose to solve the internal energy or total energy equations, the latter being conservative in energy to machine round-off error. As the simulations presented here are only mildly super-magnetosonic ($M_{\rm f}<8$; Table 5), the internal energy equation does an adequate job of conserving energy, while guaranteeing a positive-definite pressure, which is of greater importance here than strict energy conservation.

Like ZEUS-3D, AZEuS solves the ideal equations of MHD on a fully staggered mesh (zone-centred scalars, face-centred vector components) in an operator split fashion (source terms computed separately from fluxes), using directional splitting for compressive terms (scalar transport, pressure gradient, transport of the $i^{\rm th}$ component of momentum in the $i$-direction), and planar splitting for transverse terms (magnetic induction, transverse Lorentz forces, transport of the $i^{\rm th}$ component of momentum in the $j$-direction, $i\neq j$). AZEuS is upwinded in the entropy and Alfvén waves and relies on a modest amount of Von Neumann & Richtmyer (1950) artificial viscosity to stabilise compressive (fast and slow magnetosonic) waves. Interpolations are performed using the second order, monotonised scheme of van Leer (1977) and, for transverse terms, interpolations are performed implicitly in each plane using the Consistent Method of Characteristics (CMoC; Clarke, 1996b).

As for the AMR module, we have adapted the block-based method of Berger & Colella (1989) and Bell et al. (1994) for the staggered mesh of AZEuS. Significant effort was spent minimising errors caused by waves passing across grid boundaries, which is of particular importance to this work. This includes the development and implementation of third-order interpolation schemes in which mass, momentum, and energy are conserved to machine round-off error. Prolongation of magnetic field is done using a method based on Li & Li (2004), ensuring the validity of the solenoidal condition to machine round-off error regardless of how various 2-D meshes abut, overlap, and overlay each other. Indeed, we find it critical for the solenoidal condition to be valid to machine round-off error even within the boundaries. The interested reader is referred to Ramsey et al. (2012) for details.

All simulations are initialised with nine static, nested grids (including the base grid) with a refinement ratio $\nu=2$. Table 3 gives the extent (in AU) of each of the 2-D grids ($z_{\rm max}$ and $r_{\rm max}$) excluding the boundary regions, along with their resolution, $\Delta$, in each of the $z$- and $r$-directions. Thus, level 1—the coarsest ‘base’ grid—is resolved with $2,\!548\times 160$ zones (including 2 boundary zones at each edge), while each of levels 2–9 are resolved with $640\times 160$ zones.

In addition, smaller grids are added and removed dynamically based on how well the radial gradient of $B_{\varphi}$ is resolved near the symmetry axis. By definition, in axisymmetry both $v_{\varphi}$ and $B_{\varphi}$ should be zero on axis. In our simulations, we find that while $v_{\varphi}$ obliges, $B_{\varphi}$ does not always. Specifically, as a jet propagates, a hot, low-velocity ‘spine’ of strong helical field develops along the symmetry axis. With insufficient resolution, the decline of $B_{\varphi}$ from its maximum value off-axis to zero on-axis is buried within a single zone, creating an ‘inverted profile’ for $B_{\varphi}$, one whose magnitude declines away from the symmetry axis. This generates an axial current density, $J_{z}\propto\partial_{r}(rB_{\varphi})$ of opposite sign to $B_{\varphi}$ (physically, in this situation, $J_{z}$ and $B_{\varphi}$ should have the same sign) which, in turn, exerts a Lorentz force, $F_{r}\propto-J_{z}B_{\varphi}$, directed radially outward (instead of inward). Left unchecked, these unphysical forces occasionally trigger rather dramatic ‘numerical explosions’, sending vast bubbles of hot, rarefied gas expanding into the solution. As interesting as these events are to watch, they are completely numerical in origin and destroy the integrity of the simulation.

We have therefore imposed a ‘Lorentz criterion’ in which a level of refinement is added in any region near the axis where the gradient in $B_{\varphi}$ is insufficiently resolved. Specifically, we require the unitless gradient:

[TABLE]

where $\Delta$ is the zone size, $B=\sqrt{B_{\rm p}^{2}+B_{\varphi}^{2}}$ is the local magnetic field strength, and $N=6$ is the minimum number of zones we require to resolve the radial profile of $B_{\varphi}$. To avoid ‘mesh trashing’ (Khokhlov, 1998), the threshold for removing a grid is $2N$. In practise, we must also guard against ‘frivolous’ inverted $B_{\varphi}$ profiles. Frequently, noisy and dynamically inactive values of $B_{\varphi}$ can create inverted profiles near the symmetry axis and such occurrences should not trigger the insertion of a new grid. In these simulations, we do not go beyond refinement level 9.

3.2 Initial conditions

Young protostellar discs can extend for hundreds of AU, but have inner radii, $r_{\rm i}$, of 3–5 stellar radii, $R_{*}$ (Calvet et al., 2000). For a typical T Tauri star ($M_{*}=0.5\,M_{\odot}$), $R_{*}=2.5R_{\odot}$. Thus, we adopt $r_{\rm i}=0.05$ AU and use this as our length scale. Our finest static grid (level 9; Table 3) resolves $r_{\rm i}$ with eight zones which, based on test simulations of different resolutions, is sufficient for numerical convergence with respect to the physics of the jet launching mechanism.

3.2.1 The atmosphere

The atmosphere is initialised in hydrostatic equilibrium (HSE),

[TABLE]

where $\phi$ is the gravitational potential of $M_{*}$. Since the second term is not a perfect gradient, differencing it directly on a staggered-mesh commits sufficient truncation error to render the atmosphere numerically unstable. Ouyed & Pudritz (1997a, 1999), and continuing with Staff et al. (2010, 2015), address this problem by assuming a strict polytropic equation of state, $p=\kappa\rho^{\gamma}$ where $\kappa$ is constant throughout the grid even across shocks. While this has the advantage of allowing the $\rho\nabla\phi$ term to be written as a perfect gradient which eliminates the numerical truncation error and stabilises the atmosphere, it also replaces energy with entropy as the primary conserved variable. This has the unintended consequences of forbidding the formation of contact/tangential discontinuities and generating isentropic shocks, which we find adversely affects the global solution. Thus, in all of our work, we have retained the adiabatic equation of state ($p\propto\rho^{\gamma}$, but where the proportionality constant remains a function of entropy) to allow the correct entropy jump across shocks and the spontaneous formation of contact discontinuities with their required discontinuities in entropy.

This leaves, however, the numerical instability of the HSE atmosphere unsettled. We address this problem by replacing $\nabla\phi$ in Eq. (11) with the corresponding poloidal gravitational acceleration vector,

[TABLE]

where $\rho_{\rm h}$ and $p_{\rm h}$ are the hydrostatic density and pressure:

[TABLE]

Here, $\rho_{\rm i}$ and $p_{\rm i}$ are the initial density and pressure777The factor of $1/\gamma$ appearing in Eq. (3) of Ramsey & Clarke (2011) is in error. at $r=r_{\rm i}$, and $p\propto\rho^{\gamma}$ ($\gamma=5/3$) is assumed throughout the atmosphere at $t=0$. In this way, differencing Eq. (11) maintains HSE to machine round-off error indefinitely. While it is true that g determined from Eq. (12) is not numerically irrotational (no scalar function $\phi$ exists such that $\textbf{{g}}=-\nabla\phi$ to machine round-off error), this turns out to be an unnecessary requirement on g.

Still, Eq. (12) alone is insufficient to guarantee the numerical integrity of the atmosphere. Regardless of resolution, the singular nature of Eqs. (13) generates sufficient truncation errors at the origin to produce a steady, outwardly directed pressure gradient that drives a supersonic, narrow jet along the symmetry axis, destroying the integrity of the solution. This numerical effect is fixed by replacing the point mass at the origin with a uniform sphere of the same mass and a radius $R_{\rm sph}$, thus modifying the first of Eqs. (13) to:

[TABLE]

If $R_{\rm sph}$ is sufficiently resolved (e.g. four zones), the numerical jet is eliminated. The resulting ‘rounded potential’ is superior to a ‘softened potential’ since the former has no measurable effects beyond $R_{\rm sph}$. Here, we use $R_{\rm sph}=r_{\rm i}$.

As a final comment on the initialisation of pressure, unlike, e.g., Ouyed & Pudritz (1997a, b, 1999), we find no need to attribute a portion of the thermal pressure to Alfvénic turbulent pressure. This was used as a mechanism to reduce what Ouyed & Pudritz felt was an unrealistically high temperature in their outflow. In local simulations, as soon as the confining bow shock leaves the grid, expansion of jet material is free and thus isothermal and the temperature of the fluid does not fall. In our case, the jet expansion is always confined and the gas continuously does $P\,dV$ work to expand thereby reducing the temperature in the inner portions of the jet. Indeed, without additional physics such as radiative heating, the temperature within the inner portions of our jets may be unrealistically low, and certainly we have no need to assign a portion of the thermal pressure to ‘Alfvénic turbulence’.

Similar to Ouyed & Pudritz (1997a), the atmosphere is initialised with a force-free ‘hour-glass’ magnetic field distribution. The toroidal component of the magnetic field, $B_{\varphi}$, is initially zero whereas the poloidal components are set by specifying the toroidal component of the vector potential:

[TABLE]

where $z_{\rm d}$ is the ‘disc thickness’ (which we set to $r_{\rm i}$ for convenience) and where, on AZEuS’ staggered mesh, $A_{\varphi}$ is edge-centred. The $r$- and $z$-components of the magnetic field are then given by:

[TABLE]

which, when differenced, locates the poloidal magnetic field components at the face-centres and ensures $\nabla\cdot\textbf{{B}}=0$ to machine round-off error. Thus, at $(z,r)=(0,r_{\rm i})$ where $B_{\rm p}=B_{\rm i}$,

[TABLE]

where, in terms of the initial plasma beta at $(0,r_{\rm i})$, $\beta_{\rm i}$, used to characterise the simulations (e.g. Table 4),

[TABLE]

We have adopted an ‘hour-glass’ initial magnetic field distribution, in part, because of its simplicity and the availability of a closed analytical form, but we do acknowledge that the true magnetic field distribution in protostellar systems is generally unknown, and indeed will vary from system to system due to environmental differences. Other authors have studied the effect of different field distributions on outflow launching, and find that it significantly affects the collimation of the outflow (e.g. Pudritz et al., 2006; Fendt, 2006). While this aspect is worth exploring, for the current effort, we choose instead to focus only on varying the initial magnetic field strength.

Finally, to ensure the declining density and magnetic field profiles do not fall below observational limits, we add floor values $\rho_{\rm floor}\sim 10^{-6}\rho_{\rm i}$ and $B_{z,{\rm floor}}\sim 10^{-5}B_{\rm i}$ (c.f. Bergin & Tafalla 2007, Vallée 2003) to Eqs. (14) and (16). Thus, the atmosphere attains its asymptotic values by $z\sim 500\,$AU. By imposing HSE and the adiabatic gas law at $t=0$, a floor value on $\rho$ also imposes effective floor values on g, $p$, $T$, etc.

A schematic of the base (coarsest) grid is shown in Fig. 2, where contours of $\rho_{\rm h}$, the ‘hour-glass’ $B_{\rm p}$, and every tenth grid line are plotted.

3.2.2 The accretion disc, and other boundary conditions

The accretion disc, maintained in $z\leq 0$ as a boundary condition, is initially assumed to be in gravito-centrifugal balance (i.e. Keplerian), and to have a force-free magnetic field. Thus, for $z\leq 0$ and $r\geq r_{\rm i}$, $v_{\varphi}=v_{\rm K}=\sqrt{GM_{*}/r}$. We also assume a gentle ‘evaporation speed’, $v_{z}=\zeta v_{\rm K}$ with $\zeta=10^{-3}$ to transfer mass from the disc surface to the atmosphere, preventing any outflow from being ‘starved’ of material. The disc and atmosphere are initially in pressure balance with a density contrast $\eta=\rho_{\rm disc}/\rho_{\rm atm}=100$, while B is initialised from Eqs. (16).

Following Krasnopolsky et al. (1999), $\rho$, $p$, and $v_{z}$ are held to their initial conditions, while $v_{r}$ and $v_{\varphi}$ are allowed to ‘evolve’ in time according to:

[TABLE]

Magnetic boundary conditions are maintained by imposing conditions on the edge-centred induced electric field, $\textbf{{E}}=-\textbf{{v}}\times\textbf{{B}}$:

[TABLE]

where $E_{z}(0)$ is allowed to ‘float’. Since $v_{z}$ is sub-slow, these conditions are formally over-determined and, in principle, $p$ should be allowed to evolve as well. Testing this conjecture, we find that, since $\nabla p$ is $\lesssim 1\%$ of the net Lorentz force at the disc surface, a floating $p$ has only the slightest quantitative effects in the computational domain, yet rather severe consequences within the boundary. Owing to the incomplete dynamics, unphysically high temperatures develop inside the ‘disc’, forcing unnecessarily small time steps on the rest of the simulation. Thus, we fix $p$ to its initial value (matching the initial atmospheric profile) as a numerical convenience throughout the simulations.

Ideally, one would perform a full characteristic analysis at the boundary, setting amplitudes of the outwardly directed waves to zero and using the inwardly directed waves to determine properly upwinded boundary values (e.g. App. A in Del Zanna et al. 2001). Such a capacity has not yet been implemented in AZEuS.

Within the inner radius of the accretion disc ($z\leq 0$ and $r<r_{\rm i}$), we apply reflecting and conducting boundary conditions ($\textbf{{J}}=\frac{1}{4\pi}\nabla\times\textbf{{B}}\neq 0$). Thus, $\rho$, $p$, and v are reflected across $z=0$, and magnetic boundary conditions are set according to $E_{z}(-z)=-E_{z}(z)$, $E_{r}(-z)=E_{r}(z)$, and $E_{\varphi}(-z)=E_{\varphi}(z)$. At $z=0$, $E_{r}$ and $E_{\varphi}$ are evolved using the full MHD equations.

Finally, we use reflecting boundary conditions along the $r=0$ symmetry axis with inversion of $v_{\varphi}$ and $B_{\varphi}$, and outflow conditions along the outermost $z$- and $r$-boundaries, neither of which are ever crossed by anything significant to the simulations.

3.3 Scaling relations

All simulations are performed in units where $\rho_{\rm i}=r_{\rm i}=c_{\rm s,i}=1$, and where $c_{\rm s,i}$ is the sound speed at $r=r_{\rm i}$. Physical units can be restored as follows. First, from Eq. (11) and the adiabatic gas law, one can show that:

[TABLE]

Then, from Eqs. (17), (19), and the ideal gas law ($p=\rho{kT}/\langle{m}\rangle$, where $\langle{m}\rangle$ is $\sim$half a proton mass), the following scaling relations to convert from unitless to physical quantities may be derived:

[TABLE]

where $\gamma=5/3$ and a nominal magnetic field strength of 10 G at $(z,r)=(0,r_{\rm i})$ have been used. Evidently, $T_{\rm i}$ is the temperature at $(z,r)=(0,r_{\rm i})$ and $\tau_{\rm i}$ is the time scale which, for the chosen parameters, is slightly more than a day. As a representative example, simulation E required $\sim\!4.2\times 10^{7}$ time steps on the finest grid ($\sim\!164,\!000$ on the coarsest grid) over a span of $121$ yr to reach the end of the computational domain (4070 AU). The time step in these simulations is typically controlled by the Alfvén speed within several $r_{\rm i}$ of the disc, close to the symmetry axis.

4 Description of the simulations

Table 4 lists the values of $\beta_{\rm i}$ and $B_{\rm i}$ for the eight simulations, A–H, as well as the problem time and jet length at simulation end, assuming the scaling parameters in Sect. 3.3. The simulations were stopped when the tip of the leading bow shock reached the end of the coarsest grid (4,070 AU), or after $t=50,\!000\,t_{i}\simeq 153$ yr, whichever came first. For reference, after 153 yr, the inner and outer edges of the disc (at $r_{\rm i}$ and $4,\!992\,r_{\rm i}$) have undergone $\sim$9,750 and $\sim$0.03 Keplerian orbits, respectively. Time-lapse animations of the simulations described below can be found at http://people.virginia.edu/j̃pr8yu/azeus/proto_jets.html.

4.1 Overview

The simulations listed in Table 4 can be divided into three categories based on the strength of the initial magnetic field: strong (A–D); moderate (E and F); and weak (G and H). Discussion in this subsection on the origins of the jet and its principle morphological features applies mostly to the strong and moderate-field cases, and less so to the weak-field cases. While simulations G and H do generate sustained outflow, they are much more turbulent with far fewer distinctive features than in the stronger field runs.

With this in mind, when any of the simulations begin, a torsion Alfvén wave is launched into the initially stationary atmosphere at $r\geq r_{\rm i}$ by the rotating disc. The wave propagates outward at the local Alfvén speed, $a_{\rm p}$, leaving in its wake atmospheric material rotating in the same sense as the disc, and a toroidal magnetic field, $B_{\varphi}$, is twisted out of the initial $B_{\rm p}$ in a direction opposite to the rotation. Note that $v_{\varphi}$ and $B_{\varphi}$ remain nearly zero within the inner disc radius ($r<r_{\rm i}$) for most of the simulations (with the notable exceptions of the weak-field cases). The torsion wave is a transient feature, borne from the unrealistic initial conditions in which the stationary atmospheric magnetic field threads the rotating disc. That said, it is quickly overcome and absorbed by the leading bow shock of a super-fast jet launched almost immediately from the disc surface, and the torsion wave thus plays a negligible role in the overall appearance of the simulations.

With the passage of the torsion wave, the magnetic field distribution is no longer force-free, and the radial Lorentz force888$\textbf{{J}}\times\textbf{{B}}$ is actually a force density which needs to be integrated over a volume to get an actual force. becomes:

[TABLE]

where $J_{\varphi}$ (as determined from Eqs. 16) remains approximately zero, at least early in the simulations and for the stronger field cases.

The radial profile of $rB_{\varphi}$ is necessarily zero on axis, remains (essentially) zero inside $0<r<r_{\rm i}$ where the torsion wave does not pass, deviates strongly from zero beyond $r=r_{\rm i}$ reaching a global minimum at $r=r_{\rm m}$ ($\sim 1\,{\rm AU}=20r_{\rm i}$), then for the most part returns monotonically and asymptotically to zero as $r\rightarrow\infty$ (e.g., bottom middle panel of Fig. 4). This profile is directly related to the three distinct regions of the jet that develop in the simulations, as depicted in the top panel of Fig. 3 and described below.

First, inside $0<r<r_{\rm i}$, where no material is driven onto the grid from the $z=0$ boundary and where the axial field remains nearly force-free, a narrow, cold and relatively quiescent magnetic ‘spine’ develops along the symmetry axis. Its integrity is maintained throughout the simulations in all but the weakest-field cases. In 3-D, however, it is unlikely this feature could survive given the higher mode instabilities that tend to disrupt the axisymmetry of the jet (e.g. Hardee & Clarke 1995).

Second, within $r_{\rm i}<r<r_{\rm m}$, $B_{\varphi}$ and $\partial_{r}(rB_{\varphi})$ have the same sign and $F_{r}$ in Eq. (25) is directed inward, compressing material toward the axis. This sets up what becomes a dense, hot, strongly magnetised, super-fast ‘jet core’ threaded by the spine.

Third, for $r>r_{\rm m}$, $B_{\varphi}$ and $\partial_{r}(rB_{\varphi})$ have opposite sign, and $F_{r}$ in Eq. (25) is directed outward. This drives material away from the axis and opens up a wide cavity terminating radially at the tangential discontinuity (TD), located at $r=r_{\rm TD}$ which, in Fig. 4), is located at $r\sim 100$ AU. This cavity is cold, highly magnetised and, at least initially, relatively evacuated, thus allowing material driven from the disc easy passage. This outflow forms an exceedingly cold, strongly magnetised, trans-fast, and dense (but less dense than the core), ‘sheath’ that surrounds the jet core.

The torsion wave sets material spinning, and material near the disc surface with its frozen-in $B_{\rm p}$ experiences a centrifugal acceleration outward. If the angle between the disc and a magnetic field line, $\theta_{0}$ (Fig. 1), falls below a critical angle $\theta_{\rm c}=60^{\circ}$, the magnetocentrifugal ‘bead-on-a-wire’ mechanism (BWM; Henriksen & Rayburn 1971; Blandford & Payne 1982) is triggered and accelerates material away from the disc. From Eqs. (15) and (16), one can show that initially at the disc surface, $\theta_{0}=\tan^{-1}B_{z}/B_{r}<\theta_{\rm c}$ for $r\geq\sqrt{3}r_{\rm i}$, and it is only in the region $r_{\rm i}<r<\sqrt{3}r_{\rm i}$ where $\theta_{0}$ must be reduced before outflow can begin. As discussed in Sect. 5.6, this turns out to be key in the persistent generation of ‘knots’ (plasmoids) observed in one of our simulations.

The BWM requires a ‘rigid’ $\textbf{{B}}_{\rm p}$, which can occur only if $\beta_{\rm p}=8\pi p/B_{\rm p}^{2}\lesssim 1$. If $\beta_{\rm p}\gg 1$, $B_{\rm p}$ cannot provide the magnetic tension necessary to act as the ‘wire’ to guide the ‘beads’ of plasma, and a different mechanism to drive the jet must be invoked. As discussed further in Sect. 5.5, even a weak $\textbf{{B}}_{\rm p}$ in a persistently rotating environment will generate a dynamically important $B_{\varphi}$ ($\beta_{\varphi}=8\pi p/B^{2}_{\varphi}\sim 1$), and this can provide sufficient magnetic pressure to accelerate an outflow. This is sometimes referred to as the ‘magnetic tower mechanism’ (MTM; Lynden-Bell 1996).

Whether driven magneto-centrifugally and/or via a sustained magnetic pressure, outflow is strongest near $r_{\rm i}$ where the rotation is the most rapid and the magnetic field strongest. While the driving force gets progressively weaker with increasing $r$, the density in the disc and sheath also fall, and material can still be accelerated to super-fast speeds. Beyond $r_{\rm TD}$ (outer limit of the sheath), the weak Lorentz force cannot accelerate the much denser ambient material effectively. In this way, a fairly distinct boundary is visible between fast, outwardly moving material—the actual ‘jet’ consisting of a core and sheath—and the more slowly moving, shocked and much denser atmosphere behind the jet bow shock.

Thus, as annotated in the top panel of Fig. 3, the general outflow established in our simulations consists of three outwardly moving components. First, the super-fast, dense, hot, strongly and helically-magnetised narrow jet core—what we suggest corresponds to the main portion of the observable jet—is threaded by a narrow, quiescent spine of strong axial magnetic field. Second, the jet core is surrounded by a somewhat slower-moving trans-fast, under-dense, cold, highly magnetised wide sheath, which may, in part, be observable as a ‘cavity’ of emission surrounding the jet. As noted below, this region is subject to strong reflection shocks triggered by the Kelvin-Helmholtz unstable TD, that can heat sheath material significantly. As a result, some of the sheath may contribute to the observable outflow. We note in passing that the development of an under-dense sheath could account for the great distances to which jets can propagate stably in 3-D (Hardee et al. 1997). Third, between $r_{\rm TD}$ and the bow shock, the sheath is surrounded by a trans-slow, hot, weakly-magnetised, shocked ambient medium, whose density is less than that of the jet core, but greater than that of the undisturbed ambient medium. This may correspond to the ‘second wind’, first described by Stocke et al. (1988) and currently interpreted as ‘molecular winds’ (e.g., Frank et al. 2014), whose forward motion is a result of entrainment by the leading bow shock driven by the advancing jet rather than magnetic stresses near the disc.

4.2 Strong-field simulations: A–D

Differences among the four simulations in the strong-field category are relatively minor, and we use as an exemplar the final epoch of simulation A ($\beta_{\rm i}=0.1$, $B_{\rm i}=200\,$G) shown in Fig. 3, in which colour contours of temperature are superposed with magnetic field lines (white contours) and the slow surface (where $v_{\rm p}=a_{\rm p,s}$, the poloidal slow speed; black contours). Here, the magnetic field is strong enough to enforce $\textbf{{v}}_{\rm p}\parallel\textbf{{B}}_{\rm p}$ virtually everywhere inside the TD, resulting in what appears to be a largely self-similar, steady state solution. Indeed, the bottom panel in Fig. 3 could be taken from any time within the final $\sim$60% of the run and the relative simplicity of this run makes it prototypical of the description given in Sect. 4.1.

The bottom panel of Fig. 3 bears a resemblance to the early, local simulations of Uchida & Shibata (1985). Near the disc, outflow is driven almost exclusively by the BWM, with the MTM gradually contributing additional thrust as the poloidal field becomes wound up even before the Alfvén point (Sect. 5.5). Outflow is robust, long-lived, steady, and there is no reason to believe it would ever cease so long as the disc continues to provide mass. Magneto-centrifugal wind launching theory (e.g. Blandford & Payne, 1982) predicts that acceleration of the flow ceases beyond the Alfvén point, while some local simulations (e.g. Pudritz et al. 2006) show acceleration of the outflow continues until the fast point, with steady flow thereafter. However, we find that because of the MTM, a gentle acceleration persists well beyond the fast point, accelerating the advance speed of the jet to $\sim 420\,{\rm km\,s}^{-1}$ for simulation A ($v_{\rm jet}$ in Table 5), rendering the jet essentially ballistic (since $v_{\rm p,max}\sim v_{\rm jet}$). Further, the sound speed in the asymptotic ambient medium is,

[TABLE]

using the asymptotic density in Sect. 3.2.1. Thus, the external sonic Mach number of the jet in simulation A is $M_{\rm ext}\sim 540$. Of course, the jet core is much hotter than the ambient medium (by a factor of $\sim 10^{4}$ and as high as $10^{6}$ K), and the internal sonic Mach number is therefore considerably less; $M_{\rm int}\sim 22$. The internal fast magnetosonic Mach number, that which actually governs the ‘supersonic character’ of the jet, is a rather modest 2.1. See Table 5 for comparative values for the other simulations.

Hugging the axis, but not particularly apparent in even the lower panel of Fig. 3, is the relatively quiescent, cold spine with negligible $B_{\varphi}$, strong $B_{\rm p}$ and a radius which remains nearly constant ($\sim\!2r_{\rm i}$) (but not resolved by levels $<$ 6). This spine threads the dense, hot, highly magnetised, super-fast jet core which, as seen in the bottom panel of Fig. 3, reaches a radius of $\sim 0.3\,$AU (outer extent of the hot region) at $z=4\,$AU from the origin. Fig. 4 shows a radial slice of numerous variables across the computational domain at $z=1,\!200\,$AU. In these slices, the jet core is demarcated by the drop in density from $10^{8}$ to $10^{6}$ cm*-3* in the top left panel, and the global minimum of $B_{\varphi}$ in the bottom middle panel, and thus has reached a radius of about 1 AU. At $z\sim$ 3,500 AU, the jet core reaches its maximum radius of a few AU.

Surrounding the jet core is the magnetic sheath, and the boundary between the sheath and the hotter, denser shocked ambient medium is the tangential discontinuity, best visualised in the top panel of Fig. 3. The distinguishing feature of a TD (as opposed to an ordinary contact discontinuity) is the lack of a normal component of the magnetic field; B is everywhere parallel to the feature we have identified in Fig. 3 as the TD. The TD is also apparent in Fig. 4 at $r\sim 100\,$AU as a sudden drop in $v_{\rm p}$ (left centre panel), a jump in temperature (middle top panel), the end of a gradual drop in $v_{\varphi}$ (middle centre panel), a sharp rise in the plasma-$\beta$ (right centre panel), and a drop from super to sub-fast-magnetosonic speed (right bottom panel) all while the total pressure (thermal + magnetic) remains more or less continuous (top right panel).

Dynamically, the most important of these characteristics is the drop in $v_{\rm p}$, which means the TD is a shear layer and subject to the Kelvin-Helmholtz instability, manifest in the top panel of Fig. 3 as gentle undulations along its length. Between $500\,{\rm AU}\lesssim z\lesssim 900\,{\rm AU}$, these undulations have a wavelength of about $100\,$AU and gradually grow in amplitude. By $z\sim 900\,$AU, the severity of the undulation triggers a fairly strong reflection (criss-cross) shock in the sheath that thermalises enough of the kinetic energy to warm the sheath significantly (from several to a few hundred kelvins), and redirect flow parallel to the jet axis. This ceases the expansion of the jet sheath, whose outer radius at this epoch of simulation A is about 100 AU. After $z\sim 900$ AU, the K-H undulations continue with lesser amplitude and a wavelength of 250–300 AU, and trigger a series of gentler reflection shocks and rarefaction fans similar to those described for a hydrodynamical jet by Norman et al. (1982). These are all visible as discontinuities and gradations in temperature in Fig. 3. At $z\sim 2,\!400$ AU, one final, relatively strong reflection shock is triggered, after which the flow remains rather laminar and featureless, bearing a strong resemblance to the ‘nose-cone’ described by Clarke et al. (1986) for a propagating jet dominated by a toroidal magnetic field. We note in passing that the strong poloidal magnetic field, which provides some stability against the $m=0$ ‘pinch mode’ apparent in these axisymmetric simulations, would also provide some stability against higher mode instabilities in 3-D (e.g. Hardee et al. 1997), and possibly enough to preserve the nearly axisymmetric appearance of jets such as HH 34 (e.g. Devine et al., 1997).

The shocked ambient medium (sometimes referred to as the ‘second wind’) lies between the TD and the bow shock (top panel of Fig. 3) and is characterised as a warm, trans-slow, dense, relatively weakly magnetised medium. It owes its forward motion entirely to entrainment by the leading bow shock, and has virtually zero rotation and toroidal field. The sub-slow ‘islands’ (closed black contours in the top panel of Fig. 3) indicate the forward motion is trans-slow ($\sim\!10\,{\rm km\,s}^{-1}$). Where the flow is sub-slow, streamlines diverge (slightly) from the jet axis, and where the flow is super-slow, streamlines converge. These transitions are precisely coupled to the K-H undulations described above. The temperature ranges from $\sim 30,\!000\,$K just above the TD, to a few hundred K just inside the bow shock (top centre panel of Fig. 4), and thus permits the survival of molecules such as CO, which are frequently used to observe entrained outflow material (e.g. Zhang et al., 2016). Finally, the plasma-$\beta$ jumps 4.5 orders of magnitude across the TD (thermal pressure suffers a sudden increase while magnetic pressure undergoes a commensurate decrease to maintain a near-continuous total pressure), and the shocked ambient medium is dominated by thermal pressure ($\beta\sim 10$) just above the TD. The plasma-$\beta$ then drops continuously between $r_{\rm TD}$ and the bow shock where $\beta<0.1$, and the flow is once again magnetically dominated, but not to the extent (by a factor of a few) observed in the jet core and sheath (right centre panel of Fig. 4).

Last, what remains of the undisturbed ambient medium is visible above the bow shock in the top panel of Fig. 3, where the asymptotic levels for $\rho$ and $B_{z}$ (Sect. 3.2.1) are reached (to within 10%) by $z\sim 2$,300 AU.

While for convenience, the simulations are identified by the single value $\beta_{\rm i}$, this does not represent the average magnetic field strength in the resulting outflow. Figure 5 shows the evolution of the average plasma beta, defined as,

[TABLE]

where quantities in angle brackets on the right hand side are volume averages over the outflow identified as regions where $M>5$ and $v_{\varphi}>10^{-3}$. The $M=5$ contour was chosen as it was found to hug tightly inside the TD for all jets thus eliminating the shocked ambient medium and less organised outflow near the accretion disc from the average. The limit on rotational speed eliminated material from the jet spine. Mass-weighted averages were also computed (not shown), giving qualitatively similar loci as in the figure, with values for $\langle\beta_{\rm tot}\rangle$ consistently about 0.2 lower.

What is striking about Fig. 5 is that all loci (with the exception of simulation H) seem to be converging on $\langle\beta_{\rm tot}\rangle\rightarrow 0.2$–0.4, characteristic of a magnetically dominated outflow. Even simulation H, which of all the simulations had the most difficult time organising itself into an outflow, seems to be converging toward magnetic equipartition ($\langle\beta_{\rm tot}\rangle\rightarrow 1$). From these observations, we speculate that however weak or strong the initial magnetic field in the corona may be that launches the outflow, the outflow itself ends up being magnetically dominated, with an average plasma beta asymptoting to $\lesssim 1$. If true, an immediate consequence is that the magnetic properties of an observed jet may not be useful in determining what the magnetic environment may be near the protostar.

The nature of the magnetic field within the jet also evolves with time. At first, most magnetic field in the outflow is poloidal, reflecting the initial conditions. As the jet advances and rotates, $B_{\rm p}$ is twisted to produce a significant $B_{\varphi}$ which, by the end of the simulation, accounts for $>90\%$ of the transported magnetic energy density even for simulation A with the strongest initial poloidal field. This is why we find in our simulations that the MTM is the dominant acceleration mechanism even in simulation A (Sect. 5.5), and why our jets continue to accelerate well beyond the fast point.

The higher-$\beta_{\rm i}$ runs (B–D) are qualitatively identical to simulation A, but with significant quantitative differences. As $\beta_{\rm i}$ increases, jet speeds decrease as do the sonic Mach numbers, although the fast magnetosonic Mach number actually increases monotonically from $\sim 2.1$ for simulation A to $\sim 4.0$ for simulation D. These trends continue into the moderate and weak field runs (Table 5). The temperature, density, and strength of reflection shocks within the outflow also decrease with increasing $\beta_{\rm i}$, while the time taken for a strong outflow to be organised increases. This trend is apparent by simulation F (where outflow doesn’t really begin until $t\sim 0.2\,$yr), and continues on through simulation H which, in many respects, shows signs of being a ‘frustrated jet’ (Sect. 4.4). As part of this trend, the inner-most regions of the jet become less steady for higher $\beta_{\rm i}$ to the point where periodic knots start to form. These are sporadic for simulation E but, by simulation F, the knots are steady, long-lived, and dominate the inner-jet structure (Sect. 4.3 and 5.6).

Finally, near the bottom of the third panel in Fig. 3 is a prominent ‘streak’ originating from the level 8 grid boundary at $(z,r)\sim(8.0,1.3)$, and stretching to $(z,r)\sim(23,2)$ in the middle of the level 6 grid. These transient features, more prominent in the stronger field simulations than the weaker ones, are entirely numerical in origin, and triggered where the Alfvén surface intersects a grid-boundary. At such points, and then only rarely, a truncation error occurs in the momentum interpolation which results in a sudden pinch and subsequent local spike in internal energy that is advected downwind with the flow, resulting in the streak seen. We are not entirely certain why these streaks occur when they do, though they are reminiscent of the ‘magnetic field explosions’ suffered in earlier versions of ZEUS (Clarke, 1996b). We also know that if one interpolates on velocity in the grid boundaries rather than momentum density (and violate conservation of momentum between grids), these features disappear.

Despite their ominous appearance (in the streak in Fig 3, $T\sim$ a couple thousand kelvins within a cold sheath of a few kelvins), the streaks are not apparent in the other MHD variables, and this region is so magnetically dominated that such a thermal pressure anomaly has negligible dynamical consequences. We also note that the streaks are entirely absent in the weaker field runs (F–H) and, at this stage, we regard them as cosmetic. Still, understanding their origin remains an area of current investigation.

4.3 Medium-field simulations (E, F)

As the initial magnetic field strength ($B_{\rm i}$) is weakened, all measures of outflow speed decrease (Table 5) and the outflow itself becomes less steady. While jets on the observational scale look largely the same (e.g., compare top panels of Figs. 3 and 6), the smaller scale structures in the bottom panels are strikingly different.

Figure 6 shows simulation F at $t=153$ yr, having reached a length of $\sim 3,\!800$ AU999A similar simulation taken to 100 yr was described in Ramsey & Clarke (2011).. At the jet tip, $v_{\rm jet}\sim 116$ km s*-1*, $M_{\rm ext}\sim 150$, $M_{\rm int}\sim 10$, and $M_{\rm f,int}\sim 6$ (Table 5). Roughly 99% of the magnetic energy density within the magnetic sheath is in $B_{\varphi}$ from which we conclude the MTM is the dominant driver along most of the jet. Indeed, for simulation F (and even more so for weaker $B_{\rm i}$), the Alfvén surface lies close to the disk surface ($\lesssim 30$ AU) for material interior to the TD, and the BWM is ineffective at any significant height above the disc.

We have chosen simulation F as the exemplar for the medium-field runs because of the propensity and regularity of knots, seen in the lower panel of Fig. 6. Knots, which are virtually absent in simulations A–D, are present in simulation E, but not to the degree seen in simulation F. While the production of knots in simulation F is often ‘steady’, they do not represent a steady state; none of the ‘constants’ in Eqs. 1–4 are constant along field lines passing through them. In this simulation and where they are present in simulation E, only the portion of the sheath devoid of knots and away from the TD is in a quasi-steady state, at least as measured by the constancy of the WD constants (Eqs. 1–4; Figure 12).

The knots are launched from $(0,r_{\rm i})<(z,r)<(r_{\rm i},2r_{\rm i})$ where gas is both dense and hot. Knots are generated nearly from the beginning and, after a few short periods of intermittency and variability, become steady after $t\sim 9$ yr with a period ${\cal T}_{\rm knot}\sim 0.026$ yr and a wavelength $\lambda_{\rm knot}\sim 0.25$ AU.

Once launched, the knots follow a nearly axial trajectory through the sheath, never venturing further from the axis than a few AU. As can be seen in the lower panels of Fig. 6, they are best described as hot ‘towers’ (tori in 3-D) of plasma following the local magnetic field lines; truly ‘beads (tori) on a wire’. As they move downstream, they lengthen, merge and, at $z\sim 10$ AU above the disc, coalesce into the continuous hot jet core of radius $\sim 2\,$AU, which continues to expand gradually to several AU towards the head of the jet. Never drifting more than a few AU from the axis and losing their identity long before reaching observational scales, we rule these features out as precursors of HH objects.

Along their length, the knots exhibit one, two, and sometimes three extrema in $T$, likely a result of uneven magnetic confinement. The reader is encouraged to examine animations of this simulation (available on-line101010http://people.virginia.edu/j̃pr8yu/azeus/proto_jets.html), where the propagation of the knots is seen to be extremely dynamic.

Figure 7 is an axial slice of several variables at $r=0.4$ AU ($8r_{\rm i}$) at the final epoch of simulation F. This plot, however, could have been taken at any time $t>10$ yr, so regular are the knots. Relative to the cold sheath material, the knots are 10–20 times denser, $\sim\!10^{3}$ times hotter, and thus have $10^{4}$ times the thermal pressure of their immediate surroundings! Were it not for the fact that $\beta$ is still $<1$ within these ‘plasmoids’, they would explode into the ambient gas upon creation. As it is, the knots are magnetically confined (the total pressure in Fig. 7 is nearly continuous) by strong poloidal flux loops that effectively contain them as they propagate intact $\sim 50$ knot radii along the jet axis. Other differences in the knot material compared to their immediate surroundings include: $v_{\rm p}$ is 10–20% higher, $v_{\varphi}$ is $\sim 50\%$ lower, $B_{\rm p}$ is $\sim 50$% lower, and $|B_{\varphi}|$ is $\sim 20$% lower.

The fact that the knots move slightly faster than the surrounding sheath material is interesting. In this simulation, the difference of a few km/s is still super-slow, and thus the knots excite leading slow shocks—much like the jet itself excites a bow shock in the ambient medium—explaining their discontinuous leading edge (Fig. 7) and minimal diffusion.

We resume discussion on the knots in Sect. 5.6, where the physics of knot generation is addressed.

4.4 Weak-field simulations (G, H)

Figure 8 shows simulation G at $t=153$ yr having reached a length of $\sim 2,\!800$ AU. Evidently, $\beta_{\rm i}$ has passed a critical value as the qualitative appearance of the inner jet is dramatically different even from simulation F with $\beta_{\rm i}$ just a factor of four higher. The jet speeds continue to diminish with $B_{\rm i}$ ($v_{\rm jet}\sim 86$ km s*-1*, $M_{\rm ext}\sim 110$, $M_{\rm int}\sim 12$) while the fast magnetosonic Mach number increases ($M_{\rm f,int}\sim 7.5$; Table 5). However, the most striking difference is the nearly complete replacement of organised knots with ‘turbulent’ outflow and displacement of the initial poloidal field lines within the jet core and sheath (vestiges of the knots—the ‘bases’ identified in Sect. 4.3—are, however, still apparent near the jet axis). This trend is even more evident in simulation H (not shown), where the turbulent nature of the outflow extends right to the TD, all but eliminating it as a discernible feature. The spine—with its nearly straight axial field—remains prominent and hot next to the axis, and while much of the outflow is still launched from within $r_{\rm i}<r<2r_{\rm i}$, the organisation and steadiness observed in the strong and medium-field simulations is not apparent in the weak-field cases.

Still, on the observational scale (top panel of Fig. 8), the jet looks much the same as the stronger-field simulations. Regardless of how weak the magnetic field is at the base of the jet, the gravitational and rotational effects organise to amplify the magnetic field enough to generate a large-scale, supersonic outflow transporting significant magnetic field energy. This is true even for simulation H with $\beta_{\rm i}=640$, what ought to be considered an essentially hydrodynamical environment.

As will be shown in Sect. 5.5, the weak $B_{\rm p}$ is continuously wound up into $B_{\varphi}$ until such time as $a_{\varphi}\sim c_{\rm s}$. At this point, the outwardly-directed gradient in toroidal magnetic energy density is comparable to thermal pressure gradients and gravitational forces, and material can be launched and accelerated outward. This is true even for the weakest initial magnetic field (simulation H) which, with very little poloidal field to contribute to jet confinement, gives rise to a turbulent jet that entirely fills the magnetic sheath.

As seen in Fig. 5, once the outflow becomes organised in simulation G ($t\sim 1$ yr), $\langle\beta_{\rm tot}\rangle$ quickly falls below $\beta_{\rm i}$ to a surprisingly low value of $\sim 0.02$, then rises steadily to $\sim 0.6$ after which it declines asymptotically to $\sim 0.4$. Even in the weaker-field simulation H, when outflow first begins, its value of $\langle\beta_{\rm tot}\rangle$ is well under unity, and then rises—rather sporadically—toward unity as the run progresses.

The spikes in $\langle\beta_{\rm tot}\rangle$ seen for simulation H attest to the marginality with which steady outflow is established in this very weak magnetic environment. Speculating that $\beta_{\rm i}\sim$ several hundred may represent an ‘end-of-the-line’ for a successful jet launch, we ran another simulation with $\beta_{\rm i}=2,\!560$ to see if the MTM would prevail in such a weak magnetic environment. Sure enough, even in this extreme limit, a jet is launched, albeit even more turbulent than simulation H and requiring even more time to organise itself into an outflow.

Of course, 2-D axisymmetry provides an artificially favourable geometry for the MTM in which a weak poloidal field can do nothing but wind up under the relentless rotation of the fluid. Furthermore, axisymmetry is immune to all but the $m=0$ mode of the Kelvin-Helmholtz instability, which not only allows the toroidal field to develop unimpeded, but provides a perfectly rigid axis along which flow can be directed.

In 3-D, the physics is not so accommodating. With virtually no poloidal field to stabilise the flow, mixing of the shocked ambient medium with the sheath occurs along the K-H unstable TD while the turbulent flow encourages mixing of the jet core with the sheath. With less contrast in density and temperature between the core and sheath, the core is at risk of higher order MHD instabilities that may cause it to break up into tangled filaments after propagating only several jet radii (Hardee et al., 1997). Further, an advancing, toroidally confined, magnetic column is unstable to the kink instability (e.g. Jackson 1975), exacerbating its ability to propagate to observable length scales. We speculate, then, that in 3-D and in cases where $\beta_{\rm i}\gtrsim$ a few hundred, a ‘frustrated jet’ could result in which outflow never manages to organise itself to advance very far from the disc. Other than indicating when 2-D turbulence fills the jet sheath ($160<\beta_{\rm i}<640$), our present simulations are unable to determine what this critical value for $\beta_{\rm i}$ may be.

5 Analysis

5.1 Morphology

To the best of our knowledge, these are the first jet simulations to link the physics of magnetically driven outflows from Keplerian discs, as first reported by Uchida & Shibata (1985), and the physics of magnetically collimated supersonic outflows, as first reported by Clarke et al. (1986) (hereafter, CNB).

Regardless of $\beta_{\rm i}$, all simulations form a magnetically dominated outflow ($\langle\beta_{\rm tot}\rangle\rightarrow 0.2$–0.4, with the exception of simulation H) which come to resemble CNB jets (who, coincidentally, used $\beta_{\varphi}=0.2$) at observational scales. Thus, they all possess a hot, supersonic core with an advance speed of 80–420 km s*-1* and a number of oblique shocks triggered along their length. The core terminates with a ‘jet shock’ (occasionally a Mach stem), with most of the post shock jet material collected in front of the jet shock forming a ‘nose-cone’.

As discussed in CNB, the nose-cone consists of trans-fast material—hotter than in the jet core—which, without the confining magnetic field, would form the ‘back-flowing cocoon’ reported by Norman et al. (1982). Since the nose-cone is more pointed and denser than the hydrodynamical cocoon created in the absence of $B_{\varphi}$, a magnetically confined jet pushes into the ambient medium more ballistically than a hydrodynamical jet, consistent with the present simulations at the largest scale.

Surrounding the jet core and as reported by CNB, a largely evacuated ‘cocoon’ filled with cold, highly magnetised rarefied material with a strong toroidal magnetic field confines and stabilises the jet core. In the present simulations, we refer to this feature as the ‘sheath’, prominent in all simulations except H, the most weakly magnetised jet in our sample.

Surrounding the cocoon (sheath) and the leading nose-cone is a bow shock excited in the quiescent ambient medium by the passage of the supersonic jet. The shocked ambient material is accelerated forward ($\lesssim$​ 10 km s*-1*) forming what we identify as the ‘second wind’, consisting of material not launched from the disc but entrained in situ by the bow shock from the ambient atmosphere. As such, it is warm (300–30,000 K), weakly magnetised, has virtually zero toroidal velocity and field, with a magnetic field strength directly related to the magnetic conditions in the primordial atmosphere and thus proportional to $B_{\rm i}$.

Because the ambient atmosphere is so cold, the Mach number of the jet relative to the external atmosphere is very high, even for simulation H where $M_{\rm ext}\sim 100$. Thus, the bow shock is very narrow and the radius of the jet depends very weakly on $B_{\rm i}$ making the jet radius a poor indicator of the magnetic conditions near the disc surface. Indeed, the morphology of the observable jet as a whole—being so similar for all values of $B_{\rm i}$—cannot be used as an indicator of this most elusive of physical quantities.

It is worthwhile noting that our choice of magnetic field distribution (Sect. 3.2.1) plays a role in determining the outflow morphology. A steeper radial field distribution ($B(z=0,r)\propto r^{-1}$ in our case) is known to result in less collimated jets (Pudritz et al., 2006; Fendt, 2006), and could lead to jets with larger opening angles and larger radii. We also recognise that, in 3-D, with the availability of additional K-H modes, some of the jet kinetic energy would be converted to turbulent or thermal energy, and we subsequently expect the jet bow shock to not only be blunter and wider, but also propagate more slowly. Finally, on very large scales ($\sim$0.1 pc; Frank et al. 1999), jets are expected to be weakly ionised and ambipolar diffusion will be important, resulting in less efficient magnetic confinement and an additional source of heating (e.g. Pinto et al., 2008; Panoglou et al., 2012). Like the magnetic field distribution, any of these effects could produce a jet that is less ‘knife’-like than in the current simulations.

5.2 Qualitative trends

Certain quantities, such as the advance speed of the jet ($v_{\rm jet}$), reach an asymptotic limit early in the simulation (e.g. after $\sim 10$ yr), making comparisons among the simulations easy and straight-forward. Other quantities such as the jet radius ($r_{\rm jet}$) continue to grow with time, and thus one must decide how such quantities are to be compared. In particular, one could make comparisons at the same chronological time or at the same dynamical time, the latter defined as when a jet has reached a specified length.

Fig. 9 shows the poloidal velocity, $v_{\rm p}$, and the fast magnetosonic surface (white contours) for all eight simulations at the same chronological time, $t\sim 47$ yr, when simulation A reaches the end of the domain. This figure illustrates how the jet length and radius at a fixed time depend rather strongly upon $B_{\rm i}$.

However, observationally, one can be more certain of comparing jets of the same length than of the same age, and thus Fig. 10, where the jets are shown at the same dynamical time $l=2$,400 AU, is more practical. In this case, there is little to distinguish the eight simulations geometrically. At a given length, the radii of the sheath (indicated, for the most part, by the TD) and jet core increase slightly with decreasing $B_{\rm i}$, but the dependence is so weak as to make these impractical observational comparators for inferring the magnetic field properties at the base of a jet.

5.3 Weber-Davis constants revisited

We now return to the WD constants defined in Sect. 2 (Eqs. 1–4), as well as Eq. (8) for the fast speed along a steady-state field line. Since most of our global simulations do not reach a true steady state, these expressions will not be valid along all field lines. However, as seen in Figs. 3, 6 and, to a much lesser extent, Fig. 8, the first few hundred AU along field lines anchored in the disc at $0.5\lesssim r_{0}\lesssim 10$ AU are smooth and appear to reach a quasi-steady state. Within these regions the steady state constants can be examined, both to provide an in situ check on our numerical methods, and for what they can tell us about the flow.

For convenience, we have chosen the field line anchored in the disc at $r_{0}=1$ AU to perform this analysis as it exemplifies steady state behaviour for all but the weak-field simulations. However, we emphasise that these results apply equally for portions of any field line anchored between 0.5 AU and 10 AU, with some simulations showing better steady state behaviour than others. For the 1 AU field line, Fig. 11 shows plots of various speeds as a function of distance along the field line, $s$, including $v_{\rm p}$, $v_{\varphi}$, $a_{\rm p}$, $a_{\varphi}$, and $a_{\rm f}$, all defined in Sect. 2. The dashed blue line shows $r\omega_{0}/2$ (where $\omega_{0}=v_{\rm K,0}/r_{0}$ is the angular speed at the anchor point of the field line), an important quantity in understanding the BWM (Sect. 5.5).

The Alfvén point ($v_{\rm p,A}$, where $v_{\rm p}$ and $a_{\rm p}$ intersect in Fig. 11) and fast point ($v_{\rm p,f}$, where $v_{\rm p}$ and $a_{\rm f}$ intersect) are indicated in the figure with triangles and squares, respectively. The quantities $v_{\rm p,A}$, $v_{\rm p,f}$, as well as the asymptotic poloidal speed attained along the 1 AU field line, $v_{\rm p,max}$, are included in Table 5. For the 1 AU field line and most any other steady state field lines we examine, acceleration of jet material continues well beyond the fast point (via the MTM), less so for weaker $B_{\rm i}$.

Figure 12 shows profiles of the fractional differences of each WD constant (Eqs. 1–4) relative to their expected values (Eqs. 5 and 6) in each of simulations A–F. As can be seen, the WD constants remain—for the most part—constant to within 3% along the 1 AU field line. Notable exceptions include simulation F ($\pm 5\%$), which is really a transitional simulation between those which exhibit strong steady state regions (A–E), and those that do not (G and H), and the specific energy constant, $\varepsilon$, for simulations A and B, which is dominated by the difference between two large and nearly equal numbers ($v_{\varphi}^{2}/2$ and $rv_{\varphi}\Omega$ in Eq. 4). This dominance decreases with $B_{\rm i}$ and, as such, $\varepsilon$ is constant to within 3% for simulations D–F. Simulation D also exhibits a transient feature in $\varepsilon$ at $\sim$9 AU from when the field line passes through one of the streaks mentioned in Sect. 4.2; this feature is not visible in the other steady state constants and is only barely visible in Fig. 11. We note that the field lines embedded in quasi-steady state regions pass through several nested grids, and take this as evidence that our adaptive mesh and MHD algorithms maintain conserved quantities satisfactorily.

Table 6 shows the expected values of $v_{\rm p,f}$ along the 1 AU field line (Eq. 8) for simulations A–G compared to the values shown in Fig. 11. Notwithstanding simulation G, all values agree to within 1.2%, once again indicating portions of the jets do attain a quasi-steady state, and that the code is able to recognise and maintain these regions. In as much as there is a ‘numerical test’ for these simulations, this would be it.

As can be seen from Figure 13, the poloidal speed at the fast point, $v_{\rm p,f}$, along the 1 AU field line follows a rather tight ‘2/3 power law’ with $B_{\rm i}$, contrary to the prediction made by equation (9) and more in line with Eq. (10). Indeed, the fast point along all field lines passing through a quasi-steady state region show equally tight power laws $v_{\rm p,f}\sim B_{\rm i}^{\alpha}$, with $\alpha=0.67\pm 0.03$ (Table 5). This is one of the most robust power-law relationships gleaned from the simulations, and we regard this result as firm.

The discrepancy between the numerically determined power-law index and that predicted by Eq. (9) can only be caused by the factor,

[TABLE]

in Eq. (8), since it was assuming it to be independent of $B_{\rm i}$ that led to the ‘1/2 power law’ prediction in Eq. (9). Evidently, $f\sim B_{\rm i}^{1/6}$, a weak but significant dependence that accounts for the difference between the measured and predicted power law indices for $v_{\rm p,f}$. What was not anticipated in Sect. 2 was the degree to which the fast point is pushed away from the disc for higher $B_{\rm i}$ ($z_{\rm f}$ in Table 5). Along the 1 AU field line, $z_{\rm f}\sim 300\,$AU for simulation A, whereas $z_{\rm f}\sim 7\,$AU for simulation F. Since density is highly dependent upon $z$, this is where the dependence of $f$ upon $B_{\rm i}$ arises.

We continue this discussion in the next subsection.

5.4 Dependence of velocity on $B_{\rm i}$

While both $v_{\rm p,f}$ and $v_{\rm p,A}$ follow a power law in $B_{\rm i}$ along field lines in quasi-steady state, they do so with different indices, and the tightness of fit is somewhat stronger for $v_{\rm p,f}$ (Fig. 13; Table 5). This is true for other field lines in quasi-steady state as well, with similar power law indices measured.

In addition to speeds specific to the 1 AU field line, Table 5 includes various observationally-accessible speeds for each jet, including the mass-weighted averages (exclusive of the second wind) of the rotational velocity, $\langle v_{\phi}\rangle$, the axial speed along the last 500 AU of the jet core, $\langle v_{z}\rangle$, and the advance speed of the jet tip into the quiescent atmosphere, $v_{\rm jet}$. Each velocity varies with $B_{\rm i}$ reasonably well as a power law (with power-law index $\alpha$ given in Table 5), indicating that the observable kinematics of the jet is highly and simply dependent upon the initial magnetic field strength.

For the toroidal velocity, we performed the mass-weighted average,

[TABLE]

where $V$ is the volume bounded by the TD. As defined, this quantity should correspond to observationally-determinable rotation speeds of well-resolved jets such as that reported in Woitas et al. (2005). All integrations are performed from data taken at the same dynamical time, $l=2$,400 AU, which is comfortably beyond $l\sim 1,\!000$ AU where $\langle v_{\varphi}\rangle$ seems to reach its asymptotic and steady state value for all simulations. As defined, the data in Table 5 and plotted in Fig. 13 show that $\langle v_{\varphi}\rangle\sim B_{\rm i}^{0.67\pm 0.01}$.

Similarly, the quantity $\langle v_{z}\rangle$ is a mass-weighted average of the axial outflow speed given by:

[TABLE]

where here, $V$ is the volume contained by the last 500 AU of the TD and within the jet core in which the outflow speed reaches its asymptotic and steady state value for all simulations. Note that $\langle v_{z}\rangle$ should correspond to observed outflow speeds within jets.

Next, the advance speed of the jet tip into the ambient medium, $v_{\rm jet}$, is measured by fitting the asymptotic slope of the position of the tip of the jet as a function of time, and is plotted in Fig. 13. Such a velocity would be measured using multi-epoch observations of jet bow shocks. In our simulations, we find that $v_{\rm jet}\sim\langle v_{z}\rangle$ (Table 5), and thus the jets are largely ballistic. Given the greater mass density of the nose-cone relative to the ambient medium, and the confinement of the flow by strong toroidal fields, this is not a surprising result. Both $\langle v_{z}\rangle$ and $v_{\rm jet}$ follow power-laws in $B_{\rm i}$ with index $\alpha\sim 0.43\pm 0.01$ which, when combined with the $0.67$ power law for $\langle v_{\varphi}\rangle$, is consistent with the result,

[TABLE]

Finally, $v_{\rm entr}$ is a mass-weighted average of the forward speed of the entrained material (i.e. shocked, ambient material in the simulations identified as having a non-negligible forward velocity but very low $B_{\varphi}$) behind the forward 500 AU of the bow shock (the ‘second wind’). Such a velocity could be measured from CO molecular line observations, for example. It is worth noting, however, that the rotational velocity of the entrained material is negligible. In our simulations, $v_{\rm entr}$ ranges from 11 km s*-1* for simulation A to $1.1$ km s*-1* for simulation H, with a power law index of $\sim 0.56$ (Table 5). This combined with the $0.43$ power law index for $v_{\rm jet}$ is consistent with,

[TABLE]

Should either or both the ‘2/3 law’ in Eq. (26) or the ‘4/3 law’ in Eq. (27) be observed, this could be interpreted as indirect confirmation of the role played by $B_{\rm i}$ in forming and driving protostellar jets. Unfortunately, the proportionality constants in these power-law relationships may be and probably are different from jet to jet, which may make such measures challenging.

5.5 The driving mechanism

The jets are driven by the Lorentz force both directly, as a poloidal gradient in the toroidal magnetic pressure (the ‘magnetic tower mechanism’; MTM), and indirectly whereby the ability of B to exert a substantial ‘normal force’ when rotated allows it to be the ‘rigid wire’ for the ‘bead-on-a-wire mechanism’ (BWM). Both mechanisms can be identified mathematically by a suitable analysis of the relevant forces.

In axisymmetric cylindrical coordinates, the Lorentz force is given by:

[TABLE]

where $\nabla_{\rm p}=\hat{\textbf{{r}}}\partial_{r}+\hat{\textbf{{z}}}\partial_{z}$ is the poloidal gradient, $\varpi_{\varphi}=\frac{1}{8\pi}(rB_{\varphi})^{2}$ is a radially weighted toroidal magnetic pressure, $J_{\varphi}=\frac{1}{4\pi}(\partial_{z}B_{r}-\partial_{r}B_{z})$ is the $\varphi$-component of the current density, $\textbf{{B}}_{\perp}=-B_{r}\hat{\textbf{{z}}}+B_{z}\hat{\textbf{{r}}}$ is a vector perpendicular to and with the same magnitude as the poloidal magnetic field, $\textbf{{B}}_{\rm p}=B_{z}\hat{\textbf{{z}}}+B_{r}\hat{\textbf{{r}}}$, and $b_{\varphi}=\frac{1}{4\pi}rB_{\varphi}$. The last two terms in Eq. (28) are both ‘normal forces’ exerted perpendicular to $\textbf{{B}}_{\rm p}$ ($F_{\perp}$ and $F_{\varphi}$ in Fig. 14), whereas the first term, being the gradient of a function of $B_{\varphi}$ twisted out from the poloidal field, will lie along the general direction of $\textbf{{B}}_{\rm p}$ ($F_{\rm s}$ in Fig. 14). Note further that at $t=0$, $\textbf{{F}}_{\rm L}$ is identically zero since the initial ‘hour-glass’ magnetic field configuration is force-free (Eqs. 15 and 16). It is only after the disc begins to rotate and the magnetic field is distorted from its initial conditions that $\textbf{{F}}_{\rm L}\neq 0$.

For convenience, consider the problem in the co-rotating reference frame of a point on the disc at distance $r_{0}$ from the origin (where the gravitating point mass $M_{*}$ is located; see Fig. 14), whose angular speed is given by $\omega_{0}^{2}=GM_{*}/r_{0}^{3}$. Next, consider the poloidal field line $\psi$ anchored at this point, emerging from the disc at an angle $\theta_{0}$. If we think of the poloidal field line as the ‘rigid wire’ for a ‘bead’ of plasma on the surface of the disc then, as is widely known (e.g. Fig. 1 in Blandford & Payne 1982), the ‘bead’, when nudged, will accelerate out along the ‘wire’ if $\theta_{0}<\theta_{\rm c}=60\degr$.

Now, unlike the classic ‘bead-on-a-wire’ problem found in most sophomore mechanics texts, the poloidal field line is not truly ‘rigid’, regardless of its strength. When the disc starts to rotate, an Alfvén wave with speed $B_{\rm p}/\!\sqrt{4\pi\rho}$ is launched, twisting $\textbf{{B}}_{\rm p}$ in its wake. The stronger $\textbf{{B}}_{\rm p}$ is, the faster the Alfvén wave propagates and the fewer number of turns per unit length suffered by the poloidal field. Still, and regardless of its strength, $\textbf{{B}}_{\rm p}$ will be twisted by numerous full turns over a long enough distance, resulting in a restructuring of the field (rather than a perturbation) that ultimately shuts down and even reverses the effect of the BWM.

To see this, we must include the inertial (centrifugal and Coriolis) forces as observed in the co-rotating frame where the rotational speed of jet material is:

[TABLE]

and the inertial forces are given by:

[TABLE]

where $\mn@boldsymbol{\omega}=\omega_{0}\hat{\textbf{{z}}}$ and $\textbf{{v}}=v_{z}\hat{\textbf{{z}}}+v_{r}\hat{\textbf{{r}}}+v_{\varphi}^{\prime}\hat{\mn@boldsymbol{\varphi}}$.

Combining Eqs. (28) and (29), and resolving the radial inertial force into the two poloidal components, $\hat{\textbf{{s}}}$ and $\hat{\mn@boldsymbol{\perp}}$ as depicted in Fig. 14, we get:

[TABLE]

where $\theta>\theta_{0}$ is the angle between the local field line and $\hat{\textbf{{r}}}$. Other than $-\nabla p-\rho\nabla\phi$ which, in hydrostatic equilibrium, is zero, Eq. (30) represents all the forces acting on a ‘bead’ of matter transported along a poloidal field line ‘wire’, as observed in the co-rotating frame. The omission of the pressure and gravitational gradients from Eq. (30) means any insight gained will be qualitative in nature; to do the full problem with all the forces properly accounted for is why we do the simulations.

The first two terms in Eq. (30) represent the MTM and BWM respectively, and are what drive the jet. The third term is perpendicular to the poloidal magnetic field line within the poloidal plane and, in the steady state, should be zero. If not, the poloidal field line would move, contrary to the assumption (for this analysis) of steady state. The fourth term is a normal force in the symmetry direction, and is responsible for the field line acquiring a toroidal component.

From this it is evident that both the MTM and BWM contribute to the acceleration of the jet, regardless of poloidal field strength; it is simply a matter of which term dominates where. In fact, while the MTM only requires the presence of an outwardly-pointing gradient in the toroidal magnetic pressure—a condition we find true for much of the jet length regardless of the strength of $\textbf{{B}}_{\rm i}$—the conditions for the BWM are much more limiting.

Most importantly, Eq. (30) shows that the BWM term falls to zero once $v_{\varphi}=\omega_{0}r/2$, where $\omega_{0}$ is the Keplerian angular speed at the anchor point. As seen in Fig. 11, this ‘cross-over’ point—where the dashed blue and solid red lines cross—is located before the Alfvén point in all simulations, indicating that the BWM is only effective close to the disc. Table 5 includes the cross-over distance, $s_{\times}$, along the 1 AU field line where the BWM mechanism is shut down. These range from $\sim\!20\,$AU for simulation A to $<1\,$AU for simulation G, with similar values for other field lines in quasi-steady state anchored in the disc between 0.5 and 10 AU. Thus, while the BWM may be the dominant acceleration mechanism on and within short distances from the surface of the disc, the MTM is the dominant acceleration mechanism for most of the jet length in all our simulations.

Indeed, beyond $s_{\times}$, $2v_{\varphi}-\omega_{0}r<0$ and the BWM actually retards outflow, pulling matter back toward the disc! As we see in the simulations, however, there are two reasons why this doesn’t happen. First, poloidal field lines asymptote towards the $\hat{\textbf{{z}}}$-direction and $\cos\theta\rightarrow 0$ in Eq. (30), minimising the BMW term. Second, the MTM term is relentless, counteracting the negative but increasingly feeble BWM as one pulls away from $s_{\times}$. Thus, along just about the entire length of the jet and regardless of $B_{\rm i}$, the gradual, negative gradient in $B_{\varphi}$ gives rise to a net outward acceleration—however modest—well beyond the fast point.

As for the MTM, regardless of the strength of the poloidal field, the rotating disc twists the poloidal field into a dynamically active toroidal field with the passage of the Alfvén wave launched from the disc. To see this, in a time $\Delta t$, the Alfvén wave travels a distance $l=B_{\rm p}\Delta t/\!\sqrt{4\pi\rho}$, during which time the disc has wound up the field into $n=v_{\rm K}(r)\Delta t/(2\pi r)$ coils of radius $r$. Thus,

[TABLE]

since, for quasi-hydrostatic equilibrium, $v_{\rm K}\sim c_{\rm s}/\!\sqrt{\gamma-1}$ (Eq. 19). Thus,

[TABLE]

and $-\nabla_{\rm p}\varpi_{\varphi}/r^{2}$ is a dynamically important outward-pointing force, regardless of the initial poloidal field. This means that in principle, even a trace poloidal magnetic field at the disc surface is sufficient to launch an outflow into the ambient medium, however slow that outflow may be, an observation borne out by our simulation H.

We remind the reader, however, that this analysis is strictly for 2-D axisymmetry which, as has been pointed out, provides the ideal environment for the MTM. In 3-D, the rotation of the fluid which encourages $B_{\varphi}$ is accompanied by numerous modes of instability which discourages its development. We therefore anticipate a much more complicated picture in 3-D, particularly for weaker values of $B_{\rm i}$.

5.6 The knot generator

One of the most striking features of simulation F is the regularity with which ‘knots’ or ‘towers’ (rings or discs in 3-D) are launched from the inner disc, roughly in the region $r_{\rm i}<r<2r_{\rm i}$ (lowermost panel of Fig. 6). As observed in Sect. 4.3, the knots take a little while to establish themselves as various early transients occur in the simulation but, after $\sim 9$ yr, knot production in simulation F remains steady for the remainder of the run.

The only other simulation in which knots of any significance are observed is simulation E, and then only sporadically. Thus, we consider the production of knots in simulation F as ‘transitionary’ between runs with stronger $B_{\rm i}$ characterised by more steady, even laminar flow (e.g., Fig. 3), and runs with weaker $B_{\rm i}$ characterised by much more chaotic flow within the sheath and jet core (e.g., Fig. 8). We note that in the region $0.1<r<0.2\,$AU in the lower panel of Fig. 8, several distinct ‘knot bases’ are present where knots might have formed were the confinement by a poloidal magnetic field more effective.

When produced, the temperature and density within the knots is commensurate with $\rho$ and $T$ near the centre of the gravitational potential well, and these values are maintained as they venture away from the disc. Thus, by the time they reach a distance of $\sim 50$ AU from the disc, their temperature and density can be $10^{3}$ and ten times higher, respectively, than that of the surrounding medium (e.g. Fig. 7). Such knots, therefore, have a thermal pressure excess of four orders of magnitude compared to their immediate surroundings, and it is only the enshrouding poloidal magnetic flux loops that maintain equilibrium. These knots, therefore, may be considered plasmoids (Bostick, 1956) confined by naturally occurring magnetic bottles reminiscent of hydrogen pellet confinement in a thermonuclear reactor.

The mechanism by which knots are launched is very simple and quite unlike that described by Ouyed & Pudritz (1997b). In our simulations, knots are generated directly from the disc, and from poloidal field lines emerging from the disc at or very near the critical angle for the BWM ($\theta_{\rm c}=60^{\circ}$; Sects. 4.1 and 5.5). Although the knots form from warm material, they occur at a transition to a cold medium (bottom panel of Fig. 6), and we indeed observe that the magnetic field angle at this location is much closer to the dynamically cold critical value of $60^{\circ}$ (Blandford & Payne, 1982) than a dynamically warm value of $70^{\circ}$ (Pelletier & Pudritz, 1992). For strong $B_{\rm i}$, field lines are rigid at the disc surface, and plasma is guided subserviently by the BWM. For weak $B_{\rm i}$, field lines are completely at the mercy of the inertia of the plasma, and it is up to the MTM to organise a sufficiently strong toroidal field to launch the jet. However, for simulation F where neither the rigidity of the field lines nor the inertia of the plasma dominate the dynamics, there is much ‘give and take’ at the critical angle.

For field lines of marginal strength, perturbations in the outflow cause the critical field lines to wiggle back and forth across $\theta_{\rm c}$. When $\theta_{0}$ is slightly greater than $\theta_{\rm c}$, hot, dense material near the depth of the gravitational potential accumulates at the base of the field line without moving outward. As the field line is ‘loaded up’, the centrifugal inertia of the growing plasmoid bends the field line outward so that $\theta_{0}$ falls below $\theta_{\rm c}$, and the plasmoid is launched. With the field line now relieved of its mass load, it bends back to something greater than $\theta_{\rm c}$, outflow is quashed, and the cycle repeats. As in any oscillator, the system ‘overshoots’ equilibrium, establishing the regular periodicity observed.

The regularity of knot spacing is indicative of a simple harmonic oscillator whose dependence on $B_{\rm i}$ can be understood as follows. If a force-free magnetic field, $(B_{z_{0}},B_{r_{0}})$, is perturbed by changing the angle at which it emerges from the disc, $\theta_{0}$, by a small angle $\varepsilon\ll\theta_{0}$ (Fig. 15), then the components of the perturbed field become:

[TABLE]

giving rise to a non-zero toroidal current density,

[TABLE]

using $\nabla\cdot\textbf{{B}}=0$. Thus, the restoring force, $\textbf{{F}}_{\perp}$ (Eq. 28), is given by:

[TABLE]

since $\hat{\mn@boldsymbol{\theta}}_{0}$, a unit vector in the direction of increasing $\theta_{0}$, is antiparallel to $\textbf{{F}}_{\perp}$, and $B_{r_{0}}=B_{\rm p}\cos\theta_{0}=B_{\rm p}/2$ for $\theta_{0}=\theta_{\rm c}=60\degr$.

To find the equation of motion for this simple harmonic oscillator, consider the torque density, $\tau$, generated by $\textbf{{F}}_{\perp}$ about point O on the mass within the shaded wedge subtended by angle $2\varepsilon$ in Fig. 15 which, in axisymmetry, is a hollow truncated cone of length $l$ (a portion of this mass is what gets launched to form a knot). For $l\ll r_{0}$, the moment of inertia per unit volume of this wedge is $I\sim\rho l^{2}/2$, and we have using Eq. (31):

[TABLE]

which has the classic form of a simple harmonic oscillator in $\varepsilon$, with a frequency of oscillation given by $\omega\propto a_{\rm p}/\!\sqrt{lr_{0}}\sim B_{\rm i}$.

Since no simulation other than F generated a steady stream of knots, we are unable to test the predicted dependency of $\omega$. Thus, the analysis is included here only as an attempt to identify the physics responsible for the simple harmonic behaviour of the knots in simulation F. We also remark that the generation of knots, although occurring well inside the simulation domain, is a physical result of our prescribed boundary conditions; were the disc self-consistently included in these simulations, it is far from guaranteed that the knots would persist in their observed form.

5.7 Fluxes of mass, momentum, and energy

We now turn to the transport of mass, linear and angular momentum, and kinetic energy within the jet and consider these separately from that transported by material entrained by the bow shock (i.e., the ‘second wind’). In addition to the velocities discussed above, fluxes are among the few concrete physical values that can be determined from observations of jets. All fluxes reported here are axial fluxes measured at a height of $z=1,\!000$ AU above the disc using:

[TABLE]

where the integral is performed as a sum over grid points in the $r$-direction. For fluxes transported by the jet core and sheath, we mask the data by requiring $M>5$ (similar to the masking used to generate Fig. 5), which effectively identifies all material between the jet axis and the TD. For fluxes transported between the TD and bow shock (the second wind), the data are masked by requiring $M<5$. While this includes all points beyond the bow shock as well, there the velocities are zero and thus do not contribute to the fluxes.

Applying the scaling relations in Sect. 3.3, one can convert the fluxes from code to physical units via:

[TABLE]

Figure 16 shows the fluxes calculated from each simulation as a function of time. Notably, each flux reaches (or nearly reaches) an asymptotic value by the end of each simulation, including, to some extent, simulation H. This is confirmed by calculating the fluxes at heights $z=750$ and $1,\!500$ AU which give the same results, short of an appropriate offset in time. We also note that similar loci for fluxes within the second wind (not shown) also converge to asymptotic values, although these data are rather noisier. Typically, second wind fluxes are $\sim\!1$% of the corresponding jet fluxes.

The asymptotic jet fluxes (averaged over the last ten years of each simulation) are included in Table 5. Because the fluxes depend upon the velocities, they inherit the dependence the velocities have on $B_{\rm i}$ (Sect. 5.4), and thus the power law relationships indicated in the Table.

In Table 7, we summarise ranges of velocities ($v_{\rm jet}$, $\langle v_{\varphi}\rangle\sim\langle v_{\rm rot}\rangle$, and $v_{\rm entr}$) from Table 5 along with ranges for the fluxes measured from the last 10 yr of each simulation. The upper (lower) half of the table includes jet (second wind) fluxes where $M\geq 5$ ($M<5$). On comparing Tables 1 and 7, it is evident that these aspects of our simulations agree comfortably with the observations, including both ‘jet’ and ‘entrained’ fluxes. The overlap is not perfect, of course. Our mass fluxes fall within the observational range while the observed angular momentum fluxes fall within the range of our simulations. As extensive as these simulations are, they are still limited by initial conditions, geometry, evolution time, and even some physics. With the relaxation of any of these constraints, specific comparisons may well improve. Furthermore, the results presented in Table 7 depend somewhat on the values adopted in the scaling relations (Eqs. 20–24) for the stellar mass ($M_{*}$), inner disc radius ($r_{\rm i}$), plus initial magnetic field and plasma-$\beta$ values at $r_{\rm i}$ ($B_{\rm i}$, $\beta_{\rm i}$), and should be adjusted appropriately for specific comparisons. Nevertheless, we take the agreement between Tables 1 and 7 as encouragement that these global simulations have captured the essence of both the production and propagation of protostellar outflows, and are among the first to do so.

We note that observational measurements of the linear momentum flux in jets are rare. Still, that the few values of which we are aware (see Podio et al. 2006; Hartigan et al. 1994; $10^{-6}$ – $1.4\times 10^{-4}\,M_{\odot}\,{\rm yr}^{-1}\,{\rm km}\,{\rm s}^{-1}$) fall within the range of our calculated fluxes is encouraging. Meanwhile, to the best of our knowledge, there are no observational estimates of kinetic energy fluxes in protostellar outflows, so our values in Table 7 serve as a prediction should such measurements ever be made.

6 Summary and Conclusions

We present the first simulations to resolve the inner launching region of a protostellar (non-relativistic) jet, and then follow the jet to observational length scales. The simulations were performed with our adaptive grid MHD code, AZEuS (Ramsey et al., 2012), with an effective dynamic range in resolution of $\sim 6.5\times 10^{5}$.

With a putative protostellar mass of $0.5\,M_{\odot}$, results from eight axisymmetric simulations are discussed, each identical except for the value of the plasma-beta at the inner radius of the accretion disc, $\beta_{\rm i}$, ranging from $0.1$ (simulation A) to 640 (simulation H). In each simulation, a jet is launched from the inner portion of a Keplerian accretion disc maintained as boundary conditions by magneto-rotational dynamics on the grid, and outflow is followed until the leading tip of the bow shock excited in the primordial atmosphere reaches the end of the coarsest grid. One of the primary features of these simulations that sets them apart from others is at no time does any part of the outflow leave the computational domain.

While our jets are still extremely ‘young’ (lengths range from $\sim\!2,\!300$ to $\sim\!4,\!100$ AU; ages range from $\sim\!50$ to $\sim\!150$ yr), each has developed sufficiently to be directly comparable to observed jets, and to have established certain key observational properties such as asymptotic speeds, mass, momentum, and energy fluxes, rotation rates, etc.

Based on their characteristics, we have divided our eight simulations (Table 4) into three sub-groups: ‘strong-field’ (simulations A–D), ‘medium-field’ (simulations E, F), and ‘weak-field’ (simulations G, H). With this distinction in mind, we draw the following conclusions:

All eight simulations generate an organised, supersonic, magnetically confined outflow. 2. 2.

In regions of quasi-steady state, the ‘Weber-Davis constants’ (equations 1–4) remain constant for the most part to within 3%, attesting to the numerical integrity of our AMR-MHD scheme. 3. 3.

On the observational scale ($>\!1,\!000\,$AU), each jet resembles the MHD simulations of Clarke et al. (1986). Thus, they are confined by an internal toroidal magnetic field twisted out of the poloidal magnetic field in the atmosphere, and develop: (i) a hot, dense, super-fast ‘jet core’ we identify with the observed jet; (ii) a rarefied, cold, magnetised, trans-fast ‘sheath’ which we suggest may show up as an ‘emission cavity’ surrounding the jet; and (iii) a narrow, laminar, magnetically confined ‘nose-cone’ leading the jet. Surrounding the jet is a tangential discontinuity (TD) that separates jet material from the warm, trans-slow, shocked ambient medium entrained by the bow shock excited by the passage of the jet (‘second wind’; top panel of Fig. 3). Due to the very large dynamic range and adaptive resolution, to the best of our knowledge, this is the first study of its kind that self-consistently includes and clearly differentiates between these different outflow components. 4. 4.

In agreement with local studies of outflow launching, strong-field jets (e.g. Fig. 3) are initially launched by the ‘bead-on-a-wire mechanism’ (BWM; Blandford & Payne 1982) but is shut down even before the Alfvén point. Beyond that, the ‘magnetic tower mechanism’ (MTM; Lynden-Bell 1996) takes over, accelerating outflow well beyond the fast point. They are characterised by quasi-steady state flow even near the jet launching region. 5. 5.

Also in agreement with local studies, weak-field jets (e.g. Fig. 8) are launched and accelerated by the MTM, in which relentless twisting of the weak poloidal field builds up sufficient toroidal magnetic pressure to drive the outflow. These are characterised by highly turbulent flow near the launching region, with turbulence encroaching the TD as $\beta_{\rm i}$ increases. For very high $\beta_{\rm i}$, we speculate that the encroachment of turbulence on the TD will cause enough mixing between the shocked jet and ambient media to disrupt the jet in 3-D (Hardee et al. 1997). 6. 6.

Medium-field jets (e.g. Fig. 6) are transitional between the strong and weak field cases. They are launched by the BWM but accelerated soon thereafter by the MTM. Neither laminar nor turbulent in the jet launching region, medium-field jets can exhibit ‘knot-like’ structures which are generated periodically from the inner disc and propagate with the flow. 7. 7.

With the exception of simulation H (our weakest-field simulation), this study confirms that the interior dynamics of jets are dominated by a strong toroidal magnetic field. Indeed, we find that plasma-$\beta$ averaged over rotating regions where $M>5$, $\langle\beta_{\rm tot}\rangle$, asymptotes to 0.2–0.4 regardless of $\beta_{\rm i}$. Simulation H struggles to maintain outflow, but even still $\langle\beta_{\rm tot}\rangle\rightarrow 1$. Thus, this study implies that measures of magnetic field strength within the jet may not reveal much about the magnetic environment near the protostar. 8. 8.

The global nature of these simulations reveal that jets at similar ‘dynamical times’ (e.g., of equal length) have similar radii and bow shock shape. Thus, these cannot be used as a measure of magnetic field strength near the protostar. 9. 9.

Since no part of the outflow leaves the domain, we are able to determine that the advance speed of the jet (as would be measured by time-lapse images), $v_{\rm jet}$, and the average flow speed within the leading portion of the jet (as could be measured from line emission observations), $\langle v_{z}\rangle$, both vary as $B_{\rm i}^{\sim 4/9}$, where $B_{\rm i}$ is the initial magnetic field strength at the inner radius of the disc. The fact that $v_{\rm jet}\sim\langle v_{z}\rangle$ indicate the jets are essentially ballistic. 10. 10.

The average rotation speed along the latter portion of the jet varies as $\langle v_{\varphi}\rangle\sim B_{\rm i}^{2/3}$. This, along with our other results, leads to possible observable evidence for the magnetic character of protostellar jets: $v_{\rm jet}\sim\langle v_{\varphi}\rangle^{2/3}$ (Eq. 26). 11. 11.

The jet advance speed, $v_{\rm jet}$, and the advance speed of material entrained by the bow shock, $v_{\rm entr}$, all fall nicely within the realm of observational constraints (Table 1). Together, they indicate $v_{\rm entr}\sim v_{\rm jet}^{4/3}$ (Eq. 27). 12. 12.

Because no jet leaves the computational domain, we have been able to make estimates of fluxes for mass, momentum (linear and angular), and kinetic energy in our jets. Based in part on the consistency of these fluxes and velocities with those measured from observational data (Tables 1 and 7), we conclude that our global simulations are able to make the link to observations where local simulations cannot. In particular, we have shown that numerical models based solely on gravito-magneto-rotational fluid dynamics are capable of launching and driving jets that are consistent morphologically and quantitatively with the observations.

Acknowledgements

We thank ACEnet technicians Phil Romkey and Sergiy Khan for their assistance in optimising AZEuS for the ACEnet facilities. This work was supported, in part, by an NSERC Discovery Grant to DAC. The Centre for Star and Planet Formation is funded by the Danish National Research Foundation (DNRF97). JPR was supported, in part, by the Virginia Initiative on Cosmic Origins (VICO). Computing resources were provided, in part, by ComputeCanada via ACEnet and Calcul Quebéc. Additional computing facilities were provided by the University of Copenhagen HPC centre, funded in part by Villum Fonden (VKR023406), and the University of Copenhagen Electronic Research Data Archive (ERDA). This work made use of the SAO/NASA Astrophysics Data System and the MPFIT least-squares fitting package (Markwardt, 2009). We also thank the anonymous referee for a timely and constructive report.