Cyclotron resonant scattering feature simulations. II. Description of the CRSF simulation process
F. -W. Schwarm R. Ballhausen S. Falkner G. Schönherr K. Pottschmidt M. T. Wolff P. A. Becker F. Fürst D. M. Marcu-Cheatham P. B. Hemphill E. Sokolova-Lapa T. Dauser D. Klochkov C. Ferrigno J. Wilms




II. Description of the CRSF simulation process


F.-W. Schwarm et al.

Dr. Karl Remeis-Sternwarte and Erlangen Centre for Astroparticle Physics, Sternwartstrasse 7, 96049 Bamberg, Germany
Leibniz-Institut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
CRESST, Department of Physics, and Center for Space Science and Technology, UMBC, Baltimore, MD 21250, USA
NASA Goddard Space Flight Center, Code 661, Greenbelt, MD 20771, USA
Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA
Department of Physics & Astronomy, George Mason University, Fairfax, VA 22030-4444, USA
European Space Astronomy Centre (ESA/ESAC), Science Operations Department, Villanueva de la Cañada (Madrid), Spain
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Faculty of Physics, M. V. Lomonosov Moscow State University, Leninskie Gory, Moscow 119991, Russia
Sternberg Astronomical Institute, Moscow M. V. Lomonosov State University, Universitetskij pr., 13, Moscow 119992, Russia
Institut für Astronomie und Astrophysik, Universität Tübingen (IAAT), Sand 1, 72076 Tübingen, Germany
ISDC Data Center for Astrophysics, Université de Genève, chemin d’Écogia 16, 1290 Versoix, Switzerland





	Cyclotron resonant scattering features (CRSFs) are formed by scattering
of X-ray photons off quantized plasma electrons in the strong magnetic field
(of the order  10^12 G) close to the surface of 
an accreting X-ray pulsar. Due to the complex scattering cross sections the line
profiles of CRSFs cannot be described by an analytic expression. Numerical
methods such as Monte Carlo (MC) simulations of the scattering processes
are required in order to predict precise line shapes for
a given physical setup, which can be compared to observations to gain information about the
underlying physics in these systems.


	A versatile simulation code is needed for the generation of synthetic
cyclotron lines. Sophisticated geometries
should be investigatable by making their simulation
possible for the first time.


	The simulation utilizes the mean free path tables described in
the first paper of this series for the fast interpolation of propagation
lengths. The code is parallelized to make the very time consuming simulations
possible on convenient time scales. Furthermore, it can generate responses to
mono-energetic photon injections, producing Green's functions,
which can be used later to generate spectra for arbitrary continua.


	We develop a new simulation code to generate synthetic cyclotron
lines
for
complex scenarios, allowing for unprecedented
physical interpretation of the observed data. An associated XSPEC 
model implementation is used to fit synthetic line profiles
to NuSTAR data of Cep X-4.
The code has
been developed with the main goal of overcoming previous geometrical
constraints in MC simulations of CRSFs. By applying this code also to
more simple, classic geometries used in previous works, we furthermore address
issues of code verification and cross-comparison of various models. The
XSPEC model and the Green's function tables are available online
(see link in footnote, page 1).




Cyclotron resonant scattering feature simulations
    F.-W. Schwarm1,
	R. Ballhausen1,
	S. Falkner1,
	G. Schönherr2,
	K. Pottschmidt3,4,
	M. T. Wolff5,
	P. A. Becker6,
	F. Fürst7,
	D. M. Marcu-Cheatham3,4,
	P. B. Hemphill8,
	E. Sokolova-Lapa9,10,
	T. Dauser1,
	D. Klochkov11,
	C. Ferrigno12,
	J. Wilms1

    Received December 14, 2016; accepted January 23, 2017
=============================================================================================================================================================================================================================================================





§ INTRODUCTION

<http://www.sternwarte.uni-erlangen.de/research/cyclo>
Cyclotron resonant scattering features (CRSFs, or “cyclotron lines”)
result from the interaction of photons with electrons in a strong
magnetic field (B≳ 10^12 G), for example, in the accretion
column of the magnetized neutron star in accreting X-ray binaries.
Their complex shape <cit.> reveals information about
the environment of the line forming region, which is typically close
to the neutron star.
A physical model for the line formation is needed in order to translate
observed line characteristics — like the line's depth, width, and shape —
into physically meaningful parameters. 

This paper is the second in a series in which we describe calculations of
resonant cyclotron scattering for various configurations of the
accreting matter. In the first paper of the series <cit.>, we described the calculation of the
relevant cross section for the range of conditions typically expected
for accreting neutron stars in X-ray binary systems. In this work we
discuss the application of these cross sections to the simulation of synthetic cyclotron lines.
Simulations of CRSFs have been performed by various groups in the
past. They can be divided into three classes: Monte Carlo
<cit.>,
Feautrier <cit.>,
and semianalytic methods <cit.>. <cit.> provide an
overview over the methods used in previous works. 
<cit.> discuss the features of the fundamental approaches and
use multiple methods for the generation
of synthetic cyclotron lines in various parameter regimes.
For optical depths on the order of
10^-4–10^-3τ_T, where τ_T is the Thomson optical
depth, Monte Carlo (MC) simulations are the suitable method
<cit.>.
From all of the calculations mentioned above only <cit.>
provided a method for direct comparison between data and observations with
the corresponding large grids of parameters.

