The Evolution of Particle Energy Spectra in Spherically-Symmetric Systems

The Evolution of Particle Energy Spectra in Spherically-Symmetric Systems

The previous chapter introduced a spatially independent model that can be used to calculate the evolution of a particle spectrum in a pulsar wind nebula. While this model includes the effects of adiabatic energy changes and non-thermal radiation losses, it relies on using the velo- city of the outer boundary to approximate adiabatic energy changes, and an average magnetic field to calculate the synchrotron loss rate. Another limitation of this model is that the effect of diffusion can only be included in a very limited fashion, as this process is only used to calcu- late the time needed for the particles to escape from the nebula. Lastly, a spatially independent model implies that the particle source is distributed throughout the nebula, contrary to the nature of pulsar wind nebulae where the particle source is located at the centre of the system.

Pulsar wind nebulae are just one example of central source systems characterised by the ac- celeration of particles in the innermost regions of the system, with this accelerated, energetic particle population transported outward by means of various processes. Other examples of central source systems include globular clusters (see, e.g., Venter et al., 2009), and systems where cosmic rays accelerated at supernova remnants interact with molecular clouds located some distance away from said remnants (Gabici and Aharonian, 2007). Although the model to be presented in this chapter and the next is developed within the context of pulsar wind neb- ulae, the main features of the model are equally applicable to all central source systems. For a pulsar wind nebula the particles in question are typically electrons and positrons, and some of the processes included in the model, specifically synchrotron and inverse Compton losses, are only important for these particles. However, the other processes included in the model are equally important for nuclei, or even partially ionised atoms.

A specific focus of this chapter is to investigate the effect that diffusion has on the evolution of the particle spectrum. Particle evolution models developed for pulsar wind nebulae, like those of Zhang et al. (2008) and Tanaka and Takahara (2010), are in general spatially independent, and diffusion is only used to estimate the time needed for particles to escape from the system. As

61

far as is known, only two spatially dependent models that allow for an improved treatment of diffusion have thus far been presented for pulsar wind nebulae. The first is the model of Lerche and Schlickeiser (1981), where a Fokker-Planck transport equation is solved analytically in cylindrical coordinates. Although originally developed for any central source system, this model requires that the convection velocity must either increase as a function of distance from the source or remain constant, in order to make the analytical approach tractable. Apart from this limitation, the model by Lerche and Schlickeiser (1981) has an additional characteristic that follows from the chosen geometry. As this model has only one spatial dimension, specifically z, particles will only suffer adiabatic losses when the velocity has the aforementioned character.

As shown in the previous chapter, this is contrary to the nature of pulsar wind nebulae where the convection velocity decreases as a function of distance. Due to the fact that adiabatic losses play a role in the evolution of the particle spectra (see, e.g., Del Zanna et al., 2006; Bucciantini, 2008), more realistic velocity profiles must be investigated.

The second spatially dependent model presented by Van Etten and Romani (2011) places no such restriction on the convection velocity, where these authors calculated the evolution of the particle spectrum using a method similar to the one discussed in Section 4.2. In their model, the time required to transport particles with a given energy between two regions in the nebula by convection and/or diffusion is first calculated, thereby supplying the spatial evolution of the spectrum. With this transport time known, it is possible to determine the energy losses that the particles have suffered during the propagation time, and the spectrum is accordingly evolved in energy. As these authors focused specifically on modelling the non- thermal emission observed from the nebula HESS J1825–137, they did not explicitly discuss the effect of diffusion on the particle spectrum.

The main focus of the work presented in this chapter and the next is thus to develop a spatially dependent model based on a Fokker-Planck transport equation, as well as investigating the ef- fect of diffusion on the evolution of the particle spectrum. The transport equation used in this work is similar to the transport equation used by Ginzburg and Syrovatskii (1964) to describe the propagation of cosmic rays in the Galaxy, and the transport equation used by Parker (1965) to describe the modulation of cosmic rays in the heliosphere. The spatial transport processes that will be treated in this and the following chapter are those of convection, diffusion, and gradi- ent/curvature drifts in the system’s magnetic field and its irregularities, while the momentum transport processes include adiabatic energy losses in the expanding system, synchrotron ra- diation, and inverse Compton scattering. As the transport equation is solved numerically, the model does not require any limitation on the convection velocity. The models presented in this chapter were, to a large degree, specifically developed for this study, and are based on finite difference techniques. An introduction to these numerical techniques, as well the derivation of the numerical schemes used to solve the transport equation, can be found in Appendix C.

In this chapter the transport equation is first solved in a one-dimensional, spherically-symmetric

system. As a result of having only one spatial dimension, the geometry of the nebular magnetic

field cannot be included in such a system. This is left for the next chapter, where a second spa- tial dimension, that of latitude, is added in order to study the effects of gradient and curvature drifts in the magnetic field. The results of this chapter initially focus on steady-state solutions that are built up in several steps to individually demonstrate the effect that the various trans- port processes have on the evolution of the spectrum. The first of these sections describes a system that only contains convection and adiabatic energy losses, as well as a system where the processes of synchrotron and inverse Compton energy losses are included. These results are primarily used to test the model, as the effects of the aforementioned energy losses are well known, and described in many textbooks on high-energy astrophysics (see, e.g., Longair, 2011).

The next section investigates the evolution of the particle spectrum in a system where convec- tion, diffusion, and adiabatic losses are important. Although the results in this section are less well-known, the same results have also been previously found by, e.g., Skilling (1971), Jokipii and Higdon (1979), Lerche and Schlickeiser (1981), and Atoyan et al. (1995), and are employed as an additional test of the model. The last of the time-independent sections form the main con- tribution of this chapter. Here the transport equation is solved with all of the above-mentioned processes included. It is shown that diffusion not only significantly moderates the high-energy cut-off resulting from synchotron losses, but can also be used to explain the spatial evolution of the X-ray synchrotron spectra observed from pulsar wind nebulae.

