and, at the same time, be consistent with estimates for the system being
studied, which can be obtained from the emission lines. In addition, a more
complex upgrade would correspond to estimating the ionization fraction in the
accretion structure.

Another ingredient not included in the current version of the cyclops code is
the reflection from the WD surface. This process is responsible for the iron
K$\alpha$ line at ${\approx 6.4}$ keV and the Compton hump between ${\sim 10}$
and ${\sim 50}$ keV (George & Fabian, 1991). In addition, based on _NuSTAR_
observations, Mukai et al. (2015) and Luna et al. (2018) showed that
reflection is an important ingredient to be included in magnetic CV parameter
estimation, especially for IPs, which can potentially provide further
constraints. A reasonable model for the reflection was developed by Hayashi et
al. (2018), who also investigated the influence of the parameters on the iron
complex and Compton hump. A model like that could be incorporated in future
versions of the cyclops code and complementary analysis be performed (e.g.
inferring the influence of the magnetic field vector).

Finally, the thermo-hydrodynamic differential equations for post-shock flows
in the cyclops code are not solved simultaneously with the radiative transfer
for bremsstrahlung and cyclotron radiations, which is a more realistic
approach (e.g., Fischer & Beuermann, 2001). Indeed, we first model the PSR
based on the 1D approach described in Section A. Then, we perform the 3D
phase-dependent and line-of-sight-dependent radiative transfer throughout the
PSR as described in Costa & Rodrigues (2009) and Silva et al. (2013).

##  8 Summary and Conclusions

We have used an updated version of cyclops code to model X-ray spectra and
multiwavelength light curves of magnetic CVs, which includes accurate PSR
modeling, by solving the stationary one-dimensional hydro-thermodynamic
equations describing the accreting plasma, where bremsstrahlung and cyclotron
radiative processes play a role in cooling the gas from the shock until the WD
surface. We have adopted a quasi-dipolar geometry, have taken into account the
WD gravitational potential and a finite size of magnetosphere radius, and have
allowed the magnetic field to decay as the distance from the WD surface
increases. While comparing our approach with previous works, we showed that
our model agrees reasonably well with previous investigations and showed how
the temperature profiles are affected by geometrical and physical parameters.
In particular, we discussed a four-dimensional fitting scheme in modeling
magnetic CVs, whose parameters are the WD mass, the WD magnetic field, the
specific accretion rate and the threading region radius.

Our main results can be summarized as follows.

i) We showed that a balance between cyclotron and bremsstrahlung radiative
processes always exists in the presence of non-negligible cyclotron cooling,
as expected. In other words, there is always a combination of parameters such
that the importance of these emission processes in shaping the PSR profiles
are comparable.

ii) Due to the balance between bremsstrahlung and cyclotron processes, there
is always a critical specific accretion rate above which bremsstrahlung
radiation governs the cooling, and below which cyclotron radiation is the main
mechanism driving the cooling. In particular, the greater the WD magnetic
field at the PSR bottom, the greater such a critical specific accretion rate,
i.e. influence of cyclotron radiation takes place for a larger range of
specific accretion rates.

iii) Due to the balance between bremsstrahlung and cyclotron processes, the
previously found anti-correlation between the shock height and the specific
accretion rate is broken in the presence of non-negligible cyclotron cooling.
Additionally, the above-mentioned critical specific accretion rate sets the
region in the parameter space where such an anti-correlation is broken.

iv) We found that, for intermediate polars with extremely low accretion rates
(e.g. due to instability in the accretion disk or low state), the cyclotron
cooling cannot usually be neglected. This is due to the above-mentioned
balance between bremsstrahlung and cyclotron processes, which occurs at low
specific accretion rates in systems with weak magnetic fields.

v) With respect to X-ray spectra, we found that the greater the WD mass and/or
the smaller the WD magnetic field and/or the greater the specific accretion
rate and/or the greater the threading region radius, the harder the spectrum,
as expected. This is intrinsically connected with the role of such parameters
in shaping the PSR temperature and density profiles.

vi) While adopting PSRs with uniform temperature and density distributions, we
showed that the cyclops X-ray spectra are in very good agreement with those
generated by the xspec code, which is a widely used X-ray package to fit X-ray
spectra, not only with respect to distance-dependent bremsstrahlung emission
but also regarding the impact of the interstellar extinction.

vii) We showed that there exist thresholds in the parameter space below/above
which X-ray spectra become indistinguishable, resulting in a degeneracy among
the parameters. In addition, such thresholds strongly depend on the
combination of parameters. In other words, due to the balance between
bremsstrahlung and cyclotron radiative processes, the interplay between the WD
mass, the WD magnetic field and the specific accretion rate strongly affects
the thresholds mentioned above.

viii) We argued that in our four-dimensional modeling scheme, there is always
a degeneracy, regardless of the region of the parameter space one explores,
which needs to be addressed in any meaningful fitting scheme.

ix) We further argued that any attempt of solving the degeneracy problem
requires at least more constraints, i.e. having only the X-ray continuum
spectrum for a given source is not enough to look for a global minimum used to
uniquely estimate the main magnetic CV parameters.

x) We showed that even though several methods seem to some extent useful to
break the degeneracy in the parameter space, some of them alone are most
likely not enough to unambiguously solve the degeneracy problem. In
particular, despite the fact that incorporating either consistent pre-shock
region modeling or accurate distance estimates may help, most likely in a
fitting scheme one would still end up with many different models having very
similar X-ray spectra. Moreover, some very specific constraints are not always
available, such as the PSR height or the break frequency. However,
incorporating light curves to the fitting scheme looks like a promising way
out of the degeneracy problem and, based on that, we proposed a method that
seems able to solve this problem.

xi) The method we proposed corresponds to incorporating X-ray light curves in
different energy ranges to the fitting scheme, in case there is enough signal-
to-noise ratio for that, potentially together with accurate _Gaia_ distance
estimates. This is an appealing method since X-ray data already comprise
spectra and light curves, and no extra observational effort is required.
However, this method is expected to work only if the pre-shock region is
properly modeled, since this region plays a key role in shaping X-ray light
curves.

xii) While investigating the energy-dependent pulsed fraction in the X-ray
light curves modeling, we have shown that the partial self-eclipse of the PSR
should not be neglected, since it can produce pulsed fractions that smoothly
decrease with energy. This implies that the observed decrease of the pulsed
fraction with energy in some magnetic CVs is not necessarily due to
photoeletric absorption of the pre-shock region alone, and could perhaps be at
least partially explained by self-eclipse, depending on the geometrical
properties of system. In addition, this combination of self-eclipse and pre-
shock region absorption could lead to energy-dependent light curves, in the
sense that the location or existence of the minima present in the X-ray light
curves could depend on the energy.

xiii) Finally, we showed that the cyclops code, which is a 3D emission
modeling code, is a sophisticated tool that can and should be used by the
community in future magnetic CV emission modeling, not only with respect to
X-ray data, but also regarding data in virtually any frequency range,
including polarized emission. In addition, to the best of our knowledge,
cyclops is the only code available to the community able to consistently model
the 3D pre-shock region, which is an essential ingredient to properly deal
with X-ray emission from magnetic CVs.

