6.1 The magnetic field and transport equation

(1)

Transport in an Axisymmetric System

The previous chapter investigated the evolution of the particle spectrum in a spherically- symmetric central source system. The assumption inherent to these systems is that the particles only diffuse in the radial direction. However, in a more realistic system diffusion is determined by the geometry of the magnetic field. As the particles can diffuse in directions both perpendic- ular and parallel to the field lines, the particles will be transported not only in a radial direction but also in the polar and azimuthal directions. Furthermore, the geometry can also lead to the additional transport processes of gradient and curvature drift.

In the present chapter, the evolution of the particle spectrum is investigated for a system where parallel and perpendicular diffusion as well as drift can occur. This investigation requires an additional spatial dimension, and the transport equation is therefore solved in an axisymmetric geometry. The magnetic field that will be used for the simulations is the Archimedean spiral derived in Section 2.2.2. The magnetic field frozen into the plasma wind is believed to have this geometry in the region between the pulsar’s light cylinder and the pulsar wind termination shock, and based on this many PWN models assume that the magnetic field in the nebula has a similar structure (e.g. Rees and Gunn, 1974; Kennel and Coroniti, 1984a,b). This is supported by polarisation measurements of, e.g., PWNe G21.5-0.9, from which Zajczyk et al. (2012) were able to infer a toroidal structure for the inner regions of the nebula. As discussed in Section 2.2.3, the dipole nature of the pulsar’s magnetic field also leads to the formation of a neutral sheet that separates the magnetic field lines that are directed in opposite directions in the two hemispheres. The simulations therefore also take into account drifts along this neutral sheet.

This chapter commences with a brief introduction to the physics of diffusion and drift, before presenting modelling results. Apart from investigating the role of drift, the results will also show how the sign of the particle’s electric charge, as well as the orientation of the pulsar’s rotation and magnetic axes relative to each other, influence the evolution of the spectra. The aim of this chapter is not to present a detailed study, but rather an initial investigation into the possible effects that drift and polar diffusion may have on spectral evolution. As far as is known, this is the first time that such an investigation has been done for a central source system.

93

(2)

94 6.1. THE MAGNETIC FIELD AND TRANSPORT EQUATION

6.1 The magnetic field and transport equation

The following expressions and definitions can be found in Sections 2.2.2 and 2.2.3, but are repeated for the sake of convenience. The angle between the rotational Ω Ω Ω and magnetic µ µ µ axes is defined as

A ≡ cos α = Ω Ω Ω · µ µ µ

Ωµ , (2.18)

while the angle between the magnetic field and the radial direction is given by ψ = tan

⁻¹

( Ωr sin θ

V

_r

)

. (2.15)

Lastly, the magnetic field can be described by the equation B = B

∗

( r

∗

r )

2

(e

r

− tan ψe

ϕ

) [1 − 2H (θ − θ

^ns

)] , (2.20) where H represents the Heaviside step function and θ

ns

the polar extent of the wavy neutral sheet.

In Section 5.1 of the previous chapter, the transport equation

∂f

∂t + ∇ · S − 1 p

²

∂

∂p (

p

²

[

⟨ ˙p⟩f + D

p

∂f

∂p ])

= Q(r, p, t), (5.1)

was introduced, with the differential current density

S = 4πp

²

(CVf − K · ∇f) (5.2)

describing the convective and diffusive flux of particles in space. Here K(r, p, t) denotes the diffusion tensor.

6.2 Diffusion in the magnetic field

It is well-known that the Lorentz force will cause charged particles in a magnetic field B to gyrate around a point G, known as the guiding centre. Suppose that at any location in the system the magnetic field can be described by B = B

0

+ δB, where B

0

is the average magnetic field and δB variations resulting from turbulence. As a result of these variations, particles will diffuse in directions both parallel and perpendicular to the magnetic field, and in a coordinate system that is aligned with B, the diffusion tensor has the simplified form

K

_S

=







κ

_∥

0 0

0 κ

_⊥

0 0 0 κ

_⊥





 , (6.1)

where κ

_∥

describes diffusion parallel to the magnetic field, and κ

⊥

diffusion perpendicular to

the magnetic field. Note that (6.1) is based on the assumption of axisymmetric perpendicular

diffusion.

(3)

The following two sections provide the briefest of introductions to parallel and perpendicular diffusion: As such, the interested reader is referred to the discussions of, e.g., Shalchi (2009) and Engelbrecht (2013) for a more detailed treatment of these topics.

6.2.1 Parallel diffusion

Different particles will encounter an inhomogeneity in the magnetic field at different phases in their gyration. If the particles had a similar pitch-angle θ

p

before the inhomogeneity, the difference in gyro-phase leads to the particles having different pitch-angles with respect to the inhomogeneity. The result is that particles will either scatter in a forward direction, or will be reflected along the field line, depending on the value of their phase angle. The most effective scattering occurs when the scale length of the inhomogeneity is comparable to the gyro-radius of the particle (see, e.g., Dr¨oge, 2003), while particles with much larger or much smaller gyroradii will not be affected by the inhomogeneity. The process of parallel diffusion is fairly well-understood in terms of Quasilinear Theory, first introduced by Jokipii (1966).

