Mass-loading of bow shock pulsar wind nebulae
aa r X i v : . [ a s t r o - ph . H E ] S e p Mon. Not. R. Astron. Soc. , 000–000 (0000) Printed 9 November 2018 (MN L A TEX style file v2.2)
Mass-loading of bow shock pulsar wind nebulae
G. Morlino, , ⋆ M. Lyutikov, † and M. Vorster ‡ Department of Physics, Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47907-2036, USA. INFN – Gran Sasso Science Institute, viale F. Crispi 7, 67100 L’Aquila, Italy.
Accepted 2015 September 21. Received 2015 September 21; in original form 2015 May 6.
ABSTRACT
We investigate the dynamics of bow shock nebulae created by pulsars moving super-sonically through a partially ionized interstellar medium. A fraction of interstellarneutral hydrogen atoms penetrating into the tail region of a pulsar wind will undergophoto-ionization due to the UV light emitted by the nebula, with the resulting massloading dramatically changing the flow dynamics of the light leptonic pulsar wind. Us-ing a quasi 1-D hydrodynamic model of both non-relativistic and relativistic flow, andfocusing on scales much larger than the stand-off distance, we find that if a relativelysmall density of neutral hydrogen, as low as 10 − cm − , penetrate inside the pulsarwind, this is sufficient to strongly affect the tail flow. Mass loading leads to the fastexpansion of the pulsar wind tail, making the tail flow intrinsically non-stationary.The shapes predicted for the bow shock nebulae compare well with observations, bothin H α and X-rays. Key words: stars: pulsars: general – stars: winds, outflow – ISM: general
It has been estimated that between 10% and 50% ofpulsars are born with kick velocities V NS &
500 km s − (Cordes & Chernoff 1998; Arzoumanian et al. 2002). Thesepulsars will escape from their associated supernova remnantsinto the cooler, external interstellar medium (ISM) in lessthan 20 kyr (Arzoumanian et al. 2002). As this time scale issufficiently short, the pulsars are still capable of producingpowerful relativistic winds. Furthermore, comparison withtypical sound speeds in the ISM, c s, ISM = 10 −
100 km s − ,shows that the pulsars are moving with highly supersonicvelocities. The interaction of the pulsar’s wind with theISM produces a bow shock nebula with an extended tail.If a pulsar is moving through a partially ionized medium,the bow shock nebula can be detected by the characteristicH α emission resulting from the collisional and/or charge-exchange excitation of neutral hydrogen atoms in the post-shock flows and the subsequent emission via bound-boundtransitions (Chevalier et al. 1980). To date, nine such bowshock nebulae have been discovered, including three around γ -ray pulsars (Brownsberger & Romani 2014).Hydrodynamic (and hydromagnetic) models ( e.g. Bucciantini 2002a) of bow shock nebulae predict the forma-tion of a smooth two shock structure schematically shownin Fig. 1: a forward shock in the ISM separated by a con- ⋆ E-mail: [email protected] † E-mail: [email protected] ‡ E-mail: [email protected] tact discontinuity from a termination shock in the pulsarwind. In the head of the nebula the shock and the contactdiscontinuity are situated at a distance given by Eq. (1),corresponding to the position where the ram pressure of theISM balances the pulsar wind pressure. The flow structurein the head of the nebula is reasonably well understood, es-pecially in the limit of strong shocks ( i.e. , when the pulsarvelocity is much larger than the ISM sound speed) and ne-glecting the internal structure of the shocked layer (Wilkin1996) .These models further predict that the pulsar wind ter-minates at a Mach disk located approximately at a distance d back = M NS d behind the pulsar (where M NS = V NS /c s, ISM is the Mach number of the pulsar moving through theISM). At approximately the same distance the oblique for-ward shock turns into a Mach cone with an opening angle ∼ /M NS . For distances larger than d back the flow in thetail of the nebula is smooth and nearly cylindrical, althoughsome models predict the development of shear flow instabil-ities ( e.g. , Kelvin-Helmholtz instabilities).In contrast to these numerical models, H α , radio and X-rays observations show that the morphologies of bow shocknebulae are significantly more complicated. More specifi-cally, observations reveal that the tails of bow shock nebulaehave a highly irregular morphology, with Fig.2 showing four This model is very realistic when the system cools efficiently,otherwise the pressure of the shocked ISM needs to be taken intoaccount.c (cid:13)
G. Morlino, M. Lyutikov and M. J. Vorster NS A P v ISM = V ISM shoched ISMContact discontinuityUnshocked windShocked wind v wind P wind = P d d back D Figure 1.
Schematic illustration of a bow shock nebula prop-agating through fully ionized
ISM, as seen in the rest frame ofthe pulsar. The dot-dashed rectangle shows the region zoomed inFig.3. such examples. All these nebulae have a characteristic “head-and-shoulder” structure , with the smooth bow shock in thehead not evolving into a quasi-conical or quasi-cylindricalshape, but instead showing a sudden sideways expansion(s).Arguably the most famous example is the
Guitar nebula powered by the pulsar PSR B2224+65 (top left panel inFig.2). As the name suggests, this nebula has a guitar-likeshape with a bright head, a faint neck, and a body consistingof several larger bubbles.Although the morphologies of these nebulae vary fromsource to source, there are a number of common featureswhich, in our view, not only reflect the intrinsic dynamicalproperties of the flows , but which are also independent ofthe subtle details of both the pulsar winds ( e.g. , the relativeorientation of the velocity and spin axis) and of the ISM.We stress the fact that all bow shock nebulae show qualita-tively similar morphological features not expected from sim-ple fluid models. In the X-ray and radio bands the tails showhighly non-trivial morphologies with quasi-periodic varia-tions in the intensity ( e.g.
Kargaltsev & Pavlov 2008). Forexample, in the case of the Guitar nebula the tail shows quasi-periodic bubble-like structures (van Kerkwijk & Ingle2008).These peculiar tail shapes have been interpreted as theresult of density variations in the ISM (Romani et al. 1997;Vigelius et al. 2007). However, based on the following con-siderations, we find this explanation unsatisfactory: ( i ) alltails show similar morphological variations (see Fig.2); ( ii )a common characteristic of these bow shock nebulae is thatthey are all highly symmetric with respect to the directionof motion of the pulsar - this is not expected if variations aredue to the external medium; ( iii ) morphological features inH α , radio and X-rays are quasi-periodic - this is also not ex-pected from random ISM density variations. From these ob-servations we conclude that the peculiar morphological fea-tures result from the internal dynamics of the pulsar wind ,rather than through inhomogeneities in the ISM. van Kerkwijk & Ingle (2008) have also previously pro-posed that the morphology of the Guitar nebula couldbe explained by (unidentified) instabilities in the jet-likeflow of pulsar material away from the bow shock. Al-ternatively, Bucciantini & Bandiera (2001) and Bucciantini(2002b) have suggested that the mass loading of pulsar windnebulae may strongly affect their dynamics. These authorshave shown that a non-negligible fraction of neutral atomscan cross the shocked ISM behind the bow shock withoutundergoing any interaction, thereby enabling these atomsto propagate into the pulsar wind region. Once inside thewind, neutral hydrogen can be ionized by UV or X pho-tons emitted by the nebula, and possibly by collisions withrelativistic electrons and positrons, resulting in a net massloading of the wind.In order to study this scenario, Bucciantini & Bandiera(2001) and Bucciantini (2002b) extended the thin-layer ap-proximation used to model cometary nebulae (Bandiera1993; Wilkin 1996, 2000). The thin-layer approximation isconceptually analogous to a 1-D model as it neglects thethickness of the nebula, while all quantities depend only onthe distance from the apex. Despite the above-mentionedsimplifications, these models provide a good description ofthe head region of the nebulae in terms of shape, hydrogenpenetration length scale, and H α luminosity, as was laterconfirmed by more accurate 2-D axisymmetric simulations,both in the hydrodynamic (HD) regime (Bucciantini 2002a;Gaensler et al. 2004) and in the relativistic magnetohydro-dynamic (MHD) regime (Bucciantini et al. 2005). Using a 3-D model, Vigelius et al. (2007) was able to extend the studyof these systems by also taking into account either a non-uniform ambient medium, or the anisotropy of the pulsarwind energy flux. However, none of these models are able toexplain the peculiar morphology of the H α emission oftenobserved in the tail regions of bow shock nebulae.While the above-mentioned studies focused primarilyon the head of bow shock nebulae, the aim of the presentpaper is to investigate the effect of neutral hydrogen on thetail region of these nebulae . The question we would like toinvestigate is whether the mass loading of neutral hydro-gen in the pulsar wind can explain the peculiar morphologyobserved at H α , radio and X-ray energies. In order to fo-cus on the effect of mass loading on the evolution of bowshock nebulae, complications introduced by magnetic fieldpressure (and topology) are neglected in the present paper.These aspects are indeed necessary for a comprehensive andrealistic treatment of the problem, and will be the subjectof a future study.At this point the question arises as to whether one canuse observations of the heliopshere to understand the prob-lem formulated in the previous paragraph. Although massloading plays an important role in the dynamics of the so-lar wind (Baranov et al. 1971; Baranov 1990; Zank 1999),there are a number of key differences between this scenarioand the pulsar wind scenario. Firstly, the velocity of theSun through the ISM is, most likely, weakly sub-fast magne-tosonic (McComas et al. 2012), whereas the pulsar’s motionis highly supersonic; secondly, the pulsar wind is very light- composed of lepton pairs - and one would therefore expectmass loading to have a greater effect; thirdly, most of theresearch related to mass loading in the solar wind concen- c (cid:13) , 000–000 ass-loading of bow shock PWN Figure 2.
Montage of H α images of optical bow shocks associatedwith pulsar wind nebulae. Shown are J2224+65, the so-called Gui-tar nebula (Chatterjee & Cordes 2002), J0742-2822, J2030+4415,and J2124-3348 (Brownsberger & Romani 2014). trated on flow in the head region, while we are interested inthe large scale dynamics of the tail flows.The outline of the rest of the paper is as follows: in § § § §
4, whilein § § § § Fig.(1) shows the typical structure produced by a pulsar(or a normal star) moving supersonically through the ISM,as seen in the rest frame of the pulsar, when the effect ofneutrals is not taken into account. Such bow shock struc-tures are preferentially produced when the pulsar propa-gates through the warm phase of the ISM, whose typicaltemperature is T ∼ · -10 K. In this case the soundspeed is c s ∼
10 km/s, hence the pulsar’s Mach number is ≫
1. Conversely, the hot phase of the ISM has T ∼ K and c s ∼
100 km/s, implying a Mach number close toone. In the hot phase the ISM is totally ionised, while inthe worm one the presence of neutral Hydrogen cannot be neglected. In fact the typical density of the warm ISM is0.2-0.5 cm − and the ionized fraction is estimated to rangebetween 0.007-0.05 for the warm neutral medium (WNM),and 0.6-0.9 for the ionised neutral medium (WIM) (see, e.g.,Jean et al. 2009, Table 1).The distance, d , between the pulsar and the contactdiscontinuity (CD) (formed between the shocked ISM andthe shocked wind) is obtained by equating the wind pressurewith the bulk pressure of the ISM: d = (cid:18) L w πV ρ ISM c (cid:19) = 1 . × L / w, V − n − / , − cm , (1)where L w = 10 L w, erg s − is the pulsar luminosity, V NS = 300 V km s − is its peculiar velocity, and ρ ISM = m p n ISM is the density of the dragged component of the ambi-ent medium, expressed in units of n ISM = 0 . n ISM , − cm − .In order for mass loading to play a role in the dynamicevolution of bow shock nebulae, neutrals are required tocross ∆ (the distance between the bow shock and the contactdiscontinuity) and penetrate into the wind region. Hence theinteraction length inside the shocked ISM, λ , must be largerthan ∆. Chen et al. (1996) estimated that ∆ = 5 / d , avalue that was later confirmed by Bucciantini et al. (2005)using numerical simulations (note that ∆ is smaller than d ). At this stage it is important to note that different phys-ical processes will lead to the interaction of neutrals in theshocked ISM than in the pulsar wind. In the next sectionit will be shown that a significant amount of neutrals canundergo charge exchange (CE) with protons in the shockedISM and collisional ionization with electrons (specifically inthe head of the bow shock system), while the collisionalionization of neutrals by protons can be neglected when V NS .
300 km s − . Despite the neutrals undergoing CE andionization, it will further be shown that a dynamically im-portant fraction of neutrals can still penetrate into the pul-sar wind in their original state.The pulsar wind most likely consists of only electronsand positrons (Ruderman & Sutherland 1975), and CE istherefore not possible. Rather, ionization of neutrals canonly occur through photo-ionization or through the colli-sion with relativistic electrons and positrons. The formerprocess is discussed in § § § § c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster cient, D B = vr L /
3, as an upper limit for the perpendiculardiffusion, we can estimate the time needed for a thermalproton to cross the typical distance d , which is given by t cross = d /D B . Using B = 1 µ G, d = 10 cm and v = 100km s − , we get t cross ≈ yr. Hence the diffusion of ions istoo slow and can be neglected. An alternative way for ions topenetrate inside the wind is through shear flow instabilitiesthat can develop at the contact discontinuity. Nevertheless,the effects of instabilities will not be limited to the simpleinjection of ions, but will mix the ISM materials with thepulsar wind resulting in a complex tale structure and prob-ably to its disruption. In our knowledge this scenario hasnever been studied in details, and it is far beyond the aimof the present paper. After crossing the bow shock, neutral atoms can interactwith the shocked protons. The interaction length is given by λ = V NS X ion n ISM r c h σ ( v rel ) v rel i , (2)where X ion is the ionization fraction of the ISM, r c is thecompression ratio of the bow shock, σ ( v rel ) is the rele-vant cross section of the process under consideration, and h σv rel i is the collision rate averaged over the ion distribu-tion function. When the ion distribution is a Maxwellian, h σv rel i is well approximated (within 20%) by the expression(Zank et al. 1996; Blasi et al. 2012) h σv rel i ≈ σ ( U ∗ ) U ∗ , (3)where U ∗ = s π k B Tm p (4)is the average, relative speed between the incoming hydro-gen atom and ions ( T is the temperature of the shockedISM determined by the Rankine-Hugoniot jump conditions,assumed to be ≫ T ISM ). Using the fiducial values n ISM =0 . − and V NS = 300 km s − , together with an ioniza-tion fraction of 90% and r c = 4 (the typical value for strongshocks), leads to the following estimates for the mean freepaths λ ion ,p ≈ . × cm (5) λ ion ,e ≈ . × cm (6) λ CE ≈ . × cm (7)for the ionization due to collisions with protons, electronsand CE, respectively. Note that λ ion ,e has been calculatedunder the assumption that the electrons downstream of theshock equilibrate rapidly with protons, thereby acquiringthe same temperature. If this assumption does not hold, thecollisional length scale for ionization due to electrons canbecome much larger than the value reported in Eq. (6). Inaddition, the values (5)-(7) are to be taken as lower limitsas they are valid just ahead of the nebula, where the com-pression ratio and the temperature obtain their maximumvalues.From these estimates it follows that only a negligiblefraction of the neutral hydrogen will be collisionally ionizedby electrons, whereas a significant fraction of neutrals will undergo CE. The neutrals resulting from a CE event willhave a bulk speed and a temperature that are close to thatof the protons in the shocked ISM. This implies that thenewly formed neutrals tend to be dragged with the shockedprotons along a direction parallel to the contact disconti-nuity. Nevertheless, the CE process produces a diffusion ofneutrals in the nose of the nebula and it may still be possi-ble for the newly formed neutrals to enter the wind region,provided that their diffusion velocity perpendicular to thecontact discontinuity is of the same order or larger thantheir velocity parallel to the contact discontinuity. This is acomplication that will not be addressed in the present paperbut is essential to estimate the correct amount of neutralsthat can penetrate into the wind.Although a large number of neutrals will be lost dueto CE, a minimum fraction of neutrals proportional toexp[ − ∆ / ( λ CE + λ ion )] will cross the shocked ISM regionwithout suffering any interaction and will enter the pulsarwind in their original state. These neutrals will not influencethe wind structure until they are ionized, either through col-lisions with relativistic electrons and positrons, or throughphoto-ionization with photons emitted by the nebula or bythe pulsar. In § § There are three different sources of photons to account for:thermal and non-thermal radiation from the pulsar, andnon-thermal radiation from the nebula. The thermal ra-diation from the pulsar has been shown to be negligible(see Bucciantini & Bandiera 2001, Eq.(25) and discussionbelow). On the other hand, non-thermal emission from boththe pulsar and the nebula can play a role as the non-thermalpulsar luminosity is generally comparable to the luminosityof the nebula.When the radiation emitted by the nebula has a suf-ficient luminosity, incoming neutrals can be ionized beforecrossing the bow shock, and one would consequently not ex-pect any effect from mass loading. Conversely, the require-ment that hydrogen atoms are not fully ionized at the bowshock imposes an upper limit on the total ionizing luminos-ity. We derive this limit closely following a similar derivationgiven by van Kerkwijk & Kulkarni (2001).Assuming a spherical symmetry for the emission emit-ted by the nebula, the ionization fraction ξ + at a distance r from the center of the nebula is dξ + ( r ) dt = (1 − ξ + )¯ σ ph N ph e − τ ( r ) πr − ξ + n e α rec , (8)where N ph is the number of ionizing photons emitted by thenebula per unit time, and ¯ σ ph = N − R ∞ I σ ph ( ν ) n ph ( ν ) dν isthe photo-ionization cross section averaged over the photondistribution. The photo-ionization cross section is σ ph ( ν ) = 64 α − σ T ( I/hν ) / = 10 − ( hν/ Ryd) − / cm , (9)with σ T the Thompson cross section and I = 1 Ryd theionization potential. As σ ph decreases rapidly with photonenergy, the only relevant photons are those with energy & I . c (cid:13) , 000–000 ass-loading of bow shock PWN Table 1.
Tabulated quantities, in order of appearance: pulsar name, spin-down power, non-thermal X-ray emission from the pulsar, non-thermal X-ray flux from nebula, photon index of nebula X-ray emission, ISM density. References: (1) Brownsberger & Romani (2014), (2)Abdo et al. (2013), (3) Hui & Becker (2007), (4) Romani et al. (2010), (5) Stappers et al. (2003), (6) Hui & Becker (2006). The valuesof L X marked with an ∗ are calculated using the relation log L X , pwn = 1 .
51 log ˙ E − . Chandra observations.Pulsar ˙ E L X , pul L X , pwn Γ X , pwn n ISM
Reference[erg s − ] [erg s − ] [erg s − ] [cm − ]J0437–4715 0 .
55 2 . × . × ∗ – 0 .
21 1, 2J0742–2822 19 . < . × . × ∗ – 0 .
28 1, 2J1509–5850 68 . . × (1 . − . × . +0 . − . .
14 1, 2, 3J1741–2054 12 . . × . × ∗ . ± . .
44 1, 2, 4J1856–3754 3 × − – 4 . × ∗ – 0 .
05 1, 2J1959+2048 21 . . × . × . ± . .
02 1, 2, 5J2030+4415 2 .
90 2 . × . × ∗ – 2 .
69 1, 2J2124–3358 0 .
68 8 . × . ± . .
47 1, 2, 6J2225+6535 0 .
16 – 5 . × ∗ – 1 .
43 1, 2
Solving Eq.(8) in the rest frame of the pulsar where theplasma is moving in the positive direction along the x -axiswith velocity V NS , one may write dx = V NS dt . If the ex-tinction and recombination of protons are neglected, Eq.(8),written in cylindrical coordinates ( x, ρ ), simplifies to dξ + − ξ + = ¯ σ ph N ph πV NS dxρ + x . (10)The solution is straightforward and reads ξ + ( x ) = 1 − (1 − ξ + , ) exp (cid:26) − r ρ (cid:20) π − arctan (cid:18) xρ (cid:19)(cid:21)(cid:27) , (11)where ξ + , is the original ionization fraction of the ISMfar from the nebula and r ≡ ¯ σ ph N ph / (4 πV NS ) defines thetypical distance where ionization is effective. As pointedout by van Kerkwijk & Kulkarni (2001), Eq.(11) describesan ionization fraction which, for fixed impact parameter ρ ,smoothly increases as one approaches the nebula from theleft in Fig.(1), strongly increases for ρ ∼ r , before approach-ing the asymptotic value 1 − (1 − ξ + , ) e − πr /ρ as one movesaway from the pulsar towards the right.In order to estimate the value of r , it is useful to expressthe product ¯ σ ph · N ph using the X-ray luminosity of the neb-ula, L X , which we define as the luminosity in the Chandraenergy range, i.e. , between ǫ = 0 . ǫ = 8 keV.Assuming that the non-thermal emission from the nebularanging between UV and X-ray energies is a single powerlaw with a photon index Γ, the value of r can be written as r = ¯ σ ph N ph πV NS = 1 . · F (Γ) L X, V − cm , (12)where L X, is the X-ray luminosity in units of 10 erg s − and V is the pulsar speed in units of 300 km s − , while F is a dimensionless function that depends only on the photonindex F (Γ) = I / −
42Γ + 5 [ ǫ − / − Γ ] ǫ I [ ǫ − Γ ] ǫ ǫ . (13)The typical values of Γ inferred for PWNe detected in X-rays are Γ > . (see also, e.g. Kargaltsev & Pavlov 2008). It should be kept in mind that a value of Γ ∼ . It follows that F > × − , with F reaching unity whenΓ = 2 .
55. Eq.(12) shows that r can easily be smaller thanthe stand-off distance, d , implying that the majority of neu-trals reach the bow shock without having been ionized. Itis useful to express the condition r < d in terms of anupper limit for the X-ray efficiency of the nebula, defined as η X ≡ L X / L w : r < d ⇒ η X < − F (Γ) − L − / w, n − / , − . (14)Typical measured values for PWNe are η X ≈ − , althoughlarge deviations from this value are observed. For H α emit-ting bow shock nebulae it is known that hydrogen atomsare present at the position of the bow shock, therefore theinequality (14) has to be satisfied. In fact, for PWNe emit-ting both H α bow shock and X-ray radiation, the estimatedvalue of η X is closer to 10 − (see values reported in Table 2),and the inequality (14) is thus satisfied. The photo-ionization of atoms inside the wind determinesthe rate of mass loading. In order to estimate this rate, thevalue of the photon density inside the nebula is required.For simplicity it is assumed that the wind is a cylinder witha cross section πd and a length R , emitting radiation uni-formly. The photon number density inside the wind is thengiven by n ph , w = N ph h l i πd Rc . (15)The value of R can vary greatly from one bow shock nebulato another, but in general R ≫ d . The quantity h l i is themean length traversed by a photon before escaping the neb-ula. When the presence of neutrals is neglected, PWNe aretransparent to UV radiation because the e + /e − pair den-sity is very small and the Compton scattering mean free the particles responsible for the non-thermal X-rays are producedthrough Fermi acceleration.c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster path is greater than the size of the nebula . Hence for geo-metrical reasons we can assume h l i ∼ d . If we account forthe presence of neutrals inside the wind, the mean free pathof photons due to ionization is l mfp = 1 / ( σ ph n N ), where σ ph is given by Eq.(9) and n N is the neutral density insidethe wind. We do expect n N < n N, ISM ∼ . − , hence l mfp > cm. In conclusion, also when we account forthe presence of neutrals, the assumption h l i ∼ d remains agood estimate and the ratio R/ h l i is ≫
1. Using Eq. (15),the ionization length of neutrals inside the wind is estimatedto be λ ph = V NS n ph , w ¯ σ ph c = 3 . · R h l i (cid:16) η X − (cid:17) − V − n − , − F (Γ) − cm , (16)where the second equality has been obtained using Eqs. (1)and (12). One important comment is in order: the ionizationlength scale, λ ph , is ultimately not the length scale thatdetermines whether neutrals will influence the dynamics ofthe wind. λ ph is an estimate of the length scale required toionize the largest fraction of neutrals. The quantity morepertinent to the dynamics of the wind is the relative changein the density of the wind induced by mass loading. Sincethe pulsar wind is light, consisting of electron-positron pairs,only a small fraction of neutrals are required to be ionized inorder for these neutrals to become dynamically important.The more important parameter is therefore the length scalewhere the dynamics of the wind is dominated by the ionizedprotons rather than by the electron-positron pairs. Thus, wedefine λ ML , ph as the length scale where the loaded mass isequal to the initial mass of the wind, i.e. , λ ML , ph = ρ e ρ N V wind n ph , w ¯ σ ph c = n e m e n N m p V wind V NS λ ph . (17)A numerical estimate of λ ML , ph can be found using Eqs.(16)to evaluate λ ph , and estimating the electron density of thewind, n e , by equating the luminosity of the pulsar with theenergy flux in the shocked wind, i.e. , ˙ M ≡ L w / ( γc ) = n e m e πd c . The result reads λ ML , ph = 1 . · V wind c R h l i (cid:16) η X − (cid:17) − n − , − F (Γ) γ cm . (18)It is seen that λ ML , ph can be much smaller than d due tothe fact that its value is inversely proportional to the elec-tron Lorentz factor, whose typical value is γ ≈ − .Even for Lorentz factor as low as 10 , a value estimated fromthe modelling of some pulsars (see, e.g., Dubus et al. 2015), λ ML , ph remains smaller than d . It will be shown in § non-relativistic wind. On the other hand, for a rela-tivistic wind the definition of λ ML has to be modified as itbecomes necessary to take into account the relativistic iner-tia of the wind. We will discuss this issue in § This conclusion could be different only for very young PWNewhich have a larger e ± density (see Atoyan & Aharonian 1996,for details). HD model. Moreover, we will show that there is no signifi-cant difference between the relativistic and non-relativisticresults, apart from the different definition of λ ML . The collisional ionization length scale depends on the Bethecross section (Kim et al. 2000) σ Bethe ( γ ) = 4 πa α β (cid:2) M (cid:0) ln( γ β ) − β (cid:1) + C R (cid:3) , (19)where γ is the Lorentz factor of the electrons (and positrons)and β = (1 − γ − ) / . The two constants M and C R arerelated to the atomic form factors of the target, and areindependent of the incoming particle energy. For hydrogenatoms M = 0 . C R = 4 . γ = γ , furnishes the collisional ionization mean free pathin the pulsar wind: λ col = V NS n e σ Bethe ( γ ) c = 3 . × V n ISM , − γ cm . (20)As discussed in the previous section, the electron densityof the wind is estimated by equating the luminosity of thepulsar with the energy flux in the shocked wind, i.e. , ˙ M ≡L w / ( γc ) = n e m e πd c (thereby ensuring that the dynamicsof the pulsar wind is not dominated by the magnetic field).Using the same values for the parameters as those used inEq.(1) leads to the estimate n e ≈ · − γ − n ISM , − cm − .Note that in the relativistic regime ( γ ≫ σ Bethe increasesonly logarithmically with energy, and that increasing γ from10 to 10 implies a decrease in λ col of only ∼ λ ML , col as the typical mass-loading scale where the loaded mass density is equal to theinitial mass density of the wind, i.e. , λ ML , col = ρ e ρ N V wind n e σ Bethe ( γ ) c = 2 . × V wind c n N , − cm . (21)It is again stressed that the above expression is for a non-relativistic wind and that the correct expression for the en-thalpy should be used for a relativistic wind. Comparisonof Eq.(21) and Eq.(18) shows that photo-ionization is con-siderably more effective than collisional ionization, even forvery low luminosity nebulae. Collisional ionization can thusbe safely neglected in all realistic cases. There is a finite probability that neutrals can be ionizedinside the region of the unshocked wind (shaded region inFig.(1)). These ions could change the wind dynamics but,as we argue below, this process can be neglected. In fact, forparameters typical to bow shock nebulae, the Larmor radiusof these ions, calculated in the rest frame of the wind, iscomparable or larger than the size of the unshocked windregion.To show that this is, indeed, the case, an expression re-lating the termination shock radius and the wind luminosityis required. Assuming that the energy fluxes of the magneticfield and particles are similar, the wind luminosity can be c (cid:13) , 000–000 ass-loading of bow shock PWN external flow V V incomingneutrals Ρ N , V x Wind Ρ H x L u H x L P = P A = A c s = c s1 u = u < R = d Figure 3.
Sketch of the simplified model used to study the effectof mass loading in the tail of the bow shock nebula. This Figurerepresents a zoom in of the dot-dashed rectangle of Fig. (1). Notethat the presence of the bow shock is neglected. The dash-dottedarrows show the regions where the neutrals in scenario B pen-etrate into the wind from the side of the contact discontinuity.This is in contrast to scenario A where neutrals only penetrateinto the wind in a cylindrical region with cross section A . estimated as L w = B R t c , where R t is the radius of thetermination shock and B is the magnetic field strength atthe termination shock. In the rest frame of the wind thenewly formed ion have a bulk Lorentz factor equal to thewind Lorentz factor, γ w , hence the Larmor radius is: r L = γ w m p c eB = γ w m p c / R t √L w , (22)and the ratio between r L and the size of the free wind regioncan be expressed as r L R t = m p c / γ w √L w ≈ γ L − / w, . (23)The ionization of neutrals in the unshocked wind can pro-duce ions with r L > R t that will escape the free wind re-gion, suffering only a small deflection. The requirement forthe production of ions in the free wind regions is a high bulkLorentz factor and a low spin-down power L w . Consequentlywe neglect the role of these ions and we will only accountfor neutrals ionized in the shocked wind. To study the effect of mass loading in the tails of bow shocknebulae, the steady-state conservation equations for mass,momentum and energy are solved in a quasi 1-D approxi-mation. Fig. (3) shows a sketch of the model and representsan idealization of the region enclose in the dot-dashed rect-angle of Fig. (1). The quasi 1-D approximation implies thatthe transverse cross section, A , of the flow can change, butthat all the characteristic quantities of the wind, i.e. , thevelocity, u , density, ρ and pressure, P , are assumed to be afunction of the position x only. This approach further impliesthat any internal structures are neglected, in particular thefree wind region and the termination shock shown in Fig.(1). For the sake of simplicity the presence of the bow shock isalso neglected, but note that possible effects of mass load-ing on the shape of the bow shock will be discussed in § V , pressure, P , and density, ρ .Lastly, it is assumed that both the ISM and the wind arenon-relativistic with an adiabatic index γ g = 5 /
3. In thenext section we will develop the model for a wind with arelativistic temperature moving with a non relativistic bulkspeed. The latter case is the more realistic scenario for apulsar wind, whereas the former case applies more to stel-lar winds . We anticipate that the main difference betweenthese two cases lies solely in the different values used for theenthalpy. While this changes the dynamical length-scale sig-nificantly, it does not change the qualitative behavior of thesolution. Finally, note that the steady-state nature of thesystem is not guaranteed. In fact, we will show that scenar-ios exist where the steady-state assumption is most likelyviolated.The conservation equations for mass, momentum andenergy for a quasi 1-D system, written in the rest frame ofthe pulsar, are ∂ x [ ρuA ] = qA ′ , (24) ∂ x (cid:2) ρu A (cid:3) + A∂ x P = qA ′ V , (25) ∂ x (cid:20)(cid:18) ρu + γ g γ g − P (cid:19) uA (cid:21) = qA ′ V / . (26)Note that ρ = ρ e + ρ p is the total density of the wind andthat the quantities on the right-hand side of Eqs.(24)-(26)represent the mass, momentum and energy flux due to theionization of neutrals, respectively. The mass loading termis given by q = ˙ n ( m e + m p ) . (27)As we showed in § n is determined predominantly byphoto-ionization, allowing one to write˙ n = n N n ph ¯ σ ph c , (28)where n N is the local density of neutrals, n ph is the pho-ton density as seen in the rest frame of neutrals, and ¯ σ ph is the photo-ionization cross section averaged over the pho-ton distribution. For the sake of simplicity we assume thatboth n N and n ph are constant along x . The spatial indepen-dence of n N is a good approximation when x ≪ λ ph ( i.e. ,the photo-ionization length scale given in Eq.(16)), while itbecomes necessary to include the spatial dependence of n N when x & λ ph . This can be done with a simple change ofvariables as we show in Appendix A. Conversely, a correctevaluation of the spatial dependence of the photon density isnon-trivial and requires a detailed analysis of the specific ob-ject one wants to study. Here we avoid such a complicationby assuming a constant value.It is also assumed that the incoming neutrals have thesame temperature as the ISM ( ≈ K), and they can thusbe considered as being cold with respect to the wind in thenebula. As a result of this assumption, the momentum and For stellar winds it is necessary to also take into account thecharge-exchange that occurs between neutral hydrogen comingfrom the ISM and protons inside the wind.c (cid:13) , 000–000
G. Morlino, M. Lyutikov and M. J. Vorster energy injection terms on the right-hand side of Eq.(25) and(26) only contain the contribution due to the bulk motion.The variable A ′ represents the effective area crossed byneutrals. In § A ′ ≈ A = πd .However, in reality the situation is more complicated. Wealso noted that the CE process produces a diffusion of neu-trals in the shocked ISM and that this can lead to neutralspenetrating into the wind from the side of the contact dis-continuity. A rough estimate of this process, considering aconstant diffusion coefficient, would result in a mass load-ing that is proportional to A ( x ) / . Moreover, it will subse-quently be shown that the effect of mass loading is to expandthe cross section of the wind. When this happens the dis-tance between the bow shock and the contact discontinuityis reduced, thereby further increasing the probability of neu-trals penetrating from the side of the contact discontinuity,specifically in the tail region of the pulsar wind. This re-sults in a mass loading that is proportional to the total area A ( x ). A realistic solution therefore requires a description ofthe bow shock geometry and how the neutrals interact inthe shocked ISM. Such a complication necessitates the useof a 2-D simulation, and is beyond the scope of the presentwork. In this study the two opposite situations, A ′ = A (scenario A) and A ′ ∝ A (scenario B), are investigated aswe expect a realistic situation to be bracketed between thesetwo scenarios.A consequence of the steady-state assumption is thatthe pressure is required to be constant everywhere, i.e. , P = P and ∂ x P = 0, and it is thus possible to simplify Eqs.(25)and (26).As a first step to finding the solution, Eq. (24) can besubstituted into the right-hand side of Eqs. (25) and (26),leading to ∂ x [ ρuA ( u − V )] = 0 , (29) ∂ x (cid:20) ρuA (cid:0) u − V (cid:1) + γ g P uAγ g − (cid:21) = 0 . (30)These equations define two constants of the system whichcan be used to write down two expressions, one for ρ and theother for A , as functions of the velocity u only. Evaluatingthe constants of the system at the position x = x = 0,where the boundary conditions are defined as u = u , A = A and ρ = ρ , the following expression for the area isobtained: A ( u ) A = u u + ( u − u )( u − V )( γ g − M uu , (31)and the following expression for the density: ρ ( u ) ρ = (cid:20) u − V u − V + ( u − u )( u − V )( γ g − M u (cid:21) − , (32)where M = u /c s is the initial Mach number of the windand c s = p γ g P/ρ is the sound speed. Note that all quan-tities with the subscript refer to values at the boundary x . The next step requires finding an expression for the ve-locity. Differentiating Eq. (25) by parts, and using Eq. (24)leads to ∂ x u = − q A ′ ( u − V ) ρuA . (33) This differential equation can easily be integrated by partsfrom x = 0 to x , after the quantity ρuA has been sub-stituted using Eq.(29). The solution for scenario A, where A ′ = A , is straightforward and reads: u ( x ) = λu + xV λ + x , (34)where we have introduced the length scale λ = u ρ /q . If q is calculated using the photo-ionization cross section, λ corresponds exactly to the definition provided in Eq.(17).An important and noteworthy result is that the solution(34) does not depend on the initial Mach number of thewind. Eq.(34) further shows that u ( x ) decreases monoton-ically from u (0) = u to u ( ∞ ) = V , and that the typicallength scale for this transition is λu /V . The solution forscenario B ( i.e. , A ′ = A ) is slightly more complicated butcan be expressed in an implicit form as follows: xλ = F u − uu − V + F ln (cid:20) u − V u − V u + α ( u − u )2 u (cid:21) , (35)where F = 2 u V u + α ( u − V ) , (36) F = 2 u ( u − V )(2 u + α )[2 u + α ( u − V )] , (37)and α = M ( u − V )( γ g − . (38)The solution (35) also shows that u ( ∞ ) = V . This canbe readily understood: when the amount of mass loaded ismuch larger then the initial mass of the wind, the final stateof the wind will be the same as the initial state of the neu-trals. This is confirmed by the behaviour of c s , which is alsoa monotonic decreasing function of x , as well as by the result u → V ⇒ c s ( u ) →
0, which one expects from the assump-tion of cold neutrals. The continual process of mass loadingtherefore decelerates the wind from u to V . The main dif-ference between solutions (34) and (35) is the typical lengthscale where this transition occurs: in the former case it is λ A ∼ λu /V , whereas in the latter case it is λ B ∼ λ . Thisresult once again has a simple explanation: the amount ofmass loaded in scenario A is a factor A /A smaller comparedto the amount of mass loaded in scenario B. For x → ∞ thesolution (31) shows that the cross section of the flow in-creases asymptotically to the valuelim x →∞ A ( u ) = A u V (cid:18) γ g − M (cid:19) , (39)where this asymptotic value is the same for both scenariosA and B. As we expect the transition length scale to beproportional to the loaded mass, we have λ A ≃ λ B A ∞ /A ≃ λ B u /V .The solutions for u ( x ), M ( x ), and A ( x ) correspondingto scenario A are shown in Fig. (4), and the solutions corre-sponding to scenario B in Fig. (5). For the latter case threepanels are presented showing the results for three differentvalues of the Mach number, while for scenario A only thecase M = 1 is shown as the solution depends only weaklyon the value of M . All these solutions are obtained by usingthe benchmark values summarised in Table 2. These valueslead to an asymptotic expansion of A ∞ /A ≈ ∼ c (cid:13) , 000–000 ass-loading of bow shock PWN Scenario A H M = L A (cid:144) A ® u (cid:144) V ® H u L = u r (cid:144) u r M c s0 (cid:144) V - - x (cid:144) Λ Figure 4.
Structure of the wind when mass loading occurs inscenario A, as calculated using the non-relativistic model. Theplot shows, as a function of the position, the wind velocity dividedby V (solid line), the Mach number, M (dashed thin line), andthe expansion velocity u r normalised to u r , as given by Eq.(44)(dot-dashed line). The plot also shows the normalised area, A/A (dotted line) and the ISM sound speed normalized to V (thindot-dashed line). The initial wind velocity is u /V = 1000 andthe initial Mach number of the wind is M = 1. Note that thespatial coordinate, x , is normalised to the mass loading lengthscale λ . Comparison of Figs. (4) and (5) confirms that the tran-sition in scenario B occurs much faster than in scenario A.However, apart from the different scale lengths of the tran-sition, the two scenarios are very similar. Fig. (5) furthershows that when the initial Mach number increases from 0.3to 3, the transition occurs faster by a factor ∼ . e.g. Szeg¨o et al. 2000, §
2, for a review). More specifically, the solution in the 1-Dmodel generally predicts that mass loading in a supersonicflow leads to the deceleration and heating of the flow, whilemass loading in a subsonic flow causes acceleration and cool-ing. The reason why this behaviour is not observed in thepresent model is due to the fact that the flow in a quasi 1-Dmodel can expand in the radial direction while the pressureremains constant, and therefore no acceleration in the flowdirection is possible. On the other hand, if the radial expan-sion were limited for some reason ( e.g. due to hoop stressescaused by to a toroidal magnetic field), then an accelera-tion of the flow would be possible. This specific point willbe developed in more detail in a subsequent paper.
Scenario B H M = L A (cid:144) A u (cid:144) V F H u L = u r (cid:144) u r M c s0 (cid:144) V - - x (cid:144) Λ Scenario B H M = L A (cid:144) A u (cid:144) V F H u L = u r (cid:144) u r M c s0 (cid:144) V - - x (cid:144) Λ Scenario B H M = L A (cid:144) A u (cid:144) V F H u L = u r (cid:144) u r M c s0 (cid:144) V - - x (cid:144) Λ Figure 5.
The same as Fig.(4) but for scenario B. From topto bottom: solutions are calculated for the initial Mach numberequal to 3, 1 and 0.3, respectively, while the initial wind velocityis identical for all panels, i.e. , u /V = 1000. Table 2.
Summary of the parameter values used for scenarios Aand B. n ISM X ion T ISM V u M L w [cm − ] [ K ] [km s − ] [erg s − ]0.1 0.9 10 c c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster
One of the limitations of the present, quasi 1-D model is thatthe bow shock which separates the unshocked ISM from theshocked ISM is neglected. It is therefore useful to comparethe results of the present model with the expected structureof the bow shock when mass loading is absent. In Fig. (6)the expansion profiles for scenarios A and B presented in theprevious section are compared with the profile of an idealbow shock. The bow shock profile is calculated using thethin-shock approximation (see Wilkin 1996, Eq.(9)) close tothe head of the nebula. When the distance from the headbecomes large, the thin-shock approximation is no longervalid, and we therefore replace this solution with the Machcone, where the inclination angle, θ , between the bow shockand the pulsar velocity is such that sin θ = 1 /M NS . For thebow shock profile the benchmark values summarised in Ta-ble 2 are again used. For the mass loaded wind profiles threedifferent cases are plotted: scenario A with λ/d = 1, andscenario B with λ/d = 5 and λ/d = 30. Fig. (6) shows thatfor all mass loaded cases the pulsar wind in the tail of thenebula expands beyond the position of the unperturbed bowshock, and one would thus expect the geometry of the bowshock to be affected. However, while the expanding windprofile in scenario A closely follows the profile of the bowshock, the expansion of the wind profile in scenario B is muchfaster, and is capable of producing the head-shoulder shapeobserved in some H α bow shock nebulae. It should also benoted that the wind profile in scenario A never crosses thebow shock when λ & d , and one would therefore not expectsuch a configuration to affect the bow shock profile.We note that neglecting the primary bow shock ( i.e., ne-glecting the ram pressure of the ISM) leads to an incorrectprediction for the wind profile when this profile increasesbeyond the Mach cone. We show in Appendix B that whenthis occurs the formation of a secondary shock in the regionbetween the primary bow shock and the contact disconti-nuity is expected. Such a situation is schematically shownin Fig.(7). Behind the secondary shock the pressure will in-crease, resulting in a bending of the contact discontinuitytowards the axes of the wind. Consequently the rapid ex-pansion predicted in scenario B will probably be attenuatedby the presence of the secondary shock.An interesting consequence related to the presence ofthe secondary shock is that the increase in temperaturewill lead to an enhancement of Balmer emission just behindthis shock. However, the enhancement of the emission willstrongly depend on the inclination angle of the secondaryshock: if the inclination angle is close to the Mach cone, thetemperature and Balmer emission will increase only slightly,while for a larger inclination angle the temperature will in-crease significantly, resulting in a strong enhancement of theBalmer emission.Although the quasi 1-D model has limitations, it is nev-ertheless interesting to compare the expansion speed of thepulsar wind with the sound speed of the external ISM. Thisexpansion speed along the radial coordinate r , as seen by anobserver co-moving with the external medium, is u r ( x ) = drdt = V √ πA dAdx = V √ πA dAdu dudx , (40)where we have used A = πr and t = x/V . The last equality can be used to express u r in a compact form as a functionof u only. Using du/dx from Eq.(33) and calculating dA/du from Eq.(31), leads to u r = V r A πλ F ( u ) , (41)where the dimensionless function F ( u ) reads F ( u ) = A ′ A u ( u − V ) u ( u − V ) (cid:18) γ g − M u − V u (cid:19) . (42)The value of F at x = 0 is F ( u ) = u − V u (cid:18) γ g − M u − V u (cid:19) (43)for both scenarios A and B. In the limit M ≪ u ≫ V one has F ( u ) = 1, and consequently the value of u r atthe origin is u r = V r A πλ = V d λ . (44)By contrast, for very large distances one has lim x →∞ F ( u ) =0. The behavior of the radial expansion speed can be seen inFigs.(4) and (5) where the normalised speed u r ( x ) /u r = F ( u ) is shown for all plots. With a simple study of thefunction F ( u ) it is easy to demonstrate that u r is a mono-tonically decreasing function of x for scenario A, while forscenario B the function F ( u ) has a peak where u equals u peak = 5 V u (2 u + α )4 V α + u (2 u + α ) , (45)where α is given by Eq.(38). In the limit u ≫ V one has u peak → V , which corresponds to a value of the expansionspeed that is equal to u r, peak ≡ u r ( u peak ) = u r r u V √ (cid:18) γ g − M (cid:19) / . (46)In addition, always in the limit u ≫ V , the peak of theexpansion speed is located at x peak = λ M ( γ g − ×× (
12 + 2 ln " u (cid:0) M ( γ g − (cid:1) V . (47)Using parameter values that are typical for a bow shockPWN leads to u r, peak ≈ (10 − u r and x peak ≈ (2 − λ .Comparison of Fig. (5) and Fig. (6) shows that for λ . d the expansion speed in scenario B is larger than the soundspeed in the ISM when the wind crosses the unperturbedbow shock profile. In other words, when the pulsar windexpands beyond the unperturbed bow shock, the expansionis supersonic and one expects a strong modification of thebow shock. By contrast, for λ & d the expansion of thewind in the ISM is subsonic, and one therefore expects a lesspronounced deformation of the bow shock profile.From Fig. (5) it should also be noted that for λ . d , thevelocity u r can be larger than the sound speed in the windas well as the wind velocity u . When the former conditionis realised the stationary approach is no longer valid, whilethe second condition leads to a break down of the quasi 1-Dapproximation. As a result our model can no longer be used.Based on the above arguments one may speculate that c (cid:13) , 000–000 ass-loading of bow shock PWN Scen.B H Λ= d L Scen.B H Λ= d L Scen.A H Λ= d L Bow - shock mass loading region for scenario A (cid:144) d r (cid:144) d Figure 6.
Comparison between three mass loaded wind profiles,calculated using the non-relativistic model (solid lines for sce-nario B and dashed line for scenario A), and the profile of theunperturbed bow shock, calculated using the thin shock approx-imation (dot-dashed line) (Wilkin 1996). The results for scenarioB are shown for two different values of the mass loading distance, λ = 5 d (thin lines) and λ = 30 d (thick lines), while scenario Ais plotted only for λ = d . All cases have an initial Mach number M = 1. The shaded region, with radius r = d , shows the massloading region for scenario A. when λ decreases to a value close to d , the non-stationarityof the problem could give rise to a periodic structure of ex-panding bubbles, similar to those observed in the GuitarNebula. However, addressing this scenario using only ana-lytical models is a complicated matter, requiring the use offull 2-D numerical simulations.From the present investigation it is difficult to predictwhich of the scenarios, A or B, is the more realistic one.Based on the 1-D approach, one can state that when theexpansion occurs inside the bow shock, the behaviour is mostlikely well described by scenario A. In this scenario neutralsrepeatedly undergo charge exchange in the shocked ISM,acquiring the same bulk speed, while also flowing parallel tothe contact discontinuity between the shocked ISM and therelativistic wind. On the other hand, when the wind is closeto the position of the unperturbed bow shock (or crossesit) the effective distance between the new bow shock andthe contact discontinuity could be small enough to allow arelevant fraction of neutral particles to penetrate into thewind. If this is indeed the case, the expansion speed shouldincrease notably, reaching values close to the value predictedby scenario B. It should also be noted that far from thehead of the nebula the neutral fraction of the ISM shouldbe larger as the ionization due to the UV radiation emittedby the nebula is less effective, and that this should producea faster transition to scenario B. In this section we present the solution of a mass loaded windaccounting for the fact that the pulsar wind is relativistic,hence we will use a relativistic expression for the enthalpy.This is the only difference with respect to the non-relativistictreatment presented in §
3. We will show that the solution is NS Primary bow shockSecondary shock
Contactdiscontinuity
ISM TS Mach disk
Figure 7.
Sketch of the generation of a secondary shock insidethe first bow shock due to the mass-loading induced expansion ofthe wind. The shaded region behind the secondary shock showswhere one would expect an enhancement of the H α emission. very similar to the non-relativistic one, apart from a differentvalue of the transition length scale λ .The relativistic equations for the conservation of par-ticle number, energy and momentum, written in the restframe of the neutron star and for a 1-D system, are ∂ x [ n p,e uA ] = ˙ nA ′ , (48) ∂ x [ wγ w uA ] = qc γ A ′ , (49) ∂ x (cid:2) wu A (cid:3) + c A∂ x P = qc γ A ′ V . (50)Here n e,p is the numerical density of electrons (protons), w is the total wind enthalpy and γ w and γ are the Lorentzfactors of the wind and the neutrals, respectively. The quan-tities q and ˙ n are given by Eqs.(27) and (28), respectively.We notice that in a relativistic treatment it is more conve-nient to start with the conservation of the particle number,Eq.(48), rather than with the conservation of mass. For thereader’s convenience, the derivation of Eqs.(48)-(50) is re-ported in Appendix C.As was the case for the non-relativistic model, we againneglect the free wind region and the termination shock. Weconsider only the pulsar wind after the termination shockwhere it becomes marginally relativistic with a bulk Lorentzfactor γ w ≈
1. Furthermore, the neutrals are non-relativistic,hence γ ≈
1, while the steady-state assumption requiresthe pressure to be constant everywhere, i.e. , P = P and ∂ x P = 0. Using these assumptions, it is thus possible tosimplify the system (48)–(50) as follows: ∂ x [ ρ e uA ] = ˙ nm e A ′ ≈ q ( m e /m p ) A ′ , (51) ∂ x [ ρ p uA ] = ˙ nm p A ′ ≈ qA ′ , (52) ∂ x [ wuA ] = qc A ′ , (53) ∂ x (cid:2) wu A (cid:3) = qc A ′ V . (54)Note that in Eqs.(51) and (52) we neglect the contributionof the electrons to mass loading, and that the approximation q = ˙ n ( m e + m p ) ≈ ˙ nm p is used. In order to close the system(51)-(54), an expression for the enthalpy is required. It isgenerally believed that the pulsar wind predominantly con-sists of electron-position pairs with highly relativistic tem-peratures, and we assume that when mass loading occursthere is no energy transfer between electron and protons. Inother words, electrons and protons do not reach a thermalequilibrium, but evolve independently with different temper-atures. This implies that the electron gas is always highly c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster relativistic, hence the rest mass contribution to enthalpycan be neglected, i.e. , w e = ǫ e + P e = 4 P e = 4 P . On theother hand, protons are non-relativistic, hence their thermalenergy is always negligible with respect to their rest massenergy, and their enthalpy is w p = ρ p c . The total enthalpyof the wind is thus w = w p + w e = ρ p c + 4 P . (55)Using these simplifications allows one to obtain an analyti-cal solution for the wind dynamics by following a proceduresimilar to the one outlined for the non-relativistic case. Com-bining Eq.(53) and (54) furnishes a first constant of motion: uAw ( u − V ) = uA ( ρ p c + 4 P )( u − V ) = const . (56)Similarly, combining Eq.(52) and Eq.(53) furnishes a secondconstant of motion: uA ( w − ρ p c ) = 4 uAP = const . (57)Evaluating the constants at the initial position x = x ,where u = u and A = A , Eq.(57) gives the solution forthe cross section A as a function of the velocity: A ( u ) = u A /u , (58)while dividing Eq.(56) by Eq.(57) gives the solution for theproton density as a function of u : ρ p ( u ) = 4 P c u − uu − V . (59)The electron density can be obtained in a similar way usingEq.(51) rather than Eq.(52), giving the following result ρ e ( u ) = ρ e, + m e m p P c u − uu − V , (60)where ρ e, is the initial electron density. One noteworthypoint is that the solutions (58), (59) and (60) do not dependon the specific assumption regarding the injection cross sec-tion A ′ (this result is also identical to the non-relativisticcase). The solution for the velocity u ( x ) can be obtained byderiving Eq.(54) by parts, and using Eq.(53). This leads tothe expression ∂u∂x = − qc A ′ A u − V wu . (61)As was done for the non-relativistic case, we again dis-tinguish between scenario A ( A ′ = A ) and scenario B( A ′ = A ). For scenario B the integration of Eq.(61) leadsto the following implicit solution: xλ rel = uu − V − u u − V + ln (cid:20) u − V u − V (cid:21) , (62)where we have introduced the relativistic length scale λ rel = 4 P ( u − V ) qc ≃ P ρ N c u n ph ¯ σ ph c . (63)The second equality results from assuming V ≪ u ∼ c andusing the photo-ionization rate n ph ¯ σ ph c in the calculation ofthe mass loading rate, i.e. , q = m p n N n ph ¯ σ ph c . Comparisonof Eq.(63) with the definition of the mass loading lengthscale for the non-relativistic case, Eq.(17), shows that theyare identical, with the exception that the quantity 4 P /c replaces the electron mass density. In the steady-state modelthe internal pressure of the wind is equal to the externalpressure, hence 4 P corresponds to the specific enthalpy of the electron-positron wind plasma. Using Eq.(16) from § λ rel : λ rel ≈ . · T R h l i u c η − X, − V − n − N, − F (Γ) − cm , (64)where P = n ISM K B T and T = 10 T K have been used.Choosing realistic values for the parameters associated withbow shock nebulae emitting H α shows that λ rel can be aslarge as 10 cm, but we stress that this value can vary byorders of magnitude, essentially due to the fact that the val-ues of the neutral density inside the wind and the luminosityof the PWN tail are difficult to estimate. Nevertheless, forany scenario λ rel represents the lower limit for the typicalexpansion length scale. For example, in scenario A a largerexpansion length scale is predicted.Substituting A ′ = A in Eq.(61) leads to the followingequation ∂u∂x = − qc P ( u − V ) ( u − V ) u , (65)which, once integrated by parts, gives the solution xλ rel = u u − V − u u − V . (66)In this case the typical expansion length scale is λ rel , A = u V λ rel = 4 P ( u − V ) u qc V . (67)As u /V can be as large as 10 , the mass loading lengthscale in scenario A can reach 10 cm. In a realistic case oneexpects the transition to occur between λ rel and λ rel , A .Eq.(64) predicts that a density of neutrals as small as10 − cm − is sufficient to have λ rel comparable to the sizeof the nebula. Thus, for a relativistic pair-plasma wind, arelatively small amount of neutrals can strongly affect thetail flow.Fig.(8) shows the solutions for u and A correspondingto scenarios A and B, where the parameter values given inTable 2 have been used. Comparison of Fig.(8) with Figs.(4)and (5) shows that the relativistic and non-relativistic solu-tions are very similar.As was done for the non-relativistic case, the profile ofa relativistic wind undergoing mass loading is compared tothe typical profile of the unperturbed bow shock. Fig. (9)is analogous to Fig. (6), but the wind profile is now calcu-lated using the relativistic expression for scenarios A andB obtained in Eq.(58). One can see that the relativistic andnon-relativistic model give very similar profiles, with the im-portant difference that λ is replaced by λ rel . Additionally, asimilar radial expansion is predicted for the correspondingscenarios of the non-relativistic and relativistic cases. Start-ing from Eq.(40), using du/dx from Eq.(65), and calculating dA/dx from Eq.(58), leads to the following expression for theexpansion speed: u r = V s A πλ u / ( u − V ) u / (cid:18) uu (cid:19) ± , (68)where the upper and lower signs refer to scenarios A andB, respectively. Fig. (8) also shows the normalised velocity u r /u r , where u r = V (cid:18) A πλ (cid:19) / = d λ rel V , (69) c (cid:13) , 000–000 ass-loading of bow shock PWN Scenario A A (cid:144) A u (cid:144) V u r (cid:144) u r c s0 (cid:144) V - x (cid:144) Λ rel Scenario B A (cid:144) A u (cid:144) V u r (cid:144) u r c s0 (cid:144) V - x (cid:144) Λ rel Figure 8.
Structure of a pulsar wind when mass loading occurs,assuming that the wind consists of a hot relativistic electron-positron pair plasma. Top and bottom panels report results forscenario A and B, respectively, using the same values summarisedin Table 2. The various lines represent the wind speed divided by V (solid line), wind cross section (dot-dashed line) and expansionspeed (dotted line) both normalised to their initial values. Thesound speed of the ISM, c s , also normalised to V , is shown forcomparison (thin dot-dashed line). Scenario A H Λ rel = d L Scenario B H Λ rel = d L Scenario B H Λ rel = d L Bow - shock mass loading region for scenario A (cid:144) d r (cid:144) d Figure 9.
Same as Fig.(6), but for a relativistic model. Noticethat the dot-dotted line is the same as in Fig.(6) and representsthe bow shock profile in the thin shock approximation for non-relativistic wind. is the expansion speed at x = 0, with this speed beingthe same for both scenarios A and B. Similar to the non-relativistic case, u r ( x ) is also a monotonically decreasingfunction of x for scenario A, while for scenario B it has apeak at u = u peak = 5 V , corresponding to u r, peak ≡ u r ( u peak ) = u r r u V √ . (70)The peak is located at the position x peak , which, in the limit u ≫ V , is equal to x peak = λ rel (cid:18)
14 + ln (cid:20) u V (cid:21)(cid:19) . (71)It is interesting to note that this result is identical to thenon-relativistic case for the limit M →
0. In other words,all discussions presented for the non-relativistic solutions atthe end of § λ rel . d the expansion velocity can becomelarger than the wind velocity, thereby causing the quasi 1-Dapproximation to break down. In this section we provide an example of a comparison be-tween the predicted profile of the wind due to mass loadingand the observed shape of an H α bow shock. Among thepulsars showing an anomalous H α bow shock, we chose PSRJ0742-2822 for two reasons: not only has several quantitiesbeen measured with sufficient accuracy, but the pulsar isbelieved to move almost perpendicular to the line of sight,a fact which, to some extent, simplifies the comparison be-tween the model prediction and observation.We use the parameter values reported byBrownsberger & Romani (2014). The pulsar luminos-ity is L w = 1 . · erg s − , while the measured propermotion is 29 . − . The most recent distance estimategives 2 . ± . d = 2 kpc, Brownsberger & Romani (2014)estimated n ISM = 0 .
28 cm − for the total ISM density and V ⊥ = 275 km s − for the projected component of the pulsarvelocity. We assume that the velocity of the pulsar is almostperpendicular to the line of sight, hence V ≈ V ⊥ .Using the above values, the stand-off distance is esti-mated to be d = 3 . · cm. Note that this value iscompatible with the measured distance between the pulsarand the apex, measured as 1 . ′′ , corresponding to a phys-ical distance of 4 . · d/ (2kpc) cm. Using this value of d , the profile of an ideal bow shock ( i.e., no mass load-ing) is calculated according to the procedure outlined at thebeginning of §
4. The resulting solution is shown in Fig. 10(dot-dashed line), while the small dot indicates the positionof the pulsar. The unperturbed bow shock profile fits thehead of the nebula reasonably well, but fails to reproducethe “fan” structure emerging behind the bow shock head.For the mass loaded wind, profiles are calculated usingthe relativistic model developed in §
5, with Fig. 10 showingthe results for both scenarios A (dashed line) and B (solidline). In order to ensure that the wind profile in scenarioB crosses the bow shock at the precise location where the“fan” structure emerges, the profiles are calculated using c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster the value λ rel = 2 d . Note that the wind profiles start at theposition d back = 4 d behind the pulsar, which correspondsto the location of the backward termination shock accordingto numerical simulations (Bucciantini et al. 2005). As previ-ously discussed in § §
5, scenario A cannot account fora rapid expanding structure emerging from the unperturbedbow shock as the resulting profile is very smooth, while sce-nario B does predict such a structure. On the other hand,from Fig. 10 it can be seen that the expansion in scenario Bis larger than the one detected in H α .Our aim is not to exactly reproduce the shape of thebow shock as our model only predicts the shape of the rel-ativistic tail wind. The more important point that we wantto emphasise is that the location where the “fan” structureemerges is compatible with our predictions. In particular,one can ask whether the ratio λ rel /d = 2 is compatiblewith observations. This can be checked by estimating theneutral density inside the wind.Analysing the H α flux, Brownsberger & Romani (2014)found that the neutral fraction of the ISM is comparable to ∼
1. However, as the uncertainty is quite large, we use afiducial value of 0.5. It was shown in § V .
300 km s − , therefore one need only take into accountthe effect of photo-ionization. As was shown in § F X < · − erg cm − s − (Abdo et al. 2013), whichtranslates into an upper limit for the X-ray luminosity of L X < . · ( d/ erg s − . For the non-thermal spec-trum we assume a power law with an index Γ = 2. Further-more, using Eq.(11) allows one estimate the neutral densityat the location of the backward termination shock, d back .The value of this density averaged over the cross section ofthe wind is found to be ∼ . n ISM = 0 .
065 cm − . Finally,using this value in Eq.(64) and dividing the result by Eq.(1),one finds that λ rel /d > . R/ h l i . A value of λ rel /d = 2thus implies R/ h l i .
60, which is fully compatible with thewind geometry.
In this paper we have neglected the role of the magneticfield inside the wind of the nebula. Nevertheless, using theconservation of magnetic flux, we now show that generalconclusions can be drawn from the knowledge of the winddynamics only. We consider two different configurations: apurely poloidal and a purely toroidal magnetic field, withboth cases assuming that the initial magnetic field is dy-namically unimportant. In the poloidal case the flux conser-vation is applied through the cross section A , which leadsto ¯ B x A = const , where ¯ B x is the average poloidal magneticfield strength. Using the solution for the cross section fromEq.(58), one finds that ¯ B x ∝ u . As the fluid speed decreases,the poloidal magnetic field strength will also decreases, andone would not expect the wind dynamics to be affected.For the toroidal configuration the opposite situation oc-curs. The conservation of magnetic flux is applied througha poloidal surface ( i.e. , a surface parallel to the x direc-tion) which implies ¯ B φ u √ A = const , where ¯ B φ is the av-erage poloidal magnetic field inside the wind. Again using Figure 10.
Comparison of the observed H α bow shock producedby the PSR J0742-2822 (Brownsberger & Romani 2014) with thewind profile obtained using the relativistic calculations for sce-nario A (dotted line) and scenario B (solid line). The dot-dashedline shows the predicted position of the unperturbed bow shock.The small dot in the head of the nebula indicates the position ofthe pulsar. The profile of the wind starts at the location of thebackward termination shock, which is 4 d behind the pulsar. Eq.(58), one finds ¯ B φ ∝ u − / . Consequently the strengthof the toroidal component increases along x . As the toroidalmagnetic field exerts a hoop stress on the wind ∝ ¯ B φ , themagnetic pressure can easily become comparable to the in-ternal pressure, thereby reducing the expansion of the wind.These conclusions have implications for the non-thermalemission. The synchrotron emissivity is proportional to j syn ∝ n e ¯ B , where n e is the density of relativistic elec-trons and ¯ B is the average field. The density of relativisticelectrons is constant along x , and as a result the expansioninduced by mass loading will enhance j syn if the magneticfield is mainly toroidal. Conversely, for a poloidal configura-tion at the location of expansion the synchrotron emissionshould decrease.Remarkably, this prediction is compatible with the ob-servations of at least one object. The bow shock nebula pow-ered by the pulsar J1509-5850 has been detected both in X-rays and radio, and shows an anti-correlation in these twobands (Hui & Becker 2007; Kargaltsev et al. 2008): whilethe head of the nebula is mainly bright in X-rays, the bulkof the radio emission is observed from the tail of the wind ata distance roughly 4’ far away from the pulsar, at a locationwhere the X-ray emission becomes negligible. Moreover, ra-dio polarimetry measurements show that the magnetic fieldin the tail is mainly toroidal (Ng et al. 2010).In our model these observations can be explained byassuming that the enhancement of the radio emission oc-curs where the tail expands due to mass loading, while thedecrease in X-ray emission should be the consequence of ra-diative cooling, which is efficient for the X-ray emitting elec-trons but is negligible for the radio emitting electrons. Thispicture is confirmed by the fact that J1509-5850 has alsobeen observed in H α , showing an anomalous bow shock ex-pansion right at the location where the enhancement of the c (cid:13) , 000–000 ass-loading of bow shock PWN radio emission begins. A further intriguing detail is that thecomparison between radio and X-ray data suggests a signifi-cant deceleration of the flow when moving downstream, withthe speed decreasing by 1-2 orders of magnitude (Ng et al.2010). From our model it follows that this decreases can bethe result of mass loading. However, a more quantitativemodel is needed to better compare this picture with obser-vations.Noticeably, an opposite situation has been observedfor the Mouse nebula (J1747-2958), where radio polarime-try shows a poloidal magnetic field structure along the tail(Gaensler et al. 2004). In this case X-ray and radio emis-sion are both peaked in the head of the nebula and decreasesmoothly along the tail. The effect of mass loading for theMouse nebula is probably negligible, in fact the H α emissionis not detected, and this is probably due to its exceptionalluminosity able to photo-ionise the ISM around it. The X-ray luminosity measured by Chandra is L X = 5 · ergs − for a distance of 5 kpc (Gaensler et al. 2004), while thepulsar spin down luminosity is L w = 2 . × erg s − ,giving an X-ray efficiency η x = 0 .
02. Gaensler et al. (2004)also measured the photon index as 1 . < Γ X < .
5, and pro-vide an estimate of the ISM density of ≈ . − . Usingthese values we can estimate the upper limit for the X-rayefficiency from Eq. (14) to be 1 . · − . Because η X is muchlarger than this upper limit, we infer that the vast majorityof H atoms have been ionized before to reach the bow-shock. In this work we have shown that the structure of a bow shocknebula produced by a neutron star propagating through theISM can be significantly affected by the presence of neutralhydrogen in the ISM. In many cases a non-negligible fractionof neutral atoms can penetrate into the relativistic windwhere they will be ionized by UV photons emitted by thenebula. Once ionized, the new protons and electrons interactwith the wind, leading to a net mass loading of the wind inthe tail region of the nebula. More specifically, this massloading is important in the shocked part of the tail wind,while it is most likely negligible in the region of the free,unshocked wind.To investigate the effect of mass loading, we have de-veloped a steady-state hydrodynamic model using a quasi1-D approximation, with the study focusing on two specificscenarios. In the first a wind consisting of a non-relativisticplasma with an adiabatic index γ g = 5 / §
3) was in-vestigated, while the second situation investigated a windthat consists of a hot relativistic electron-positron plasma,as is the case for a pulsar wind (see § A ∞ /A = u /V , where u ∼ c is the initial ve-locity of the pulsar wind after the termination shock, and V is the speed of the neutron star measured in the restframe of the ISM, which is of the order of few hundredskm s − . The model therefore predicts A ∞ /A ≈ λ ML . We further showed that thislatter length scale is determined by the distance where theenthalpy of the loaded mass is comparable to the enthalpyof the wind (compare Eq. (17) vs. Eq. (63)).Considering a relativistic wind consisting of an electron-positron plasma, and using parameters values observed forpulsar bow shock nebulae, we showed that the mass loadingeffect could be responsible for the head-shoulder shape ob-served in many bow shock nebulae. Unfortunately it is noteasy to make a quantitative prediction for specific objects asthe amount of mass loading depends on two important pa-rameters, specifically the density of neutrals inside the windand the density of UV photons, both of which are difficultto estimate.We showed that a relatively small density of neutralsinside the wind (as small as 10 − cm − ) is sufficient to af-fect the wind, decelerating the tail flow and producing afast expansion in the transverse direction. For comparisonwe remember that the typical number density of neutralhydrogen in the warm interstellar medium (where H α bowshock nebulae are though to propagate) is ∼ . − . − (Jean et al. 2009). In order for the mass loading effect to benegligible, either the ISM should be completely ionized, orthe photo-ionization due to the nebula should be so effectivethat all neutrals are ionized before they reach the termina-tion shock in the tail. Neither of these two possibilities canbe applied to H α bow shock nebulae as the presence of H α lines is an indication that the photo-ionization is incapableof fully ionizing the ISM. If this were the case, it would notbe possible to observe the H α emission in the first place.We also made a visual comparison of the wind profilepredicted from our model with the H α bow shock observedaround PSR J0742-2822. We showed that a density of neu-trals in the wind comparable to ∼ .
06 cm − is required toexplain the presence and the location of the “fan” structureobserved behind the head of the nebula. Such a density iscompatible with the estimated neutral density of the localISM, taking into account the photo-ionization resulting fromthe UV radiation emitted by the nebula, which, in turn, wasestimated from the upper limit of the X-ray flux.Finally, even though the magnetic field is neglected inthe present model, we discussed the qualitative behaviour ofa magnetised wind flow that is being loaded with mass. Thepoloidal component of the field decreases with distances,thereby becoming dynamically unimportant. On the otherhand, the toroidal component, if present, will be amplifieddue to the compression of the wind flow. Such an ampli-fication will have two important consequences: the first isto limit the transverse expansion of the wind due to theincrease of the magnetic hoop stresses; the second is to pro-duce an enhancement of the synchrotron emission in thelocation where the expansion occurs. Remarkably, this pre-diction is confirmed by the observation of the bow shocknebula associated with the pulsar J1509-5850, where an en-hancement of the radio emission is observed far from thehead of the nebula, and at the same location where the H α emission shows an anomalous expansion of the bow-shock.The quasi 1-D approximation presented in this work hasthe advantage of having a clear and simple analytical solu-tion, thereby making it very usable. However, it also hasseveral limitations. Firstly, we neglected the presence of any c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster internal structure, as well as the bow shock. Moreover, thequasi 1-D approximation is based on the assumption thatthe transverse expansion velocity is much smaller than thewind velocity, a condition that can be violated when λ ML is of the order of, or smaller than the typical stand-off dis-tance of the nebula. In this case the quasi 1-D approximationbreaks down. Additionally, the transverse expansion speedcan be faster than the sound speed of the wind. When thisoccurs the steady-state condition is violated and the solutionpresented is no longer applicable. It is therefore necessary toimplement a time-dependent solution. This last situation isespecially interesting as it points toward to possibility ofhaving a tail wind with periodic expanding bubbles simi-lar to what is observed in the Guitar nebula. Hopefully allthis limitations can be overcome using time-dependent, 2-D simulations, which will form the topic of research of aforthcoming paper. ACKNOWLEDGMENTS
GM is grateful to Niccol´o Bucciantini for many fruitful dis-cussions on PWNe and relativistic winds. This work wasfunded through the NSF grant num. 1306672.
APPENDIX A: DEPLETION OF NEUTRALSINSIDE THE WIND
For the solutions presented in the paper we have assumedthat the density of neutrals inside the wind remains con-stant, i.e. , that the number of ionized atoms is negligiblewith respect to the total number of neutrals. This is trueonly on a length scale much smaller than the ionizationlength scale, λ ph . Here we show how to modify the model,should this assumption no longer be valid, by taking intoaccount the fact that the density of neutrals is no longerconstant, but decreases as a function of distance.When the depletion of neutrals becomes important, sce-nario B (where neutrals can also penetrate into the windregion through the side of the contact discontinuity) can nolonger be described using a 1-D model as the neutral densityinside the wind will depend on both x and r . The discussionwill therefore be limited to scenario A, where the neutralspenetrate only through the head of the nebula.It is useful to introduce a new spatial variable, ξ , definedas dξ = ¯ n N ( x ) dx, (A1)where ¯ n N ≡ n N ( x ) /n N is the numerical density of neu-trals normalised to its initial value. Substituting ξ for n N ( x ),Eqs.(24)-(26) (or Eqs.(48)-(50) for the relativistic case) canbe rewritten by replacing ∂ x → ∂ ξ , and using˙ n = n N n ph ¯ σ ph c (A2)rather than ˙ n on the right-hand side. This new set of equa-tions is formally identical to Eqs.(24)-(26), where a constantdensity is assumed, with the only exception that x is re-placed by ξ . The solutions presented in § § x is substituted by the dependence on ξ . To find the full newsolutions therefore only requires that one finds the function ξ ( x ), relating the new variable ξ to the physical distance x .This can be easily done using Eq.(A1) along with the equa-tion for the evolution of the neutral density, which, for thescenario A, is: V dn N dx = − ˙ n = − n N n ph ¯ σ ph c , (A3)leading to a simple exponential solution n N ( x ) = n N e − x/λ ph , (A4)where λ ph is defined in Eq.(16). Finally, substitutingEq.(A4) into Eq.(A1), and integrating one obtains ξ ( x ) = λ ph (cid:16) − e − x/λ ph (cid:17) . (A5)It is easy to check that when x ≪ λ ph one has ξ = x , andthat the solutions presented in § § APPENDIX B: FORMATION OF SECONDARYSHOCKS
One of the main results of this study is that the mass loadingof a pulsar wind leads to a sideways expansion of the tail.As the contact discontinuity between the tail flow of a PWNand the ISM acts as impenetrable wall for the ISM flow, asecondary shock can be created. As we demonstrate below,this secondary shock will appear when the contact disconti-nuity makes an angle with the flow that is larger than theMach cone.As a model problem, consider a flow with velocity, V ,and internal sound speed, c s , bounded by a wall with theparabolic profile, r = x /L . A shock will appear when thecharacteristics first intersect (Landau & Lifshitz 1959). Thecharacteristics originating at point { x , r } are r ( t ) = r + c s t,x ( t ) = x + ( c s + V ) t. (B1)Eliminating t leads to r = x /L + c s c s + V ( x − x ) . (B2)The first intersection occurs where ∂r/∂x = 0: x = 12(1 + M ) L,M = V /c s . (B3)At this point the angle that the wall makes with the flow istan θ = 1 / (1 + M ) . (B4)Thus, for a highly supersonic flow a shock forms when theangle of attack becomes larger than the Mach angle. APPENDIX C: HYDRODYNAMICAL MODELFOR RELATIVISTIC WIND
To derive Eqs. (48)-(50) we start by writing a relativisticformulation of the mass-loading problem, and then we spe-cialise our equations to describe a typical pulsar wind. Ithas been shown by Komissarov (1994) (see also Lyutikov2003) that the covariant relativistic equation for a perfect c (cid:13) , 000–000 ass-loading of bow shock PWN fluid with the inclusion of mass loading can be written inthe following form: ∂T νµ ∂x ν = qcτ µ , (C1)where T νµ = wu ν u µ + P g νµ (C2)is the energy momentum tensor, w = ǫ + P is the totalenthalpy, and u µ = γ w (1 , u /c ) is the four-velocity of theplasma, with γ w denoting the Lorentz factor of the wind.The quantity qcτ µ represents the mass loading term, where τ µ = γ (1 , V /c ) is the four-velocity of the neutrals movingwith a Lorentz factor γ = (1 − V /c ) − / , and q = X i ˙ n i γ i m i (C3)is the mass injected per unit time per unit volume calculatedin the frame moving with velocity V , i.e. , the rest frame ofthe neutrals. It is assumed that the incoming neutrals havethe same temperature as the ISM ( ≈ K), and they canthus be considered as being cold with respect to the wind inthe nebula, i.e , γ i = 1. As the rate of injected electrons andprotons is the same, ˙ n e = ˙ n p = ˙ n , the equations describingthe evolution of the electron and the proton densities aresimilar, i.e. , ∂∂x ν [ n e,p cu ν ] = ˙ n , (C4)where n e,p is the numerical density of electrons (protons),with an expression for ˙ n given by Eq.(28). As was donefor the non-relativistic case, we assume that both n N and n ph are constant along x . A method to obtain the solutionwhen the depletion of neutrals inside the wind is taken intoaccount is described in Appendix A.The next step is to rewrite the conservation equations(C1) and (C4) for a 1-D system, integrating over a small vol-ume with a cross section A and a length dx (see Komissarov1994). Additionally assuming that the system is in a steady-state, the conservation equations for the particles’ number,energy and momentum in the x direction become: ∂ x [ n p,e uA ] = ˙ nA ′ , (C5) ∂ x [ wγ w uA ] = qc γ A ′ , (C6) ∂ x (cid:2) wu A (cid:3) + c A∂ x P = qc γ A ′ V . (C7)Note that the non-relativistic equations (24)-(26) can be re-covered from Eqs.(C5)-(C7) using the first order approxi-mation for the Lorentz factors, i.e. , γ w ≈ u /c and γ ≈ V /c . REFERENCES
Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS,208, 17Arzoumanian, Z., Chernoff, D. F., & Cordes, J. M. 2002,ApJ, 568, 289Atoyan, A. M. & Aharonian, F. A. 1996, A&A suppl., 120,C453Bandiera, R. 1993, A&A, 276, 648Baranov, V. B. 1990, Space Science Reviews, 52, 89 Baranov, V. B., Krasnobaev, K. V., & Kilikovskii, A. G.1971, Soviet Physics - Doklady, 15, 791Blasi, P., Morlino, G., Bandiera, R., Amato, E., & Caprioli,D. 2012, ApJ, 755, 121Brownsberger, S. & Romani, R. W. 2014, ApJ, 784, 154Bucciantini, N. 2002a, A&A, 387, 1066Bucciantini, N. 2002b, A&A, 393, 629Bucciantini, N., Amato, E., & Del Zanna, L. 2005, A&A,434, 189Bucciantini, N. & Bandiera, R. 2001, A&A, 375, 1032Chatterjee, S. & Cordes, J. M. 2002, ApJ, 575, 407Chen, Y., Bandiera, R., & Wang, Z.-R. 1996, ApJ, 469, 715Chevalier, R. A., Kirshner, R. P., & Raymond, J. C. 1980,ApJ, 235, 186Cordes, J. M. & Chernoff, D. F. 1998, ApJ, 505, 315Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581,A27Gaensler, B. M., van der Swaluw, E., Camilo, F., et al.2004, ApJ, 616, 383Hales, C. A., Gaensler, B. M., Chatterjee, S., van derSwaluw, E., & Camilo, F. 2009, ApJ, 706, 1316Hui, C. Y. & Becker, W. 2006, A&A, 448, L13Hui, C. Y. & Becker, W. 2007, A&A, 470, 965Janssen, G. H. & Stappers, B. W. 2006, A&A, 457, 611Jean, P., Gillard, W., Marcowith, A., & Ferri`ere, K. 2009,A&A, 508, 1099Kargaltsev, O., Durant, M., Pavlov, G. G., & Garmire, G.2012, ApJS, 201, 37Kargaltsev, O., Misanovic, Z., Pavlov, G. G., Wong, J. A.,& Garmire, G. P. 2008, ApJ, 684, 542Kargaltsev, O. & Pavlov, G. G. 2008, in American Instituteof Physics Conference Series, Vol. 983, 40 Years of Pulsars:Millisecond Pulsars, Magnetars and More, ed. C. Bassa,Z. Wang, A. Cumming, & V. M. Kaspi, 171–185Kim, Y.-K., Santos, J. P., & Parente, F. 2000, Phys. Rev.A, 62, 052710Komissarov, S. S. 1994, MNRAS, 269, 394Landau, L. D. & Lifshitz, E. M. 1959, Fluid mechanics(Oxford: Pergamon Press, 1959)Lyutikov, M. 2003, MNRAS, 339, 623McComas, D. J., Alexashov, D., Bzowski, M., et al. 2012,Science, 336, 1291Ng, C.-Y., Gaensler, B. M., Chatterjee, S., & Johnston, S.2010, ApJ, 712, 596Romani, R. W., Cordes, J. M., & Yadigaroglu, I.-A. 1997,ApJL, 484, L137Romani, R. W., Shaw, M. S., Camilo, F., Cotter, G., &Sivakoff, G. R. 2010, ApJ, 724, 908Ruderman, M. A. & Sutherland, P. G. 1975, ApJ, 196, 51Stappers, B. W., Gaensler, B. M., Kaspi, V. M., van derKlis, M., & Lewin, W. H. G. 2003, Science, 299, 1372Szeg¨o, K., Glassmeier, K.-H., Bingham, R., et al. 2000,SSRv, 94, 429van Kerkwijk, M. H. & Ingle, A. 2008, ApJL, 683, L159van Kerkwijk, M. H. & Kulkarni, S. R. 2001, A&A, 380,221Vigelius, M., Melatos, A., Chatterjee, S., Gaensler, B. M.,& Ghavamian, P. 2007, MNRAS, 374, 793Wilkin, F. P. 1996, ApJL, 459, L31Wilkin, F. P. 2000, ApJ, 532, 400Zank, G. P. 1999, Space Science Reviews, 89, 413 c (cid:13) , 000–000 G. Morlino, M. Lyutikov and M. J. Vorster
Zank, G. P., Pauls, H. L., Williams, L. L., & Hall, D. T.1996, JGR, 101, 21639 c (cid:13)000