Advances in Space Research 37 (2006) 1975–1978 www.elsevier.com/locate/asr
Numerical modeling of the pulsar wind interaction with ISM S.V. Bogovalov a
a,*
, V.M. Chechetkin
a,b
, A.V. Koldoba
a,b
, G.V. Ustyugova
b
Moscow Engineering Physics Institute (State University), Kashirskoje shosse 31, Moscow 115409, Russia b The Keldysh Institute of Applied Mathematics RAS, Miusskaya sq. 4, Moscow 125047, Russia Received 29 October 2004; received in revised form 24 October 2005; accepted 25 October 2005
Abstract Time dependent numerical simulation of relativistic wind interaction with interstellar medium was performed. The winds are ejected from magnetosphere of rotation powered pulsars. The particle flux in the winds is assumed to be isotropic. The energy flux is taken as strongly anisotropic in accordance with prediction of the MHD theory of the relativistic winds. The modeling has been performed for the wind magnetization in the range 3 · 103–101. The numerical solutions reproduce the most spectacular features observed in the central part of plerions: toroidal structure and jet-like features. Increase of the wind’s magnetization results in decrease of the size of the synchrotron nebula. 2006 Published by Elsevier Ltd on behalf of COSPAR. Keywords: Pulsars; Plerions; Relativistic plasma; Computer simulations
1. Introduction Rotation-powered pulsars (RPS) lose their rotational energy in the form of a wind. The wind carries out energy and angular momentum of the pulsars. Particles of the wind move relativistically together with being frozen into the wind magnetic field. They do not emit synchrotron radiation which could bring us direct information about the winds. Thus, the winds would be unobservable to us unless they interact with the ISM. The pulsar winds are terminated by a shock wave at the interaction with the ISM. The shock redistributes the energy of particles accelerating a part of them up to 115 eV (Aharonian and Atoyan, 1998) and randomizes their pitch angle. They start to emit synchrotron and inverse Compton radiation in the post-shock region forming a plerion (or synchrotron nebula). In the last few years, telescopes on the Hubble, Chandra, and XMM observatories bring us new observational data about the central part of plerions. Observations in X-rays (Brinkmann et al., 1985; Weisskopf, 2000) and in *
Corresponding author. Tel. +7 095 323 90 45. E-mail address:
[email protected] (S.V. Bogovalov).
0273-1177/$30 2006 Published by Elsevier Ltd on behalf of COSPAR. doi:10.1016/j.asr.2005.10.047
optics (Hester et al., 1995) discovered spectacular torus and jet-like structures in the central part of the Crab Nebula. The same structures were found by Chandra around Vela pulsar (Helfand et al., 2001; Pavlov et al., 2000, 2001). Along with rather close similarity between the Crab and Vela plerions they have differences. The bolometric Crab Nebula luminosity is equal to dozens of percentages of the pulsar rotational losses. The Vela Nebula is rather compact and has low luminosity compared with the rotational losses of the pulsar, 5 · 105 (Helfand et al., 2001). The toroidal and jet-like structures were found in several other plerions: around PSR 1509-58 (Kaspi et al., 2000), plerions in the supernova remnant G0.9+1 (Gaensler et al., 2001), and G54.1+0.3 (Lu and Wang, 2001). Apparently, formation of these structures is a rather general phenomenon for plerions. Integral characteristics of the Crab Nebula are well described by the shock model (Kennel and Coroniti, 1984; Rees and Gunn, 1974). This theory explains the spectra and luminosity of the Crab Nebula in the energetic band from optics up to TeV gamma-rays pretty well (Atoyan and Aharonian, 1996, de Jager, 1992). However, the physics of the interaction of the wind with the interstellar medium is essentially simplified in this theory. The wind
1976
S.V. Bogovalov et al. / Advances in Space Research 37 (2006) 1975–1978
is considered as isotropic. An attempt to study the elongation of the Crab Nebula along the rotational axis has been carried out by Begelman and Li (1992). However, the theory in its original form is not able to explain the observable structures in the center of the Crab Nebula. It has been shown in the series of works that the formation of the toroidal and jet-like structures in the central part of the plerion is provided by two factors – anisotropy of the energy flux in the pulsar winds and by the affect of the toroidal magnetic field on the dynamics of plasma in the post-shock region (Lyubarsky, 2002; Bogovalov and Khangoulian, 2002a,b; Khangoulian and Bogovalov, 2003). All the attempts to reproduce the plasma flow in the post-shock region analytically (Begelman and Li, 1992; Khangoulian and Bogovalov, 2003) can be applied only to some limiting cases. The realistic calculation of the plasma flow in the post-shock region can be done only numerically. Recently, several groups have performed numerical modeling of the post-shock flow which confirmed predictions regarding the formation of the toroidal and jet-like structures in the central part of plerions (Komissarov and Lyubarsky, 2003, 2004; Del Zanna et al., 2004). Here, we present basic results of the Russian team. 2. Pulsar winds Energy flux density in the wind is a sum of the electromagnetic energy flux 4pc E B and the kinetic energy flux c2mc2nv densities, where c is the Lorentz factor of the wind bulk motion, n and v are the plasma density in the comoving coordinate system and the velocity, and E and B are the electric and magnetic fields, respectively. The relativistic winds from pulsars flow radially. Their magnetic collimation is infinitesimally small (Beskin et al., 1998; Bogovalov, 2001). Toroidal magnetic field Bu generated at the rotation of the pulsar is proportional to Bpsinh, where h is the polar angle and Bp is the poloidal magnetic field (Bogovalov, 1999). Electric field in the wind is connected with the magnetic field through ideality condition E ¼ cv Bu provided that r Xc , where X is the angular velocity of the pulsar. Therefore, the electric field is proportional to Bpsin h as well. Thus, the Poynting flux appears proportional to B2p sin2 h. The ratio of the Poynting flux over the kinetic energy flux is close to 104 at the light cylinder of the Crab pulsar (Coroniti, 1990). Practically all the energy flux is concentrated in the electromagnetic field. According to Kennel and Coroniti (1984) the situation at the shock is reversed. The ratio of the Poynting flux over the kinetic energy flux becomes equal to r = 3 · 103 right upstream the shock (Kennel and Coroniti, 1984). The kinetic energy flux here is a sum of the initial kinetic energy flux plus the initial Poynting flux. Therefore, the Lorentz factor of the wind takes the form (Bogovalov and Khangoulian, 2002a) c ¼ c0 þ cm sin2 h;
ð1Þ
where c0 200 is the initial Lorentz factor of the wind near the light cylinder. According to (Daugherty and Harding, 1982) c0 200 or even less (Hibshmann and Arons, B20 2 6 7 2001). cm ¼ ðXrv 0 Þ 4pn0 mc in conventional mod2 c 10 –10 0 els of the pulsar winds. Here n0, B0 are the initial density and the magnetic field at the light cylinder. It was assumed for the simulations that the particle flux c0n0c and the poloidal magnetic field B0 are isotropic. cm is constant in the case under the consideration and c is proportional to sin2 h. 3. Method The computer simulation includes the solution of the system of equations under boundary conditions with some initial state. The equations of relativistic MHD flows were considered in detail by Lichnerowicz (1967) and Anile (1989). The numerical method for the solution of these equations is presented in Kodoba et al. (2002). Direct numerical simulation of the plasma flow takes place only in the plerion limited by the termination shock at the inner side and by the contact discontinuity at the outer side, where the plerion comes into contact with the shocked ISM material. Location of these boundaries is unknown a priori. The dynamics of the pulsar wind in the preshock zone and of the ISM motion are not simulated, and their influence on the plasma flow in the plerion is taken into account by an effective fashion, determining boundaries of the corresponding regions. The position and propagation of the termination shock is calculated by an approximate solution of a problem about decay of a discontinuity between the ultrarelativistic wind and the plasma of the plerion. For the outer boundary (contact discontinuity) two types of boundary conditions were used. In the first type (type I), the pressure pext was specified as the outer pressure and was defined from the initial spherically symmetric solution. Thus, the total pressure in the flow on the outer sides is prescribed to the external pressure, qu2 + P = Pext. In the second boundary condition (type II), the pressure of the external medium is considered equal to zero and the external pressure is defined by the inertial pressure of the shell of the external medium. This pressure is proportional to the value qu2 at the boundary of the simulation box. These two types of boundary conditions on the outer boundary correspond to different types of confinement taking place in different plerions. The first one corresponds to the outer pressure confinement taking place, for example, in the Vela Nebula. The second one corresponds to the inertial confinement taking place in the Crab Nebula. Comparison of simulations for these two boundary conditions has shown that the flow in the central part appears essentially identical for both boundary conditions. Therefore, below we present the results only for the boundary condition of Type I.
S.V. Bogovalov et al. / Advances in Space Research 37 (2006) 1975–1978
1977 0.938 0.891 0.844
0.928
6
0.882 0.835 0.789 0.743 0.696 0.650
5
4
0.798 0.751 0.704 0.657 0.610
3
0.604 0.557
z
0.511 0.465 0.418 0.372 0.325
3
2
0.279 0.233
0.564 0.517 0.470
z
4
2
0.423 0.376 0.330 0.283 0.236 0.189
1
0.186
1
0.143
0.140 0.094
0.096 0.049
0.047
0
2
r
4
6
0
2
4
r
Fig. 1. The velocity distribution in the post-shock region for r = 3 · 103 (left panel), and for r = 101 (right panel). Solid lines – flow lines, color represents v/c. Pulsar is located at the left lower corner. Coordinates are scaled by the shock radius at the equator.
The plasma flow is considered as axisymmetric. There is symmetry in relation to the equatorial plane as well. Therefore, all the computational domain is located in the quarter of the total flow region. A spherically symmetric hydrodynamical solution (Kennel and Coroniti, 1984) is taken as the initial state in all cases. 4. Results We performed simulation for the pulsar winds with different magnetization starting with r = 3 · 103 and finishing with r = 0.1. Basically, our results agree well with the numerical simulations performed by Komissarov and Lyubarsky (2004) and Del Zanna et al. (2004). In Fig. 1, we present the distribution of velocity in units of the light velocity in the post-shock region. The distribution in Fig. 1(left panel) is produced by the wind with r = 3 · 103. White space in this figure is filled by the relativistic wind from the pulsar which is terminated by the shock front. The terminating shock front is already not sphere like it was assumed in the (Kennel and Coroniti, 1984) theory. In 3D space, the terminating shock takes the form of a torus. In the post-shock region, the shocked flow forms a toroidal structure in 3D space due to energy anisotropy of the pulsar wind. The jets are formed due to magnetic collimation of the post-shock flow. It is remarkable that all the flow from the pulsar is finally collimated along the rotational axis. Although it seems to contradict to observations, it is necessary to keep in mind that we observe apparently only a small and the brightest part of the jets. In addition to the terminating shock front, the secondary shocks are formed in the region of the flow which remains supersonic after passage of the terminating shock. Surprisingly, that the velocity of the plasma at the distance >2rsh is close to 0.5c, the value agrees well with that derived
from observation of the Crab Nebula. This is in strong contrast with the (Kennel and Coroniti, 1984) theory which predicts that downstream the shock the velocity of plasma goes down from the value 0.33c to zero with distance. In Fig. 1 (right panel), we present the results of simulation of the post-shock flow produced by the wind with r = 0.1. The increase of the wind magnetization results in the formation of more compact nebula and more effective collimation of the flow. This behavior agrees well with the observation of the compact nebula around Vela pulsar. This witnesses in favor of a higher magnetization of the wind from Vela. One of the most unusual results is the motion of the plasma towards the pulsar along the jet. This happens due to our assumption about constancy of the wind magnetization. The increase of wind magnetization to the axis in comparison with the equatorial region could eliminate this discrepancy. 5. Conclusion Our simulations show that the combined effect of energy anisotropy in the pulsar winds together with the wind magnetization results in the formation of toroidal and jet-like structures in the central part of the plerions. Increase of magnetization of the winds naturally results in decrease of the size of the nebula which could explain the fact that some pulsars are surrounded by rather compact nebula like it happens in the case of the Vela pulsar. Along with the general agreement between observations and calculation there are discrepancies. Nevertheless, it is evident that the detailed comparison with observations and final conclusions about wind magnetization could be possible after an accurate calculation of the synchrotron brightness distribution of the plerions for a wide set of magnetization parameters.
1978
S.V. Bogovalov et al. / Advances in Space Research 37 (2006) 1975–1978
Acknowledgments The work is performed under support of collaborative INTAS-ESA Grant N 120-99 and RFBR Grants N 0302-17089 and N 03-02-16548. References Anile, A.M. Relativistic Fluids and Magnetofluids. Cambridge Univ. Press, Cambridge, 1989. Atoyan, A., Aharonian, F.A. On the mechanism of gamma radiation in the Crab Nebula. MNRAS 278, 525–541, 1996. Aharonian, F.A., Atoyan, A. Nonthermal radiation of the Crab Nebula, in: Proceeding of International Conference on Neutron Stars and Pulsars. Univ. Acad Press Inc., p. 439, 1998. Begelman, M.C., Li, Z.-Y. An axisymmetric magnetohydrodynamic model for the Crab pulsar wind bubble. ApJ 397, 187–195, 1992. Beskin, V.S., Kuznetsova, I.V., Rafikov, R.R. On the MHD effects on the force-free monopole outflow. MNRAS 299, 341–348, 1998. Bogovalov, S.V. On the physics of cold MHD winds from oblique rotators. A& A 349, 1017–1026, 1999. Bogovalov, S.V. Acceleration and collimation of relativistic plasmas ejected by fast rotators. A& A 371, 1155–1168, 2001. Bogovalov, S.V., Khangoulian, D.V. The Crab Nebula: interpretation of Chandra observations. Astronomy Lett. 28, 373–385, 2002a. Bogovalov, S.V., Khangoulian, D.V. On the origin of the torus and jetlike structures in the centre of the Crab Nebula. MNRAS 336, L53– L55, 2002b. Brinkmann, W., Aschenbach, B., Langmeier, A. Xray morphology of the Crab nebula. Nature 313, 662–664, 1985. Coroniti, F.V. Magnetically striped relativistic magnetohydrodynamic winds – The Crab Nebula revisited. ApJ 349, 538–545, 1990. Del Zanna, L., Amato, E., Bucciantini, N. Axially symmetric relativistic MHD simulations of Pulsar Wind Nebulae in Supernova Remnants. On the origin of torus and jet-like features. A & A 421, 1063–1073, 2004. de Jager, O.C., Harding, A.K. The expected high-energy to ultra-highenergy gamma-ray spectrum of the Crab Nebula. ApJ 396, 161–172, 1992. Daugherty, J.K., Harding, A.K. Electromagnetic cascades in pulsars. ApJ 252, 337–347, 1982.
Gaensler, B.M., Pivovarof, M.J., Garmire, G.P. Chandra Observations of the Pulsar Wind Nebula in Supernova Remnant G0.9+0.1. ApJ 556, L107–L111, 2001. Helfand, D.J., Gotthelf, E.V., Halpern, J.P. Vela Pulsar and Its Synchrotron Nebula. ApJ 556, 380–391, 2001. Hester, J.J., Scowen, P.A., Sankrit, R., et al. WFPC2 Studies of the Crab Nebula. I. HST and ROSAT Imaging of the Synchrotron Nebula. ApJ 448, 240–263, 1995. Hibshmann, J.A., Arons, J. Pair production multiplicities in rotationpowered pulsars. ApJ 560, 871–884, 2001. Kaspi, V.M., Gaensler, B.M., Arons, J., Pivovaroff, M.J., Kawai, N., Tamura, K. Chandra observations of the Young, Energetic Pulsar PSR B1509-582001. AAS 197, 8312, 2000. Kennel, C.F., Coroniti, F.V. Confinement of the Crab pulsar’s wind by its supernova remnant. ApJ 283, 694–709, 1984. Khangoulian, D.V., Bogovalov, S.V. The role of a magnetic field in the formation of Jet-like features in the Crab Nebula. Astronomy Lett. 29, 495–501, 2003. Kodoba, A.V., Kuznetsov, O.A., Ustyugova, G.V. An approximate Riemann solver for relativistic magnetohydrodynamics. MNRAS 333, 932–942, 2002. Komissarov, S., Lyubarsky, Y.E. The origin of peculiar jet-torus structure in the Crab nebula. MNRAS 344, L93–L96, 2003. Komissarov, S., Lyubarsky, Y.E. Synchrotron nebulae created by anisotropic magnetized pulsar winds. MNRAS 349, 779–792, 2004. Lichnerowicz, A. Relativistic Hydrodynamics and Magnetohydrodynamics. Benjamin Press, New York, 1967. Lyubarsky, Y.E. On the structure of the inner Crab Nebula. MNRAS 329, L34–L36, 2002. Lu, F., Wang, D. Chandra Observation of SNR G54.1+0.3 – a Close Cousin of the Crab Nebula. AAS 33, 1491, 2001. Pavlov, G.G., Sanwal, D., Garmire, G.P., Zavlin, V.E., Burwitz, V., Dodson, R.G. Observations of the Vela Pulsar and its Compact Nebula with the Chandra high resolution camera. AAS 32, 733, 2000. Pavlov, G.G., Zavlin, V.E., Sanwal, D., Burwitz, V., Garmire, G.P. The X-Ray spectrum of the Vela Pulsar resolved with the Chandra X-Ray Observatory. ApJ 552, L129–L133, 2001. Rees, M.J., Gunn, J.E. The origin of the magnetic field and relativistic particles in the Crab Nebula. MNRAS 167, 1–12, 1974. Weisskopf, M.C. et al. Discovery of spatial and spectral structure in the X-ray emission from the Crab Nebula. ApJ 536, L81–L84, 2000.