<cit.> performed a first comparison of simulated MC line
profiles to data from the X-ray pulsar A0535+26. By convolving
precalculated Green's functions, describing the response of a CRSF
medium to monoenergetic photon injection, with arbitrary
continua, <cit.> disentangled the time consuming MC
simulations from the choice of continuum. This approach enabled
us to develop a code which can be used to compare MC simulated spectra
to observational data using standard X-ray astronomical data analysis
packages such as XSPEC <cit.> or ISIS <cit.>. A
disadvantage of these earlier calculations, however, is that they are
limited to very simple predefined geometries. In addition, while the
general properties of the line behavior are correct, unfortunately
there was an error in the integration of the scattering cross section
code by <cit.>. 

Building on previous results and experience we have developed a new
simulation tool to calculate synthetic line profiles for arbitrarily
complex cylindrically symmetric geometries. This generalization naturally
requires considerable computational effort, which must be met with new
strategies. The mean free path (MFP) tables used for the interpolation
of the thermally averaged scattering cross sections have been
described in paper I. Here, we will provide the description and
applications of the full MC scattering code, which has been written
with the prime goal of imprinting cyclotron lines on the continuum
emission of accreting X-ray pulsars, and which includes a working
fit model. The XSPEC model code and the Green's tables necessary to
generate synthetic cyclotron lines are provided, as well as instructions
on how to use them. The discussion of more sophisticated physical
scenarios and the application to observational data will be the subject of successive papers in this series.
Here, we restrict ourselves to presenting examples of selected
synthetic spectra for illustration and for comparison to other works.

The structure of this paper is as follows. In Sect. <ref>
the physical assumptions are laid out and the simulation process is explained
in detail, including a description of the scattering process and commonly used
geometries. The duration of an exemplary simulation run is analyzed with
respect to the used number of CPUs and computer systems. This is followed by a
description of the Green's function scheme, which is utilized to enable
the fitting of the synthetic cyclotron lines to observational data, taking
a step toward the application of the simulation in Sect. <ref>.
In the same section, a comparison of synthetic CRSFs to previous works is performed.
Furthermore, the Green's functions are applied to a physical continuum
for several combinations of line parameters. Also, the availability of the
XSPEC model and the associated Green's function tables is discussed.
The XSPEC model is fitted to NuSTAR data of Cen X-4
and the resulting magnetic field strength and cyclotron line shape
is compared to results from empirical models.
Additional technical details about the model usage and its parameters are
given in the Appendix.



§ SIMULATION OF CYCLOTRON LINES

In this section we first describe the physical setup and explain our
Monte Carlo approach in detail before discussing how these simulations
can be sped up. We also introduce various classical accretion
geometries which we will be using in Sect. <ref> in
order to validate our results with respect to the ones from earlier works using
the MC method.



 §.§ Physical setup

We consider the interactions of photons with electrons in strong
magnetic fields via the cyclotron scattering process. Other particles,
such as protons which might form “proton cyclotron lines”
<cit.> in the spectra at energies below the
electron cyclotron lines, are neglected. There are several other
processes altering the properties and paths of X-ray photons in the
vicinity of magnetic fields near the critical field strength <cit.>,

    B_crit=
    m_e^2 c^3/eħ = 4.413× 10^13 G .

Photon splitting, which describes the process of splitting one photon
into two photons, dominates over Compton scattering for very high
magnetic field strengths and densities <cit.>. Pair
production becomes possible if the photon energy exceeds 2m_ec^2
<cit.>. Apart from photon-electron interactions,
<cit.> also analyzed the properties of Møller and Bhabha
scattering where electrons are scattered off electrons or positrons,
respectively. In the regime we are interested in, cyclotron scattering
is the dominant process and therefore we will neglect all the other
possible interactions.

As further discussed in paper I, the momenta of the CRSF forming
electrons perpendicular to the direction of the B-field are
quantized to discrete values corresponding to the Landau energy
levels, .
The electrons can move freely parallel to the field
lines. Inelastic scattering of photons off these electrons leads to
the formation of CRSFs which appear as absorption-like lines at the
fundamental Landau energy E_1 and its integer multiples <cit.>.

The cyclotron scattering cross section strongly depends on the
incident photon angle and energy <cit.>.
The angle ϑ_in is measured with respect to the magnetic
field axis and mostly specified by its cosine, μ_in = cosϑ_in.
In the rest frame of
an electron, the photons' energies and propagation angles are Lorentz
boosted due to the electrons' motion parallel to the B-field. This
makes the accurate analytic treatment difficult <cit.>
because a photon's energy and angle are relativistically coupled in
the electron rest frame.

Cyclotron absorption and emission processes as well as resonant
scattering generate highly complex line profiles. In contrast to Compton
scattering which does not absorb or create photons, the number of
photons that contribute to the spectrum is not conserved
<cit.>. A photon can be absorbed by an electron if it has
the right energy and angle to excite the electron to a higher Landau
level. Almost immediately, the electron will emit “spawned” photons
during its successive de-excitation to the ground state. Since
de-excitation preferentially takes place to the next lower Landau
level <cit.>, the majority of spawned photons will have
energies corresponding to the energy difference between neighboring
Landau levels. This is close to the fundamental energy E_1
regardless of the initial Landau level of the electron. The transition
process continues until the ground state is reached, effectively
increasing the number of photons within the medium.
Therefore complex computational
methods are necessary to determine the exact shape of cyclotron
resonance features.



 §.§ Monte Carlo simulation


