Radiative transfer inside accretion columns of neutron stars at extremely high mass accretion rates

Anton Pannekoek Institute for Astronomy

B

T

Radiative transfer inside accretion

columns of neutron stars at extremely high

mass accretion rates

Author: Boris WOLVERS Student ID: 10801936 Supervisor: Dr. Alexander MUSHTUKOV Examiner: Prof. dr. Michiel van der Klis Second reader: Dr. Anna Watts

Credits: 15 EC 03-04-2017 - 10-07-2017

Abstract

In this bachelor project we have studied the radiation transfer inside an accretion column above highly magnetized neutron stars. Accreting material is confined in to a narrow wall of magnetic funnel in the column. The problem of radiative transfer in an accretion column can be considered as a problem of radiative transfer in a layer. We have constructed a numerical model of radiative transfer of polarized X-ray radiation of a layer with fixed geometrical thickness. We assumed that the major process shaping the radiative transfer is Compton scattering, which is strongly modified in case of the high magnetic fields typical for X-ray pulsars (XRPs). The photon polarization state in highly magnetized plasma can be described by 2 normal modes, which significantly simplifies the problem. Using our code (written in C) we can get the distribution of photon mode intensities over the directions in the center and at the border of a layer. The distribution in the center is crucially important for future estimations of radiative pressure inside the column, which shapes the whole column structure. The distribution of polarized intensities at the border is directly related to the pulse profiles of XRPs and will be important for upcoming X-ray polarimeters. Furthermore, a better understanding of accretion columns in general may provide an explanation of recently discovered pulsating ultra-luminous (L∼ 1040erg/s) X-ray sources.

Nederlandse Samenvatting

Neutronensterren zijn enorm fascinerende objecten. Moet je eens voorstellen, een ster die twee keer zoveel weegt als onze zon maar dan ter grootte van Amsterdam! Hierdoor heeft een neutronenster een enorm hoge dichtheid. Een theelepel neutro-nenster zou maar liefst 5 miljard ton wegen! Doordat neutroneutro-nensterren zo compact zijn hebben ze ook een heel sterk zwaartekrachtsveld. Maar daar blijft het niet bij, ze hebben ook hele sterke magnetische velden. Deze zijn zodanig sterk, dat ze de fun-damentele eigenschappen van materie zelfs kunnen veranderen. Het is dan ook niet verwonderlijk dat in de laatste decennia het onderzoek naar deze objecten flink is toegenomen, aangezien de extreme fysica die deze sterren kenmerkt niet na te boot-sen is in laboratoria hier op aarde. Zorgvuldig onderzoek doen naar deze objecten kan er voor zorgen dat we fundamentele natuurkunde beter begrijpen.

In dit bachelor project doen we theoretisch onderzoek naar neutronensterren met een extreem hoog magnetisch veld, in de nabijheid van een normale ster (zie figuur

1).

FIGURE1: Een neutronenster met een extreem hoog magnetisch veld dat materie opslokt van een buurtster. Hierbij komt r¨ontgenstraling

vrij.

Doordat de neutronenster een sterke aantrekkingskracht heeft kan de neutronen-ster materie van de nabije neutronen-ster opslokken, dit proces wordt accretie genoemd. Bij dit proces zal er een zogenaamde accretieschijf ontstaan rond de neutronenster en va-nuit deze schijf kan er materie met hoge snelheid naar de magnetische polen van de neutronenster stromen. Wanneer de materie het oppervlak van de neutronenster nadert zal er r¨ontgenstraling vrijkomen die dan uiteindelijk op aarde gedecteerd kan worden en naderhand bestudeerd. De r¨ontgenstraling die uitgezonden wordt zorgt overigens ook voor een stralingsdruk met een gradi¨ent tegengesteld aan de richting van de invallende materie. In het geval er weinig materie op de neutronenster valt zal de stralingsdruk ook niet zo hoog zijn. Maar in het geval er veel materie op de ster valt dan zal de stralingsdruk erg hoog worden. Op een gegeven moment wordt de stralingsdruk dermate hoog dat de invallende materie het oppervlak van de ster niet eens zal bereiken. Vlak voordat dit gebeurt zal de straling dan ook de meeste lichtkracht opleveren. De maximale lichtkracht die uitgezonden kan worden wordt de Eddingtonlichtkracht genoemd.

Alleen er is iets geks aan de hand in de huidige astronomie. Uit observationele waarnemingen blijkt dat er tegenwoordig al drie neutronensterren zijn die een licht-kracht kunnen uitzenden hoger dan de Eddingtonlichtlicht-kracht! Dit is zo bizar dat er tot op heden nog geen goed model is om dat te kunnen verklaren. Echter, een mo-gelijk model dat het zou kunnen verklaren is het accretie kolom model. Dit model beschrijft hoe bij een hoge massaoverdracht vanuit de accretieschijf zogenaamde cilinders bovenop de magnetische polen kunnen verschijnen. Deze cilinders bestaan uit heet plasma waaruit mogelijk de hoge lichtkracht uitgezonden kan worden. Echter, dit model is ongelofelijk ingewikkeld en om dit model beter te kunnen begrijpen moeten we weten hoe de straling in zo’n kolom/cilinder zich gedraagt. Het voor-naamste proces van de straling in zo’n kolom is Compton verstrooiing. Echter, in een extreem sterk magnetisch veld zal de straling zich heel anders gedragen. Om deze straling zo nauwkeurig mogelijk te beschrijven worden er in dit project dif-ferentiaalvergelijkingen opgesteld, rekening houdend met het extreem sterk mag-netische veld en wat dit veld teweeg kan brengen op hele kleine schaal. Het bli-jkt namelijk dat fotonen met verschillende polarisaties zich door het magnetische veld verschillend gedragen. De opgestelde vergelijkingen worden daarna met be-hulp van zelfgeschreven code met de computer opgelost omdat deze vergelijkingen behoorlijk complex zijn. Niet alleen zijn ze complex maar het duurt ook erg lang om ze te berekenen, in sommige gevallen zelfs een paar weken. Desalniettemin hebben we resultaten kunnen verkrijgen waardoor we het accretie kolom model beter kunnen begrijpen en zijn we een stap dichterbij om de té hoge lichtkracht van de waargenomen neutronensterren te kunnen beschrijven!

Contents

Abstract 1

Nederlandse Samenvatting 2

Contents 4

1 Introduction 5

1.1 Physics in a strong magnetic field . . . 7

1.1.1 Charged particles and cyclotron energy . . . 7

1.1.2 Photon polarization in magnetized plasma . . . 8

1.1.3 Cross sections of linearly polarized modes . . . 8

1.2 Physical parameters inside the column. . . 10

1.2.1 Optical thickness of the accretion column . . . 10

1.2.2 Typical escape time of photons from the accretion column . . . 12

1.2.3 Length of advection . . . 13

1.2.4 Initial sources of emission . . . 13

2 Radiative transfer equations 15 2.1 One-dimensional case . . . 15

2.1.1 Numerical method . . . 16

2.1.2 Numerical results . . . 17

2.2 Three-dimensional: non-polarized . . . 19

2.3 Three-dimensional: polarized . . . 21

2.3.1 Scattering cross sections . . . 24

3 Results 28 3.1 Non-polarized . . . 28

3.1.1 Average number of scatterings . . . 30

3.2 Polarized . . . 31

3.2.1 Initial emissivity distributions of O-and X-mode photons . . . 33

3.2.2 Ratio of the photon energy density in the center . . . 34

3.2.3 Comparison of the photon energy densities . . . 35

4 Conclusions 37

1

Introduction

Neutron stars (NSs) are among the most dense objects in the universe. They can weigh as twice as much as the sun whilst their radius is ∼ 10 − 15 km. This leads to huge gravitational fields. On top of that they possess extremely high magnetic fields (B & 1011 G) which can alter the behaviour of matter down to the smallest scales. The extreme physics which occurs at NSs makes these objects of great ob-servational interest because such physical conditions can not be obtained here on Earth. NSs, therefore, provide us with information that can be used to test and mod-ify theoretical predictions of physics at extreme conditions (extreme densities, huge gravitational and magnetic fields).