6.2.2 Perpendicular diffusion

Fluctuations δB that are perpendicular to B

0

will lead to the diffusion of magnetic field lines.

As the guiding centres of the particles will follow these lines, the particles will subsequently diffuse perpendicular to B

0

. This picture is known as the Field Line Random Walk limit (see, e.g., Matthaeus et al., 1995) of Quasilinear Theory (Jokipii, 1966). It has however been found that this description is only accurate at higher particle energies, with the Nonlinear Guiding Centre Theory developed by Matthaeus et al. (2003) providing a more accurate picture over a wider energy range. Nevertheless, a precise description of perpendicular diffusion is still an outstanding problem of astrophysics, and a full discussion thereof is beyond the scope of the present study. For an excellent introduction to this topic, the interested reader is referred to Shalchi (2009).

6.3 Drift in the magnetic field

6.3.1 Gradient and curvature drift

If the pitch angle of the particle is θ

p

< 90

^◦

, the guiding centre G will have a velocity component V

_∥

that is parallel to B. Additionally, gradients and curvatures in the magnetic field will cause the guiding centre to drift with a velocity V

D

in a direction that is perpendicular to B.

Averaging over both pitch and phase angle, the total velocity of the guiding centre can be written as

⟨V

G

⟩ = ⟨V

∥

⟩ +

⟨ pv 3qB

( B × ∇B

B

²

+ ∇ × B B

)⟩

, (6.2)

(4)

96 6.3. DRIFT IN THE MAGNETIC FIELD

D D

B f

d

³

x <v

G

>

<v

_f

>

Figure 6.1: The right-hand side of the figure shows a negatively charged particle drifting as a result of a gradient in the magnetic field. The left hand-side shows the trajectories of the particles passing through the volume element d³x. It follows from this illustration that ⟨Vf⟩ and ⟨VG⟩ are not equivalent.

where the second and third terms on the right-hand side describe gradient and curvature drifts.

The derivation of ⟨V

G

⟩ is a standard result of plasma physics, and can be found in textbooks such as those of Rossi and Olbert (1970) and Chen (1985). Note that this derivation requires that the scale length of the variations in the magnetic field must be significantly larger than the gyroradius of a particle.

In a system where particles are transported by convection, diffusion, and drift, the differential streaming density is given as

S = 4πp

²

(CVf − K · ∇f − ⟨V

f

⟩f) , (6.3)

where the last term describes the flux of particles through the volume element d

³

x as a result of drift. While (6.2) is an exact description of the motion of G for an average particle, it is important to note that ⟨V

f

⟩ and ⟨V

G

⟩ are not equivalent. This is illustrated in Figure 6.1, where it can be seen that the gradient in B causes an electron to drift towards the right. However, as a result of a gradient in the distribution function f (x, p, t), the net flux of particles in the volume d

³

x may be in a completely different direction.

It was first shown by Jokipii and Levy (1977), and later by Burger (1987), that the pitch- and phase-angle averaged velocity of the particles located in the momentum interval p + dp is given by

⟨V

f

⟩ = ⟨V

f,∥

⟩ + ⟨ κ

_⊤

f

B

B × ∇f ⟩

, (6.4)

where

κ

_⊤

= βpc

3qB . (6.5)

This expression is exact even if the distribution function has first-order anisotropies f = f

0

(1 + cos θ

_p

).

However, for the purpose of this study f is taken as isotropic, and consequently ⟨V

f,∥

⟩ = 0.

For an isotropic distribution of particles one has

∇ · (⟨V

f

⟩f) = V

D

· ∇f, (6.6)

where

V

_D

= ∇ × κ

⊤

B

B (6.7)

(5)

is the average drift velocity as a result of gradient and curvature drifts, as discussed by Jokipii and Levy (1977) and Burger (1987). Using Liouville’s theorem, Fisk (1976) demonstrated that

∇ · V

D

= 0. The corollary of this result is that drift will only transport particles between different parts of the system, and that the streamlines of V

D

must be closed.

In an axisymmetric system described using the spherical coordinates (r, θ, ϕ), the components of (6.7) are

V

_D,r

= κ

_⊤

rB

(

cot θB

_ϕ

+ ∂B

_ϕ

∂θ − 2 B

_ϕ

B

∂B

∂θ )

, (6.8a)

and

V

_D,θ

= κ

_⊤

B

( 2 B

_ϕ

B

∂B

∂r − B

_ϕ

r − ∂B

_ϕ

∂r )

. (6.8b)

When the the magnetic field has an Archimedean spiral geometry (2.20), these components reduce to

V

_D,r

= κ

_⊤

r

(

− sin ψ

tan θ + cos ψ (2 sin

²

ψ − 1 ) ∂ tan ψ

∂θ )

, (6.9a)

and

V

_D,θ

= κ

_⊤

(

cos ψ (1 − 2 sin

²

ψ ) ∂ tan ψ

∂r + 3 sin ψ r

)

. (6.9b)

