# Polarization

This is a short tutorial for modeling polarized X-ray flux from a rapidly rotating neutron star.

Let us first make the initial setup (see the Modeling tutorial for more detailed explanations):

[1]:

%matplotlib inline

from __future__ import print_function, division

import os
import numpy as np
import math
import time

from matplotlib import pyplot as plt
from matplotlib import rcParams
from matplotlib.ticker import MultipleLocator, AutoLocator, AutoMinorLocator
from matplotlib import gridspec
from matplotlib import cm

import xpsi

from xpsi.global_imports import _c, _G, _dpr, gravradius, _csq, _km, _2pi

freq = 600.0

bounds = dict(distance = (0.1, 1.0),                     # (Earth) distance
mass = (1.0, 3.0),                       # mass
cos_inclination = (0.0, 1.0))            # (Earth) inclination to rotation axis

spacetime = xpsi.Spacetime(bounds=bounds, values=dict(frequency=freq))

ceding = False

bounds = dict(super_colatitude = (None, None),
phase_shift = (0.0, 0.1),
super_temperature = (None, None))

if ceding:

bounds = dict(super_colatitude=(None,None),
phase_shift = (0.0, 1.0),
super_temperature = (None, None),
cede_colatitude = (None, None),
cede_azimuth = (None, None),
cede_temperature = (None, None))

# a simple circular, simply-connected spot
primary = xpsi.HotRegion(bounds=bounds,
values={}, # no initial values and no derived/fixed
symmetry=True,
#symmetry=False,
omit=False,
cede=ceding,
concentric=False,
sqrt_num_cells=32,
min_sqrt_num_cells=10,
max_sqrt_num_cells=64,
#num_leaves=100,
num_leaves=121,
num_rays=200,
do_fast=False,
atm_ext="Pol_BB_Burst",
prefix='p')

class derive(xpsi.Derive):
def __init__(self):
"""
We can pass a reference to the primary here instead
and store it as an attribute if there is risk of
the global variable changing.

This callable can for this simple case also be
achieved merely with a function instead of a magic
method associated with a class.
"""
pass

def __call__(self, boundto, caller = None):
# one way to get the required reference
global primary # unnecessary, but for clarity
return primary['super_temperature'] - 0.2

bounds['super_temperature'] = None # declare fixed/derived variable

secondary = xpsi.HotRegion(bounds=bounds, # can otherwise use same bounds
values={'super_temperature': derive()}, # create a callable value
symmetry=True,
omit=False,
cede=ceding,
concentric=False,
sqrt_num_cells=32,
min_sqrt_num_cells=10,
max_sqrt_num_cells=100,
num_leaves=100,
num_rays=200,
do_fast=False,
atm_ext="Pol_BB_Burst", #the simplest polarized atmosphere extension
is_antiphased=True,
prefix='s') # unique prefix needed because >1 instance

from xpsi import HotRegions

hot = HotRegions((primary, secondary))

h = hot.objects[0]

hot['p__super_temperature'] = 6.0 # equivalent to primary['super_temperature'] = 6.0

class CustomPhotosphere(xpsi.Photosphere):
""" Implement method for imaging."""

@property
def global_variables(self):
""" This method is needed if we also want to ivoke the image-plane signal simulator. """

return np.array([self['p__super_colatitude'],
self['p__phase_shift'] * _2pi,
self['p__super_temperature'],
self['s__super_colatitude'],
(self['s__phase_shift'] + 0.5) * _2pi,
self.hot.objects[1]['s__super_temperature']])

bounds = dict(spin_axis_position_angle = (None, None))
#Set here stokes=True to calculate the polarized X-ray signal.
photosphere = CustomPhotosphere(hot = hot, elsewhere = None, stokes=True, bounds=bounds,
values=dict(mode_frequency = spacetime['frequency']))

print(photosphere['mode_frequency'] == spacetime['frequency'])

star = xpsi.Star(spacetime = spacetime, photospheres = photosphere)