Subsequently, the focus shifts to time-dependent solutions that are again divided into three subsections. In the first of these, the evolution of an initially empty system is investigated, with this scenario again used as a test of the numerical scheme. In the next section, the effect of a time-dependent magnetic field on the evolution of the system is investigated (see Section 3.3.3).

As diffusion is inversely proportional and synchrotron losses proportional to the magnetic field, it is important to study any temporal variations that the magnetic field may be subjected to. The inclusion of a temporal component in the model was specifically motivated by the fact that the expanding outer boundary of the nebula may have an influence on the evolution of the particle spectra. This is investigated in the last of the sections dealing with time-dependent solutions.

The research presented in this chapter has been published in Vorster and Moraal (2013).

5.1 The transport equation

The general transport equation, expressed in terms of the directionally averaged distribution function f (r, p, t), is given by (see, e.g., Ginzburg and Syrovatskii, 1964; Parker, 1965)

∂f

∂t + ∇ · S − 1 p ²

∂

∂p (

p ² [

⟨ ˙p⟩f + D p

∂f

∂p ])

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

where p represents momentum, r the spatial coordinates, and t time. The differential current density

S = 4πp ² (CVf − K · ∇f) (5.2)

describes the convective and diffusive flux of particles in space, with V(r, t) denoting the con- vection velocity and K(r, p, t) the diffusion tensor. The Compton-Getting coefficient (see, e.g., Gleeson and Axford, 1968; Gleeson and Urch, 1973)

C = − 1 3

∂ ln f

∂ ln p (5.3)

takes into account the fact that an observer moving relative to the rest frame in which f is ho- mogeneous and isotropic will measure a modified current density. This modification produces a spectral effect similar in nature to the Doppler effect on photons. For a power-law spectrum f ∝ p ^−α , this coefficient has the value C = α/3.

The third term in (5.1) describes the transport of particles in momentum space as a result of continuous energy losses. In a spherically expanding medium, energy or momentum losses include the adiabatic loss rate

⟨ ˙p⟩ ^ad p = − 1

3 ∇ · V. (5.4)

However, Gleeson and Webb (1978) showed that when the particles reside in a wind with velo- city V, this is the momentum loss rate in the frame moving with the wind, while the loss rate in the stationary frame of reference is

⟨ ˙p⟩ ad

p = 1

3 V · ∇f

f . (5.5)

The modification of the convective streaming from Vf to CVf , and the adiabatic cooling rate from (5.4) to (5.5), produces effects in the transport equation that cancel each other. Gleeson and Webb (1978) showed that when C is dropped in S, the momentum loss rate in (5.1) transforms from (5.5) back to its more familiar form given by (5.4). It should be noted at this point that the process of first-order Fermi acceleration in shocks is formally included in the transport equation if it is treated as a discontinuous jump in the flow velocity.

For the non-thermal processes of synchrotron radiation and inverse Compton scattering in the Thomson regime, the momentum loss rate is given by

⟨ ˙p⟩ ^n−t = z p p ² (

1 + U IC

U _B )

, (5.6)

where

z _p = 4

3 σ _T 1 (m e c) ²

B ²

8π . (5.7)

It will be recalled from Section 4.2 that U IC and U B represent the energy density of the photon and magnetic fields, respectively. The transformation of the more familiar synchrotron and in- verse Compton energy loss rate (4.8) to the momentum loss rate (5.6) is discussed in Appendix B.1. It will be further recalled from Section 4.2 that the inverse Compton loss rate given in (5.6) is valid when U IC /U _B ≲ 3. If this condition is violated, the inverse Compton loss rate has to be modified, as discussed in Appendix B.2.

From (5.6) it can be seen that synchrotron radiation and inverse Compton scattering have the

same momentum dependence, implying that both processes effect the same evolution of the

particle spectrum. It is therefore possible to take inverse Compton scattering into account by suitably redefining the spatial and temporal variation of B(r, t). If the target photons are provided by the CMBR, the modification amounts to simply adding a constant value to B(r, t).

Explicit reference to inverse Compton losses will consequently be omitted from all subsequent discussions.

The fourth term in (5.1) describes stochastic or second-order Fermi acceleration in magneto- hydrodynamic turbulence, with D p acting as the momentum diffusion coefficient. Lastly, any sources or sinks are included in the transport equation through Q(r, p, t).

Inserting the expressions for S and ⟨ ˙p⟩ total = ⟨ ˙p⟩ ad + ⟨ ˙p⟩ n−t into (5.1), imposing spherical sym- metry (∂f /∂ϕ, ∂f /∂θ = 0) on the system, and neglecting momentum diffusion, D p ∂f /∂p = 0, leads to

∂f

∂t = κ ∂ ² f

∂r ² + [ 1 r ²

∂

∂r (r ² κ ) − V ] ∂f

∂r + [ 1

3r ²

∂

∂r (r ² V ) + z p p ] ∂f

∂ ln p + 4z _p pf + Q, (5.8) expressed in the spherical coordinate r. In this case K reduces to a single effective isotropic diffusion coefficient κ(p, r, t). For a more complete derivation of the various terms present in (5.8), Appendix B can be consulted.