While (6.4) can be used to formally include gradient and curvature drift into the transport equation, Axford (1965) and Parker (1965) have shown that it is also possible to take this trans- port process into account by modifying the diffusion tensor in the reference frame of the mag- netic field to become K = K

_S

+ K

_A

, where the symmetric (diffusion) tensor is given by (6.1), and the asymmetric (drift) tensor by

K

A

=







0 0 0

0 0 κ

_⊤

0 −κ

⊤

0 



 . (6.10)

It can be shown (see, e.g., Steenkamp, 1995) that

∇ · (K

A

· ∇f) = V

D

· ∇f (6.11)

if K

_A

has the form prescribed by (6.10), and the differential streaming density (6.3) can there- fore be alternatively formulated as

S = 4πp

²

(CVf − [K

S

+ K

_A

] · ∇f) . (6.12) It is this formulation that will be used for the numerical simulations.

6.3.2 Neutral sheet drift

Apart from gradient and curvature drift, particles within two gyroradii 2R

0

of the neutral

sheet will also be subjected to an additional drift motion. Figure 6.2 shows the trajectory of

two positively charged particles with the same momenta. It can be seen in the figure that

(6)

98 6.3. DRIFT IN THE MAGNETIC FIELD

Figure 6.2: Illustration showing two positively charged particles with the same momenta drifting along the neutral sheet (solid horizontal line). The figure is taken from Ferreira (2002).

the direction of gyration changes when the particles cross the sheet, leading to the particles drifting along the sheet. The particle close to the neutral sheet completes only a small fraction of its gyro-orbit before the gyration direction changes. On the other hand, the particle located further away completes a larger fraction of its gyro-orbit.

Figure 6.3: The neutral sheet drift velocity (VD,ns) profile calculated by Burger (1987). The neutral sheet is located along the y-axis, while the drift velocity has been normalised to its maximum value of VD,ns= 0.463v.

Burger (1987) calculated that the maximum drift velocity along the neutral sheet is V

D,ns

= 0.463v, where v is the velocity of the particle. Furthermore, this author also calculated that V

_D,ns

decreases as a function of distance from the neutral sheet, as shown in Figure 6.3. The area under the curve represents the average neutral sheet drift velocity, with Burger (1987) finding that

⟨V

D,ns

⟩ = 4R

0

v

6 = 2κ

_⊤

, (6.13)

with κ

⊤

defined in (6.5). For the simulations the velocity curve is approximated using a tri-

angle, and is sufficient provided that the areas under the curve and approximation are the

same. To include neutral sheet drift in the problem, the components V

D,ns,r

and V

D,ns,θ

at a

given r and θ position are calculated using this approximation, with the results added to the

corresponding velocity components (6.9a) and (6.9b) describing gradient and curvature drift.

(7)

Figure 6.4: Average drift velocity streamlines in the meridional plane for the qA > 0 scenario, as illus- trated by Pesses et al. (1981). The streamlines are the same for the qA < 0 scenario, but in the opposite direction.

6.3.3 Drift motion in an axisymmetric system

Research focusing on the modulation of cosmic rays in the heliosphere has shown that apart from the electric charge of the particle q, the orientation of Ω relative to µ also plays a role in determining the drift motion in the system, and it follows from (2.18) that one can distinguish between an A < 0 and an A > 0 orientation. For the former scenario, the angle between Ω and µ is α > 90

^◦

, while the latter scenario has α < 90

^◦

. In an analogous manner, it then becomes possible to distinguish between two types of drift scenarios, depending on whether qA < 0, or qA > 0.

Figure 6.4 illustrates the average drift velocity streamlines in the meridional plane for a qA > 0 heliosphere. As an example, a positively charged particle in an A > 0 orientation will drift in towards the equatorial region, as well as out along the neutral sheet. By contrast, a negatively charged particle in an A > 0 orientation will drift in along the neutral sheet, and then towards the poles. As the primary difference between the system studied here and the heliosphere is related to the location of the particle source, the drift motion for the two systems should be similar in nature.

6.4 The two-dimensional, steady-state transport model

For an an axisymmetric steady-state system where the particles are convected outward by a

purely radial wind V

r

, while simultaneously diffusing, drifting, and losing energy as a result

(8)

100 6.4. THE TWO-DIMENSIONAL, STEADY-STATE TRANSPORT MODEL

of adiabatic losses, the general transport equation (5.1) becomes κ

_rr

∂

²

f

∂r

²

+ κ

_θθ

r

²

∂

²

f

∂θ

²

+ [ 1

r

²

∂

∂r (r

²

κ

_rr

) + 1 r sin θ

∂

∂θ (sin θκ

_θr

) − V

r

] ∂f

∂r + [ 1

r

²

∂

∂r (rκ

_rθ

) + 1 r

²

sin θ

∂

∂θ (sin θκ

_θθ

) ] ∂f

∂θ +

[ 1 3r

²

∂

∂r (r

²

V

r

) ] ∂f

∂ ln p + Q (r, θ, p, t)

= 0,

(6.14)

with the expressions for the diffusion coefficients supplied by the transformation of K = K

_S

+ K

_A