#Let us first initialize the star with these values but change them later:
p = [1.4,
12.5,
0.2,
math.cos(1.25),
0.0,
0.0,
1.0,
0.075,
6.2,
0.025,
math.pi - 1.0,
0.2]
if ceding:
p = [1.4, #mass
0.2, #distance
math.cos(1.25), #cos_inclination
0.0, # spin_axis_position_angle
0.0, #p__phase_shift
1.0, #p__super_colatitude
6.2, #p__super_temperature
0.1, #p__cede_colatitude
0.0, #p__cede_azimuth
6.2, #p__cede_temperature
0.025, #s__phase_shift
math.pi - 1.0, #s__super_colatitude
math.pi-1.0, #s__cede_colatitude ..
0.0, #s__cede_azimuth
6.2] #s__cede_temperature

star(p)

#These are the values used in the light curve computation:
star['mass'] = 1.4 #2.7088795
star['distance'] = 0.2
incl_deg = 40.0 #90.0 #40.0
star['cos_inclination'] = math.cos(math.pi*incl_deg/180.0)#math.cos(2.0)
theta_deg = 60.0
star['p__super_colatitude'] = math.pi*theta_deg/180.0 #0.0 #2.0
rho_deg = 10.0 #10.0 #0.001 #10.0
tplanck = 1.0219978 #1.0 #in keV #1 keV -> log10(T[K]) = 7.06 (out of bounds originally)
star['p__super_temperature'] = np.log10(tplanck*11604525.0061657)

if ceding:
star['p__cede_colatitude'] = math.pi*theta_deg/180.0
star['p__cede_azimuth'] = 0.0 #0.0
star['p__cede_temperature'] = np.log10(tplanck*11604525.0061657)
star['s__super_colatitude'] = math.pi-math.pi*theta_deg/180.0
star['s__cede_colatitude'] = math.pi-math.pi*theta_deg/180.0
star['s__cede_azimuth'] = 0.1 #math.pi/2.0
star['s__cede_temperature'] = np.log10(tplanck*11604525.0061657)

rcParams['text.usetex'] = False
rcParams['font.size'] = 14.0

def veneer(x, y, axes, lw=1.0, length=8):
""" Make the plots a little more aesthetically pleasing. """
if x is not None:
if x[1] is not None:
axes.xaxis.set_major_locator(MultipleLocator(x[1]))
if x[0] is not None:
axes.xaxis.set_minor_locator(MultipleLocator(x[0]))
else:
axes.xaxis.set_major_locator(AutoLocator())
axes.xaxis.set_minor_locator(AutoMinorLocator())

if y is not None:
if y[1] is not None:
axes.yaxis.set_major_locator(MultipleLocator(y[1]))
if y[0] is not None:
axes.yaxis.set_minor_locator(MultipleLocator(y[0]))
else:
axes.yaxis.set_major_locator(AutoLocator())
axes.yaxis.set_minor_locator(AutoMinorLocator())

axes.tick_params(which='major', colors='black', length=length, width=lw)
axes.tick_params(which='minor', colors='black', length=int(length/2), width=lw)
plt.setp(axes.spines.values(), linewidth=lw, color='black')

def plot_pulse(stokes="I"):
""" Plot hot region signals before telescope operation. """
fig = plt.figure(figsize=(7,7))

ax.set_ylabel('Signal [arbitrary normalisation]')
ax.set_xlabel('Phase [cycles]')

if stokes=="I":
temp1 = np.sum(photosphere.signal[0][0], axis=0)
temp2 = np.sum(photosphere.signal[1][0], axis=0)
elif stokes=="Q":
temp1 = np.sum(photosphere.signalQ[0][0], axis=0)
temp2 = np.sum(photosphere.signalQ[1][0], axis=0)
elif stokes=="U":
temp1 = np.sum(photosphere.signalU[0][0], axis=0)
temp2 = np.sum(photosphere.signalU[1][0], axis=0)
else:
print("ERROR: Stokes option must be 'I', 'Q', or 'U'")
exit()

ax.plot(hot.phases_in_cycles[0], temp1/np.max(temp1), 'o-', color='k', lw=0.5, markersize=2)
ax.plot(hot.phases_in_cycles[1], temp2/np.max(temp2), 'o-', color='r', lw=0.5, markersize=2)

