Cyclotron resonant scattering feature simulations. II. Description of the CRSF simulation process

TL;DR

This paper introduces a new, versatile Monte Carlo simulation code for modeling cyclotron resonant scattering features in X-ray pulsars, enabling complex geometry investigations and improved data interpretation.

Contribution

The paper presents a novel, flexible simulation tool for CRSF modeling that overcomes previous geometric limitations and facilitates detailed spectral analysis.

Findings

The code successfully simulates complex CRSF geometries.

It produces Green's functions for arbitrary continuum spectra.

Application to NuSTAR data demonstrates effective spectral fitting.

Abstract

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. 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 at http://www.sternwarte.uni-erlangen.de/research/cyclo .

Tables

Table 1. Table 1: Cyclotron emission rates for a magnetic field B = 0.1 ​ B crit 𝐵 0.1 subscript 𝐵 crit B=0.1\,B_{\mathrm{crit}} , averaged over initial and final electron spins and summed over all possible polarization modes, in units of ω B = e ​ B / m subscript 𝜔 B 𝑒 𝐵 𝑚 \omega_{\mathrm{B}}=eB/m (see, e.g., Herold et al. , 1982 ) . The rates are given for transitions from the initial Landau level n i ≤ 7 subscript 𝑛 𝑖 7 n_{i}\leq 7 (first column) to all final Landau levels n f subscript 𝑛 𝑓 n_{f} for which n f < n i − 1 subscript 𝑛 𝑓 subscript 𝑛 𝑖 1 n_{f}<n_{i}-1 (first row).

$n$	0	1	2	3	4	5	6
1	$7.8708\times {10}^{-4}$
2	$2.4432\times {10}^{-4}$	$1.7366\times {10}^{-3}$
3	$1.1313\times {10}^{-4}$	$5.8499\times {10}^{-4}$	$2.2486\times {10}^{-3}$
4	$6.3951\times {10}^{-5}$	$2.8415\times {10}^{-4}$	$8.6338\times {10}^{-4}$	$2.5362\times {10}^{-3}$
5	$4.0800\times {10}^{-5}$	$1.6579\times {10}^{-4}$	$4.4945\times {10}^{-4}$	$1.0698\times {10}^{-3}$	$2.6992\times {10}^{-3}$
6	$2.8234\times {10}^{-5}$	$1.0818\times {10}^{-4}$	$2.7355\times {10}^{-4}$	$5.8940\times {10}^{-4}$	$1.2188\times {10}^{-3}$	$2.7889\times {10}^{-3}$
7	$2.0710\times {10}^{-5}$	$7.6123\times {10}^{-5}$	$1.8357\times {10}^{-4}$	$3.7232\times {10}^{-4}$	$7.0200\times {10}^{-4}$	$1.3252\times {10}^{-3}$	$2.8336\times {10}^{-3}$

Table 2. Table 2: Parameter ranges covered by CRSF Green’s function tables

Parameter	Minimum	Maximum
Magnetic field [$B_{\mathrm{crit}}$]	0.01	0.12
Electron temperature $k_{B}T$ [keV]	03	15
Cosine to $B$-field, $\mu =\cos \vartheta$	0.001	0.999
Optical depth [${\tau }_T$]	$1\cdot {10}^{-4}$	$3\cdot {10}^{-3}$

Table 3. Table 3: Best-fit parameters for the cyclofs model in combination with a Fermi-Dirac cutoff continuum model.

Parameter	Value	Unit
$C_{\mathrm{FPMB}}$	$1.0341\pm 0.0020$
$N_H$	${2.1}_{-0.5}^{+0.6}$	${10}^{22}$ cm-2
$A_{\mathrm{cont}}$a𝑎aa𝑎a$a$In photons keV-1s-1cm-2 at 1 keV	${0.230}_{-0.025}^{+0.036}$
$\Gamma$	${1.30}_{-0.07}^{+0.09}$
$E_{\mathrm{cut}}$	${29.1}_{-1.9}^{+2.2}$	keV
$E_{\mathrm{fold}}$	${9.8}_{-0.4}^{+0.6}$	keV
$A_{10keV}$a𝑎aa𝑎a$a$In photons keV-1s-1cm-2 at 1 keV	${0.025}_{-0.006}^{+0.009}$
$E_{10\mathrm{keV}}$	${16.1}_{-0.7}^{+0.5}$	keV
${\sigma }_{10\mathrm{keV}}$	${6.0}_{-0.6}^{+0.7}$	keV
${\tau }_{\parallel }$	$\left({2.59}_{-0.24}^{+0.16}\right)\times {10}^{-4}$	${\tau }_{\mathrm{Th}}$
B	${3.627}_{-0.031}^{+0.030}$	${10}^{12}G$
Tb𝑏bb𝑏b$b$The electron temperature in this fit is only slightly below the current upper boundary of $15\mathrm{keV}$ covered by the table. The corresponding upper error is therefore given at the $68\%$ confidence level.	${14.4}_{-0.9}^{+0.5}$	keV
$\mu =\cos \vartheta$	${0.82}_{-0.05}^{+0.04}$
$A_{\mathrm{iron}}$a𝑎aa𝑎a$a$In photons keV-1s-1cm-2 at 1 keV	$\left({9.6}_{-1.0}^{+1.2}\right)\times {10}^{-4}$
$E_{\mathrm{iron}}$	${6.491}_{-0.028}^{+0.027}$	keV
${\sigma }_{\mathrm{iron}}$	$0.30\pm 0.05$	keV
red. ${\chi }^2$	1.17
d.o.f.	862