While it is often more convenient to express transport equations in terms of the omni-directional distribution function f (r, p, t), it may be argued that the particle density N (r, p, t) is a more natural quantity to work with. This is due to the fact that the expressions for synchrotron and inverse Compton emission spectra are usually stated in terms of N (r, p, t) (see, e.g., Blu- menthal and Gould, 1970; Longair, 2011). The relation between these two quantities is given by N (r, p, t) = 4πp ² f (r, p, t), and the power-law spectrum f ∝ p ^−(α+2) therefore translates to N ∝ p ^−α .

5.1.1 Scaling of the transport equation

To keep the model as general as possible, it is useful to introduce a dimensionless scaling for the variables in (5.8). For the radial distance, the scaling

r = r ˜

˜

r _S (5.9)

is used, where ˜ r is the unscaled variable and ˜ r _S a characteristic length, chosen to be the size of the system. This leads to a system with a scaled radial dimension of 0 < r ≤ 1. A similar scaling is used for momentum

p = p ˜

˜

p _S , (5.10)

with the difference that ˜ p _S is not chosen to be the maximum particle momentum, but rather

an intermediate value. For the numerical simulations, the scaled momentum range is chosen

as 10 ⁻³ ≤ p ≤ 10 ⁴ . It is also possible to scale the particle energy in a similar fashion, such

that E = ˜ E/ ˜ E _S . As ˜ E ≈ ˜ pc for relativistic particles, it then follows that E = p. For the

numerical calculations the scale energy is chosen as ˜ E S = 1 TeV, with the scaled momentum range translating to the energy range 1 GeV ≤ ˜ E ≤ 10 PeV. This energy range includes the electron energy range needed to produce the non-thermal emission observed from PWNe.

The convection velocity in the system is expressed as a fraction of the speed of light V = V ˜

c , (5.11)

while the diffusion coefficient is scaled as

κ = ˜ κ

˜

r _S c . (5.12)

It is useful to define a scaling factor for time, chosen as the light transit time through the system t = c˜ t

˜ r S

. (5.13)

Lastly, the non-thermal energy loss coefficient (5.7) is scaled as z _p = z ˜ p r ˜ _S p ˜ _S

c . (5.14)

5.1.2 Useful time scales

The energy lost by a particle is directly related to its residence time in the system. In the absence of diffusion, the convection time scale is given by

τ con =

∫ r r

0

dr

V (r) , (5.15)

while the diffusion time scale was given in Section 4.2 as τ _dif = r ²

6κ . (4.17)

It should be remembered that diffusion is a stochastic process, and that τ dif represents an aver- age propagation time. Furthermore, this diffusion time scale was derived by Parker (1965) un- der the assumption that κ is independent of r, and that the particles undergo negligible energy losses during their propagation through the system. As both these assumptions are generally unrealistic for PWNe, it should be kept in mind that (4.17) is at best an approximation.

The synchrotron cooling time is estimated to be (see, e.g., De Jager and Djannati-Ata¨ı, 2009)

˜

τ syn ∼ 8.4 × 10 ³

B ² _µG E _TeV kyr. (5.16)

5.2 One-dimensional, steady-state solutions

5.2.1 Parameter values in the steady-state model

In the hydrodynamic modelling of Chapter 3, the evolution of the fluid quantities in the region

between the pulsar’s light cylinder and the termination shock r 0 formed part of the simula-

tions. However, in this chapter and the next, the calculation of the evolution of the particle

spectrum will be confined to the PWN, defined as the region between r 0 and the outer bound- ary of the PWN. As such, it would not make sense to place the inner boundary of the present simulations upstream of r 0 . As discussed in Sections 2.3.1 and 2.3.3, not only is the accelera- tion of particles at the relativistic MHD shock of the PWN still far from well understood, but synchrotron observations of the inner regions of a number of PWNe (see, e.g. G21.5-0.9: Boc- chino et al., 2005; Camilo et al., 2006; Vela X: Mangano et al., 2005; MSH 15-52: Sch¨ock et al., 2010) provide a direct measurement of the particle spectrum injected into the PWN at r 0 .

The above-mentioned synchrotron observations indicate that a two-component power-law spectrum is required to account for the radio and X-ray synchrotron emission from PWNe. For simplicity, the spectrum injected at the inner boundary of the present simulations is chosen to be the single power law f (p, r 0 ) ∝ p ⁻⁴ , or equivalently, N (p, r 0 ) ∝ p ⁻² . This is comparable to the values of N (p, r 0 ) ∝ p ^−1.8±0.2 − N(p, r 0 ) ∝ p ^−2.5±0.1 observed for G21.5-0.9 (Safi-Harb et al., 2001; Bocchino et al., 2005), N (p, r 0 ) ∝ p ^{−2.00±0.04} observed for Vela X (Mangano et al., 2005), N (p, r ₀ ) ∝ p ^−1.6±0.2 observed for 3C 58 (Slane et al., 2004), and N (p, r 0 ) ∝ p ^{−2.32±0.04} observed for MSH 15-52 (Sch¨ock et al., 2010). The specific choice for the power-law index used in the model is not only motivated by X-ray synchrotron observations taken in the vicinity of r 0 , but it is also the hardest spectrum that one would expect from diffusive shock acceleration.

Two-dimensional, axisymmetric simulations show a complicated flow and magnetic field struc- ture in PWNe (see, e.g., Komissarov and Lyubarsky, 2003; Del Zanna et al., 2004). However, to ob- tain the one-dimensional solutions presented in this section, a radial flow and azimuthal mag- netic field is assumed, identical to the conditions implemented in the hydrodynamic model of Chapter 3. The assumption of an azimuthal magnetic field is also based on the derivation presented in Section 2.2.2. It is further assumed that the steady-state, ideal MHD limit (3.2) holds, allowing one to relate the radial profiles of the velocity and magnetic field. With the chosen flow and magnetic field structure, (3.2) reduces to