In the last few decades the research on to these objects has increased vastly. Mean-while there are many (sub)-types of NSs estabilished (Harding, 2013), all of them important to our understandig of extreme physics. However, some of them are more easily detectable than others, such as X-ray pulsars (XRPs), which are highly magne-tized NSs in close binary systems. In these systems the NS can accrete matter of the companion star which in turn can release high energy X-ray radiation. The radiation can be studied here on Earth in order to help us understand the internal structure of NSs. The maximum luminosity that can be achieved in the case of spherical accre-tion is limited by the Eddington luminosity, which is LEdd≈ 2 × 1038_{erg s}−1_{for the} case of NSs. However, in case of highly magnetized NSs the luminosity can reach or even much exceed the Eddington limit (Basko and Sunyaev,1976; Mushtukov et al.,

2015a).

Nowadays there are three known pulsating ultra-luminous X-ray sources (ULXs) which have lumonisities exceeding the Eddington limit by more than two orders of magnitude: M82 X-2 (Bachetti et al.,2014), NGC 5907 ULX (Israel et al., 2017a) and NGC 7793 P13 (Israel et al.,2017b), which are all accreting NSs. The ULX in NGC 5907 is the brightest one among these three sources. Its maximum luminosity is detected to be ≈ 2 × 1041erg s−1 (Israel et al.,2017a), three orders of magnitude higher than the Eddington limit. At the moment there is no accurate model, which can describe all observational properties of ULXs powered by accreting NSs. How-ever, the model of accretion columns above highly magnetized NSs can provide an explanation of the basic features of NSs at super-Eddington luminosity (Basko and Sunyaev,1976; Mushtukov et al.,2015a; Mushtukov et al.,2015b). According to the model the accretion disk around the NS is interrupted at the magnetospheric radius and then the matter flows along the magnetic field lines on to the magnetic poles. For low mass accretion rates the matter reaches the NSs surface and creates hot spots of area S ∼ 1010cm2. In case of high mass accretion rates (> 1017g s−1)an accretion column is formed above the NS’s magnetic poles (see figure 1.1). The luminosity, above which the accretion column arises, is called the critical luminosity. It is of or-der of 1037erg s−1 and depends on surface magnetic field strength and geometry of

accretion channel (Mushtukov et al.,2015a).

FIGURE1.1: On the left: A highly magnetized NS with accretion columns above the mag-netic poles. The accretion column is a cylinder-shaped column confined by the magmag-netic field lines and its structure is defined by gravity and radiation pressure due to Compton scattering. On the right: Detailed overview of the column with a geometrical thickness of d, where H is the height of the column, x is the position inside the column, perpendicular to

the magnetic field line and ∆h is the advection length. (Mushtukov et al.,in prep)

At the top of the column the matter is halted due to a radiation dominated shock (Lyubarskii and Syunyaev,1982). The matter then slowly settles down through the magnetically confined column. The structure of the column is defined by gravity and the radiation pressure associated with Compton scattering, which is the main contributor of opacity inside the accretion column. The height of the column is de-pendent on the mass accretion rate and magnetic field strength and could be as large as the NS’s radius (Poutanen et al.,2013). It is suspected that the high luminosities originate from the accretion column (Lyubarskii and Syunyaev,1988; Mushtukov et al.,2015b).

Understanding the internal structure of accretion column is crucially important for estimations of the accretion column height and understanding of the physical pro-cesses which are going on inside the column. It is also important for our understand-ing of observational manifestation of accretunderstand-ing NSs at high mass accretion rates. However, the problem of the accretion column is incredibly complex: the magneto-hydrodynamical processes in extreme gravitational and magnetic fields have to be considered simultaneously with the radiative transfer. The radiative transfer calcu-lations are complicated by themselves because the basic processes of radiation and matter interaction are strongly dependent on photon energy, photon momentum and polarization (Harding and Lai,2006).

In this project the radiation inside the accretion column is described by using a sim-plified model, where we consider a layer with a fixed geometrical thickness of d and

is aligned with the magnetic field lines. We construct and discuss the equations of radiative transfer in this layer. We take into account the most important features of the radiative transfer in a strongly magnetized medium. We assume that advection is negligible in the accretion column, which is indeed the case, up to some luminos-ity value ((Mushtukov et al.,2015b), see also section1.2.2). The transfer equations are solved numerically and the results give us an idea about the distribution of the radiation of the layer and therefore inside the column. The results obtained will be also very important for theoretical predictions in the future. In a few years time the Imaging X-ray Polarimetry Explorer (IXPE) will be launched into an orbit around our Earth. The IXPE will exploit the polarization state of light of astrophysical ob-jects such as X-ray pulsars. The results of this project could be helpful in the inter-pretation of polarized radiation detected with IXPE.

This thesis has the following structure. First we consider the physics in high mag-netic fields and explain some decisions we made to simplify the radiative trans-fer equations. Then we continue with the physical parameters suitable to accretion columns and out of that we can already conclude some important features. In section

2, first the radiative transfer equations for the one - dimensional problem are solved to get used to the physical quantities associated with transfer equations in general. After that we expand the equations to the three - dimensional case and finally we take into account the photon polarizations. At the end of the thesis, in section3, the obtained results are presented and discussed.

1.1

Physics in a strong magnetic field

High magnetic fields alter the basic properties of matter and radiation (Harding and Lai,2006), which in an astrophysical setting could lead to various observable phe-nomena. Explaining every detail of the processes caused by high magnetic fields would be beyond the scope of this thesis. However, we consider the most impor-tant aspects that underlie some of the decisions made in this project. The behaviour of charged particles, photon polarization and photon propagation in high magnetic fields are such aspects which are important to consider.

1.1.1 Charged particles and cyclotron energy

The motion of charged particles, for example, electrons, is quantized in circular orbits when a magnetic field is applied. Each circular orbit corresponds to a dis-crete quantum state and these seperate states are called Landau levels. In the non-relativistic case the energy difference between the states for an electron equals the cyclotron energy, which is given by:

Ecyc= ~eB

mec = 11.6 B12keV, (1.1)

where ~ is the reduced Planck constant, e is the electron charge, B is the magnetic field strength (B12 = B/1012G), meis the electron mass and c is the speed of light. From this it can be seen that the cyclotron energy is dependent on the field strength. The cyclotron energy provides a useful scaling for the particle energies in a mag-netized plasma. A higher magnetic field results in a higher cyclotron energy. In the vicinity of highly magnetized NSs, the cyclotron energy would be extremely high, in case of XRPs it is in the X-ray or even gamma-ray energy range. If the magnetic field strength is extremely high, the photon energy is well below the cyclotron energy

even inside the accretion column. In this thesis we assume that the typical photon energies are below the cyclotron energy and most of the electrons are occupying the ground Landau level.

1.1.2 Photon polarization in magnetized plasma

Polarized radiation is commonly described by four Stokes parameters (Zheleznyakov,

1996). The radiative transfer in this case is given by a system of four integro-differential equations. However, the strongly magnetized plasma in the vicinity of NSs is anisotropic and birefringent (Harding and Lai,2006). The phase and group velocities of photons depend on their polarization state. As a result, the radiation can be described as a su-perpositions of two normal modes, which independently propagate in a magnetized plasma. The mixing of polarization states is possible only due to the state switches in scattering events. Generally, the normal modes in magnetized plasma are ellip-tically polarized. However, the parameters of elliptical polarization depend on the angle between photon momentum and B-field direction and local magnetic field strength. If the magnetic field strength is sufficiently high, the normal modes are almost linearly polarized in the directions across the strong magnetic field (Gnedin and Pavlov,1974). Along the field lines the modes are circularly polarized. The ex-act ellipticity of the normal modes depends both on magnetic field strength and local mass density of plasma. It is interesting that a magnetized vacuum is a birefringent medium by itself, but the normal modes there are linearly polarized.