from field-aligned to spherical coordinates,

K (r, θ, ϕ) =







κ

_rr

κ

_rθ

κ

_rϕ

κ

_θr

κ

_θθ

κ

_θϕ

κ

_ϕr

κ

_ϕθ

κ

_ϕϕ







=







cos ψ 0 sin ψ

0 1 0

− sin ψ 0 cos ψ













κ

_∥

0 0

0 κ

_⊥

κ

_⊤

0 −κ

⊤

κ

_⊥













cos ψ 0 − sin ψ

0 1 0

sin ψ 0 cos ψ







=







κ

_∥

cos

²

ψ + κ

_⊥

sin

²

ψ −κ

⊤

sin ψ (κ

_⊥

− κ

∥

) cos ψ sin ψ

κ

_⊤

sin ψ κ

_⊥

κ

_⊤

cos ψ

(κ

⊥

− κ

∥

) sin ψ cos ψ −κ

⊤

cos ψ κ

⊥

cos

²

ψ + κ

_∥

sin

²

ψ





 .

(6.15)

For a more detailed discussion on the derivation of (6.14) from (5.1), Appendix B can be con- sulted.

Comparing (6.6) with the coefficients of the first derivatives appearing in (6.14) shows that the components of V

D

are given by

V

_D,r

= 1 r sin θ

∂

∂θ (sin θκ

_θr

) , (6.16a)

and

V

_D,θ

= 1 r

²

∂

∂r (rκ

_rθ

) . (6.16b)

Note that these expressions are equivalent to those given by (6.9a) and (6.9b), respectively.

As was the case for the time-dependent one-dimensional transport equation (5.8), the steady-

state axisymmetric transport equation (6.14) is again solved using the Douglas ADI method,

with the numerical scheme derived in Appendix C.6.3. The numerical model is adapted from

the one originally developed by Moraal (personal communication), where the latter model has

been used extensively by, e.g., Moraal and Gleeson (1975), Potgieter and Moraal (1985), and Rei-

necke et al. (1993) to model the modulation of cosmic rays in the heliosphere.

(9)

6.4.1 Parameter values and radial profiles used

For the axisymmetric simulations, parameters similar to those described in Section 5.2.1 are used. The unscaled convection velocity at the termination shock r

0

is again chosen to be V ˜

₀

= 0.3c, while the magnetic field is chosen to be ˜ B

₀

= 35 µG. This smaller magnetic field, compared to the value of ˜ B

₀

= 350 µG used in the previous chapter, is chosen in order to better illustrate the effect of drift. The smaller value chosen is nevertheless still reasonable within the context of pulsar wind nebulae. In contrast to the one-dimensional simulations of the previ- ous chapter, where it is implicitly assumed that diffusion only takes place perpendicular to the magnetic field lines, the present simulations require the specification of both a parallel (κ

_∥

) and perpendicular (κ

⊥

) diffusion coefficient. In order to have comparable results, the magnitude of the perpendicular diffusion coefficient for the axisymmetric simulations is chosen to be similar to that of the one-dimensional simulations, i.e. κ

⊥,0

= 4 × 10

²⁴

cm

²

s

⁻¹

. For the parallel diffu- sion coefficient a value that is a hundred times larger (see, e.g., Jokipii and Kota, 1995; Giacalone and Jokipii, 1999) is used in the simulations, i.e. κ

_∥,0

= 4 × 10

²⁶

cm

²

s

⁻¹

.

For the simulations a purely radial convection velocity with the spatial dependence V (r) ∝ 1/r is used, and it follows from (2.15) and (2.20) that this leads to a radially independent magnetic field B(r) = B

0

. It has been demonstrated in Chapter 3, as well as by Kennel and Coron- iti (1984a,b), that the magnetic field can either increase or decrease as a function of position, depending on the ratio of electromagnetic to particle energy σ. A radially independent mag- netic field is therefore useful as this represents an intermediate scenario. It will, however, be recalled from Section 5.2.5 that the radial dependence of the parameters do not qualitat- ively influence the spectral evolution. Note that the large angular velocity of the pulsar Ω in (2.15) effectively implies that ψ = 90

^◦

for all values of r and θ. The only exception is at the poles (θ = 0

^◦

, 180

^◦

) where ψ = 0

^◦

. Here the magnetic field has only a radial component that decreases as B

r

(r) ∝ 1/r

²

. Following Section 5.2.1, it is again assumed that the spatial and mo- mentum dependence of the diffusion and drift coefficients appearing in the diffusion tensor (6.15) can be separated.

Apart from choosing values for the convection velocity and diffusion coefficient that are sim- ilar to those used for the simulations of Chapter 5, the radial profiles used for the convection velocity, magnetic field, and by implication the diffusion coefficient, are also similar to those used for the scenario shown in Figure 5.3(a) of Section 5.2.4. The only difference between the two scenarios is that the present simulations use a magnetic field that is ten times smaller.

However, as synchrotron losses are neglected in the present simulations, while the diffusion

coefficient is scaled arbitrarily, the value of B

0

plays no role in the evolution of the spectra

