Variability of Bright X-ray
Pulsars
Report Bachelor Project Physics and Astronomy, size 15 EC,
conducted between 29-03-2016 and 07-07-2016
Author:
Patrick Verhagen
Daily supervisor: dr. A. Mushtukov Examiner: prof. dr. M.B.M. van der Klis Second reviewer: dr. A.L. Watts
Scientific abstract
In this thesis I used the theory of accretion onto highly magnetized neutron stars (NSs) in their bright states to create a model in Python which simulates the relation between the pulse fraction, defined as Fmax−Fmin
Fmax+Fmin, where Fmin and Fmax are the minimum and maximum observed flux within the pulse period at a given luminosity respectively, and the total accretion luminosity of a NS. Effects from special relativity are taken into account while general relativity is ignored. The predicted flux consists of reflected photons from the NS surface and of photons which are directly observed from the columns according to recent results of Poutanen et al. (2013). Theoretical pulse profiles were also constructed for different accretion column heights and therefore different luminosites. The predicted pulse fractions were then compared to observe pulse profile variations of X-ray pulsar V0332 + 53 during its outburst in 2004 - 2005, observed by the RXTE observatory (Lutovinov et al. 2015). The results from the model indicate that the pulse fraction firstly becomes smaller and then larger as luminosity and mass accretion rate increase. The simulated results are in good agreement with observed variations in the pulse fraction (Lutovinov et al. 2015; Tsygankov et al. 2010), but for high luminosities (L > 2.5 · 1038 erg · s−1) the model seems to underestimate the pulse fractions. This might be a consequence of excluding general relativity in the model. The model also predicts eclipsings of accretion columns by the NS itself which manifest themselves as very special features (steep drops) in the pulse profile. Since eclipses are caused by geometrical effects, their detection might be important for future measurements of NS radii.
Popular abstract (Dutch)
Neutronensterren zijn, op zwarte gaten na, de meest dichte waargenomen objecten in het heelal. Ze kunnen anderhalf keer zo zwaar zijn als de Zon en een straal van 70 000 keer kleiner hebben. Als een zware ster van minimaal ongeveer 10 keer de massa van de zon explodeert, kan er een ontstaan. Bij zo’n explosie stort de kern van de ster ineen en versmelten de protonen en elektronen met elkaar tot dicht verpakte neutronen en mogelijk andere exotische deeltjes. Doordat de massa van de ineenstortende kern nagenoeg hetzelfde blijft, maar het object veel kleiner wordt, zal de nieuw ontstane neutronenster snel gaan roteren.
Door een magnetische behoudswet in de natuurkunde zal het magneetveld zeer sterk worden. Bij massaoverdracht van een nabije ster (accretie) wordt materie door het magneetveld naar de magnetische polen van de neutronenster geleid. Dit gebeurt bij zulke hoge temperaturen dat de materie verandert in een plasma: een gas dat vooral bestaat uit ionen en vrije elektronen. Dit plasma kan met snelheden van een half keer de lichtsnelheid bewegen vlakbij de neutronenster. Bij deze hoge energieën kan er röntgenstraling ontstaan, die waargenomen kan worden vanaf de Aarde. Doordat de accretiematerie kan botsen op fotonen, die van de neutronenster af bewegen, kan de materie afgeremd worden en kunnen er zogenaamde accretiekolommen, soort cilinders met materie boven de magneetpolen, ontstaan.
In mijn model heb ik deze accretiekolommen gesimuleerd en bekeken hoe wij, hier op Aarde, het licht ervan kunnen waarnemen. Vanaf de Aarde wordt de röntgeninten-siteit als functie van de tijd waargenomen als een golfpatroon door het “vuurtoreneffect”: door de rotatie zullen de magnetische polen de ene keer richting de Aarde staan en de andere keer niet. Hierdoor varieert de lichtintensiteit als een cirkelbeweging en krijg je een periodiek verband.
Door verschillende modellen te bouwen met respectievelijk 1 puntvormige lichtbron boven de ster, dan 2 puntvormige lichtbronnen, dan 2 cilindervormige lichtbronnen en vervolgens nog de rotatie van de ster te implementeren, heb ik theoretische data kunnen creëren over hoe het eruit zou zien als je naar een neutronenster zou kijken met deze emissiegeometrie. De waargenomen fotonen kunnen in 2 delen gesplitst worden: een deel, dat gereflecteerd wordt vanaf het oppervlak van een neutronenster naar ons toe, en een deel dat direct vanaf de kolommen waargenomen wordt. Het doel is de relatie tussen de variatie in lichtintensiteit, als het verschil tussen de maximum- en minimumwaarden van het golfpatroon, en de hoogte van accretiekolommen te vinden. Voor lage kolommen (hoogte van de kolom gelijk aan ongeveer 0.2 keer de straal van de ster) lijkt de variatie kleiner te worden naarmate de kolommen hoger worden, doordat het licht gelijkmatiger wordt verdeeld over het oppervlak van de ster. Bij rotatie worden namelijk afwisselend verschillende delen van de neutronenster vanaf de Aarde waargenomen en zal het verschil tussen de maximale en minimale lichtintensiteit kleiner worden naarmate een groter deel van de neutronenster belicht wordt door de kolommen. Hierbij domineert dus het aandeel van gereflecteerde fotonen. Voor hoge kolommen (hoogte van de kolom gelijk
kolom. Een huidig probleem in de astronomie is dat het vaststellen van de straal van een neutronenster erg moeilijk gaat. Dit model kan aan de hand van verduisteringen van de accretiekolommen door de neutronenster zelf een schatting geven voor de straal door te kijken hoe sterk deze verduisteringen voorkomen als steile dalingen in het golfpatroon. Hierdoor zou het mysterie van de precieze inwendige structuur van de neutronenster misschien opgelost kunnen worden.
Popular abstract (English)
Neutron stars (NSs) are, except for black holes, the most compact objects in the Universe. They can be 1.5 times as heavy as the Sun and have a radius of 70 000 times smaller. If a star of minimally approximately 10 times the mass of the Sun explodes, a NS can arise. After such an explosion the core of star collapses and the protons and electrons melt together into compactly packed neutrons and possibly other exotic particles. Because the mass of the collapsing core basically stays the same, the new NS will rotate faster.
Figure 2: An artist’s impression of the accretion process near a rotating neutron star.
Because of a magnetic conservation law in physics, the new magnetic field will be very strong. During mass transfer from a nearby star to the NS (accretion) the matter is channeled to the magnetic poles of the NS by the magnetic field. This happens at such high temperatures that the matter turns into a plasma: a gas primarily consisting of ions and free electrons. This plasma can reach velocities of half the speed of light in the vicinity of the NS. At these high energies X-ray radiation can emerge, which can be observed from Earth. Because this accretion matter can collide with photons which move away from the NS, the matter can be decelerated and so-called accretion columns, kind of cylinders consisting of matter above the magnetic poles, can arise.
In my model I simulated these accretion columns and looked how we, here on Earth, can observe the light from them. From the Earth the X-ray intensity as a function of time is observed as a wave-pattern because of the “lighthouse-effect": because of the rotation, the magnetic poles will one time be directed towards the Earth and another time not. Be-cause of this the light intensity varies as a circular motion and a periodic relation will arise.
this emission geometry. The observed photons can be split in 2 parts: one part, that is being reflected from the surface of the NS towards us, and a part which can be directly observed from the columns. The goal is to find the relation between the variation in light intensity, as the difference between the maximum and minimum values of the wavepattern, and the height of the accretion columns. For low columns (with heights like 0.2 times the radius of the star) the variation seems to become smaller as the columns become higher, because the light is more homogeneously distributed over the NS surface. During rotation different parts of the NS are alternatingly observed from Earth and the difference between the maximum and minimum light intensity will become smaller as a larger part of the NS surface is illuminated. So for this the contribution of reflected photons dominates. For high columns (with heights of circa the radius of the star) the contribution of directly observed photons dominates, and therefore the variation in light intensity becomes larger as the columns become higher. This happens because a higher column can emit more photons than a lower column. A current problem within astronomy is that determining the radius of a NS is very hard. This model could give an estimation for the radius with the help of column eclipses by the neutron star itself by looking at how strong these eclipses appear as steep drops in the wave pattern. Because of this the mystery of the exact internal structure of the NS could maybe be solved.
1 Introduction 6
1.1 Motivation . . . 6
1.2 Theoretical background . . . 6
1.2.1 Magnetic field of Neutron Stars . . . 6
1.2.2 Introduction to XRPs . . . 7
1.2.3 General relativistic effects . . . 7
1.3 Approach . . . 11
1.3.1 General approach . . . 11
1.3.2 Physical steps taken in the model . . . 11
2 Results 19 2.1 Point sources . . . 19 2.2 Columns . . . 22 2.3 Pulse fractions . . . 25 2.4 Pulse profiles . . . 27 3 Discussion 28 3.1 1 point source & 2 point sources . . . 28
3.2 Comparison isotropic with anisotropic emission . . . 28
3.3 Pulse fractions . . . 29
3.4 Pulse profiles . . . 29
4 Conclusions 31
A Part of the code 32
Introduction
1.1
Motivation
Since the prediction of neutron stars (NSs) in 1934 by Baade and Zwicky (Baade & Zwicky 1934) new laboratories had been discovered where extreme physics could be tested. Not only do NSs test general relativistic effects with their strong graviational fields, but some of them might also have very strong magnetic fields in the order of 1012 G and even higher. The combination of both of these effects results in the accumulation of matter onto the stars, which can be observed from Earth in the form of light. Trying to understand what we observe is an essential part of science and therefore understanding how the accretion process works, can teach us more about these extreme beasts of nature.
1.2
Theoretical background
1.2.1 Magnetic field of Neutron Stars
In this thesis we are looking at NSs with strong magnetic fields in the order of ≥ 1012 G. A possible explanation for strong magnetic fields given by Spruit (2007) is the so-called fossil field hypothesis. The fossil field hypothesis primarily uses flux conservation to explain strong magnetic fields. Let us define the flux of the progenitor star Fprog:
Fprog=
Z
dSprog B⊥≈ 4πR2progB⊥ (1.1) Then the flux of the NS FNS will be equal to the flux of the progenitor star Fprog and the following relation can be determined:
Fprog= 4πR2progB⊥,prog= FNS= 4πR2NSB⊥,NS, (1.2) and therefore B⊥,NS B⊥,prog = R 2 prog R2NS (1.3)
If the radius of the progenitor star is taken as 1011 cm and the radius of the NS as 106 cm, then the magnetic field of the NS is 1010 times stronger than the magnetic field of the progenitor star. Since the strongest fields in main sequence stars are around 104 G (Spruit 2007), a magnetic field strength of 1014 G can be reached by the NS using flux conservation. In this sense NSs are the strongest magnets in the Universe.
1.2.2 Introduction to XRPs
X-ray pulsars (XRPs) are NSs which have strong magnetic fields (in the order of 1012 G) and accrete matter from a companion. The strong magnetic fields channel the accretion flow towards the magnetic poles of XRPs. A NS is covered by a geometrically thin, but optically thick atmosphere which mostly consists of free electrons and ions (Potekhin 2010).
If the mass accretion rate ˙M is low enough, then the matter reaches the NS surface
and converts its kinetic energy to photons. The photons are then emitted from parts of the NS surface which are called hot spots. This happens, since a low ˙M causes a low
radiation pressure, therefore the gravitational attraction is stronger and hot spots at the magnetic poles are created. For high ˙M the radiation pressure becomes so strong
that matter, during the accretion process, stops above the NS surface and accretion columns emerge (Basko & Sunyaev 1976; Mushtukov et al. 2015a). The height of an accretion column depends on the mass accretion rate and can be as high as the NS radius. If the mass accretion rate is higher, the luminosity is also higher (L = GM ˙M /R)
and this results in a higher radiation pressure. Since the height of accretion columns is dependent on the balance between the radiation pressure and graviational attraction, higher accretion columns are possible at higher mass accretion rates.
The theory on the structure of these accretion columns is still under debate (Lyubarskii & Syunyaev 1988; Mushtukov et al. 2015b). It is also expected that the flux distribution over the directions differs between XRPs with hot spots and XRPs with accretion columns (Gnedin & Sunyaev 1973). The theory on this is not robust enough yet to explain the differences between these 2 situations, but there is a relevant amount of data on the flux and spectral variability which could be helpful in the clarification of the theory (Tsygankov et al. 2006; Lutovinov et al. 2015).
1.2.3 General relativistic effects
Typical NS masses and radii are about 1.4 M and 10 - 15 km respectively (Haensel et al. 2007). For 1.4 M the Schwarzschild radius Rs = 2GMc2 ≈ 4.2 km, which means that the typical radius of a NS is 2-3 times the Schwarzschild radius. This means that effects of general relativity have to be taken into account. Also, since the dimensionless free-fall velocity at the surface of a NS is βff =p
Rs/R ≈ 0.5, effects of special relativity have to be considered.
According to Lightman (1975) the differential equation which describes the motion of photons around a gravitational point source is:
d2u
dφ2 + u = 3u
2, (1.4)
where u = Rs
2r and φ is the angle between photon momentum and the vector which is drawn from the NS center to the photon. The numerical solution to eq. (1.4) is shown in fig. 1.1 on the next page for α0 = π/2 and R = 2 · RS and in fig. 1.2 on the following page, where α0 is the angle between the normal to surface and the photon momentum during emission from the surface (see fig. 1.5 on page 10).
Figure 1.1: Numerical solution to eq. (1.4) on page 7 with φ0= 0, α0 = π/2 and R = 2 · RS,
where M = 1.4M .
Figure 1.2: Numerical solution to eq. (1.4) on page 7 with R = 106 cm, scaled to R = 1.0, and
M = 1.4M , which results in an RS ≈ 0.42. α0is the angle between the normal to the surface
and the photon momentum during the moment of emission. The green circle represents the NS surface, while the dashed lines represent the photon trajectories for different α0angles.
In case of flat space-time u = RS
2r << 1, the differential equation could be approxi-mated as ddφ2u2 + u = 0. The numerical solution to this case is displayed in fig. 1.3 on the following page.
Figure 1.3: Numerical solution to ddφ2u2 + u = 0 with φ0= 0, α0= π/2 and the radius of the NS
equal to R = 2 · RS, with RS = 2GM/c2, where M = 1.4M .
Figure 1.4: Numerical solution to d2u
dφ2 + u = 0 with R = 10
6 cm, scaled to R = 1.0, and
M = 1.4M , which results in an RS ≈ 0.42. α0is the angle between the normal to the surface
and the photon momentum during the moment of emission. The green circle represents the NS surface, while the dashed lines represent the photon trajectories for different α0angles.
For numerical solutions, this second order ODE was split into a system of two first order ODEs by taking v = dφdu:
du dφ = v dv dφ = 3u 2− u .
If the initial distance between the photon and the NS center is larger than 2RS, then the following approximation can be used (Poutanen et al. 2013):
cos α = RS r + 1 −RS r cos ψ, (1.5)
where α is the angle between the photon momentum and the normal to the surface,
RS = 2GMc2 is the Schwarzschild radius, r is the distance between the photon and the NS center and ψ is the angle between the trajectory into infinity and the normal to the surface at the starting point of the trajectory (see fig. 1.5).
Figure 1.5: Illustration of the situation, in which a photon is emitted from the surface travelling through curved spacetime (Poutanen et al. 2013).
When the photon is emitted from the NS surface and α0= π/2, we will have α = α0 and r = R. Equation (1.5) can be rewritten as the following equation and used to find ψ for given NS radius R with mass M = 1.4M :
cos α0 = RS R + 1 −RS R cos ψ (1.6) RS 2R = R S R − 1 cos ψ (1.7) ψ = sec−1 1 − R RS . (1.8)
The light bending has a number of consequences. One of them is that a distant observer can always see more than half of the NS surface.
The observable fraction of the surface Ω = 2π
ψ R
0 dθ sin θ
4π depends on the NS compact-nessRR
S
. In fig. 1.6 on the next page the observable fraction Ω is plotted as a function of the NS compactnessRR
S
Figure 1.6: The observable fraction of the NS surface Ω as a function of R/RS, where R is
the NS radius and RS is the Schwarzschild radius 2GM/c2, using eq. (1.5) on page 10. The
approximation holds for R/RS > 2.0. The horizontal dashed lines are Ω = 1 and Ω = 0.5. In
the limit of very high R/RS, Ω goes to 0.5, which is in agreement with flat space-time, since the
mass is fixed.
1.3
Approach
1.3.1 General approach
In this thesis I would like to predict certain physical properties of highly magnetized accreting NSs: namely the PF dependence on L. A bright XRP will be considered as a sphere illuminated by 2 accretion columns. The height of an accretion column depends on the mass accretion rate ˙M . A model is created in Python. The model simulates the
PF as a function of L and simulates pulse profiles for certain conditions. Effects from special relativity are taken into account, but the effects of general relativity are ignored. The final results will be in the form of a theoretical pulse profile and the pulse fraction as a function of luminosity of the source. The theoretical pulse fractions are compared to the observed pulse fractions (Lutovinov et al. 2015) in order to test the validity of the model. The used observational results give pulse fractions in the energy range of 25-45 keV versus the bolometric luminosity during an outburst of XRP V0332 + 53 in 2004 -2005.
1.3.2 Physical steps taken in the model
First a model was set up where the accretion columns were represented by point sources located on the magnetic axis which isotropically emit photons at height h above the NS surface. The photons collide with the NS surface. The consequent flux Fin incident on the NS surface at latitude θB, which is the angle between the magnetic axis and the normal to the surface where the photon reaches the NS surface:
Fin(θB) = L ps 4πD(θB, h)2 cos αin(θB, h), (1.9)
where Lps is the luminosity of a point source and D(θB, h) = q h2+ 2R(R + h)(1 − cos θ B), cos αin= h cos θB− R(1 − cos θB) D , (1.10) where D(θB, h) is the distance from the point source to a point on the NS surface and
αin describes the angle between the initial photon momentum and the normal to the surface at latitude θB. The angles for the magnetic reference frame are illustrated in fig. 1.7.
Figure 1.7: Illustration of the NS with the angles from the magnetic and the observer’s reference frame.
In order to calculate the observed flux, one can move from the magnetic reference frame to the observer reference frame, where θ0 and φ0 are the spherical angles, as illustrated in fig. 1.7, and α0 is the angle between the normal to the surface and the momentum of the reflected photon.
Let ~r0 be the unit vector from the center of the NS to the observer and let ξ be the angle between the magnetic axis and the observer’s line of sight (see fig. 1.7). Then
ˆ
mx is the rotation matrix which rotates a vector with angle ξ around the x-axis:
~ r0 = r sin θ0cos φ0 r sin θ0 sin φ0 r cos θ0 , mˆx = 1 0 0 0 cos ξ − sin ξ 0 sin ξ cos ξ . (1.11)
Using eq. (1.11) on page 12 the coordinates of vector ~r0 can be calculated in the magnetic reference frame:
~ rB= ˆmxr~0 = r sin θ0cos φ0
r sin θ0sin φ0cos ξ + r cos θ0sin ξ
r sin θ0sin φ0sin ξ + r cos θ0cos ξ
= r sin θBcos φB r sin θBsin φB r cos θB . (1.12)
If energy conservation is considered the assumption can be made that the total flux re-emitted by a unit area of the NS equals to the total absorbed flux by the NS surface. Therefore, the flux which is absorbed by the NS surface at latitude θB satisfies the following relation:
Fin(θB) = 2π π/2
Z
0
dα0 Iout(θB, α0) cos α0sin α0, (1.13)
where Iout is the intensity of the outgoing radiation. If the intensity of the outgoing radiation is isotropically distributed, i.e. it is independent on α0, eq. (1.13) can then be rewritten as:
Fin(θB) = πIout(θB) (1.14) One could also assume that the reflection of photons happens anisotropically. If the NS atmosphere is optically thick, then photons which travel along the normal to the NS surface have a higher probability than other reflected photons to exit the atmosphere (Mihalas 1978). According to Poutanen et al. (2013) the outwards intensity Iout can be
expressed as follows: Iout(θB, α0) = 1 +π 2 cos α0 3F in(θB) π(π + 3). (1.15)
The flux which the observer sees, is radiated from half of the sphere in case of flat geometry. This flux can be calculated by integration over the visible part of a NS:
Fobs,sur = 2π Z 0 dϕ0 π/2 Z 0
dθ0sin θ0cos θ0I(θ0, ϕ0, α0), (1.16)
where θ0 and φ0 are the polar angles in the observer’s reference frame (see fig. 1.7 on page 12). If the observer is very far away (r >> R), then α0 ≈ θ0. This results in the following 2 equations using the z-component of eq. (1.12) :
I(θ0, ϕ0, α0) = Iout(θB, α0), cos θB = sin ξ sin ϕ0sin θ0+ cos ξ cos θ0. (1.17) So far, the assumption has been made that the point sources emit radiation isotrop-ically. In reality, the emitted photons are subject to relativistic Compton scattering, because an accretion column is surrounded by a layer of fast moving plasma with ve-locities as high as 0.5c (Lyubarskii & Syunyaev 1988). As a result, according to special relativity, the photons from the accretion column are beamed towards the NS surface and its distribution over the directions is far from the isotropical case. The anisotropy can be described by f (φ) which describes the dependence of the outgoing flux on angle
φ, defined as the angle between the magnetic axis and the photon momentum. In order
to use a normalized function f (φ) over all angles, the following requirement is made:
2π π
Z
0
dφ sin φf (φ) = 1. (1.18)
The flux distribution over the NS surface can be recalculated then:
Fin(θB) = L psf (φ) 4πD(θB, h)2 cos α(θB, h). (1.19)
Let dL(φ)dφ be the angular distribution of the luminosity, then this function is normalized in the following way:
π
Z
0
dφdL(φ)
dφ = Ltot, (1.20)
where Ltot is the total luminosity of 1 accretion column.
Using eq. (1.18) and eq. (1.20), eq. (1.19) can be rewritten as follows:
Fin(θB) = dL(φ)/dφ 2π sin φD(θB, h)2 cos α(θB, h), (1.21) where φ = α − θB.
Before integrating over the surface, it is useful to introduce a dimensionless angular luminosity distributiondL(φ)dφ ∗: dL(φ) dφ ∗ = dL(φ)/dφ Ltot , π Z 0 dφ dL(φ) dφ ∗ = 1 (1.22)
According to Poutanen et al. (2013) dLdφ can be well described as:
dL dφ = I0sin 2φ γ5(1 − β cos φ)4 1 +π 2 sin φ γ(1 − β cos φ) , (1.23)
where I0 is the normalization constant, β is the dimensionless plasma velocity in the layer of free-falling plasma and γ is the corresponding Lorentz factor. One can see that the luminosity distribution over the directions depends on the velocity of the plasma, as illustrated in fig. 1.8 on the following page, and on the projected solid angle which is proportional to sin φ. So the intensity stays constant over the different
φ angles, but the flux does not, because it is defined by the solid angle where the
intensity is concentrated. For high velocities the photons are beamed towards the NS surface and for low velocities the radiation distribution over the directions is mostly affected by the geometrical projection of the emitting area for a given direction. If the dimensionless velocity β = v/c = 0, the luminosity distribution takes the form of a symmetric distribution around φ = π/2. This corresponds to the classical result of the theory on stellar atmospheres (Mihalas 1978).
Figure 1.8: The angular distribution of the luminosity dLdφ for different dimensionless electron free-falling velocities β = v/c. At high β the radiation is beamed strongly towards the NS surface, which is located in the downwards direction, and at low β the radiation is more symmetric around
φ = π/2.
One can perform normalization with eq. (1.22) on page 14, and then the dimensionless angular luminosity distributiondL(φ)dφ ∗ can be calculated:
π
Z
0 dφI0
sin2φ1 +π2γ(1−β cos φ)sin φ
γ5(1 − β cos φ)4 = 7πI0 6 (1.24) =⇒ dL(φ) dφ ∗ =
6 sin2φ1 +π2γ(1−β cos φ)sin φ
7πγ5(1 − β cos φ)4 . (1.25) In case of extended accretion columns the luminosity might be nonhomogeneously distributed over the column height. Let g(h) be the function which describes the distribution of the luminosity over the height of the accretion column with the following normalization: g(h) = dL(h) dh , H Z 0 dh g(h) = Ltot. (1.26) The column is approximated as an infinite amount of infinitesimally thick cylinder slabs with luminosity g(h)dh consisting of point sources. The flux distribution over the NS surface in case of accretion columns can be obtained by integration over a part of the emitting point sources in the accretion column according to the following equation:
Fin(θB) = H Z hmin dh dL(φ(θB,h)) dφ ∗ g(h) 2πD(θB, h)2 sin φ(θ B, h) cos αin(θB, h), (1.27)
where hmin is the minimum height of the accretion column which contributes to the flux at latitude θB:
hmin =
R
The observed flux consists of two parts: the flux Fobs,col, which is directly observed by the observer, without intervention of the NS surface, is calculated by integrating the angular luminosity distribution times the height distribution over the part of the accretion column which can be seen by the observer:
Fobs,col(ξ) = 1 2π sin φ0 hmax Z hmin dh dL(φ0) dφ ∗ g(h), (1.29)
where φ0 = ξ or φ0= π − ξ and ξ is the angle between the observer’s line of sight and the magnetic axis. φ0 is chosen depending on which accretion column is being considered. It is clear that Fobs,colis completely defined by the angle ξ for a system of given parameters. In order to construct a rotation matrix ˆm∗ which transforms coordinates from the magnetic frame to the observer’s frame, one has to take the inverse of the matrix ˆmx (see eq. (1.11) on page 12) which describes a rotation of angle −ξ around the x-axis.
ˆ m∗ = ˆm−1x = 1 0 0 0 cos ξ sin ξ 0 − sin ξ cos ξ (1.30) Then: ~ r0 = ˆm∗r~B = r sin θBcos φB
r sin θBsin φBcos ξ + r cos θBcos ξ − sin φBsin θBsin ξ + cos θBcos ξ
= r sin θ0cos φ0 r sin θ0sin φ0 cos θ0 (1.31)
The second part of the observed flux which consists of the radiation reprocessed by the NS atmosphere is defined by the flux distribution over the NS surface and the orientation of a NS in the observer’s reference frame:
Fobs,in(ξ) = 2π Z 0 dφB π Z 0 dθBsin θBFin(θB) π × max {0, cos θ0} , (1.32)
where cos θ0= − sin φBsin ξ + cos θBcos ξ. The maximum of the set {0, cos θ0} is taken in order to prevent integration over the part of the NS surface which is invisible to a distant observer. The integration was performed in the magnetic reference frame in order to avoid numerical problems caused by possibly high incident fluxes Fin(θB) in the vicinity of the magnetic poles.
According to Poutanen et al. (2013) the reprocessed flux dominates the total observed flux in the case of accretion columns with low H and the direct flux dominates in the case of accretion columns with high H.
According to Mushtukov et al. (2015b) the total luminosity of an accretion column can be approximated as:
L ≈ 38 l 0/d0 50 κ T κ⊥ f H R LEdd, f H R ≡ log 1 +H R − H R + H, (1.33)
where l0 and d0 define the base geometry of the accretion column, κT and κ⊥ are the Thomson scattering opacity and magnetic Compton scattering opacity respectively,
Figure 1.9: A visual representation of the used model for the optical thickness of an accretion column. h is the height of the accretion column, d0 is the geometrical thickness and l0 is the
geometrical length/circumference of an accretion column (Mushtukov et al. 2015b).
and LEdd is the Eddington luminosity for a given NS mass. The assumption that l0/do
50 ≈ κT
κ⊥ ≈ 1 is made in the initial model. After fitting the model on observational data, these parameters can be determined with certain errors. The attentive reader can see that the luminosity of such accretion columns can exceed 40 LEdd and these luminosi-ties belong to so-called ultra-luminous X-ray sources (Feng & Soria 2011). The height where the matter is stopped, varies inside the accretion column. Because the radiation pressure is lower at the boundary of the column than at the center of the column, the mat-ter stops higher at the cenmat-ter of the column (Lyubarskii & Syunyaev 1988). As a result, the photons have to penetrate a layer of free-falling particles. The scattering by these particles causes the special relativistic beaming effect as described in eq. (1.24) on page 15. The geometrical thickness of a matter sinking region is (Mushtukov et al. 2015b):
x d/2 = 1 − h H R + H R + h 1/2 . (1.34)
Then the optical thickness of the free-fall region τff can be estimated as follows:
τff = 1.7 × 105 l0 L37 βff 1 − x d/2 , (1.35)
where βff is the dimensionless freefall velocity and the dimensionless luminosity L37=
L/1037 erg · s−1. Looking at eq. (1.35) one can see that τff > 1 for the expected set of
parameters. It means that most of the photons from a sinking region will be scattered in the free-falling region and, therefore beamed towards the NS surface.
In order to implement the optical thickness τff in the code, the height dependence function g(h) is split in two functions:
where g1(τff, h) corresponds to photons which penetrate the free-falling layer without
scattering there and have a momentum perpendicular to the column surface and g2(τff, h)
corresponds to photons which are scattered towards the NS surface. To construct both of these functions, the normalization
H
R
0
dh g(h) = 1 has to be taken into account:
g1(τff, h) = g(h)e−τff, g2(τff, h) = g(h)(1 − e−τff) (1.37) However, it was already mentioned that radiation from the column is dominated by beamed photons, therefore g2(τff, h) ≈ g(h).
In order to construct the pulse profile three angles have to be defined, because the orientation of the rotational axis now has a significance: η is the angle between the rotational axis and the magnetic axis, µ is the angle between the rotational axis and the observer line of sight and δ is the angle between the projection of the rotational axis on the magnetic axis and the line between and the rotational axis and observer as illustrated in fig. 1.10.
Figure 1.10: Sketch with η as the angle between the rotational axis and the magnetic axis, µ as the angle between the rotational axis and the observer line of sight and δ as the angle between the projection of the rotational axis on the magnetic axis and the line between the rotational axis and the observer, also called the phase-angle.
Take ~B as the unit vector in the direction of the magnetic axis from the center of the
NS and ~O as the unit vector in the direction of the observer from the center of the NS.
These vectors have the following coordinate representations:
~ O = sin µ cos µ 0 , ~ B = sin η cos δ cos η sin η sin δ . (1.38)
ξ is still defined as the angle between the observer’s line of sight and the magnetic axis.
Therefore one can get ξ as a function of µ, η and δ:
Results
For all of the following numerical results the flux distribution over the height of the accretion columns was chosen as a constant:
H
Z
0
dh g(h) = 1 ⇒ g(h) = 1
H, (2.1)
however one can choose any g(h) in the model.
2.1
Point sources
We started with 1 point source isotropically emitting photons. The first result of the model was the flux absorbed by the NS Finas a function of the latitude at the NS surface
θB. To calculate this eq. (1.9) on page 11 was used. In fig. 2.1 on the next page the point source was placed at a distance h = R and at a distance h = 100R from the NS surface.
For 2 point sources, the total flux absorbed by the NS surface changes. Using ro-tational symmetry of π radians, the total flux absorbed as a function of θB becomes:
Fin,2ps(θB) = Fin(θB) + Fin(π − θB). This resulted in the produced data by the model for h = R and h = 100R as seen in fig. 2.2 on the following page.
Figure 2.1: The flux distribution over the NS surface if 1 point source at height above the NS surface along the magnetic axis h = R and h = 100R is considered. The flux for h = 100R is multiplied by 104. The higher the point source above the NS surface, the larger the illuminated
area. The flux from the point source at h = 100R becomes 0 at θB≈ π/2, while the flux from
the point source at h = R is 0 at approximately θB ≈ 1.04. So the point source at h = 100R
is able to provide flux for halve of the sphere, while the point source at h = R is eclipsed for latitudes θB between 1.04 and π/2.
Figure 2.2: The flux distribution over the NS surface if 2 point sources at h = R and h = 100R are considered. The flux for h = 100R is multiplied with 104.
Using eq. (1.21) on page 14 and eq. (1.23) on page 14 we can take effects of relativistic photon-electron scattering into account and this results in emission from the point sources beamed towards the NS surface. In fig. 2.3 the inwards flux as a function of latitude θB is calculated with the radiation beaming.
Figure 2.3: The flux distribution over the NS surface in case of 2 anisotropically emitting point sources at h = R and h = 100R are considered. The flux for h = 100R is multiplied with a factor of 2.9316 · 106. As a consequence of the anisotropic emission the magnetic poles themselves do not absorb any flux by the point sources. This is caused by the projection of the solid angle as seen in fig. 1.8 on page 15.
2.2
Columns
Since accretion columns are more accurately approximated as cylinders than as point sources, the point sources are converted to radiating cylinders with infinitely small radii and a total height of H. In fig. 2.4 column heights of H = 0.2R and H = R are taken. The incident flux as a function of angle θB is then calculated using eq. (1.27) on page 15.
Figure 2.4: The flux distribution over the NS surface in case of 2 anisotropically emitting accretion cylinders with height H = 0.2R and height H = R are considered. In contrary to the 2 anisotropically emitting point sources, flux is now absorbed maximally at the magnetic poles.
The flux from the NS surface which is detected by the observer depends on the flux distribution over the NS surface and on the orientation of the NS in the observer’s reference frame. The NS orientation is defined by the angle ξ (see fig. 1.7 on page 12). The column heights H = 0.1R, H = 0.5R and H = R were taken.
Figure 2.5: The observed flux, which is caused by re-emitted radiation, as a function of the angle between the magnetic reference frame and observer’s line of sight ξ (see fig. 1.7 on page 12). The used accretion column heights are H = 0.1R, H = 0.5R and H = R.
In order to calculate the flux which is directly observed from the accretion columns eq. (1.29) on page 16 was used. Column 1 is chosen as the column at θB= 0 and column 2 is chosen as the column at θB = π. In fig. 2.6 the directly observed flux is shown as a function of the angle between the magnetic axis and the observer’s line of sight ξ for
H = R and in fig. 2.7 on the next page for H = 100R. In fig. 2.8 on the following page
the directly observed flux for 2 columns is plotted for H = 0.1R, H = 0.5R and H = R.
Figure 2.6: Distribution of the directly observed flux from the accretion columns over the directions ξ for H = R. The blue dashed line represents column 1, the green dash-dotted line represents the flux of column 2, and the red solid line represents the sum of the fluxes. It is interesting that for most of the angles ξ the directly observed flux is dominated by photons from the column on the opposite side of the NS in flat space-time, because of special relativistic radiation beaming.
Figure 2.7: Directly from the column with H = 100R observed flux as a function of angle ξ. The blue dashed line represents column 1, the green dash-dotted line represents the flux of column 2, and the red solid line represents the sum of the fluxes. For physically unreasonable high accretion columns, the directly observed flux distributes in a "fan-beam" form, which was assumed to be the general case by Gnedin & Sunyaev (1973).
Figure 2.8: Direct flux from the accretion columns with different heights as a function of angle ξ.
The sum of the reflected flux (see fig. 2.5 on page 23) and the directly observed flux from the columns (see fig. 2.8) gives the total observed flux. This flux is presented in fig. 2.9 on the next page.
Figure 2.9: Total observed flux as a function the angle between the observer’s line of sight and magnetic axis ξ for different accretion column heights. The total observed flux consists of the flux which is absorbed and re-emitted by the NS, and of the flux which is directly observed from the accretion columns.
2.3
Pulse fractions
In order to simulate the pulse fraction, defined as P F (H, ξmin, ξmax) = Fmax−FFmax+Fminmin, where
Fmaxis the maximum flux in the interval [ξmin, ξmax] and Fmin is the minimum flux in the interval [ξmin, ξmax]. The total flux from the NS with accretion columns above the surface was calculated. Considering the effects of optical thickness g(h) from eq. (1.37) on page 18 was used. Since the model simulates the pulse fraction as a function of the accretion column height, eq. (1.33) on page 16 was used to find the luminosities associated with the column heights. The pulse fractions were calculated for accretion column heights between H = 0 and H = 1.5R. These data were then plotted against the observational data by Lutovinov et al. (2015) in fig. 2.10 on the following page. In order to find the luminosity out of the accretion colum height eq. (1.33) on page 16 was used with l0d050 κTκ⊥ = 0.4 for the pulse fractions in fig. 2.10 on the following page.
Figure 2.10: The pulse fraction as a function of luminosity simulated by the model with an isotropic reflection from the NS surface and observational data by Lutovinov et al. (2015) are being shown in the figure. For the data produced by the model l0d0
50
κT
κ⊥ = 0.4 was taken. Until now an isotropic reflection was considered in the model. By using eq. (1.15) on page 13 anisotropic reflection was assumed and other pulse fractions were calculated and fitted on the observational data using the increase in pulse fraction around L = 2 · 1038erg · s−1. This resulted in taking l050d0κTκ⊥ = 0.25. The results are shown in fig. 2.11.
Figure 2.11: The pulse fraction as a function of luminosity simulated by the model with an anisotropic reflection from the NS surface and observational data by Lutovinov et al. (2015) are being shown in the figure. For the data produced by the model l0d0
50
κT
2.4
Pulse profiles
Several pulse profiles, which are defined as the total observed flux as a function of the phase-angle δ, were simulated using eq. (1.39) on page 18 and for given parameter µ and
η. In fig. 1.10 on page 18 this situation is illustrated. First, the angles η and µ were
taken both as π/2 for fig. 2.12 for the accretion column heights H = {0.1R, 0.5R, R}.
Figure 2.12: Flux as a function of the phase-angle δ. For the angle between the rotational axis and magnetic axis, η, and the angle between the rotational axis and observer line of sight, µ, the value π
2 was chosen. The situation is illustrated in fig. 1.10 on page 18. In this figure eclipsings
can be seen as steep drops in the pulse profile at δ = {0, π, 2π}.
In fig. 2.13 the pulse profiles for the same accretion column heights as in fig. 2.12 were simulated for η = 3π4 and µ = π4.
Figure 2.13: Flux as a function of the phase-angle δ. For the angle between the rotational axis and magnetic axis, η, 3π4 was chosen and for the angle between the rotational axis and observer line of sight, µ, the value π4 was chosen. The situation is illustrated in fig. 1.10 on page 18.
Discussion
3.1
1 point source & 2 point sources
The results for the inwards flux from 1 point source (see fig. 2.1 on page 20) seem physically reasonable. If the point source is located very high above the NS surface (h >> R), then the component of the momentum perpendicular to the surface, is 0 at
θB= π/2 and is maximum at θB= 0, since at that angle all of the photon momentum contributes towards the flux. If the point source is closer to the NS surface, for example
h = R, then a fraction of the emitted photons is eclipsed by the NS surface. This would
result in a smaller illuminated area on the NS surface. For example it can clearly be seen in fig. 2.1 on page 20 that the NS surface at θB ≥ π/3 is not illuminated by the point source at h = R.
Since the symmetry argument can be made that the incident flux as a function of the latitude θB for 2 point sources Fin,2ps(θB) = Fin(θB) + Fin(π − θB), it should be no surprise that with 2 point sources, the diagram can be reflected in θB = π/2 and still have the same features.
3.2
Comparison isotropic with anisotropic emission
When anisotropic emission is considered (see eq. (1.21) on page 14) the distribution of flux over the surface changes significantly. The flux is at a minimum of 0 at θB= {0, π} instead of at a maximum if fig. 2.2 on page 20 and fig. 2.3 on page 21 are compared to each other. In order to understand this, one has to look at eq. (1.23) on page 14: if
θB= 0, then φ = α. Using the definitions of α and D described in eq. (1.10) on page 12: cos α = h
D(θB, h) = h
h = 1
α = arccos 1 = 2κπ : κ N ∪ {0}
Since φ = α and eq. (1.21) on page 14 ∝ sin φ, the inwards flux equals 0 at θB = 0. Because the situation is rotationally symmetric at π radians, the same argument can be made for θB = π.
However, this results concerns the situation of point sources above the NS surface. In case of extended accretion columns the flux distribution over the NS surface has its maximum at the poles if g(h) = constant.
3.3
Pulse fractions
According to the model presented here the total observed flux consists of the direct flux from the accretion column and the flux re-emitted from the NS surface. For relatively low luminosities and therefore low accretion column heights the total flux is dominated by re-emitted photons. If this would also be the case for relatively high luminosities, the pulse fraction would only decrease with higher luminosity, since the area of illumination on the NS surface would also be larger, see fig. 2.5 on page 23. The larger illuminated area would result in a more homogoneous flux distribution over the whole NS surface. However, at relatively high mass accretion rates the total observed flux is dominated by the direct flux from the accretion column if fig. 2.5 on page 23 and fig. 2.8 on page 24 are compared to eachother. In this case the pulse fraction increases with higher luminosity, since the flux can take on more varied values. An interesting observation is that for high mass accretion rates the flux is actually dominated by photons from the column on the opposite side of the star from the observer’s perspective, due to radiation beaming, see fig. 2.6 on page 23 and fig. 2.7 on page 24.
A combination of both of these modes of behavior results in a model which provides results in good agreement with observational data (Lutovinov et al. 2015) as seen in fig. 2.10 on page 26 and fig. 2.11 on page 26. This result shows a pencil beam diagram for the distribution of flux over the directions rather than the fan beam diagram earlier assumed by Gnedin & Sunyaev (1973).
The calculated pulse fractions which correspond to isotropic reflection are smaller in comparison with Lutovinov et al. (2015) for low luminosities (1037 erg · s−1 < L <
1.0 · 1038 erg · s−1) and high luminosities (L > 2.5 · 1038erg · s−1). The results calculated for anisotropic reflection seem to be in better agreement with the observational results for the low and high luminosities than in the isotropic reflection case. However, they also do not fit with all of the observations within the errors.
Another sidenote to comparing the simulated pulse fractions with the observed pulse fractions is that the limited resolution in time from observations could result in an underestimation of the pulse fraction. The fact is that the maximum observed flux could be underestimated and the minimum flux overestimated, which would result in an underestimated pulse fraction.
It is uncertain how GR effects will exactly influence the results, but since the ob-servable part of the NS surface is > 1/2 for typical NS compactnesses (R/RS), as a result of light bending (see fig. 1.6 on page 11), the observer will see more reflected flux.
3.4
Pulse profiles
The visibility of eclipsing in the pulse profile is expected to be dependent on the luminos-ity. In fig. 2.12 on page 27 one can see eclipses as steep drops at θB = {0, π, 2π} for the column height H = R. The luminosity should be sufficiently high, because only in this case the observed flux would be dominated by the direct flux contribution. In fig. 2.12 on page 27 one can also see that for H = 0.1R the eclipses are less severe than for H = R,
which confirms the previous statement. At the same time it is easier to eclipse small accretion columns, which correspond to relatively small luminosities. This can be seen in fig. 2.12 on page 27 as longer lasting eclipses in the pulse profile for H = 0.1R than for
H = R. As a result, some luminosity range is expected where eclipses of the accretion
columns by the NS surface are observable.
One would expect steep drops of the flux at its maximum during the pulse period. Such drops were already observed in pulse profiles of XRP V0332+53 at the highest luminosity during the outburst (Lutovinov et al. 2015). It is important to note that the orientation of a NS in the observer’s reference frame is not the only factor which determines the accretion column eclipsing. The eclipsing is also determined by the NS radius. This is the case, because a larger NS would more likely eclipse its accretion columns than a smaller NS. Since the height of accretion columns is expected to be compatible with the NS radius in the brightest state of a given XRP (Poutanen et al. 2013), one could set constraints on the mass-radius relation of NSs at high mass accretion rates if an accurate model for the accretion column is available. It is principally possible to distinguish between the direct flux from the column and the reflected flux (Postnov et al. 2015). Since it is expected that during eclipsing a part of the total observed flux is dominated by the reflected signal instead of the directly observed one, it would be possible to verify the hypothesis of eclipsing by phase-resolved spectroscopy.
Conclusions
In order to explain the relation between the pulse fraction and the luminosity of bright X-ray pulsars a model was constructed based on the idea of accretion columns above the NS surface.
According to the constructed model the total observed flux can be split in two dif-ferent contributions: the flux directly observed from the accretion columns and the reflected photons from the NS surface. For relatively low luminosities and therefore small accretion column heights, the total observed flux mostly consists of the reflected photons, while for high luminosities, the total observed flux is dominated by the direct flux from the accretion columns. As expected and seen in fig. 2.8 on page 24 the directly observed flux from the accretion columns is higher for higher accretion columns and higher luminosities. It is interesting that the direct flux is dominated by photons from the accretion column on the opposide side of the star with respect to the observer (see fig. 2.6 on page 23). This effect is caused by special relativistic radiation beaming. For the direct contribution towards the total flux the pulse fraction increases with luminosity and accretion column height.
For the contribution of re-emitted photons from the NS surface, it looks like the pulse fraction decreases with luminosity, see fig. 2.5 on page 23. This is the case, because if the luminosity is higher, the accretion column is higher and this results in a larger illuminated area on the NS surface. Combining both of these behaviours in fig. 2.10 on page 26 and fig. 2.11 on page 26 ensues a pulse fraction - luminosity relation, which is in good agreement with recent observational results (Lutovinov et al. 2015).
By simulating a theoretical pulse profile in fig. 2.12 on page 27 and fig. 2.13 on page 27 eclipsing could be observed as unexpected steep drops in the pulse profile. It is unknown how GR effects will exactly change the results, but it is certain that the reflected contri-bution to the total observed flux will change, because then more than half of the NS surface is observable for typical NS compactnesses (see fig. 1.6 on page 11).
Part of the code
from __future__ import division # Python 2 compatibility import math as m
import numpy as np
import matplotlib.pyplot as plt import scipy.integrate as integrate import matplotlib
'''
The angle zeta in the code corresponds with angle xi in the thesis. '''
# Initial conditions using CGS L = 1.0 # Convenience
R = 1.0 # Convenience c = 3.0*10**10
G = 6.674 * 10**(-8)
M = 1.4 * 1.98855 * 10**(33) # 1.4 solar masses in g def ksi_to_delta(ksi,mu,eta):
'''
Calculates the phase-angle for a given orientation '''
return m.acos( (m.fabs(m.cos(ksi))- m.fabs(m.cos(mu)*m.cos(eta)) )
/ ( m.sin(mu)*m.sin(eta) ) )
,→
def ksi_to_delta_original(ksi,mu,eta):
return m.acos( (m.cos(ksi)- m.cos(mu)*m.cos(eta) ) / ( m.sin(mu)*m.sin(eta) ))
,→
def distance_travelled(h,R,theta_b): '''
Calculates the travelled distance for reflected photons between the point of emission and NS surface
,→
'''
def perp_factor(theta_b,h,R): '''
Calculates the perpendicular (to the NS surface) component for the photon momentum.
,→
'''
return (h*m.cos(theta_b) - R*(1.0-m.cos(theta_b))) / distance_travelled(h,R,theta_b)
,→
def hmin_theta_b(theta_b,R): '''
Calculates the minimum accretion column height from which the column is observable for a given latitude.
,→
'''
return R*(m.cos(theta_b)**(-1) - 1) def hmin(zeta,R):
'''
Calculates the minimum accretion column height from which the column is observable for a given orientation.
,→
'''
return R*( (1.0/m.sin(zeta)) - 1.0) def phi(theta_b,h,R):
'''
Calculates the photon momentum angle for a given latitude. '''
phi = m.acos(perp_factor(theta_b,h,R)) - theta_b # alpha - theta_b return phi
def beta(M,h): '''
Calculates the free-fall velocity of a particle towards the NS surface. ,→ ''' R_s = 0.42 beta = ( R_s/(R+h))**(0.5) return beta def gamma(M,h): gamma = (1.0-beta(M,h)**2.0)**(-0.5) return gamma
def luminosity_change_star(theta_b,h,R,M): '''
dl/d(phi)* as a function of phi(theta_b,h,R) '''
function = (7*m.pi)**(-1) * 6.0 * I_0 *
m.sin(phi(theta_b,h,R))**(2.0) * (1.0+ (0.5*m.pi *
m.sin(phi(theta_b,h,R))/(gamma(M,h)*(1.0-beta(M,h)*m.cos(phi(thc eta_b,h,R)))))) / ( gamma(M,h)**(5.0) * ( 1.0-(beta(M,h)*m.cos(phi(theta_b,h,R)) ))**4.0 ) ,→ ,→ ,→ ,→ return function
def luminosity_change_star_beta_0(theta_b,h,R): '''
dl/d(phi)* as a function of phi(theta_b,h,R) with beta = 0 '''
I_0 = 1.0
function = (7*m.pi)**(-1) * 6.0 * I_0 *
m.sin(phi(theta_b,h,R))**(2.0) * (1.0 + (0.5 * m.pi *
m.sin(phi(theta_b,h,R) ) ) )
,→ ,→
return function
def luminosity_change_star_not_normalized_phi(phi,M,h): I_0 = 1.0
function = I_0 * m.sin(phi)**(2.0) * (1.0 + (0.5*m.pi *
m.sin(phi)/(gamma(M,h)*(1.0-beta(M,h)*m.cos(phi))))) / ( gamma(M,h)**(5.0) * (1.0- (beta(M,h)*m.cos(phi) ))**4.0 )
,→ ,→
return function
def luminosity_change_star_not_normalized_phi_beta(phi,beta): I_0 = 1.0
gamma = (1.0-beta**(2.0))**(-0.5)
function = I_0 * m.sin(phi)**(2.0) * (1.0 + (0.5*m.pi *
m.sin(phi)/(gamma*(1.0-beta*m.cos(phi))))) / ( gamma**(5.0) *
(1.0- (beta*m.cos(phi) ))**4.0 )
,→ ,→
return function
def plot_dldphi_phi_beta(beta1,beta2,beta3): phi_list = np.linspace(0.0,m.pi,1000) dldphi_list_1 = []
dldphi_list_2 = [] dldphi_list_3 = [] for phi in phi_list:
dldphi_list_1.append( luminosity_change_star_not_normalized_phi_beta(phi,beta1) ) ,→ dldphi_list_2.append( luminosity_change_star_not_normalized_phi_beta(phi,beta2) ) ,→ dldphi_list_3.append( luminosity_change_star_not_normalized_phi_beta(phi,beta3) ) ,→
matplotlib.rcParams.update({'font.size': 20}) plt.axvline(x=m.pi*0.5, linestyle=":",lw=3)
plt.plot(phi_list,dldphi_list_2,label=r'$\beta=0.2$',lw=3,ls='--') plt.plot(phi_list,dldphi_list_3,label=r'$\beta=0.5$',lw=3,ls='-.') plt.ylabel(r'$\frac{\mathrm{d}L}{\mathrm{d}\phi}$')
plt.xlabel(r'$\phi$ in radians',labelpad=20) plt.legend(loc="best")
plt.show()
def plot_dldphi_phi_beta_polar(beta1,beta2,beta3): phi_list = np.linspace(0.0,m.pi,1000)
dldphi_list_1 = [] dldphi_list_2 = [] dldphi_list_3 = [] for phi in phi_list:
dldphi_list_1.append( luminosity_change_star_not_normalized_phi_beta(phi,beta1) ) ,→ dldphi_list_2.append( luminosity_change_star_not_normalized_phi_beta(phi,beta2) ) ,→ dldphi_list_3.append( luminosity_change_star_not_normalized_phi_beta(phi,beta3) ) ,→ dldphi_list_1.extend(dldphi_list_1) dldphi_list_2.extend(dldphi_list_2) dldphi_list_3.extend(dldphi_list_3)
phi_list_1 = np.append(phi_list,np.linspace(2.0*m.pi-phi_list[0],2.c
0*m.pi-phi_list[999],1000))
,→
phi_list_2 = np.append(phi_list,np.linspace(2.0*m.pi-phi_list[0],2.c
0*m.pi-phi_list[999],1000))
,→
phi_list_3 = np.append(phi_list,np.linspace(2.0*m.pi-phi_list[0],2.c
0*m.pi-phi_list[999],1000))
,→
matplotlib.rcParams.update({'font.size': 20})
ax = plt.subplot(111,polar=True) # Introduces polar axes ax.set_theta_direction(1)
ax.set_rlabel_position(180)
ax.set_xticklabels([r'0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$', r'$5\pi/4$', r'$3\pi/2$', r'$7\pi/4$'])
,→ ax.set_theta_zero_location('S') ax.grid(True) ax.arrow(0,0,0,3,head_width=0.05) ax.annotate(xy=(0.05,3),s=r'$\vec{v}$') ax.plot(phi_list_1,dldphi_list_1,label=r'$\beta=0$',lw=3) ax.plot(phi_list_2,dldphi_list_2,label=r'$\beta=0.2$',lw=3,ls='--') ax.plot(phi_list_3,dldphi_list_3,label=r'$\beta=0.5$',lw=3,ls='-.') plt.legend(loc='best') plt.show() def L_tot(M,h): '''
'''
integral = simpson(lambda phi:luminosity_change_star_not_normalizedc _phi(phi,M,h),0,m.pi,1000)
,→
return integral
def luminosity_change_star2(zeta,h,M): '''
dl/d(phi)* as a function of zeta '''
I_0 = 1.0
function = (gamma(M,h)**(5.0) * (1.0- (beta(M,h)*m.cos(m.pi-zeta) ))**4.0 * (7*m.pi))**(-1) *6.0*I_0 * m.sin(m.pi-zeta)**(2.0) *
(1.0+ (0.5*m.pi * m.sin(m.pi
-zeta)/(gamma(M,h)*(1.0-beta(M,h)*m.cos(m.pi-zeta)))))
,→ ,→ ,→
return function
def luminosity_change_star2_beta_0(zeta): '''
dl/d(phi)* as a function of zeta. '''
I_0 = 1.0
function = ((7*m.pi))**(-1) *6.0*I_0 * m.sin(m.pi-zeta)**(2.0) *
(1.0+ (0.5*m.pi * m.sin(m.pi - zeta)))
,→
return function def L_tot_2(M,h):
'''
Normalization of function should give 1.0. '''
integral = simpson(lambda
zeta:luminosity_change_star2(zeta,h,M),0.0,m.pi,1000)
,→
return integral
def heightdependence(H,h,function_number): '''
Describes the distribution of flux over the height of the column g(h). ,→ ''' if function_number == 0: # g(h) = 1/H return H**(-1.0) elif function_number == 1: # g(h) = 2h/H^2 return 2.0*h*H**(-2.0) else: # g(h) = (2/H)(1-h/H) return (2.0/H) - ( (2.0*h) / (H*H))
def flux_observer_col_1(zeta,H,M,R,function_number): '''
Calculates the directly observed flux from one accretion column. '''
integral_no_eclipse_col_1 = (0.5/m.pi)*simpson(lambda h:
(1.0/m.sin(m.pi-zeta))*luminosity_change_star2(zeta,h,M)*heightc dependence(H,h,function_number), 0.0 , H , 1000) # Description without eclipse, phi is pi-zeta
,→ ,→ ,→
integral_eclipse_col_1 = (0.5/m.pi)*simpson(lambda h:
(1.0/m.sin(m.pi-zeta))*luminosity_change_star2(zeta,h,M)*heightc dependence(H,h,function_number), hmin(zeta,R),H,1000) #
Description eclipse, phi is pi-zeta
,→ ,→ ,→
if zeta <= 0.5 * m.pi: # if accretion column is not eclipsed return max(0.0,integral_no_eclipse_col_1)
else: # if accretion column is partially eclipsed if hmin(zeta,R) > H:
return 0.0 else:
return max(0.0,integral_eclipse_col_1) def plot_flux_observer_col_1(M):
'''
Plots the directly observed flux for different accretion column heights.
,→
'''
zeta_list = np.arange(0.001,m.pi,(m.pi-0.001)/1000.0) flux_list = []
flux_list_2 = [] flux_list_3 = [] flux_list_4 = [] flux_list_5 = []
for i in range(len(zeta_list)):
flux_list.append(flux_observer_col_1(zeta_list[i],0.1,M,1.0,0) + flux_observer_col_1(m.pi-zeta_list[i],0.1,M,1.0,0)) ,→ flux_list_2.append(flux_observer_col_1(zeta_list[i],0.5,M,1.0,0c ) + flux_observer_col_1(m.pi-zeta_list[i],0.5,M,1.0,0)) ,→ ,→ flux_list_3.append(flux_observer_col_1(zeta_list[i],1.0,M,1.0,0c ) + flux_observer_col_1(m.pi-zeta_list[i],1.0,M,1.0,0)) ,→ ,→ plt.plot(zeta_list,flux_list,label=r'$H=0.1R$',lw=3,ls='--') plt.plot(zeta_list,flux_list_2,label=r'$H=0.5R$',lw=3,ls='-.') plt.plot(zeta_list,flux_list_3,label=r'$H=R$',lw=3,ls='-') plt.axvline(x=0.5*m.pi,lw=3,ls=':') plt.xlim(0.0,m.pi)
matplotlib.rcParams.update({'font.size': 20}) plt.xlabel(r'$\xi$ in radians',labelpad=20) plt.legend(loc="upper right")
plt.ylabel(r'Relative flux') plt.show()
factor = max(0.0,perp_factor(theta_b,h,R)) # Prevents the perpendicular component to be negative
,→
function = factor * (luminosity_change_star(theta_b,h,R,M) *
heightdependence(H,h,function_number) ) / (2.0*m.pi*distance_travelled(h,R,theta_b)**(2.0) * m.sin(phi(theta_b,h,R))) ,→ ,→ ,→ return function
def flux_inwards(theta_b,H,R,M,function_number): if theta_b >= m.pi*0.5:
theta_b = m.pi-theta_b
integral = simpson(lambda h: flux_inwards_function_general(thetc a_b,h,R,M,H,function_number), hmin_theta_b(theta_b,R) , H,1000)
,→ ,→
else:
integral = simpson(lambda h: flux_inwards_function_general(thetc a_b,h,R,M,H,function_number), hmin_theta_b(theta_b,R) , H,1000)
,→ ,→
return max(0.0,integral)
def flux_inwards_1ps_iso(h,R,theta_b): '''
Calculates the incident flux for a latitude theta_b considering 1 isotropically emitting point source.
,→
'''
flux = (1.0 / (4.0 * m.pi * distance_travelled(h,R,theta_b)**2.0) )*perp_factor(theta_b,h,R)
,→
return max(0.0,flux)
def flux_inwards_1ps_aniso(h,R,theta_b): '''
Calculates the incident flux for a latitude theta_b considering 1 anisotropically emitting point source.
,→
'''
factor = max(0.0,perp_factor(theta_b,h,R))
function = factor * luminosity_change_star(theta_b,h,R,M) /
(2.0*m.pi*distance_travelled(h,R,theta_b)**(2.0) * m.sin(phi(theta_b,h,R))) ,→ ,→ return function def plot_flux_inwards_1ps_iso(h,R): '''
Plots the incident flux as a function of latitude theta_b considering 1 isotropically emitting point source.
,→
'''
theta_b_list = np.linspace(0.0,0.5*m.pi,1000) flux_list = []
for theta_b in theta_b_list:
flux_list.append(10.0**(4.0)*flux_inwards_1ps_iso(h,R,theta_b)) flux_list_2.append(flux_inwards_1ps_iso(1.0,1.0,theta_b))
matplotlib.rcParams.update({'font.size': 20}) plt.xlim(0.0,m.pi*0.5)
plt.plot(theta_b_list,flux_list_2,label=r'$h=R$',lw=3,ls='-') plt.plot(theta_b_list,flux_list,label=r'$h=100R$',lw=3,ls='--') plt.legend(loc="best")
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0)) plt.xlabel(r'$\theta_B$ in radians',labelpad=20)
plt.ylabel(r'Relative flux') plt.show()
def plot_flux_inwards_2ps_iso(h,R): '''
Plots the incident flux as a function of latitude theta_b considering 2 isotropically emitting point sources.
,→
'''
theta_b_list = np.linspace(0.0,m.pi,1000) flux_list = []
flux_list_2 = []
for theta_b in theta_b_list:
flux_list.append(
10.0**(4.0)*(flux_inwards_1ps_iso(h,R,theta_b)
+flux_inwards_1ps_iso(h,R,m.pi-theta_b)) )
,→ ,→
flux_list_2.append(flux_inwards_1ps_iso(1.0,1.0,theta_b)+flux_ic nwards_1ps_iso(1.0,1.0,m.pi-theta_b))
,→
matplotlib.rcParams.update({'font.size': 20})
plt.plot(theta_b_list,flux_list_2,label=r'$h=R$',lw=3)
plt.plot(theta_b_list,flux_list,label=r'$h=100R$',ls='--',lw=3) plt.axvline(x=0.5*m.pi,linestyle=":",lw=3)
plt.xlim(0.0,m.pi)
plt.legend(loc="lower right")
plt.xlabel(r'$\theta_B$ in radians',labelpad=20) plt.ylabel(r'Relative flux')
plt.show()
def plot_flux_inwards_2ps_aniso(h,R): '''
Plots the incident flux as a function of latitude theta_b considering 2 anisotropically emitting point sources.
,→
'''
theta_b_list = np.linspace(0.001,m.pi-0.001,1000) flux_list = []
flux_list_2 = []
