© The Authors 2024

1. Introduction

Low-mass X-ray binaries hosting a weakly magnetized neutron star (NS-LMXBs) are among the brightest sources of the X-ray sky and have been studied since the dawning of X-ray astronomy. In the Milky Way, about 150 NS-LMXBs are known (Avakyan et al. 2023). They have historically been divided into two main classes – Z and atoll sources – on the basis of the shape they track in color–color or color–intensity diagrams (Hasinger & van der Klis 1989). The main characteristics of Z sources are a persistently bright luminosity (L_X ≳ L_Edd) and a relatively stable spectrum that changes slightly while the source moves along the Z track. The X-ray spectrum is usually modeled in terms of a soft thermal component, plus an unsaturated Comptonization spectrum of ∼1 keV black-body-like photons off a thermal electron distribution with kT_e  ∼  3–5 keV in a high optical depth environment and τ  ∼  5–10 depending on the assumed geometry, slab, or spherical, respectively (Di Salvo et al. 2000, 2001; D’Amico et al. 2001; Lavagetto et al. 2004; Farinelli et al. 2008, 2009; Homan et al. 2018). Additional features are the presence of an iron fluorescence line, as well as a transient, hard X-ray tail above 30 keV, which appears correlated to the radio activity and in turn the formation of a jet (Paizis et al. 2006). On the other side, atoll sources cover a wider range of luminosities (from ≲10⁻³ up to ∼L_Edd), sharing in this fashion a common property with black-hole X-ray binaries (Done et al. 2007). The changes in luminosity are linked to variations of the spectral state, with the brightest luminosity achieved when the X-ray spectrum is Z-like (soft state), and dimmer fluxes (hard state) in which the spectrum is characterized by kT_e ∼ 30 keV and τ∼ a few (Gierliński & Done 2002; Maccarone & Coppi 2003; Falanga et al. 2006; Cocchi et al. 2011). There is also a subclass of bright atoll sources (GX 9+9, GX 9+1, GX 3+1, and GX 13+1) that have never shown spectral transition, with a rather stable X-ray spectrum typical of Z sources (Mainardi et al. 2010; Iaria et al. 2020).

For a long time, understanding the shape of the accretion process in NS-LMXBs has been a challenge: in particular, it has long been discussed about the hydrodynamical configuration of the inner region responsible for the strong, dominating, thermal Comptonization bump. It has been proposed as either a diffuse accretion disk corona (White & Holt 1982) or the region between the inner accretion disk and the NS surface (Frank et al. 2002). In the latter case, in the literature it is often defined as the boundary layer (BL); albeit, in some cases it is also called spreading layer (Inogamov & Sunyaev 1999; Suleimanov & Poutanen 2006).

The origin of the direct thermal emission at lower energies has been a matter of debate as well. Historically, this has led to spectral modeling degeneracy, in which the observed X-ray spectrum could be equally well described by an unsaturated Comptonization component plus either a direct black-body (BB) emission, attributed to the NS surface (Western Model, White et al. 1988a), or a direct contribution from the accretion disk (Eastern Model, Mitsuda et al. 1989). Somewhat surprisingly, even missions with unprecedented broad-band coverage such as BeppoSAX (Boella et al. 1997) have not been able to remove the degeneracy. A strong boost in favor of the Eastern Model scenario has been obtained from the analysis of a sample of atoll and Z sources using Fourier frequency-resolved spectroscopy of archival data of the Rossi X-Ray Timing Explorer (Gilfanov et al. 2003; Revnivtsev & Gilfanov 2006; Revnivtsev et al. 2013). In these works, the authors have shown that the unsaturated Comptonization region can be reasonably described in terms of a spreading layer covering the NS surface in a belt-like fashion, with a power spectrum resembling the Fourier frequency-resolved spectrum at the period of quasi-periodic oscillations.

Significant expectations for a major leap in knowledge were put on the Imaging X-ray Polarimetry Explorer (IXPE, Weisskopf et al. 2022), because the energy range of the spacecraft (2–8 keV) covers a substantial fraction of the emission bulk in NS-LMXBs, in particular in the region where the soft and hard components overlap. Indeed, X-ray polarimetry does constitute the effective third axis of an ideal vector in which the other two components are the spectral and temporal analysis, respectively, and it can be fundamental in removing model degeneracy. Joint observations of IXPE with other X-ray facilities covering a wider energy range have shown how spectropolarimetric information can provide unprecedented constraints on the inner region of accreting X-ray binaries (see e.g., Krawczynski et al. 2022, for the black hole Cyg X–1). A sample of NS-LMXBs in soft state have been observed so far by IXPE; however, despite some very intriguing results, the data statistics have not allowed us to put unambiguously tight constraints on the polarization degree (PD) and polarization angle (PA) of the different components that contribute to the total observed spectrum.

Different values of the PD of the BL component have been observed when performing model-dependent spectropolarimetric analysis. Concerning Z sources, for Cyg X-2 Farinelli et al. (2023) found a value of ≈4 ± 2% during a possible normal branch (NB) phase, while when studying GX 5–1 Fabiani et al. (2024) obtained 5.7 ± 1.4% in the horizontal branch (HB) and 4.3 ± 2.0% in the state time-integrated over the NB and flaring branch (FB). However, in the second case the source showed a peculiar behavior with the presence of a direct BB component plus a Comptonization feature. For Sco X-1, the PD obtained by La Monaca et al. (2024) was 1.3 ± 0.4% when the source was in the soft apex joining the NB with the FB. Studying the outburst of the peculiar source XTE J1701–462, Cocchi et al. (2023) modeled the IXPE-alone X-ray spectrum with a simple BB + multicolor BB disk model and derived a PD of ≈6 ± 0.5% for the hard BB component when the source was in the HB, the value dropping to ≈1.2 ± 0.7% in the NB.

