Issue 
A&A
Volume 601, May 2017



Article Number  A99  
Number of page(s)  11  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201630250  
Published online  10 May 2017 
Cyclotron resonant scattering feature simulations^{⋆}
II. Description of the CRSF simulation process
^{1} Dr. Karl RemeisSternwarte and Erlangen Centre for Astroparticle Physics, Sternwartstrasse 7, 96049 Bamberg, Germany
email: fritz.schwarm@sternwarte.unierlangen.de
^{2} LeibnizInstitut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
^{3} CRESST, Department of Physics, and Center for Space Science and Technology, UMBC, Baltimore, MD 21250, USA
^{4} NASA Goddard Space Flight Center, Code 661, Greenbelt, MD 20771, USA
^{5} Space Science Division, Naval Research Laboratory, Washington, DC 203755352, USA
^{6} Department of Physics & Astronomy, George Mason University, Fairfax, VA 220304444, USA
^{7} European Space Astronomy Centre (ESA/ESAC), Science Operations Department, Villanueva de la Cañada, 28692 Madrid, Spain
^{8} Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
^{9} Faculty of Physics, M. V. Lomonosov Moscow State University, Leninskie Gory, 119991 Moscow, Russia
^{10} Sternberg Astronomical Institute, Moscow M. V. Lomonosov State University, Universitetskij pr., 13, 119992 Moscow, Russia
^{11} Institut für Astronomie und Astrophysik, Universität Tübingen (IAAT), Sand 1, 72076 Tübingen, Germany
^{12} ISDC Data Center for Astrophysics, Université de Genève, chemin d’Écogia 16, 1290 Versoix, Switzerland
Received: 14 December 2016
Accepted: 23 January 2017
Context. Cyclotron resonant scattering features (CRSFs) are formed by scattering of Xray photons off quantized plasma electrons in the strong magnetic field (of the order 10^{12} G) close to the surface of an accreting Xray pulsar. Due to the complex scattering crosssections, 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.
Aims. 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.
Methods. 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 timeconsuming simulations possible on convenient time scales. Furthermore, it can generate responses to monoenergetic photon injections, producing Green’s functions, which can be used later to generate spectra for arbitrary continua.
Results. 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 X4. 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 crosscomparison of various models. The XSPEC model and the Green’s function tables are available online (see link in footnote, page 1).
Key words: Xrays: binaries / stars: neutron / methods: numerical
© ESO, 2017
1. Introduction
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 Xray binaries. Their complex shape (Isenberg et al. 1998; Schönherr et al. 2007; Nishimura 2008; Fürst et al. 2015; Schwarm et al. 2017, and references therein) reveals information about the environment of the lineforming 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 (Schwarm et al. 2017, Paper I in the following), we described the calculation of the relevant crosssection for the range of conditions typically expected for accreting neutron stars in Xray binary systems. In this work, we discuss the application of these crosssections 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 (e.g., Wang et al. 1988, 1989; Araya & Harding 1996, 1999; Schönherr et al. 2007; Nobili et al. 2008a,b), Feautrier (e.g., Meszaros & Nagel 1985; Alexander et al. 1989; Nishimura 2008), and (semi)analytic methods (e.g., Wang et al. 1993). Wang et al. (1988) provide an overview of the methods used in previous works. Isenberg et al. (1998) 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 (Isenberg et al. 1998). From all of the calculations mentioned above, only Schönherr et al. (2007) provided a method for direct comparison between data and observations with the corresponding large grids of parameters.
Araya & Harding (1996) performed a first comparison of simulated MC line profiles to data from the Xray pulsar A0535+26. By convolving precalculated Green’s functions, describing the response of a CRSF medium to monoenergetic photon injection, with arbitrary continua, Schönherr et al. (2007) disentangled the timeconsuming MC simulations from the choice of continuum. This approach enabled us to develop a code that can be used to compare MC simulated spectra to observational data using standard Xray astronomical data analysis packages such as XSPEC (Arnaud 1996) or ISIS (Houck & Denicola 2000). 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 crosssection code by Sina (1996, see Paper I for details).
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 crosssections have been described in Paper I. Here, we 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 Xray 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. 2, 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 number of CPUs and computer systems used. 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. 3. 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 X4 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.
2. 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 that we use in Sect. 3 in order to validate our results with respect to the ones from earlier works using the MC method.
2.1. 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” (Ibrahim et al. 2002, 2003) in the spectra at energies below the electron cyclotron lines, are neglected. There are several other processes altering the properties and paths of Xray photons in the vicinity of magnetic fields near the critical field strength (Canuto & Ventura 1977), (1)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 (Adler 1971). Pair production becomes possible if the photon energy exceeds 2m_{e}c^{2} (Daugherty & Harding 1983). Apart from photonelectron interactions, Sina (1996) 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 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 Bfield are quantized to discrete values corresponding to the Landau energy levels, E_{n} ~ n·12 keV·B/ 10^{12} G. 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 absorptionlike lines at the fundamental Landau energy E_{1} and its integer multiples (Gnedin & Sunyaev 1974; Canuto & Ventura 1977; Araya & Harding 1996; Araya 1997; Schönherr et al. 2007, and references therein).
The cyclotron scattering crosssection strongly depends on the incident photon angle and energy (Paper I and Canuto et al. 1971; Canuto & Ventura 1977; Ventura 1979; Mészáros & Ventura 1979). 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 Bfield. This makes the accurate analytic treatment difficult (Canuto et al. 1971) 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 (Bonazzola et al. 1979). 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 deexcitation to the ground state. Since deexcitation preferentially takes place to the next lower Landau level (Latal 1986), 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.
2.2. Monte Carlo simulation
The advantage of a Monte Carlo simulation is that its photontracing 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. 1 and Sect. 2.3 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 homogeneous density is usually calculated from the optical depth into a given direction and the corresponding extension of the cylinder. For the classical “slab 10” and “slab 11” geometries, the direction parallel to the Bfield axis is chosen, while in the “cylinder” geometry, the optical depth is defined perpendicular to this axis (see Fig. 2). The physical properties inside the accretion column are given by the density, the magnetic field, the electron temperature, and the 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. 2 together with their implied seed emission patterns, for comparison of our results to Isenberg et al. (1998) in Sect. 3.1.
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; and a plane source describes an emitting plane perpendicular to the Bfield axis at a given height. Variants emitting only photons upwards, that is, in the direction of the Bfield 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 deexciting electron to the point where the spawned photons are interacting with or escaping the medium ensure selfconsistent treatment.
The most timeintensive parts of processing are parallelized using the Message Passing Interface (1994, MPI), decreasing the required CPU time efficiently as shown in Fig. 3. Because available computing power has increased significantly since earlier approaches to this problem, we were also able to introduce additional features to the simulation that 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, Bfield gradients, and an inhomogeneous density stratification.
2.3. 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. 1. Seed photons, generated by the configured photon sources, must be processed along with 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.
Fig. 1 Flowchart of the complete Monte Carlo process. The mean free path calculation and the electron momentum sampling block are described in detail in Paper I. Solid lines correspond to photonrelated steps, dashed lines show electronrelated processes. Diamond shapes depict decisions. 
Each photon is propagated according to its current path length, which is sampled from the mean free path tables provided by Schwarm et al. (2017). 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 processes other than resonant cyclotron scattering, which is justified by the comparably small time scales of the electron decay compared to the collisional excitation rate (Bonazzola et al. 1979; Latal 1986; Arons et al. 1987; Araya 1997; Isenberg et al. 1998). Nevertheless, they can be excited to higher Landau levels in the course of resonant cyclotron scattering. The final Landau level for each interaction is sampled by comparing the scattering crosssections 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 crosssection 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 deexcitation 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. Deexcitation 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 1 for initial Landau levels up to n_{i} = 7. The sixth and seventh levels are provided as a reference and in order to allow for estimations to be made 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, , and the corresponding transition to the ground state, , becomes in agreement with the corresponding calculation by Latal (1986, Table 5 therein). The ratio is becoming smaller for larger fields as shown by Latal (1986), 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 deexcitation 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 reemits 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 Schönherr et al. (2014) for the calculation of phasedependent cyclotron line spectra and emission patterns from a cylindrically symmetric accretion column throughout which such a flow with relativistic velocities is expected (see, e.g., Becker & Wolff 2007). 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 X4 the omission of bulk velocity is justified by the chosen geometry; in “slab” geometries, the lineforming region is assumed to reside at the neutron star hot spot at which the matter has been decelerated to nonrelativistic 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 angledependent decay rates. All spawned photons are further propagated in the same way as the seed photons until they leave the medium.
2.4. Classic geometries
Figure 2 depicts various geometries that differ in radius to height ratio and the position and type of the seed photon source.
Cyclotron emission rates for a magnetic field B = 0.1 B_{crit}, averaged over initial and final electron spins and summed over all possible polarization modes, in units of ω_{B} = eB/m (see, e.g., Herold et al. 1982).
Fig. 2 Accretion column geometries used in earlier cyclotron calculations. The slab 10 geometry is a bottomilluminated slab. The slab 11 geometry corresponds to a slab illuminated by a seed photon plane in the middle of the slab. In cylinder geometry, the seed photons are injected isotropically from a line in the center of the cylinder. Spectra emerging from the cylinder geometry are characterized by the optical depth in the radial direction, while the total optical depth along the vertical axis is used to define the scale of slab geometries. 
Slab geometries mimic a thin scattering volume of infinite radius. We approximate this in our numerical simulations by setting the radius to a value that is very large in comparison 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 Bfield 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 slab10/slab11 or cylinder 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 Isenberg et al. (1998), in the slab 10 geometry, the medium is solely illuminated from the bottom, while in the slab 11 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 Bfield axis, that is, the cylinder radius in units of optical depth 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 slab 10 and slab 11 geometry in Fig. 4 and spectra for cylinder geometry in Fig. 5. The CRSF shape is, in general, highly sensitive to the simulated geometry. The formation of line wings, for example, is especially pronounced in slab 11 geometry for viewing angles almost parallel to the magnetic field, which can be seen in Fig. 4.
2.5. CPU time
We reduce the calculation time by using MPI to distribute the work among multiple processors. Figure 3 shows that our simulation is parallelized efficiently. We performed tests on different computer systems and found that using approximately eight CPUs for each simulation provides a convenient compromise between simulation time and efficiency.
Fig. 3 Absolute CPU time (filled symbols) on the left axis and multicore efficiency (lines) on the right axis for the simulation of synthetic spectra. The efficiency is calculated as η = t_{1}/ (nt_{n}), where t_{n} is the execution time on n CPUs. The simulation speed tests have been performed on an AMD Opteron 2.2 GHz system (dashed line) as well as on an Intel Xeon 2690 2.6 GHz system (solid line and filled symbols). The simulations have been performed for the same parameters as in Fig. 4 below, but with only 100 input photons per energy bin instead of 10 000. Absolute CPU times depend on the chosen parameter values. 
2.6. Green’s table approach
Performing a timeconsuming 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 (Schönherr et al. 2007), 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 monoenergetic 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 (Becker & Wolff 2007), can then be obtained by convolving it with the Green’s functions.
In general, our approach requires the interpolation of 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 collected. 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 regrid 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 Bfield value not covered by the grid of precalculated values. Since the Green’s functions are selfsimilar in E/B, we interpolate the Green’s functions for different B in this “energyshifted” 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 Xray telescopes for the energy range where cyclotron lines are observed. This takes approximately 300 CPU hours for one geometry on relatively coarsegrained magneticfield 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 Bfields, the fully relativistic quantumelectrodynamic calculations are unnecessarily complex, while for much higher fields cyclotron scattering may not be the dominant process (see, e.g., Sina 1996). Much higher temperatures are not expected to occur in accreting Xray binaries. For the testing of different geometries and other parameters, we currently concentrate on an optimized subset in order to save CPU time. Table 2 lists the parameters for which we have calculated Green’s function tables. 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 A.1. The Appendix also describes the extrapolation beyond the covered parameter regime. However, this extrapolation should be taken with extreme caution.
3. 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 endtoend 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 X4 is performed to show the applicability of the model to real data and as a demonstration of the new CRSF fitting model.
3.1. Code verification
Parameter ranges covered by CRSF Green’s function tables.
Fig. 4 Comparison of this work (solid lines) to Isenberg et al. (1998, their Fig. 7; dashed lines). Red lines denote the slab 11 geometry while green lines show the slab 10 geometry. The emission in different angular bins is labeled by the cosine of the viewing angle to the magnetic field axis, μ = cosϑ. Temperatures of 7 keV and 8 keV were used for the slab 10 and slab 11 geometries, respectively. 
In order to validate the general simulation, Fig. 4 compares spectra obtained using our simulation with earlier calculations by Isenberg et al. (1998) for two different accretioncolumn 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 Bfield 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 highenergy cutoff, as less photons are spawned from the higher harmonics (Isenberg et al. 1998; Araya 1997; Schönherr et al. 2007, and references therein):
Previous works based on the numerically calculated cyclotron crosssections obtained by Sina (1996), such as Araya (1997) and Schönherr et al. (2007), do not show such a good agreement with the spectra from Isenberg et al. (1998) due to their usage of erroneous integrated crosssections in the simulation code (Schwarm et al. 2012, and Paper I).
Some slight deviations remain in the line wings. These are formed by spawned photons. In the spectra of Isenberg et al. (1998), 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 crosssection used by Isenberg et al. (1998, Eq. (11) therein) for modeling the resonant scattering at energies above the first harmonic line. Isenberg et al. (1998) use this approximation (see also Daugherty & Ventura 1977; Fenimore et al. 1988), only strictly valid 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) (see also Canuto et al. 1971; Herold 1979; Ventura 1979; Wasserman & Salpeter 1980; Harding & Daugherty 1991; Graziani 1993) 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 4 therefore assumes temperatures of 7 keV and 8 keV for the slab 10 and slab 11 geometries, respectively, following Isenberg et al. (1998). These temperatures correspond to the Compton temperature where the energy transfer is minimal. Lamb et al. (1990) 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.
3.2. Cyclotron lines for the Becker & Wolff (2007) continuum
Fig. 5 Continuum from Becker & Wolff (2007) with imprinted cyclotron lines for several angles μ = cosϑ to the magnetic field axis and optical depths τ. The continuum parameters are the same as in the spectrum from Becker & Wolff (2007, Fig. 6 therein) for Her X1. The different colors show three angles to the magnetic field axis. The optical depth is increasing from the bottom to the top as shown by the axis on the righthand side. 
Figure 5 shows synthetic cyclotron lines imprinted on the continuum from Becker & Wolff (2007) for the purpose of illustrating the influence of the continuum on the CRSFs and as an example of the cylinder geometry. The continuum parameters are the same ones as in the theoretical calculation from Becker & Wolff (2007) for Her X1. This calculation agrees very well with the BeppoSAX data reported by dal Fiume et al. (1998). A spectral fit of a recent XSPEC implementation of the same continuum model to NuSTAR data of Her X1 can be found in the work by Wolff et al. (2016). Figure 5 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 emissionlike 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 lineforming 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 lineforming region with a magnetic field gradient. Though the accurate handling of angle mixing requires the inclusion of relativistic effects (Schönherr et al. 2014), 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 A.2 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 angledependent continua is straight forward.
3.3. 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, cyclofs, to imprint cyclotron lines on arbitrary continua. See Appendix A for a description of the model. Table 2 shows the parameter ranges covered by these tables. Each table corresponds to one geometry. The tables have been calculated using both thermally averaged and polarizationaveraged crosssections.
The cyclofs 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 2, 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.
3.4. Fitting Cep X4
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 X4, also known as GS 2138+56, was discovered in June and July observations in 1972 with the OSO7 Xray telescope (Ulmer et al. 1973). It is an accreting highmass Xray binary (HMXB) with a pulse period of ~66 s as found in its 1988 outburst with Ginga (Makino & GINGA Team 1988a,b) and later confirmed by Koyama et al. (1991). The optical counterpart is a Be star (Roche et al. 1997) at a distance of 3.8 ± 0.6 kpc (BonnetBidaud & Mouchet 1998). Koyama et al. (1991) 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 Xray pulsars. Mihara et al. (1991) 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. McBride et al. (2007) performed a spectral and timing analysis of the latter outburst and found a cyclotron line at . NuSTAR (Harrison et al. 2013) observed Cep X4 close to the maximum of the outburst and during its decline on June 18/19, 2014, and July 1/2, 2014, respectively (Fürst et al. 2015). The cyclotron line was measured at an energy of ~30 keV in both observations. Fürst et al. (2015) found that its shape deviates from a simple Gaussian line profile (see also Jaisawal & Naik 2015), which makes this observation a good candidate for an application of the physical CRSF model described in this work.
We reextracted the NuSTAR data from ObsID 80002016002 (June 18/19, 2014, exposure time 40.5 ks) near the maximum of the outburst using the same settings and procedure as described by Fürst et al. (2015) 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 modules A and B (FPMA and FPMB) were extracted separately, using circular regions with radii of 120′′. We used data in the 3.6–55 keV band and rebinned the data to a signaltonoise ratio (S/N) of 10 below 45 keV and a S/N of 5 above that.
We use the same empirical continuum model as Fürst et al. (2015), which was already found by McBride et al. (2007) to describe the continuum well, that is, a powerlaw with a FermiDirac cutoff (Tanaka 1986; Makishima & Mihara 1992), (2)with normalization constant A, photon index Γ, cutoff energy E_{cut}, and folding energy E_{fold}.
An improved version of the tbabs model, namely tbnew^{1}, is used with abundances by Wilms et al. (2000) and crosssections by Verner et al. (1996) in order to account for photoelectric absorption of the continuum. A narrow iron Kα line has been used as in the analysis by Fürst et al. (2015). Their lowenergy black body was dismissed because it did not improve the fit. Also, contrary to the analysis by Fürst et al., but in agreement with the continuum model used by McBride et al. (2007), 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 (Mihara 1995; Coburn et al. 2002). We found the width of this broad emission component to be in agreement with, but much better constrained than, the corresponding component used by McBride et al. (2007). The centroid energy of 16.1/1.3 = 12.4 keV is below the value of 14.4 keV found by McBride et al. (2007). 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 apart from the iron line with a fixed redshift of z = 0.3 (see, e.g., Schönherr et al. 2007). The XSPEC convolution model zashift is used for that purpose.
The luminosity of 1–6 × 10^{36} erg s^{1} (Fürst et al. 2015) is low compared to other accreting Xray pulsars exhibiting cyclotron lines (Fürst et al. 2015) and well below the theoretical limit where a shock in the accretion column is expected to form (Becker et al. 2012). This suggests the usage of a slab geometry for the CRSF model, cyclofs, 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: C_{FPMA}× tbabs × (redshift(cyclofs(powerlaw × fdcut + gauss_{10 keV})) + gauss_{iron}), where C_{FPMA} is the crosscalibration of focal plane module FPMB relative to FPMA.
Starting with the parameters obtained from the fits by McBride et al. (2007) and Fürst et al. (2015), the physical CRSF model for a slab 10 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. 6a together with the model for both detectors, FPMA and FPMB, and the corresponding residuals are shown in Fig. 6b. The best fit parameters can be found in Table 3, with uncertainties given at the 90% confidence level.
Fig. 6 a) Unfolded spectrum, bestfit model, and individual continuum components for an observation (ObsID 80002016002) of Cep X4. The data and residuals from NuSTAR FPMA and FPMB are shown in red and blue, respectively. b) Residuals to the bestfit model. c) Residuals to the bestfit model with a Gaussian absorption line centered at the cyclotron line position found by the bestfit model. Only the width and depth of the line have been fitted. d) Ratio between the physical and the empirical model. In the empirical model, all line parameters were allowed to vary including the centroid energy. The data have been regrouped for improved clarity. 
Bestfit parameters for the cyclofs model in combination with a FermiDirac cutoff continuum model.
As discussed by Müller et al. (2013), 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. Koyama et al. (1991) noticed an absorption feature at ~31 keV using a powerlaw plus exponential cutoff. Mihara et al. (1991) used a powerlaw times cyclotron scattering cutoff (Tanaka 1986) to describe the continuum and found the feature to reside at 30.5 ± 0.4 keV. Using a model consisting of negative and positive powerlaws with a common exponential cutoff (Mihara 1995, NPEX), Makishima et al. (1999) detected the cyclotron line at 28.8 ± 0.4 keV. Jaisawal & Naik (2015) included a fundamental cyclotron line at 27.5 ± 0.4 keV, 27.7 ± 0.4 keV, or 29.6 ± 0.5 keV in their bestfit models for an NPEX, CompTT, or FermiDirac plus black body continuum model, respectively, in their analysis of the Suzaku observation of Cep X4 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. McBride et al. (2007) and Fürst et al. (2015), both using a FermiDirac continuum model, found cyclotron lines at and (and during the decline of the 2014 outburst), respectively. The analyses of all these works differ by more than the continuum model; some use additional gabs 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 pseudoLorentzian profiles such as cyclabs; different instruments might be responsible for systematic deviations; and the physical behaviour of the source itself, such as variations of the height of the lineforming region with luminosity (see, e.g., Becker et al. 2012), 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 cyclofs model by a multiplicative absorption model, namely gabs, at the energy where such a simple absorption line would be expected for the bestfit magnetic field strength, that is, 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 gabs model is still within the redshift model and therefore shifted to ~32.3 keV by zashift. 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. 6c, which clearly illustrates that the centroid energy of the Gaussian absorption line model is off the cyclotron line. Only when all parameters of the gabs model are left free, can the profile of the cyclotron line be represented by a Gaussian absorption line with a centroid energy of before redshift (red. χ^{2} of 1.17 for 875 d.o.f.). The value of after redshift reduction is almost consistent with the value found by Fürst et al. (2015) for the higher energetic Gaussian absorption line used to model the asymmetric line shape. The ratio of this model to the bestfit model using the physical CRSF model is shown in Fig. 6d. Evidently, the shapes of the line models differ significantly as expected.
The model for Cep X4 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 cyclofs model provides the possibility to average over multiple angles in order to overcome this inaccuracy, albeit in an approximative way (see Sect. A.2 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 which may be tied together, might help to further constrain the highly degenerate parameter regime.
4. 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 Isenberg et al. (1998) and find a good agreement overall for both the slab 10 and the slab 11 geometry.
In order to show the influence of the continuum, the Her X1 model spectrum from Becker & Wolff (2007) has been used to generate synthetic spectra for several optical depths and angles to the magnetic field.
A new XSPEC fitting model, cyclofs, 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 describes the cyclotron line shape on physical grounds, we successfully fitted NuSTAR spectra of Cep X4. 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 (Isenberg et al. 1998; Schönherr et al. 2007; Nishimura 2008) 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 selfconsistent spectrum of the accretion column. The simulation code presented here provides the flexibility to address these challenges. The associated XSPEC model mechanism allows the corresponding results to be made available in an easily usable and familiar way.
Other applications of the new code include the simulation of observable (Ferrigno et al. 2011) phase lags at the CRSF energy (Schönherr et al. 2014), a combination with models for relativistic light bending to obtain selfconsistent pulse profiles of the phase dependent CRSF behavior (Falkner et al., in prep.), comparisons to observational data from other sources, and a study of the dependence of the CRSF profile on the magnetic field geometry.
Acknowledgments
This work has been partially funded by the Deutsche Forschungsgemeinschaft under DFG grant number WI 1860/111 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 140291345. M.T.W. 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 slxfig package by Davis (2014).
References
 Adler, S. L. 1971, Annals of Physics, 67, 599 [NASA ADS] [CrossRef] [Google Scholar]
 Alexander, S. G., Meszaros, P., & Bussard, R. W. 1989, ApJ, 342, 928 [NASA ADS] [CrossRef] [Google Scholar]
 Araya, R. A. 1997, Ph.D. Thesis, Johns Hopkins University, Baltimore, USA [Google Scholar]
 Araya, R. A., & Harding, A. K. 1996, A&AS, 120, C183 [NASA ADS] [Google Scholar]
 Araya, R. A., & Harding, A. K. 1999, ApJ, 517, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes (San Francisco: ASP), ASP Conf. Ser., 101, 17 [Google Scholar]
 Arons, J., Klein, R. I., & Lea, S. M. 1987, ApJ, 312, 666 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, P. A., & Wolff, M. T. 2007, ApJ, 654, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, P. A., Klochkov, D., Schönherr, G., et al. 2012, A&A, 544, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonazzola, S., Heyvaerts, J., & Puget, J. L. 1979, A&A, 78, 53 [NASA ADS] [Google Scholar]
 BonnetBidaud, J. M., & Mouchet, M. 1998, A&A, 332, L9 [NASA ADS] [Google Scholar]
 Canuto, V., & Ventura, J. 1977, Fund. Cosmic Phys., 2, 203 [Google Scholar]
 Canuto, V., Lodenquai, J., & Ruderman, M. 1971, Phys. Rev. D, 3, 2303 [NASA ADS] [CrossRef] [Google Scholar]
 Coburn, W., Heindl, W. A., Rothschild, R. E., et al. 2002, ApJ, 580, 394 [NASA ADS] [CrossRef] [Google Scholar]
 dal Fiume, D., Orlandini, M., Cusumano, G., et al. 1998, A&A, 329, L41 [NASA ADS] [Google Scholar]
 Daugherty, J. K., & Harding, A. K. 1983, ApJ, 273, 761 [Google Scholar]
 Daugherty, J. K., & Ventura, J. 1977, A&A, 61, 723 [NASA ADS] [Google Scholar]
 Davis, J. E. 2014, http://www.jedsoft.org/fun/slxfig/ [Google Scholar]
 Fenimore, E. E., Conner, J. P., Epstein, R. I., et al. 1988, ApJ, 335, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrigno, C., Falanga, M., Bozzo, E., et al. 2011, A&A, 532, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fürst, F., Pottschmidt, K., Miyasaka, H., et al. 2015, ApJ, 806, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, I. N., & Sunyaev, R. A. 1974, A&A, 36, 379 [NASA ADS] [Google Scholar]
 Graziani, C. 1993, ApJ, 412, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Harding, A. K., & Daugherty, J. K. 1991, ApJ, 374, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Herold, H. 1979, Phys. Rev. D, 19, 2868 [NASA ADS] [CrossRef] [Google Scholar]
 Herold, H., Ruder, H., & Wunner, G. 1982, A&A, 115, 90 [NASA ADS] [Google Scholar]
 Houck, J. C., & Denicola, L. A. 2000, in Astronomical Data Analysis Software and Systems IX, eds. N. Manset, C. Veillet, & D. Crabtree (San Francisco: ASP), ASP Conf. Ser., 216, 591 [Google Scholar]
 Ibrahim, A. I., SafiHarb, S., Swank, J. H., et al. 2002, ApJ, 574, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Ibrahim, A. I., Swank, J. H., & Parke, W. 2003, ApJ, 584, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Isenberg, M., Lamb, D. Q., & Wang, J. C. L. 1998, ApJ, 505, 688 [NASA ADS] [CrossRef] [Google Scholar]
 Jaisawal, G. K., & Naik, S. 2015, MNRAS, 453, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K., Kawada, M., Tawara, Y., et al. 1991, ApJ, 366, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Lamb, D. Q., Wang, J. C. L., & Wasserman, I. M. 1990, ApJ, 363, 670 [NASA ADS] [CrossRef] [Google Scholar]
 Latal, H. G. 1986, ApJ, 309, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, F., & GINGA Team 1988a, IAU Circ., 4575 [Google Scholar]
 Makino, F., & GINGA Team 1988b, IAU Circ., 4577 [Google Scholar]
 Makishima, K., & Mihara, T. 1992, in Frontiers Science Series, eds. Y. Tanaka, & K. Koyama, 23 [Google Scholar]
 Makishima, K., Mihara, T., Nagase, F., & Tanaka, Y. 1999, ApJ, 525, 978 [NASA ADS] [CrossRef] [Google Scholar]
 McBride, V. A., Wilms, J., Kreykenbohm, I., et al. 2007, A&A, 470, 1065 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Message Passing Interface 1994, MPI: A MessagePassing Interface Standard, Tech. Rep. FUTCS94230, University of Tennessee, Knoxville [Google Scholar]
 Meszaros, P., & Nagel, W. 1985, ApJ, 298, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Mészáros, P., & Ventura, J. 1979, Phys. Rev. D, 19, 3565 [NASA ADS] [CrossRef] [Google Scholar]
 Mihara, T. 1995, Ph.D. Thesis, Dept. of Physics, Univ. of Tokyo [Google Scholar]
 Mihara, T., Makishima, K., Kamijo, S., et al. 1991, ApJ, 379, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, S., Ferrigno, C., Kühnel, M., et al. 2013, A&A, 551, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nishimura, O. 2008, ApJ, 672, 1127 [NASA ADS] [CrossRef] [Google Scholar]
 Nobili, L., Turolla, R., & Zane, S. 2008a, MNRAS, 386, 1527 [NASA ADS] [CrossRef] [Google Scholar]
 Nobili, L., Turolla, R., & Zane, S. 2008b, MNRAS, 389, 989 [NASA ADS] [CrossRef] [Google Scholar]
 Roche, P., Green, L., & Hoenig, M. 1997, IAU Circ., 6698 [Google Scholar]
 Schönherr, G., Wilms, J., Kretschmar, P., et al. 2007, A&A, 472, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schönherr, G., Schwarm, F.W., Falkner, S., et al. 2014, A&A, 564, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarm, F., Schönherr, G., Wilms, J., & Kretschmar, P. 2012, PoS, INTEGRAL 2012, 153 [Google Scholar]
 Schwarm, F.W., Schönherr, G., Falkner, S., et al. 2017, A&A, 597, A3 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sina, R. 1996, Ph.D. Thesis, University of Maryland, College Park, MD [Google Scholar]
 Tanaka, Y. 1986, in Radiation Hydrodynamics in Stars and Compact Objects, eds. D. Mihalas, & K.H. A. Winkler (Berlin: Springer Verlag), IAU Colloq., 89, Lect. Notes Phys., 255, 198 [Google Scholar]
 Ulmer, M. P., Baity, W. A., Wheaton, W. A., & Peterson, L. E. 1973, ApJ, 184, L117 [NASA ADS] [CrossRef] [Google Scholar]
 Ventura, J. 1979, Phys. Rev. D, 19, 1684 [NASA ADS] [CrossRef] [Google Scholar]
 Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, J. C. L., Wasserman, I. M., & Salpeter, E. E. 1988, ApJS, 68, 735 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, J. C. L., Lamb, D. Q., Loredo, T. J., Wasserman, I. M., & Salpeter, E. E. 1989, Phys. Rev. Lett., 63, 1550 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wang, J. C. L., Wasserman, I., & Lamb, D. Q. 1993, ApJ, 414, 815 [NASA ADS] [CrossRef] [Google Scholar]
 Wasserman, I., & Salpeter, E. 1980, ApJ, 241, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Wolff, M. T., Becker, P. A., Gottlieb, A. M., et al. 2016, ApJ, 831, 194 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: 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 either setting the environment variable CYCLOFS_TABLE to the table’s location or by defining the table location as initialization string in the model.dat file.
Table A.1 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.
Appendix A.1: 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 extrapolationviasuccessiveconvolution 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 (Schönherr et al. 2014).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.
Appendix A.2: 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 Bfield 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 mu_N 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.
XSPEC model parameters.
All Tables
Cyclotron emission rates for a magnetic field B = 0.1 B_{crit}, averaged over initial and final electron spins and summed over all possible polarization modes, in units of ω_{B} = eB/m (see, e.g., Herold et al. 1982).
Bestfit parameters for the cyclofs model in combination with a FermiDirac cutoff continuum model.
All Figures
Fig. 1 Flowchart of the complete Monte Carlo process. The mean free path calculation and the electron momentum sampling block are described in detail in Paper I. Solid lines correspond to photonrelated steps, dashed lines show electronrelated processes. Diamond shapes depict decisions. 

In the text 
Fig. 2 Accretion column geometries used in earlier cyclotron calculations. The slab 10 geometry is a bottomilluminated slab. The slab 11 geometry corresponds to a slab illuminated by a seed photon plane in the middle of the slab. In cylinder geometry, the seed photons are injected isotropically from a line in the center of the cylinder. Spectra emerging from the cylinder geometry are characterized by the optical depth in the radial direction, while the total optical depth along the vertical axis is used to define the scale of slab geometries. 

In the text 
Fig. 3 Absolute CPU time (filled symbols) on the left axis and multicore efficiency (lines) on the right axis for the simulation of synthetic spectra. The efficiency is calculated as η = t_{1}/ (nt_{n}), where t_{n} is the execution time on n CPUs. The simulation speed tests have been performed on an AMD Opteron 2.2 GHz system (dashed line) as well as on an Intel Xeon 2690 2.6 GHz system (solid line and filled symbols). The simulations have been performed for the same parameters as in Fig. 4 below, but with only 100 input photons per energy bin instead of 10 000. Absolute CPU times depend on the chosen parameter values. 

In the text 
Fig. 4 Comparison of this work (solid lines) to Isenberg et al. (1998, their Fig. 7; dashed lines). Red lines denote the slab 11 geometry while green lines show the slab 10 geometry. The emission in different angular bins is labeled by the cosine of the viewing angle to the magnetic field axis, μ = cosϑ. Temperatures of 7 keV and 8 keV were used for the slab 10 and slab 11 geometries, respectively. 

In the text 
Fig. 5 Continuum from Becker & Wolff (2007) with imprinted cyclotron lines for several angles μ = cosϑ to the magnetic field axis and optical depths τ. The continuum parameters are the same as in the spectrum from Becker & Wolff (2007, Fig. 6 therein) for Her X1. The different colors show three angles to the magnetic field axis. The optical depth is increasing from the bottom to the top as shown by the axis on the righthand side. 

In the text 
Fig. 6 a) Unfolded spectrum, bestfit model, and individual continuum components for an observation (ObsID 80002016002) of Cep X4. The data and residuals from NuSTAR FPMA and FPMB are shown in red and blue, respectively. b) Residuals to the bestfit model. c) Residuals to the bestfit model with a Gaussian absorption line centered at the cyclotron line position found by the bestfit model. Only the width and depth of the line have been fitted. d) Ratio between the physical and the empirical model. In the empirical model, all line parameters were allowed to vary including the centroid energy. The data have been regrouped for improved clarity. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.