In this bachelor thesis we consider two linearly-polarized normal modes only, O-and X-mode, where the electric vector of the O-mode is parallel to the plane formed by the photon momentum and magnetic field. The electric vector of the X-mode is perpendicular to the plane formed by the photon momentum and the field. The cross sections of Compton scattering in strong magnetic field are well approximated by the cross sections obtained for linearly-polarized normal modes if the photon en-ergy is taken to be below the cyclotron enen-ergy. For the directions along the field lines the modes are circularly polarized. However, the cross sections of the two normal modes are exactly the same in the directions along the field lines. If we are con-cerned about the size of the cross sections along the field lines, then we can obtain that by the two linearly-polarized normal modes. Therefore, the appearance of nor-mal modes in a strongly magnetized plasma allows us to simplify the description of radiative transfer: a system written for four Stokes parameters reduces to a system of two radiative transfer equations written in terms of the linearly-polarized normal modes.

It has to be mentioned that under certain conditions the birefringence of the magne-tized plasma and the magnemagne-tized vacuum compensate each other and a description of radiative transfer in terms of normal modes is not applicable anymore, this ef-fect is called “the vacuum resonance” (Gnedin, Pavlov, and Shibanov,1978). In this case the radiative transfer has to be described using the four Stokes parameters. This special situation, which might be important in some range of magnetic field strength and plasma mass density, is beyond the scope of this thesis.

1.1.3 Cross sections of linearly polarized modes

In the non-magnetized vacuum the photon scattering by electron is described by the Thomson cross section:

σT= 8π 3  e2 mec2 2 = 6.65 × 10−25cm2. (1.2) A strong magnetic field modifies the scattering cross section. Non-relativistic cross sections for the linearly-polarized normal modes strongly depend on photon energy and momentum and can be described by (Blandford and Scharlemann,1976):

σO= σT  sin2θB+1 2cos 2_θB  E2 (E + Ecyc)2 + E2 (E − Ecyc)2  (1.3) σX= σT 2  E2 (E + Ecyc)2 + E 2 (E − Ecyc)2  , (1.4)

where σO and σX are the cross sections for O-and X-mode respectively, σT is the cross section due to non-magnetized scattering and θBis the angle between the pho-ton momentum and the magnetic field direction. Both equations (1.3 and 1.4) de-pend on photon energy and cyclotron energy. However, the cross section of the X-mode does not depend on the angle between photon momentum and magnetic field direction. When scattering, the photons change their momentum and polariza-tion state (see secpolariza-tion2.3.1for detailed information).

In figure1.2 the total cross section is given as a function of photon energy for O-and X-mode photons which propagate initially across the magnetic field direction. This plot however is obtained with the full relativistic expression for O-and X-mode cross sections. In the plot some peaks can be seen, where the scattering cross section exceeds Thompson cross section by several orders of magnitude. The peaks are re-lated to electron transition between Landau levels and arise at certain photon ener-gies. The peaks are referred as the resonant features of the scattering. The difference between the full relativistic equation to obtain this plot and the equations (1.3 and

1.4) we are using are the number of peaks that can be shown. Nonetheless, in our case the photon energy is assumed to be lower than the cyclotron energy. That is done for simplicity, but it is reasonable assumption for the case of extremely high magnetic field strength. Therefore, we are concerned about the shape of the two graphs below the first resonant peak at E/mec2 _{∼ 0.09 in figure1.2. As can be seen} from the figure, for low energy photons the cross section for the X-mode is much smaller than the Thomson cross section. At the same time the scattering cross sec-tion for O-mode photons is almost equal to σT(at least for this particular direction of initial photon motion). This dramatic difference of cross sections is very important for the radiation field inside the accretion column.

FIGURE1.2: Cross section as a function of photon energy for the two polarized modes. The angle of direction of both modes is across the field ( θB = π₂) (Mushtukov, I.Nagirner, and

Poutanen,2016).

1.2

Physical parameters inside the column

Until now we have introduced the theory of physics in a strong magnetic field which has led to a simplification of the radiative transfer in terms of two linearly-polarized modes. Before we construct and discuss the transfer equations by them selves, we should consider the physical parameters appropiate for the theory of accretion col-umn, which then can be related to the transfer equations.

1.2.1 Optical thickness of the accretion column

The optical thickness τ is related to the absorption of the intensity, it is defined as a length measured in the length of photon free path. An optical thickness of τ = 0 means that there is no reduction in the intensity, whereas τ  1 means that the intensity would be significantly reduced due to absorption or scattering processes. The main mechanism of opacity in accretion columns is Compton scattering. To de-termine the optical thickness inside the accretion column we should first dede-termine the number density of electrons inside the column. The number density of the elec-trons is given by the local mass density, chemical composition and ionization degree of material. Because of extreme temperature inside the accretion column, we can assume that the matter is completely ionized. The local mass density is defined by the mass accretion rate and velocity according to the law of mass conservation. Let us consider that matter is accreted from the companion star, the speed of the matter is given by the freefall velocity:

vff(r) = r

2GM

r , (1.5)

where G is the gravitational constant, M is the mass of the NS and r is the distance between the infalling matter and the center of the NS. The freefall velocity of the matter thus increases towards the NS as the distance r decreases. However, the maximum velocity that can be attained is not at the NS surface, it would be obtained at a certain height above the NS:

vff(R + H) =

r 2GM

R + H, (1.6)

where R is the NS radius and H the height above the NS. In the scenario of high mass accretion rates, the accreting matter would not reach its maximum velocity at the NS’s surface because of the occurrence of a radiation-dominated shock at the top of the accretion column (Basko and Sunyaev,1976). The speed of the matter at the top of the column will be deaccelerated to roughly v(R + H) = vff(R + H)/7(Lyubarskii and Syunyaev,1988). Then the matter would slowly settle down, through the col-umn, to the surface of the NS where it reaches a speed of v(R) = 0. Now, consider a typical NS with mass M = 1.4M, R = 10 km and a height above the NS of 5 km, which is a typical height of the accretion column (Mushtukov et al.,2015b). Plug-ging these values into equation1.6gives the free-fall velocity at the top of accretion column vff(10 + 5) = 0.53c. Therefore, the velocity behind the radiation dominated shock can be roughly estimated as v(10 + 5) = v_ff(10 + 5)/7 = 0.075c. As a result, we know that the velocity in accretion column decreases from ∼ 0.075c to 0 at the NS surface. The exact velocity dependence is unknown and should be provided by an accurate model of accretion column. In this thesis we assume that the velocity of accreting material in accretion column is a linear function of height.

Now, we can relate the local velocity in the column to the number density of the electrons inside the column by considering the mass conservation law:

˙ M

2S = ρ v, (1.7)

where ˙Mis the total mass accretion rate, S is the area of cross section of the column, ρ is the mass density of the plasma, v is the speed of the plasma and the factor 2 is due to two columns at the NS. If the accreting plasma consists of fully ionized hydrogen, then: ne ≈ np. The number density of the protons is related to the mass density: np= ρ/mp, so substituting equation1.7gives:

ne(h) ≈ np= ρ mp

= M˙

2S v(h) mp

, (1.8)

but it could be also written in terms of β. As we know v ∈ [0, 0.075c] and β = v/c, thus β ∈ [0, 0.075]. Then:

ne( ˙M , β) ≈ ˙ M

2 S β c mp. (1.9)

The optical thickness is related to the electron number density: τ = σned, where d is the geometrical thickness of the column and σ is the scattering cross section. The cross section is given by σ = κemp, where κe is the opacity due to Compton scat-tering, combining these two expressions along with the expression of the number density1.9gives: τ ≈ M˙ 2 S β c mpκempd ≈ ˙ M κed 2 S β c ≈ ˙M17d4S −1 10β −1_, (1.10) where ˙M17 = ˙M /1017g s−1, S10 = S/1010cm2 and d4 = d/104cm. Because we are interested in supercritical mass accretion rates in a range 1017_{− 10}21_{g s}−1_{, the optical}

thickness due to non-magnetized Compton scattering is typically in the range τ ∈ [10, 1000].

1.2.2 Typical escape time of photons from the accretion column

We can give an estimation of the time of escape tescof photons emitted from inside the accretion column. If this time is significantly smaller than the time of plasma to settle down to the NS surface, we can ignore the effects of advection, otherwise the radiative transfer is strongly affected by hydrodynamics and cannot be considered separately. In order to determine the tesc, we have to get an expression for the aver-age number of scatterings which photons have undergone to get to the border. The expression for the average number of scatterings of the photons in the case of a layer of given optical thickness was obtained analytically (Nagirner,1972):

hni(τ, τ0) = 1

2ψ(τ )ψ(τ0− τ ), (1.11)

where τ is the optical depth where the photon was initially emitted, τ0 is the total optical thickness of the layer and ψ(x) is Sobolev’s function for conservative scat-tering. In the limit where τ and τ0 are both  1, we can also use the following approximation:

ψ(x) =√3(x + 0.71). (1.12)

In order to get the average number of scatterings for photons to get to the border of the accretion column we can take the average of the number of scatterings emitted from different optical depth:

hnib(τ0) = hhni(τ, τ0)iτ. (1.13) where the subscript b stands for border and hnib(τ0)is the average number of scat-terings for photons to get to the border for a fixed optical depth τ0 of the layer, the subscript τ stands for the different optical depth and hhni(τ, τ0)iτ is the average of the number of scatterings emitted from different optical depth. To solve equation

1.13we first get the expression for hni(τ, τ0)by substituting equation1.12into equa-tion1.11which gives:

hni(τ, τ0) = 3

2(τ + 0.71)(τ0− τ + 0.71), (1.14) now we can use equation1.13to rewrite an integral:

hni_b(τ0) = 3 2 Z τ0 0 (τ + 0.71)(τ0− τ + 0.71) τ0 dτ, (1.15)

solving the integral gives:

hnib(τ0) = 3 2  τ2 0 6 + 0.71τ0+ (0.71) 2  , (1.16)

and in the limit for big optical thickness (τ >> 1), this reduces to: hnib(τ0) ≈ 3 2  τ2 0 6  = 1 4  τ₀2  , (1.17)

Now that we have an expression for the average amount of scatterings to get to the border, we can determine the time of escape. Consider first the mean free path of the photon, that is the distance a photon travels before it collides with an electron. According to the definition of the optical thickness, it is simply given by:

l = d τ0

. (1.18)

We further approximate the speed of the photon in the plasma equal to the speed of light in vacuum c. Then the tesc is the distance divided by the speed:

tesc= s v = hni_b(τ0) · l c = hn_bi(τ₀) · d c · τ0 , (1.19)

using the exppresion for the average amount of scatterings using equation1.17gives: tesc = 1 4  τ₀2  · d c τ0 = 1 4  τ0d c  . (1.20)

A typical value for d is ∼ 104_{cm, which along with an optical thickness of τ = 1000} gives a time of escape of tesc ≈ 10−5_{s. The time for the matter to settle down from a} height of 5 km on to the NS surface is 10−4s, so tesc  tsettleand we can ignore the effects of advection.

1.2.3 Length of advection

Although we can ignore the effects of advection, we have found an interesting fea-ture which simplifies a lot of models where advection is taken into account. Let us denote as h1 the height of the matter above the NS where it initially emits a photon and denote h2 as the height of the matter when the initially emitted photon has es-caped the column. Then we can denote the height difference as ∆h = h1− h2, which we refer as the advection length. The advection length in the accretion column is independent of the speed of the plasma inside the column. We can show that by using the expression for the time of escape1.20:

∆h = tescc β = 1 4  τ0d c  c β = 1 4  ˙ M17d4S10−1β−1d c  c β = 1 4  ˙ M17d4S10−1d  , (1.21)

so in order to estimate the advection length ∆h we do not have to know what the speed is inside the column. ∆h can be obtained by the physical parameters of the column it self and the mass accretion rate.

1.2.4 Initial sources of emission

The initial sources of emission which produce the radiation inside the column are as-sociated to the mass accretion. As the matter sinks through the column, its potential energy is converted to radiation. The local energy release (Q+_{) is defined by height} derivative of the kinetic and potential energy flux. Because the potential energy flux through the unit area is ˙mghand the kinetic energy flux is ˙mv2_{/2, we get:}

Q+(h) = ˙m  GM (R + h)2 + c 2_βdβ dh  = 1024m˙10  1.328 m (R + h)2₆ + 9β dβ dh6  erg cm−3s−1, (1.22) where ˙m10= ˙m/1010g cm−2s−1is a mass accretion rate onto a unit area ( ˙m = ˙M /S). The initial emission coefficient is directly related to the local energy release:

ε(0)(h) = Q +_(h)

4π . (1.23)

In case of ULXs powered by NSs the mass accretion rate per unit area ˙m10∼ 1 and, therefore, the initial emission coefficient is of order of ε(0)∼ 1023_{erg cm}−3_s−1_ster−1_. In this thesis calculations are performed in dimensionless form, when the initial emissivity ε(0) is taken to be 1 everywhere within the layer. Actual initial emissivity depends on mass accretion rate and dynamical parameters within the column (see equation1.22). As a result, the final radiation energy density is obtained in units of

1 3 10 17_d4ε(0) 23 ' 2.6 × 1016d4m10˙  1.328 m (R + h)2₆ + 9β dβ dh6  erg cm−3s−1, (1.24) which depend on physical parameters in the system: mass accretion rate per unit area, geometrical thickness of accretion channel, NS mass and radius, local velocity of matter and its derivative.

2

Radiative transfer equations

2.1

One-dimensional case

Before solving the radiative transfer equations for the three-dimensional polarised case, we take a look at the one-dimensional case to introduce the problem and under-stand the method of numerical solutions. In the one-dimensional case photons are assumed to be able to propagate only in to opposite directions perpendicular to the edge of the accretion channel, the radiation field depends only on the x-coordinate. We denote the direction of photon propagation as "positive" if the photons are prop-agating from smaller to larger x-coordinates and "negative" if they are propagat-ing from larger to smaller x-coordinates. Therefore, we have two seperate radiative transfer equations which describe the variations of the intensity accordingly:

dI+(x) dx = −αI+(x) + ε+(x) + ε (0) + (x) (2.1) dI−(x) dx = −αI−(x) + ε−(x) + ε (0) − (x) (2.2)

where I+(x), I−(x)are the intensities in the positive and negative direction corre-spondingly, ε+(x), ε−(x)are the emission coefficients due to the scattering, ε(0)+ (x), ε(0)₋ (x)are the emission coefficients due to initial sources of radiation and α is the ab-sorption coefficient. In our calculations of the radiative transfer the initial emission coefficient is caused by the mass accretion, as mentioned in the previous chapter (see section1.2.4). The absorption coefficient will be a consequence of altering the direction of the intensity and no real absorption will take place during scattering pro-cesses. In the one-dimensional case the absorption coefficient would be a constant for each scattering, so α = constant. In the polarized case the absorption coefficient αwould depend on the angle between photon momentum and magnetic field direc-tion.

Because we can take the absorption coefficient as a constant, there is a relation be-tween the emission coefficients, absorption coefficients due to the scattering and lo-cal intensities (the energy conservation law):

ε+(x) + ε−(x) = α(I+(x) + I−(x)). (2.3) If the direction of scattered photons does not depend on its initial direction then ε+= ε−, and we get:

ε+(x) = ε−(x) = α

2.1.1 Numerical method