The advantage of a Monte Carlo simulation is that its photon tracing approach
allows for the simulation of the radiative transfer in very complex setups of
the scattering medium. In general, in a MC simulation a seed photon is
generated at a place where the primary photons originate. The photon is
assigned an energy, position, and direction. A random number is then drawn from
an exponential distribution, which depends on the photon mean free path. This
random number is the optical depth that the photon will travel.  As discussed
in paper I, in order to speed this step up we use precalculated mean free path
tables (see Fig. <ref> and Sect. <ref> below for
more details).

The scattering geometry is realized in our simulation by describing the
medium with a list of cylinders with arbitrary dimensions, which may be combined
to model all kinds of cylindrically symmetric shapes of accretion columns. The
geometry of each cylinder is parameterized by its inner — to allow for hollow
cylinders — and outer radii and the distances of its bottom and top to the
neutron star surface. The physical properties inside are given by its
homogeneous density[The density is usually calculated from the optical
depth into a given direction and the corresponding extension of the cylinder.
For the classical “” and “” geometries the
direction parallel to the B-field axis is chosen, while in the “cylinder”
geometry the optical depth is defined perpendicular to this axis (see
Fig. <ref>)], magnetic field, electron temperature, and velocity
towards the neutron star. Multiple cylinders can be stacked on top of, or
inside of each other in order to properly simulate parameter gradients.  Here,
we still use the simple geometries, which are explained in
Fig. <ref> together with their implied seed emission patterns,
for comparison of our results to <cit.> in Sec. <ref>.

The seed photons in a simulation run result from the
configured photon sources. Different emission patterns are available
for each individual source: a point source emits
photons isotropically from its static origin, a line source
corresponds to photon emission from a line along part of
the magnetic field axis, meaning that the height of each photon above
the hypothetical neutron star is sampled individually, a plane source
describes an emitting plane perpendicular to the B-field axis at a
given height. Variants emitting only photons upwards, that is, in the
direction of the B-field axis, or downwards are available for the
point and plane source types to provide a convenient way of preventing
unprocessed photons from showing up in the resulting data. 

The implementation of photon spawning is straightforward: successive photon
generation and propagation of the spawned photons from the coordinates of the
de-exiting electron to the point where the spawned photons are interacting with
or escaping the medium ensure self-consisted treatment.

The most time-intensive parts of processing are parallelized using
the <cit.>, decreasing the required CPU time efficiently as shown in
Fig. <ref>.
Because available computing power has increased significantly since
earlier approaches to this problem, we were also able to introduce
additional features to the simulation which allow us to specify the
conditions to be simulated in a more flexible way. For example, more
complex physical settings for the accretion column, including velocity
gradients, B-field gradients, and an inhomogeneous density
stratification.



 §.§ Description of the scattering process


The core of the MC simulation is the propagation of individual photons
through a given CRSF medium, the process of which is illustrated by
Fig. <ref>. Seed photons, generated by the configured
photon sources, must be processed as well as spawned photons,
originating from the transition of an electron to a lower Landau
level. For this purpose the seed photons are injected into the
simulation and photons spawned during their interactions are
added to the current list of photons and processed the same way.

Each photon is propagated according to its current path length, which
is sampled from the mean free path tables provided by
<cit.>. The photon is stored to a file in FITS format if
it escapes the medium and the simulation starts with the next seed
photon. Otherwise the momentum of the interacting electron is sampled,
again using precalculated tables (paper I). We assume that the
electrons are in the ground state, that is, we neglect the
excitation of electrons through other processes than resonant cyclotron
scattering, which is justified by the comparably small time scales of
the electron decay compared to the collisional excitation rate <cit.>.
They can be excited to higher
Landau levels in the course of resonant cyclotron scattering, though. The final Landau level
for each interaction is sampled by comparing the scattering cross
sections for all possible final levels below the maximum Landau level
n_f,max, which has been set to 5 for the simulations
performed here, in compliance with the number of Landau levels taken
into account for the MFP table calculations. In a similar manner, the
scattering angle is sampled by assigning a random number to the
cumulative angular distribution gained by integrating the cross
section over all possible final photon angles. The kinematic
calculations can be carried out once the initial properties of the
interacting electron and photon have been sampled, gaining the final
photon energy and electron momentum. The photon altered by the
scattering process is now further propagated, starting again with the
mean free path interpolation.

The excited electron will produce spawned photons during its
de-excitation to the ground state if it has been excited to a higher
Landau level during the scattering process. The number of spawned
photons depends on the Landau level the electron is excited to, which
is mainly determined by the photon energy in the electron rest frame.
De-excitation preferentially takes place to the next lower Landau
level, producing a photon with an energy corresponding to the energy
difference of the involved Landau levels. 
In order to show this quantitatively, we calculated the transition rates for a magnetic field of 0.1 B_crit, averaged over initial and final spin and
summed over final polarization. They are shown in Table <ref> for initial Landau levels
up to n_i = 7. The sixths and seventh levels are provided as a reference and in order to allow
for estimations of the probability that a photon is spawned at the resonance energy of the maximum
Landau level taken into account in the simulation, for instance, via the transition 7 → 2.
The ratio between the rate for a transition from Landau level n_i = 2 to Landau level
n_f = 1, that is, Γ_21≈ 1.74 × 10^-3, and the corresponding
transition to the ground state, Γ_20≈ 2.44 × 10^-4, becomes
Γ_21/Γ_20≈ 7.11 in agreement with the
corresponding calculation by <cit.>. The ratio
is becoming smaller for larger fields as shown by <cit.>, as well.



