Linear and nonlinear evolution of current-carrying highly magnetized jets
aa r X i v : . [ a s t r o - ph . H E ] J u l Mon. Not. R. Astron. Soc. , 1–14 (2014) Printed 28 September 2018 (MN L A TEX style file v2.2)
Linear and nonlinear evolution of current-carrying highlymagnetized jets
M. Anjiri , ⋆ , A. Mignone , G. Bodo and P. Rossi Dipartimento di Fisica Generale “Amedeo Avogadro” Universit`a degli Studi di Torino, Via Pietro Giuria 1, I-10125 Torino, Italy Department of Physics, School of Sciences, Ferdowsi University of Mashhad, Mashhad, 91775-1436, Iran INAF/Osservatorio Astronomico di Torino, Strada Osservatorio 20, I-10025 Pino Torinese, Italy
28 September 2018
ABSTRACT
We investigate the linear and nonlinear evolution of current-carrying jets in a periodicconfiguration by means of high resolution three-dimensional numerical simulations.The jets under consideration are strongly magnetized with a variable pitch profile andinitially in equilibrium under the action of a force-free magnetic field. The growthof current-driven (CDI) and Kelvin-Helmholtz (KHI) instabilities is quantified usingthree selected cases corresponding to static, Alfv´enic and super-Alfv´enic jets.During the early stages, we observe large-scale helical deformations of the jetcorresponding to the growth of the initially excited CDI mode. A direct comparisonbetween our simulation results and the analytical growth rates obtained from lineartheory reveals good agreement on condition that high-resolution and accurate dis-cretization algorithms are employed.After the initial linear phase, the jet structure is significantly altered and, whileslowly-moving jets show increasing helical deformations, larger velocity shear are vi-olently disrupted on a few Alfv´en crossing time leaving a turbulent flow structure.Overall, kinetic and magnetic energies are quickly dissipated into heat and duringthe saturated regime the jet momentum is redistributed on a larger surface area withmost of the jet mass travelling at smaller velocities. The effectiveness of this processis regulated by the onset of KHI instabilities taking place at the jet/ambient interfaceand can be held responsible for vigorous jet braking and entrainment.
Key words: instabilities - ISM: jets and outflows - stars: jets - MHD - methods:numerical
The investigation of instabilities on the propagation ofcollimated magnetized flows has been an outstandingresearch subject during the last 30 years (Cohn 1983;Ferrari & Trussoni 1983; Hardee et al. 1992; Todo et al.1993; Appl 1996). These instabilities have a substantialimportance in characterizing the mechanism of formationand evolution of various observed structures in such flows(Bahcall et al. 1995; Biretta 1996; Raga & Noriega-Crespo1998; Rosado et al. 1999; Nakamura & Meier 2004;Nakamura et al. 2007; Mignone et al. 2010). Broadly speak-ing, instabilities in astrophysical jets fall in two classes:external instabilities (caused by the relative motion betweenjet and external medium) and intrinsic instabilities (relatedto the toroidal component of magnetic fields). ⋆ E-mail:[email protected] (MA)
One of the most important external instabilities isthe Kelvin-Helmholtz instability (KHI), responsible forthe interaction and mixing between the jet and ambi-ent medium as well as the transfer of linear momen-tum and jet braking (Bodo et al. 1995; Hardee et al. 1995;Bodo et al. 1998; Rossi et al. 2008). The KHI generallyleads to distortions of the interface between the jet andthe ambient fluids, eventually producing shocks and tur-bulent mixing with ambient material. Analyses of theKHI have been extensively accomplished by many re-searchers either in a linear regime (Turland & Scheuer 1976;Payne & Cohn 1985; Hardee et al. 1992; Bodo et al. 1996;Hardee 2000; Perucho et al. 2004a; Perucho & Lobanov2007; Osmanov et al. 2008) or a nonlinear regime(Bodo et al. 1994; Hardee et al. 1997; Bodo et al. 1998;Rosen et al. 1999; Rossi et al. 2004; Perucho et al. 2004b,2005; Rossi et al. 2008). For example, the linear analysisby Perucho et al. (2004a) showed that the linear growth of c (cid:13) Anjiri et al.
KH modes is smaller for faster and colder relativistic jets.Perucho & Lobanov (2007) indicated that the growth ratesof KHI modes are significantly reduced by the presence ofa thick shear-layer. This important result was confirmed bynumerical simulations as well (see, e.g, Perucho et al. 2005,2007). Some results have been derived for cylindrical rela-tivistic jets that clarified the importance of the nonlinearevolution of KHI in the dichotomy of FRI/FRII extragalac-tic radio sources (Rossi et al. 2004, 2008). In this respect,the hydrodynamic simulations performed by Rossi et al.(2008) on relativistic light jets have shown that the non-linear growth of KHI promotes a strong interaction betweenthe jet and the external medium with a consequent mixingand remarkable deceleration.On the other hand, intrinsic instabilities are gener-ally related to the magnetohydrodynamic (MHD) struc-ture of the jet. In this sense, a key role is played bythe relative strength between the poloidal and toroidal (orazimuthal) components of magnetic field. This ratio, re-ferred to as the magnetic pitch parameter, plays an impor-tant role in triggering the so-called current-driven instabil-ities or CDI. According to the mechanisms responsible forjet acceleration and collimation (Blandford & Payne 1982;Romanova & Lovelace 1992), the toroidal magnetic field is asignificant component that is expected to dominate far fromthe central engine. Its destabilizing action is at the base ofthe CDI (Bateman 1980) and can deeply affect the morpho-logical structure of jets (Ferrari et al. 2011). In this respect,the observations and analytical models presented for the jetin M87 (Sikora et al. 2005; Hardee 2006, 2011; Walker et al.2008, 2009) confirm the role of the CDI in conversion Poynt-ing flux to kinetic flux flow within a few hundred gravita-tional radii. Nevertheless, there it has been suggested overthe years that this flux conversion process may be accompa-nied by other efficient mechanisms such as matter entrain-ment and jet expansion that can significantly change the na-ture of jets from magnetically (sub-Alfv´enic) to kinetically(super-Alfv´enic) dominated (Hardee 2006, 2011).Among CDI, the | m | = 1 (or kink) mode is themost important one and grows faster than the other modes(Appl et al. 2000; Nakamura & Meier 2004). In the kink in-stability, magnetic field lines are compressed on the innerside of a deformed cylindrical flux tube and the magneticpressure exerted by toroidal component becomes larger thanthe net magnetic tension. Thus, the jet curvature magni-fies so that some of the magnetic energy accumulated bythe twisting is released by a kink (Spruit 1996). Eventually,the hoop stress provided by the toroidal field is reduced(Eichler 1993) and leads to a whole helical deformation jet(Mizuno et al. 2009, 2011) and/or even jet complete disrup-tion (Nakamura & Meier 2004).Linear perturbative analysis of the CDI in MHD jetshas been largely performed under different physical con-ditions by several groups (Cohn 1983; Appl & Camenzind1992; Eichler 1993; Appl 1996; Appl et al. 2000; Wanex2005; Bonanno & Urpin 2011). In most of them the force-free limit is considered for helical magnetic field. For in-stance, Appl et al. (2000) showed that for cold supermagne-tosonic jets with dominant azimuthal fields, the CD modesare confined to the jet interior and develop quickly on atime scales of the order of Alfv´en crossing time in a frameof reference co-moving with the jet. Wanex (2005) demon- strated that jets with axial sheared flow and sheared mag-netic fields reduced the linear growth of CD modes at con-stant magnetic pitch. In addition, Bonanno & Urpin (2011)considered the case with both azimuthal and axial mag-netic fields and showed that the length scale of CD kinkmodes is strongly sensitive to the magnetic pitch parame-ter. In the relativistic MHD jets, the linear analysis of CDIby Istomin & Pariev (1994, 1996) demonstrated that cylin-drical jets with constant axial (poloidal) magnetic fields arestable against kink CD modes. More recently, the studiesby Bodo et al. (2013) on the non-rotating magnetized jetsrevealed the splitting of the CD kink modes into an innerand outer modes at high magnetization.On the other side, several numerical studies haveconcentrated on the nonlinear development of CD in-stabilities in MHD jets under various physical assump-tions (Todo et al. 1993; Lery & Frank 2000; Lery et al.2000; Baty & Keppens 2002; Nakamura & Meier 2004;Nakamura et al. 2007; Carey & Sovinecy 2009). In the workby Lery et al. (2000), the non-linear analyses of CDI for coldsuper fast magneto-sonic jets demonstrates that the CDIplay a significant role in re-distributing the current densityin the inner parts of the jet. Baty & Keppens (2002) consid-ered the nonlinear interaction between CDI and KHI sur-face modes on the propagation of supersonic jets showingthat CDI prevent the development of KH vortices at the jetsurface. According to their results, the magnetic field de-formations induced by the nonlinear growth of CDI modesprovide a stabilizing factor on the KHI-driven vortical struc-tures leading to a substantial decrease in the mixing processbetween jets and ambient medium and preventing disruptiveeffects.More recently, Mizuno et al. (2009) studied the devel-opment of kink instabilities on helically magnetized rela-tivistic static columns showing that, for small pitch val-ues, the growth rate of instability rapidly increases duringthe linear stages and that the nonlinear evolution featuresa continually growing helically twisted column. These au-thors further outstretched this work by considering a sub-Alfv´enic velocity shear surface (Mizuno et al. 2011), show-ing that the temporal evolution of the CDI is dramaticallyreduced with respect to the static column. Additionally,O’Neill et al. (2012) assumed a similar approximation tothat of Mizuno et al. (2009) and studied various local mod-els of co-moving magnetized plasma columns in force-free,rotational, pressure-confined equilibrium configurations andfound that the details of initial force balance strongly affectthe resulting column morphology.In this work, we investigate the stability of stronglymagnetized jets with initial equilibrium structure describedby a force-free helical magnetic field. We adopt a periodicconfiguration representative of a jet section far from thelaunching region and consider radially sheared, axial flowswith different velocities (Section 2). The configurations aredestabilized using an exact eigenfunction corresponding thefastest growing m = 1 CDI mode of the linearized MHDequations and the evolution is followed through the lin-ear and nonlinear phases using three-dimensional numericalsimulations (Section 3). We first assess the impact of gridresolution on the growth of the CDI mode during the lin-ear phase by performing a close comparison with the resultsfrom normal mode analysis. Hence we explore the nonlinear c (cid:13) , 1–14 inear and nonlinear evolution of MHD jets behaviour in terms of jet morphology, shear-induced effectstriggered by the onset of KHI modes with particular at-tention to jet braking and momentum transfer from the jetto ambient medium as the system approaches the saturatedregime. Contrary to previous studies, our results are relevantto jets with a smaller plasma- β ( β ≈ − ) and are basedon much larger numerical resolutions ( &
20 zones on the jetradius). Finally, our findings are summarized in Section 4.
In the following, we will investigate the dynamical evolutionof current-carrying jets by means of three-dimensional nu-merical simulations. The simulations are performed by solv-ing the time-dependent ideal MHD equations in Cartesiancoordinates ∂ρ∂t + ∇ · ( ρ v ) = 0 , (1) ∂ ( ρ v ) ∂t + ∇ · h ρ vv − BB i + ∇ p t = 0 , (2) ∂ E ∂t + ∇ · [( E + p t ) v − ( v · B ) B ] = 0 , (3) ∂ B ∂t + ∇ · ( vB − Bv ) = 0 , (4) ∂ ( ρ T ) ∂t + ∇ · ( ρ T v ) = 0 , (5)where ρ, v and B denote the mass density, bulk velocityand magnetic field, respectively. Note that a factor 1 / √ π has been absorbed in the definition of B . The total (gas+ magnetic) pressure p t is denoted with p t = p + B / E = pΓ − ρ v B , (6)where the internal energy obeys the perfect gas law withspecific heat ratio Γ = 5 /
3. Equation (5) sets the evolutionof a dynamically passive scalar field that is used to trackfluid elements initially residing in the jet ( T jet = 1 for R < T jet = 0 otherwise) or to mark surfaces of constantmagnetic flux which, for our case, are initially concentriccylinders ( T mag = R ).Alternatively to equation (3) we also solve, during theinitial stages of evolution for t <
35, the entropy equationaway from shocks ∂s∂t + v · ∇ s = 0 , (7)where s = p/ρ Γ . In a highly magnetized plasma, solvingequation (7) instead of equation (3) has the advantage ofbeing more accurate and preventing the occurrence of neg-ative pressure values which could otherwise be triggered bythe truncation error of the scheme when retrieving p fromequation (6).The numerical computations are carried using the MHDmodule of the PLUTO code for astrophysical gas dynamics(Mignone et al. 2007). In order to solve equation (1)-(4), we use a third-order algorithm based on a second order Runge-Kutta time stepping, piecewise linear/parabolic reconstruc-tion and the HLL Riemann solver. The employment of amore accurate Riemann solver such the Roe or HLLD Rie-mann solver (see the discussion in O’Neill et al. 2012) leads,unfortunately, to severe numerical difficulties when dealingwith such low-beta plasma configurations. Our initial condition consists of an infinitely long axisym-metric jet moving in the vertical direction. The jet has aninitial equilibrium structure that depends on the cylindri-cal radius r only and characterized by a force-free magneticfield B ( r ) = B φ ( r )ˆ φ + B z ( r )ˆ z , constant gas pressure andabsence of rotations ( v φ = v r = 0). In such a way, thetoroidal and poloidal magnetic field components obey thetime-independent radial component of the momentum equa-tion (equation 2) 1 r d ( r B φ ) dr + dB z dr = 0 , (8)while density and vertical velocity can be chosen arbitrar-ily as they do not explicitly appear in the radial balanceequation. Here, we set ρ = ρ h η + 1 − η cosh( r/r j ) i , v z = v z cosh( r/r j ) , (9)where ρ is the on-axis density, r j is the jet radius, η isthe ambient to jet density contrast and v z is the jet ax-ial velocity. Gas pressure is initially constant and equal to p ( r ) = ρ c s where c s is the speed of sound.Following Bodo et al. (2013) we prescribe the azimuthalcomponent of magnetic field to be B φ = − H c r j r s − exp (cid:18) − r a (cid:19) , (10)where a = 0 . r j is the magnetization radius and H c de-termines the maximum field strength. The chosen profilecorresponds to a current mainly distributed inside the jetand peaked on the axis. Equation(10) yields a current-freefield at large distances as B φ ∼ /r for r → ∞ and a linearprofile, B φ ∼ r , for r → B z = H c r j a s P c a − √ π erf (cid:18) r a (cid:19) , (11)where erf() is the error function while P c = lim r → (cid:12)(cid:12)(cid:12)(cid:12) rB z B φ (cid:12)(cid:12)(cid:12)(cid:12) (12)is the magnetic pitch parameter on the axis.The strength of the magnetic field is controlled by thevalue of H c which, without loss of generality, is fixed by thecondition that the average Alfv´en speed over the jet beamis always unity:¯ v A ≡ s r j ρ Z r j ( B φ + B z ) r dr = 1 . (13) c (cid:13) , 1–14 Anjiri et al.
Figure 1.
Radial profile of the jet axial velocity (left-hand panel), poloidal and toroidal magnetic field B z and B φ (middle panel) andmagnetic pitch P = rB z /B φ (right-hand panel). The solid, dotted and dashed lines in the left-hand panel refer to A10, A1, A0 cases,respectively. Table 1.
Simulation cases and parameters describing the initial jet configuration. Here M s ≡ v z /c s and M A ≡ v z / ¯ v A are, respectively,the sonic and Alfv´enic Mach numbers, L h and L max define the horizontal extent of the computational domain (see the text), L z = 2 π/k is the vertical domain extent (equal to one perturbation wavelength). The number of points in the three directions together with thenumber of zones per jet radius ( N j ) are given in columns 8 and 9, and columns 10 and 11 for low and high resolution computations,respectively. Finally, k and − Im( ω ) represents the wavenumber and growth rate of the fastest-growing CDI mode shown in Fig. 2.Low resolution High resolutionCase v z M s M A L h L max L z N x × N y × N z N j N x × N y × N z N j k − Im( ω )A0 0 0 0 8 . . × × ≈
12 576 × × ≈
24 1 .
164 0 . . . × × ≈
12 640 × × ≈
25 0 .
806 0 . .
85 30 29 . × × ≈
11 640 × × ≈
22 0 .
213 0 . Likewise, we take the jet density and radius to be our refer-ence density and length, i.e., ρ = 1 and r j = 1.Our equilibrium jet model is thus written in terms ofthe pitch parameter P c , the jet velocity v z (in units of theAlfv´en speed), the jet density contrast η and the sound speed c s . In the remainder of this paper, we will set η = 1 and c s = 1 /
10 corresponding to highly magnetized jets with β = 2 p/ B ≈ − and density equal to that of the am-bient medium. Also, we take P c = 0 . B z ( r ) and B φ ( r ) (and the pitch) are the same in all cases.We consider three simulation cases with different val-ues of the axial velocity v z = 0 , ,
10 corresponding, re-spectively, to the static jet case, a trans-Alfv´enic jet and toa super-Alfv´enic flow. The three cases will be referred toas A0, A1 and A10 and are listed, together with specificsimulation parameters in Table 1. Each simulation case isrepeated twice using low and high resolution.The computational domain is the Cartesian box definedby x, y ∈ [ − L max , L max ] and z ∈ [0 , L z ], where L z = 2 π/k corresponds to one wavelength of the chosen eigenmode (see § N x , N y and N z , respectively, and is given in Table1 for low and high resolutions runs. The grid has uniformspacing ∆ x = ∆ y = ∆ z = L z /N z in the vertical directionand in the region x, y ∈ [ − L h , L h ], where L h < L max (see Table 1). Outside of this region, for | x | , | y | ∈ [ L h , L max ], themesh spacing increases in geometrical progression. In orderto obtain cubic cells in the central region of the domain,the number of zones in the x and y directions is chosento be 2 N z L h /L z while the stretched portions of the gridare discretized symmetrically with N x / − N z L h /L z and N y / − N z L h /L z zones on each segment. The number ofzones used to resolve the jet radius, N j , is reported in Table1 and is, to the extent of our knowledge, the largest one em-ployed so far in simulations of low plasma- β in periodic jetconfigurations. This provides finer resolution in the neigh-bourhood of the jet where most of the dynamics is expectedto take place. Finally, we use periodic boundary conditionsin the vertical direction while outflow conditions hold at theremaining boundaries. In order to ease up the comparison with the predictions fromglobal normal mode stability analysis, we perturb the equi-librium configuration using an exact eigenmode of the lin-earized MHD equations (Appl & Camenzind 1992). Setting V = ( ρ, v , B , p ) the array of fluid variables, we perturb theinitial condition as V = V + δ V , where V are the equi-librium profiles described in Section 2.2 while δ V = ǫ Re h Y ( r ) exp (i ωt − i kz − i mφ ) i (14)is the perturbation. Here Re[ ] denotes the real part, Y ( r )is a complex eigenfunction, ω ( k ) is the corresponding (com- c (cid:13) , 1–14 inear and nonlinear evolution of MHD jets Figure 2.
Linear growth rates for m = 1 for different jet axial velocities as predicted from normal mode analysis. From left to right,we plot the CD and KH growth rates as a function of wavenumber k for different jet axial velocity v z = 0 , , ,
10. The CD mode isplotted using the solid blue line while the fundamental KH mode shows as a solid red line. The purple curves result from the merging andsplitting of the CD and KH modes through the ’X’ point. The dot-dashed and dashed lines represent reflected modes while the verticalline gives the perturbation wavenumber used in the numerical simulation. plex) eigenvalue, k and m are the (real) axial and azimuthalwavenumbers and ǫ is a small real number chosen in sucha way that the transverse velocity perturbation amplitudeequals 0 .
01. The eigenfunction Y ( r ) is determined by solv-ing a boundary value problem with appropriate conditionsat small and large radii see, for instance, Appl & Camenzind(1992); Appl et al. (2000) or Bodo et al. (2013) for the rel-ativistic treatment. Instability occurs when Im[ ω ( k )] < k and m are realnumbers while the growth rate of instability is given by theimaginary part of the complex eigenvalue ω . A plot of thegrowth rate for m = 1 as a function of k is given in Fig. 2for different velocities. For v z = 0 (blue curve, first panel),only the CDI mode is present while the appearance of avelocity shear (equation 9) triggers additional KHI modes.When v z = 1 a new KH mode with comparable growthrate appears at larger wavenumbers (red curve in the secondpanel of Fig. 2) and partially overlaps with the CD mode.The CD and KH growth rates curves intersect through an’X’ point at k ≈ t ≈ . v z = 4) consists of a mode with larger growth rate anda second mode with smaller amplitude with mixed CD/KHproperties. By further increasing the velocity they reach theconfiguration shown in the right-hand panel of Fig. 2.In addition we note that reflected KH modes(Bodo et al. 1989) appear as well for v z = 10. Overall, weexpect both KH and CD modes to play a role although itmay not be possible to unequivocally isolate their contribu-tions. We choose the eigenfunction corresponding to thewavenumber k = k at which the growth rate is approxi-mately maximum (vertical dotted line in Fig. 2) and reportthe precise value in Table 1. Beware that only modes withwavenumbers given by an integer multiple of k can actuallydevelop in our computational box.Finally, since the perturbation Y ( r ) in equation (14) isknown as a function of the cylindrical radius, we use bilinearinterpolation to obtain the values on our discrete Cartesiandomain at t = 0. In the following sections we present our simulation resultsseparately according to their linear and nonlinear evolu-tion. Several diagnostic are computed in a way similar toMignone et al. (2013) by introducing the horizontal averageoperator ¯ Q ( z, t ) = h Q, χ i = Z Q ( x , t ) χ ( x , t ) dx dy Z χ ( x , t ) dx dy , (15)where Q ( x , t ) is any flow quantity, χ ( x , t ) is a weight func-tion and x = ( x, y, z ) is the position vector. Likewise, wedefine the volume average operator as¯ Q ( t ) = hh Q, χ ii = Z Q ( x , t ) χ ( x , t ) dx dy dz Z χ ( x , t ) dx dy dz . (16)When averaging with χ = 1 we use the short-hand notation hh Q ii ≡ hh Q, ii . During the initial stages of the evolution, the perturbationgrows and the system slowly departs from equilibrium. Aswe shall see, this phase is characterized by the dominanceof the m = 1 CDI mode. c (cid:13) , 1–14 Anjiri et al.
Figure 3.
The evolution of transverse velocity ¯ v tr as a function of time for Case A0 (left), Case A1 (middle) and Case A10 (right).Dashed and solid lines represent the low-resolution and the high-resolution runs, respectively. Dotted line is the reference slope accordingto the linear theory. Velocity is in units of the Alfv´en speed v A while time is units of r j /v A . Figure 4.
Power spectra of the radial displacement for the Case A0 (left-hand panel), Case A1 (middle panels) and Case A10 (right-handpanel). The plot shows P k |F mag ( m, z k ) | /N z and the discrete Fourier transform is taken using the perturbations of T mag . As an indicator of the growth of the instability, we use thevolume-average of the transverse velocity, defined by¯ v tr ( t ) = DDq v x + v y EE , (17)where v x and v y are the horizontal components of velocity.Equation (17) allows us to perform a direct comparison withthe result of linear stability analysis, as shown in Fig. 3where we plot the time history of the averaged transversevelocity ¯ v tr for both low and high resolution computations.Complementary, we also quantify the growth of the in-stability by performing a Fourier decomposition of the quan-tity ∆ T mag ( x , t ) = T mag ( x , t ) − T mag ( x ,
0) which, in thelimit of small departures from equilibrium, represents es-sentially the linear radial displacement. In order to computethe discrete Fourier transform in the φ direction, we con-sider all the zones lying in a small annular region satisfying | p x + y − | < ∆ x and then interpolate the resulting se-quence of values on an azimuthal 1D grid using N φ regularlyspaced points. The discrete Fourier transform is then com-puted as F mag ( m, z k ) = 1 N φ N φ − X j =0 ∆ T mag ( r , φ j , z k ) e − i mφ j (18)where r = 1 is a fiducial radius, ∆ T mag ( r , φ j , z k ) rep- resents the regularly gridded data in the φ direction and m is the azimuthal wave number. The power spectrumis then computed by averaging in the vertical direction: P k |F mag ( m, z k ) | /N z and it is shown in the three panelsof Fig 4.In the static column case most of the main source of en-ergy is magnetic and only the CD mode is present. For t . t ≈ r j / ¯ v A ) and ( t ≈ r j / ¯ v A ) for thelow- and high-resolution runs, respectively and the measuredgrowth rate (averaged between 6 < t <
10) is ω Lo ≈ . ω Hi ≈ . ω A = 0 . t = 24,the prevalence of the m = 1 mode is evident.The evolution of the Alfv´enic jet (Case A1) takes placeon a time scale similar to the static column case and theduration of the linear phase is approximately the same ( t . ω Lo ≈ .
247 and ω Hi ≈ .
246 for (8 < t <
12) showing convergence for finermesh spacing ( ω A = 0 . t = 24 clearlyshowing that the kink mode dominates.In the super-Alfv´enic jet, the grid resolution plays a cru- c (cid:13) , 1–14 inear and nonlinear evolution of MHD jets Figure 5.
Three-dimensional structures of the jets for Case A0 (left-hand panels), Case A1 (middle panels) and Case A10 (right-handpanels) for the high resolution runs. The top panels show magnetic pressure isosurfaces (Pm, blue-green) together with density isosurfaces(orange-red). Bottom panels display magnetic flux surfaces (i.e. surface with constant T mag ) coloured by field intensity and superimposedon magnetic field lines. cial role in determining the numerical value of the growthrate (right-hand panel in Fig. 3). Indeed, from several nu-merical experiments not reported here, we found this par-ticular case to be the most challenging one owing to thelarge discretization noise arising from the truncation errorof the scheme. The noise triggers grid-sized additional per-turbations with sufficiently large amplitudes on top of theinitially chosen CDI eigenmode. We note that this spuriouseffect could be reduced by the combined action of a higherorder interpolation scheme (PPM) and by solving for theentropy equation rather than the total energy equation. Inaddition, the results obtained at higher resolution improveappreciably over the low-resolution computations. A linearphase is observed for t .
20 and the measured growth ratesare ω Lo = 0 .
322 and ω Hi = 0 .
321 for (3 < t <
9) at lowand high resolutions, respectively, to be compared with thetheoretical value of ω A = 0 . m = 1 mode is clear fromthe right-hand panel of Fig. 4. The three-dimensional structures of selected jet cases arevisible in the panels of Fig. 5 where density and magneticpressure isosurfaces (top) and magnetic flux surfaces definedby the condition T mag = const (bottom) are shown at t = 24(for case A0 and A1) and t = 15 (for case A10).In the static column case (left-hand panel), the ini-tial displacement grows into a twisted helical deformationof the column with density enhancements corresponding toregions of smaller magnetic perturbation and viceversa. Thisgives rise to a double-stranded spiral structure in whichthe outward magnetic flux tube stretches and its diameterwiden progressively. Owing to the flux-freezing condition,mass and internal energy inside the flux tube are depleted.A similar structure has been also found by other authors(Baty & Keppens 2002; Mizuno et al. 2009; O’Neill et al.2012) using different setups. The end of the linear phase ismarked by a density build-up on the axis with flow velocitiesof the order of the Alfv´en speed.In the Alfv´enic jet case, the three-dimensional structureat t = 24 (middle panel in Fig. 5) reveals some similaritieswith the static column case although we observe the for- c (cid:13) , 1–14 Anjiri et al.
Figure 6.
Top panels: volume-integrated jet energy fractions as a function of time for the three simulation cases. Kinetic, magnetic andthermal contributions are shown using blue, green and red coloured lines, respectively, while total jet energy is plotted in black. Dashedand solid line styles refer to low and high resolution computations. Bottom panels: energy gain by the ambient normalized to the initialtotal jet energy. The thin vertical dotted line marks when fluid material begins to leave the lateral side of the computational domain. mation of fast magnetosonic disturbances propagating intothe ambient medium. At later stages these fronts steepeninto weak shock waves. This configuration makes the doublehelix pattern more difficult to form owing to the presenceof the velocity shear. However, similarly to Mizuno et al.(2011), a helical density structure stills persists surroundingthe central magnetic helix and propagates along the jet.For the super-Alfv´enic jet (case A10), the linear growthof the perturbation is accompanied by the formation ofstrong waves propagating obliquely away from the axis andlater steepening into shocks (see Section 3.2). The three-dimensional structure, rendered in the right-hand panel ofFig. 5 at t = 15, shows a large wavelength non-axisymmetricdeformation as well as the presence of small scale surfacemodes not seen in the previous two cases. This can be in-spected from Fig. 6 for 15 . t .
20. Note also, that duringthis phase, both low and high resolution computations showsimilar trends.
The transition from the linear to the nonlinear phase leadsto an overall change of morphology characterized by largescale energy and momentum redistribution. In general, theoriginal equilibrium configuration is destroyed after a fewAlfv´en crossing time following the end of the linear phasealthough the different growth of CD or KH modes deeply affects the evolutionary stages of the jets in the three casesconsidered here.
The energy initially carried by the jet is subject to consid-erable variation as a result of the instability and the ensu-ing interaction with the external medium. In the top panelsof Fig. 6, we plot the relative contributions of the volume-integrated kinetic, thermal and magnetic energies of the jet E X , jet ( t ) = Z E X T jet dx dy dz , (19)normalized to the initial total energy of the jet, E jet (0) = P X E X , jet (0). Here, X = { kin , mag , th } is a short-handnotation used to label the different contributions, i.e. E kin = ρ v / E mag = B / E th = p/ (Γ − E X , amb = E X , amb ( t ) − E X , amb (0) normalized to E jet (0).Here E X , amb is defined as in equation (19) with T jet → (1 − T jet ).Although the three contributions remain approximatelyconstant during the initial evolutionary stages, the end of thelinear phase is marked by a rapid transfer of energy from thejet to the ambient medium and, to a lesser degree, by a par-tial dissipation of magnetic and kinetic energies into heat. Inthe static column case, a fraction ≈ −
30% of the initial c (cid:13) , 1–14 inear and nonlinear evolution of MHD jets Figure 7.
Three-dimensional structure of the A10 jet (high resolution case) at three different instants t = 17 (left), t = 19 (middle) and t = 22 (right) showing magnetic pressure isosurfaces with density isosurfaces (top panels) and gas pressure (bottom panels). Steepeningof waves into shocks is evident at t ≈
19 in the pressure structure. jet magnetic energy is transferred to the ambient mediumthereby accelerating and ultimately heating the surround-ing gas. This transition occurs more violently when the jetvelocity increases and the beam becomes more kineticallydominated. KH-driven fast magnetosonic disturbances aresheared by the flow and later steepen forming shock wavesthat provide an efficient dissipation mechanism of mechani-cal energy. This situation is best depicted in Fig. 7 showingdensity, magnetic and thermal pressure isosurfaces at differ-ent times for the super-Alfv´enic jet which becomes forcefullydisrupted on a very rapid time scale around t ≈
20. Here thejet loses up to ≈ −
90% of its initial energy which becomesthen available to the ambient medium mostly in the form ofheat and, secondly, kinetic energy.The employment of high resolution introduces modestvariations during the onset of the nonlinear phase and re-sults, for sheared jets, in a somewhat more efficient heat gen-eration. Computations performed at higher numerical reso-lution yield a larger energy gain by the ambient medium and,most remarkably in Case A10, result in an earlier onset ofthe nonlinear phase. However, its effect is less recognizableat later times.Note that the total energy is not strictly conserved sincethe entropy equation is selectively used to evolve the internalenergy during the early stages ( t <
As demonstrated in Section 3.1, the different jet configura-tions become linearly unstable to the growth of the m = 1CDI mode with wavelength equal to the vertical box size.This induces large-scale helical displacements of the jet asalready shown in Fig. 5. In order to quantify the amount ofjet deformation and distortion we compute the barycentercoordinates as in Mignone et al. (2010, 2013)¯ x ( z, t ) = h x, χ i , ¯ y ( z, t ) = h y, χ i , (20)where χ = ρ T jet is the jet mass density.In Fig. 8 we plot the radial distance ¯ R ( z, t ) = p ¯ x ( z, t ) + ¯ y ( z, t ) as a function of the vertical coordinatefor the selected cases at different simulation times for highresolution (low-resolution computations behave similarly).For the static column (left-hand panel), the initial magneticsurface expands and stretches sideways while maintaininga simple helical structure with constant growing radius up c (cid:13) , 1–14 Anjiri et al.
Figure 8.
Radial distance of the jet barycenter from the axis as a function of the vertical distance. From left to right: A0, A1 andA10. Different colours and symbols refer to the simulation times described in the legend. Note that a perfect helical displacement isrepresented by a horizontal line. to t &
80 where ¯ R ≈ .
8. At later times, the helical de-formation stretches up ¯ R ≈ . t &
24, middle panel inFig. 8). Here, for t &
24, we see the development of an addi-tional perturbation mode with wavenumber k = 4 k where k = 2 π/L z is the initial perturbation wavenumber (see Sec-tion 2.3). This may be explained by the fact that the initialequilibrium for the A1 jet is liable not only to the CDI butalso to the KHI mode that has comparable growth rate (sec-ond panel of Fig. 2). This mode grows on top of the CDImode and it reaches a maximum amplitude at t ≈
80 whileit is later modified by nonlinear interactions. The same ef-fect is also found in the super-fast jet, where large amplitudeperturbations with smaller wavelength begin to develop for t &
20 (see the right-hand panel in Fig. 8).The final configuration approaching the saturated stateis shown in Fig. 9 where we display both the fluid density ρ (red colour) and jet density ρ T jet (green isosurfaces) in thehigh resolution simulations. At the centre of the domain, acentral region with larger magnetic energy and gas density isformed. This region, on the other hand, contains very little ofthe initial jet material which has been wrapped and twistedinside the helical magnetic flux tubes stretching sideways.These magnetic surfaces enclose a progressively larger vol-ume thereby leading to a substantial dilution of mass andmagnetic energy contained therein, see Fig. 9. The samebehaviour has been described in the force-free simulationcases presented by O’Neill et al. (2012) with the only dif-ference that, in our case, the nonlinear growth is not solelydominated by the m = 1 mode.A different evolution is observed in super-Alfv´enic jet(case A10), where the growth of instability goes along withthe onset of small-scale surface perturbations developing ontop of the m = 1 CDI mode. While the large-scale columndeformation takes place on a few Alfv´en time scale, these ad-ditional modes are triggered by the velocity shear and growas small ripples on the jet surface eventually dominating thesmall-scale structure and leading to a quick violent disrup-tion of the initial cylindrical configuration. This transition is mediated by the onset of KHI acting at the jet/ambientinterface and becoming very efficient in promoting entrain-ment, momentum transfer and mixing. Indeed, by t > Since the three components of momentum are conserved,it is instructive to compute how the total jet mass is re-distributed as a function of the velocity during the evolution.To this end we partition the velocity value range in smallintervals of width ∆ v and compute the mass-velocity densityfunction as δmδv ( v, t ) = 1∆ v hh Θ( v ) , ρ T jet ii , (21)where v ∈ [ v x , v y , v z ] is a velocity component while Θ( v ) =1 when the local zone velocity falls inside the given veloc-ity bin: | v i,j,k − v | < ∆ v/ v ) = 0 otherwise. Inother words, the numerator of equation (21) picks out thosecomputational cells having velocity between v − ∆ v/ v + ∆ v/ v is the width of the velocity bin. Notealso that P v ∆ vδm/δv = 1. Using the mass distributionfunction defined by equation (21), we compute the averageand variance as ¯ v ( t ) = X v v δmδv ( v, t )∆ v , (22) σ v ( t ) = X v ( v − ¯ v ) δmδv ( v, t )∆ v , (23)respectively.Fig. 10 shows the distributions of mass as a function of v ≡ v z at different times for Cases A0, A1 and A10 in thehigh resolution simulations. At t = 0 the distributions arestrongly peaked around the initial jet velocity with smalldispersion around the mean value meaning that most of thejet mass moves at the same speed. As the system evolves, thenet effect of the instabilities taking place inside the jet is thatof spreading and re-distributing the initial jet material and c (cid:13)000
20 (see the right-hand panel in Fig. 8).The final configuration approaching the saturated stateis shown in Fig. 9 where we display both the fluid density ρ (red colour) and jet density ρ T jet (green isosurfaces) in thehigh resolution simulations. At the centre of the domain, acentral region with larger magnetic energy and gas density isformed. This region, on the other hand, contains very little ofthe initial jet material which has been wrapped and twistedinside the helical magnetic flux tubes stretching sideways.These magnetic surfaces enclose a progressively larger vol-ume thereby leading to a substantial dilution of mass andmagnetic energy contained therein, see Fig. 9. The samebehaviour has been described in the force-free simulationcases presented by O’Neill et al. (2012) with the only dif-ference that, in our case, the nonlinear growth is not solelydominated by the m = 1 mode.A different evolution is observed in super-Alfv´enic jet(case A10), where the growth of instability goes along withthe onset of small-scale surface perturbations developing ontop of the m = 1 CDI mode. While the large-scale columndeformation takes place on a few Alfv´en time scale, these ad-ditional modes are triggered by the velocity shear and growas small ripples on the jet surface eventually dominating thesmall-scale structure and leading to a quick violent disrup-tion of the initial cylindrical configuration. This transition is mediated by the onset of KHI acting at the jet/ambientinterface and becoming very efficient in promoting entrain-ment, momentum transfer and mixing. Indeed, by t > Since the three components of momentum are conserved,it is instructive to compute how the total jet mass is re-distributed as a function of the velocity during the evolution.To this end we partition the velocity value range in smallintervals of width ∆ v and compute the mass-velocity densityfunction as δmδv ( v, t ) = 1∆ v hh Θ( v ) , ρ T jet ii , (21)where v ∈ [ v x , v y , v z ] is a velocity component while Θ( v ) =1 when the local zone velocity falls inside the given veloc-ity bin: | v i,j,k − v | < ∆ v/ v ) = 0 otherwise. Inother words, the numerator of equation (21) picks out thosecomputational cells having velocity between v − ∆ v/ v + ∆ v/ v is the width of the velocity bin. Notealso that P v ∆ vδm/δv = 1. Using the mass distributionfunction defined by equation (21), we compute the averageand variance as ¯ v ( t ) = X v v δmδv ( v, t )∆ v , (22) σ v ( t ) = X v ( v − ¯ v ) δmδv ( v, t )∆ v , (23)respectively.Fig. 10 shows the distributions of mass as a function of v ≡ v z at different times for Cases A0, A1 and A10 in thehigh resolution simulations. At t = 0 the distributions arestrongly peaked around the initial jet velocity with smalldispersion around the mean value meaning that most of thejet mass moves at the same speed. As the system evolves, thenet effect of the instabilities taking place inside the jet is thatof spreading and re-distributing the initial jet material and c (cid:13)000 , 1–14 inear and nonlinear evolution of MHD jets Figure 9.
Three dimensional renderings for Case A0 (left-hand panels), Case A1 (middle) and Case A10 (right-hand panels) at the endof the simulations. In the top panels, we display fluid density (sliced plane, red colour) and jet density ρ T jet (green isosurfaces). In thebottom panels, we show magnetic energy in the sliced planes (blue-green) while velocity vector field is given by the arrows. Figure 10.
Mass-velocity density distribution as function of the vertical velocity v z for the three different configurations. In each panel,we plot δm/δv at different times using the coloured solid lines shown in the legend. momentum over a much wider surface area with consequentreduction of the average jet velocity.For Case A0, the mass-velocity density function re-mains symmetric around the initial central value (¯ v z ≈ . . σ v z . . t . . t . c (cid:13) , 1–14 Anjiri et al.
Figure 11.
Relative contributions of the volume-integrated jetmomentum (solid lines) and ambient momentum (dashed lines)for case A1 (black) and A10 (red) as functions of time. tive within the range 10 . t .
40 and gives rise to a sequenceof non-symmetric and highly irregular distribution profiles(blue and green lines in Fig. 10). During the final steadystate the jet mass is again characterized by a symmetricdistribution peaked around ¯ v z ≈ .
012 with σ v z ≈ . t ≈
20, in particular, the jet averagevelocity has already halved and, by t ≈
30, most of the jetmass moves at a much smaller velocity (¯ v z < . v z ≈ .
086 and small dispersion(variance σ v z ≈ . q tot = ρv z .In Fig. 11, we plot ¯ q jet / ¯ q tot and ¯ q amb / ¯ q tot as a function oftime for Case A1 and A10 in the high resolution simulationswhere the volume-integrated jet and ambient momenta are,respectively, defined as¯ q jet ( t ) = hh ρv z T jet ii , ¯ q amb ( t ) = hh ρv z (1 − T jet ) ii , (24)while ¯ q tot = ¯ q jet + ¯ q amb . As expected, the end of the lin-ear phase is distinguished by a net transfer of momentumfrom the jet to the ambient medium. The time scale of thistransition closely reflects the changes observed in the mass-velocity distribution profiles and occurs faster in the super-fast jet case. Indeed, for Case A1 and t & t & In this work we have presented 3D numerical simulationsof magnetized jets evolving from an initial cylindrical equi-librium configuration described by uniform density, smallmagnetic pitch and highly magnetized plasma ( β ≈ − ).Three cases have been considered corresponding to a static,trans-Alfv´enic and a super-Alfv´enic jet with Alfv´enic Machnumber M A = 0 , ,
10, respectively. Simulations have beenperformed by perturbing the initial equilibrium state withthe exact eigenfunction of linear perturbative theory corre-sponding to one wavelength of the fastest-growing CDI modewith m = 1.Our results demonstrate that the predicted lineargrowth rate is well reproduced using a high order reconstruc-tion method and a grid resolution of at least ≈
10 zones perjet radius. Higher resolution may be needed in the case ofsuper-fast jets to reduce grid-induced numerical noise.Overall, the linear evolution is characterized by large-scale growing helical deformations of the plasma columntriggered by the excitation of the CDI mode. The insta-bility breaks the initial axial symmetry and develops on afew tens of Alfv´en crossing times while proceeding faster forthe super-fast jet. Density and magnetic field tend to forma double-helix pattern featuring regions of alternating en-hanced density and magnetic field. The magnetic flux tubestretches sideways and wrap regions of depleted mass andinternal energy. As the flow velocity is increased, the jet de-formation is accompanied by the propagation of strong fast-magnetosonic waves later steepening into shocks. Duringthis phase, the energy budget remains essentially constantand little energy and momentum exchange are observed.After the initial transient phase, the equilibrium is con-siderably altered and the final structure strongly depends onthe velocity shear layer which inevitably introduces couplingbetween CDI and KHI modes. When a velocity-shear is notpresent, the growth of the CDI mode result in large-scale he-lical deformations that do not lead to complete disruption.In this sense, our findings favourably compare to the force-free configuration of O’Neill et al. (2012) who accomplished3D simulations of local comoving jets. However, despite thefact that all magnetized columns under consideration areunstable to CDI modes, the presence of KH surface modesleads to the formation of small-scale distortions that mayeventually crumble and destroy the helical structure. Thiseffect becomes more pronounced at larger velocities and, inthe most severe case, leads to the complete jet disruptionand the formation of a chaotic turbulent flow on a veryshort time scale. These results are in contrast with the find-ings of Baty & Keppens (2002) who considered periodic jetconfigurations threaded by weaker magnetic fields than theones considered here and found that the presence of CDImodes provides a stabilizing nonlinear interaction mecha-nism weakening the disruptive effect of KHI perturbations.The onset of the nonlinear phase is marked by a net dis-sipation of energy into heat, a process that proceeds gradu-ally for slowly moving jets but becomes violently amplifiedby the formation of strong magnetosonic shocks for fast jets.Concurrently, as the system evolves towards the saturatedregime, the instabilities act so as to spread and re-distributethe initial jet material and momentum over a larger surfacearea with consequent reduction of the jet average velocity c (cid:13) , 1–14 inear and nonlinear evolution of MHD jets and hence favouring jet braking in few Alfv´en time scales.We have verified that, for configurations with non-vanishinginitial axial velocity, a large fraction ( & percent ) of theinitial jet momentum is transferred to the ambient medium.The re-distribution process is efficiently regulated by thepresence of KHI operating at the jet/ambient interface andcan thus be held responsible for promoting entrainment, mo-mentum transfer and mixing.Future extension of this work will enlarge this analysisto global jet simulations in both classical and relativisticregimes. ACKNOWLEDGEMENTS
A.M. wish to thank G. Mamatsashvili for very valuable sup-port on linear perturbative analysis. We acknowledge theCINECA Award no. HP10BCP4GU, 2013 for the availabil-ity of high performance computing resources and support.
REFERENCES
Appl S., & Camenzind M., 1992, A&A, 256, 354.Appl S., 1996, A&A, 314, 995.Appl S., Lery T., Baty H., 2000, A&A, 355, 818.Bahcall J. N., Kirhakos S., Schneider D. P., Davis R. J.,Muxlow, T. W. B., Garrington S. T., & Conway R. G.,1995, ApJ, 452, L91.Bateman G., 1980, MHD Instabilities, MIT Press, Cam-bridge.Baty H., Keppens R., 2002, ApJ, 580, 800.Biretta J. A., 1996, in Hardee P., Bridle A., Zensus J. A.,eds, ASP Conf. Ser. Vol. 100, Energy Transport in RadioGalaxies & Quasars. Astron. Soc. Pac., San Francisco, p.187Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883.Bodo G., Rosner R., Ferrari A., & Knobloch E., 1989, ApJ,341, 631Bodo G., Massaglia S., Ferrari A., & Trussoni E., 1994,A&A, 283, 655.Bodo, G., Massaglia, S., Rossi P., Rosner R., Malagoli A.,Ferrari Aet al., 1995, A&A, 303, 281.Bodo G., Rosner R., Ferrari A., & Knobloch E., 1996, ApJ,470, 797.Bodo G., Rossi P., Massaglia S., Ferrari A., Malagoli A.,Rosner R, 1998, A&A, 333, 1117.Bodo G., Mamatsashvili G., Rossi P., & Mignone A. , 2013,MNRAS, 434, 3030.Bonanno A., Urpin V., 2011, A&A, 525, A100.Carey C. S., & Sovinec C. R., 2009, ApJ, 699, 362.Cohn H., 1983, ApJ, 269, 500.Eichler D., 1993, ApJ, 419, 111.Ferrari A., Trussoni E., 1983, MNRAS, 205, 515.Ferrari A., Mignone A., Campigotto M., 2011, in BonannoA., Kosovichev A., eds, Proc. IAU Symp. 274, Advances inPlasma Astrophysics. Cambridge Univ. Press, Cambridge,p. 410.Hardee P. E., 2000, ApJ, 533, 176.Hardee P. E., 2006, in Hughes P. A., Bregman J. N., eds,AIP Conf. Proc. Vol.856, Relativistic Jets: The Common Physics of AGN, Microquasars, and Gamma-Ray Bursts.Am. Inst. Phys., New York, p. 57.Hardee P. E., 2011, in Romero G. E., Sunyaev R. A.,Belloni T., eds, Proc. IAU Symp. 275, Jets at all Scales.Cambridge University Press, Cambridge, p. 41.Hardee P. E., Cooper M. A., Norman M. L., Stone J. M.,1992, ApJ, 399, 478.Hardee P. E., Clarke D. A., Howell D. A., 1995, ApJ, 441,644.Hardee P. E., Clarke D. A., Rosen A., 1997, ApJ, 485, 533.Istomin Y. N., & Pariev V. I., 1994, MNRAS, 267, 629.Istomin Y. N., & Pariev V. I., 1996, MNRAS, 281, 1.Lery T., Frank A., 2000, ApJ, 533, 897.Lery T., Baty H., Appl S., 2000, A&A, 355, 120.Mignone A., Bodo G., Massaglia S., Matsakos T., TesileanuO., Zanni C., & Ferrari A., 2007, ApJS, 170, 228.Mignone A., Rossi P., Bodo G., Ferrari A., Massaglia S.,2010, MNRAS, 402, 7M.Mignone A., Striani E., Tavani M., & Ferrari A. 2013, MN-RAS, 436, 1102Mizuno Y., Lyubarsky Y., Nishikawa K. I., Hardee P. E.,2009, ApJ, 700, 684.Mizuno Y., Hardee P. E., Nishikawa K. I., 2011, ApJ, 734,19.Nakamura M., Meier D. L., 2004,in Bertin G., Farina D.,Pozzolieds R., eds, AIP Conf. Proc. Vol.703, Plasma in theLaboratory and in the Universe: New Insights and NewChallenges. Am. Inst. Phys., New York, p.308.Nakamura M., Li H., Li Sh., 2007, ApJ, 656, 721.Osmanov Z., Mignone A., Massaglia S., Bodo G., FerrariA., 2008, A&A, 490, 493.O’Neill S. M., Beckwith K., & Begelman M. C., 2012, MN-RAS, 422, 1436.Payne D. G., Cohn H., 1985, ApJ, 291, 655.Perucho M., Hanasz M., Mart´ı J. M., Sol H., 2004a, A&A,427, 415.Perucho M., Mart´ı J. M., Hanasz M., 2004b, A&A, 427,431.Perucho M., Mart´ı J. M., Hanasz M., 2005, A&A, 443, 863.Perucho M., Lobanov A. P., 2007, A&A, 469, L23.Perucho M., Hanasz M., Mart´ı J. M., & Miralles J. A.,2007, Phys.Rev.E, 75, 6312.Raga A.; Noriega-Crespo A., 1998, AJ, 116, 2943.Romanova M. M., Lovelace R. V. E., 1992, A&A, 262, 26.Rosado M., Raga A. C., & Arias L., 1999, AJ, 117, 462.Rosen A., Hardee P. E., Clarke D. A., Johnson A., 1999,ApJ, 510, 136.Rossi P., Bodo G., Massaglia S., Ferrari A., & Mignone A.,2004, Ap & SS, 293, 149.Rossi P., Mignone A., Bodo G., Massaglia S., & Ferrari A.,2008, A&A, 488, 795.Sikora M., Begelman M., Msejski G., & Lasota J.-P., 2005,ApJ, 625, 72.Spruit H. C., 1996, Astrophysical Journal, 2022S.Todo Y., Uchida Y., Sato T., & Rosner R., 1993, ApJ, 403,164.Turland B. D., Scheuer P. A. G., 1976, MNRAS, 176, 421.Walker R. C., Ly C., Junor W., Hardee P. J., 2008, J.Phys.: Conf. Ser., 131a, 2053, 2053.Walker R. C., Ly C., Junor W., Hardee P. J., 2009, inHagiwara Y., Fomalont E., Tsuboi M., Murata Y., eds,ASP Conf. Ser, Vol. 402, Approaching Micro-Arcsecond c (cid:13) , 1–14 Anjiri et al.
Resolution with VSOP-2: Astrophysics and Technologies.Astron. Soc. Pac., San Francisco, p. 227 227.Wanex Lucas F., 2005, Ap&SS., 298, 337. c (cid:13)000