Atoll sources in a soft state have been investigated as well. In the case of the bright source GX 9+9 (Ursini et al. 2023), the presence of a reflection component, necessary from the spectral point of view, added further degeneracy in the spectropolarimetric analysis. Among different combinations, the authors fixed the PD of the reflection to 10%, obtaining 3 ± 1% for the BL. For the same component in 4U 1820–303, Di Marco et al. (2023) obtained 5.3 ± 0.2%, while in the case of GS 1826−238, because of the intrinsically low polarization level of the source (≲1.3%), it was not possible to perform a detailed spectropolarimetric study (Capitanio et al. 2023).

2. Modeling the polarization of the boundary layer around NS

The continuously increasing amount of data coming from IXPE, as well as future satellites such as the enhanced X-ray Timing and Polarimetry mission (Zhang et al. 2019), needs a parallel development of theoretical and numerical simulations, which are fundamental for data interpretation. Several codes for computing polarization in the X-ray band have been developed both in the past and more recently. The cases ignoring general relativity effects were treated, for instance, by Sunyaev & Titarchuk (1985, hereafter ST85) and Poutanen & Svensson (1996, hereafter PS96), who considered a slab geometry. Both algorithms were based on an iterative procedure, the main difference consisting of the use of the Thomson (ST85) or Klein-Nishina (PS96) scattering kernel. Most efforts in the theoretical investigation have so far been focused on the polarization properties of accretion disks around black holes (Schnittman & Krolik 2009, 2013; Dovčiak et al. 2011; Abarr & Krawczynski 2020; Taverna et al. 2021; Loktev et al. 2022; Podgorný et al. 2022).

Instead, X-ray polarimetry of XRBs hosting an NS has received much less attention, and in this framework the seminal paper of Lapidus & Sunyaev (1985, hereafter LS85) deserves to be mentioned, as the authors considered the polarization properties of X-ray bursters modeled with an accretion disk and a BL approximately at a right angle, both during the burst and the out-of-burst phases. Despite being a work dating back to the mid-1980s, to our knowledge this is the most detailed treatment of the problem so far. For our work, it is important to point out that LS85 modeled the BL as an equatorial belt around the NS of geometrical thickness ΔR ≪ R_ns and variable height-scale ratio H/R_ns. In particular, LS85 did not consider the direct disk emission and showed that the total PD of radiation directly coming from the BL plus the fraction intercepted by the disk and reflected toward the observer may achieve 6% for a viewing angle around 70°. The results of LS85 do not consider light bending of radiation from the emission region to the observer frame. The problem was instead considered to study spectropolarimetric properties of X-ray pulsars (XRPs), in which the local emission comes from the hot spots at the magnetic poles (e.g., Pihajoki et al. 2018; Bogdanov et al. 2019).

The aim of our work is to investigate the polarization properties of the BL using an approach similar to that used for XRPs, but considering the phase-independent problem in which the BL covers the NS surface in a belt-like manner. We adopted a two-phase approach: in Sect. 2.1, we computed the PD and PA of radiation in the BL comoving frame through the solution of the radiative transfer equation (RTE) for polarized radiation in plane-parallel geometry. As long as the radial thickness of the BL is significantly less than the NS radius, this is a good approximation. Subsequently, we performed angular integration over the BL surface and calculated the total received signal in the observer frame as a function of inclination with respect to the NS spin axis (see Sect. 2). We present the numeral results in Sect. 3 and provide a discussion in Sect. 4. Finally, we draw our conclusions in Sect. 5.

2.1. The radiative transfer equation of polarized radiation in plane-parallel geometry

We considered plane-parallel geometry and define the total specific intensity as I^tot = I_l + I_r, where I_l and I_r are the intensities with the electric field contained in the meridional plane and perpendicular to it, respectively (Chandrasekhar 1960, hereafter C60). The meridional plane is defined as the plane formed by the normal to the slab and the propagation direction of photons.

The coupled RTEs for an electron-scattering atmosphere in the Thomson limit are (Pomraning 1973)