Although the spacing of the
cyclotron resonances is not perfectly harmonic these photons have
energies very close to the fundamental energy in the electron's rest
frame. For Landau levels above the first excited state an electron
emits multiple photons during its successive de-excitation to the next
lower level until the ground state is reached. The spawned photons
have to be boosted to the neutron star frame of reference before they are further
processed, since the electron that re-emits a photon has some velocity
component parallel to the magnetic field due to the electrons' thermal momentum
distribution (see paper I).
The simulation code features the possibility to configure an additional velocity
component of the simulated medium. This leads to an additional boosting factor
for all photons entering the medium and can be used to simulate the bulk velocity
of the accreting matter. This has been
used by <cit.> for the calculation of phase dependent cyclotron
line spectra and emission patterns from a cylindrically symmetric accretion
column throughout which such a flow with relativistic velocities is expected
<cit.>.
We neglect bulk motion in this work because it is unnecessary for the
purpose of code verification. For the application of the CRSF model to Cep X-4
the omission of bulk velocity is justified by the chosen geometry: in “slab”
geometries the line forming region is assumed to reside at the neutron star hot
spot at which the matter has been decelerated to non-relativistic velocities.
The boosting caused by the thermal momentum distribution of the electrons, together with
the slightly anharmonic spacing of the Landau levels, is responsible
for the formation of line wings around the fundamental cyclotron line.
This apparent excess of photons occurs at energies slightly
above or below the first resonance E_1. The intermediate Landau
levels occupied by the electron during this process are sampled by
making use of the corresponding cyclotron decay rates. Each transition
produces a photon with an energy equal to the energy difference of the
Landau levels. Its angle is sampled according to the angle
dependent decay rates. All spawned photons are further propagated in
the same way as the seed photons until they leave the medium.



 §.§ Classic geometries



Figure <ref> depicts various geometries which differ in
radius to height ratio and the position, and type, of the seed photon source.

Slab geometries mimic a thin scattering volume of infinite radius. We
approximate this in our numerical simulations by setting the radius to a value
very large compared to the height of the medium. We find empirically that a
factor of 1000 between the radius and the height is a good choice in the
simulations. The output spectrum does not change significantly for a larger
extension in the radial direction, that is, only few photons, if any,
travel that far perpendicular to the B-field axis without escaping through
the bottom or top of the slab. The height of the slab is specified in terms of
the Thomson optical depth. The density in the slab can therefore be calculated in
the simulation from the slab height or radius, for /
or  geometry, respectively, and the corresponding optical depth. For slab
geometries the optical depth parallel to the magnetic field is used to describe
the medium, because the slab radius is infinite or at least much larger than
its height. For a cylinder geometry the optical depth perpendicular to the
magnetic field axis is used.

Using the notation of <cit.>, in the 
geometry the medium is solely illuminated from the bottom, while in
the  geometry a source plane in the middle produces seed
photons within the medium.

In the cylinder geometry the medium's radius is much smaller than its
height. Therefore, the optical depth perpendicular to the B-field
axis, that is, the cylinder radius in units of optical depths is used
as a parameter to describe the medium instead. Here, the seed photons
are also produced within the column but only along the cylinder's
axis.

In this work we show spectra for  and 
geometry in Fig. <ref> and spectra for cylinder
geometry in Fig. <ref>.
The CRSF shape is, in general, highly
sensitive to the simulated geometry. The formation of line wings, for
example, is especially pronounced in  geometry for
viewing angles almost parallel to the magnetic field, which can be
seen in Fig. <ref>.



 §.§ CPU time


We reduce the calculation time by using MPI to
distribute the work among multiple processors.
Figure <ref> shows that our simulation is parallelized
efficiently. We performed tests on different computer systems and
found that using ∼8 CPUs for each simulation provides a
convenient compromise between simulation time and efficiency.



 §.§ Green's table approach

Performing a time consuming simulation for each generation of a
spectrum is not practicable in anticipation of our goal to fit
simulated spectra to observational data. A large number of simulation
runs will be required where the shape of the emerging spectrum is
simulated as a function of the input parameters. In addition to the
magnetic field strength and the plasma temperature, the spectral shape
of the incident photons is also variable. This significantly increases
the computing time as potentially a large parameter space needs to be
covered. As it was shown previously <cit.>,
this large expense of computing time can be avoided by calculating the
Green's function of the radiative transfer problem. We propagate
photons through a medium with a given geometry, magnetic field
strength, and temperature for a well sampled grid of mono-energetic
seed photon energies. The photons escaping the medium are collected
and binned to spectra for a grid of output angles. These emerging
spectra, normalized to the seed photon flux in a given
angular bin, correspond to the Green's function of the radiative
transfer problem. The emerging spectrum for an arbitrary seed photon
spectrum, like a cutoff power law continuum or a continuum due to
Comptonization in a radiative shock <cit.>, can then be
obtained by convolving it with the Green's functions.

In general, our approach requires to interpolate the Green's
functions for the energies on which the seed photon continuum is
defined. This step is necessary because these energies are typically
defined by the response matrix of the detector with which the data
were taken. This energy grid is outside of our control. We have
calculated grids with a logarithmic energy spacing that is fine enough
that linear interpolation is sufficient to re-grid the Green's
functions in energy.