when drift is neglected. Therefore, when drift is not taken into account it is possible to directly

compare the present results with those of Section 5.2.4, and in particular with those shown in

Figure 5.3(a).

(10)

Figure 6.5: The radial (a) and angular dependence (b) of the diffusion coefficient κrr.

6.4.2 Numerical considerations

Similar to the steady-state model of Section 5.2, momentum is again used as the stepping para- meter in the present scheme, along with an identical initial condition. The momentum range is chosen to be 10

⁻³

≤ p ≤ 10

⁴

with a logarithmic step size of ∆ ln p = −0.02.

In the first half of the momentum step, the numerical scheme is solved implicitly in the radial direction using the same boundary conditions described in Section 5.2.2, with the radial grid ranging between r

0

≤ r ≤ 1, and ∆r = 0.002.

In the second half of the momentum step the scheme is solved implicitly in the polar direction over the interval 0

^◦

≤ θ ≤ 180

^◦

using a step size ∆θ = 0.5

^◦

, with a Neumann condition

∂f

∂θ = 0

θ=0^◦,180^◦

(6.17)

imposed at both boundaries. This boundary condition is based on the assumption that there is a symmetry between the eastern and western hemispheres. Note that the boundaries in the θ-direction coincide with the positions of the poles, and that the equator is located at θ = 90

^◦

.

6.4.3 Radial and angular dependence of the transport coefficients

Before turning to the evolution of the particle spectra, it is worthwhile to first investigate the radial and angular dependence of the coefficients present in (6.14). When interpreting these results, it should be kept in mind that the diffusion and drift velocities scale with momentum p, and that the graphs in Figures 6.5 and 6.6 are specifically for the case p = 1.

The first important coefficient is κ

rr

, as defined in (6.15). Figure 6.5(a) shows the radial depend- ence at three polar angles, and Figure 6.5(b) the angular dependence at three radial distances.

As ψ ≈ 90

^◦

for 0

^◦

< θ < 180

^◦

, it follows from (6.15) that diffusion in the radial direction is

essentially perpendicular to the magnetic field. Figure 6.5(a) shows that at these polar angles

(11)

κ

_rr

is independent of radial position, reflecting the fact that the components of (6.15) scale as κ

_i,j

∝ 1/B. At the poles where ψ = 0

^◦

, only parallel diffusion is important, and the magnetic field B

r

∝ 1/r

²

leads to the radial dependence κ

rr

∝ r

²

, as shown in Figure 6.5(a). Due to the choice of parameters κ

∥

= 100κ

_⊥

, Figure 6.5(b) shows that κ

rr

decreases rapidly when θ is small, before approaching a constant value at higher polar angles. It can be seen from (6.15) that κ

θθ

= κ

_⊥

, and the polar diffusion coefficient will therefore have a radial and angular dependence that is similar to that of κ

rr

when 0

^◦

< θ < 180

^◦

.

The next set of quantities to be investigated are the drift velocities V

D,r

and V

D,θ

, as described by (6.16a) and (6.16b) respectively. These velocities are a part of the coefficients associated with the first derivatives that appear in the transport equation (6.14), along with terms that contain the derivatives ∂κ

θθ

/∂θ and ∂κ

rr

/∂r. It was found that the ∂κ

θθ

/∂θ-term is generally signific- antly smaller than the other terms, and will therefore not be investigated. As the ∂κ

rr

/∂r-term appears in a coefficient together with the drift and convection velocities V

D,r

and V

r

, this term can be interpreted as a ”diffusion velocity” V

κ,r

. This coefficient determines the rate at which a density feature is transported in the radial direction without a dispersion in shape.

The left-hand panels of Figure 6.6 show the radial dependence of the velocities V

D,r

, V

D,θ

, and V

_κ,r

at three angular positions, while the right-hand panels show the angular dependence of these velocities at three radial positions. Shown for comparison is V

r

at the same radial and angular positions.

Figure 6.6(a) shows that at θ = 0

^◦

both V

D,r

and V

κ,r

increase as function of radial distance, while V

D,θ

= 0. This is again related to the fact that B ∝ 1/r

²

at the poles. By contrast, Figures 6.6(c) and (e) show that at higher polar angles V

D,θ

and V

D,r

decrease as a function of position.

Due to the purely azimuthal structure of the magnetic field, drift at higher polar angles is essentially in the polar direction, and consequently V

D,r

becomes small, as seen in Figure 6.6(c), and qualitatively shown in Figure 6.4. The exception is in the equatorial region, as shown in Figure 6.6(e), where neutral sheet drift becomes important. This can also be seen from the angular profiles shown in Figures 6.6(b), (d), and (f), where a peak related to neutral sheet drift appears at θ ∼ 85

^◦

. This implies that neutral sheet drift (6.13) is numerically distributed over a region that extends ±5

^◦

from the neutral sheet. Lastly, the angular profiles also show that all velocities are essentially independent of polar angle.

For most radial distances and polar angles the drift velocities are smaller than the convective

velocity, implying that drift will have only a moderate effect on spectral evolution. However, at

the poles and near the equator the parameters are chosen such that drift becomes the dominant