$$\begin{array}{cc}\hfill \mu \frac{\mathrm{d}{I}_l(\tau ,\mu )}{\mathrm{d}\tau }& =-I_l(\tau ,\mu )+\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{1,1}(\mu ,{\mu }^{\prime })I_l(\tau ,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{1,2}(\mu ,{\mu }^{\prime })I_r(\tau ,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime }+S_l(\tau ,\mu ),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \mu \frac{\mathrm{d}{I}_r(\tau ,\mu )}{\mathrm{d}\tau }& =-I_r(\tau ,\mu )+\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{2,1}(\mu ,{\mu }^{\prime })I_l(\tau ,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{2,2}(\mu ,{\mu }^{\prime })I_r({\mu }^{\prime },\tau )\mathrm{d}{\mu }^{\prime }+S_r(\tau ,\mu ),\hfill \end{array}$$(2)

where μ is the cosine of the angle of radiation beam with respect to the normal and S_l, r(τ) is the source term describing the distribution over τ of the seed photons. We considered initially unpolarized seed photons; thus, $S_l(\tau )=S_r(\tau )=\frac{1}{2}S(\tau )$.

The Chandrasekhar scattering matrix 𝒫_i, j(μ, μ′) is defined as (ST85) follows:

$$\begin{array}{cc}& {\mathcal{P}}_{1,1}(\mu ,{\mu }^{\prime })=2(1-{\mu }^2)(1-{\mu }^{\prime 2}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\mathcal{P}}_{1,2}(\mu ,{\mu }^{\prime })={\mu }^2,\hfill \end{array}$$(4)

$$\begin{array}{cc}& {\mathcal{P}}_{2,1}(\mu ,{\mu }^{\prime })={\mu }^{\prime 2},\hfill \end{array}$$(5)

$$\begin{array}{cc}& {\mathcal{P}}_{2,2}(\mu ,{\mu }^{\prime })=1.\hfill \end{array}$$(6)

The system of two integro-differential equations – Eqs. (1) and (2) – for the unknowns I_l(τ, μ),I_r(τ, μ) is solved with an iterative procedure with appropriate boundary conditions, similarly to ST85. First, we consider the distribution over τ of photons which escape without interaction with matter. Mathematically, this is equivalent to ignoring the scattering terms (i.e., the integrals), leaving only the source function S(τ) and the absorption terms in the right sides of the system. We note that this is not a true absorption in terms of photon loss due to bound-bound or bound-free processes, but rather a decrease of the photon intensity in the given direction, μ, as Thomson scattering diffuses a fraction of the incoming intensity away from the μ direction.

For the unscattered radiation, Eqs. (1) and (2) then become

$$\begin{array}{cc}& \mu \frac{\mathrm{d}{I}_l(\tau ,\mu )}{\mathrm{d}\tau }=-I_l(\tau ,\mu )+\frac{S(\tau )}{2},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \mu \frac{\mathrm{d}{I}_r(\tau ,\mu )}{\mathrm{d}\tau }=-I_r(\tau ,\mu )+\frac{S(\tau )}{2},\hfill \end{array}$$(8)

for which the general solution is

$$\begin{array}{c}\hfill I_{l,r}^0(\tau ,\mu )=C_1{\mathrm{e}}^{-\tau /\mu }+\frac{1}{2\mu }{\displaystyle {\int }_0^{\tau }}{\mathrm{e}}^{-(\tau -{\tau }^{\prime })/\mu }S({\tau }^{\prime })\mathrm{d}{\tau }^{\prime }.\end{array}$$(9)

We define the bottom of the slab at τ = 0 and the top at τ = 2τ₀. We note that we conventionally define the total optical depth of the slab as 2τ₀ to keep the same definition of ST85. The condition of no incident radiation at the base of the slab reads I_l, r(0, μ)=0, which implies C₁ = 0. The explicit form of Eq. (9) depends on the initial seed photon distribution S(τ).

We consider two cases here. In the first one, which is labeled Case 1 and is important for our further discussion, seed photons have isotropic angular distribution for μ ≥ 0 and are located at the base of the atmosphere, namely S(τ)=S₀δ(τ). We note that with this configuration S(τ)=0 for μ <  0, and in this case Eq. (9) becomes

$$\begin{array}{cc}& I_{l,r}^0(\tau ,\mu )=\frac{S_0}{2\mu }{\mathrm{e}}^{-\tau /\mu }\mu \ge 0,\hfill \end{array}$$(10)

$$\begin{array}{cc}& I_{l,r}^0(\tau ,\mu )=0\mu <0.\hfill \end{array}$$(11)

For the second case, labeled Case 2, we consider isotropic seed photons uniformly distributed over the slab, namely S(τ)=S₀, and in this case the solution of the RTE gives

$$\begin{array}{cc}& I_{l,r}^0(\tau ,\mu )=\frac{S_0}{2}(1-{\mathrm{e}}^{-\tau /\mu })\mu \ge 0,\hfill \end{array}$$(12)

$$\begin{array}{cc}& I_{l,r}^0(\tau ,\mu )=\frac{S_0}{2}[1-{\mathrm{e}}^{-(2{\tau }_0-\tau )/\mu }]\mu <0.\hfill \end{array}$$(13)

2.2. Solution for scattering order k >  0

When considering the specific intensity of radiation subjected to k scatterings, the source term S(τ) disappears, and the contribution to the specific intensity across the μ direction comes from photons of scattering order $I_{l,r}^{k-1}$ diffused into the μ direction from all other directions. This contribution arises from both photons, which remain in the same polarization mode (e.g., $I_l^{k-1}\to I_l^k$), as well as those switching from one mode to the other (e.g., $I_r^{k-1}\to I_l^k$).

The system of Eqs. (1) and (2) then becomes

$$\begin{array}{cc}\hfill \mu \frac{\mathrm{d}{I}_l^k(\tau ,\mu )}{\mathrm{d}\tau }& =-I_l^k(\tau ,\mu )+\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{1,1}(\mu ,{\mu }^{\prime })I_l^{k-1}({\mu }^{\prime })\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{1,2}(\mu ,{\mu }^{\prime })I_r^{k-1}({\mu }^{\prime })\mathrm{d}{\mu }^{\prime },\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill \mu \frac{\mathrm{d}{I}_r^k(\tau ,\mu )}{\mathrm{d}\tau }& =-I_r^k(\tau ,\mu )+\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{2,1}(\mu ,{\mu }^{\prime })I_l^{k-1}({\mu }^{\prime })\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{2,2}(\mu ,{\mu }^{\prime })I_r^{k-1}({\mu }^{\prime })\mathrm{d}{\mu }^{\prime }.\hfill \end{array}$$(15)

As the angular integration is performed over the domain μ′∈[−1 : 1], we define the upstream intensity I^( + )(τ, μ) as the one directed toward the top layer of the slab, and downstream intensity I^( − )(τ, μ) as that directed toward the bottom layer. In this way, the variable μ is positively defined for both I^( + )(τ, μ) and I^( − )(τ, μ), from bottom to top in the first case, and the opposite in the second one (see Fig. 1).

Fig. 1. Sketch of radiative transfer problem in plane-parallel geometry. Seed photons at the base of the slab are shown here. At any location τ, a positive angle μ relative to the slab normal can be defined for an upward and downward intensity component, respectively.

Considering, for example, the l-mode, the integral terms can thus be written as

$$\begin{array}{cc}& {\displaystyle {\int }_{-1}^1}{\mathcal{P}}_{i,j}(\mu ,{\mu }^{\prime })I_l^{k-1}(\tau ,{\mu }^{\prime })\mathrm{d}{\mu }^{\prime }=\hfill \\ \hfill & {\displaystyle {\int }_0^1}{\mathcal{P}}_{i,j}(\mu ,{\mu }^{\prime })[I_l^{(+)k-1}(\tau ,{\mu }^{\prime })+I_l^{(-)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime }.\hfill \end{array}$$(16)

Given that the intensities I_r, l(τ, μ) are now expressed as the contribution of the upstream and downstream factors, the coupled RTEs switch from two to four integro-differential equations according to

$$\begin{array}{cc}& \mu \frac{\mathrm{d}{I}_l^{(+)k}(\tau ,\mu )}{\mathrm{d}\tau }=-I_l^{(+)k}(\tau ,\mu )\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{1,1}(\mu ,{\mu }^{\prime })[I_l^{(+)k-1}(\tau ,{\mu }^{\prime })+I_l^{(-)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{1,2}(\mu ,{\mu }^{\prime })[I_r^{(+)k-1}(\tau ,{\mu }^{\prime })+I_r^{(-)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime },\hfill \end{array}$$(17)

$$\begin{array}{cc}& \mu \frac{\mathrm{d}{I}_r^{(+)k}(\tau ,\mu )}{\mathrm{d}\tau }=-I_r^{(+)k}(\tau ,\mu )\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{2,1}(\mu ,{\mu }^{\prime })[I_l^{(+)k-1}(\tau ,{\mu }^{\prime })+I_l^{(-)}(2{\tau }_0-\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime }\hfill & \hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{2,2}(\mu ,{\mu }^{\prime })[I_r^{(+)k-1}(\tau ,{\mu }^{\prime })+I_r^{(-)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime },\hfill \end{array}$$(18)

$$\begin{array}{cc}& \mu \frac{\mathrm{d}{I}_l^{(-)k}(\tau ,\mu )}{\mathrm{d}\tau }=-I_l^{(-)k}(\tau ,\mu )\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{1,1}(\mu ,{\mu }^{\prime })[I_l^{(+)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })+I_l^{(-)k-1}(\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{1,2}(\mu ,{\mu }^{\prime })[I_r^{(+)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })+I_r^{(-)k-1}(\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime },\hfill \end{array}$$(19)

$$\begin{array}{cc}& \mu \frac{\mathrm{d}{I}_r^{(-)k}(\tau ,\mu )}{\mathrm{d}\tau }=-I_r^{(-)k}(\tau ,\mu )\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{2,1}(\mu ,{\mu }^{\prime })[I_l^{(+)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })+I_l^{(-)k-1}(\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime }\hfill \\ \hfill & +\frac{3}{8}{\displaystyle {\int }_0^1}{\mathcal{P}}_{2,2}(\mu ,{\mu }^{\prime })[I_r^{(+)k-1}(2{\tau }_0-\tau ,{\mu }^{\prime })+I_r^{(-)k-1}(\tau ,{\mu }^{\prime })]\mathrm{d}{\mu }^{\prime }.\hfill \end{array}$$(20)

The first two equations of the system solve for the upstream intensities $I_l^{(+)k}$ and $I_r^{(+)k}$ with τ computed from bottom to top, while the other two solve for the downstream components $I_l^{(-)k}$ and $I_r^{(-)k}$ with τ computed from top to bottom. Once defined the seed photons distribution S(τ), the algorithm first seeks for the solution $I_{l,r}^0(\tau ,\mu )$ as defined in Eqs. (10)–(11) or (12)–(13), and then iteratively integrates over τ equations up to a user-defined maximum number of scatterings k^max.

The values of the specific intensity in Eqs. (17)–(20) are computed at fixed values of μ_i, which are then preliminarily defined over a given grid, and the angular integration over μ′ is performed using the Gauss-Legendre (GL) quadrature rule. We remind the reader that GL quadrature is used for computing the integral of functions bounded over the domain x ∈ [a : b] according to

$$\begin{array}{c}\hfill {\displaystyle {\int }_a^b}f(x)\mathrm{d}x\approx {\displaystyle \sum _{i=0}^{N-1}}f(x_i)w_i,\end{array}$$(21)

where a <  x_i <  b are the roots of the Nth Legendre polynomial, and w(x_i) are the quadrature weights. Of course, the higher the value of N, the greater the precision in approximating the true integral value. We checked that for N ≳ 20 there are no further appreciable changes in numerical results, so we fixed N = 30. We eventually found that it was convenient to compute the solution at the μ_i values coincident with the roots of the GL quadrature. The numerical integration over τ was performed using an explicit embedded Runge-Kutta Prince-Dormand method with adaptive step-size control provided by the Gnu Scientific Library¹.

By definition, the iterative method provides the solution for each scattering order $I_{r,l}^k(\tau ,\mu )$. We are interested in the following quantities at the top of the slab (τ = 2τ₀), where it is assumed that an observer is located: angular distribution of the total specific intensity (limb-darkening law); angular distribution of the total polarization.

Defining the total intensities of the two polarization modes (r, l) at the top of the slab by summing over all the orders of scattering, we obtain

$$\begin{array}{c}\hfill I_{r,l}^{(+)\mathrm{tot}}(\mu )={\displaystyle \sum _{k=0}^n}{I}_{\mathrm{r},\mathrm{l}}^{(+)k}(\mu ,2{\tau }_0),\end{array}$$(22)

and the total polarization comes directly to

$$\begin{array}{c}\hfill P^{\mathrm{tot}}(\mu )=\frac{I_r^{(+)\mathrm{tot}}(\mu )-I_l^{(+)\mathrm{tot}}(\mu )}{I_r^{(+)\mathrm{tot}}(\mu )+I_l^{(+)\mathrm{tot}}(\mu )}.\end{array}$$(23)

Dropping the “tot” superscript for clarity, one may see that P(μ) can be positive or negative. As in planar and spherical configurations, the Stokes parameter U = 0; here, P = Q. Keeping in mind the definition of the PA,

$$\begin{array}{c}\hfill \chi =\frac{1}{2}\mathrm{Atan}\frac{U}{Q},\end{array}$$(24)

it results that for Q >  0, χ = 0; otherwise, χ = π/2. The PA is then parallel or perpendicular to the slab plane for positive and negative values of the PD, respectively.

We also considered the percentage of the flux emerging from the top of the slab, where the comoving observer is located. From the definition of flux and using the GL quadrature rule, one obtains

$$\begin{array}{cc}& F^{(+)}={\displaystyle \sum _{i=0}^{N-1}}{I}_{r,l}^{(+)\mathrm{tot}}({\mu }_i){\mu }_{i}w({\mu }_i),\hfill \end{array}$$(25)

$$\begin{array}{cc}& F^{(-)}={\displaystyle \sum _{i=0}^{N-1}}{I}_{r,l}^{(-)\mathrm{tot}}({\mu }_i){\mu }_{i}w({\mu }_i),\hfill \end{array}$$(26)

$$\begin{array}{cc}& ϵ(2{\tau }_0)=\frac{F^{(+)}}{F^{(+)}+F^{(-)}},\hfill \end{array}$$(27)

where F^( + ) and F^( − ) are the upstream (τ = 2τ₀) and downstream fluxes (τ = 0), respectively.

The parameter ϵ(2τ₀) is by definition equal to 0.5 for seed photons distribution symmetric with respect to the center of the slab (as in ST85), while it depends on the optical depth for seed photons located at the base.

2.3. Polarimetry of the boundary layer

The PD and PA in Eqs. (23) and (24) are defined at the top of the slab, which is identified as the surface of the BL (i.e., comoving frame). As already mentioned in Sect. 1, the mathematical formalism for passing from the comoving to the observer frame has been developed to investigate spectropolarimetric properties of accreting XRPs. We consider the method defined in Poutanen (2020, hereafter P20) to which we refer the reader for details. The most significant difference between our work and results obtained for XRPs is the change in the contribution to the total observed PD and PA from just two small spots at the magnetic poles of the NS to an equatorial belt covering the NS surface up to a given latitude θ^max. This can be achieved thanks to the additive property of the Stokes parameters, so the computed quantities can be derived as the contribution from each elementary area dΣ′ of the BL surface integrated over the domain θ ∈ [π/2 − θ^max : π/2] and ϕ ∈ [0 : 2π] (see Fig. 2 and Eqs. (17)–(21) in P20). The proper area element using angular coordinates is

$$\begin{array}{c}\hfill d{\mathrm{\Sigma }}^{\prime }=\gamma (\theta )R\sin \theta \mathrm{d}\theta \mathrm{d}\varphi ,\end{array}$$(28)

Fig. 2. Proposed geometry for soft state of NS-LMXBs; the disk extends very close to the NS surface, and then a BL forms with ΔR ≪ Rns. The BB photons from the NS photosphere are located at the base of the BL, which extends up to some latitude θmax computed from the orbital plane, and has a given temperature kTe and radial optical depth 2τ0. Only photons in the northern hemisphere (θ <  π/2) reach the observer. The total signal is obtained by integrating the contribution from each elementary area dΣ′ shown by the blue box over the BL surface.

where γ(θ) is the Lorentz factor of the matter velocity at given latitude with

$$\begin{array}{c}\hfill \beta (\theta )=\frac{2\pi R_{\mathrm{bl}}}{c}\frac{{\nu }^{\mathrm{bl}}}{\sqrt{1-\eta }}\sin \theta .\end{array}$$(29)

The impact parameter in Eq. (29) is η = R_g/R_bl, where R_bl ≳ R_ns, and for which we set the value η = 1/2.5, while we consider a BL co-rotating with the NS at frequency v^bl = 500 Hz.

For each elementary area dΣ′ at coordinates θ and ϕ on the BL surface, one defines the Stokes vector in the observer frame as

$$\begin{array}{c}\hfill S=I^{\mathrm{tot}}d\mathrm{\Omega }[1,P\cos 2\chi ,P\sin 2\chi ].\end{array}$$(30)

For the relation between dΩ and dΣ′ as well as I^tot and P between the observer and comoving frame, we refer the reader to Eqs. (17)–(20) of P20, while χ is given in Eq. (58) of the same paper and is the PA with respect to the projection of the NS spin onto the sky.

By defining [BOLD]𝒮^tot as the Stokes vector obtained integrating [BOLD]𝒮^tot over dΩ, the total PD and PA are then derived using the standard definitions, with Q and U as the second and third component of [BOLD]𝒮^tot, respectively (see Eq. (24)). It is important to note that the upper boundary for the polar angle is θ = π/2, because only radiation from the BL northern hemisphere is received, the southern part being totally intercepted by the inner disk (and vice versa for an observer located in the opposite position).

While performing two-dimensional integration over the BL surface, the comoving frame angular dependence of the specific intensity and PD are obtained by a linear interpolation of the array of values stored in a text file and preliminarily computed for plane-parallel geometry. Numerical results were obtained with a C code optimized with a third-party adaptive integrator².

3. Results

We first discuss the results obtained from the solution of the RTE in plane-parallel geometry (see Eqs. (1) and (2)). The solutions presented by ST85 considered three different configurations for the initial seed photons distribution, namely (i) sources at the center of the slab, (ii) uniform photon distribution, (iii) and photons at the boundaries. These configurations are by definition totally symmetric with respect to the center of the slab, and there is no difference while considering physical quantities in the upward or downward direction. We consider a uniform seed photon distribution as well, but unlike the ST85 case we set up a geometry where seed photons are located at the bottom of the slab. In this case, the upward and downward symmetry breaks and different spectropolarimetric values are obtained at τ = 0 or τ = 2τ₀.

We focus on the results at the top of the slab, where we assume that the comoving observer is located. A configuration with a plane-parallel electron-scattering atmosphere located above a source of photons can be representative of a geometrically thin BL partially covering the surface of the NS, with ΔR ≪ R_ns, or the atmosphere above an accretion disk. The latter case was also recently investigated by Taverna et al. (2021) for both the cases of pure electron scattering and a combination of scattering plus absorption; although, the authors reported their results only in the rest frame of the atmosphere. The case of uniform seed photons distribution, on the other hand, is in our opinion more suitable to be considered if they originate from free-free emission processes (Narayan & Yi 1995). However, we present results concerning observational differences between the two configurations here.

In Fig. 3, we report the angular behavior of the total specific intensity and PD for seed photons located at the base of the slab (Case 1) and for different values of the optical depth. The specific intensity values are normalized to the case μ = 1. When 2τ₀ ≲ 1, the PD is always negative, which implies a PA aligned with the slab normal. For the intermediate value 2τ₀ = 2, PA swapping occurs around μ ∼ 0.3, while for 2τ₀ = 5 the solution becomes virtually indistinguishable from the Chandrasekhar limit. When seed photons are uniformly distributed, on the other hand, except for a small interval at μ ≲ 0.1 for 2τ₀ = 5, PD is always negative, as shown in Fig. 4. For low optical depths, with both configurations the total intensity I^tot(μ) does not show monotonic behavior (see also Fig. 4 of ST85); however, it should be kept in mind that the observed flux at large distances is F^tot(μ)∝I^tot(μ)μ, and in the case of an accretion disk atmosphere this mitigates the first-glance impression that at some more edge-on viewing angles the flux is higher than the face-on case (neglecting relativistic light bending effects).

Fig. 3. Angular distribution of total PD (upper panel) and normalized specific intensity (lower panel) emerging from the upper layer of the plane-parallel atmosphere and different values of the optical depth. Seed photons are located at the base of the slab (Case 1). For positive values of the PD, the PA is coplanar to the slab; for negative values it is parallel to the normal.

Fig. 4. Same as Fig. 3, but for seed photons uniformly distributed over τ (Case 2).

Fig. 5. Total PA (upper panel), PD (middle panel), and normalized intensity (lower panel) of the BL as a function of the viewing angle with seed photons at the base of the slab (Case 1) for an optical depth of 2τ0 = 5, spin frequency of νbl = 500 Hz, and two values of the BL latitude. The PA is measured with respect to the projection on the sky of the NS spin axis.

With spectropolarimetric quantities at hand in the comoving frame (i.e., the top of the slab), we computed the polarimetric parameters as well as the flux in the observer frame as a function of the viewing angle. The results with high optical depth deserve particular attention because this is a typical value of the Comptonization component of NS-LMXBs in the soft state (see references in Sect. 1).

We considered two maximum latitudes of the BL, namely θ^max = 45° and θ^max = 30°, with results for Case 1 and optical depth 2τ₀ = 5 shown in Fig. 5. As expected, the PD for θ^max = 45° is lower than that for θ^max = 30°, as a consequence of the qualitative fact that a less extended BL is intrinsically “spherically less symmetric” than a more extended one – actually a BL covering the whole NS surface would have, by definition, null polarization. It is also interesting to note that when θ^max = 45°, the observed flux monotonically increases by about 30%, while passing from an edge-on to a face-on view of the system, while for θ^max = 30° it achieves its maximum value at i ∼ 50°. The BL properties for lower values of 2τ₀ are instead shown in Fig. 6 for θ^max = 45°. Here, a PD value of less than 1% and thus comparable to that of Fig. 5 is obtained for 2τ₀ = 2, but the PA is 90° apart, thus making the two cases easily distinguishable from the observational point of view. The angular dependence of the observed flux is, on the other hand, qualitatively similar to the case of higher optical depth and θ^max = 45°.

Fig. 6. Same as in Fig. 5, but for different values of the optical depth and fixed BL latitude θmax = 45° and spin νbl = 500 Hz.

We then considered the configuration with seed photons uniformly distributed over the slab (Case 2), with results presented in Fig. 7. Here, the PD for 2τ₀ = 5 is comparable to that of Case 1, but again the PA is rotated by 90°, avoiding degeneracy in the geometrical setup of the input radiation. The flux hence increases monotonically while progressively decreasing the observer viewing angle. A shared property in all cases is that the PA is not strictly aligned to the NS spin, or at a right angle to it, due to polarization plane rotation being the effect of the BL angular velocity.

Fig. 7. Same as Fig. 5, but with seed photons uniformly distributed (Case 2) and different values of the optical depth of the BL, with a latitude of θmax = 45° and spin of νbl = 500 Hz.

As we assumed a BL corotating with the NS, the effect depends on the spin value that we set to ν^bl = 500 Hz. Lower values of ν^bl of course decrease the deviation of PA from the two orthogonal directions, leaving a shift of ≲2° due only to the light-bending effect for the case of nonrotating NS/BL. Regardless of PA rotation being the result of special and general relativistic effects, the corresponding PD for high optical depth of the BL is ≲0.5%.

Fig. 8. Fraction of flux emerging from the top of the slab in the comoving frame as a function of the optical depth (see Eqs. (25)–(27)) when seed photons are located at the base.

4. Discussion

The behavior of the physical quantities presented in this work must be considered in their integral form; namely, we consider the angular dependence of intensity and PD in the NS comoving frame (i.e., the top layer of the slab) obtained by summing over all orders of scatterings (see Eqs. (22) and (23)). Our choice is motivated by two main reasons, which we explain below. The first one is that linking the photon energy of the emerging spectrum to the number of scatterings experienced in the medium is not as trivial as it could appear: for an electron thermal distribution, the energy of a photon that escapes the medium after k scattering is (Rybicki & Lightman 1979)

$$\begin{array}{c}\hfill E_k\approx E_0{(1+\frac{4k{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2})}^k,\end{array}$$(31)

where E₀ is the initial energy. This is true as long as E₀ ≲ E_k ≪ 4kT_e; otherwise, the recoil effect is to be taken into account. For the typical electron temperature of NS-LMXBs in the soft state (kT_e ∼ 3–5 keV), the validity of Eq. (31) fails in the energy interval that partially overlaps with the IXPE bandwidth (2–8 keV) so that the identification $I_{l,r}^k(\mu )\equiv I_{l,r}(E_k)$ appears too rough an approximation. We note that ST85 presented results both for the energy-integrated spectrum, as in Eq. (22), and for what they call qualitatively “hard radiation spectrum”, namely for $I_{l,r}^k(\mu )$ with $k>\overline{k}$, the average number of scatterings. When discussing the relation between the number of scatterings and the energy, they were only able to quantify this in the power-law-like regime of the emerging spectrum. It is also worth pointing out that Eq. (31) holds on average, namely there may be photons that are subjected to k >  0 scattering and leave the medium with energy lower than the initial one. The contributions of photons leaving the medium with lower or higher energies than input one is dictated by the shape of Green’s function (GF), for which I(E)∝e^α + 3 for E <  E₀ and I(E)∝e^−α (up to E ∼ kT_e) for E >  E₀, and where the spectral index α depends on the temperature and optical depth (Sunyaev & Titarchuk 1980, hereafter ST80). The red wing of the GF is steeper than the blue wing; however, for an electron temperature of the order of a few kilo-electron-volts, the power-law regime of the GF (defined by Eq. (31)) is confined to a narrow band. In turn, the exact form of I(E) and PD(E) requires either Monte Carlo methods or an explicit dependence on energy of the RTE (see Eqs. (1) and (2)). The second reason is the fact that despite its challenging results, the data statistics of IXPE observations of NS-LMXBs do not allow thorough investigations of detailed energy-dependent trends for PD and PA for each spectral component. When the behavior of PD and PA as a function of energy is presented, it indeed refers to a model-independent data analysis (see references in Sect. 1), in which the polarimetric quantities are the results of contribution from all components of the systems (Comptonization region, accretion disk, as well as disk reflection). For model-dependent spectropolarimetric analysis instead, the PD and PA of the single components are given in the whole IXPE band (2–8 keV). We thus claim that presenting our energy-integrated theoretical results (i.e., summed over all scattering orders k) is sufficient to provide theoretical support to IXPE model-dependent data interpretation.

One potential objection to our results could arise from the fact that we used the RTE of polarized radiation in the Thomson approximation (i.e., the Chandrasekhar scattering matrix) to compute the polarization properties of radiation emerging from the BL (which we identify as the region responsible for the Comptonization spectrum of NS-LMXBs). The difference from an energy-dependent treatment of the electron scattering cross-section is marginal in this context, however. Indeed, the electron temperature of the Comptonization medium in the NS-LMXBs soft state never exceeds kT_e ∼ 3–5 keV (e.g., Done et al. 2007), and at this temperature the Klein-Nishina corrections are marginal.

The same considerations can be of course applied to polarimetry as well. The post-scattering PD of unpolarized and polarized radiation is (Matt et al. 1996)

$$\begin{array}{cc}& {\mathrm{\Pi }}_{\mathrm{U}}=\frac{1-{\cos }^2\mathrm{\Theta }}{1+{\cos }^2\mathrm{\Theta }-2+\eta +1/\eta },\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\mathrm{\Pi }}_{\mathrm{P}}=2\frac{1-{\sin }^2\mathrm{\Theta }{\cos }^2\mathrm{\Phi }}{\eta +{\eta }^{-1}-2{\sin }^2\mathrm{\Theta }{\cos }^2\mathrm{\Phi }},\hfill \end{array}$$(33)

where η is the ratio of the photon energies after and before scattering in the electron rest frame. For kT_e ≲ 5 keV, it is easy to see that η ∼ 1, and Π_U − Π_P have a tiny dependence on energy. Considering different optical depths of the slab and two different distributions of the seed photons, we find the following main results.

• A PA almost aligned with the NS spin can be obtained only for optically thick configurations and seed photons located at the base of the medium. Nevertheless, with this physical and geometrical setup, the total PD of the BL is ≲0.5%. Up to now, measurements of absolute PA alignment have been possible only for the Z sources Cyg X-2 (Farinelli et al. 2023) and Sco X-1; while in the former case the PA resulted to be parallel with the direction of the radio jet, in the second a polarization rotation was found by La Monaca et al. (2024) with respect to previous observations performed with PolarLight in which the PA was approximatively aligned with the radio jet as well (Long et al. 2022).

• For PA in the direction of the NS spin axis, the value is slightly rotated by ∼4°–7° for two fiducial values of the BL latitude (30° and 45°, respectively) and BL rotational frequency ν^bl = 500 Hz as a result of both special and general relativistic effects. For the extreme case of a nonrotating NS/BL system, we checked a deviation due to light bending effect ≲2°. This PA misalignment is lower than the accuracy with which it has been measured in Cyg X-2 and Sco X-1, but it could represent a powerful scientific case for future higher sensitivity facilities.

The upper limit we put on the PD of the BL at high optical depths additionally shows that whenever model-dependent spectropolarimetry attributes PD values on the order of a few percent or more to the Comptonized component (see Sect. 1), the exceeding part of the polarization signal must be attributed to an extra-component, which provides a strong contribution in terms of polarimetry, but a less significant one to the observed flux. Indeed, if the energy coverage is narrow, or the data statistics is not of high quality, a simple two-component continuum model is often enough to describe the observed spectrasatisfactorily, with the contribution of other physical processes somewhat embedded by the Comptonization component. The most natural candidate for the excess in the net polarized signal is reflection from the inner accretion disk illuminated by the radiation of the BL (LS85; Coughenour et al. 2018; Iaria et al. 2020; García et al. 2022), with the possibility of achieving PD values as high as 50% (Matt 1993); albeit, contribution from an outflowing subrelativistic wind cannot be excluded (Di Marco et al. 2023; Poutanen et al. 2023). On the other hand, if disk reflection is detected but spectropolarimetry does not allow us to split contribution to total polarization over the single components unambiguously, our theoretical results provide strong theoretical support to put tight constraints on at least one of them (the BL), helping to reduce model degeneracy.

It is also important to point out that the results presented here implicitly assume albedo A = 0 at the bottom of the slab (see Fig. 1). Physically, this resembles the case of a neutral or poorly ionized material absorbing the back-scattered radiation of a surrounding ionized atmosphere in an accretion disk, with the effect being dominant for photon energies ≲10 keV (Morrison & McCammon 1983; Haardt & Maraschi 1991). For an NS photosphere illuminating from the bottom of the BL, at the characteristic temperature of the BB-like seed photon’s (∼1–2 keV) matter is ionized and partially reflective; thus, polarimetric and angular properties of reflected radiation need to be computed either numerically (C60; Basko et al. 1974; Magdziarz & Zdziarski 1995; Poutanen & Svensson 1996) or via Monte Carlo methods (White et al. 1988b). In Fig. 8, we report the fraction of flux at the top of the layer as a function of the optical depth of the slab. The behavior can be described by the functional form

$$\begin{array}{c}\hfill \epsilon (2{\tau }_0)={\mathrm{e}}^{-a{(2{\tau }_0)}^b},\end{array}$$(34)

with a = 0.72 ± 0.013 and b = 0.54 ± 0.014. For progressively higher values of 2τ₀, an increasing fraction 1 − ε(2τ₀) of radiation is back-scattered to the bottom of the layer, where it can be partially reflected.

Mathematically, the solution of Eqs. (1) and (2) for radiation reflected at the inner boundary can be obtained by substituting the source functions S_l(τ, μ) and S_r(τ, μ) with some terms that are functionals of the incident radiation, which we may express as $\delta (\tau )I_{l,r}^{\mathrm{refl}}(\mu )=\mathcal{F}[I_{l,r}^{(-)}(2{\tau }_0,\mu )]$. The presence of a partially reflecting surface at the bottom of a scattering medium increases the average number of scattering of photons before escaping, with a hardening of the Comptonization spectral slope (Titarchuk et al. 1997). Polarization properties of emerging radiation can be affected as well. However, for high optical depths, where diffusion approximation holds, photons lose information about their initial angular distribution and polarization state (ST80; ST85). This is particularly true for photons initially located at the base of the scattering atmosphere, such as Case 1 treated in the present work. We performed a numerical check by using the crude approximation of seed photons initially polarized at a constant 20% level (see Figs. 24–26 of C60 for a more precise solution). For optical depths 2τ₀ ≳ 4, the PD angular distribution resulted to have a behavior qualitatively comparable to that shown in Fig. 3, with deviations in terms of PD value on the order of 15% for 2τ₀ = 4, progressively decreasing for increasing 2τ₀. The deviations when passing to the observer frame integrating over the BL area of course resulted to be on the same order, in turn preserving the upper limit of 1% obtained for unpolarized seed photons Thus, for albedo A = 0, results are valid for any value of the optical depth of the slab and the polarimetric properties of emerging radiation are fully described by the solution of Eqs. (1) and (2).

On the other hand, when the albedo A >  0, two things can be deduced. Firstly For high optical depths (2τ₀ ≳ 4), the photon diffusion regime is applicable; the total flux emerging from the slab is higher than the case A = 0, but polarization properties and angular distribution of emerging radiation do not vary significantly, independently of the polarization status of the initial source function. Secondly, as this is the typical optical depth of the Comptonization spectra of NS-LMXBs in a soft state, we conclude that results for the case shown in Fig. 5 are reliable. For low optical depths (less than a few), polarimetric properties of radiation emerging at the top of the slab after reflection depends on the initial state of the source function (i.e., reflected intensity at τ = 0), which then must be computed carefully.

5. Conclusions

In this paper, we presented results concerning the PD and PA of radiation scattered in a BL around an NS. These quantities were obtained by first solving the coupled system of RTEs for the two polarization modes in the comoving frame for an electron-scattering atmosphere with plane-parallel geometry and Thomson scattering cross-section, and then passing to the observer frame with back-tracing light propagation in the Schwarzschild metric. For seed photons located at the base of the slab, and typical optical depths of the Comptonization spectra of NS-LMXBs in a soft state, the PD is less than 0.5%, while the PA is in the same direction as the NS spin axis, but it is slightly rotated by a few degrees; the amount depends on the latitudinal extent of the BL as well as its rotational frequency. For lower optical depths or seed photons uniformly distributed over the slab, the PA is instead at about a right angle with the same amount of rotation.

The significantly higher values of the PD for the Comptonized component (BL) obtained with IXPE spectropolarimetric analysys of NS-LMXBs in the soft state so far are indicative that another process, presumably inner disk reflection, is responsible for it. To our knowledge, this is the first time that an upper limit on the BL polarization is provided with a thorough numerical approach, and the results presented here provide strong theoretical support for the interpretation of IXPE observations of NS-LMXB in the soft state, where the contribution of the single components to the total polarization signal cannot be completely split when the model-dependent spectropolarimetric analysis is performed.

Acknowledgments

This work is carried out as part of the research funded by INAF Grant 2023 titled “Spin and Geometry in accreting X-ray binaries: The first multi frequency spectro-polarimetric campaign”.

References

All Figures

	Fig. 1. Sketch of radiative transfer problem in plane-parallel geometry. Seed photons at the base of the slab are shown here. At any location τ, a positive angle μ relative to the slab normal can be defined for an upward and downward intensity component, respectively.
In the text

Fig. 2. Proposed geometry for soft state of NS-LMXBs; the disk extends very close to the NS surface, and then a BL forms with ΔR ≪ Rns. The BB photons from the NS photosphere are located at the base of the BL, which extends up to some latitude θmax computed from the orbital plane, and has a given temperature kTe and radial optical depth 2τ0. Only photons in the northern hemisphere (θ <  π/2) reach the observer. The total signal is obtained by integrating the contribution from each elementary area dΣ′ shown by the blue box over the BL surface.

In the text

	Fig. 3. Angular distribution of total PD (upper panel) and normalized specific intensity (lower panel) emerging from the upper layer of the plane-parallel atmosphere and different values of the optical depth. Seed photons are located at the base of the slab (Case 1). For positive values of the PD, the PA is coplanar to the slab; for negative values it is parallel to the normal.
In the text

	Fig. 4. Same as Fig. 3, but for seed photons uniformly distributed over τ (Case 2).
In the text

	Fig. 5. Total PA (upper panel), PD (middle panel), and normalized intensity (lower panel) of the BL as a function of the viewing angle with seed photons at the base of the slab (Case 1) for an optical depth of 2τ0 = 5, spin frequency of νbl = 500 Hz, and two values of the BL latitude. The PA is measured with respect to the projection on the sky of the NS spin axis.
In the text

	Fig. 6. Same as in Fig. 5, but for different values of the optical depth and fixed BL latitude θmax = 45° and spin νbl = 500 Hz.
In the text

	Fig. 7. Same as Fig. 5, but with seed photons uniformly distributed (Case 2) and different values of the optical depth of the BL, with a latitude of θmax = 45° and spin of νbl = 500 Hz.
In the text

	Fig. 8. Fraction of flux emerging from the top of the slab in the comoving frame as a function of the optical depth (see Eqs. (25)–(27)) when seed photons are located at the base.
In the text