A second interpolation step in the magnetic field strength is required
to obtain a Green's function for a B-field value not covered by the
grid of precalculated values. Since the Green's functions are
self-similar in E/B we interpolate the Green's functions for
different B in this “energy shifted” system. In a similar way
further interpolations in temperature, output angle, and optical depth
are used to approximate the final spectrum for parameter combinations
off the grid points.

Using MFP interpolation tables, we are able to produce Green's
functions for a large parameter range and with a resolution better
than the resolution of all currently flying X-ray telescopes for the
energy range where cyclotron lines are observed. This takes about
300 CPU hours for one geometry on rather coarse grained magnetic
field and temperature grids. In principle the MFP table mechanism, and
therefore the creation of Green's tables, is reasonable for magnetic
fields 0.01 ≲ B/B_crit≲ 0.12 and
temperatures 0 keV≤ T ≲ 50 keV. For
much smaller B-fields the fully relativistic quantum-electro-dynamic
calculations are unnecessarily complex, while for much higher fields
cyclotron scattering may not be the dominant process <cit.>. Much higher
temperatures are not expected to occur in accreting X-ray
binaries. For the testing of different geometries and other parameters
we currently concentrate on an optimized subset in order to save CPU
time. Table <ref> lists the parameters we have
calculated Green's function tables for. The Green's functions needed
for the evaluation of model spectra for arbitrary parameters within
the precalculated ranges are interpolated from these tables using the
methods given in Appendix <ref>. The Appendix also
describes the extrapolation beyond the covered parameter regime.
However, this extrapolation should be taken with extreme caution.



§ APPLICATION


We now compare
our results to those of previous works and discuss the generation of
synthetic CRSF spectra for a physical continuum. We use end-to-end
comparison to previous works since we do not have access to their
intermediate data products such as mean free path or momentum sampling
information. A comparison to a NuSTAR observation of Cep X-4
is performed to show the applicability of the model to real data and
as a demonstration of the new CRSF fitting model.



 §.§ Code verification



In order to validate the general simulation,
Fig. <ref> compares spectra obtained using our
simulation with earlier calculations by <cit.> for two
different accretion column geometries. Four different angular
regimes are shown. The upper left plot shows the simulated
spectra for a viewing angle almost parallel to the magnetic field
axis. The lower right plot shows spectra for viewing angles
almost perpendicular to the B-field axis. The corresponding angle
bins are defined in terms of the cosine of the viewing angle
μ = cosϑ.

The calculations agree well: the line width decreases with increasing
viewing angle, while the depth of the line decreases with decreasing viewing
angle. Indications for a third harmonic line can be seen in spectra
emerging perpendicular to the magnetic field, while for viewing angles
parallel to the magnetic field only a complex fundamental line is
observed. The strong line wings become less pronounced for
continua with a high energy cutoff, as less photons are spawned from
the higher harmonics <cit.>:

Previous works based on the numerically calculated cyclotron cross
sections obtained by <cit.>, such as <cit.> and
<cit.>, do not show such a good agreement with the
spectra from <cit.> due to their usage of erroneous
integrated cross sections in the simulation code <cit.>.

Some slight deviations remain in the line wings. These are formed by
spawned photons. In the spectra of <cit.> a larger number
of spawned photons escape at small angles to the magnetic field. These
extra photons appear at larger angles in our simulation. The reason
for this disagreement is probably a different angular redistribution
scheme or the approximation of the total
scattering cross section used by <cit.>
for modeling the resonant scattering at harmonics above the first one.
<cit.> use this approximation <cit.>, which is strictly valid only near
the line center, for resonant scattering involving the excitations of
higher harmonics. They resort to a numerical integration over the
scattering angle of the more correct form given in their Eq. 10
<cit.>
for transitions involving only the ground state and the first excited
state.

Although the discussion of cyclotron resonant scattering as a cooling
process for the electrons in the accretion column is beyond the scope
of this paper, we briefly give reasoning for the two different
temperatures used in the comparison. The electron temperature depends
on the photon interactions which in turn depend on the electron
temperature. In the regime where cyclotron resonant scattering is the
dominant cooling process, a convenient choice for an equilibrium
temperature is the temperature at which the energy transfer is zero.
The point of equilibrium depends on the geometry.
Figure <ref> therefore assumes temperatures of
7 keV and 8 keV for the  and 
geometries, respectively, following <cit.>. These
temperatures correspond to the Compton temperature where the energy
transfer is minimal. <cit.> used this technique to show that
the ratio of this equilibrium temperature to the magnetic field is
fairly constant if the temperature is determined solely due to
cyclotron resonant scattering.



 §.§ Cyclotron lines for the <cit.> continuum



Figure <ref> shows synthetic cyclotron lines imprinted on the
continuum from <cit.> for the purpose of illustrating
the influence of the continuum on the CRSFs and as an example
for the cylinder geometry.
The continuum parameters are the same ones as in the theoretical calculation
from <cit.> for Her X-1. This calculation agrees very well
with the BeppoSAX data reported by <cit.>. A spectral fit
of a recent XSPEC implementation of the same continuum model to
NuSTAR data of Her X-1 can be found in the work by
<cit.>.
Figure <ref> shows that the lines become deeper with increasing
optical depth. The line width decreases with increasing viewing angle
to the magnetic field, because of the smaller influence of Lorentz
boosting in this regime. A cylinder geometry Green's table has been
used for imprinting the cyclotron lines on the continuum. This
approach leads to emission like behavior for very small optical depths
and large angles to the magnetic field. Two higher harmonics can
clearly be seen and are especially pronounced for large angles and
relatively high optical depths of τ = 3 ·
10^-3 τ_T.