transport process.

(12)

Figure 6.6: The radial, (a),(c),(e), and angular dependence, (b),(d),(f), of the drift (VD,r and VD,θ) and diffusion velocities (Vκ,r), along with the radial convection velocity (Vr).

(13)

6.5 Results

6.5.1 Solutions to the transport equation

Solutions to (6.14) were calculated for a no-drift scenario, as well as for the qA < 0 and qA > 0 scenarios defined by (2.18). The results are shown in Figure 6.7, where spectra are plotted at selected r and θ values. These solutions are specifically for the case where the magnetic and rotation axes of the pulsar are aligned, while the spectra for an oblique rotator will be discussed in Section 6.5.2. Comparison of the no-drift spectra in Figures 6.7(c) and (i) with those of Figure 5.3(a) shows that both the one- and two-dimensional models predict the same evolution of the particle spectrum, thereby indicating that the two-dimensional model works correctly (see last paragraph in Section 6.4.1).

The qA < 0 scenario

It follows from Figures 6.5 and 6.6 that, apart from the polar and equatorial regions, the trans- port processes are independent of polar angle, so that the spectra at θ = 45

^◦

in Figures 6.7(b), (e), and (h) are the spectra that one would typically expect in the largest part of the system.

Figure 6.7(b) shows that the no-drift and qA < 0 spectra in the inner part of the system (r = 0.1) are indistinguishable at lower energies ( ˜ E ≲ 0.7 TeV), implying that convection is the dominant transport processes in this energy regime. At higher energies it can be seen that in the energy range 0.7 TeV ≲ ˜ E ≲ 200 TeV the qA < 0 spectrum is initially softer than that of the no-drift scenario, before hardening again, leading to the formation of a ”trough”-like feature (it should be noted that this nomenclature is not entirely correct as the spectra have been multiplied by E ˜

²

to enhance spectral features). At ˜ E ≳ 200 TeV the qA < 0 spectrum has a similar index to that of the no-drift spectrum, although the intensity of the former is marginally lower than that of the latter. Similar spectral features are also seen at larger radial distances, as shown in Figures 6.7(e) and (h). A closer inspection of these figures does however show that the trough-like spectral feature grows with radial distance.

Spectral features similar to those described above also develop at the poles (θ = 0

^◦

), as shown

in Figures 6.7(a), (d), and (g). The main difference is that the qA < 0 spectrum at the lowest

energies ( ˜ E ≲ 0.2 TeV) has a higher intensity and larger spectral index compared to the no-drift

scenario. As drift scales with momentum p, one would not expect this transport process to lead

to new spectral features at the lowest energies, and it is therefore not clear if this behaviour can

be ascribed to a physical process, or whether it results from the choice of boundary conditions

for the r and/or θ-grid. Lastly, Figures 6.7(c), (f), and (i) show that the trough-like feature

is especially pronounced in the equatorial region (θ = 90

^◦

), and it is only at the very lowest

energies ( ˜ E ≲ 0.02 TeV) that the qA < 0 spectra resemble the spectra of the no-drift scenario.

(14)

106 6.5. RESULTS

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

θ = 0°

(a)

r=0.1

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

θ = 0°

(a)

r=0.1

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

θ = 0°

(a)

r=0.1

θ = 45°

(b) no drift

θ = 45°

(b) no drift

qA < 0

θ = 45°

(b) no drift

qA < 0 qA > 0

θ = 90°

(c)

θ = 90°

(c)

θ = 90°

(c)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

r=0.3

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

r=0.3

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

r=0.3

(e) (e)

(e) (f)(f)(f)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

r=0.9

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

r=0.9

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

r=0.9

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

Figure 6.7: The evolution of the spectra for a steady-state, axisymmetric system. The spectra are shown at selected r and θ values, as indicated on the figures. In every figure three spectra are shown, corresponding to the no-drift, qA < 0, and qA > 0 scenarios.

(15)

The qA > 0 scenario

Figure 6.7(b) shows that the qA > 0 spectrum in the inner part of the system (r = 0.1) becomes softer than that of the no-drift scenario at ˜ E ∼ 1 TeV, leading to a larger discrepancy in intensity between the qA > 0 and no-drift spectra at higher energies. At ˜ E ≳ 200 TeV the spectral index of the qA > 0 scenario is again similar to that of the no-drift scenario, as well as to the spectral index of the qA < 0 scenario. It is also seen from Figure 6.7(a) that the intensity of the qA > 0 spectrum is lower than that of both the no-drift and qA < 0 scenarios. This spectral behaviour can also be seen at larger radial distances, as shown in Figures 6.7(e) and (h), as well as in the polar regions, as shown in Figures 6.7(a), (d), and (g). A closer inspection of these last three figures shows that the intensity of the qA > 0 spectra for energies ˜ E ≲ 0.2 TeV is marginally lower than those of the no-drift spectra. It is again not clear whether this discrepancy has a physical or numerical origin.