V Br = constant = V 0 B ₀ r ₀ . (5.17)

Figure 3.4 shows that the convective flow downstream of the pulsar wind termination shock is subsonic, with the post-shock velocity in the immediate vicinity of the termination shock being of the order ˜ V ₀ = 0.3c. This value is therefore adopted in the present simulations. Additionally, the figure also shows that the magnetic field increases with distance, so that ˜ B ∝ r, while the velocity decreases as ˜ V ∝ 1/r ² . As discussed in Section 3.3.5, these are the expected profiles when the ratio of electromagnetic to particle energy is σ < 0.01. On the other hand, it was also discussed in the same section that the velocity remains almost independent of r when σ = 1, leading to the dependence ˜ B ∝ 1/r.

To illustrate the effect of the radial dependence on the solutions, three scenarios are investig- ated: V ∝ 1/r, V ∝ 1/r ^0.5 , and V = V 0 , with these profiles chosen for a number of reasons.

The first is to emphasise that the present transport model has applications beyond the study of

PWNe. The second reason is motivated by previous modelling of PWNe. From the literature

only four PWNe could be found where the radial profile of the velocity has been inferred. This includes the Crab nebula (Kennel and Coroniti, 1984a,b), the inner region of the Vela PWN (Vor- ster, 2010), and MSH 15-52 (Sch¨ock et al., 2010). For the Crab nebula a profile similar to the one in Figure 3.4 was derived, while a decreasing velocity profile was derived for the other two sources. It is important to note that these three models neglected diffusion. With the effect of diffusion included, Van Etten and Romani (2011) derived a velocity profile ˜ V ∝ 1/r ^0.5 for HESS J1825–137, motivating the second profile that will be used in the simulations. The final reason follows from the results themselves, from which it will be seen that different radial profiles do not lead to new spectral features. This is addressed in more detail at the end of Section 5.2.5.

The chosen velocity profiles, together with (5.17), respectively lead to the magnetic profiles, B = B ₀ , B ∝ 1/r ^0.5 , and B ∝ 1/r. For HESS J1825–137, Van Etten and Romani (2011) derived a magnetic field profile B ∝ 1/r ^0.7 , similar to the second scenario. For the magnetic field it is difficult to specify a characteristic value at the termination shock as very few estimates of B ˜ 0 have been made. Additionally, ˜ B 0 is also expected to vary significantly between different nebulae. One need only consider the variation of the average magnetic fields ¯ B in various nebulae, as ¯ B depends on ˜ B ₀ . For example, the Crab Nebula has an average magnetic field strength of ¯ B ∼ 300 µG (see, e.g., Trimble, 1982), whereas it is estimated that the Vela PWN has a lower value of ¯ B ∼ 3 µG (see, e.g., De Jager and Djannati-Ata¨ı, 2009). One estimate of ˜ B ₀ that can be found in the literature is that of Van Etten and Romani (2011), where these authors estimated ˜ B ₀ ≈ 400 µG for HESS J1825–137. For the simulations, a similar value of ˜ B ₀ = 350 µG is chosen as the magnetic field strength at the inner boundary.

Following, e.g., Lerche and Schlickeiser (1981), it is assumed that the radial (κ r ) and momentum (κ p ) dependence of the diffusion coefficient are separable,

κ(r, p) = κ ₀ κ _r κ _p , (5.18)

with κ 0 a constant. Diffusion occurs as a result of particle interaction with irregularities in the magnetic field, and experience with the propagation of cosmic rays in the heliosphere shows that the radial dependence of the diffusion coefficient can be modelled as κ r ∝ 1/B, while the momentum dependence scales roughly as κ p ∝ p (see, e.g., Caballero-Lopez and Moraal, 2004, and references therein). With the chosen magnetic field profiles discussed in the previous paragraph, this leads to the radial dependencies κ r = κ ₀ , κ r ∝ r ^0.5 , and κ r ∝ r in the model, while the momentum dependence is chosen to be κ p = ˜ p/˜ p _S .

To account for the escape of particles with ˜ E > 100 GeV from the gamma-ray emission region

in the Vela PWN, a spatially independent diffusion coefficient of ˜ κ ∼ 10 ²⁷ E _TeV cm ² s ⁻¹ has

been estimated by Hinton et al. (2011). A similar ˜ κ has also been derived by Van Etten and

Romani (2011) to explain TeV gamma-ray emission from the source HESS J1825–137. If a PWN

has a size of ˜ r _S = 10 pc and ˜ E _S = ˜ p _S c = 1 TeV is chosen as the scaling for the energy, then

it follows from (5.12) and (5.18) that the normalisation κ 0 = 10 ⁻³ will lead to a value for κ

that is similar to the values derived from observations. For the purposes of illustration, the

smaller value κ 0 = 4 × 10 ⁻⁵ is chosen at the inner boundary, and is motivated by the scaling κ ∝ 1/B. As the average magnetic field in the model is larger than ¯ B in the Vela PWN, a smaller diffusion coefficient is appropriate.

5.2.2 Numerical considerations of the steady-state scheme

The steady-state (∂f /∂t = 0) version of (5.8) is solved using the second-order accurate Crank- Nicolson solution method, with the numerical scheme derived in Appendix C.6.1. In this steady-state version, momentum plays the role usually assigned to time, and is thus used as the stepping parameter of the parabolic equation. Because particles migrate downwards in energy or momentum space as a result of energy losses, the numerical scheme steps logarith- mically from high to low momentum values. For the solutions presented in this section, the radial grid ranges between r 0 ≤ r ≤ 1, with ∆r = 0.001, while the stepping parameter ranges between 10 ⁻³ ≤ p ≤ 10 ⁴ , with step size ∆ ln p = −0.02.