In order to calculate the intensity distribution throughout the line we assume a line with geometrical length d, which is located within the x-coordinates: x ∈ [0, d]. The intensity distribution in the whole line changes after each scattering. At each point x in the line the intensity will decrease after each consecutive scattering until the intensity would be a constant. Numerically we have to calculate the intensity in the whole line of each scattering seperately. The initial emission coefficients are given by ε(0)₊ (x) = ε(0)₋ (x) = C1. With a suitable integrating factor the differential equations

2.1and2.2can be solved. Which will lead to:

I₊(0)(x) = x Z 0 dx0ε(0)₊ (x0) exp[−α(x − x0)], I₋(0)(x) = d Z x dx0ε(0)₋ (x0) exp[−α(x0− x)], (2.5) which are the intensities in both directions of photons which have not been scattered at all. From these equations the intensity at the arbitrary point x can be calculated. Notice that we are integrating over the emission coefficients, because the emission coefficients at each point of the line contribute to the intensity at the arbitrary point x. Of course, in the calculations one should not calculate just one arbitrary point x, but every x in the entire line. If I₊(0)(x)and I−(0)(x)are calculated for the entire line we can get the emission coefficients after the first scattering:

ε(1)₊ (x) = ε(1)− (x) = α 2  I₊(0)(x) + I−(0)(x)  , (2.6)

and if the emission coefficients ε(1)₊ (x) and ε(1)− (x) are known, then the intensities after first scattering can be calculated:

I₊(1)(x) = x Z 0 dx0ε(1)₊ (x0) exp[−α(x − x0)], I−(1)(x) = d Z x dx0ε(1)− (x0) exp[−α(x0− x)]. (2.7) and emission coefficients which have been scattered twice:

ε(2)₊ (x) = ε(2)₋ (x) = α 2



I₊(1)(x) + I₋(1)(x), (2.8) and so on. These iterations will take place until the intensity doesn’t change any-more, so after the nth-scattering the total local intensity is given by

I±(x) = N X

n=0

I_±(n)(x), (2.9)

while the total emission coefficient is given by ε±(x) =

N X

n=0

ε(n)_± (x), (2.10)

2.1.2 Numerical results

To calculate the intensity and emission distributions a code in C has been constructed, the integrals in the equations of2.5were solved with the trapezoid method. The plots were made using matplotlib. The reason to use C is that the calculations of the radia-tive transfer equations takes a lot of computational time, so efficient code is neces-sary.

The results of our calculations are illustrated in figure 2.1. In figure 2.1a the dis-tributions of the total emission (using equation 2.10) are presented, and in figure

2.1bthe total intensity distributions (using equation2.9) in the positive direction for different values of the optical depth τ of the line (figure 2.1b). The geometrical length d = 10, which is taken arbitrary. In figure2.1ait can be seen that the total emission gets higher for higher optical thickness. This is reasonable because for higher opti-cal thickness more scatterings occur in the line. If there are more scatterings it takes more time for the intensity to leave the system and, thus, resulting in a higher emis-sion. The same explanation applies for the intensities which can be seen in figure

2.1b. The intensities are higher for higher optical thickness because for higher op-tical thickness it takes more time for photons to leave the system. Notice however that for I+ the intensities are 0 at the left border. This is reasonable because at the border the intensity is leaving the system. And at the left border only I− would be contributing for the intensity leaving the system. At the right border we can see that the intensity is the same for all optical thickness, due to conservation of energy: the energy initially emitted in the line should be the same as the total intensity leaving the system.

(A) Distributions of the total emission coeffi-cients in the line for different optical thickness. For higher optical thickness more scatterings oc-cur which leads to higher emission coefficients.

(B) Distributions of the total intensities in the positive direction for different optical thickness. For higher optical thickness there are more scat-terings which results in higher peaks. Also for I+it can be seen that there is no contribution at

the left border. At the right border the intensity is the same for all values of the optical thickness.

FIGURE2.1: Distribution of total emission and intensities.

In figure2.2the total intensity in the center as a function of the optical thickness is shown taking into account only a limited number of scatterings. It can be seen that

for high optical thickness one should consider more scatterings. From the previous figures it is known that higher optical thickness results in higher intensities.

FIGURE2.2: Total intensities in the center as a function of the optical thickness for different number of scatterings taking into account. It can be seen that more scatterings are needed for higher optical thickness, and if the number of scatterings is insufficient, then the intensity

is strongly underestimated.

If the number of scatterings is insufficient in the numerical calculations, then the intensity is strongly underestimated. In the subsequent chapter it will be explained more thoroughly how a sufficient number of scatterings is taken into account for various optical depth.

2.2

Three-dimensional: non-polarized

In this chapter we expand the radiative transfer equation from one-dimensional to three-dimensional. For the three-dimensional case we consider a slab of length d which not only radiates into the x-direction but as well in the θ and ϕ direction. We have chosen a right-handed coordinate system as the reference frame of the layer. The x-coordinate is perpendicular to the layer. The θ angle is the angle between the photon momentum and the positive x-axis, where the slab with thickness d is aligned with the x-axis. The angle ϕ is the angle between the positive y-axis and the projection of the photon momentum in the yz-plane. The reference frame of the layer is illustrated in figure2.3:

FIGURE2.3: The reference frame of the layer with the photon angles θ and ϕ. The magnetic field is directed in the positive z-direction.

where the magnetic field is in the positive z-direction. With this choice of coordinates the radiative transfer equation for isotropic scattering (where the absorption and emission coefficients do not depend on the direction) takes the following form:

cos θdI(x, θ)

dx = −αI(x, θ) + ε(x) + ε (0)_(x),

(2.11) Note, that the intensity depends on x and θ only because of plane symmetry of the problem. This equation can be solved again with a suitable integration factor to obtain the intensity for a specific scattering:

