Issue 
A&A
Volume 684, April 2024



Article Number  A62  
Number of page(s)  9  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202348915  
Published online  04 April 2024 
The polarization of the boundary layer around weakly magnetized neutron stars in Xray binaries
^{1}
INAF – Osservatorio di Astrofisica e Scienza dello Spazio, Via P. Gobetti 101, 40129 Bologna, Italy
email: ruben.farinelli@inaf.it
^{2}
Institut für Astronomie und Astrophysik, Kepler Center for Astro and Particle Physics, University of Tübingen, Sand 1, 72076 Tübingen, Germany
^{3}
Department of Astronomy, University of Geneva, Chemin d’Ecogia 16, 1290 Versoix, Switzerland
Received:
12
December 2023
Accepted:
13
January 2024
Context. Xray binaries hosting a compact object have been among the main targets of the Imaging Xray Polarimetry Explorer (IXPE) since its launch, due to their high brightness in the 2–8 keV energy band. The spectropolarimetric analysis performed so far has proved to be of great importance in providing constraints on the accretion geometry of these systems. However, the data statistics is not enough to unambiguously disentangle the contribution of the single components to the net observed polarimetric signal.
Aims. In this work, we aim to present a model for computing the polarization degree and polarization angle of the boundary layer around weakly magnetized neutron stars in lowmass Xray binaries in the soft state. The main motivation is to provide strong theoretical support to data interpretation of observations performed by IXPE or future satellites for Xray polarimetry.
Methods. The results were obtained by modeling the boundary layer as an equatorial belt around the compact object and locally approximating it as a planeparallel scattering atmosphere, for which the associated radiative transfer equation for polarized radiation in the Thomson limit was solved. The polarimetric quantities were then transformed from the comoving frame to the observer frame using the numerical methods formerly developed for Xray pulsars.
Results. For typical values of the optical depth and electron temperature of the boundary layer of these systems in a soft state, the polarization degree was less then 0.5%, while the polarization angle was rotated by ≲5° with respect to the neutron star spin axis due to special and general relativistic effects for fast rotation, the amount progressively decreasing for lower spin frequencies. The derived quantities can be used to remove degeneracy when multicomponent spectropolarimetry is performed.
Key words: polarization / radiative transfer / relativistic processes / scattering
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Lowmass Xray binaries hosting a weakly magnetized neutron star (NSLMXBs) are among the brightest sources of the Xray sky and have been studied since the dawning of Xray astronomy. In the Milky Way, about 150 NSLMXBs 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 Xray spectrum is usually modeled in terms of a soft thermal component, plus an unsaturated Comptonization spectrum of ∼1 keV blackbodylike 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 Xray 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^{−3} up to ∼L_{Edd}), sharing in this fashion a common property with blackhole Xray binaries (Done et al. 2007). The changes in luminosity are linked to variations of the spectral state, with the brightest luminosity achieved when the Xray spectrum is Zlike (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 Xray 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 NSLMXBs 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 Xray spectrum could be equally well described by an unsaturated Comptonization component plus either a direct blackbody (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 broadband 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 frequencyresolved spectroscopy of archival data of the Rossi XRay 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 beltlike fashion, with a power spectrum resembling the Fourier frequencyresolved spectrum at the period of quasiperiodic oscillations.
Significant expectations for a major leap in knowledge were put on the Imaging Xray 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 NSLMXBs, in particular in the region where the soft and hard components overlap. Indeed, Xray 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 Xray facilities covering a wider energy range have shown how spectropolarimetric information can provide unprecedented constraints on the inner region of accreting Xray binaries (see e.g., Krawczynski et al. 2022, for the black hole Cyg X–1). A sample of NSLMXBs 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 modeldependent spectropolarimetric analysis. Concerning Z sources, for Cyg X2 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 timeintegrated 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 X1, 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 IXPEalone Xray 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 Xray 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 Xray 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 KleinNishina (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, Xray 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 Xray bursters modeled with an accretion disk and a BL approximately at a right angle, both during the burst and the outofburst phases. Despite being a work dating back to the mid1980s, 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 heightscale 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 Xray 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 phaseindependent problem in which the BL covers the NS surface in a beltlike manner. We adopted a twophase 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 planeparallel 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 planeparallel geometry
We considered planeparallel 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 electronscattering atmosphere in the Thomson limit are (Pomraning 1973)
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, .
The Chandrasekhar scattering matrix 𝒫_{i, j}(μ, μ′) is defined as (ST85) follows:
The system of two integrodifferential 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 boundbound or boundfree 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
for which the general solution is
We define the bottom of the slab at τ = 0 and the top at τ = 2τ_{0}. We note that we conventionally define the total optical depth of the slab as 2τ_{0} 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_{1} = 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_{0}δ(τ). We note that with this configuration S(τ)=0 for μ < 0, and in this case Eq. (9) becomes
For the second case, labeled Case 2, we consider isotropic seed photons uniformly distributed over the slab, namely S(τ)=S_{0}, and in this case the solution of the RTE gives
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 diffused into the μ direction from all other directions. This contribution arises from both photons, which remain in the same polarization mode (e.g., ), as well as those switching from one mode to the other (e.g., ).
The system of Eqs. (1) and (2) then becomes
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 planeparallel 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 lmode, the integral terms can thus be written as
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 integrodifferential equations according to
The first two equations of the system solve for the upstream intensities and with τ computed from bottom to top, while the other two solve for the downstream components and with τ computed from top to bottom. Once defined the seed photons distribution S(τ), the algorithm first seeks for the solution as defined in Eqs. (10)–(11) or (12)–(13), and then iteratively integrates over τ equations up to a userdefined 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 GaussLegendre (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
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 RungeKutta PrinceDormand method with adaptive stepsize control provided by the Gnu Scientific Library^{1}.
By definition, the iterative method provides the solution for each scattering order . We are interested in the following quantities at the top of the slab (τ = 2τ_{0}), where it is assumed that an observer is located: angular distribution of the total specific intensity (limbdarkening 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
and the total polarization comes directly to
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,
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
where F^{( + )} and F^{( − )} are the upstream (τ = 2τ_{0}) and downstream fluxes (τ = 0), respectively.
The parameter ϵ(2τ_{0}) 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
Fig. 2. Proposed geometry for soft state of NSLMXBs; the disk extends very close to the NS surface, and then a BL forms with ΔR ≪ R_{ns}. 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 kT_{e} 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
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 corotating 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
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 twodimensional 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 planeparallel geometry. Numerical results were obtained with a C code optimized with a thirdparty adaptive integrator^{2}.
3. Results
We first discuss the results obtained from the solution of the RTE in planeparallel 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τ_{0}.
We focus on the results at the top of the slab, where we assume that the comoving observer is located. A configuration with a planeparallel electronscattering 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 freefree 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τ_{0} ≲ 1, the PD is always negative, which implies a PA aligned with the slab normal. For the intermediate value 2τ_{0} = 2, PA swapping occurs around μ ∼ 0.3, while for 2τ_{0} = 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τ_{0} = 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 firstglance impression that at some more edgeon viewing angles the flux is higher than the faceon 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 planeparallel 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. 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 NSLMXBs 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τ_{0} = 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 edgeon to a faceon view of the system, while for θ^{max} = 30° it achieves its maximum value at i ∼ 50°. The BL properties for lower values of 2τ_{0} 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τ_{0} = 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τ_{0} = 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 lightbending 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)
where E_{0} is the initial energy. This is true as long as E_{0} ≲ E_{k} ≪ 4kT_{e}; otherwise, the recoil effect is to be taken into account. For the typical electron temperature of NSLMXBs 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 appears too rough an approximation. We note that ST85 presented results both for the energyintegrated spectrum, as in Eq. (22), and for what they call qualitatively “hard radiation spectrum”, namely for with , 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 powerlawlike 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_{0} and I(E)∝e^{−α} (up to E ∼ kT_{e}) for E > E_{0}, 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 kiloelectronvolts, the powerlaw 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 NSLMXBs do not allow thorough investigations of detailed energydependent 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 modelindependent 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 modeldependent 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 energyintegrated theoretical results (i.e., summed over all scattering orders k) is sufficient to provide theoretical support to IXPE modeldependent 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 NSLMXBs). The difference from an energydependent treatment of the electron scattering crosssection is marginal in this context, however. Indeed, the electron temperature of the Comptonization medium in the NSLMXBs soft state never exceeds kT_{e} ∼ 3–5 keV (e.g., Done et al. 2007), and at this temperature the KleinNishina corrections are marginal.
The same considerations can be of course applied to polarimetry as well. The postscattering PD of unpolarized and polarized radiation is (Matt et al. 1996)
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 X2 (Farinelli et al. 2023) and Sco X1; 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 X2 and Sco X1, 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 modeldependent 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 extracomponent, 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 twocomponent 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 backscattered 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 BBlike 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
with a = 0.72 ± 0.013 and b = 0.54 ± 0.014. For progressively higher values of 2τ_{0}, an increasing fraction 1 − ε(2τ_{0}) of radiation is backscattered 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 . 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τ_{0} ≳ 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τ_{0} = 4, progressively decreasing for increasing 2τ_{0}. 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τ_{0} ≳ 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 NSLMXBs 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 electronscattering atmosphere with planeparallel geometry and Thomson scattering crosssection, and then passing to the observer frame with backtracing 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 NSLMXBs 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 NSLMXBs 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 NSLMXB in the soft state, where the contribution of the single components to the total polarization signal cannot be completely split when the modeldependent 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 Xray binaries: The first multi frequency spectropolarimetric campaign”.
References
 Abarr, Q., & Krawczynski, H. 2020, ApJ, 889, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Avakyan, A., Neumann, M., Zainab, A., et al. 2023, A&A, 675, A199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Basko, M. M., Sunyaev, R. A., & Titarchuk, L. G. 1974, A&A, 31, 249 [NASA ADS] [Google Scholar]
 Boella, G., Butler, R. C., Perola, G. C., et al. 1997, A&AS, 122, 299 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Capitanio, F., Fabiani, S., Gnarini, A., et al. 2023, ApJ, 943, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1960, Radiative Transfer (New York: Dover) [Google Scholar]
 Cocchi, M., Farinelli, R., & Paizis, A. 2011, A&A, 529, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cocchi, M., Gnarini, A., Fabiani, S., et al. 2023, A&A, 674, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Coughenour, B. M., Cackett, E. M., Miller, J. M., & Ludlam, R. M. 2018, ApJ, 867, 64 [CrossRef] [Google Scholar]
 D’Amico, F., Heindl, W. A., Rothschild, R. E., & Gruber, D. E. 2001, ApJ, 547, L147 [CrossRef] [Google Scholar]
 Di Marco, A., La Monaca, F., Poutanen, J., et al. 2023, ApJ, 953, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Di Salvo, T., Stella, L., Robba, N. R., et al. 2000, ApJ, 544, L119 [CrossRef] [Google Scholar]
 Di Salvo, T., Robba, N. R., Iaria, R., et al. 2001, ApJ, 554, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Done, C., Gierliński, M., & Kubota, A. 2007, A&ARv, 15, 1 [Google Scholar]
 Dovčiak, M., Muleri, F., Goosmann, R. W., Karas, V., & Matt, G. 2011, ApJ, 731, 75 [CrossRef] [Google Scholar]
 Fabiani, S., Capitanio, F., Iaria, R., et al. 2024, A&A, in press, https://doi.org/10.1051/00046361/202347374 [Google Scholar]
 Falanga, M., Götz, D., Goldoni, P., et al. 2006, A&A, 458, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Farinelli, R., Titarchuk, L., Paizis, A., & Frontera, F. 2008, ApJ, 680, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Farinelli, R., Paizis, A., Landi, R., & Titarchuk, L. 2009, A&A, 498, 509 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Farinelli, R., Fabiani, S., Poutanen, J., et al. 2023, MNRAS, 519, 3681 [NASA ADS] [CrossRef] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, 3rd edn. (Cambridge: Cambridge University Press) [Google Scholar]
 García, J. A., Dauser, T., Ludlam, R., et al. 2022, ApJ, 926, 13 [CrossRef] [Google Scholar]
 Gierliński, M., & Done, C. 2002, MNRAS, 337, 1373 [CrossRef] [Google Scholar]
 Gilfanov, M., Revnivtsev, M., & Molkov, S. 2003, A&A, 410, 217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haardt, F., & Maraschi, L. 1991, ApJ, 380, L51 [Google Scholar]
 Hasinger, G., & van der Klis, M. 1989, A&A, 225, 79 [NASA ADS] [Google Scholar]
 Homan, J., Steiner, J. F., Lin, D., et al. 2018, ApJ, 853, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Iaria, R., Mazzola, S. M., Di Salvo, T., et al. 2020, A&A, 635, A209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Inogamov, N. A., & Sunyaev, R. A. 1999, Astron. Lett., 25, 269 [NASA ADS] [Google Scholar]
 Krawczynski, H., Muleri, F., Dovčiak, M., et al. 2022, Science, 378, 650 [NASA ADS] [CrossRef] [Google Scholar]
 La Monaca, F., Di Marco, A., Poutanen, J., et al. 2024, ApJ, 960, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Lapidus, I. I., & Sunyaev, R. A. 1985, MNRAS, 217, 291 [NASA ADS] [CrossRef] [Google Scholar]
 Lavagetto, G., Iaria, R., di Salvo, T., et al. 2004, Nucl. Phys. B Proc. Suppl., 132, 616 [Google Scholar]
 Loktev, V., Veledina, A., & Poutanen, J. 2022, A&A, 660, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Long, X., Feng, H., Li, H., et al. 2022, ApJ, 924, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Maccarone, T. J., & Coppi, P. S. 2003, A&A, 399, 1151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Magdziarz, P., & Zdziarski, A. A. 1995, MNRAS, 273, 837 [Google Scholar]
 Mainardi, L. I., Paizis, A., Farinelli, R., et al. 2010, A&A, 512, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matt, G. 1993, MNRAS, 260, 663 [NASA ADS] [CrossRef] [Google Scholar]
 Matt, G., Feroci, M., Rapisarda, M., & Costa, E. 1996, Radiat. Phys. Chem., 48, 403 [Google Scholar]
 Mitsuda, K., Inoue, H., Nakamura, N., & Tanaka, Y. 1989, PASJ, 41, 97 [NASA ADS] [Google Scholar]
 Morrison, R., & McCammon, D. 1983, ApJ, 270, 119 [Google Scholar]
 Narayan, R., & Yi, I. 1995, ApJ, 452, 710 [NASA ADS] [CrossRef] [Google Scholar]
 Paizis, A., Farinelli, R., Titarchuk, L., et al. 2006, A&A, 459, 187 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pihajoki, P., Mannerkoski, M., Nättilä, J., & Johansson, P. H. 2018, ApJ, 863, 8 [Google Scholar]
 Podgorný, J., Dovčiak, M., Marin, F., Goosmann, R., & Różańska, A. 2022, MNRAS, 510, 4723 [CrossRef] [Google Scholar]
 Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (New York: Pergamon Press) [Google Scholar]
 Poutanen, J. 2020, A&A, 641, A166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [CrossRef] [Google Scholar]
 Poutanen, J., Veledina, A., & Beloborodov, A. M. 2023, ApJ, 949, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Revnivtsev, M. G., & Gilfanov, M. R. 2006, A&A, 453, 253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Revnivtsev, M. G., Suleimanov, V. F., & Poutanen, J. 2013, MNRAS, 434, 2355 [NASA ADS] [CrossRef] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: Wiley) [Google Scholar]
 Schnittman, J. D., & Krolik, J. H. 2009, ApJ, 701, 1175 [CrossRef] [Google Scholar]
 Schnittman, J. D., & Krolik, J. H. 2013, ApJ, 777, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Suleimanov, V., & Poutanen, J. 2006, MNRAS, 369, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Sunyaev, R. A., & Titarchuk, L. G. 1980, A&A, 86, 121 [NASA ADS] [Google Scholar]
 Sunyaev, R. A., & Titarchuk, L. G. 1985, A&A, 143, 374 [NASA ADS] [Google Scholar]
 Taverna, R., Marra, L., Bianchi, S., et al. 2021, MNRAS, 501, 3393 [NASA ADS] [Google Scholar]
 Titarchuk, L., Mastichiadis, A., & Kylafis, N. D. 1997, ApJ, 487, 834 [NASA ADS] [CrossRef] [Google Scholar]
 Ursini, F., Farinelli, R., Gnarini, A., et al. 2023, A&A, 676, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weisskopf, M. C., Soffitta, P., Baldini, L., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 026002 [NASA ADS] [CrossRef] [Google Scholar]
 White, N. E., & Holt, S. S. 1982, ApJ, 257, 318 [NASA ADS] [CrossRef] [Google Scholar]
 White, N. E., Stella, L., & Parmar, A. N. 1988a, ApJ, 324, 363 [NASA ADS] [CrossRef] [Google Scholar]
 White, T. R., Lightman, A. P., & Zdziarski, A. A. 1988b, ApJ, 331, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, S., Santangelo, A., Feroci, M., et al. 2019, Sci. China Phys. Mech. Astron., 62, 29502 [Google Scholar]
All Figures
Fig. 1. Sketch of radiative transfer problem in planeparallel 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 NSLMXBs; the disk extends very close to the NS surface, and then a BL forms with ΔR ≪ R_{ns}. 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 kT_{e} 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 planeparallel 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.