Particles propagating away from the source immediately lose energy, and as a result particles with momentum p max can only exist at the source located at r 0 = 0.01, leading to the ”initial”

condition f (p max , r) = 0 for all r > r 0 . As the source of particles Q can be described by a δ-function at the inner boundary of the system, it is possible to include the injected particles in the problem by specifying the inner boundary condition described below.

The number of particles per momentum interval that flow through the inner radial boundary must be equal to the total number of particles produced per time and momentum interval, Q ₀ = Q ^∗ (p/p ^∗ ) ^−(α+2) , with Q ^∗ and p ^∗ normalisation constants. This implies that

I

S ₀ · dA 0 = Q ₀ , (5.19)

where dA 0 = r ² ₀ sin θdθdϕe r is the surface element on the inner boundary, and S 0 is given by (5.2). Integrating (5.19) leads to the inner boundary condition

CV ₀ f − κ 0

∂f

∂r = Q ^∗ 1 4πr ₀ ²

( p 0

p ^∗

) −(α+2)

(5.20) that is solved simultaneously with (5.8). A similar boundary condition was derived by Ng and Gleeson (1975) to describe the evolution of the energy spectrum of cosmic rays produced by solar flares. In the absence of diffusion, (5.20) reduces to

f (r 0 , p 0 ) = Q ^∗ CV ₀

1 4πr ² ₀

( p 0

p ^∗

) −(α+2)

. (5.21)

To simulate particles escaping from the system, the free-escape condition f (r n , p) = 0 is im-

posed at the outer boundary. The effect that the choice of outer boundary condition has on

the solutions will be discussed in more detail in Section 5.2.4, where additional solutions are

presented for a no-escape boundary, defined in the model as S = 0.

5.2.3 Convection and convection-synchrotron

This section presents steady-state solutions to (5.8) for the convection (conv) and convection- synchrotron (conv-sync) scenarios. The former includes only the processes of convection and adiabatic losses, while the latter scenario additionally includes the process of synchrotron ra- diation. The solutions for the two systems are divided into three subsets corresponding to the three different radial velocity/magnetic field profiles used, with these solutions shown in Figures 5.1(a)–(c). Shown in Figures 5.1(d)–(f) are the radial profiles for the particle intensities corresponding to the convection scenario.

Figures 5.1(a)–(c) show that in the absence of diffusion and synchrotron losses, the intensity of the spectra decreases with increasing radial distance, while retaining the spectral shape of the source spectrum. As this decrease is a result of adiabatic losses, it follows from the energy loss rate (5.4) that the particles suffer larger losses when the velocity profile changes from V ∝ 1/r to a constant radial profile V = V 0 .

Turning to the radial profiles, Figures 5.1(d)–(f) show that the particle intensity can be de- scribed by simple power laws: (d) N ∝ r ^−4/3 ; (e) N ∝ r ⁻² ; and (f) N ∝ r ^−8/3 . It can also be seen that the numerical scheme becomes unstable in the outer parts of the system, leading to the intensities of the numerical solutions oscillating about these power laws. To illustrate these oscillations more clearly, the ˜ E = 10 ⁻³ TeV radial profiles shown in Figures 5.1(d)–(f) have been enlarged for the region extending between 0.88 ≤ r ≤ 0.9, and are shown in Figures 5.2(a)–(c).

In the absence of diffusion, the character of (5.8) changes from parabolic to hyperbolic, and the Crank-Nicolson method may no longer be the most suitable scheme for solving this equation.

As the instabilities in Figures 5.1(d)–(f) are independent of energy, the correct shape of the spectrum is nevertheless obtained from the numerical scheme.

Apart from such possible stability issues, it is also important to investigate the accuracy of the numerical scheme. The simplest way to do this is to compare the results from the numerical scheme with analytical solutions. For the convection scenario, the transport equation (5.8) is amenable to a simple analytical solution when the velocity profile is described by

V (r) = V ₀ ( r ₀ r

) β

. (5.22)

With this prescription, (5.8) reduces to

∂f

∂r − β − 2 3r

∂f

∂ ln p = 0. (5.23)

Using the method of characteristics, (5.23) can be reduced to three first-order, ordinary differ- ential equations, with the first of these equations given by

dr

dp = − 3r

(β − 2)p . (5.24)

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) V ∝ 1/r, B = const.

(a)

r=0.1 r=0.9 source

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) V ∝ 1/r, B = const.

(a)

r=0.1 r=0.9 source

conv

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) V ∝ 1/r, B = const.

(a)

r=0.1 r=0.9 source

conv

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) V ∝ 1/r, B = const.

(a)

r=0.1 r=0.9 source

conv conv-sync

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) V ∝ 1/r, B = const.

(a)

r=0.1 r=0.9 source

conv conv-sync

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5

(b) r=0.1

r=0.9 source

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5

(b) r=0.1

r=0.9 source

conv

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5

(b) r=0.1

r=0.9 source

conv

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5

(b) r=0.1

r=0.9 source

conv conv-sync

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5

(b) r=0.1

r=0.9 source

conv conv-sync

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., B ∝ 1/r

(c) r=0.1

r=0.9 source

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., B ∝ 1/r

(c) r=0.1

r=0.9 source

conv

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., B ∝ 1/r

(c) r=0.1

r=0.9 source

conv

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., B ∝ 1/r

(c) r=0.1

r=0.9 source

conv conv-sync

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., B ∝ 1/r

(c) r=0.1

r=0.9 source