I(i)(x, θ) = x2 Z x1 dx0 ε (i)_(x0₎ cos θ exp  −α|x − x 0_| cos θ  , (2.12)

where x1 and x2 depend on the angle θ and therefore take different values for par-ticular cases (see under for further explanation of the values of x1 and x2). For the emission coefficient we have to take a double integral in the polar angles θ and ϕ. Just like the one-dimensional case the emission coefficient is only dependent on x, so for the three-dimensional case we get:

ε(i+1)(x) = α 4π Z 2π 0 dϕ0 Z π 0 dθ0sin θ0I(i)(x, θ0), (2.13) and because the absorption α does not depend on the angles it can be taken out of the integral. The intensity does not depend on ϕ because of plane symmetry, so the integral over dϕ yields 2π and therefore the emission coefficient for a specific scattering is given by:

ε(i+1)(x) = α 2

Z π 0

dθ0sin θ0I(i)(x, θ0), (2.14) The numerical approach is very similar to the one-dimensional problem. From the initial emission coefficients the initial intensity can be calculated using equation2.12. Then we calculate the emission (with equation2.14) after the first scattering using the initial intensity which is already obtained, this computational scheme could just be repeated until nth-scattering. However, for the three-dimensional case we will obtain the nth-scattering differently than for the one-dimensional case. In order to take into account a sufficient number of scatterings we will use the energy conser-vation law. Based on the fact that the initial amount of energy emitted in a layer should be equal to the amount of energy emitted from the borders of a layer. This only happens after a sufficient number of scatterings. The equation is given by:

2 Z d 0 dx ε(0)(x) = Z π/2 0

dθ sin θ| cos θ| I(d, θ) + Z π

π/2

dθ sin θ| cos θ| I(0, θ), (2.15)

where I(d, θ) = Pn

i=0I(i)(d, θ) refers to the total intensity at the right border and I(0, θ) =Pn_i=0I(i)_{(0, θ)refers to the total intensity at the left border. When the} con-servation law is satisfied with sufficient accuracy we can stop the calculations in the code. Furthermore, notice that in this particular problem there is no such thing as I+ or I−, however numerically there is some resemblance to these quantities. To calcu-late the emission2.14at some x an integral over θ has to be taken. For θ ∈ [0, π/2) the limits of integration of2.12should be from 0 to x and for θ ∈ (π/2, π] the limits in

2.12should be from x to d. For the angle θ = π/2 one can see that we would divide by zero due to the cos-term in equation2.12, but the intensity for this particular case can be obtained analytically.

After the total intensity has been calculated at each point we can get the temper-ature at each point inside the column. First we can relate the total intensity to the energy density, which is:

u(x) = 1 c Z 2π 0 dϕ0 Z π 0 dθ0sin θ0I(x, θ0), (2.16) where c is the speed of light, after that we can get the temperature from the energy density: u(x) = aT4_{ef f}(x) → Tef f(x) =  u(x) a 1₄ , (2.17)

2.3

Three-dimensional: polarized

The radiative transfer equations for the polarized case are quite similar to the non-polarized case in three-dimensions, but we have to take into account the two differ-ent polarization modes. The normal modes have differdiffer-ent group and phase veloci-ties and, as a result, must be described separately. Thus, we consider two seperate sets of radiative transfer equations. One corresponds to O-mode and the other one to X-mode. The polarized radiative transfer equation in terms of vectors and matrices is:

cos θdI(x, θ, ϕ)

dx = −α∗(θ, ϕ)b I(x, θ, ϕ) + ε(x) + ε

(0)_{(x), (2.18)}

where the notations for the intensities, emission and absorption coefficients are merged like so: I =  I1 I2  , ε =  ε1 ε2  , α∗_b =  α11+ α21 0 0 α12+ α22  , (2.19) and where I1 and I2 are the intensities in the two polarization modes and ε1and ε2 are the emission coefficients in the two modes. Subscript 1 corresponds to O-mode and subscript 2 corresponds to X-mode. α∗b is a matrix containing total absorption coefficients for both modes respectively. If we want to know the absorption coeffi-cient corresponding to mode 1, then we have to consider what the contribution is to this absorption coefficient. This contribution may come from the possibility of switching from one mode to another. For example, consider the intensity of mode 1. The photons of mode 1 may scatter, and after scattering the photon can either stay in mode 1 or switch to mode 2. In both cases the scatterings are reducing the intensity of mode 1. The general way of describing these absorption coefficients is: αi = α1i+ α2i, where αiis the total absorption coefficient for mode i, α1iis the ab-sorption coefficient where a photon in mode i scatters into mode 1 and α2i is the absorption coefficient where a photon of mode i scatters into mode 2. We can also define the probability of switching from mode "j" to mode "i" due to the scattering:

Pij = αij α1j+ α2j

≤ 1. (2.20)

Notice however that the absorption coefficients depend on the photon angles θ and ϕand hence so does the probability.

Numerically the intensity of each separate scattering can be described by the inten-sity vector: I(i)(x, θ, ϕ) = x2 Z x1 dx0 exp  − bα∗(θ, ϕ)|x − x 0_| cos θ  ε(i)_(x0₎ cos θ , (2.21)

where x1 and x2 are the positions inside the layer and depends on the θ angle, the same as in non-polarized case. Notice that the argument of the exponent under the integral in equation 2.21is a matrix. In this case the exponent is also a matrix of the same dimension, which can be obtained as a corresponding power series. The expression for the emission vector of each separate scattering is given by:

ε(n+1)(x) = 1 4π 2π Z 0 dϕ0 π Z 0 dθ0 sin θ0α(θ_b 0, ϕ0)I(n)(x, θ0, ϕ0), (2.22) whereαbis: b α =  α11 α12 α21 α22  . (2.23)

The intensities are dependent on three variables: x, θ and ϕ, whereas the emission coefficient is still dependent only on the x-coordinate. If we take a look at equation

2.22then we see the expressionα. The top row ofb αbconsists of the absorption coef-ficients α11and α12. Meaning that these coefficients would contribute to ε1, as these coefficients both describes photon modes scattered into mode 1. The second row of b

αon the other hand, consists of absorption coefficients contributing to ε2.

The computational scheme for the polarized case is the same as the previous cases of one-dimensional and three-dimensional non-polarized. From the initial emission coefficients we get the intensity of non-scattered radiation all over the layer, then we can get the emission coefficients due to the first scattering and intensity of radiation scattered only once, then we get emissivity after the second scattering and so on. However, what is substantially different from the previous cases are the following parameters {α1, α2, P11, P22}, which describe absorption coefficients for different modes and probabilities of polarization switch due to the scattering. We have to know these four parameters in order to calculate the radiative proccesses for the po-larized case. If these parameters are known we can also deduce the other absorption parameters which are contained in the absorption matrices α∗b and α. α1b and α2 can be obtained from the expressions of the O-and X-mode scattering cross sections, already discussed in the introduction (see section1.1.3). The absorption coefficients are related to the scattering cross sections as αi= σine, where neis the electron num-ber density. Therefore, the absorption coefficients can be represented as follows:

α1 = αT  sin2θB+ 1 2cos 2_θ B  E2 (E + Ecyc)2 + E 2 (E − Ecyc)2  , (2.24) α2 = αT 2  E2 (E + Ecyc)2 + E2 (E − Ecyc)2  , (2.25)

where αT = σTneis the absorption coefficient due to Thomson scattering and σT = 6.65 × 10−25cm2is the Thomson cross section, E is the photon energy and Ecycis the cyclotron energy. Note that α1, the absorption coefficient of the O-mode, depends both on photon energy E and the angle between photon momentum and magnetic field direction θB, while α2, the absorption coefficient of the X-mode, depends on the photon energy only.

The other parameters P11and P22are approximately given by (Miller,1995):

P11≈ 1 −  E Ecyc 2 , P22≈ 3 4. (2.26)

However, before we can perform the numerical calculations one should notice that the absorption coefficients contains the magnetic field angle θB. The equations for radiative transfer contain photon angles θ and ϕ, so in order to perform the calcula-tions we should express the magnetic field angles in terms of photon angles. Both

reference frames are given in figure2.4 but let us visualize the magnetic reference frame first in figure2.4a:

(A) The right-handed magnetic reference frame where the magnetic field line is in the positive z-direction. θBis refered as the angle between

the magnetic field direction and the photon mo-mentum. The geometrical length of the layer is

aligned with the x-axis.

(B) The right-handed reference frame of the layer where the magnetic field line is in the pos-itive z-direction. θ is the angle between the posi-tive x-axis with the photon momentum. The ge-ometrical length of the layer is aligned with the

x-axis.

FIGURE 2.4: The reference frames of the magnetic field and of the layer.

where we have used a right-handed coordinate frame with the x-axis perpendicular to the layer and the z-axis alined with the magnetic field direction. The Cartesian coordinates of the unit vector of the magnetic reference frame is given by:

     x = sin θBcos ϕB y = sin θBsin ϕB z = cos θB.

Figure2.4brepresents the reference frame related to layer. The Cartesian coordinates of the unit vector could be expressed in the polar angles θ and ϕ as follows

     x = cos θ y = sin θ cos ϕ z = sin θ sin ϕ,

so we have obtained expressions of the Cartesian coordinates in terms of the mag-netic field angles θB, ϕB and in terms of the photon angles θ,ϕ. With these expres-sions we can change from the magnetic reference frame to the reference frame related to the layer by equating the expressions of the two coordinate frames to each other. For example, notice that in the expressions of the absorption coefficients2.24there is a cos2θBterm. Looking at the expressions obtained from the reference frames we can equate the Cartesian z-coordinate of both frames to each other which gives the

following: cos θB = sin θ sin ϕ, so that would give us cos2_{θB = sin}2_{θ sin}2_{ϕ. For the} sin2θBterm in equation2.24we get: sin2θB= 1 − sin2θ sin2ϕ.

2.3.1 Scattering cross sections

In order to get a detailed description of the scattering one can use the differential cross section. The differential cross section describes not only the probability of scat-tering of a photon of given initial parameters, but also the probability of the final parameters of scattered photons (final photon momentum and polarization state). In this thesis we use the cross sections derived by (Herold,1979), who used methods of quantum electrodynamics and the assumption of electrons occupying the ground Landau level. The differential cross section dσij/dΩ0_B is related to the complex am-plitude of scattering aij as follows (Herold,1979):

dσij dΩ0_B =

3

32πσT|aij|

2_{, (2.27)}

where σij is the cross section of a photon which initially was in state j and after scattering it is in state i and ΩBgives the final direction of photon momentum. The complex amplitudes are dependent on the photon momentum, photon energy and the polarization. For the case of non-relativistic scattering the complex amplitudes are given by:

a11= 2 sin θBsin θ0B+ cos θBcos θ0B  E E + Ecyc ei(ϕB−ϕ0B)+ E E − Ecyc e−i(ϕB−ϕ0B)  a22= E E + Ecyce i(ϕB−ϕ0B)+ E E − Ecyce −i(ϕB−ϕ0B)  a21 a12  =  E E + Ecyc ei(ϕB−ϕ0B)− E E − Ecyc e−i(ϕB−ϕ0B)   −i cos θB i cos θ_B0  , where the angles with an accent ϕ0_B,θ0_Bare the final angles in the scattering process and the angles without an accent ϕB,θB are the initial photon angles. In figures

2.5and2.6 the differential cross sections (dσ11/dΩ0_B) σ₁−1 and (dσ21/dΩ0_B) σ−1₁ , are displayed as a function of the final θ_B0 angle, where we have divided by the total cross section. The plots contain information about probabilities of possible changes in polarization state of photons due to the scattering. These plots are polar plots where the radial and polar coordinate represents the differential cross section and final magnetic field angle θBrespectively. The differential cross sections for photons with different photon energies, in terms of the cyclotron energy, and different initial magnetic field angle θBare given in both of the figures2.5and2.6. The initial angle ϕB = 0is fixed for all cases.

FIGURE2.5: The distribution of the differential cross section (dσ11/dΩ0B)σ −1

1 for several

dif-ferent photon energies and initial θBangles.

In figure2.5the distribution of the differential cross section (dσ11/dΩ0_B)σ−1₁ for sev-eral different photon energies and initial θB angles is shown, the direction of the magnetic field is as well illustrated and pointed to the right. It could be seen that photons propagating along the magnetic field would likely still propagate in the same direction after scattering. The probability distribution for the photons prop-agating along the magnetic field angle is the same for the two different energies of E = 0.01Ecycand E = 0.5Ecyc. For the inital angle of π/4 however, the probability of scattering for low energy photons is higher than for high energy photons. But the final direction for photons with initial angle of π/4 is most likely perpendicu-lar to the magnetic field. Notice further that in figure2.5three lines are coinciding: θB = π/2 E = 0.01Ecyc, θB= π/2 E = 0.5Ecycand θB= π/4 E = 0.01Ecyc.

FIGURE 2.6: The distribution of the differential cross section (dσ21/dΩ0B) σ −1

1 for several

In figure2.6the distribution of the differential cross section (dσ21/dΩ0_B)σ₁−1is shown, where we can obtain the probability of scattering from mode 1 to mode 2. We con-sider first the photons propagating along the magnetic field line. It could be seen that the probability of scattering into another direction is equal for all angles. If we would compare the photons propagating along the field line from this figure2.6

with figure2.5, it can be seen that the probability of switching is equal in the final direction along the field line. However, if the final direction is not along the mag-netic field line, it would most likely switch to mode 2. If we compare the photon with initial angle π/4 and a energy of E = 0.5Ecyc from figure 2.6with figure 2.5, it can be seen that this photon would most likely stay in mode 1 but the direction would be altered perpendicular to the field line. Notice that in figure2.6there are less plots compared with figure2.5. These plots are left out because the probability was almost equal to zero. Thus for the other plots of different angles and different energies, they would stay in mode 1 after scattering.

From the expression of the differential cross section2.27we can obtain an expres-sion for the total cross section by integration over parameters of the final photon (i.e. direction of its momentum and polarization state). Remember from the expression of the absorption coefficient we had to take a sum of the contributing parts to the absorption coefficient, for example: αi = α1i+ α2i. If we want to calculate the total cross section for mode i, the sum of the contributing parts of the differential cross sections has to be taken. The total scattering cross section is then given by:

σj(θB, ϕB, E/Ecyc) = 2 X i=1 π Z 0 dθ0_B 2π Z 0 dϕ0_Bsin θ0_Bdσij dΩ0_B, (2.28) where dσij/dΩ0_Bis the differential cross section. The total cross section is a function of the initial angles, polarization state and photon energy. In figure2.7the total cross section is given for mode 1 (O-mode) and 2 (X-mode) as a function of the photon energy for various initial angles θB. The total cross section does not depend on the angle ϕB.

From figure2.7it can be seen that the total cross sections for O-and X-mode pho-tons depend on their photon energy. In both plots (left and right) of the figure the first resonant feauture (at E/Ecyc = 1) of the scattering can be seen where the cross section exceeds Thompson cross section by several orders of magnitude. In the in-troduction (see section1.1.3) it was already shown that the cross sections for O-and X-mode photons strongly differ in case of low energy photons. From figure2.7it can be additionally seen that the cross sections of O-mode photons strongly depends on the initial propagation angle with the magnetic field as well. If O-mode photons propagate in the direction of the magnetic field, the cross section, in case of low en-ergy, can be as low as the cross section for X-mode photons with the same photon energy. For initial directions across the magnetic field, the cross sections of O-mode photons are almost equal to the Thompson cross section. The cross sections for X-mode photons however do not depend on the initial propagation angle. Which is reasonable, based on the physical properties of X-mode photons (see equation1.4). The radiation inside the accretion column is strongly affected by the properties of cross sections for O-and X-mode photons.

FIGURE2.7: On the left: Total cross sections for O-mode as a function of photon energy for various initial direction of photon momentum. The cross sections not only depend on the photon energy but as well on the initial direction. On the right: Total cross sections for X-mode as a function of photon energy for various initial direction of photon momentum. The cross sections only depend on the photon energy and not on the initial direction. Both plots

3

Results

In the following results the initial emission coefficient is taken as a constant which is the same throughout the layer. For the non-polarized case the initial emission coefficient is ε0 = _4πd1 , where d is the geometrical length of the layer and taken to be 1. The initial emission coefficient for the polarized case depends on the initial distribution of O-and X-mode. In the case of equal initial emissivity for the two modes then the initial emission coefficient is given by: ε0 = _8πd1 , with d taken to be 1. If we are concerned about real physical values in the final result of any of the plots then we would simply multiply the dimensionless results with the constant obtained in section1.2.4. For now the plots contains dimensionless results, because at this point we are concerned about the relationship among quantities.

3.1

Non-polarized

In figure3.1the total intensity distributions of the non-polarized case are illustrated for various optical thickness. According to equation2.12the intensity depends on the position inside the layer and on the propagation angle θ of the photon. The presented plots are heatmaps where the colour corresponds to the strength of the di-mensionless intensity, the intensity is plotted in logarithmic scale and as a function of x and θ.

From the plots in figure3.1 it can be seen that for higher optical thickness the in-tensity is higher, which was also the case for the one-dimensional problem. The intensity is higher for higher optical thickness because for higher optical thickness it takes more time for photons to leave the system. Due to the vast increase of the intensity for higher optical thickness, the θ dependence becomes less visible as well. Notice furthermore that the intensity leaving the system at the right border (which can be seen best in the plot for optical thickness τ = 10) is higher for the angles θ ∈ [0, π/2]than it is for θ ∈ [π/2, π]. This should be the case because the intensity at the border should be higher for photons leaving the system than for photons going into the system.

FIGURE3.1: Total intensity distribution of the non-polarized case as a function of the photon angle θ and position x inside the layer for various optical thickness, τ = {10, 100, 200, 500}. The total intensity strongly depends on the position inside the layer, where the intensity is at it highest in the center and drops towards the borders. It depends on the angle as well, which can be seen best in the plot with the optical thickness of τ = 10. For higher optical

thickness the θ dependence becomes less visible.

Thus for higher optical thickness the intensity gets higher due to the increase of scatterings inside the system. The increase of scatterings leads to a higher photon energy density as well. In figure3.2the photon energy density for various optical thickness is illustrated. In order to obtain the photon energy density we have used equation2.16. As can be seen out of the figure, the photon energy density is higer for higher optical thickness.

FIGURE 3.2: Photon energy density distribution as a function of the position x inside the layer for various optical thickness. For higher optical thickness the photon energy density is higher due to the increase of scatterings and thus photons are locked up in the layer for

3.1.1 Average number of scatterings

To check the validity of our numerical calculations we can compare the average number of scatterings, obtained numerically, with the expression of the analytical approximation1.17. As a reminder, the average number of scatterings are the scat-terings which a photon would have undergone since initially emitted. The expres-sion of the average number of scaterings obtained numerically is given by:

hn(x, θ)i = N X i=0 iI(i)(x, θ) Itotal(x, θ), (3.1)

where i is the number of each separate scattering and capital N is the N th-scattering. In figure3.3 the average number of scatterings inside the layer is illustrated for an optical thickness τ = 10 and τ = 500. The plots in this figure are heatmaps where the average number of scatterings are plotted as a function of the photon angle θ and position x inside the layer in logarithmic scale. For an optical thickness τ = 500the θ dependence becomes less visible due to the rapid increase of the number of scatterings. In figure3.4 the average number of scatterings as a function of the optical thickness is shown. Where we compared the numerical solutions with the analytical results using equation1.17. It can be seen that our numerical solutions corresponds well to the analytical solutions.

FIGURE3.3: The average number of scatterings as a function of the photon angle θ and the position x inside the layer for various optical thickness, τ = {10, 500}. The average number of scatterings are the scatterings a photon would have undergone since initally emitted. The dependence on θ becomes less visible for higher optical thickness due to the rapid increase

FIGURE3.4: The average number of scatterings as a function of the optical thickness. The numerical results (red dotted line) are compared with the analytical approximation (dashed

blue lines).

3.2

Polarized

Next we consider the radiative transfer problem for the case of polarized radiation which can be described in terms of two normal modes. We consider the total in-tensities of O-and X-mode and they both depend on three variables: x, θ and ϕ coordinates, which gives a direction of photon momentum. From the expressions of the cross sections (equation1.3 and 1.4) it can be seen that they depend on the photon energy (in units of the cyclotron energy) as well. Let us first make a distinc-tion in the energy, we consider photons which have an energy of E = 0.1Ecyc or E = 0.5Ecyc. For both of these energies we take a look at the total intensities in both polarization modes, at the border and in the center of the layer. The presented plots are heatmaps where the colour corresponds to the strength of the dimensionless in-tensity, the intensity is plotted in logarithmic scale and as a function of the photon angle ϕ and θ. All of the plots shown in this section are for fixed optical thickness of τ = 100, due to non-magnetized Compton scattering.

In figure 3.5 the total intensity distributions of O-and X-mode at left the border and in the center with an energy of E = 0.1Ecycis illustrated. The rows of the figure represents O-and X-mode respectively and the columns of the figure represents the border and center respectively.

FIGURE3.5: The total intensity distributions of O-and X-mode photons at the left border and in the center as a function of the photon angles ϕ and θ. The energy of photons is E = 0.1Ecycwithin a layer of optical thickness τ = 100. (a) O-mode at the left border. (b)

O-mode in the center. (c) X-mode at the left border. (d) X-mode in the center.

From figure3.5(a and b) it can be seen that for O-mode photons the intensity strongly depends on both photon angles ϕ and θ. This is due to the specific features of the scattering cross section in strong magnetic field. The intensity for O-mode photons leaving the system at the left border (figure3.5a) is higher for the directions close to the magnetic field lines. In the center, the intensity for O-mode photons (figure3.5

b) is sharply peaked in the directions along the magnetic field lines. This is because the cross sections for O-mode photons along the magnetic field lines is at its smallest value (see figure2.7in section2.3.1). For X-mode photons however (figure3.5c and d), the intensity does not depend on the ϕ angle, but it still does depend on the θ angle. The direction of the intensity for X-mode photons leaving the system is not strongly peaked in the directions along the magnetic field lines, as it is for O-mode. Furthermore, it can be seen that the intensity in the center is in overall higher than at the borders and that the intensity of O-mode photons is higher than that for X-mode photons, both at the border and in the center.

In figure 3.6 the total intensity distributions of O-and X-mode at left the border and in the center with an energy of E = 0.5Ecycis illustrated. The rows of the figure represents O-and X-mode respectively and the columns of the figure represents the border and center respectively.

FIGURE3.6: The total intensity distributions of O-and X-mode photons at the left border and in the center as a function of the photon angles ϕ and θ. The energy of photons is E = 0.5Ecycwithin a layer of optical thickness τ = 100. (a) O-mode at the left border. (b)

O-mode in the center. (c) X-mode at the left border. (d) X-mode in the center.

From figure3.6it can be seen that the general shape of all of the plots are the same as in figure3.5. The main difference however is that the intensity distributions of the plots in figure3.6are less peaked (e.g. the difference in the highest and lowest intensities are smaller) compared with the intensity distributions in figure3.5. 3.2.1 Initial emissivity distributions of O-and X-mode photons

Until now we assumed that the initial emissivity of O-and X-mode photons was the same throughout the layer. In a real physical setting we do not know whether the distributions are equal to each other or not, unless we have an accurate model de-scribing the mechanisms of initial radiation in details. In figure3.7the total intensity distributions of final O-mode photons in the center of the layer for various initial dis-tributions of photon modes is shown, to determine whether the initial disdis-tributions would make a difference in our results. Therefore, in figure3.7, we have assumed how the final distributions of O-mode photons would be if initially the layer con-sisted of only O-mode photons or of X-mode photons. We have done this for two different energies as well: E = 0.1Ecycand E = 0.5Ecyc.

FIGURE3.7: The total intensity distributions of O-mode photons in the center for various initial distributions of photon modes as a function of the photon angles ϕ and θ. (a) Initially only O-mode photons with an energy of E = 0.1Ecyc. (b) Initially only O-mode photons with

an energy of E = 0.5Ecyc. (c) Initially only X-mode photons with an energy of E = 0.1Ecyc.

(d) Initially only X-mode photons with an energy of E = 0.5Ecyc.

If the plots (b) and (d) from figure3.7 are compared with plot (b) from figure 3.6, it can be seen that there are no differences in the total intensity. If plot (a) from figure3.7is compared with plot (b) from figure 3.5no big differences can be seen. However, if we compare plot (c) from figure3.7with plot (b) from figure3.5we see a huge difference in the intensity, the intensity is less high in the directions along the magnetic field lines but the intensity is still sharply peaked.

3.2.2 Ratio of the photon energy density in the center

To determine the ratio of photon modes in the center we have determined the ratio of the photon energy density of O-and X-mode in the center of the layer. In fig-ure3.8the ratio of the photon energy density for two different photon energies as a function of the optical thickness is illustrated. It can be seen that the ratio is strongly dependent on the photon energy. In case of E = 0.1Ecycthe photon energy density is dominated by O-mode photons. However, one has to take into account that the layer of optical thickness τ = 100, due to non-magnetized scattering, is still opti-cally thin for photons of relatively low energy. From figure2.7in section2.3.1it can be seen that for relatively low photon energy the total cross sections for X- mode photons are about two orders of magnitude lower than the total cross sections for O-mode photons (across the field). As a result, X-mode photons leaves the system earlier than O-mode photons and hence the photon energy density is dominated by O-mode photons in case of E = 0.1Ecyc.

In case of E = 0.5Ecyc the photon energy density is dominated by X-mode pho-tons. For this particular energy the total cross sections for X-mode photons are still lower than total cross sections for O-mode photons, but not by two orders of mag-nitude. As a result, X-mode photons stay longer inside the layer. At the same time,