As a very prominent trough-like spectral feature develops in the qA < 0 scenario, one would expect that a prominent spectral feature should also develop for the qA > 0 scenario. Such a feature does indeed develop, but is limited to the equatorial region. Compared to the no-drift scenario, Figures 6.7(c), (f) and (i) show that at lower energies ( ˜ E ≲ 10 TeV) the qA > 0 spectra are harder, while being softer at larger energies. This leads to the formation of a ”peak”-like structure that appears to be the inverse of the trough-like feature of the qA < 0 scenario.

Figure 6.6(e) shows that the velocity associated with neutral sheet drift is larger than the other velocities, and as both the trough- and peak-like features are most prominent in the ecliptic plane where the neutral sheet is located, it is natural to associate these spectral features with neutral sheet drift. Therefore, in the qA < 0 scenario the particles will drift inward along the neutral sheet, making it more difficult for the particles to propagate to larger radial distances, while the exact opposite occurs for the qA > 0 scenario.

6.5.2 The effect of neutral sheet drift and an oblique rotator

To illustrate that the trough- and peak-like spectral features develop as a result of neutral sheet drift, simulations similar to those of Figure 6.7 were performed, but with the neutral sheet drift velocity set to zero. It should be kept in mind that this represents a non-physical scenario as it leads to ∇ · V

D

̸= 0 (see Section 6.3.1). This scenario is therefore used purely for the purposes of illustration. The results are shown in the middle rows of Figures 6.8 and 6.9. The former figure shows the spectra at various angular positions and at r = 0.1, while the latter figure shows the spectra at the same angular positions, but for r = 0.9.

These figures show that when neutral sheet drift is absent, the spectra of the qA < 0 scenario

are very similar to those of the no-drift scenario. The only significant difference occurs at the

poles and for energies ˜ E ≲ 0.2 TeV. For the qA > 0 scenario the spectra at lower energies

( ˜ E ≲ 1 TeV) are similar to those of the no-drift scenario, both in spectral index and intensity,

(16)

108 6.5. RESULTS

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

r=0.1, θ=0°

(a)

Fig 6.7 (α=0°)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

r=0.1, θ=0°

(a)

Fig 6.7 (α=0°)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

r=0.1, θ=0°

(a)

Fig 6.7 (α=0°)

r=0.1, θ=45°

(b) no drift

r=0.1, θ=45°

(b) no drift

qA < 0

r=0.1, θ=45°

(b) no drift

qA < 0 qA > 0

r=0.1, θ=90°

(c)

r=0.1, θ=90°

(c)

r=0.1, θ=90°

(c)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

V

_sheet

=0

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

V

_sheet

=0

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

V

_sheet

=0

(e) (e)

(e) (f)(f)(f)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

α =30 °

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

α =30 °

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

α =30 °

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

Figure 6.8: The spectra at r = 0.1 for the following three scenarios: a system with a flat neutral sheet (top row), a system without a neutral sheet (middle row), and a system that again has a neutral sheet, but where the angle of inclination between the rotation and magnetic axes of the pulsar is α = 30^◦(bottom row). The top row of figures are the same as those in Figure 6.7, but have been repeated here for ease of reference.

(17)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

r=0.9, θ=0°

(a)

Fig 6.7 (α=0°)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

r=0.9, θ=0°

(a)

Fig 6.7 (α=0°)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

r=0.9, θ=0°

(a)

Fig 6.7 (α=0°)

r=0.9, θ=45°

(b) no drift

r=0.9, θ=45°

(b) no drift

qA < 0

r=0.9, θ=45°

(b) no drift

qA < 0 qA > 0

r=0.9, θ = 90°

(c)

r=0.9, θ = 90°

(c)

r=0.9, θ = 90°

(c)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

V

_sheet

=0

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

V

_sheet

=0

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

E ~

2

N (norm. units)

(d)

V

_sheet

=0

(e) (e)

(e) (f)(f)(f)

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

α =30 °

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

α =30 °

10

^-10

10

^-8

10

^-6

10

^-4

10

^-2

10

⁰

10

²

10

⁴

10

⁶

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

2

N (norm. units)

E ~

=p ~c (TeV)

(g)

α =30 °

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(h)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

10

^-3

10

^-2

10

^-1

10

⁰

10

¹

10

²

10

³

E ~

=p ~c (TeV)

(i)

Figure 6.9: The same as Figure 6.8, but for r = 0.9.

(18)

110 6.6. SUMMARY

while having a similar index, but lower intensity, at higher energies. As neutral sheet drift is absent, it thus follows that this deviation from the no-drift spectra is a result of gradient and curvature drift.

Additionally, the effect that a larger tilt angle α of the neutral sheet has on the evolution of the spectra also has to be investigated. In a more realistic scenario the rotation and magnetic axes of the pulsar will not be aligned, resulting in a wavy neutral sheet. Therefore, simulations using α = 30

^◦

were also performed, with the results plotted in the bottom rows of Figures 6.8 and 6.9.

Figure 6.8 shows that at r = 0.1 the α = 0

^◦

and α = 30

^◦

scenarios have similar qA < 0 spectra.

Figure 6.9 shows that this is also true for larger radial distances (r = 0.9), with the exception that the intensity of the α = 30

^◦