Table 4. Table 4: XSPEC model parameters.

Name	Description
tau_para	Parallel optical depth [${\tau }_T$]
tau_perp	Perpendicular optical depth [${\tau }_T$]
B	Magnetic field [${10}^{12}$G]
T	Electron temperature $k_{B}T$ [keV]
mu	Cosine of the viewing angle to the $B$-field, $\mu =\cos \vartheta$
mu_N	Number of mu-values to average over (0 for exact value of mu)
mu_min	Minimum value of mu for averaging
mu_max	Maximum value of mu for averaging

Equations

$$B_{\mathrm{crit}}=\frac{m_{\mathrm{e}}^2{c}^3}{e\hslash }=4.413\times 10^{13}\mathrm{G}.$$

$$F(E)=A{E}^{-\Gamma }\frac{1}{1+e^{-(E-E_{\mathrm{cut}})/E_{\mathrm{fold}}}},$$

Full text

*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 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.

*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 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).

1 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\gtrsim 10^{12}$ G), for example, in the accretion column of the magnetized neutron star in accreting X-ray 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 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 (Schwarm et al., 2017, paper I in the following), 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 (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 over 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}\tau_{\mathrm{T}}$, where $\tau_{\mathrm{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 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, Schönherr et al. (2007) 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 (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 cross section 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 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. 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 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. 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 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.

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 which we will be using 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 X-ray photons in the vicinity of magnetic fields near the critical field strength (Canuto & Ventura, 1977),

[TABLE]

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_{\mathrm{e}}c^{2}$ (Daugherty & Harding, 1983). Apart from photon-electron 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 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, $E_{n}\sim n\cdot 12\,\mathrm{keV}\cdot 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 absorption-like 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 cross section 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 $\vartheta_{\mathrm{in}}$ is measured with respect to the magnetic field axis and mostly specified by its cosine, $\mu_{\mathrm{in}}=\cos\vartheta_{\mathrm{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 (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 de-excitation to the ground state. Since de-excitation 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 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. 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 physical properties inside are given by its homogeneous density111The density is usually calculated from the optical depth into a given direction and the corresponding extension of the cylinder. For the classical “slab 1-0” and “slab 1-1” 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. 2), 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. 2 together with their implied seed emission patterns, for comparison of our results to Isenberg et al. (1998) in Sec. 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, 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 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 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.

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 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 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 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 (Bonazzola et al., 1979; Latal, 1986; Arons et al., 1987; Araya, 1997; Isenberg et al., 1998). 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,\mathrm{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_{\mathrm{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 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\rightarrow 2$. The ratio between the rate for a transition from Landau level $n_{i}=2$ to Landau level $n_{f}=1$, that is, $\overline{\Gamma}_{21}\approx 1.74\times 10^{-3}$, and the corresponding transition to the ground state, $\overline{\Gamma}_{20}\approx 2.44\times 10^{-4}$, becomes $\overline{\Gamma}_{21}/\overline{\Gamma}_{20}\approx 7.11$ 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 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 Schönherr et al. (2014) 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 (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 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.

2.4 Classic geometries

Figure 2 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 slab1-0/slab1-1 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 1-0 geometry the medium is solely illuminated from the bottom, while in the slab 1-1 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 slab 1-0 and slab 1-1 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 1-1 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 $\sim$8 CPUs for each simulation provides a convenient compromise between simulation time and efficiency.

2.6 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 (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 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 (Becker & Wolff, 2007), 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\lesssim B/B_{\mathrm{crit}}\lesssim 0.12$ and temperatures $0\,\mathrm{keV}\leq T\lesssim 50\,\mathrm{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 (see, e.g., Sina, 1996). 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 2 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 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 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.

3.1 Code verification

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 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 $\mu=\cos\vartheta$.

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 (Isenberg et al., 1998; Araya, 1997; Schönherr et al., 2007, and references therein):

Previous works based on the numerically calculated cyclotron cross sections 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 cross sections 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 cross section used by Isenberg et al. (1998, Eq. 11 therein) for modeling the resonant scattering at harmonics above the first one. Isenberg et al. (1998) use this approximation (see also Daugherty & Ventura, 1977; Fenimore et al., 1988), 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 (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 1-0 and slab 1-1 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

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 for the cylinder geometry. The continuum parameters are the same ones as in the theoretical calculation from Becker & Wolff (2007) for Her X-1. 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 X-1 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 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 $\tau=3\cdot 10^{-3}\,\tau_{\mathrm{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 (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 angle dependent 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 thermally and polarization averaged cross sections.

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 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 (Ulmer et al., 1973). It is an accreting high mass X-ray binary (HMXB) with a pulse period of $\sim 66\,\mathrm{s}$ as found in its 1988 outburst with Ginga (Makino & GINGA Team, 1988a, b) and later confirmed by Koyama et al. (1991). The optical counter part is a Be star (Roche et al., 1997) at a distance of $3.8\pm 0.6\,\mathrm{kpc}$ (Bonnet-Bidaud & Mouchet, 1998). Koyama et al. (1991) found that the addition of a cyclotron line at $31\,\mathrm{keV}$ improved the fit but they did not include it in their discussion for the sake of comparison to other cataloged X-ray pulsars. Mihara et al. (1991) proposed a cyclotron line at an energy of $30.5\,\pm\,0.4\,\mathrm{keV}$ and deduced a magnetic field of $2.6\times 10^{12}\,(1+z)\,\mathrm{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 $30.7^{+1.8}_{-1.9}\,\mathrm{keV}$. NuSTAR (Harrison et al., 2013) 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 (Fürst et al., 2015). The cyclotron line was measured at an energy of $\sim 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 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 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 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 Fürst et al. (2015), which was already found by McBride et al. (2007) to describe the continuum well, that is, a power-law with a Fermi-Dirac cutoff (Tanaka, 1986; Makishima & Mihara, 1992),

[TABLE]

with normalization constant $A$, photon index $\Gamma$, cutoff energy $E_{\mathrm{cut}}$, and folding energy $E_{\mathrm{fold}}$.

An improved version of the tbabs model, namely tbnew222 http://pulsar.sternwarte.uni-erlangen.de/wilms/research/tbabs/, is used with abundances by Wilms et al. (2000) and cross sections by Verner et al. (1996) in order to account for photoelectric absorption of the continuum. A narrow iron K$\alpha$ line has been used as in the analysis by Fürst et al. (2015). Their low energy 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 $\sim$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\,\mathrm{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 but 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\times 10^{36}\,\mathrm{erg}\,\mathrm{s}^{-1}$ (Fürst et al., 2015) is low compared to other accreting X-ray 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: CFPMA $\times$ tbabs $\times$ (redshift(cyclofs(powerlaw $\times$ fdcut + gauss${}_{\mathrm{10\,keV}}$)) + gaussiron), where CFPMA is the cross-calibration 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 1-0 geometry converged towards an acceptable fit with a reduced $\chi^{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.

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 $\sim 31\,\mathrm{keV}$ using a power-law plus exponential cutoff. Mihara et al. (1991) used a power-law times cyclotron scattering cutoff (Tanaka, 1986) to describe the continuum and found the feature to reside at $30.5\pm 0.4\,\mathrm{keV}$. Using a model consisting of negative and positive power-laws with a common exponential cutoff (Mihara, 1995, NPEX), Makishima et al. (1999) detected the cyclotron line at $28.8\pm 0.4\,\mathrm{keV}$. Jaisawal & Naik (2015) included a fundamental cyclotron line at $27.5\pm 0.4\,\mathrm{keV}$, $27.7\pm 0.4\,\mathrm{keV}$, or $29.6\pm 0.5\,\mathrm{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. McBride et al. (2007) and Fürst et al. (2015), both using a Fermi-Dirac continuum model, found cyclotron lines at $30.7_{-1.9}^{+1.8}\,\mathrm{keV}$ and $30.39_{-0.14}^{+0.17}\,\mathrm{keV}$ (and $29.42_{-0.24}^{+0.27}\,\mathrm{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 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 pseudo-Lorentzian 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 line forming 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 $\sim 28\,\mathrm{keV}$ to $\sim 31\,\mathrm{keV}$, and differences of $\sim 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 best-fit magnetic field strength, that is, $E_{\mathrm{gabs}}=B_{\mathrm{\texttt{cyclofs}}}\times 11.57\,\mathrm{keV}\approx 41.96^{+0.35}_{-0.36}\,\mathrm{keV}$ in the frame of the neutron star (i.e., before redshifting) 444Note 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 $\sim 32.3\,\mathrm{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 $\chi^{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 $40.01_{-0.14}^{+0.15}\,\mathrm{keV}$ before redshift (red. $\chi^{2}$ of 1.17 for 875 d.o.f.). The value of $E=40.01\,\mathrm{keV}/(1+z)=30.78_{-0.11}^{+0.12}\,\mathrm{keV}$ 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 best-fit 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 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 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 these parameters might 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 an overall good agreement for both, the slab 1-0 and the slab 1-1 geometry.

In order to show the influence of the continuum, the Her X-1 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 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\times 10^{12}\,\mathrm{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 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 (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 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.

Acknowledgements.

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 slxfig package by Davis (2014).

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

• •

setting the environment variable CYCLOFS_TABLE to the table’s location

• •

defining the table location as initialization string in the model.dat file

Table 4 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.

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 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 (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.

A.2 Angular averaging

The model is designed to calculate spectra for exactly one viewing angle, defined by its cosine mu$=\cos\vartheta$, 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 $\mu=0.5(\mu_{\mathrm{min}}+\mu_{\mathrm{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 $\mathrm{\textit{mu\_N}}=0$.