conv conv-sync

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r (d)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 (e)

E ~ 1 E ~

2 E ~

3 E ~

4

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const.

(f)

E ~ 1 E ~

2 E ~

3 E ~

4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

Figure 5.1: Numerical solutions to the one-dimensional, steady-state version of the transport equa-

tion (5.8), using various radial profiles for the velocity and magnetic field. Figures (a)–(c) show the

spectra calculated in the inner and outer parts of the system for the convection (conv) and convection-

synchrotron (conv-sync) scenarios. The term ”convection” is used as short-hand for the phrase ”convec-

tion with adiabatic losses”. Note that the differential number density N = 4πp ² f is plotted, instead of

the distribution function f . To highlight characteristics of the spectrum, N is multiplied by ˜ E ² . Figures

(d)–(f) show the corresponding radial profiles in particle intensity for the convection scenario at various

energy values. The dotted lines represent the analytical solution (5.26).

Figure 5.2: Numerical instabilities that arise when diffusion is neglected. The radial profiles shown in this figure are expansions of the ˜ E = 10 ⁻³ TeV profiles shown in Figures 5.1(d)–(f).

Integrating (5.24) gives

p ₀ = p ( r ₀ r

) (β−2)/3

, (5.25)

where p 0 is the initial momentum (energy) of the particle at the inner boundary of the system, located at r 0 . Here p 0 > p, due to the fact that the particles are subjected to energy losses. As f does not appear explicitly in (5.23), the solution to the partial differential equation is constant f (r ₀ , p ₀ ) = f (r, p) along the curves described by (5.25). Inserting (5.25) into the boundary condition (5.20) leads to the analytical solution

f (r, p) = Q ^∗ 4πr ₀ ² CV 0

[ p ( r ₀

r

) (β−2)/3 ] −(α+2)

H (r − r 0 ) , (5.26) where H is the Heaviside function.

The analytical solution (5.26) shows that the spectra at any position in the system have the same p-dependence as the source spectrum, with the only difference being that the intensities are reduced. This behaviour can be ascribed to the fact that the fractional momentum change as a result of adiabatic losses (5.4) is independent of momentum. The same result has also been previously found by, e.g., Jokipii and Higdon (1979) and Lerche and Schlickeiser (1981). The analy- tical solution (5.26) further shows the well-known result that particles in a spherical geometry will not be subjected to adiabatic energy changes when α = −2, i.e., when the convection velocity is described by the profile V (r) ∝ 1/r ² .

Figures 5.1(d)–(f) also show the analytical solutions calculated using (5.26), where it can be

seen that the numerical and analytical solutions not only predict the same radial dependence of the particle intensity, but it also follows that the particle intensity predicted by the numerical scheme is within 10% of the analytical solutions. The numerical scheme can thus be used with confidence.

The next step is to extend the model by including the effect of synchrotron losses, along with the processes of convection and adiabatic losses. The inclusion of this process leads to the de- velopment of a sharp cut-off in the spectra, as shown in Figures 5.1(a)–(c). The behaviour of the spectra in Figures 5.1(a)–(c) is a well-known effect of synchrotron losses, and is discussed in many textbooks on astrophysics such as Longair (2011). The synchrotron losses suffered by a particle are directly related to its confinement time in a system, and in the case of a convection- synchrotron system, (5.15) shows that the confinement time τ ad is independent of E. On the other hand, the synchrotron time scale τ syn (E) given by (5.16) is an estimate of the time re- quired by a particle to lose the largest fraction of its energy as a result of synchrotron radiation.

Therefore, the particles in a convection-synchrotron system that have τ syn ≲ τ ad will have lost all their energy, leading to the cut-off visible in Figures 5.1(a)–(c). As synchrotron losses are proportional to E ² , as shown in (5.6), in contrast to the E dependence of adiabatic losses de- scribed by (5.4), it is the higher-energy particles that are affected by synchrotron losses, while the low-energy particles are only subjected to adiabatic losses.

Figures 5.1(a)–(c) also show that the cut-off shifts to lower energies with increasing r. Com- pared to r = 0.1, the confinement time increases as one moves to r = 0.9. Coupled with the fact that τ syn ∝ 1/E, one consequently has τ syn ≲ τ _ad for smaller energy values at r = 0.9, leading to the aforementioned shift. A comparison of Figures 5.1(a)–(c) shows that, for a given value of r, the cut-off energy increases as the velocity profile changes from V ∝ 1/r to V = constant. Not only does τ syn increase from Figure 5.1(a) to Figure 5.1(c), but the synchrotron losses suffered by the particles also decrease. This last statement follows from the fact that the synchrotron loss rate is proportional to B ² , as shown in (5.6).

Additionally, a small peak is visible just before the cut-off in the spectra at r = 0.9. The most in- tuitive explanation for this peak would be that synchrotron losses lead to a pile-up of particles at the edge of the cut-off. Such an effect is indeed expected for the source spectrum N ∝ p ^−α when α < 2 (see, e.g., Longair, 2011). However, the peak in Figure 5.1 (where the source spec- trum has α = 2) is not a physical effect, but directly related to the instabilities visible in Figures 5.1(d)–(f). The spectra at r = 0.9 shown in Figures 5.1(a)–(c) all correspond to the case where the numerical scheme over-predicts the intensities. It follows from Figures 5.2(a),(b) that if one were to plot the spectra at the neighbouring grid points, then one would find an under- prediction of the intensities.

Following the discussion presented thus far, the energy ˜ E _syn where the synchrotron break

should appear in Figures 5.1(a)–(c) can be estimated by setting τ syn = τ _ad . The value for the

magnetic field used in the estimate is the average value in the region ranging from r 0 to r. The

estimates for the three scenarios, together with the values predicted by the numerical scheme,

Table 5.1: The synchrotron break energy estimated using (5.15) and (5.16), along with the values predicted by the numerical scheme.

V ∝ 1/r, B = const. V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5 V = const., B ∝ 1/r Estimate Figure 5.1(a) Estimate Figure 5.1(b) Estimate Figure 5.1(c)

r = 0.1 1.2 ∼ 0.7 12.6 ∼ 6 102 ∼ 20

r = 0.9 0.015 ∼ 0.01 2.9 ∼ 1.4 267 ∼ 4

are summarised in Table 5.1.

A comparison of the results derived from Figures 5.1(a),(b) with the order-of-magnitude es- timates shows a factor ≲ 2 difference in the break energies. Note that this difference does not reflect negatively on the numerical simulations. Rather, this highlights the necessity of using more rigorous calculations, such as solving the Fokker-Planck transport equation, in- stead of only relying on the estimates that are often used. However, these estimates are useful in determining whether the numerical scheme predicts the correct values within an order-of- magnitude.

For the scenario shown in Figure 5.1(c), the estimate predicts that the break energy at r = 0.9 should be larger than the value at r = 0.1. Using the average magnetic field to calculate the estimates in Table 5.1 is a better approximation in the inner part of the system where the magnetic field is large, and the particles consequently lose the most energy. However, in the last scenario B ∝ 1/r, leading to a magnetic field that differs significantly between the inner and outer boundary of the PWN. As τ syn ∝ ¯ B ² , the average magnetic field used for the region r ₀ ≤ r ≤ 0.9 will lead to a significant under-prediction of the synchrotron losses, leading to the discrepancy in the predicted value of ˜ E _syn .

5.2.4 Convection-diffusion

To determine the effect of diffusion on the evolution of the particle spectrum, this section con-

siders steady-state solutions to (5.8) with the processes of convection, diffusion and adiabatic

losses included, where it will be recalled from Section 5.2.1 that the radial dependence of the

diffusion coefficient scales as κ r ∝ 1/B(r). A comparison of the spectra in the convection-

diffusion system, Figures 5.3(a)–(c), with those of the convection system, Figures 5.1(a)–(c),

shows that these spectra can be divided into two regimes. At lower energies convection

and adiabatic losses dominate and the spectra retain the shape of the source spectrum. At

higher energies diffusion is the more important transport process, leading to an evolution in

the spectral index. With the free-escape (Dirichlet) condition imposed at the outer boundary,

the diffusion-dominated parts of the spectra become softer by one power-law index so that

N ∝ p ⁻³ , or expressed in terms of energy, N ∝ ˜ E ⁻³ . This last form follows from the fact that

the particles are relativistic, i.e. E = pc.

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) V ∝ 1/r, κ = const.

(a) r=0.1 r=0.9

source

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) V ∝ 1/r, κ = const.

(a) r=0.1 r=0.9

source free esc

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) V ∝ 1/r, κ = const.

(a) r=0.1 r=0.9

source free esc

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) V ∝ 1/r, κ = const.

(a) r=0.1 r=0.9

source free esc no esc

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) V ∝ 1/r, κ = const.

(a) r=0.1 r=0.9

source free esc no esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , κ ∝ r ^0.5

(b) r=0.1 r=0.9

source

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , κ ∝ r ^0.5

(b) r=0.1 r=0.9

source free esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , κ ∝ r ^0.5

(b) r=0.1 r=0.9

source free esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , κ ∝ r ^0.5

(b) r=0.1 r=0.9

source free esc no esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V ∝ 1/r ^0.5 , κ ∝ r ^0.5

(b) r=0.1 r=0.9

source free esc no esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., κ ∝ r

(c) r=0.1

r=0.9 source

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., κ ∝ r

(c) r=0.1

r=0.9 source free esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., κ ∝ r

(c) r=0.1

r=0.9 source free esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., κ ∝ r

(c) r=0.1

r=0.9 source free esc no esc

10 ^-3 10 ^-2 10 ^-1 10 ⁰ 10 ¹ 10 ² 10 ³ E ~

=p ~c (TeV) V = const., κ ∝ r

(c) r=0.1

r=0.9 source free esc no esc

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

10 ^-6 10 ^-5 10 ^-4 10 ^-3 10 ^-2 10 ^-1

0.01 0.1

r 2 E ~ 7/4 N (norm. units)

r (norm. units) V ∝ 1/r, κ = const.

(d)

E ~ 1 E ~

2 E ~

3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V ∝ 1/r ^0.5 , κ ∝ r ^0.5 (e)

E ~ 1 E ~

2

E ~ 3

E ~ 4

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

0.01 0.1

r (norm. units) V = const., κ ∝ r (f)

E ~ 1 E ~

2

E ~ 3

E ~ 4

E ~

1 =10 ^-3 E ~

2 =10 ^-1 E ~

3 =10 ¹ E ~

4 =10 ³

Figure 5.3: The same as Figure 5.1, but for a convection-diffusion system. Also shown are the predicted

spectra for a system with a no-escape outer boundary. The radial profiles of Figure (d)–(f) are for a

convection-diffusion system with a free-escape outer boundary, while the dotted lines represent the

analytical solution calculated for the convection scenarios of Figure 5.1.

Table 5.2: The energy values where convection and diffusion are equally important, as estimated from (5.28) and the numerical calculations.

V ∝ 1/r, B = const. V ∝ 1/r ^0.5 , B ∝ 1/r ^0.5 V = const., B ∝ 1/r

κ r = const. κ ∝ r ^0.5 κ ∝ r