veneer((0.05,0.2), (0.05,0.2), ax)
#fig.savefig("figs/pulse_profile"+stokes+"X.pdf")

nene = 281 #128

from numpy import logspace
from numpy import log
evere=0.5109989e6 # electron volts in elecron rest energy
x_l, x_u = -3.7 , .3 # lower and upper bounds of the log_10 energy span
IntEnergy = logspace(x_l,x_u,nene), log(1e1)*(x_u-x_l)/(nene-1.)
x, x_weight = IntEnergy
energies = (x*evere)/1e3

star.update() #Calculating the space-time integrals etc.

import time

primary.image_order_limit = 1

/=============================================\
| X-PSI: X-ray Pulse Simulation and Inference |
|---------------------------------------------|
|                Version: 2.2.3               |
|---------------------------------------------|
|      https://xpsi-group.github.io/xpsi      |
\=============================================/

Imported GetDist version: 1.4.7
Imported nestcheck version: 0.2.1
Creating parameter:
> Named "frequency" with fixed value 6.000e+02.
> Spin frequency [Hz].
Creating parameter:
> Named "mass" with bounds [1.000e+00, 3.000e+00].
> Gravitational mass [solar masses].
Creating parameter:
> Named "radius" with bounds [4.430e+00, 1.600e+01].
Creating parameter:
> Named "distance" with bounds [1.000e-01, 1.000e+00].
> Earth distance [kpc].
Creating parameter:
> Named "cos_inclination" with bounds [0.000e+00, 1.000e+00].
> Cosine of Earth inclination to rotation axis.
Creating parameter:
> Named "super_colatitude" with bounds [0.000e+00, 3.142e+00].
> The colatitude of the centre of the superseding region [radians].
Creating parameter:
> Named "super_radius" with bounds [0.000e+00, 1.571e+00].
Creating parameter:
> Named "phase_shift" with bounds [0.000e+00, 1.000e-01].
> The phase of the hot region, a periodic parameter [cycles].
Creating parameter:
> Named "super_temperature" with bounds [3.000e+00, 7.600e+00].
> log10(superseding region effective temperature [K]).
Creating parameter:
> Named "super_colatitude" with bounds [0.000e+00, 3.142e+00].
> The colatitude of the centre of the superseding region [radians].
Creating parameter:
> Named "super_radius" with bounds [0.000e+00, 1.571e+00].
Creating parameter:
> Named "phase_shift" with bounds [0.000e+00, 1.000e-01].
> The phase of the hot region, a periodic parameter [cycles].
Creating parameter:
> Named "super_temperature" that is derived from ulterior variables.
> log10(superseding region effective temperature [K]).
Creating parameter:
> Named "mode_frequency" with fixed value 6.000e+02.
> Coordinate frequency of the mode of radiative asymmetry in the
photosphere that is assumed to generate the pulsed signal [Hz].
Creating parameter:
> Named "spin_axis_position_angle" with bounds [-1.571e+00, 1.571e+00].
> Spin axis position angle measured from the north counterclock-
wise to the projection of the rotation axis on the plane of the
True


And then let us make the final integration and plot the results for all Stokes parameters.

We first check Stokes I profiles (integrated over all energies) for the primary and secondary hot spots:

[2]:

start = time.time()
end = time.time()
print("Time spent in integration (unpolarized):",end - start)

plot_pulse(stokes="I")

Time spent in integration (unpolarized): 3.1238365173339844


Let us also check Stokes Q and U profiles:

[3]:

plot_pulse(stokes="Q")
plot_pulse(stokes="U")


## Creating synthetic data

To produce polarized synthetic data, we can couple X-PSI to ixpeobssim as an input model. The input model file needs to be placed in the config folder of the ixpeobssim repository . An example of the input model file is provided here. To prodcue the synthetic Stokes profiles and their associated errors one should execute source setup.sh in the main directory of ixpeobssim and then run a script like simulate_amsp_xpsi.py found here. For more details, check the ixpeobssim documentation pages.

## Fitting the data

Modeling X-ray polarized data using Bayesian inference is not yet part of this tutorial. However, a preliminary example of a such procedure can be found in this example script: TestRun_PolNum_split_inference.py.

[ ]: