Emission from Millisecond Pulsar Binaries
CJT van der Merwe
orcid.org 0000-0002-6363-1139
Dissertation accepted in partial fulfilment of the requirements for
the degree
Master of Science in Astrophysical Sciences
at the
North-West University
Supervisor:
Prof C Venter
Co-Supervisor:
Dr AK Harding
Co-Supervisor:
Dr Z Wadiasingh
Graduation December 2020
30871964
I. Kant
“The Lord who made the earth, the Lord who forms it to establish it, the Lord is his name, says this: Call to me and I will answer you and tell you great and incomprehensible things you do not know.”
Jeremiah 33:2-3 (CSB)
“When I observe your heavens, the work of your fingers, the moon and the stars, which you set in place, what is a human being that you remember him, a son of man that you look after him?”
Abstract
Natural and Agricultural Sciences Centre for Space Research
Magister Scientae in Astrophysical Sciences
Orbitally-modulated X-ray and Gamma-ray Emission from Millisecond Pulsar Binaries
by C.J.T. van der Merwe
In the 11 years since its launch, the Fermi Large Area Telescope has yielded a large number of new detections. Among these new detections are about 30 black widow and redback systems. These sources are compact binaries in which a millisecond pulsar heats and may even ablate its low-mass companion by its intense wind of relativistic particles and radiation. In such systems, an intrabinary shock can form as a site of particle ac-celeration and associated non-thermal emission. We model the X-ray and gamma-ray synchrotron and inverse-Compton spectral components for select spider binaries, includ-ing diffusion, convection, and radiative energy losses in an axially-symmetric, steady-state approach. Our new multi-zone code simultaneously yields energy-dependent light curves and orbital-phase-resolved spectra. Using parameter studies and matching the observed X-ray spectra and light curves, and Fermi Large Area Telescope spectra where available, with a synchrotron component, we can constrain certain model parameters. For PSR J1723–2837 these are notably the magnetic field and bulk flow speed of plasma moving along the shock tangent, the shock acceleration efficiency, and the multiplicity and spectrum of pairs accelerated by the pulsar. This affords a more robust predic-tion of the expected high-energy and very-high-energy gamma-ray flux. We find that nearby pulsars with hot or flaring companions may be promising targets for the fu-ture Cherenkov Telescope Array. Moreover, many spiders are likely to be of significant interest to future MeV-band missions such as AMEGO and e-ASTROGAM.
Keywords: binaries: close – pulsars: individual (PSR B1957+20, PSR J1723−2837, PSR J1311−3430, PSR J2339−0533) – radiation mechanisms: non-thermal – X-rays: binaries – gamma rays: stars – stars: winds, outflows.
Abstract ii List of Figures vi List of Tables ix Abbreviations x 1 Introduction 1 1.1 Problem Statement . . . 1 1.2 Research Aims . . . 2 1.3 Research Outcomes . . . 3 2 Theoretical Background 4 2.1 Supernovae . . . 4 2.1.1 Type Ia supernovae . . . 5 2.1.2 Core-collapse supernovae . . . 5 2.2 Pulsars . . . 6 2.2.1 Observations . . . 6 2.2.2 Theoretical discussion . . . 9 2.2.3 Millisecond Pulsars . . . 10 2.2.4 Pulsar magnetosphere . . . 12
2.2.5 Pair multiplicity and pair cascades . . . 14
2.2.6 Striped wind structure . . . 17
2.3 Shocks . . . 18
2.3.1 Non-relativistic shocks . . . 18
2.3.2 Relativistic shocks . . . 21
2.3.3 Shock acceleration particle transport . . . 22
2.4 Radiation Mechanisms . . . 25
2.4.1 Synchrotron Radiation . . . 25
2.4.2 Inverse Compton Scattering . . . 28
2.4.3 Synchrotron Self-Compton In The Binary Pulsar Case . . . 30
2.4.4 γγ-absorption . . . 31
2.5 Doppler Boosting . . . 33
2.6 Binary Systems . . . 36
2.6.1 Classification by observational properties . . . 36
2.6.2 Classification by physical properties . . . 37 iii
2.7 Spider Binaries . . . 38
2.7.1 Observational properties . . . 38
2.7.2 Theoretical discussion . . . 42
3 UMBRELA - The Numerical Model 49 3.1 Formalism and Assumptions . . . 49
3.2 Particle Injection and Acceleration . . . 51
3.3 Shock Orientation and Division into Spatial Zones . . . 54
3.4 Particle Transport Equation . . . 55
3.5 Radiation Losses . . . 58
3.6 Beaming and Emission . . . 60
3.7 Shadowing . . . 67
3.8 Convergence and Energy Conservation . . . 69
4 Results 72 4.1 Case Study: Parameter Variation for PSR J1723−2837 . . . 72
4.1.1 Magnetic field . . . 73
4.1.2 Acceleration efficiency . . . 73
4.1.3 Particle injection spectral index . . . 73
4.1.4 Pair multiplicity . . . 74
4.1.5 IBS radius . . . 75
4.1.6 Companion temperature . . . 77
4.1.7 Mass ratio . . . 77
4.1.8 Bulk-flow momentum . . . 78
4.1.9 Minimum energy fits . . . 79
4.2 Application of the Model to Individual Sources . . . 80
4.2.1 PSR J1723−2837 (RB) – ICDP . . . 82 4.2.2 PSR J1311−3430 (BW) – ICDP . . . 83 4.2.3 PSR J2339−0533 (RB) – ICDP . . . 85 4.2.4 PSR B1957+20 (BW) – SCDP . . . 86 5 Radio Observations 89 5.1 Motivation . . . 89 5.2 MeerKAT . . . 90
5.3 FRB Open Time Proposal . . . 91
5.4 Method . . . 91
5.5 Preliminary Results, Interpretation, and Future Aspirations . . . 91
6 Conclusion and Future Work 93 6.1 Previous Work: Existing Voids . . . 94
6.2 UMBRELA . . . 95
6.3 Future Work . . . 96
A Timescales 97
2.1 Light curves of a number of Type Ia SNe after the luminosity-width rela-tion has been taken into account. Taken from Longair(2011). . . 5 2.2 This Table depicts the typical timescales for each stage of evolution for a
15M star. This Table was taken fromWoosley & Janka (2005). . . 6
2.3 A schematic model of a pulsar showing the misaligned rotation and
mag-netic axes. Taken from www.cv.nrao.edu . . . 8 2.4 P-Pdot diagram taken from www.cv.nrao.edu. . . 11 2.5 A schematic representation of the magnetic field and charge distribution
about a rotating magnetised NS according to the analysis of Goldreich & Julian(1969). Taken fronLongair(2011). . . 12
2.6 Schematic representation of the e± cascade above the PC of a young
pulsar. Obtained for illustrative purposes fromTimokhin & Harding(2019). 15 2.7 Schematic diagram representing an electron gyrating in a magnetic field
~
B. The parallel and perpendicular components of the electron’s velocity are also shown. . . 25 2.8 The spectrum of the SR of a single electron shown linear (left) and
loga-rithmic (right) scales . . . 27 2.9 γγ-absorption cross section as a function of x = 12(1 − µ), taken from
B¨ottcher et al. (2012) . . . 32 2.10 Frequency-dependent radio eclipses observed from the BW PSR B1957+20.
Image obtained from Ryba & Taylor(1991). . . 40 2.11 Optical light curve for PSR B1957+20 showing the dip in magnitude when
the companion comes between the pulsar and observer (pulsar superior conjunction). Obtained from Reynolds et al. (2007) . . . 41 2.12 X-ray light curves for (A) BW PSR 1957+20 and (B) RB PSR J1723-2837
showing the double-peaked structure. The shadow regions indicate where radio eclipses were observed for the source. Obtained from Huang et al.
(2012) and Hui et al.(2014), respectively. . . 42 3.1 Schematic diagram indicating a cross section of the shock wrapping around
the pulsar, with parameters defined as indicated. The pulsar wind is em-anating from the pulsar, indicated by green lines, and the particles are captured by and flow along the intrabinary shock. In this work, we ap-proximate the shock as a 3D thin hemispherical shell, rather than one of finite thickness and multi-layer structure (as alluded to in the inset). The particles are accelerated at the shock locale to ultrarelativistic energies and acquire a slight anisotropy in their steady-state distribution function (mildly relativistic “bulk motion”) along the shock tangent, as indicated by an increase in the bulk momentum βΓ and the corresponding blue colour. 50
3.2 Schematic of particle injection and transport in zones along the shock colatitude sectors. A fraction of QPSR is injected into the ith zone of
the shock (as indicated by large arrows). The small arrows between cells signify the direction of bulk motion of particles flowing away from the nose of the shock. . . 55 3.3 Schematic diagram depicting the geometry for the calculation of n∗. . . . 59
3.4 Schematic indicating the shock geometry and beamed emission from the
shock. The shock is around the pulsar. In the top panel, the inclination i = 75◦, while in the bottom panel i = 40◦. Blue indicates flow along the shock surface directed towards the observer, and red indicates flow away from the observer. The normalised orbital phase is indicated in the left corner of each panel. Obtained from van der Merwe et al. in press. For the BW case, refer to Fig. 3 inWadiasingh et al.(2017). . . 62
3.5 The plots show the integrated SR and IC luminosity per zone for an RB.
See description in text. . . 64
3.6 The plots show the integrated SR and IC luminosity per zone for a BW.
See description in text. . . 65
3.7 These plots show the effect of shadowing when θmax is varied, on the
model-predicted light curves and SEDs. The black lines on the plots rep-resent the output when shadowing is incorporated, other colours reprep-resent the output without shadowing. Here Rsh= 0.16a. . . 68
3.8 These plots show the effect of shadowing when Rshis varied, on the
model-predicted light curves and SEDs. The black lines on the plots represent the output when shadowing is incorporated, other colours represent the output without shadowing. . . 69 3.9 The Figures above show the limit where τeff,i/τmin,i 1. (A) shows the
energy input into each zone cumulatively and (B) shows the cumulative effect of energy lost due to SR, IC, and escaping particles. It can be seen that more energy is lost as one progresses through the zones along the IBS, increasing from a factor 2 to a factor 5 as the zone number increases. 71
4.1 Model SED plots for PSR J1723−2837 depicting the effect of varying the
B-field strength at the shock, Bsh. . . 73
4.2 Model SED plots for PSR J1723−2837 depicting the effect of varying
acceleration efficiency acc that limits γmax. . . 74
4.3 Model SED plots for PSR J1723−2837 depicting the effect of varying
injection spectral index p. . . 74
4.4 Model SED plots for PSR J1723−2837 depicting the effect of varying the
pulsar pair multiplicity Mpair. . . 75
4.5 Model SED plots for PSR J1723−2837 depicting the effect of varying the
shock radius Rsh. . . 76
4.6 Plot for PSR J1723−2837 depicting the zone-to-companion distance ρ
(solid lines) as well as the value of usp (dashed lines) versus zone for the
RB case. . . 76
4.7 SED plots for PSR J1723−2837 depicting the effect of varying the
com-panion temperature Tcomp. . . 77
4.8 SED plots for PSR J1723−2837 depicting the effect of varying the stellar mass ratio q. . . 78
4.9 Model SED plots for PSR J1723−2837 depicting the effect of varying the bulk-flow momentum βΓ. . . 78 4.10 Model SED plots for PSR J2339−0533 indicating two different fitting
scenarios: (1) large γminso that we fit X-ray data using the intrinsic
single-particle SR spectral slope of 4/3 (black line) and (2) using a lower γmin
so that the X-ray spectrum is matched by varying the particle spectral index p (gray line). . . 79 4.11 Plot for PSR J1723−2837 indicating the SED where we show two cases,
one matching Fermi LAT spectral data (black line) and the other not (grey line). . . 82 4.12 Plot for PSR J1723−2837 indicating the energy-dependent light curves.
The light curves correspond to the black SED model. The top panel (black line) is a model light curve for an energy close to that of the observed X-ray one. The blue (SR) and red (IC) lines in the bottom panel illustrate energy-dependent pulse shapes, for energies as indicated in the legend. . . 82 4.13 Plot for PSR J1311−3430 depicting the SED for the quiescent and flaring
states. The dashed grey line represents an SED model for the flaring state in which all parameters are identical to the quiescent state (black line), except the companion temperature and magnetic field at the shock. The solid grey line represents an alternative model to the flaring-state spectral data. . . 83 4.14 Plot for PSR J1311−3430 depicting the energy-dependent light curves
for the quiescent (blue and purple) and flaring (orange and red) states. We do not exhibit light curve data due to their low quality. The light curves should only be taken as a qualitative indication of the predicted double-peak shape. . . 84 4.15 The model SED for PSR J2339−0533. . . 85 4.16 The model energy-dependent light curves for PSR J2339−0533. . . 86 4.17 Plot for PSR B1957+20 depicting the model SED for both i = 65◦(gray)
and i = 85◦ (black). . . 87 4.18 Plot for PSR B1957+20 depicting the model energy-dependent light curves
for i = 85◦. The top panel shows a fit to X-ray data and the bottom panel shows a fit to available modulated GeV data. . . 87 5.1 The three images above show the post-processing results of the
observa-tions done for the three FRB posiobserva-tions. The current best location of the FRB is indicated by a cyan coloured dot at the centre of each image. Credit: J Chibueze. . . 92 A.1 Timescale plots for PSR J1723−2837 depicting the first and last spatial
zones. . . 97 A.2 Timescale plots for PSR J2339−0533 depicting the first and last spatial
zones. . . 98 A.3 Timescale plots for the quiescent state of PSR J1311−3430 depicting the
first and last spatial zones. . . 99 A.4 Timescale plots for PSR B1957+20 depicting the first and last spatial
zones for inclination angle i = 64◦ (top). Similarly, the bottom panels are for the alternative fit with inclination angle of i = 84◦. . . 99
4.1 Model Parameters for Illustrative Cases. . . 81
UMBRELA Unravelling the Multi-wavelength Beamed Radiation for Energetic Leptons in Arachnids
BW Black Widow
RB Redback
IBS Intrabinary Shock
PC Polar Cap
SED Spectral Energy Distribution
LC Light Cylinder MHD Magnetohydrodynamic SN Supernova NS Neutron Star BH Black Hole WD White Dwarf
LOS Line of Sight
RPP Rotation Powered Pulsar
PWN Pulsar Wind Nebula
MSP MIllisecond Pulsar
RRAT Rotating Radio Transients
SGR Soft Gamma-ray Repeaters
AXP Anomallous X-ray Pulsars
CCO Compact Central Objects
SNR Supernova Remnant
INS Isolated Neutron Stars
ANS Accreting Neutron Stars
LMXB Low-mass X-ray Binary
HMXB High-mass X-ray binary
WRXB Wolf-Rayet X-ray binary
GC Globular Cluster
CR Curvature Radiation
SR Synchrotron Radiation
SSC Synchrotron Self-Compton
ICS Inverse Compton Scattering
FFE Force-Free Electrodynamics
Introduction
1.1
Problem Statement
Black widow and redback systems are compact binaries in which a millisecond pulsar (MSP) heats and may even ablate its low-mass companion by its intense wind of rela-tivistic particles and radiation. The pulsar wind drives mass loss from the companion, and an intrabinary shock forms as a site of particle acceleration. Radio, optical, and X-ray follow-up of unidentified Fermi Large Area Telescope (LAT) sources has expanded the number of these systems from four to nearly 301.
It is believed that the main difference between the two types of binaries is that in the black widow case the intrabinary shock is formed around the companion with a very low mass (∼ 0.01M ) (Fruchter et al., 1988), whereas with the redback case the
intra-binary shock is formed around the pulsar and the companion star is more massive (∼ 0.2M ) (Roberts, 2011). These binary systems have orbital periods of a few hours
and the formation of an intrabinary shock between the pulsar and its companion is believed to be the reason for the observation of radio eclipses and double-peaked X-ray light curves that result from Doppler-boosted emission from the shock surrounding the companion (Wadiasingh et al.,2015). Observations show that these systems typically have strong non-thermal power-law components with hard photon indices of Γ ∼ 1 − 1.5, which implies a hard electron spectrum and efficient particle acceleration (Roberts et al.,
2015). Variable heating of the companion of various sources has been inferred via optical light curves as well as flux increasing during ‘flaring states’, which may be attributed to magnetic events on the rapidly rotating companion or the pulsar wind interaction with the intrabinary shock (An et al.,2017).
1
See the public Fermi LAT pulsar list.
These systems are very mysterious; many aspects of underlying physics of the system and its environment are not well-known. The mechanisms for particle emission and ac-celeration in these extreme astrophysical systems are not well understood. The structure of the intra-binary shock is unknown and is assumed to be spherical in this work, but a recent study by Kandel et al. (2019) implemented new geometries. The pulsars wind content and composition is another point of interest for the future of this research. An interesting question is whether there are some of these binary systems that may have very large TeV fluxes that we would be able to detect with telescopes such as the future Cherenkov Telescope Array (CTA) and the High Energy Stereoscopic System (H.E.S.S.) To address these questions, we model the X-ray and gamma-ray spectral components from the intrabinary shocks of nearby ‘spider binaries’, including diffusion, convection, and radiative energy losses in an axially-symmetric, steady-state approach. The code simultaneously yields energy-dependent light curves and orbital phase-resolved spectra. Using parameter studies and fitting X-ray and gamma-ray spectra and light curves, we will constrain certain model parameters and estimate the very-high-energy gamma-ray flux (> 100 GeV) for a few promising sources. We will demonstrate that nearby binaries and those in a ‘flaring state’ are promising targets for CTA, and may also be detectable by H.E.S.S. for optimistic parameter choices. Constraining the inverse Compton emission via TeV observations may be useful to probe the particle acceleration in the shock vicinity as well as the pulsar wind content.
1.2
Research Aims
The overarching aim of this dissertation is to obtain a better understanding of the underlying physics and environment of these mysterious millisecond pulsar binaries. This will be done by using and refining a numerical emission model that calculates the X-ray and gamma-ray spectral components from nearby spider binaries. The predictions for the synchrotron radiation (SR) spectrum can be fitted to the X-ray data that are found in literature to constrain certain model parameters such as the pulsar mass, companion stars temperature, magnetic field at the intrabinary shock, etc. Once the models have been calibrated using the X-ray and GeV data, we can use them to predict the TeV gamma-ray flux expected from some sources. The model also simultaneously produces energy-dependent light curves, which together with the X-ray spectra and confronting these predictions with data, will enable us to give a more robust prediction of the expected TeV gamma-ray flux. Constraining the gamma-ray flux, which is believed to be produced mostly by inverse Compton emission (ICS), will help probe the particle acceleration at the shock as well as the pulsar wind content. The detectability of these
sources by telescopes such as the future Cherenkov Telescope Array (CTA), Fermi Large Area Telescope (LAT) and H.E.S.S. will also be investigated during this dissertation. In our capacity as H.E.S.S. members, we will contribute to the observational planning of this instrument. In particular, we will attempt to monitor some sources optically, using optical indications of flaring as triggers for H.E.S.S. to observe these sources for some hours. A detection of a modulated TeV signal will strengthen our preliminary interpretation of shock acceleration and leptonic emission, and provide input for the planning of future CTA observations.
1.3
Research Outcomes
During the course of the research for this dissertation the following output was produced:
• Conference presentations:
– The 64th Annual Conference of The SA Institute of Physics (SAIP2019). – High Energy Astrophysics in Southern Africa meeting (HEASA2019). • Articles:
– van der Merwe et al. has been submitted to the Astrophysical Journal (ApJ) and will also be published on ArXiv.
• Observational proposals: Based on predictions made with the model described in this dissertation, we submitted a proposal to H.E.S.S. to observe two sources, namely PSR J1723-2837 and PSR J1311-3430, which was approved and observ-ing time was allocated to study PSR J1723-2837. Observations and analysis are ongoing.
Theoretical Background
In this chapter, I will give a brief discussion of various topics that will provide context for the work presented in the next chapters of this dissertation. I will start with dis-cussing supernovae and the formation and characteristics of pulsars in Sections 2.1 and 2.2. The discussion will then move to the pulsar magnetosphere in Section 2.2.4, the production of pair cascades therein in Section 2.2.5 and the properties of the pulsar wind in Section 2.2.6. Shocks and shock acceleration is discussed in Section 2.3. I will then look at the radiation mechanisms present in these systems in Section 2.4, followed by a discussion on Doppler boosting in Section 2.5. Lastly, I will give an overview of spider binaries - what they are, what we know and what we expect to see, in Section 2.7.
2.1
Supernovae
At the end of a very massive star’s life, the fuel burnt (mainly hydrogen and helium) to counteract the gravitational forces created by the star’s mass to maintain hydrodynamic equilibrium runs out, the star collapses, and an event called a supernova (SN) occurs. SNe are extremely violent and luminous stellar explosions, which occur when a star as a whole explodes and ejects its envelope at high velocity. Since the formation of neutron stars (NSs) and black holes (BHs) is believed to be associated with events that release extreme amounts of energy on very short timescales, SNe seem to be a natural candidate (as a progenitor) for explaining their existence (Woosley & Janka,2005). In what follows, I discuss the two main categories of SNe namely, Type Ia vs. all other SNe that are believed to result from the core collapse of massive stars.
2.1.1 Type Ia supernovae
Type Ia SNe are very unique in that they all seem to display very similar characteristics as can be seen in Fig 2.1. They are the most luminous SNe with typical absolute B magnitudes of MB = −19.5 ± 0.1 (Longair, 2011). This apparent consistency in their
characteristics enables the use of their light curves and the luminosity-width relation to determine the absolute magnitudes of very distant SNe, which has also contributed greatly to determining the redshift-distance relation out to redshifts of one or greater (Pinto & Eastman,2000). This class of SN is very rare and is associated with the nuclear explosion of carbon-oxygen white dwarfs (WDs) with masses close to the Chandrasekhar
mass of 1.4M , which is the mass above which WD become gravitationally unstable.
This limit can be reached if the WD has a binary companion from which it accretes matter onto its surface. The most favoured companion type is another WD star - the system would inspiral due to energy loss via gravitational radiation and the merger of the two WDs could then give rise to a Type Ia SN.
Figure 2.1: Light curves of a number of Type Ia SNe after the luminosity-width relation has been taken into account. Taken fromLongair(2011).
2.1.2 Core-collapse supernovae
During a star’s lifetime it will go through many different stages. As outlined inWoosley & Janka(2005), a star heavier than 8M will during its lifetime pass through successive
stages of hydrogen, helium, carbon, neon, oxygen, and silicon fusion in its centre. Once these burning processes have completed, the star will have formed a core of about 1.5M
containing iron-group elements. The typical timescale and temperature for each stage of a typical massive star’s evolution is given in Fig. 2.2 below.
Figure 2.2: This Table depicts the typical timescales for each stage of evolution for a 15M star. This Table was taken fromWoosley & Janka(2005).
The initiation of the core collapse is a result of two processes. First, at high densities (> 1010g cm−3) electrons are captured by the iron-group nuclei, increasing their neutron number and decreasing the number of unbound electrons, which reduces the pressure in the stellar interior and, thus its ability to maintain the star’s stable structure. Second, at high temperatures found in this stage of the star’s evolution, radiation begins to turn the iron nuclei into helium in a process called photodisintegration. This process occurs when a high-energy photon is absorbed by a nucleus, which is then split into lighter elements whilst releasing a neutron, proton, or alpha particle. The core then collapses to a hot, dense, and neutron-rich sphere called a proto-neutron star, where the collapse is suddenly stopped due to neutron degeneracy pressure (Podsiadlowski, 2014). This abrupt stop generates a shock wave that interacts with the still collapsing outer region of the core. Due to photodisintegration and neutrino losses, this shock wave does not eject the rest of the star. After a few milliseconds, all outward velocities are reduced to zero and the proto-neutron star starts to accrete the surrounding materials. If the accretion continues the proto-neutron star will collapse into a BH and there will not be an SN explosion. However, if the proto-neutron star does not collapse (despite a large luminosity of neutrinos), it will radiate ∼ 1053erg in a few seconds and eventually form a neutron star with a radius of approximately 10 km.
2.2
Pulsars
2.2.1 Observations
In 1967 a strange source consisting entirely of bright radio signals was discovered by
signal consisted of a series of pulses with a pulse period of around 1.33 s. The source, PSR 1919+21, was the first pulsar to be identified. Following the work ofPacini(1967),
Pacini(1968), andGold(1968), pulsars were identified as isolated, rotating, magnetised NSs. The observation of these pulses indicated that the NS’s magnetic and rotation axis need to be misaligned and that these pulses originate from beams of radio emission emitted along the magnetic axis of the NS (see Fig. 2.3). Thus, the radio emission from a pulsar is continuous, but beamed so that any observer will see a pulse if the beam sweeps across his/her line-of-sight (LOS).
Since their discovery, more than 2 800 pulsars have been catalogued1. Pulsars display a rich variety of characteristics, which has made it a daunting task to create a clearcut classification scheme. In this dissertation my focus is on MSPs and I will discuss them in more detail in Section 2.2.3. Below I will give an overview of the different types of pulsars.
Rotation-Powered Pulsars (RPPs) are NSs that are spinning down due to the torque created from the misalignment of their magnetic and rotation axes that leads to magnetic dipole radiation. They emit energy at multiple wavelengths and are detectable from
radio to gamma-ray wavebands (Harding, 2013). In pulsar wind nebulae (PWNe) the
presense of RPPs is also seen through a wind of energetic particles from a relatively young pulsar interacting with the surrounding material left behind after an SN explosion (Slane,2017). RPPs have two main subpopulations, namely normal pulsars and MSPs. The gamma-ray pulsar population has also been found to consist of radio-faint pulsars and gamma-ray MSPs (Abdo et al., 2009b,c). Another subclass of RPPs is Rotating Radio Transients (RRATs), which produce very short, bright radio pulses at irregular time intervals. For a detailed review on RRATs, see Keane et al. (2011).
Magnetars are NSs that derive most of their power from the energy stored in their extreme magnetic fields that have values of 1014− 1016G (Kaspi & Beloborodov,2017).
Historically there are two subclasses of magnetars - Soft Gamma-Ray Repeaters (SGRs) and Anomalous X-Ray Pulsars (AXPs). However, this distinction has become blurred in recent times. They respectively exhibit repeating soft gamma-ray and hard X-ray outbursts. They have longer spin periods than normal pulsars, but have very large spin-down rates (Longair, 2011). Both subclasses possess very high X-ray luminosities that cannot be explained by magnetic dipole spin down that powers RPPs (Harding,2005). The magnetic field of these sources are so strong that radiative processes need to be approached using quantum electrodynamics.
Compact Central Objects (CCOs) are X-ray sources that are detected close to the cen-tres of young supernova remnants (SNRs) that have no significant emission at other
1
wavelengths and also show no signs of having a binary companion (Harding, 2013). It is estimated that they have relatively weak magnetic fields in the range 1010− 1011 G.
They display thermal-like spectra that are usally well described using the sum of two high-temperature black bodies and small emitting radii (De Luca,2017). The formation and subsequent characteristics of this class of pulsar are still not well-understood and is a topic of ongoing research.
Isolated Neutron Stars (INSs) are thermally cooling sources that primarily emit soft X-rays, but also show faint optical and UV counterparts. These sources differ from CCOs in that they are not associated with SNRs or nebulae (Kaspi et al.,2004). They also display higher magnetic fields of ∼ 1013 G and appear to be nearby middle-aged NSs (Harding, 2013). Very few of these objects have been discovered and the limited emission of these objects makes it a challenging task to better understand these sources. Accreting Neutron Stars (ANSs) are NSs that are found in a binary system where matter is accreted via the companion star from either an accretion disc in the case of Roche-Lobe overflow, or from the stellar wind (Longair, 2011). The gravitational energy of the infalling matter is responsible for a fraction of the observed radiation seen from these systems. ANSs are found in three classes of binaries, namely low-mass X-ray binaries (LMXBs), high-mass X-ray binaries (HMXBs), and gamma-ray binaries. These subclasses are discussed in Section 2.6.
Figure 2.3: A schematic model of a pulsar showing the misaligned rotation and magnetic axes. Taken from www.cv.nrao.edu
2.2.2 Theoretical discussion
I will now summarise various characteristics that we can derive for normal RPPs. The spin period of a pulsar (P ) can be measured with high accuracy and the change in this parameter over time ( ˙P ) is very important in determining various characteristics of a pulsar. Pulsars spin down over time due to loss of rotational energy, which can be described with a braking index n defined by ˙Ω = −κΩn (Manchester & Taylor,1977), where Ω = 2π/P is the angular frequency of rotation.
As mentioned above, the magnetic and rotation axes of a pulsar are typically misaligned. This results in a varying dipole moment that causes the pulsar to lose energy due to magnetic dipole radiation. In the case of a magnetic dipole rotating at a magnetic inclination angle α with respect to the rotational axis, with a magnetic dipole moment m = 12BR3, the amount of power radiated is given by the Larmor formula
Prad= 2 3 ¨ m2⊥ Ω4 c3 = 2 3c3(BR 3sin α)2 2π P 4 , (2.1)
where Ω is the angular velocity of the NS and m⊥ is the magnetic dipole component
that is perpendicular to the rotation axis.
The rate of loss of rotational energy from the NS can be determined directly from the spin-down rate of the pulsar (Roberts et al.,2005)
− dErot dt = −d dt 1 2IΩ 2 = −IΩ ˙Ω = −I 2π P −2π P2P˙ = 4π 2I ˙P P3 . (2.2)
If the second derivative of the angular frequency is known, the braking index can be determined. For a constant braking index n, the second derivative can be written as
¨
Ω = −nκΩ(n−1)Ω = −n ˙˙ Ω2Ω−1 and then solving for n yields
n = Ω ¨Ω ˙ Ω2 = 2 − P ¨P ˙ P2 . (2.3)
For a dipolar magnetic field, n = 3. Assuming that the deceleration of an NS can be described by a braking index throughout its lifetime, integrating ˙Ω = −κΩn yields
1 (n − 1)
h
Ω−(n−1)− Ω−(n−1)0 i= κτ, (2.4) where τ is the age of the pulsar and Ω0is its initial angular velocity. If we further require
that n > 1 and Ω0 Ω, the age of a pulsar can be estimated,
τ = Ω −(n−1) κ(n − 1) = − Ω (n − 1) ˙Ω = P (n − 1) ˙P . (2.5)
This is called the characteristic age of a pulsar. Using this expression, it can be inferred that the typical lifetime for a normal pulsar is about 105− 108 years (Longair,2011).
Lastly, we can also estimate the magnetic field flux density at the surface of the NS. Assuming all the power that is radiated is produced by the loss of rotational energy, i.e. ˙E ≈ Prad, a lower limit to the magnetic field strength can be determined (Lyne &
Graham-Smith,2012). From Eq. (2.1) and (2.2), 2 3c3( 1 2BR 3sin α)2 2π P 4 = 4π 2I ˙P P3 , (2.6) ∴ B2 = 2c 3IP ˙P 3π2R6sin2α , (2.7) ∴ B > 2c3I 3π2R6 1/2 (P ˙P )1/2 , (2.8)
substituting in typical values (CGS units),
∴ B = 2 (3 × 10 10 cm s−1)3· (1045 g cm2) 3π2 (106 cm)6 1/2 (P ˙P )1/2 , (2.9) ∴ B & 4.3 × 1019 (P ˙P )1/2 G . (2.10)
The derivation above assumes a dipolar magnetic field structure, however studies of mul-tipolar fields and their effect on the properties mentioned above have been performed by
Petri(2015) andP´etri(2019). It is also good to note, that if other angular-momentum-loss mechanisms (e.g. retrograde accretion) play a role, the magnetic field could be smaller than the above estimate.
Since most of these characteristics can be calculated using P and ˙P , a large amount of information can be represented by a plot called the P − ˙P diagram. An example is shown in Fig. 2.4.
2.2.3 Millisecond Pulsars
In 1982, Backer and his colleagues discovered the first millisecond pulsar (MSP), B1931+20, which has a pulse period of 1.56 ms (Backer et al., 1982). Today, well over 300 MSPs have been discovered2. MSPs have relatively weak magnetic fields compared to pulsars and stable spin periods. Their weak magnetic fields (∼ 108− 1010G) and small ˙P values
imply that they are very old pulsars, typically 109− 1010years (Becker & Pavlov,2002)
of age. This unique group is located in the bottom left of Fig. 2.4, where those which are 2
Figure 2.4: P-Pdot diagram taken from www.cv.nrao.edu.
members in binary systems are enclosed by circles. As seen in the Figure, the majority of the MSP population are found to have a companion.
One of the popular formation scenarios for MSPs is the so-called ‘recycling scenario’ (Alpar et al., 1982). This theory argues that MSPs were originally ordinary pulsars which, after long periods of time, have lost so much rotational energy that they move down in P − ˙P space into the graveyard region (bottom right Fig. 2.4) and ‘turn off’. A pulsars ‘turns off’ when the rotational period becomes so long that the polar cap (PC) voltage falls below a critical value needed to create electron-positron pairs, thought to be vital for creating radio emission. The short periods of MSPs are naturally explained by the presence of a companion. The ‘dead’ pulsar can then again be brought to life for an encore when it starts accreting matter from a binary companion. The companion could be acquired through gravitational capture by the heavy NS (if it was originally not a binary) or has always been part of the binary system - the more likely scenario is an open question that can only be answered through extensive population synthesis studies. When mass is transferred from the companion to the NS, it is spun up. During the accretion phase, significant X-ray emission is produced and the NS is observed as a low-mass X-ray binary (Kaspi & Kramer,2016). If the pulsar has weak field, a large spin up can occur, since the spin-up rate is limited by the surface magnetic field of the given pulsar (Verbunt et al.,1987). This spin up through accretion is then believed to be the reason for the millisecond periods observed from these sources. Once the accretion phase has finished, the spun-up pulsar may become rotationally powered (again exceeding the critical polar cap potential needed to spark radio emission) and be observed as an MSP.
Some MSPs have been found to transition between rotation-powered radio MSPs and low-mass X-ray binaries (Archibald et al.,2009;Papitto et al.,2013), giving more insight, and also verifying the recycling theory in the evolutionary track of MSPs.
The presence of LMXBs in globular clusters (GC) suggests that MSPs, because of their association with close binary systems, should also be present in GCs. The discovery of transitional MSPs lends weight to the recycling scenario and the high density of stars in GCs also makes the scenario of the low-mass companion being gravitationally captured by the ‘dead’ pulsar more plausible. This seems to be a reasonable assumption, as 25
pulsars3 have been detected in the nearby GC, 47 Tucanae and more in Terzan 5, all
having pulse periods faster than 6 ms (Longair,2011).
The formation of isolated MSPs is not well-understood, but these are believed to have originally been in a binary system in which the companion has been completely evapo-rated by the intense pulsar wind and radiation or the system was tidally disrupted after the formation of the MSP (Becker & Pavlov,2002).
2.2.4 Pulsar magnetosphere
Figure 2.5: A schematic representation of the magnetic field and charge distribution about a rotating magnetised NS according to the analysis ofGoldreich & Julian(1969).
Taken fronLongair(2011).
The magnetosphere is defined as the region in which the motion of charge particles are dominated by the extreme magnetic fields produced by the pulsar. In the simplest approximation, a pulsar is considered to be a perfectly conducting sphere initially sur-rounded by vacuum, with its magnetic dipole moment ~p0 aligned with the rotation axis
of the pulsar. A uniform magnetic field is frozen into a sphere that rotates at angular frequency Ω. An induced electric field is expected, but as a result of the infinite conduc-tivity of the star, is cancelled out by electric charges rearranging themselves such that (Goldreich & Julian,1969)
~
E + (~v × ~B) = ~E + [(~Ω × ~r
c ) × ~B] = 0. (2.11)
This results in a charge distribution within the pulsar that can be calculated using ~
∇ · ~E = ρ/0. This charge distribution needs to match the external vacuum solution of
the Laplace’s equation, ∇2Φ = 0. The external electrostatic potential is of quadrupolar form (Larmor,1884) φ = −B0ΩR 5 NS 6πr3 3 cos 2θ − 1 , (2.12)
where B0 is the polar magnetic flux density, RNS is the radius of the NS, and θ is the
polar angle as measured from the magnetic axis. The above expression then implies that there is a surface charge distribution on the sphere.
It was shown by Goldreich & Julian (1969) that this region surrounding the star will not remain a vacuum, but will actually be filled with plasma. This is the result of the force exerted by enormous induced electric field, exceeding the work function of charges on the surface and dragging the charges off the surface into the magnetosphere. Thus, there is fully-conductive plasma surrounding an NS and electric currents can flow in the magnetosphere. The charge distribution within this corotating magnetosphere is given by (Holloway & Pryce,1981)
ρe= 0∇ · ~~ E = 0B0 R3 NS r3 1 − 3 cos2θ . (2.13)
Setting ρe = 0 and solving for θ, it can be shown that it displays the property of
separating positive and negative charges along so called ‘zero charge cones’ or ‘null charge surfaces’ at an angle θ = cos−1(1/3)1/2 to the magnetic axis of the NS for an aligned rotator (for a static dipolar magnetic field).
The magnetic field and charge distribution for an aligned, rotating, magnetised NS is illustrated in Fig. 2.5. A very important parameter to define is the light cylinder (LC) or corotating sphere which is located at rLC = Ωc (Longair, 2011). At rLC the speed
of material corotating with the NS is equal to the speed of light. Within this LC, closed field lines are in a dipolar configuration and particles attached to these field lines corotate with the NS and form part of the corotating magnetosphere. Field lines that extend through the LC are called open field lines and particles attached to these lines
can escape to infinity. Outside the LC radius, the magnetic field lines follow a spiral configuration when viewed form above.
The PC boundary is defined by the field lines that are tangental to the LC. Charged particles within this region are thus attached to open field lines and escape to infinity. To calculate the angle subtended by the PC, one can use the fact that a particular line satisfies the relation k = sin θ/r, with k a constant labelling that line. Thus for small values of θPC, the angular radius of the PC is (Longair,2011)
θPC= RNS rLC 1/2 = ΩRNS c 1/2 (2.14)
and the radius of the PC is
RPC= θPCRNS≈ R3/2r −1/2
LC . (2.15)
The potential difference ∆φ between the pole and radius of the LC can be found by using Eq. (2.11), setting r = RNS, requiring θPC to be small and by substituting into
Eq. (2.14), ∆φ ≈ ΩB0R 2 NSθ2PC 2 ≈ Ω2B0R3NS 2c ∝p ˙E. (2.16)
2.2.5 Pair multiplicity and pair cascades
In the previous subsection, I gave a brief overview of the pulsar magnetosphere and the important quantities used to explore the physics just outside an NS’s surface. In this subsection, I will discuss the pair multiplicity and pair cascades expected to be present around all pulsars, which play a very important role in understanding not only the emission we see from pulsars, but also aids in unravelling the mystery of the complex magnetic field structure of pulsars.
Shortly after the discovery of pulsars (Hewish et al.,1968), different models describing the surrounding area outside a rotating NS were published byGoldreich & Julian(1969) andSturrock(1971). These models proposed that the magnetosphere is filled with dense pair plasma which screens the electric field as opposed to another model byHoyle et al.
(1964) that concluded that the plasma density surrounding the pulsar is very low. This pair plasma is believed to be created by gamma-rays in the very strong magnetic field found near the PCs, since the strong magnetic and electric field near the surface of an NS create an environment that is favourable for magnetic pair cascades to be initiated. The rapidly rotating NS induces a large electric field above the PC that accelerates particles along the curved, open magnetic field lines to Lorentz factors exceeding 107
(Daugherty & Harding, 1982). This acceleration along a curved trajectory leads to radiation called curvature radiation (CR). We did not model the emission mechanisms of the pulsar and its surrounding environment explicitly in this dissertation, thus I will not describe the detail of CR - for an excellent overview of the basic theory, seeTimokhin & Harding(2019). If the resulting CR from the accelerated particles extends to energies of ∼ 1 GeV, the curvature photons have a high probability of producing an electron-positron pair (e±) in a process called magnetic photon pair production (γ + B → e±) (Arons & Scharlemann,1979).
Figure 2.6: Schematic representation of the e± cascade above the PC of a young pulsar. Obtained for illustrative purposes fromTimokhin & Harding(2019).
As primary electrons are accelerated along the magnetic field lines, they emit CR photons almost tangential to the magnetic field lines. These CR photons are absorbed by the magnetic field when the angle between their momentum vector and the magnetic field becomes large enough and each photon spawns an e±pair. These pairs have a momentum component perpendicular to the magnetic field, which is almost instantaneously radiated via SR, thus continuing on their trajectory along the magnetic field lines (Daugherty & Harding, 1982). The secondary particles remain relativistic, but have much lower energies than the primary electron and will not be able to create subsequent pairs via CR. The synchrotron photons produced by the initial e± pair are also emitted almost tangentially to the magnetic field and once the angle between their momentum and the magnetic field becomes large enough, they will generate the next generation of pairs (see Fig. 2.6). This production of pairs, initiating a cascade effect, will continue until the
produced SR photons have energies falling below the critical energy needed to create an e± pair. When the energetics of the accelerated particles drops below the minimum energy required to produce e± pairs, the photon has then crossed the so-called ‘death line’ (Sturrock,1971). A study of the conditions of pair creation and subsequent death lines was performed by Harding et al.(2002).
The total number of pairs that can be created in a cascade is determined by a combi-nation of factors. First, the number of primary electrons and their initial energy are important, as larger values for both quantities will lead to more pairs being produced. Second, the energy threshold for the formation of pairs is vital - a lower threshold will produce more pairs. Lastly, the efficiency of particle energy going into pair creation in-stead of kinetic energy of particles and photons below the formation threshold also plays a role. Typical values for the pair multiplicity (Mpair) from pulsar models is 102− 105
(Timokhin & Harding,2019). Monte Carlo simulations have revealed the spectra of the pair cascades to be hard power laws (∼ 1.5) with low- and high-energy cutoffs that are dependent on the pulsar period and its derivative Harding & Muslimov(2011).
For MSPs, the low-energy cutoff is around γpair,min ∼ 104 and the high-energy cutoff
varies between γpair,max ∼ 106−107. The pair production by MSPs is not well-understood
since the surface magnetic field of most MSPs is too low to initiate pair cascades in dipolar fields. This led to suggestions by Ruderman & Sutherland(1975) and Arons & Scharlemann (1979) that non-dipolar magnetic fields near the NS surface might be the key to producing the pair multiplicities needed to explain the radio emission seen from these sources. Studies were then initiated byHarding & Muslimov(2011) into magnetic fields that are not purely dipolar. They found that deviations from a purely dipolar field could enable pair cascades and significantly increase the pair multiplicities for MSPs. This effect is believed to be caused by the extension of the spectra to lower pair energies that results from the offset-PC configuration. Recent results on MSP J0030+0451 (Miller et al., 2019; Riley et al., 2019) from the Neutron Star Interior Composition Explorer (NICER) implied a hotspot geometry on the NS surface that was far from antipodal, strongly suggesting the presence of non-dipolar fields as well as a severely off-set dipole. This supports the suggestion that MSP surface fields are non-dipolar, giving them the ability to produce a higher than expected pair multiplicity.
Thus the production of e± pairs by pulsars plays an essential role in producing charged particles in their magnetospheres and winds, as well as the generation of the pulsed radio emission we observe. In Section 3 I will discuss how pair multiplicity is incorporated into our model.
2.2.6 Striped wind structure
The magnetisation of pulsar winds is defined as the ratio of the magnetic energy flux to the kinetic energy flux in the flow (Lyubarsky,2003)
σ = B
2
8πN mc2Γ , (2.17)
where B is the strength of the magnetic field, N is the plasma density in the emission frame and Γ is the bulk Lorentz factor of the flow.
The magnetisation of pulsar winds has been inferred to be high (σ 1) interior to the LC, where the wind is magnetically dominated by the pulsar’s strong magnetic field, and weak (σ 1) farther away from the pulsar, usually at the shock formed in binary systems, which indicates that the winds are particle dominated here (Kennel & Coroniti,
1984b). This change in σ is not well-understood and is known as the ‘σ-problem’. It has been suggested that the accelerations of the particles in the pulsar wind may be driven by magnetic reconnection, which is a very efficient acceleration mechanism and can explain X-ray and gamma-ray emission observed from sources if it occurs at a high enough rate (Kirk & Skjæraasen,2003). To incorporate magnetic reconnection, the theory of ‘striped winds’ has come to the fore. In what follows, I shall give a brief overview of the pulsar wind and formation of striped winds.
The energetics of the plasma inside the magnetosphere are mostly dominated by elec-tromagnetic fields, thus particle inertia can be neglected in a first approximation - this is called force-free electrodynamics (FFE;Mochol (2017)). It was explained in the pre-vious sections that it is now believed that the pulsar magnetosphere is mostly filled with plasma and corotates with the NS. A consequence of this is that a Poynting vector exists , where ~S = 4πc E × ~~ B, which is responsible for transportation of energy from the NS into a wind (Michel, 1973a). For an oblique rotator with a magnetosphere filled with plasma, the spin-down luminosity is given bySpitkovsky (2006)
L ≈ 3
2L0(1 + sin χ) , (2.18)
where L0 = Bs2Ω4R6NS/6c3, Bs is the magentic field at the NS equator, RNS is the NS
radius and χ is the angle between the magnetic and rotation axes.
The topology of the magnetic field lines in the FFE description exhibits the formation of a current sheet outside the LC (P´etri, 2012). In the commonly accepted pulsar model, the magnetic field is believed to have a dipolar geometry inside the LC and changes smoothly to a monopolar, radial geometry at infinity. In the FFE approach (Michel,1973b), a monopole configuration is radial in the meridional plane and creates
an Archimedean spiral in the equatorial plane regardless of the magnetic field topology inside the LC. However, one can also implement a split monopole configuration (Michel,
1969), where two half magnetic monopoles of opposite polarity are joined together in the equatorial plane (Bogovalov, 1999). Then according to Ampere’s law, a current sheet must form, within which the magnetic field vanishes and the pressure needed to support equilibrium is supplied by hot plasma (Mochol,2017).
For an oblique rotator, the current sheet oscillates around the equatorial plane as the pulsar rotates, connecting the equator with field lines of opposite polarity every half period (Mochol, 2017). The current sheet, now forming into parallel regions, can be approximated by spherical shells, separating the regions of opposite magnetic polarity -creating a ‘striped wind’ structure (Coroniti,1990). In the split-monopole configuration, the maximum Poynting flux is located at the equator, thus most energy is carried out within a wedge of the striped wind defined by θ = π/2 ± χ.
2.3
Shocks
Normally, perturbations in a gas are propagated away from the source at the speed of sound in that specific medium. However, if a disturbance propagates at velocities greater than the speed of sound in the medium, it cannot behave like a sound wave would - a discontinuity is formed between the regions before and behind the disturbance, called a shock wave. In the following subsections, I will describe basic shock physics and then also discuss relativistic shocks and jump conditions associated with them.
2.3.1 Non-relativistic shocks
In order to make this a fruitful discussion, I will first define a few quantities and the scenario we will be considering.
Consider a discontinuity between two regions of the fluid flow. Ahead of the shock (upstream), the gas has a pressure p1, density ρ1, temperature T1 and the speed of
sound in this region is c1. The region behind the discontinuity (downstream) has a
pressure, density, and temperature of p2, ρ2, and T2 respectively. Adopting a reference
frame that moves at a velocity U , in which the shock wave is stationary will aid us in this discussion. In this frame, the undisturbed gas ahead of the shock moves towards the shock at a velocity v1 = |U | and when it passes through the shock, the gas moves
The behaviour of the gas as it moves through the shock wave is described by a well-known set of conservation relations, described inLongair(2011). First, mass flux needs to be conserved in this process, hence
ρ1v1 = ρ2v2. (2.19)
Second, the energy passing through a unit area parallel to v1, per unit time (i.e., the
energy flux) must be conserved. Considering only shock waves that are perpendicular to v1 and v2, this requirement can be written as
ρ1v1 1 2v 2 1 + w1 = ρ2v2 1 2v 2 2+ w2 , (2.20)
where w = m+ pV is the enthalpy per unit mass, m is the internal energy per unit
mass, and V is the specific volume V = 1/ρ.
Finally, the momentum flux must also be conserved as the gas moves through the shock,
p1+ ρ1v21 = p2+ ρ2v22. (2.21)
These three conditions are the well-known relations referred to as the shock conditions. There are many informative results that can be obtained using these shock conditions. For simplicity, consider a perfect gas for which the enthalpy is given by w = γ0pV /(γ0−1), where here γ0 = CP/CV is the ratio between the specific heat capacities. The following
useful relations, called the Rankine-Hugoniot relations, can then be established (Landau & Lifshitz,1959) V2 V1 = p1(γ + 1) + p2(γ − 1) p1(γ − 1) − p2(γ + 1) , (2.22) T2 T1 = [2γM 2 1 − (γ − 1)][2 + (γ − 1)M12] (γ + 1)2M2 1 , (2.23) p2 p1 = 2γM 2 1 − (γ − 1) (γ + 1) , (2.24) ρ2 ρ1 = v1 v2 = (γ + 1) (γ − 1) + 2/M12 , (2.25)
where M1 is the Mach number defined by M1 = U/c1 = v1/c1 = v1/γp1V1. In the strong
shock limit, i.e., M1 1, these expressions become
T2 T1 = 2γ(γ − 1)M 2 1 (γ + 1)2 , (2.26) p2 p1 = 2γM 2 1 (γ + 1), (2.27) ρ2 ρ1 = (γ + 1) (γ − 1). (2.28)
What these relations show is that while the temperature and pressure can become arbi-trarily large, the density ratio attains a specific maximum value. For a monatomic gas that has γ = 5/3, the density ratio is four in the limit for strong shocks. This shows how effective shock waves can be at heating gas to very high temperatures, which is observed in SN explosions and their remnants.
In high-energy astrophysics there are many situations that involve an object being driven supersonically into a gas, or equivalently, supersonic gas flowing past a stationary object. A simple example to illustrate these concepts was set out by Landau & Lifshitz(1959), wherein a piston is driven supersonically into a cylinder containing stationary gas. In the reference frame of the shock front, which moves with velocity vs, the velocity of the
piston is U , the velocity of the inflowing stationary gas is v1 = |vs| and the gas moving
behind the shock (i.e., between the cylinder head and the shock front) is v2.
Using the relations in Eq. (2.22)-(2.25), it can be shown that the speed of the shock is in this scenario is given by (Longair,2011)
vs= (γ + 1) 4 U + [c 2 1+ (γ + 1)2U2 16 ] 1/2. (2.29)
This expression determines the distance between a piston that is moving at a supersonic speed and the shock front. In the case a of a very strong shock wave, i.e. U c1, this
expression becomes
vs
U =
(γ + 1)
2 . (2.30)
For a monatomic gas, vs/U = 4/3. This is significant, as it shows the existence of a
stand-off distance at which the shock is located that is intrinsically determined by the pressure balance of the interacting media.
All sound waves propagate roughly at speeds of v ≈ (/ρ)1/2, where is the energy density in the medium (Longair,2011). Thus, in plasma that is magnetically dominated ( = B2/8π), the speed at which hydrodynamic waves can propagate is v ≈ B/4π(ρ)1/2, which is known as the Alfv´en speed.
2.3.2 Relativistic shocks
Shocks appear in many astrophysical systems, and since large electromagnetic fields and relativistic particles are usually present, an MHD approach is needed to evaluate the shock conditions. For the sake of this dissertation, I will focus only on the shocks that form in systems in which highly relativistic particles are present, specifically in this case a pulsar interacting with its environment. In this system, the relativistic wind produced by the pulsar interacts with the surrounding media and will form a shock. In what follows, I will highlight the assumptions and considerations made to derive the Rankine-Hugonoit relations for a similar system (a PWN) inKennel & Coroniti(1984a). The pulsar’s spin-down luminosity, just ahead of the shock, can be divided into particle and magnetic luminosity as follows:
L = 4πnγuR2shmec3(1 + σ) , (2.31)
where n is the proper density, u is the radial speed of the flow, γ2 = 1 + u2, Rsh is the
radial distance of the shock from the pulsar, me is the electron mass, c is speed of light,
σ is the ratio of magnetic plus electric energy to particle energy flux,
σ = B
2
4πnuγmec2
, (2.32)
and B is the magnetic field measured by the observer.
To simplify the derivation of the shock conditions a few assumptions were made. It was assumed that the magnetic field has a toroidal structure far from the pulsar. The azimuthal and latitudinal components of the flow velocity were neglected and the wind’s pressure was assumed to be negligible compared to the dynamic and magnetic pressure ahead of the relativistic shock. Lastly, the shock was considered to be stationary between the pulsar and the observer and fluids were assumed to enter the shock at a 90◦ angle. The MHD conservation laws are then given by:
n1u1 = n2u2, (2.33) E = u1B1 γ1 = u2B2 γ2 , (2.34) γ1µ1+ EB1 4πn1u1 = γ2µ2+ EB2 4πn1u1 , (2.35) µ1u1+ P1 n1u1 + B 2 1 8πn1u1 = µ2u2+ P2 n1u1 + B 2 2 8πn1u1 , (2.36)
where the subscripts indicate (1) upstream and (2) downstream parameters. E and B denote the shock-frame electric and magnetic fields, µ denotes the specific enthalpy, and the number density in the shock frame can be denoted by N = γn.
Using these relations, and also assuming that P1 and u2/u1 are small, and also that
u1/γ1 ≈ 1, Kennel & Coroniti (1984a) showed that downstream flow speed, magnetic
field, pressure, temperature, and density ratios depend solely on σ. When σ is large, u22 = σ + 1 8+ 1 64σ + · · · , (2.37) γ22 = σ + 9 8+ 1 64σ + · · · , (2.38) B2 B1 = N2 N1 = 1 + 1 2σ, (2.39) T2 u1mc2 = 1 8√σ 1 − 19 64σ + · · · (2.40) For σ < 0.1, u22 = 1 + 9σ 8 , (2.41) γ22 = 9 + 9σ 8 , (2.42) T2 u1mc2 = √1 18(1 − 2σ) , (2.43) B2 B1 = N2 N1 = 3(1 − 4σ) . (2.44)
In the limit σ → 0, the downstream speed u2/γ2 approaches 1/3, B2/B1 and N2/N1
ap-proach 3. The downstream particle pressure P2 and dynamical pressure Π2, normalised
to the upstream dynamical pressure, approach P2 n1mc2u21 = 2 3, Π2 n1mc2u21 = 1 3, (2.45)
which are the expected results for a relativistic, hydrodynamic shock. Thus, shocks with large σ-values are not very effective in converting the total upstream energy flux into thermal energy (which leads to SR discussed in Section 2.4), whereas shocks with small σ-values seem to be much more efficient in converting wind luminosity into SR luminosity.
2.3.3 Shock acceleration particle transport
An important part of modelling the physics and possible emission of shocks is the accel-eration of the particles as they move not only across the shock, but also along it. Here
I will give a short overview of this concept.
Fermi’s original picture (Fermi,1949) paints a scenario where a test particle is reflected from an infinitely massive ‘magnetic mirror’ that is a result of instabilities in the Galactic magnetic field. These mirrors are moving in random directions with a velocity v = βc, thus particles can collide with the mirror when its moving towards them (gaining energy) or when the mirror is moving away (losing energy). Fermi argued that the probability for head-on collisions is greater than for head-tail collisions, so the particles should gain energy on average. The particles have a wide distribution of pitch angles, and thus the average energy gain per collision can be shown to be (Longair,2011)
∆E E
∝ β2 . (2.46)
This process, called second-order Fermi acceleration, will lead to an exponential increase in the particle’s energy, since it will increase by the same fraction after each collision it encounters.
First-order Fermi acceleration, on the other hand, considers a shock geometry where particles can diffuse and cycle through unshocked and shocked regions by crossing the velocity discontinuity. The basic idea is that after a particle moves across the shock (from upstream to downstream) it can encounter a moving change in magnetic field that can reflect it back across the shock towards the upstream region at an increased velocity. The process of a particle initially moving from downstream to upstream regions will also result in a gain of energy. The average gain in energy when crossing the shock is (Bell,1978).
∆E E
∝ β . (2.47)
However, particles are also convected downstream and some escape during this process. Taking this into account,Bell(1978) showed that the differential spectrum of the high-energy particles is
N (E) dE ∝ E−2 dE . (2.48)
There is, however, a mystery to first-order Fermi acceleration. Ionisation losses restrict the acceleration of particles starting at low energies. For particles to be able to move across the shock, the particles need to be injected with energies greater than the asso-ciated maximum energy loss rate or the initial acceleration process must be very rapid as to overcome the ionisation energy losses. This dilemma is known as the injection problem (Balogh & Treumann,2013).
The typical differential energy spectrum describing particle energies is given by (Longair,
2011)
N (E) dE = constant × E −(1+(ατesc)−1) dE , (2.49)
where τesc is the time a particle spends in the acceleration region and α is used to
describe the fractional energy gain of a particle via ˙E = −αE.
Once the injection spectral shape is known, the diffusion-loss equation can be used to evaluate the time evolution of the particle spectrum. The equation is given by
∂N (E)
∂t = D∇
2N (E) + ∂
∂E[b(E)N (E)] + Q(E, t) , (2.50)
where D is the spatial diffusion coefficient and b(E) represents energy gains/losses per unit time. The first term describes the diffusion of the particles, the second term de-scribes energy gains/losses, and the final expression dede-scribes particle injection rate (and spectrum).
Efforts can also be made towards evaluating much more complex equations that describe the transport of particles in shock environments that include many different effects. In this dissertation, we solve a Boltzmann-type convection-diffusion equation, which includes the effects of radiative losses and a spatially-independent diffusion coefficient κ(E) given by (van Rensburg et al.,2018;Ginzburg & Syrovatskii,1964)
∂N (E) ∂t = −~V · ~∇Ne + κ(E)∇2Ne+ ∂ ∂E( ˙EtotNe) − ~∇ · ~V + Q(E, t) , (2.51)
where Ne is the differential particle energy distribution, ~V is the bulk velocity of the
plasma as it moves along the shock, and ˙Etot is the total energy loss rate due to radiation
and adiabatic losses.
Shock structures are found in many astrophysical environments, and are sometimes found in binary systems. Particle acceleration in these shocks plays an important role in the current understanding of how high-energy emission arises from many different types of binary systems. Studies into the shape of these shocks in binary systems have been conducted in some detail by Canto et al. (1996) and Wilkin (1996), but in this dissertation we assume, as a first approximation, a hemispherical shock. I discuss these concepts and how they are incorporated in the work of this dissertation in Chapter 3.
Figure 2.7: Schematic diagram representing an electron gyrating in a magnetic field ~
B. The parallel and perpendicular components of the electron’s velocity are also shown.
2.4
Radiation Mechanisms
2.4.1 Synchrotron Radiation
In the non-relativistic case, a single charged particle being accelerated by some force will radiate power according to the Larmor formula (Rybicki & Lightman,1979)
− dE dt = 2q 2a2 3c3 , (2.52)
where q is the charge of the particle, a is the acceleration of the particle, and c is the speed of light in a vacuum. For particles moving at relativistic velocities this expression becomes − dE dt = 2q 2γ4 3c3 (|a⊥| 2 + γ2 ak 2 ) . (2.53)
Let us now consider what happens if we subject the electron to a magnetic field ~B. The equation of motion of the electron under the influence of a magnetic field ~B can be described by the following formula involving the Lorentz force
meγ d~v dt = e c ~ v × ~B. (2.54)
This equation can be separated into a component along the field (vk) and in the plane
~v, dvk dt = 0, dv⊥ dt = e γmec ~ v⊥× ~B . (2.55)
These equations encapsulate that the electron with follow a helical path as it follows a combination of circular motion (v⊥) and a uniform motion along the field (vk),
respec-tively. Substituting these equations back into Eq. (2.53),
− dE dt = 2e 2γ4 3c3 evB sin α γmec 2 = 2e 4γ2v2B2sin2α 3m2 ec5 , (2.56)
where α is the pitch angle between the direction of the magnetic field lines and the velocity vector of the electron. Rearranging the above expression,
− dE dt = 2 " 8π 3 e2 mec2 2# cv c 2B2 8π γ 2sin α (2.57)
The first term in square brackets on the right-hand side is the well known Thomson cross section σT. Therefore,
− dE
dt
= 2 σTc β2Umagγ2sin2α , (2.58)
where Umag = B2/8π is the energy density of the magnetic field and β = v/c. This
result applies for electrons with an arbitrary pitch angle α. However, the pitch angle is likely to be randomised by irregularities in the magnetic field distribution or streaming instabilities (Longair, 2011), thus we expect the distribution of the pitch angles for a population of high-energy electrons to be isotropic. To take this into account, we average over an isotropic distribution of pitch angles p(α)dα = 12sin α, obtaining the average energy loss rate
− dE dt = 2 σTc v c 2 Umagγ2 1 2 Z π 0 sin3αdα = 4 3σTc Umagβ 2γ2. (2.59)
We next set out to calculate an expression for the total emissivity per unit frequency by an electron. To do this, we first need a description of the single-particle spectrum of an electron, which is set by a critical frequency that is located near the maximum of the spectrum (Blumenthal & Gould,1970),
νc=
3 2γ
2ν
gsin α, (2.60)
where νg is the non-relativistic gyrofrequency of the electron in the field given by
νg =
eB 2πmec
Figure 2.8: The spectrum of the SR of a single electron shown linear (left) and logarithmic (right) scales
Now, the SR emissivity by a single electron is
j(ν) = √ 3 2π e3B sin α mc2 F (x), (2.62) where F (x) = x Z ∞ x K5/3(ϕ)dϕ, (2.63)
with x = ν/νc and K5/3 is a modified Bessel function of order 5/3. This spectrum
has a broad maximum, centred roughly around ν ≈ νc and attains a maximum at
νmax= 0.29νc.
The function F (x) has asymptotic forms for high and low frequencies which are given by F (x) ∼ √4π 3Γ(13) x 2 1/3 , x 1, (2.64) F (x) ∼ π 2 1/2 e−xx1/2, x 1. (2.65)
The high-frequency emissivity of the electron is therefore given by an expression of the form
j(ν) ∝ ν1/2exp(−ν/νc) , (2.66)
which is dominated by the exponential cutoff when ν νc. It can be seen from the
above expressions that for low frequencies, the emissivity is given by an expression of the form
Next, we work out the emission spectrum for a power-law distribution of electron
en-ergies, N (E)dE = κE−pdE, where N (E)dE is the number density of electrons in the
energy range E to E + dE. Since the spectrum of SR is so sharply peaked near the crit-ical frequency νc, one can use the approximation that an electron of energy E radiates
away its energy at νc (Longair,2011). Therefore the energy radiated in the frequency
range ν to ν + dν can be attributed to electrons with energies in the range E to E + dE and so the emissivity
J (ν)dν = − dE
dt
N (E) dE. (2.68)
One can now obtain an a expression for the emissivity per unit volume
J (ν) ∝ ν−(p−1)/2∝ ν−a, (2.69)
where p is the index of the electron energy spectrum and a is the spectral index.
2.4.2 Inverse Compton Scattering
Inverse Compton scattering (IC) is a process in which a soft (low-energy) background target photon is upscattered to high energies through interaction with a high-energy electron. In what follows, I give a brief outline of this process.
The Thomson (elastic) limit is valid when γ mec2
1, (2.70)
where γ is the electron Lorentz factor, me is the mass of the electron, and is the
soft-photon energy. The mean energy of a Compton-scattered photon 1for an isotropic
photon gas is given by (Blumenthal & Gould,1970) h1i = 4
3γ
2hi , (2.71)
where hi is the mean energy of the soft photons. The total energy loss of a single electron is (Rybicki & Lightman,1979)
−dE dt = 4 3σTcγ 2β2U ph, (2.72)
where Uphis the photon energy density. Note the similarity with the expression for ˙ESR