Estimate Figure 5.3(a) Figure 5.3(b) Figure 5.3(c)

r = 0.1 75 ∼ 90 ∼ 100 ∼ 110

r = 0.9 75 ∼ 20 ∼ 40 ∼ 90

This softening of the spectra is typical of escape losses, and is related to the momentum de- pendence of the diffusion coefficient. For the general case κ ∝ p ^λ , the diffusion-dominated parts of the spectra are described by N ∝ p ^−α−λ (see, e.g., Jokipii and Higdon, 1979; Lerche and Schlickeiser, 1981; Atoyan et al., 1995), with α = 2 and λ = 1 in the present model.

The importance of diffusion relative to convection can be estimated using the P´eclet number (see, e.g., Prandtl, 1953)

ξ = V r

κ . (5.27)

Note that ξ is essentially the ratio of the convection and diffusion time scales, (5.15) and (4.17).

When ξ ≫ 1 the system is convection dominated, while being diffusion dominated for ξ ≪ 1.

With the ideal MHD limit imposed, the scaling V = V 0 (r ₀ /r) ^β and κ r ∝ 1/B (see Section 5.2.1) implies that κ = κ 0 (r/r ₀ ) ^(1−β) κ _p . Inserting these radial profiles into (5.27) leads to

ξ = V ₀ r ₀ κ 0 κ p

, (5.28)

where it will be recalled from Section 5.2.1 that κ p = ˜ p/˜ p _S . It thus follows that the value of ξ is not uniquely defined for the system. For the spectra in Figures 5.3(a)–(c), one has ξ > 1 at low energies, while ξ < 1 at high energies. The energy ˜ E _dif = ˜ p _dif c where the two transport processes are equally important can be found by setting ξ = 1. As (5.28) only depends on the values at the inner boundary, ξ is not only independent of radial position, but also has the same value for the three different scenarios shown in Figures 5.3(a)–(c). Using the values r 0 = 0.01, V 0 = 0.3 and κ 0 = 4×10 ⁻⁵ chosen for the numerical calculations leads to the estimate E ˜ _dif ∼ 75 TeV. The values of ˜ E _dif derived from Figure 5.3 are summarised in Table 5.2.

Table 5.2 shows that the estimates derived from (5.28) and Figure 5.3 are in reasonable agree- ment. It is again emphasised that (5.28) is at best an order-of-magnitude estimate. Table 5.2 further shows that the numerical simulations predict that the value of ˜ E _dif decreases as a func- tion of radial position. Because κ increases as a function of radial distance, while V decreases, diffusion will become more important than convection in the outer parts of the system, leading to a decrease in the value of ˜ E _dif .

Figures 5.3(d)–(f) show the radial profiles for the three scenarios in the convection-diffusion

simulations. These radial profiles are similar to those of the convection system, Figures 5.1(d)–

(f), for energies smaller than ˜ E _dif . Again, this is a reflection of the fact that convection dom- inates at lower energies. At ˜ E = 10 ³ TeV, where diffusion is important, the gradients of the radial profiles are smaller than for the convection counterparts, while also having a reduced intensity. These reduced intensities are a consequence of particles escaping from the system, f (r = 1, p) = 0, while the smaller gradients can be attributed to the fact that diffusion will reduce the gradients in particle intensity. However, this decrease in the gradient is partially cancelled as particles can escape freely at the outer boundary. The free-escape condition is further reflected in the cut-off that appears in the outer part of the system at r ≥ 0.7.

To investigate the role of the outer boundary condition, spectra were also calculated using a no-escape boundary condition, S(r = 1, p) = 0. The results are shown together with the free- escape spectra in Figures 5.3(a)–(c), where it can be seen that the spectra at the highest energies converge. This is a direct result of particles not being able to escape from the system. At the highest energies diffusion will effectively cancel any radial gradients in the particle density, leading to a radially-independent intensity. At r = 0.9, adiabatic losses have reduced the intensity in the low-energy part of the spectrum below that of the converging intensity at high energies, leading to the sharp transition between the diffusion and convection regimes.

The difference between the free-escape and the no-escape conditions represents the maximum effect that the outer boundary can have on the solutions, as a more realistic system will have a partial escape boundary. It should nevertheless be kept in mind that the results presented in Figures 5.3(a)–(c) are of an artificial nature, as a no-escape boundary will lead to an increase in particle intensities with time. A correct treatment of a no-escape boundary problem therefore requires that the transport equation be solved time-dependently.

5.2.5 Convection-diffusion-synchrotron

The final set of steady-state solutions focuses on the evolution of the spectra when synchro- tron losses are important in a diffusive system. This section deals with solutions to (5.8), with the transport processes of convection and diffusion included, together with adiabatic and syn- chrotron energy losses. The results for three scenarios are plotted in Figures 5.4(a)–(c), with the corresponding radial profiles for the particle intensity shown in Figures 5.4(d)–(f).

A comparison between Figures 5.1(a)–(c) and Figures 5.3(a)–(c) shows that the particles that

lose the most energy due to synchrotron radiation are generally the particles with the largest

diffusion coefficient. This follows from the fact that both diffusion and synchrotron losses

scale with energy. However, as discussed in Section 5.2.1, the model assumes that κ ∝ p, while

(5.6) shows that synchrotron losses scale as p ² . Therefore, at the highest energies synchrotron

losses will dominate diffusion. With the exception of the highest energies, it follows from the

escape time scale (4.17) that τ dif ∝ 1/κ, and diffusion thus decreases the propagation time

of particles through the system, compared to the propagation time through a system where

only convection is important. The synchrotron losses suffered by the particles are therefore