These are theoretical spectra, which in contrast to observed spectra
correspond to a line forming region with a constant magnetic field
seen from a specific angle. Observations are expected to be smeared
out due to an angle mixing with phase and an extended line forming
region with a magnetic field gradient. Though the accurate handling of
angle mixing requires the inclusion of relativistic effects
<cit.>, its influence on the CRSF line profiles, in the
reference frame of the neutron star, can be roughly estimated by averaging over
multiple viewing angles. Appendix <ref> provides details
on how this can be done easily with the model presented in the
following. All continua used in this work are angle averaged
but the combination of the model with angle dependent continua is
straight forward.



 §.§ The CRSF model and table availability


Together with this work we distribute improved Green's function tables
for the classical geometries, which can be used together with an XSPEC
local model, , to imprint cyclotron lines on arbitrary
continua. See Appendix <ref> for a description of the
model. Table <ref> shows the parameter ranges
covered by these tables. Each table corresponds to one geometry. The
tables have been calculated using thermally and polarization averaged
cross sections.

The  model convolves the Green's table set by the user with the
given continuum. It can, in principle, extrapolate beyond the
ranges provided in Table <ref>, but this should be
used with care and currently triggers a warning message. The model and
preliminary Green's function tables are available online (see link in
footnote, page 1). The currently rather coarse grained parameter grids
will be refined successively and will be made available at the same
location.





 §.§ Fitting Cep X-4


In order to demonstrate the applicability of the model to observed data, a
comparison between empirically fitted line profiles with synthetic ones from
the simulation described above will be performed in the following. The results
further motivate the necessity of a cyclotron line model based on firm physical
grounds for the description of CRSF line profiles.

Cep X-4, also known as GS 2138+56, was discovered in
1972 June and July observations with the OSO-7 X-ray
telescope <cit.>. It is an accreting high mass X-ray binary
(HMXB) with a pulse period of ∼66 s as found
in its 1988 outburst with Ginga
<cit.> and later confirmed by <cit.>. The optical
counter part is a Be star <cit.> at a distance of 3.8 ±
0.6 kpc <cit.>. <cit.> found that the
addition of a cyclotron line at 31 keV improved the fit
but they did not include it in their discussion for the sake of
comparison to other cataloged X-ray pulsars. <cit.> proposed
a cyclotron line at an energy of 30.5 ± 0.4 keV and
deduced a magnetic field of 2.6 × 10^12 (1+z) G
for the source.
RXTE/PCA observed further outbursts in 1997 and 2002.
<cit.> performed a spectral and timing analysis of the
latter outburst and found a cyclotron line at
30.7^+1.8_-1.9 keV. NuSTAR
<cit.> observed Cep X-4 close to the maximum of
the outburst and during its decline on 2014 June 18/19 and 2014 July
1/2, respectively <cit.>. The cyclotron line was measured at
an energy of ∼30 keV in both observations. <cit.>
found that its shape deviates from a simple Gaussian line profile <cit.>,
which makes this observation a good candidate for an application of
the physical CRSF model described in this work.

We re-extracted the NuSTAR data from ObsID 80002016002 (2014
June 18/19, exposure time 40.5 ks) near the maximum of the outburst
using the same settings and
procedure as described by <cit.> but using the CalDB version
20160922 and the NuSTAR data analysis software (NuSTARDAS)
version 1.6.0 as distributed with HEASOFT 6.19. The source and
background spectra for focal plane module A and B (FPMA and FPMB) were
extracted separately, using circular regions with radii of 120.
We use data in the 3.6–55 keV band and rebinned the data to a
signal-to-noise (S/N) of 10 below 45 keV and a S/N of 5 above that.

We use the same empirical continuum model as <cit.>,
which was already found by <cit.> to describe the continuum
well, that is, a power-law with a Fermi-Dirac cutoff <cit.>,

    F(E) = A E^-Γ1/1 + e^-(E - E_cut) / E_fold ,

with normalization constant A, photon index Γ, cutoff energy E_cut,
and folding energy E_fold.