spectra are reduced for energies ˜ E ≳ 200 TeV. The consequence is that the trough-like feature is less pronounced, particularly in the equatorial region.

As was the case for qA < 0, the qA > 0 spectra for the α = 0

^◦

and α = 30

^◦

scenarios are very similar in the inner part of the system (r = 0.1), as shown in Figure 6.8. Figure 6.9 shows that this is also true at larger radial distances (r = 0.9) and for angular positions that are not in the equatorial region. At θ = 90

^◦

the peak-like feature completely disappears when α = 30

^◦

. This is in contrast to the qA < 0 spectrum, where a diminished trough-like feature is still visible when α = 30

^◦

.

When α ̸= 0, the wavy neutral sheet leads to an alternating drift of particles in the polar and equatorial directions along the neutral sheet (see Figure 6.4). This results in an increase of the polar component of the neutral sheet drift, while the radial component decreases. The effect of neutral sheet drift is thus reduced, as the same amount of particles that drift in a polar direction will be equal to the number of particles that drift in an equatorial direction. Consequently the trough- and peak-like structures associated with neutral sheet drift (cf. middle rows of Figures 6.8 and 6.9) become less pronounced.

6.6 Summary

This chapter investigated the effects that parallel and perpendicular diffusion, as well as drift, have on the evolution of the particle spectrum. These processes depend on the structure of the magnetic field, and for the simulations presented in this chapter a dipole Archimedean spiral geometry was used. Apart from the polar regions, where the magnetic field has only a radial component, the large angular velocity of the pulsar leads to a magnetic field that is almost purely azimuthal, a specific field geometry often assumed in one-dimensional PWN models (see, e.g., Kennel and Coroniti, 1984a,b; Van der Swaluw, 2003b). Due to the geometry of the field, particles will be subjected to gradient and curvature drift, while the dipole nature of the field leads to the presence of a neutral sheet, which in turn leads to drift along this structure.

As a study of this nature has never been performed for a central source system, this chapter

(19)

is intended as an initial investigation. In order to better illustrate the effects of the transport processes, synchrotron losses were neglected in the simulations.

The first set of simulations focused on a system where the magnetic and rotation axes of the pulsar are aligned (i.e., for a tilt angle α = 0

^◦

). In this specific scenario the system has a flat neutral sheet that is confined to the equatorial plane. When drift is neglected in the simula- tions, the spectral evolution is similar to that calculated in Section 5.2.4. When drift is included, it was found that the qA < 0 spectra, as defined by (2.18) in Section 6.1, develop a ”trough”-like feature at all radial distances and polar angles, with this feature being particularly pronounced in the equatorial region (see Figure 6.7). By contrast, a ”peak”-like structure that appears to be the inverse of the structure seen for the qA < 0 scenario, develops in the equatorial spectra of the qA > 0 scenario. However, for larger angular positions the evolution of the qA > 0 spectra is very similar to a scenario where only convection, adiabatic losses and diffusion occur.

Based on the result that the trough- and peak-like structures are most prominent in the equat- orial region, a natural conclusion is that they are the result of neutral sheet drift, and it was found that setting the neutral sheet drift velocity artificially to zero causes these features to disappear (see the middle rows of Figures 6.8 and 6.9). Furthermore, the qA < 0 spectra are very similar to those of the no-drift scenario, while the qA > 0 spectra have a lower intensity at higher energies compared to the no-drift spectra. These results lead one to tentatively con- clude that polar drift and diffusion are in opposite angular directions when qA < 0, causing these two processes to effectively cancel each other. On the other hand, these processes are in the same direction when qA > 0, leading to a deviation from the no-drift spectra. The effect is, however, limited.

As a more realistic PWN should have α > 0

^◦

, a scenario where α = 30

^◦

was also investigated (see the bottom rows of Figures 6.8 and 6.9). It was found that the trough-like feature still develops in the qA < 0 spectra, although it is less pronounced compared to the α = 0

^◦

scenario.

Away from the equatorial region the qA > 0 spectra in the α = 0

^◦

and α = 30

^◦

scenarios are similar: The peak-like structure completely disappears from the equatorial qA > 0 spectra when α = 30

^◦

, and they resemble those at smaller polar angles.

Due to the fact that both electrons and positrons are present in equal numbers in a PWN, the

total emission observed from a nebula will be the sum of the qA < 0 and qA > 0 spectra. Figure

6.7 shows that, apart from the equatorial region where neutral sheet drift is important, it will

be difficult to observe any spectral changes resulting from drift when α is small. Furthermore,

α is expected to be large for pulsars, and it follows from Figures 6.8 and 6.9 that any spectral

features related to neutral sheet drift will also be diminished. Nevertheless, drift may still be

important when synchrotron losses are taken into account, particularly if the magnetic field

has an angular dependence. As shown in Figure 6.4, particles with a specific electric charge

will drift to regions with a strong magnetic field, while the oppositely charged particles will

drift into a region with a weaker magnetic field. The result is that the synchrotron loss rate for

electrons and positrons will differ, potentially leading to large differences in the spectra.

(20)

112 6.6. SUMMARY