## Availability of the cyclops code

The cyclops code is written in idl and is freely available under request. In
order to obtain a copy of the cyclops code described in this paper, please
write to<EMAIL_ADDRESS>
## Acknowledgements

We would like to thank an anonymous referee for the comments and suggestions
that helped to improve this paper. We thank Alfredo N. Iusem for pleasantly
explaining many issues associated with the Newton-Raphson root finding method
to solve boundary value problems. We also thank the MCTIC/FINEP (CT-INFRA
grant 0112052700) and the Embrace Space Weather Program for the computing
facilities at the National Institute for Space Research, Brazil. D.B. was
supported by São Paulo Research Foundation (FAPESP, grant #2017/14289-3) and
partially supported by ESO/Gobierno de Chile. C.V.R. thanks the grant
303444/2018-5 from CNPq and the grant #2017/14289-3 from São Paulo Research
Foundation (FAPESP). M.R.S. acknowledges support from the Millennium Nucleus
for Planet Formation (NPF) and Fondecyt (grant 1181404). M.C. thanks the grant
#2015/25972-0 by FAPESP. G.J.M.L. is a member of the CIC-CONICET (Argentina)
and acknowledge support from grant ANPCYT-PICT 0901/2017. A.S.O acknowledges
São Paulo Research Foundation (FAPESP) for financial support under grant
#2017/20309-7. S.G.P. acknowledges the support of a Science and Technology
Facilities Council (STFC) Ernest Rutherford Fellowship. P.E.S. acknowledges
FAPESP for financial support under grant #2017/13551-6. M.Z. acknowledges
support from CONICYT PAI (Concurso Nacional de Inserción en la Academia 2017,
Folio 79170121) and CONICYT/FONDECYT (Programa de Iniciación, Folio 11170559).

## References

  * Aizu (1973) Aizu, K. 1973, Progress of Theoretical Physics, 49, 1184
  * Arnaud (1996) Arnaud, K. A. 1996, in ASP Conf. Ser., Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Balucinska-Church & McCammon (1992) Balucinska-Church, M., & McCammon, D. 1992, ApJ, 400, 699
  * Belloni et al. (2021) Belloni, D., Schreiber, M. R., Salaris, M., Maccarone, T. J., & Zorotovic, M. 2021, MNRAS, 505, L74
  * Bernardini et al. (2017) Bernardini, F., de Martino, D., Mukai, K., et al. 2017, MNRAS, 470, 4815
  * Canalle et al. (2005) Canalle, J. B. G., Saxton, C. J., Wu, K., Cropper, M., & Ramsay, G. 2005, A&A, 440, 185
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Costa & Rodrigues (2009) Costa, J. E. R., & Rodrigues, C. V. 2009, MNRAS, 398, 240
  * Coti Zelati et al. (2016) Coti Zelati, F., Rea, N., Campana, S., et al. 2016, MNRAS, 456, 1913
  * Cropper (1990) Cropper, M. 1990, Space Sci. Rev., 54, 195
  * Cropper et al. (1998) Cropper, M., Ramsay, G., & Wu, K. 1998, MNRAS, 293, 222
  * Cropper et al. (1999) Cropper, M., Wu, K., Ramsay, G., & Kocabiyik, A. 1999, MNRAS, 306, 684
  * Dorman & Arnaud (2001) Dorman, B., & Arnaud, K. A. 2001, in ASP Conf. Ser., Vol. 238, Astronomical Data Analysis Software and Systems X, ed. J. Harnden, F. R., F. A. Primini, & H. E. Payne, 415
  * Ezuka & Ishida (1999) Ezuka, H., & Ishida, M. 1999, ApJS, 120, 277
  * Ferrario et al. (2015) Ferrario, L., de Martino, D., & Gänsicke, B. T. 2015, Space Sci. Rev., 191, 111
  * Ferrario et al. (2020) Ferrario, L., Wickramasinghe, D., & Kawka, A. 2020, Advances in Space Research, 66, 1025
  * Ferrario et al. (1989) Ferrario, L., Wickramasinghe, D. T., & Tuohy, I. R. 1989, ApJ, 341, 327
  * Fischer & Beuermann (2001) Fischer, A., & Beuermann, K. 2001, A&A, 373, 211
  * Fujimoto & Ishida (1997) Fujimoto, R., & Ishida, M. 1997, ApJ, 474, 774
  * George & Fabian (1991) George, I. M., & Fabian, A. C. 1991, MNRAS, 249, 352
  * Gronenschild & Mewe (1978) Gronenschild, E. H. B. M., & Mewe, R. 1978, A&AS, 32, 283
  * Hōshi (1973) Hōshi, R. 1973, Progress of Theoretical Physics, 49, 776
  * Hameury et al. (1987) Hameury, J. M., King, A. R., Lasota, J. P., & Ritter, H. 1987, ApJ, 316, 275
  * Hayashi & Ishida (2014a) Hayashi, T., & Ishida, M. 2014a, MNRAS, 438, 2267
  * Hayashi & Ishida (2014b) —. 2014b, MNRAS, 441, 3718
  * Hayashi et al. (2011) Hayashi, T., Ishida, M., Terada, Y., Bamba, A., & Shionome, T. 2011, PASJ, 63, S739
  * Hayashi et al. (2018) Hayashi, T., Kitaguchi, T., & Ishida, M. 2018, MNRAS, 474, 1810
  * Hellier (2001) Hellier, C. 2001, Cataclysmic Variable Stars (Springer)
  * Hoogerwerf et al. (2006) Hoogerwerf, R., Brickhouse, N. S., & Mauche, C. W. 2006, ApJ, 643, L45
  * Imamura et al. (1987) Imamura, J. N., Durisen, R. H., Lamb, D. Q., & Weast, G. J. 1987, ApJ, 313, 298
  * Ishida & Ezuk (1999) Ishida, M., & Ezuk, H. 1999, in ASP Conf. Ser., Vol. 157, Annapolis Workshop on Magnetic Cataclysmic Variables, ed. C. Hellier & K. Mukai, 333
  * Kaastra (1992) Kaastra, J. 1992, Internal SRON-Leiden Report, updated version 2.0
  * Karzas & Latter (1961) Karzas, W. J., & Latter, R. 1961, ApJS, 6, 167
  * Kashyap & Drake (2000) Kashyap, V., & Drake, J. J. 2000, Bulletin of the Astronomical Society of India, 28, 475
  * Liedahl et al. (1995) Liedahl, D. A., Osterheld, A. L., & Goldstein, W. H. 1995, ApJ, 438, L115
  * Luna et al. (2018) Luna, G. J. M., Mukai, K., Orio, M., & Zemko, P. 2018, ApJ, 852, L8
  * Luna et al. (2010) Luna, G. J. M., Raymond, J. C., Brickhouse, N. S., et al. 2010, ApJ, 711, 1333
  * Luna et al. (2015) Luna, G. J. M., Raymond, J. C., Brickhouse, N. S., Mauche, C. W., & Suleimanov, V. 2015, A&A, 578, A15
  * Meggitt & Wickramasinghe (1982) Meggitt, S. M. A., & Wickramasinghe, D. T. 1982, MNRAS, 198, 71
  * Mewe et al. (1985) Mewe, R., Gronenschild, E. H. B. M., & van den Oord, G. H. J. 1985, A&AS, 62, 197
  * Mewe et al. (1986) Mewe, R., Lemen, J. R., & van den Oord, G. H. J. 1986, A&AS, 65, 511
  * Mhlahlo et al. (2007) Mhlahlo, N., Buckley, D. A. H., Dhillon, V. S., et al. 2007, MNRAS, 378, 211
  * Mukai (2017) Mukai, K. 2017, PASP, 129, 062001
  * Mukai et al. (2015) Mukai, K., Rana, V., Bernardini, F., & de Martino, D. 2015, ApJ, 807, L30
  * Nauenberg (1972) Nauenberg, M. 1972, ApJ, 175, 417
  * Norton et al. (2008) Norton, A. J., Butters, O. W., Parker, T. L., & Wynn, G. A. 2008, ApJ, 672, 524
  * Oliveira et al. (2019) Oliveira, A. S., Rodrigues, C. V., Palhares, M. S., et al. 2019, MNRAS, 489, 4032
  * Pacholczyk & Swihart (1975) Pacholczyk, A. G., & Swihart, T. L. 1975, ApJ, 196, 125
  * Patterson (1994) Patterson, J. 1994, PASP, 106, 209
  * Ramsay (2000) Ramsay, G. 2000, MNRAS, 314, 403
  * Ramsay & Cropper (2004) Ramsay, G., & Cropper, M. 2004, MNRAS, 347, 497
  * Revnivtsev et al. (2009) Revnivtsev, M., Churazov, E., Postnov, K., & Tsygankov, S. 2009, A&A, 507, 1211
  * Revnivtsev et al. (2011) Revnivtsev, M., Potter, S., Kniazev, A., et al. 2011, MNRAS, 411, 1317
  * Revnivtsev et al. (2010) Revnivtsev, M., Burenin, R., Bikmaev, I., et al. 2010, A&A, 513, A63
  * Rybicki & Lightman (1986) Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics (WILEY-VCH)
  * Sarty et al. (2008) Sarty, G. E., Saxton, C. J., & Wu, K. 2008, Ap&SS, 317, 239
  * Saxton et al. (2007) Saxton, C. J., Wu, K., Canalle, J. B. G., Cropper, M., & Ramsay, G. 2007, MNRAS, 379, 779
  * Schreiber et al. (2021) Schreiber, M. R., Belloni, D., Gänsicke, B. T., Parsons, S. G., & Zorotovic, M. 2021, Nature Astronomy, arXiv:2104.14607
  * Schure et al. (2009) Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751
  * Semena et al. (2014) Semena, A. N., Revnivtsev, M. G., Buckley, D. A. H., et al. 2014, MNRAS, 442, 1123
  * Shaw et al. (2020) Shaw, A. W., Heinke, C. O., Mukai, K., et al. 2020, MNRAS, 498, 3457
  * Silva et al. (2013) Silva, K. M. G., Rodrigues, C. V., Costa, J. E. R., et al. 2013, MNRAS, 432, 1587
  * Suleimanov et al. (2016) Suleimanov, V., Doroshenko, V., Ducci, L., Zhukov, G. V., & Werner, K. 2016, A&A, 591, A35
  * Suleimanov et al. (2008) Suleimanov, V., Poutanen, J., Falanga, M., & Werner, K. 2008, A&A, 491, 525
  * Suleimanov et al. (2005) Suleimanov, V., Revnivtsev, M., & Ritter, H. 2005, A&A, 435, 191
  * Suleimanov et al. (2019) Suleimanov, V. F., Doroshenko, V., & Werner, K. 2019, MNRAS, 482, 3622
  * Tomsick et al. (2016) Tomsick, J. A., Rahoui, F., Krivonos, R., et al. 2016, MNRAS, 460, 513
  * Van Box Som et al. (2018) Van Box Som, L., Falize, É., Bonnet-Bidaud, J.-M., et al. 2018, MNRAS, 473, 3158
  * Warner (1995) Warner, B. 1995, Cambridge Astrophysics Series, 28
  * Wickramasinghe & Ferrario (2000) Wickramasinghe, D. T., & Ferrario, L. 2000, PASP, 112, 873
  * Wu (2000) Wu, K. 2000, Space Sci. Rev., 93, 611
  * Wu et al. (1994) Wu, K., Chanmugam, G., & Shaviv, G. 1994, ApJ, 426, 664
  * Yuasa et al. (2010) Yuasa, T., Nakazawa, K., Makishima, K., et al. 2010, A&A, 520, A25
  * Zhu et al. (2017) Zhu, H., Tian, W., Li, A., & Zhang, M. 2017, MNRAS, 471, 3494

##  Appendix A Post-shock region (PSR)

###  A.1 Basic Assumptions and Geometry

In our modeling, the PSR properties are obtained from the solution of the
stationary one-dimensional hydro-thermodynamic differential equations
describing the accreting plasma. We further assume a dipole-like geometry
(i.e. cubic cross-section variation, e.g., Hayashi & Ishida, 2014a; Suleimanov
et al., 2016) and equipartition between ions and electrons (e.g., Wu et al.,
1994; Cropper et al., 1999; Van Box Som et al., 2018). We would like to draw
the readers attention to the fact that, unlike the PSR modeling, a three-
dimensional geometry is adopted in the remaining parts of the cyclops code,
e.g., the radiative transport and pre-shock region. Hereafter, we will only
discuss the PSR modeling, although in Section 2 we provide more details about
the other aspects of the cyclops code.

A stationary model with equipartition is a reasonable assumption provided that
the characteristic time-scale of energy exchanged between electrons and ions
(i.e. electron-ion equipartition time-scale) is much shorter than the cooling
time-scale. There are two particular situations in which such an assumption
does not hold. First, for very low specific accretion rates and massive WDs,
the non-equipartition area in the PSR might be well extended (${\gtrsim 80}$
%), which means that the equipartition approach is no longer valid (e.g.,
Hayashi & Ishida, 2014a). This is because of the low density within the PSR,
which makes energy exchange between electrons and ions inefficient. However,
Hayashi & Ishida (2014b) applied their model considering non-equipartition
between ions and electron to the IPs EX Hya and V1223 Sgr and found that even
for a relatively small specific accretion rate ($\sim 0.05$ g s-1 cm-2), the
equipartition between ions and electrons is quickly achieved close to the
shock. Second, for very strong magnetic fields and massive WDs, models
assuming non-equipartition predict considerably different shock heights and
qualitatively different density and temperature distributions in regions far
from the WD surface than models assuming equipartition (e.g., Wu, 2000; Saxton
et al., 2007; Sarty et al., 2008).

With respect to the cross-section variation, a dipolar geometry is a more
realistic approach than standard cylindrical models in which the cross-section
is assumed to be constant. Indeed, Saxton et al. (2007) showed that the
inclusion of a dipolar accretion geometry provides a harder continuum and in
turn lowers estimates for WD masses when compared with masses inferred from
models excluding this effect. In addition, Hayashi & Ishida (2014a) showed
that differences between both approaches become important when the shock
height is ${\gtrsim 20}$ % the WD radius (see also Canalle et al., 2005). This
is particularly true for accretion columns in IPs, in which the magnetic field
is only moderately strong and the shock position might be relatively distant
from the WD surface so that changes in the cross-section might play a
significant role.

In Fig. 2 we illustrate the geometrical and physical aspects of our PSR model.
In what follows, we describe the differential equations, boundary conditions
and the numerical scheme to solve these equations.

###  A.2 Boundary Value Problem

The steady state equations presented below are for a plasma that is restricted
to flow along the WD magnetic field lines. We assumed here that the plasma is
a fully ionized gas respecting the ideal gas law

$P\ =\ n\,k_{B}\,T\,,$ (A1)

where $P$ is the average pressure, $T$ is the average temperature,
$k_{B}=1.380658\times 10^{-16}$ erg K-1 is the Boltzmann constant, and $n$ in
the average number density given by

$n\ =\ \frac{\rho}{\mu\,m_{H}}\,,$ (A2)

where $\rho$ is the mass density, ${\mu=0.615}$ is the mean molecular mass of
the fully ionized gas with solar metallicity, and ${m_{H}=1.673525\times
10^{-24}}$ g is the atomic hydrogen mass.

We further assume that the flow is symmetric with respect to the magnetic line
passing through the PSR center, variations in the cross-section mimic a
dipolar geometry (as suggested by Hayashi & Ishida, 2014a) and that the hydro-
thermodynamic quantities are functions of only $z$ (see Fig. 2) within the
PSR. Under such assumptions, the differential equations related to the mass
continuity, the momentum conservation and the energy conservation are (e.g.,
Hayashi & Ishida, 2014a):

$\frac{{\rm d}}{{\rm d}z}\left(S\rho v\right)\ =\ 0,$ (A3)

$\frac{{\rm d}}{{\rm d}z}\left(\rho v^{2}+P\right)+\frac{\rho
v^{2}}{S}\frac{{\rm d}S}{{\rm d}z}+\rho\,g_{\rm WD}=\ 0,$ (A4)

$v\frac{{\rm d}P}{{\rm d}z}+\gamma P\frac{{\rm d}v}{{\rm
d}z}+(\gamma-1)\left(\Lambda-\frac{\rho v^{3}}{2S}\frac{{\rm d}S}{{\rm
d}z}\right)=\ 0,$ (A5)

where $\gamma=5/3$ is the adiabatic index, $v(z)$ is the flow velocity, $S(z)$
is the cross-section, $\Lambda(z)$ is the cooling function, which is described
in Section A.4, and $g_{\rm WD}(z)$ is the WD gravitational field given by

$g_{\rm WD}(z)\ =\ \frac{GM_{\rm WD}}{z^{2}},$ (A6)

where $G=6.67259\times 10^{-8}$ cm3 g-1 s-2 is the gravitational constant and
$M_{\rm WD}$ is the WD mass.

Note that Equation A3 contains an integral of the form

$S(z)\,\rho(z)\,v(z)\ =\ S(z)\,\dot{m}(z)\ =\ \dot{M}_{\rm WD}$ (A7)

where $\dot{M}_{\rm WD}$ is the accretion rate and $\dot{m}(z)$ is the
specific accretion rate, i.e. the accretion rate per unit area. Notice that
$\dot{M}_{\rm WD}$ is constant but both $S$ and $\dot{m}$ vary throughout the
PSR such that their product is constant and equals to $\dot{M}_{\rm WD}$. We
follow here the formulation by Hayashi & Ishida (2014a) and define the cross-
section as

$S(z)\ =\ S_{\rm b}\,\left(\frac{z}{R_{\rm WD}}\right)^{n},$ (A8)

where $S_{\rm b}$ is the accretion area at the bottom of the PSR (see Fig. 2),
$n=0$ corresponds to the cylindrical geometry and $n=3$ to the dipolar
geometry. The bottom of the PSR corresponds to ${z=R_{\rm WD}}$ and ${S=S_{\rm
b}}$, the latter coming directly from the adopted magnetic accretion geometry.

The cross-section derivative through the PSR is then given by

$\frac{{\rm d}S}{{\rm d}z}\ =\ n\,\left(\frac{S_{\rm b}}{R_{\rm
WD}}\right)\,\left(\frac{z}{R_{\rm WD}}\right)^{n-1}.$ (A9)

In order to facilitate the solution of the PSR structure, instead of solving
Equations A3, A4 and A5 in the variable $z$ we utilize a different independent
variable. The variable $z$ varies from $z=z_{\rm sh}=R_{\rm WD}+H_{\rm sh}$
until $z=R_{\rm WD}$, where $z_{\rm sh}$ is the shock position (i.e. $z_{\rm
sh}=R_{\rm WD}+H_{\rm sh}$), $H_{\rm sh}$ is the shock height with respect to
the WD surface (see Fig. 2), and $R_{\rm WD}$ is the WD radius given by
(Nauenberg, 1972)

$\frac{R_{\rm WD}}{{\rm cm}}\ =\ 7.8\times 10^{8}\,\sqrt{\left(\frac{M_{\rm
WD}}{1.44\,{\rm M}_{\odot}}\right)^{-2/3}\ -\ \left(\frac{M_{\rm
WD}}{1.44\,{\rm M}_{\odot}}\right)^{2/3}}.$ (A10)

We define a variable $z^{\prime}$ such that it goes from $0$ to $H_{s}$ (e.g.,
Suleimanov et al., 2005), i.e.

$z^{\prime}\ :=\ z_{\rm sh}-z.$

Equations A3, A4 and A5 can be translated into the following two-dimensional
system of ordinary differential equations (Hayashi & Ishida, 2014a), which has
$z^{\prime}$ as the independent variable:

$\frac{{\rm d}P}{{\rm
d}z^{\prime}}=\frac{\left(\gamma-1\right)\dot{m}\left(\Lambda-\frac{\rho
v^{3}}{2S}\frac{{\rm d}S}{{\rm d}z^{\prime}}\right)\ +\
\left(\dot{m}/v\right)\gamma Pg_{\rm WD}}{\gamma P\ -\ \dot{m}v},$ (A11)

$\frac{{\rm d}v}{{\rm d}z^{\prime}}=\frac{g_{\rm
WD}}{v}-\frac{1}{\dot{m}}\,\frac{{\rm d}P}{{\rm d}z^{\prime}}.$ (A12)

The system variables are the pressure $(P)$ and the velocity $(v)$ and they
depend on $z^{\prime}$, which is the independent variable. The system of
ordinary differential equations expressed by Equations A11 and A12 is
associated with boundary conditions (Section A.3), which all together
correspond to the boundary value problem defining the PSR.

###  A.3 Boundary Conditions

The system described by Equations A11 and A12 is associated with a boundary
value problem, in which the boundary conditions are given at the shock
position (i.e. $z^{\prime}=z_{\rm sh}-z=0$) and at the WD surface (i.e.
$z^{\prime}=z_{\rm sh}-z=H_{\rm sh}$). At the shock position we assume the
Rankine-Hugoniot condition and at the WD surface we impose that the velocity
is null.

The pre- and the post-shock values of the hydrodynamic variables are related
by the Rankine-Hugoniot strong-shock jump conditions. At the shock (just above
it), the plasma falls with its free-fall velocity $(v_{\rm ff})$. The
conditions for the jump are such that, at the shock, the density increases by
a factor of four and the velocity decreases by a factor of four, i.e.

$v_{\rm post}\ =\ \frac{v_{\rm pre}}{4}\ =\ \frac{v_{\rm
ff}}{4}=\frac{1}{4}\,\sqrt{2GM_{\rm WD}\,\left(\frac{1}{z_{\rm
sh}}-\frac{1}{R_{\rm th}}\right)},$ (A13)

$\rho_{\rm post}\ =\ 4\rho_{\rm pre}\ =\ 4\frac{\dot{m}}{v_{\rm pre}}\ =\
4\frac{\dot{m}}{v_{\rm ff}}\ =\ \frac{\dot{m}}{v_{\rm post}},$ (A14)

where $R_{\rm th}$ is the threading region position relative to the WD center
(see Fig. 1) and corresponds to the point where the magnetic pressure first
dominates the gas ram pressure. This position can also be understood as the
magnetosphere limit.

As already pointed out by Suleimanov et al. (2016), the assumption of a finite
magnetosphere is more accurate than assuming a free-fall velocity from the
infinity. This is rather important for systems with small magnetospheres (a
few $R_{\rm WD}$), where matter reaches considerably smaller velocities
compared to the case of accretion from infinity (see also Luna et al., 2018).

The boundary conditions at the shock (i.e. ${z^{\prime}=0}$), which are
indicated with sub-index ‘${\rm sh}$’, are then given by (e.g., Suleimanov et
al., 2005)

$\left\\{\begin{array}[]{l}v_{\rm sh}\ =\ 0.25\,\sqrt{2GM_{\rm
WD}\,\left(1/z_{\rm sh}-1/R_{\rm th}\right)},\\\ \\\ \rho_{\rm sh}\ =\
\dot{m}_{\rm sh}\,/\,v_{\rm sh},\\\ \\\ P_{\rm sh}\ =\ 3\,\dot{m}_{\rm
sh}\,v_{\rm sh},\\\ \\\ T_{\rm sh}\ =\ 3\,\mu\,m_{H}\,v^{2}_{\rm
sh}\,/\,k_{B},\\\ \end{array}\right.$ (A15)

and the condition at the bottom, i.e., at the WD surface is

$v_{\rm b}\ =\ v(z^{\prime}=H_{\rm sh})\ =\ 0.$ (A16)

Before describing the numerical scheme to solve the boundary value problem
associated with the PSR (Section A.5), we discuss the cooling function adopted
here.

###  A.4 Cooling Function

The cooling function $(\Lambda)$ adopted here includes effects of both
bremsstrahlung and the cyclotron processes, i.e.

$\Lambda\ =\ \Lambda_{\rm brem}+\Lambda_{\rm cyc}.$ (A17)

The optically thin bremsstrahlung cooling is given by (Van Box Som et al.,
2018)

$\frac{\Lambda_{\rm brem}}{{\rm erg\,cm^{-3}\,s^{-1}}}\ =\ 5.01\times
10^{6}\,\left(\frac{\rho}{10^{-9}\,{\rm
g\,cm^{-3}}}\right)^{2}\,\left(\frac{T}{10^{8}\,{\rm K}}\right)^{1/2}.$ (A18)

To evaluate which one of these two radiative processes dominates the cooling,
one usually utilizes the cooling ratio ($\epsilon_{\rm sh}$), which is defined
as the ratio between the bremsstrahlung and the cyclotron cooling time-scales
at the shock position, given by (Van Box Som et al., 2018)

$\begin{multlined}\epsilon_{\rm sh}\ =\ 4.32\,\left(\frac{B_{\rm sh}}{10\,{\rm
MG}}\right)^{2.85}\left(\frac{\dot{m}_{\rm sh}}{{\rm
g\,s^{-1}\,cm^{-2}}}\right)^{-1.85}\left(\frac{M_{\rm WD}}{0.8\,{\rm
M_{\odot}}}\right)^{2.925}\left(\frac{R_{\rm WD}}{7\times 10^{8}\,{\rm
cm}}\right)^{-2.925}\left(\frac{S_{\rm sh}}{10^{15}\,{\rm
cm^{2}}}\right)^{-0.425},\end{multlined}\epsilon_{\rm sh}\ =\
4.32\,\left(\frac{B_{\rm sh}}{10\,{\rm
MG}}\right)^{2.85}\left(\frac{\dot{m}_{\rm sh}}{{\rm
g\,s^{-1}\,cm^{-2}}}\right)^{-1.85}\left(\frac{M_{\rm WD}}{0.8\,{\rm
M_{\odot}}}\right)^{2.925}\left(\frac{R_{\rm WD}}{7\times 10^{8}\,{\rm
cm}}\right)^{-2.925}\left(\frac{S_{\rm sh}}{10^{15}\,{\rm
cm^{2}}}\right)^{-0.425},$ (A19)

where values with sub-index ‘${\rm sh}$’ are quantities at the shock position.
Bremsstrahlung dominates in cases where ${\epsilon_{\rm sh}<1}$ and cyclotron
governs otherwise (see Section 3 for more details).

The optically thick cyclotron cooling is given by (e.g., Canalle et al., 2005;
Saxton et al., 2007)

$\begin{multlined}\Lambda_{\rm cyc}\ =\ \Lambda_{\rm brem}\,\epsilon_{\rm
sh}\left(\frac{P}{P_{\rm sh}}\right)^{2}\left(\frac{\rho}{\rho_{\rm
sh}}\right)^{-3.85}\left(\frac{S}{S_{\rm
sh}}\right)^{-0.425}\left(\frac{B}{B_{\rm
sh}}\right)^{2.85},\end{multlined}\Lambda_{\rm cyc}\ =\ \Lambda_{\rm
brem}\,\epsilon_{\rm sh}\left(\frac{P}{P_{\rm
sh}}\right)^{2}\left(\frac{\rho}{\rho_{\rm
sh}}\right)^{-3.85}\left(\frac{S}{S_{\rm
sh}}\right)^{-0.425}\left(\frac{B}{B_{\rm sh}}\right)^{2.85},$ (A20)

where the quantities $P$, $\rho$, $S$ and $B$ are not fixed and vary through
the PSR and again quantities with sub-index ‘${\rm sh}$’ are given at the
shock position. In particular, the cross-section variation is given by Eq. A8
and the WD magnetic field assumed here is a dipole, which implies that its
variation with $z$ is given by

$B(z)\ =\ B_{\rm b}\,\left(\frac{z}{R_{\rm WD}}\right)^{-3},$ (A21)

where $B_{\rm b}$ is the magnetic field at the PSR bottom given by

$B_{\rm b}\ =\ \frac{B_{\rm p}}{2}\,\sqrt{4-3\cos^{2}\varrho},$ (A22)

where $B_{\rm p}$ is the magnetic field intensity at the pole, and $\varrho$
is a function of $\beta$, $B_{\rm lat}$ and $B_{\rm long}$. We would like to
stress that in the cyclops code, the parameter is $B_{\rm p}$ (i.e. the
magnetic field intensity at the pole), while $B_{\rm b}$ (i.e. the magnetic
field intensity at the PSR bottom) is computed accordingly provided the
geometry of the model (i.e. the PSR center colatitude $\beta$ and the magnetic
field axis latitude $B_{\rm lat}$ and longitude $B_{\rm long}$). See Fig. 1
for more details on the geometrical aspects of the cyclops code. Finally, for
$B_{\rm lat}=90^{o}$ and $B_{\rm long}=0^{o}$, i.e. in the case the magnetic
axis coincides with the WD rotation axis and the PSR and magnetic axis are in
the same plane (like in Fig. 1), Eq. A22 reduces to an expression consistent
with eq. 45 in Canalle et al. (2005).

###  A.5 Numerical Scheme

The boundary value problem is given by the system of ordinary differential
equations defined by the Eqs. A11 and A12, coupled with the conditions A15 and
A16, where the quantity to be determined is the shock height $H_{\rm sh}$. To
solve such a boundary value problem, we iteratively solve the ordinary
differential equations. To do so, we guess shock heights and assume the
conditions A15 as the initial conditions, until we satisfy the condition A16,
at the WD surface, i.e. null velocity at the PSR bottom.

Since the derivatives of the system tend to be huge as the velocity tends to
zero (due to the asymptotic behavior of the solutions while approaching the
PSR bottom), we solve the initial value problem with a 4–5th-order Runge-
Kutta-Fehlberg integrator with adaptive step-size control, in order to reach
very small velocities (i.e. with a precision of ${\sim 10^{-8}}$). Indeed,
given the non-smooth behavior of the system when ${v\rightarrow 0}$,
monitoring the local truncation error here is the best way to ensure accuracy,
by adjusting the step-size.

We solve the whole boundary value problem with a simple shooting algorithm
based on the Newton-Raphson root finding method. We start with a trial
integration that satisfies the boundary conditions at the shock position, from
${v=v_{\rm sh}}$ to ${v=10^{-8}}$ cm s-1. The final position (i.e. PSR bottom)
is then compared with the WD surface position (i.e. ${z^{\prime}=H_{\rm sh}}$)
and the discrepancy is used to adjust the guessing $H_{\rm sh}$, until the
boundary condition at the PSR bottom is eventually satisfied, with a precision
of ${\sim 10^{-8}}$.

##  Appendix B Consistency with previous works

Figure 20: Comparison of different cyclops solutions with previous works. In
the top and middle rows, we compare our temperature profiles with Cropper et
al. (1999, top left panel), Van Box Som et al. (2018, top right panel) and
Hayashi & Ishida (2014a, middle left and right panels). In the bottom left
panel, we show the temperature at the shock $(T_{\rm sh})$ as a function of
the WD mass $(M_{\rm WD})$ and compare with Yuasa et al. (2010). In the bottom
right panel, we provide the shock height in units of the WD radius $(H_{\rm
sh}/R_{\rm WD})$ as a function of the specific accretion rate at the WD
surface $(\dot{m}_{\rm b})$ and compare with Hayashi & Ishida (2014a). Our
profiles are in very good agreement with these previous works, especially the
temperature at the shock and the profile shape as a whole. Small differences
in the shock height ${(\lesssim\,10\,\%)}$ are likely due to different cooling
prescriptions as well as numerical techniques to solve the system of ordinary
differential equations defining the PSR structure.

With the aim of showing that our modeling provide consistent results, we
compare in this section our approach with several previous works with similar
methods. Starting with our PSR modeling, we show, in the top and middle panels
of Fig. 20, temperature profiles built with parameters consistent with Cropper
et al. (1999, top left panel), Van Box Som et al. (2018, top right panel) and
Hayashi & Ishida (2014a, middle panels). In the comparisons, a cylindrical
model is defined by setting $n=0$ in Eqs. A8 and A9, while a dipolar model is
defined by $n=3$.

In the comparison with Cropper et al. (1999), we assumed a cylindrical model
and set ${R_{\rm th}\gg R_{\rm WD}}$, ${B=0}$, ${\dot{m}=1\leavevmode\nobreak\
{\rm g\,s^{-1}\,cm^{-2}}}$, ${M_{\rm WD}=1.0\leavevmode\nobreak\ {\rm
M_{\odot}}}$ and ${\mu=0.615}$. In the comparison with Van Box Som et al.
(2018), we also assumed the cylindrical model and set ${R_{\rm th}\gg R_{\rm
WD}}$, ${B=10}$ MG constant throughout the PSR,
${\dot{m}=1\leavevmode\nobreak\ {\rm g\,s^{-1}\,cm^{-2}}}$, ${M_{\rm
WD}=0.8\leavevmode\nobreak\ {\rm M}_{\odot}}$ and ${\mu=0.5}$. Finally, in the
comparison with Hayashi & Ishida (2014a), we assumed both cylindrical and
dipolar models and set ${R_{\rm th}\gg R_{\rm WD}}$, ${B=0}$, ${M_{\rm
WD}=0.7\leavevmode\nobreak\ {\rm M}_{\odot}}$ and ${\mu=0.615}$. We then
considered two low values for $\dot{m}_{\rm b}$, namely ${10^{-4}\,{\rm
g\,s^{-1}\,cm^{-2}}}$ and ${10^{-2}\,{\rm g\,s^{-1}\,cm^{-2}}}$.

Notice that our profiles are in good agreement with results from the above-
mentioned works and the different shock heights $(\lesssim\,10\,\%)$ are
likely due to different cooling prescriptions and numerical methods. For
instance, Hayashi & Ishida (2014a) adopted the collisional ionization
equilibrium cooling function calculated by spex package (Schure et al., 2009),
in which cooling due to emission lines, which becomes more efficient for
temperatures below a few $10^{6}$ K, is included. On the other hand, cyclops
does not currently handle line cooling and only bremsstrahlung cooling
contributes to the cooling function in the middle panels of Fig. 20. This
explains why the shock heights in cyclops are slightly greater than those from
Hayashi & Ishida (2014a) since the inclusion of line contribution to the
cooling function slightly enhances the cooling, resulting in turn in slightly
shorter PSRs.

Moving forward, we compare our results with Yuasa et al. (2010) and Hayashi &
Ishida (2014a), for a set of solutions, in which we assumed parameters
consistent with those used by these authors, i.e. ${R_{\rm th}\gg R_{\rm
WD}}$, ${B=0}$ and ${\mu=0.615}$. In the comparison with Yuasa et al. (2010),
we adopted the cylindrical model, fixed ${\dot{m}=10\leavevmode\nobreak\ {\rm
g\,s^{-1}\,cm^{-2}}}$, and let $M_{\rm WD}$ vary from $0.2$ to $1.4$ M⊙. In
the comparison with Hayashi & Ishida (2014a), we adopted the cylindrical and
dipolar models, fixed ${M_{\rm WD}=1.2\leavevmode\nobreak\ {\rm M}_{\odot}}$
and let $\dot{m}_{\rm b}$ vary from $10^{-4}$ to ${10^{2}\leavevmode\nobreak\
{\rm g\,s^{-1}\,cm^{-2}}}$.

The comparisons are shown in the bottom row of Fig. 20, where in the left
panel we show the shock temperature $T_{\rm sh}$ against $M_{\rm WD}$ similar
to fig. 11 in Yuasa et al. (2010), and in the right panel we show the shock
height $H_{\rm sh}$ from the WD surface against $\dot{m}_{\rm b}$ similar to
fig. 5 in Hayashi & Ishida (2014a). Notice that we also find a good agreement
between our results and theirs. Again, small differences found in these plots
are connected with different bremsstrahlung cooling prescriptions and
numerical techniques to solve the boundary value problem associated with the
PSR structure. In particular, we predict slightly higher (lower) shock
temperatures for low-mass (high-mass) WDs, in comparison with Yuasa et al.
(2010). Additionally, we predict slightly larger shock height, especially for
high ${\dot{m}_{\rm b}\leavevmode\nobreak\ (\gtrsim 1\leavevmode\nobreak\ {\rm
g\,s^{-1}\,cm^{-2}})}$, in comparison with Hayashi & Ishida (2014a), because
line cooling is not included in the cyclops code, which makes the overall
cooling less efficient in our case.

The correlation and anti-correlation shown in the bottom row of Fig. 20 are
interesting and not difficult to understand. With regards to the correlation
between the shock temperature $(T_{\rm sh})$ and $M_{\rm WD}$, we notice that
$T_{\rm sh}$ strongly depends on the gas velocity at the shock position
$(v_{\rm sh})$, which in turn depends on $M_{\rm WD}$, as $v_{\rm sh}$ is a
result of the potential energy converted into kinetic energy. In this way, the
greater $M_{\rm WD}$, the greater $v_{\rm sh}$, an in turn the greater $T_{\rm
sh}$.

Now, the anti-correlation between $H_{\rm sh}$ and $\dot{m}_{\rm b}$ exists
for negligible cyclotron cooling cases and is directly related to the
bremsstrahlung cooling efficiency in the PSR. Indeed, the greater
$\dot{m}_{\rm b}$, the greater the density in the PSR. As the bremsstrahlung
cooling rate strongly depends on the density (see Eq. A18), the greater the
density, the greater the cooling rate, and in turn the smaller the PSR height.
This implies then that the smaller $\dot{m}_{\rm b}$, the greater $H_{\rm
sh}$. As discussed in Section 3.2, such an anti-correlation is not always
present in the case of non-negligible cyclotron radiation. In particular,
there is a critical specific accretion rate above which there is an anti-
correlation and below which there is a correlation. Such a critical specific
accretion rate strongly depends on the WD mass and the WD magnetic field.

Figure 21: Comparison of two different cyclops X-ray spectra with Suleimanov
et al. (2016). In both panels, the spectra are normalized by their values at
$1$ keV. Our spectra are in very good agreement with Suleimanov et al. (2016)
since the differences are very small for energies below ${\sim 20}$ keV
($\lesssim 10$ %) and still small for energies greater than that (${\lesssim
30}$ %). Such small difference between both works are likely due to slightly
different bremsstrahlung prescriptions as well as PSR modeling, which have
more impact on the hard part of the spectrum.

Now, moving to the X-ray modeling, we show in Fig. 21 two X-ray spectra built
with parameters consistent with Suleimanov et al. (2016). In both cases, we
assumed negligible cyclotron cooling, ${\dot{m}=1\leavevmode\nobreak\ {\rm
g\,s^{-1}\,cm^{-2}}}$, and ${\mu=0.615}$. In addition, while in the left panel
we adopted ${R_{\rm th}=2\,R_{\rm WD}}$ and ${M_{\rm
WD}=0.6\leavevmode\nobreak\ {\rm M}_{\odot}}$, in the right panel we set
${R_{\rm th}=1000\,R_{\rm WD}}$ and ${M_{\rm WD}=0.9\leavevmode\nobreak\ {\rm
M}_{\odot}}$. Notice that our X-ray spectra are in good agreement with the
results by Suleimanov et al. (2016). In particular, the spectra agree
reasonably well for energies below ${\sim 20}$ keV since differences are
usually smaller than ${\sim 10}$ %. However, for energies greater than that,
differences can be a bit higher (${\sim 20-30}$ %).

This is likely a consequence of the different PSR modeling and bremsstrahlung
prescription in both cases (compare our Section 2 and Appendix A with section
3 in Suleimanov et al. 2016), which affect more the hard part than the soft
one. For instance, the boundary value problem in Suleimanov et al. (2016, see
their eq. 21) is slightly different from ours, which follows Hayashi & Ishida
(2014a). In addition, their cooling function is computed by apec using the
AtomDB atomic database, which is different from our approach (see Equations
A17, A18 and A20). Moreover, the calculation of the X-ray spectrum in
Suleimanov et al. (2016, see their eq. 29) is different from what is done in
the cyclops code (see section 2.2 in Silva et al., 2013). Despite the fact
that there are differences in the PSR modeling and bremsstrahlung emission in
both works, provided that the differences between cyclops and Suleimanov et
al. (2016) are small in most portions of the X-ray energy range, we can
conclude that both results are consistent with one another.

##  Appendix C Comparison with xspec

Figure 22: Comparison between X-ray spectra from cyclops and xspec (model
bremss), for four different values of the temperature, namely $1$ (top left
panel), $10$ (top right panel), $50$ (bottom left panel) and $100$ keV (bottom
right panel). In all models, the electronic density is $10^{16}$ cm-3 and the
distance is $100$ pc. Clearly, the spectra from cyclops are in excellent
agreement with those from xspec, provided that the differences are usually
smaller than $10$ %, irrespective of the energy. We can then conclude that the
cyclops code properly deals with non-normalized spectra, i.e., it can take
into account the distance while fitting an object.

Despite the fact that the PSR modeling in the cyclops code includes the
solution of the stationary one-dimensional hydro-thermodynamic differential
equations describing the accreting plasma as described in Appendix A, it can
also handle PSRs with uniform temperature and density distributions, i.e.,
assuming that the temperature and densities are the same in the entire PSR.
This allows us to compare its bremsstrahlung emission modeling with that of
the xspec code (Arnaud, 1996; Dorman & Arnaud, 2001). Such a comparison is an
essential step in order to demonstrate that the distance-dependent flux in the
cyclops code is correctly calculated and the interstellar extinction is
properly taken into account.

With the goal of showing that cyclops provides consistent results, throughout
this section we compare it with xspec for several cases of homogeneous PSRs.
Before that, we shall emphasize that cyclops and xspec can only be compared in
the optically thin regime, because this is the only case currently implemented
in the xspec code, although for X-rays this assumption is not problematic. The
cyclops code, on the other hand, solves the radiative transfer without any
approximation in this regard, since it also deals with the usually optically
thick cyclotron emission, which happens in optical wavelengths.

In the following comparison, we utilized the model bremss (Karzas & Latter,
1961) in the xspec package. Even though a slightly different bremsstrahlung
prescription is adopted in this model with respect to cyclops, it corresponds
to the most suitable model to perform such a comparison. Other models such as
mekal (Mewe et al., 1985, 1986; Kaastra, 1992; Liedahl et al., 1995) have
bremsstrahlung recipes more similar to that of cyclops and could in principle
be better for our purpose. However, these models include emission lines of
many elements, which would not result in meaningful comparisons, provided that
currently there is no prescription for line emission in the cyclops code.
While discussing interstellar extinction, we adopted the phabs model in xspec,
which calculates the required photoelectric absorption cross-sections. In
addition, since our primary goal is to show that not only the shape but also
the flux of the cyclops spectra are consistent with those of xspec, we assume
in all cases a distance of $100$ pc. Moreover, we fix the electronic density
to $10^{16}$ cm-3 and focus here on the impact of the PSR temperature in
shaping the X-ray spectra as well as in determining their X-ray flux.

In order to properly compare cyclops with xspec, we need first to make sure
that the emitting source is the same in both cases. For that we utilize the
normalization factor for the bremss model adopted in xspec given by

${\rm Norm}\ =\ \frac{3.02\times
10^{-15}}{4\,\pi\,D^{2}}\,\int\,n_{e}\,n_{I}\,{\rm d}V$ (C1)

where $D$ is the distance to the source, $n_{e}$, $n_{I}$ are the electron and
ion number densities, and $V$ is the volume. For our purposes here, we assume
that $n_{e}=n_{I}$, which implies that Eq. C1 can be expressed as

${\rm Norm}\ \simeq\ 2.524\times 10^{-53}\,\left(\frac{D}{\rm
cm}\right)^{-2}\,\left(\frac{n_{e}}{\rm
cm^{-3}}\right)^{2}\,\left(\frac{S_{\rm b}}{\rm
cm^{2}}\right)\,\left(\frac{H_{\rm sh}}{\rm cm}\right)$ (C2)

where $S_{\rm b}$ is the accretion area at the PSR bottom and $H_{\rm sh}$ is
the PSR height (see Fig. 2). This normalization factor directly implies that
the ‘observed’ bremsstrahlung flux coming out of a model depends not only on
the distance but also on the PSR properties.

In order to compare the spectra from cyclops and xspec, we first calculated
the PSR of a generic model with the formalism described in Appendix A. That
procedure provides a non-homogeneous PSR, i.e., the temperature and density
vary along the PSR. However, as mentioned before in this section, we assume
here that the temperature and density are constant throughout the PSR. To
achieve that, we kept the PSR geometrical properties and simply forced the
temperature and density to be constant. Thus, having the PSR cross-section and
height from this generic model and fixing the distance, we can compute the
normalization factor with Eq. C2. This normalization factor is then used to
make sure we have the same emitting region properties, i.e., volume and
distance, in cyclops and xspec, guarantying in turn that we can compare the
fluxes. Finally, the xspec spectra have been calculated using the standard
procedure for the bremss model.

We show in Fig. 22 the X-ray spectra from cyclops and xspec, considering four
different temperatures, namely $1$, $10$, $50$ and $100$ keV. We can clearly
see from the figure that the differences between the bremss model in xspec and
cyclops are usually smaller than $10$ %, irrespective of the energy and
assumed temperature. Given that cyclops and xspec are different codes, it is
not surprising that their spectra are not identical. However, they are very
similar, being consistent with one another. We can then conclude that the
cyclops code provides X-ray spectra in excellent agreement with the widely
used xspec code. Moreover, the precision of the analytical expressions of the
free-free Gaunt factor adopted in xspec and cyclops is around $10$ % (see,
e.g., Mewe et al., 1986). So the differences between the two codes are within
the precision we have in the bremsstrahlung emissivities.

We would like to draw the readers attention to the fact that cyclops X-ray
flux is independent of the assumed spatial resolution in the plane of the sky,
which implies that the area of each and every cell in the 3D grid is correctly
calculated. This is because such a resolution defines the size of each cell in
the grid and therefore its solid angle, which also depends on the distance. In
addition, to obtain the flux in each line-of-sight that represents the PSR,
the emerging intensity from the cell closer to the observer, which results
from the radiative transport along the entire line-of-sight, is multiplied by
the cell area (in steradian). Therefore, since the flux is not affected by
variations in the spatial resolution, we can conclude that the angular areas
of the cells are correctly calculated in cyclops.

Figure 23: Comparison between X-ray spectra from cyclops and xspec
(combination of models bremss and phabs), taking into account photoelectric
absorption, for three different values of the hydrogen column density $N_{\rm
H}$, namely $10^{20}$ (left panel), $10^{21}$ (middle panel), and $10^{22}$
cm-2 (right panel). In all models, the temperature is $20$ keV, the electronic
density is $10^{16}$ cm-3, and the distance is $100$ pc. As expected,
absorption is negligible for low values of $N_{\rm H}$, but becomes relevant
for ${N_{\rm H}>10^{20}}$ cm-2, at energies below a few keV. Most importantly,
the impact of absorption in cyclops and xspec is very similar (differences are
smaller than $10$ %), except for energies below $1$ keV. We can then conclude
that both codes are consistent with one another.

The last aspect of the X-ray modeling we compare with xspec is the
interstellar extinction. We show in Fig. 23 how the photoelectric absorption
affects the spectra, by assuming a temperature of $20$ keV and three different
values of hydrogen column density $N_{\rm H}$, namely $10^{20}$, $10^{21}$,
and $10^{22}$ cm-2. We can clearly see that the impact of absorption in both
codes is rather similar, irrespective of $N_{\rm H}$. On the other hand, for
energies below $1$ keV, models in the cyclops code are clearly much more
absorbed in comparison with what the combination of phabs and bremss models
provide. In any event, like in the previous case in which we varied the
temperature, differences between both codes are smaller than $10$ %, except
for the above-mentioned discrepancy at the softer part. Therefore, as before,
we can conclude that spectra in cyclops and xspec are in excellent agreement.