An improved version of the  model, namely [
<http://pulsar.sternwarte.uni-erlangen.de/wilms/research/tbabs/>], is used
with abundances by <cit.> and cross sections by <cit.> in
order to account for photoelectric absorption of the continuum.  A narrow iron
Kα line has been used as in the analysis by <cit.>. Their low
energy black body was dismissed because it did not improve the fit. Also,
contrary to the analysis by , but in agreement with the
continuum model used by <cit.>, we added a “10 keV” feature, that
is, a Gaussian emission line that facilitates the modeling of the continuum in
the ∼8–20 keV range, which has been applied before, for various sources
and instruments <cit.>. We found the width of this broad
emission component to be in agreement with but much better constrained than the
corresponding component used by <cit.>. The centroid energy of 16.1
/ 1.3 = 12.4 keV is below the value of 14.4 keV found by
<cit.>. The additional factor of 1.3 is necessary for the comparison
because of the gravitational redshift: in contrast to previous modeling
approaches, our model includes a redshifting of all components but the iron
line with a fixed redshift of z = 0.3 <cit.>. The
XSPEC convolution model  is used for that purpose.

The luminosity of 1–6 × 10^36 erg s^-1 <cit.>
is low compared to other accreting X-ray pulsars exhibiting cyclotron
lines <cit.> and well below the theoretical limit where a
shock in the accretion column is expected to form <cit.>.
This suggests the usage of a slab geometry for the CRSF model,
, which is the only model that we use to fit the CRSF
line profile, meaning that we do not use any other absorption line
component.

The complete fit model used is thus:
,

where  is the cross-calibration of
focal plane module FPMB relative to FPMA.

Starting with the parameters obtained from the fits by
<cit.> and <cit.>, the physical CRSF model for a
slab 1-0 geometry
converged towards an acceptable fit with a reduced χ^2 of 1.17
for 862 degrees of freedom.
The unfolded spectrum is shown in Fig. <ref>a together with
the model for both detectors, FPMA and FPMB, and the corresponding
residuals are shown in Fig. <ref>b. The best fit parameters can be found in
Table <ref>, with uncertainties given at the 90% confidence
level.





As discussed by <cit.>, the position of the cyclotron line
can be affected by the continuum model. This complicates the
comparison of the magnetic field strength to previous works.
<cit.> noticed an absorption feature at ∼
31 keV using a power-law plus exponential cutoff.
<cit.> used a power-law times cyclotron scattering cutoff
<cit.> to describe the continuum and found the feature to
reside at 30.5 ± 0.4 keV. Using a model consisting of
negative and positive power-laws with a common exponential cutoff
<cit.>, <cit.> detected the cyclotron
line at 28.8 ± 0.4  keV.
<cit.> included a fundamental cyclotron line at 27.5 ± 0.4 keV,
27.7 ± 0.4 keV, or 29.6 ± 0.5 keV in their
best-fit models for an NPEX, CompTT, or Fermi-Dirac plus blackbody continuum model,
respectively, in their analysis of the Suzaku observation of Cep X-4
in 2014. The corresponding widths of the Gaussian absorption
lines used for fitting the CRSFs differ from each other as well between
the different continuum models.
<cit.> and
<cit.>, both using a Fermi-Dirac continuum model, found cyclotron lines at
30.7_-1.9^+1.8 keV and
30.39_-0.14^+0.17 keV (and
29.42_-0.24^+0.27 keV during the decline of the 2014
outburst), respectively.
The analyses of all these works differ by more than the continuum model: some
use additional  model components to model the continuum and/or asymmetries in
the fundamental line profile, different models are used to describe the line
shape including Gaussian absorption lines and pseudo-Lorentzian profiles such
as , different instruments might be responsible for systematic
deviations, and the physical behaviour of the source itself, such as variations
of the height of the line forming region with luminosity <cit.>, might lead to different magnetic field strengths — and
therefore differences of the measured cyclotron line energy — between
observations. The range of the values from the previous works listed above,
from ∼28 keV to ∼31 keV, and differences of
∼8% resulting from the application of different continuum models to
the same observation further illustrate the incomparableness of
cyclotron line energies resulting from differing analyses.

In order to compare the physical line model with an empirical model of the line
shape on the basis of the same continuum, we replaced the 
model by a multiplicative absorption model, namely , at the energy
where such a simple absorption line would be expected for the best-fit magnetic
field strength, that is, E_gabs = B_×
11.57  keV≈ 41.96^+0.35_-0.36 keV in the frame of the neutron
star (i.e., before redshifting)
[Note that the physical cyclotron line model has been
  substituted in place, that is, the  model is still within the
  redshift model and therefore shifted to ∼32.3 keV
  by .].
Leaving this energy and all other parameters frozen while fitting the width and
the depths of the empirical absorption line results in a reduced χ^2 of
1.57 for 876 degrees of freedom. The corresponding residuals are shown in
Fig. <ref>c, which clearly illustrates that the centroid energy of the
Gaussian absorption line model is off the cyclotron line. Only when all parameters of the
 model are left free, can the profile of the cyclotron line be represented
by a Gaussian absorption line with a centroid energy of 40.01_-0.14^+0.15 keV before
redshift (red. χ^2 of 1.17 for 875 d.o.f.).
The value of E = 40.01 keV / (1 + z) = 30.78_-0.11^+0.12 keV after redshift reduction is almost consistent with the value
found by <cit.> for the higher energetic Gaussian absorption line used to model
the asymmetric line shape. The ratio of this model to the best-fit model using
the physical CRSF model is shown in Fig. <ref>d. Evidently, the shapes
of the line models differ significantly as expected.

The model for Cep X-4 presented here makes no claim of
uniqueness. Instead many assumptions are made: other geometries might
fit the spectrum equally well, neglecting bulk velocity becomes
questionable if either the continuum or the cyclotron line are formed
in a region of the accretion column with a significant velocity, and
the usage of an empirical continuum with a “10 keV” feature —
the origin of which is unclear — is questionable per se,
to name just a few. Furthermore, only one viewing angle to the
magnetic field axis is taken into account, which is unrealistic considering
that the data are averaged over pulse phase. The  model
provides the possibility to average over multiple angles in order to
overcome this inaccuracy albeit in an approximative way (see Sect. <ref>
for details). The width of the cyclotron line is strongly affected by
the viewing angle and the temperature. Here, the width is mainly fitted
by the angle to the magnetic field axis. Magnetic field, temperature,
and velocity gradients are neglected, though they might largely influence
the CRSF line width as well. Studies with more complex configurations
of the CRSF medium are needed for estimating their influence quantitatively.
Physical continuum models should be combined with the physical model for the
CRSF line shape presented here. Their combined application to many observations
of diverse sources covering a large parameter space of both continuum and line
profile model parameters — some of these parameters might be tied together
— might help to further constrain the highly degenerate parameter regime.



§ SUMMARY

We have described our Monte Carlo code for the generation of synthetic
cyclotron lines. The simulated line profiles have complex shapes and show a
strong dependency on the viewing angle. We have compared our results to the
work from <cit.> and find an overall good agreement for both, the
 and the  geometry.

In order to show the influence of the continuum, the Her X-1 model
spectrum from <cit.> has been used to generate synthetic spectra for
several optical depths and angles to the magnetic field.

A new XSPEC fitting model, , which works on precalculated
tables storing the response of the Monte Carlo simulation to monoenergetic
photon injections, has been introduced and is available online (see link in
footnote, page 1).

Using this model, which is describing the cyclotron line shape on physical
grounds, we successfully fitted NuSTAR spectra of Cep X-4.
The resulting magnetic field strength, B = 3.6 × 10^12 G,
has been found to differ significantly from fits with a Gaussian absorption
line. The reason for this difference might lie in the theoretically complex
profile <cit.> of the fundamental
cyclotron line. This might lead to a different best fit value for the magnetic
field strength — even for the almost symmetrical and smooth line shapes seen
in observations — if the modeled CRSF shape is taken into account. Further
studies of more complex physical setups, including the exploration of other
geometries and parameter gradients, the inclusion of angular mixing due to
relativistic light bending, and a combination with an equally physical
continuum model are necessary in order to obtain a fully self consistent
spectrum of the accretion column. The simulation code presented here provides
the flexibility to address these challenges. The associated XSPEC model
mechanism allows for making the corresponding results available in an easily
usable and familiar way.

Other applications of the new code include the simulation of observable
<cit.> phase lags at the CRSF energy <cit.>, a
combination with models for relativistic light bending to obtain self
consistent pulse profiles of the phase dependent CRSF behavior (Falkner et al.,
2016, in prep.), comparisons to observational data from other sources, and a
study of the dependence of the CRSF profile on the magnetic field geometry.


This work has been partially funded by the Deutsche Forschungsgemeinschaft
under DFG grant number WI 1860/11-1 and by the Deutsches Zentrum für Luft-
und Raumfahrt under DLR grant numbers 50 OR 1113, 50 OR 1207, 50 OR 1410,
and 50 OR 1411. We also acknowledge the Russian Foundation for Basic Research
grant number 14-02-91345. MTW is supported by the Chief of Naval Research and
by the National Aeronautics and Space Administration Astrophysical Data
Analysis Program. We thank the International Space Science Institute in Bern
for inspiring team meetings. The fruitful discussions within the MAGNET
collaboration also had a very positive impact on this work. The figures in this
work have been produced using the  package by <cit.>.


aa






§ THE XSPEC MODEL

Our XSPEC model uses precalculated Green's function tables to imprint
synthetic CRSFs on arbitrary continua. These Green's tables have
been calculated from the response of a medium to monoenergetic photon
injection on a wide range of parameters. Therefore each geometry has
its own Green's function table, which can be selected for use by


  * setting the environment variable CYCLOFS_TABLE to the
  table's location

  * defining the table location as initialization string in the
   file


Table <ref> lists all model parameters. At least one
optical depth, the magnetic field, the temperature, and the angle or
the angular range and the number of angular points must be set for a
successful model evaluation.





 §.§ Interpolation

The model can be configured to interpolate between parameter values in different ways depending on the
parameter: currently a linear interpolation scheme is utilized for all parameters.
Model spectra are also extrapolated, if the desired parameter
combination is not covered by the Green's table. A linear scheme is used
here, as well.
The model prints out a warning message if
extrapolation is used. The extrapolation method used for the optical depth is
slightly different from the others: For a given optical depth out of the table
range, the model convolves its output iteratively with a suitable optical depth
within the table range. This extrapolation-via-successive-convolution method
turned out to be as accurate as linear extrapolation close to the boundaries
but yields much better results than any other method for extrapolation over
orders of magnitude. This is especially useful for studies in the regime of
high optical depths where calculation time increases significantly. The
accuracy of such an extrapolation over orders of magnitude is questionable,
though, as it depends on the assumption of isotropic angular redistribution,
which is normally not justified, as we have shown before <cit.>.
It is, nevertheless, very useful for studying the influence of increased line
depth on the overall combined model flux and therefore included as default
behavior with a warning message.



 §.§ Angular averaging

The model is designed to calculate spectra for exactly one viewing angle,
defined by its cosine mu= cosϑ, to the
magnetic field axis. It provides the possibility for averaging over a range
of angles to the B-field axis, by returning the mean value of
mu_N spectra between mu_min and mu_max for each
energy bin. If mu_N is set to 0, only one spectrum for the angle
mu is calculated. If it is set 1, one spectrum right in the middle of
the angular range specified by mu_min and mu_max, that is,
for an angle μ = 0.5 (μ_min + μ_max), is returned. A
mu_N value of 2 will result in a spectrum averaged from two points,
namely mu_min and mu_max. For higher values, the additional
points are equally spread over the angular range. Note that the parameter
mu is not used at all if mu_N is larger than zero and should
be frozen during fitting in order to avoid useless iterations. The same applies
to the parameters mu_min and mu_max in the case of
mu_N = 0.

