© ESO 2021

1. Introduction

XB 1916-053 is a low-mass X-ray binary system (LMXB) showing dips and type I X-ray bursts. The source was the first LMXB in which periodic absorption dips were detected (Walter et al. 1982; White & Swank 1982), the analysis of which allowed an orbital period of 50 min and its compact nature to be estimated.

The optical counterpart of the source (a V = 21 star) was discovered by Grindlay et al. (1988). Swank et al. (1984) discussed that the companion star (CS) is not hydrogen exhausted by analyzing the thermonuclear flash models of X-ray bursts, while Paczynski & Sienkiewicz (1981) showed that X-ray binary systems with orbital periods shorter than 81 min cannot contain hydrogen-rich CS. Galloway et al. (2008), by studying two consecutive type I X-ray bursts temporally separated by 6.3 h, estimated that the CS contains a 20% fraction of hydrogen under the hypothesis of a NS mass of 1.4 M_⊙. The same authors estimated a distance to the source of 8.9 kpc.

The optical light curve of the source shows a periodicity of 3027.4 ± 0.4 (Grindlay et al. 1988). Callanan et al. (1995) found that the optical period remained stable over a baseline of seven years, and refined the period estimate to 3027.551 ± 0.004 s. The most recent X-ray orbital ephemeris of the source gives an orbital period of 3000.6496 ± 0.0008 s and an orbital period derivative of 1.44(6) × 10⁻¹¹ s s⁻¹ (Iaria et al. 2015). The 1% discrepancy between the optical and X-ray periods was explained by Grindlay et al. (1988), who invoked the presence of a third body with a period of 2.5 d and a retrograde orbit. This scenario predicts that the optical period reflect the real period of the binary system, while the observed X-ray period is modified by an increase in mass transfer influenced by the presence of the third body.

An alternative scenario suggests that the real period of XB 1916-053 is actually the observed X-ray period. The SU Ursae Majoris (SU UMa) superhump scenario was initially adopted by White (1989). The superhump phenomenon is the appearance of a periodic or quasi-periodic modulation in the light curves of the SU UMa class of dwarf novae while showing superoutburst activity. The period of the superhumps is longer than the orbital period of the system in which it is observed. The commonly accepted interpretation of the phenomenon ties the longer modulations to the beat between the binary period and the apsidal precession period of the elliptical accretion disk.

Hydrodynamic simulations show that accretion disks in cataclysmic variable (CV) systems with a low-mass secondary are tidally unstable with a high probability of forming a precessing eccentric disk (Whitehurst 1988). Hirose & Osaki (1990) found the dependence of the dynamical term from the apsidal precession period of the disk. Chou et al. (2001) analyzed optical and X-ray data estimating the superhump period and the associated apsidal precession period of 3.9087 ± 0.0008 days. Adopting the SU SMa scenario and the relation discussed by Hirose & Osaki (1990), Chou et al. (2001) inferred a mass ratio q = M₂/M₁ of 0.022, where M₁ and M₂ are the NS and CS mass, respectively.

Retter et al. (2002) detected a further X-ray periodicity at 2979 s in the RXTE light curves of XB 1916-053. A similar periodicity had already been observed by Smale et al. (1989) analyzing Ginga data. Retter et al. (2002) interpreted the 2979 s period invoking the presence of a negative superhump (also called an infrahump) where the X-ray period is the orbital period of the system that beats with the nodal precession period of 4.86 days. This means that the disk is tilted with respect to the equatorial plane. Finally, Hu et al. (2008), by adopting the relation between the nodal angular frequency and the mass ratio q proposed by Larwood et al. (1996) and assuming a nodal precession period of 4.86 days, inferred q ≃ 0.045.

Iaria et al. (2015) used a baseline of 37 years to infer the accurate orbital ephemeris of the source from the dip arrival times. Although they proposed several models to describes the dip arrival times, the authors focused their attention on the ephemeris containing a quadratic term plus a long sinusoidal modulation. The quadratic term implies the presence of an orbital period derivative of 1.44(6) × 10⁻¹¹ s s⁻¹, while the long modulation of 25 years was explained with the presence of third body with mass of 0.055 M_⊙ orbiting around the binary system. Iaria et al. (2015) showed that the large orbital period derivative can be explained only assuming a high non-conservative mass transfer rate with only 8% of the mass transfer rate accreting onto the NS.

The spectrum of XB 1916-053 shows absorption lines superimposed on the continuum emission. The Fe XXV and Fe XXVI absorption lines were detected for the first time by Boirin et al. (2004) from the analysis of XMM-Newton data. Iaria et al. (2006), using Chandra data, detected prominent absorption lines in the spectrum associated with the presence of Ne X, Mg XII, Si XIV, SXVI, and Fe XXVI ions. From the plasma diagnostic the authors proposed that the lines should be produced at the outer rim of the accretion disk. Analyzing the same data, Juett & Chakrabarty (2006) suggested a thickness for the X-ray absorber of < 3.2 × 10⁹ cm, assuming that it is located at the outer edge of the accretion disk. Finally, analyzing Suzaku data, Gambino et al. (2019) set an upper limit of < 1 × 10¹⁰ cm on the distance of the ionized absorber from the NS, placing it at the innermost region of the accretion disk.

In this work we update the orbital ephemeris of XB 1916-053 by expanding the available baseline to 40 years using ten Chandra observations and one Swift/XRT observation. We revise the estimation of the NS mass suggesting that a mass of 1.4 M_⊙ can also explain the observed orbital period derivative and luminosities. From theoretical considerations we discuss that a mass ratio of 0.048 can explain both the apsidal and the nodal precession period of the accretion disk taking into account the presence of a spiral wave in the accretion disk with a pitch angle lower than 30°. From the spectroscopic analysis we study the absorption lines detected in the spectrum during the persistent emission at different source luminosities and during the dips.

We note that the same Chandra data have recently been analyzed by Trueba et al. (2020). The authors observed a redshift in the absorption lines that can be interpreted as a gravitational redshift. We reached similar results by assuming the same scenario, but adopting a different analysis.

2. Observations

We analyzed ten new observations of XB 1916-053 collected by the Chandra observatory in 2018 and the Swift/XRT observation (ObsID. 87248021) taken on 2017 July 30; furthermore, we reanalyzed the Chandra observation taken in 2004 August 7 (see Iaria et al. 2006).

The new Chandra observations were performed from 2018 June 11 to August 6 using the onboard High Energy Transmission Grating Spectrometer (HETGS) in timed graded mode; the detailed list of the observations is shown in Table 1. HETGS consists of two types of transmission gratings, the medium energy grating (MEG) and the high-energy grating (HEG). The HETGS provides high-resolution spectroscopy from 1.2 to 31 Å (0.4–10 keV) with a peak spectral resolution of λ/Δλ ∼ 1000 at 12 Å for first-order HEG. The dispersed spectra were recorded with an array of six charge-coupled devices (CCDs) that are part of the Advanced CCD Imaging Spectrometer-S.

We reprocessed the data using the FTOOLS ver. 6.26.1 and the CIAO ver. 4.11 packages. The brightness of the source required additional efforts to mitigate pileup issues. A 512 row subarray (the first row = 1) was applied during the observations, reducing the CCD frame time to 1.7 s. We ignored the zeroth-order events in our analysis focusing on the first-order HEG and MEG spectra. We reprocessed the data adopting a scale factor (width_factor_hetg) that multiplies the approximate one sigma width of the HEG/MEG mask in the cross-dispersion direction of 15, instead of the default value of 35, in order to avoid that HEG and MEG overlap in the Fe K region of the spectrum. We obtained HEG and MEG region widths of 33 and 44 pixels, respectively.

We applied barycentric correction to the events adopting the source coordinated obtained by Iaria et al. (2006) and extracted the light curve for each observation retaining only the first-order MEG and HEG data in the 0.3–10 keV energy range. We show the light curves of the ten observations in Fig. 1. During the persistent emission, the observations taken in 2018 June (ObsID. 20171, 21103, 21104, 21105) have a higher count rate than those taken in July-August (ObsID. 20172, 21662, 21663, 21664, 21106, 21666). Observations 21103, 21104, and 21105 show a count rate close to 7 c s⁻¹ (20171 shows a count rate of 5 c s⁻¹); observations 20172, 21662, 21663, and 21664 show a count rate of 4.5 c s⁻¹; and observations 21106 and 21666 show a count rate of 3 c s⁻¹. We observed five type I X-ray bursts, four during the high-flux state and one when the flux of the source dropped below 5 c s⁻¹. Furthermore, the light curves corresponding to observations 21104 and 21105 (when XB 1916-056 showed a high flux) do not show dips, while the dipping activity is intense at lower fluxes. A similar behavior was observed during a long Suzaku observation of the source (see Gambino et al. 2019).

Fig. 1. First-order MEG+HEG light curves of the observations in the 0.3–10 keV energy range. The bin time is 128 s.

XB 1916-053 was monitored by the X-Ray Telescope (XRT; Burrows et al. 2005) on board the Swift Observatory (Gehrels et al. 2004). We analyzed the Swift/XRT observation 87248021 taken on 2017 July 30 10:54:17, with elapsed and exposure times of 87 ks and 8.9 ks, respectively. The observation was performed in photon counting (PC) mode. The data were locally reprocessed by the UK Swift Science Data Center using HEASOFT v6.26.1. We applied barycentric correction to the events using the ftool barycorr¹. Using XSELECT we extracted the source events in the 1–4 keV energy range adopting a circular region of 100″; the background events in the same energy range were extracted using a same size circular region free from the source. Finally, we applied the exposure-correction to the source and background-light curves and extracted the background-subtracted light curve using the ftools xrtlccorr and lcmath².

3. Updated orbital ephemeris of XB 1916-053

First, we filtered the light curves to exclude the type I X-ray bursts; then we selected the events in the 1–4 keV energy range to estimate the dip arrival times. We estimated one dip arrival time from the Chandra observations taken in June, one from those taken in July-August, and one from the XRT observation.

The light curves were folded adopting P₀ = 3000.6511 s as the orbital period. For the first dip arrival time we folded the light curve obtained from observation 20171, where the dips are the deepest, adopting a folding time T_fold = 58281 MJD; to obtain the second we folded together the light curves of observations 20172, 21664, 21106, and 21666, adopting T_fold = 58333.65 MJD. Finally, the XRT light curve was folded assuming T_fold = 57965 MJD.

We fitted the dips with a simple model consisting of a step-and-ramp function, where the count rates before, during, and after the dip are constant and the intensity changes linearly during the dip transitions (see Iaria et al. 2015, for a description). We show the X-ray dip arrival times in Table 2.

To obtain the delays with respect to a constant period reference, we used the period P₀ = 3000.6511 s and reference epoch T₀ = 50123.00873 MJD. Hence, we added the three dip arrival times to the 27 values reported by Iaria et al. (2015).

We initially fitted the delays with a quadratic function y(t) = a + bt + ct², where t is the time in days with respect to T₀, a = ΔT₀ is the correction to T₀ in units of seconds, b = ΔP/P₀ in units of s d⁻¹ with ΔP the correction to the orbital period, and c = 1/2 Ṗ/P₀ in units of s d⁻² with Ṗ is the orbital period derivative. The quadratic form does not fit the data well, returning a χ²(d.o.f.) of 182(27). We show the delays (red points) and the quadratic curve (black curve) in the top panel of Fig. 2. The best-fit values of the parameters are shown in the second column of Table 3.

Fig. 2. Top panel: delays vs time (red points). The black curve is the best fit of the LQ function. The brown curve is the best fit of the LQS function. Bottom panel: residuals in units of σ vs time, corresponding to the LQS function

To improve the modeling of the delays we added a cubic term to the previous parabolic function, y(t) = a + bt + ct² + dt³, where d is defined as and indicates the temporal second derivative of the orbital period. Fitting with this cubic function, we obtained a χ²(d.o.f.) of 129(26) with an F-test probability of chance improvement of only 3 × 10⁻³ with respect to the quadratic form (hereafter LQ function). We also tried to fit the delays using a linear plus a sinusoidal function (the LS function in Iaria et al. 2015), where a and b are defined as above, while A and P_mod are the amplitude in seconds and the period in days of the sinusoidal function, respectively. Finally, t_ϕ is the time in days referred to T₀ at which the sinusoidal function is null. Since the best-fit value of P_mod is larger then the available time span, we fixed it to the best-fit value P_mod = 73574 d. We found a χ²(d.o.f.) of 132(26); the best-fit parameters are shown in the third column of Table 3.

We added a quadratic term to the LS function to take into account the possible presence of an orbital period derivative (hereafter LQS function). We obtained a value of χ²(d.o.f.) of 37(24) and an F-test probability of chance improvement with respect to the LQ function of 1.8 × 10⁻⁸. We show the best-fit curve of the LQS model (brown curve) and its relative residuals in units of sigma in the top and bottom panels of Fig. 2, respectively. The best-fit values are shown in the forth column of Table 3. The obtained best-fit values are compatible with the previous ones reported by Iaria et al. (2015).

The corresponding LQS ephemeris is

(1)

with P_mod = 9099 ± 302 d (24.9 ± 0.8 yr), t_ϕ = 1262 ± 133 d, and A = 130 ± 14 s. The corresponding orbital period derivative is Ṗ = 1.46(3) × 10⁻¹¹ s s⁻¹. Our analysis suggests that a quadratic or a quadratic plus a cubic term do not adequately fit the delays, as already shown by Iaria et al. (2015).

We folded the 0.3–10 keV first-order MEG+HEG light curves obtained from the observations (type I X-ray bursts excluded), adopting the orbital period inferred by the LQS ephemeris and changing the reference time so that the dip falls at the orbital phase 0.5. The epoch-folded light curves are shown in Fig. 3.

Fig. 3. First-order MEG+HEG folded light curves in the 0.3–10 keV energy range adopting the LQS ephemeris. An ad hoc phase shift of 0.5 is imposed. The label in the plot shows the correspondence between colors and observations. The blue dashed horizontal lines indicate the arbitrary count rate threshold of 3.5 c s−1 and 5.5 c s−1 adopted to select the four sets of observations (see text).

We distinguish four different sets of epoch-folded light curves at different count rates: the first set (hereafter set A) associated with observation 4584 has a persistent count rate higher than 10 c s⁻¹ (gold curve in Fig. 3); the folded curves of the second set (hereafter set B), have a persistent count rate between 6 and 8 c s⁻¹ and are associated with observations 21103, 21104, and 21105 (orange, blue, and magenta curve in Fig. 3); the third set of curves (hereafter set C) has a persistent count rate between 3.5 and 5.5 c s⁻¹ and is associated with observations 20171, 20172, 21662, 21663, and 21664 (black, red, brown, gray, and purple curve in Fig. 3); the fourth set (hereafter set D) shows a count rate close to 3.5 c s⁻¹ outside the dip associated with observations 21106 and 21666 (yellow and green curve in Fig. 3).

4. Spectral analysis

We extracted the combined spectrum from each distinguished set, calling the spectra A, B, C, and D, using the CIAO script combine_grating_spectra and combining the corresponding first-order MEG and HEG spectra. Since our aim is the spectral analysis out of the dips, we excluded the orbital phases in which the dip falls when it is present. The excluded phases are shown in Table 4. Spectrum B was obtained by excluding the phases between 0.45 and 0.6 from observation 21103; the other two observations do not show dipping activity, and we included all the phases. We rebinned the resulting first-order MEG and HEG spectra of selections A, B, C, and D to have at least 250 counts per energy channel. Such a large rebinning could affect the spectroscopic study of narrow discrete features because it risks losing energy resolution. However, we tried to fit these spectra using a rebinning so as to have 25 counts per energy channel, verifying that we obtained results that are compatible with those reported in the following. For the analysis we adopted the energy ranges 0.7–7 keV and 0.9–10 keV for the first-order MEG and HEG spectrum, respectively.

To fit the spectra we used XSPEC v12.10.1p; we adopted the cosmic abundances and the cross sections derived by Wilms et al. (2000) and Verner et al. (1996a), respectively. To take into account the interstellar absorption we adopted the Tübingen-Boulder model (TBABS in XSPEC). We fitted the spectrum using a Comptonized component (NTHCOMP in XSPEC; Zdziarski et al. 1996). Finally, we added a constant to take into account the different normalizations of the two instruments. The initial model, called Model 1, is defined as

(2)

We fixed at 0 the value of the parameter inp_type in NTHCOMP, assuming the seed-photon spectrum injected in the Comptonizing cloud to be a blackbody. The spectra A, B, C, and D were fitted simultaneously. The value of the seed-photon temperature kT_bb was tied at the same value in each spectrum and, finally, the value of the interstellar equivalent neutral hydrogen column density was tied at the same value in spectra C and D.

We show the unfolded spectra and the residuals in the left panels of Fig. 4. The best-fit parameters are shown in Table 5. Strong residuals compatible with absorption lines are visible at 1 keV, 1.47 keV, 2 keV, 2.62 keV, and 6.97 keV. Furthermore, below 0.9 keV there is a clear excess in the residuals.

Fig. 4. Unfolded spectrum and residuals adopting Model 1 (left panels) and Model 2 (right panels). The black, red, green, and blue data correspond to the first-order MEG of spectra A, B, C, and D, respectively. The cyan, purple, yellow, and orange data correspond to first-order HEG spectrum of spectra A, B, C, and D, respectively. The prominent absorption lines in residuals of Model 1 (bottom left panel) are associated with the presence of Ne X, Mg XII, Si XIV, S XVI, and Fe XXVI ions.

We fitted each observed absorption line with a GAUSSIAN component. The energies, widths, intensities, and equivalent widths of each identified line are shown in Table 6 for each spectrum. We detected the absorption lines associated with the presence of ions of Ne X, Mg XII, Si XIV, S XVI, CaXX, and Fe XXVI. The line associated with the Ca XX ions is detected only in spectrum A. Spectrum D only showed the line associated with Si XIV, probably because of the low flux of the spectrum. For the sake of clarity the sizes of the energy channels at which we identified the absorption lines are shown in Table 6.

To account for the presence of the prominent absorption lines, to the continuum of Model 1 we added the multiplicative component ZXIPCF³, which takes into account a partial covering of ionized absorbing material. The component reproduces the absorption from photoionized matter illuminated by a power-law source with spectral index Γ = 2.2, and it assumes that the photoionized absorber has a microturbulent velocity of 200 km s⁻¹ (see Reeves et al. 2008; Miller et al. 2007, for applications to AGN and Seyfert 1 galaxies). Recently, this component was adopted to fit the absorbing features in the Fe-K region associated with highly ionized matter surrounding the eclipsing NS-LMXB AX J1745.6-2901 (Ponti et al. 2015), the dipping source XB 1916-053 (Gambino et al. 2019), and the eclipsing source MXB 1659-298 (Iaria et al. 2019). The parameters of this component are N_Hξ, log(ξ), f, and z: N_Hξ describes the equivalent hydrogen column density associated with the ionized absorber, log(ξ) describes the ionization degree of the absorbing material, f is a dimensionless covering fraction of the emitting region, and z is the redshift associated with the absorption features. The parameter ξ is defined as ξ = L_x/(n_Hr²), where L_x is the X-ray luminosity incident on the absorbing material, r is the distance of the absorber from the X-ray source, and n_H is the hydrogen atom density of the ionized absorber.

To take into account that the ionized absorber can scatter the radiation out of the line of sight via Thomson or Compton scattering, we added the multiplicative component CABS to the model. The only parameter of the component CABS is N_Hcabs, which describes the equivalent hydrogen column density associated with the scattering cloud. Finally, to account for a partial covering of the scattering material, we multiplied CABS by the component PARTCOV. This is a convolution model that allows us to convert an absorption component into a partially covering absorption component, quantified by the parameter f_cabs. In order to make the model self-consistent, we tied f_cabs to the parameter f of the ionized absorber. Moreover, we linked the equivalent hydrogen column density parameter of the scattering cloud, N_Hcabs, to the value of the equivalent hydrogen column density of the ionized absorber, N_Hξ. While fitting simultaneously spectra A, B, C and D, we imposed f, log(ξ), and the seed-photon temperature to assume the same value for each spectrum. The model closely fits the absorption lines; however, the presence of an excess in the residuals below 0.9 persists. For this reason we added an absorption edge with the energy threshold fixed at 0.871 keV to account for the presence of O VIII ions along the line of sight, recently observed by Gambino et al. (2019) analyzing the Suzaku spectra of the source.

The adopted model, hereafter called Model 2, is defined as

The addition of ZXIPCF and the absorption edge improves the fits, with a χ²(d.o.f.) of 3944(4542) and a Δχ² of 1013 with respect to Model 1. We show the best-fit values in Table 5 for each spectrum. The unfolded spectra and the corresponding residuals are shown in the right panels of Fig. 4; the residuals associated with the absorption lines and the excess below 0.9 keV observed adopting Model 1 are now absent.

We found that the ionized absorber covers more than 86% of the emitting source and the ionization parameter log(ξ) is 4.33 ± 0.02. The equivalent hydrogen column density associated with the ionized absorber N_Hξ is , , , and (10 ± 7) × 10²² cm⁻² for spectra A, B, C, and D, respectively. Furthermore, we observed a similar redshift value for each spectrum (between 1.0 × 10⁻³ and 1.3 × 10⁻³).

The Γ parameter of NTHCOMP goes from 1.6 for spectrum A to 1.8 for spectrum C and D, indicating that the spectral shape softens going from spectrum A to spectrum D. The electron temperature of the Comptonizing cloud is compatible with 3–4 keV for all the spectra even if it is not well constrained for spectra C and D, while the seed-photon temperature is 0.16 ± 0.04 keV.

Finally, the best-fit value of the optical depth of the absorption edge is 0.14 ± 0.05 and the equivalent hydrogen column density of the interstellar matter is compatible with 0.5 × 10²² cm⁻² for each spectrum.

Assuming a distance to the source of 8.9 kpc (Galloway et al. 2008), the unabsorbed luminosity in the 0.1–100 keV energy range is 1.46 × 10³⁷ erg s⁻¹, 1.0 × 10³⁷ erg s⁻¹, 0.67 × 10³⁷ erg s⁻¹, and 0.41 × 10³⁷ erg s⁻¹ for A, B, C, and D, respectively.

5. Spectral analysis of the dips

To analyze the average spectrum during the dips we extracted the events from observations 20171, 20172, 21106, 21662, 21663, 21664, and 21666. The dip events were selected from observation 20171 taking into account the good times intervals (GTIs) obtained from the first-order MEG+HEG light curve in the 0.3–10 keV energy range where the count rate was lower than 3 c s⁻¹. Similarly, we used the threshold of 2 c s⁻¹ for ObsID. 20172, 2.4 c s⁻¹ for ObsID. 21662, 2 c s⁻¹ for ObsID. 21663 and 21664, 1.5 c s⁻¹ for ObsID. 21106, and 1.1 c s⁻¹ for ObsID. 21666. We extracted the first-order HEG and MEG spectrum from each observation and combined them using the CIAO script combine_grating_spectra. Each spectrum was then grouped to have at least 100 counts per energy channel; the resulting spectrum has an exposure time of 9.5 ks. The energy range adopted for the spectral analysis is 0.9–7 keV and 1–9 keV for MEG and HEG, respectively.

We fitted the spectrum adopting Model 2; however, we excluded the EDGE component because of the low flux of the spectrum. does not allow us to constrain it. We kept the photon-index value, the electron temperature, and the seed-photon temperature of the NTHCOMP component fixed to 1.8, 4.3 keV, and 0.16 keV, respectively. We did this so that the dip-spectrum was extracted from the observations belonging to set C and D, under the assumption that the spectral shape of the continuum emission is similar both inside and outside the dips.

We show the best-fit parameters in Table 7; the unfolded spectrum and the corresponding residuals are shown in Fig. 5.

Fig. 5. Unfolded spectrum and residuals of the dip spectrum adopting Model 2. The black and red data correspond to the first-order MEG and HEG spectra, respectively.

We found that during the dip the ionized absorber is characterized by log(ξ), suggesting a low ionization of the absorber. The equivalent hydrogen column density N_Hξ associated with the ionized absorber is cm⁻², which is larger than the values obtained during the persistent emission by a factor between 3 and 6.

The covering fraction is > 0.68, compatible with the values obtained for spectra A, B, C, and D. The equivalent hydrogen column density associated with the interstellar matter is (0.9 ± 0.2) × 10²² cm⁻². This value is larger than the observed value during the persistent emission. We cannot exclude that a partial contribution to this value comes from local neutral hydrogen in the hypothesis that during the dip the absorber is composed of neutral and mildly ionized matter.

Finally, the 0.1–100 keV unabsorbed luminosity, assuming a distance to the source of 9 kpc, is 0.57 × 10³⁷ erg s⁻¹.

6. Discussion

From the analysis of the Chandra and Swift/XRT observations we inferred three new dip arrival times that, added to the previous 27 reported by Iaria et al. (2015), allowed us to update the orbital ephemeris of XB 1916-053 by extending the temporal baseline from 36 to 40 years. The largest baseline excludes some solutions previously plausible with the available data, such as the LQ ephemeris shown by Iaria et al. (2015). We find that the LQS ephemeris better models the dip arrival times, improving the constraints on the orbital period derivative and the sinusoidal modulation. The orbital period derivative is Ṗ = 1.46(3) × 10⁻¹¹ s s⁻¹ and the sinusoidal modulation has a period of P_mod = 9099 ± 302 d.

Iaria et al. (2015) discussed that such a large orbital period derivative can be explained only assuming that the mass transfer rate from the CS to the NS is highly non-conservative. The authors estimated that more that 90% of the mass transfer rate from the CS is lost by the system, the mass ratio q = m₂/m₁ (m₁ and m₂ are the NS and CS mass in units of solar masses) is close to 0.013 and the NS mass should be larger than 2.1 M_⊙.

Here we show that the conclusions of Iaria et al. (2015) were strongly constrained by the hypothesis that the matter leaves the system from the inner Lagrangian point. To discuss this point we have to estimate the mass ratio q of XB 1916-053.

6.1. Mass ratio q of XB 1916-053

XB 1916-053 shows a superhump period P_sh (Chou et al. 2001) and a negative superhump period, also called infrahump period, P_ih (Retter et al. 2002; Hu et al. 2008). Chou et al. (2001) discussed that the optical modulation close to 3028 s is likely caused by the coupling of the orbital motion with a 3.9-day disk apsidal precession period, like the superhumps in SU UMa-type dwarf novae, obtaining from this assumption q ≃ 0.022, while Hu et al. (2008) discussed the infrahump period of 2979.3(1.1) s detected by Retter et al. (2002) as being due to the nodal precession period of the tilted accretion disk, and finding q ≃ 0.045. Considering that the orbital period of the system is close to 3000.66 s, the apsidal precession period of the disk should be 3.9087(8) d, while the retrograde nodal precession period 4.86 d. However, the reported values of q do not depict a self-consistent scenario. We show below that the value of q = 0.048 is consistent with the detected values of P_sh and P_ih only assuming an outer radius of the disk truncated at a 3:1 resonance.

Hirose & Osaki (1990), performing hydrodynamic simulations of accretion disks to study the superhump phenomenon in SU UMa stars, showed that the apsidal precession frequency ω_p can be written as ω_p/ω_orb = Z(r)q/(1 + q)^1/2, where

and

which is the Laplace coefficient of order 0 in celestial mechanics (see Brouwer & Clemence 1961, chapter 15, Eq. (42)) and r is the ratio of the accretion disk radius r_disk to the orbital separation a. Combining the expressions we find that

(3)

Hirose & Osaki (1990) demonstrated that the tidal instability in the accretion disks is caused by the resonance between the particle orbits in the disk and the companion star with a 3:1 period ratio in SU UMa dwarf novae. Since r = r_disk/a can be expressed by the equation

(4)

with j = 3 and k = 2 for a 3:1 resonance (see Frank et al. 2002, Eq. (5.125)), we can combine the two last equations to obtain

(5)

using the value ω_p/ω_orb ≃ 0.0089 inferred by Chou et al. (2001) we find q ≃ 0.0256.

However, Hirose & Osaki (1990) estimated the apsidal precession period of the disk taking into account only the dynamical effects of the matter. Lubow (1991a,b, 1992) showed that the apsidal precession rate ω for an eccentric disk is given by three terms ω = ω_dyn + ω_press + ω_tr, where ω_dyn is the dynamical precession frequency discussed by Hirose & Osaki (1990) and shown above, ω_press is a term due to the pressure, and ω_tr is a transient term related to the time-derivative of the mode giving rise to the dynamical precession that can be neglected in steady state. Since XB 1916-053 is a persistent source that does not show outburst, and since Callanan et al. (1995) verified the stability of the optical period over seven years, we do not consider the last term as we assume that the disk precesses in a steady state.

Lubow (1992) showed that the pressure term can be written as , where k is the radial wavenumber of the mode, c_s the sound speed of the gas, and Ω is the frequency of a particle in the disk at a given radius. The negative sign is due to the pressure term associated with a spiral arm in the disk, which acts in the opposite sense with respect to the dynamical term. For a spiral wave, the pitch angle θ is related to k by the relation tanθ = (kR)⁻¹, where R is the distance from the central source.

Pearson (2006), considering a 3:1 resonance and adopting Eq. (4), expressed the term associated with the pressure of the spiral wave ω_press as

(6)

where a is the orbital separation, j = 3, and η_A = 0.5[c_s/(ω_orb a)]²tan⁻²θ. We combined Eqs. (5) and (6) and obtained

(7)

where ω is now the apsidal precession frequency of the disk taking into account the pressure term.

Hu et al. (2008), assuming a 3:1 resonance, inferred a value of q ≃ 0.045 based on the negative superhump with a period of 2979.3 s observed by Retter et al. (2002), and interpreted as the beat period between the orbital period and the 4.86 d nodal precession period of the disk. The value of q ≃ 0.045 was obtained from the expression proposed by Larwood et al. (1996) and Montgomery (2009)

(8)

where ω_n is the angular frequency associated with the nodal precession of the disk, ω_k is the Keplerian frequency of the matter at a given radius r of the accretion disk that is defined in units of orbital separation a, and the angle δ describes the orbit of the CS with respect to the NS. For δ ≃ 0 the orbits of the two bodies are coplanar (Papaloizou & Terquem 1995). The minus sign is present because the nodal frequency is retrograde with respect to the orbital frequency (see Fig. 6.18 in Hellier 2001).

Rewriting the Keplerian frequency in the known form ω_k = (GM₁/R³)^1/2 and combining it with Kepler’s third law we find ω_k = ω_orb(1 + q)^−1/2r^−3/2. We combined the last expression with Eq. (8), obtaining

(9)

Assuming that r = r_disk/a (the same assumption made by Hu et al. 2008) and combining the last equation with Eq. (4) for a 3:1 resonance, we infer

(10)

Because ω_n/ω_orb ≃ 7.146 × 10⁻³, we adopt the nodal period of 4.86 days and find that q ≃ 0.048 for cos δ ≃ 1.

We show in Fig. 6 the Roche-lobe radius of the NS as a function of q (blue curve), the outer radius of the accretion disk for a 2:1 resonance (yellow curve) and for a 3:1 resonance (black curve). The green curve represents the outer radius of the disk adopting the value of ω/ω_orb ≃ 8.9 × 10⁻³ obtained by Chou et al. (2001), who neglect the pressure term due to the spiral wave in the apsidal precession rate (see Eq. (5)). The red curve represents the outer radius of the disk assuming a nodal precession of 4.86 days obtained from Eq. (10). All the radii are in units of orbital separation a.

Fig. 6. Characteristic radii of XB 1916-053 in units of orbital separation vs the mass ratio q. Show are the NS Roche-lobe radius (blue curve), the radius of the disk at which a 2:1 and 3:1 resonance could occur (yellow and black curves), the radius associated with the apsidal rate of the disk for a 3:1 resonance without and with the pressure term due to the presence of a spiral wave in the disk (green and brown curve), and the radius associated with the nodal precession period for a 3:1 resonance (red curve).

Initially, we note that the green and red curve do not have intersections that exclude any self-consistent scenario; in other words, the outer radius of the disk estimated from the observed apsidal precession period is not compatible with that obtained from the observed nodal precession period if we ignore the pressure term. Taking into account ω_press, the outer radius of the disk obtained from the apsidal precession period can intersect the red curve. It is possible to obtain a pair of r and q values for a 3:1 resonance assuming that the parameter η_A ≃ 0.0049. For that value we find q ≃ 0.048 and an outer disk radius r ≃ 0.474 in units of orbital separation. It could be possible a solution for a 2:1 resonance for which q ≃ 0.0315, r ≃ 0.62, and η_A ≃ 0.0089; however, since the truncation radius of the disk due to the tidal interaction with the CS is r_T = R_T/a = 0.6/(1 + q) (for q between 0.03 and 1; Paczynski 1977), we obtain r >  r_T, making the solution unrealistic.

A further indication that the mass ratio q could be close to 0.048 is given by the empirical expressions inferred for the CVs, in which the period excess ϵ depends on q. Using the values of the superhump period P_sh = 3027.5510(52) s and the orbital period P_orb = 3000.6508(9) (see Table 1 in Retter et al. 2002, and references therein) we can define the period excess as

(11)

Using the empirical form ϵ = 0.18q + 0.29² (Patterson et al. 2005) we obtain q ≃ 0.045. Goodchild & Ogilvie (2006) inferred ϵ = 0.2076(3)q − 4.1(6) × 10⁻⁴, where the small offset from the origin is due to pressure-induced retrograde precession of a stable eccentric mode in systems with very low q. Using the last relation we infer that q = 0.0452 ± 0.0004.

We obtain q ≃ 0.048 for

(12)

Lubow (1992) gave a range of 0.01–0.05 for the normalized speed of sound, c_s/(ω_orb a), and inferred that for a tightly wrapped spiral arm the pitch angle θ is between 5.7° and 31°. We show in Fig. 7 that for η_A ≃ 0.0049 this condition is verified, and we obtain that the pitch angle is between 5.8° and 26.7°. From simulations Montgomery (2001) restricted the range of the pitch angle between 13° and 21°; these values restrict the range of the normalized sound of speed between 0.023 and 0.038. Using Kepler’s third law and assuming q = 0.048 we obtain that the sound speed c_s is between km s⁻¹ and km s⁻¹. This range of values is compatible with a disk temperature in the outer region between 3 × 10⁴ K and 7 × 10⁴ K (see Eq. (2.21) in Frank et al. 2002), which is compatible with the results obtained by Nelemans et al. (2006) from the analysis of the optical band of XB 1916-053. The authors found that a LTE model consisting of pure helium plus overabundant nitrogen closely fits the observed spectrum finding that the model has a temperature of ∼3 × 10⁴ K.

Fig. 7. Pitch angle of the spiral wave vs the normalized sound of speed (black curve). The red horizontal lines correspond to the angles of pitch between 13° and 21°, as obtained by Montgomery (2001) from simulations.

6.2. Neutron star mass m₁

Assuming that the CS fills its Roche lobe, we can obtain an indication of the NS mass. The CS is a degenerate star and its radius R₂ is given by , where X is the fraction of hydrogen in the star. This equation has to be corrected for the thermal bloating factor b_f ≥ 1, which is the ratio of the CS radius to the radius of a star with the same mass and composition; this star is completely degenerate and supported only by the Fermi pressure of the electrons. Adopting the Roche-lobe prescription of Eggleton (1983), using Kepler’s third law and imposing that R_l2 = b_fR₂, where R_l2 is the Roche-lobe radius of the CS, we obtain

(13)

where P is the orbital period in seconds.

Heinke et al. (2013) inferred X ∼ 0.14 for a NS mass of 1.4 M_⊙, compatible with the measurements of Nelemans et al. (2006). We can predict how X varies with respect to m₁. Even though the CS is almost a pure helium dwarf, the small percentage of hydrogen that can be transferred onto the NS via Roche-lobe accretion can significantly change the energy released due to hydrogen’s larger energy release per nucleon (energy released per nucleon Q_nuc = 1.6 + 4.0X MeV nucleon⁻¹; Cumming 2003). Galloway et al. (2008), using RXTE data, measured the ratio α of burst to persistent flux between two consecutive type I X-ray bursts temporally spaced by 6.3 h obtaining α = 78.8 ± 0.3. We rewrite Eq. (6) shown by Galloway et al. (2008) as

(14)

in which we assumed a NS radius of 10 km. Adopting α = 78.8 ± 0.3 we find that X is 0.21, 0.3, and 0.39 for a NS mass of 1.4, 1.6, and 1.8 M_⊙ respectively, while for a NS mass of 2.1 M_⊙ we find that X ≃ 0.52, in contrast with the pure helium dwarf nature of the CS. Furthermore, assuming an orbital period value P = 3000.66 s and a mass ratio q = 0.048, and combining Eqs. (13) and (14), we find

(15)

We obtain that b_f is 1.97, 1.92, and 1.86 for a NS mass of 1.4, 1.6, and 1.8 M_⊙, respectively. In the following we explore the NS mass between 1.4 and 1.8 M_⊙ in order to have a fraction of hydrogen X in the CS lower than 0.4.

Finally, we can express m₂ from Eq. (13) remembering that m₁ = m₂q⁻¹. We find

(16)

6.3. Orbital period derivative and luminosity of XB 1916-053

Since XB 1916-053 has an orbital period close to 50 min we assume that the angular momentum lost via magnetic braking is negligible with respect to the angular momentum lost via gravitational waves. Combining Eqs. (1)–(3) shown by Rappaport et al. (1987) with the assumption that the stellar adiabatic index ξ_ad is −1/3, as suggested by the authors, we obtain

(17)

where P_h is the orbital period in hours, β is the fraction of the mass that transferred from the CS accretes onto the NS and α is linked to the rate of specific angular momentum lost via mass loss by the binary system. Defining , where a is the orbital separation and P the orbital period, we can define α as the ratio (d/a)², where d is the distance from the center of mass (CM) of the binary system from which the matter leaves. Finally, the function Λ(q, α, β) is defined as

(18)

Iaria et al. (2015) assumed that a large fraction of the mass transferred by the CS leaves the system at the inner Lagrangian point. Since the masses of the two bodies are unknown, we estimate the ratio b₁/a as a function of the mass ratio q between 0.005 and 0.1, where b₁ is the distance of the inner Lagrangian point from the NS. We note that Plavec & Kratochvil (1964) found an expression of b₁/a for values of q higher than 0.1 and we verified a posteriori that it does not work for the range of q of interest to us. We derived a more generic expression,

(19)

imposing that the sum of the gravitational forces of the two bodies and the centrifugal force is null at the inner Lagrangian point (L₁); the accuracy of this relation is 0.5% for q between 0.005 and 0.1. As reported above the free parameter α is defined as (d/a)², where d indicates the distance to the CM from where the binary system loses specific angular momentum due to mass loss. The minimum value of d/a corresponds to the ejection point coinciding with the inner Lagrangian point (i.e., d/a = b₁/a − q/(1 + q), which gives 0.726 for q = 0.048), while it is reasonable to assume that the maximum value of d/a corresponds to the ejection point coinciding with the CS position (i.e., d/a = 1 − q/(1 + q), which gives 0.954 for q = 0.048). Consequently, in the following we explore the range of d/a between 0.726 and 0.954.

Since we expect q <  0.1 we can approximate the term q/(1 + q) as ≃q in Eq. (16). By differentiating Eq. (16) and some algebraic manipulation we obtain

(20)

As suggested by Rappaport et al. (1987), the term is much smaller than ṁ₂/m₂; neglecting that term we can rewrite Eq. (20) in the compact form

(21)

where

Combining Eqs. (17) and (21) we obtain

(22)

where α only depends on q and Ṗ₋₁₁ is the orbital period derivative in units of 10⁻¹¹.

On the other hand, the source luminosity can be written as L_x = −β(GM₁Ṁ₂)/R_NS, where R_NS is the NS radius. This relation can be opportunely rewritten as erg s⁻¹ for a NS radius of 10 km; combined with Eq. (21) it becomes

(23)

where L₃₇ is the luminosity in units of 10³⁷.

From the study of the dip arrival times we infer an orbital period derivative Ṗ = 1.46 ± 0.03) × 10⁻¹¹ s s⁻¹ and an orbital period P = 3000.662(1) s. Studying Eq. (22) for q = 0.048, we explored the pairs of β and d/a for which it is possible to obtain the Ṗ value assuming NS masses of 1.4, 1.6, and 1.8 M_⊙.

To explain a Ṗ value of (1.46 ± 0.03) × 10⁻¹¹ s s⁻¹, a non-conservative mass transfer is required; we obtain that the parameter β has to be smaller than 0.35 for each value of the NS mass studied. We find that the closer the ejection point is to the Lagrangian point, the smaller the fraction of accreted matter is. The possible pairs of β and d/a are indicated with the red curve in Fig. 8.

Fig. 8. Constraint of the matter-ejection point position and β for q = 0.048. Left, middle, and right panels describe the constraints for m1 values 1.4, 1.6, and 1.8 M⊙. The red curve represents the orbital period derivative of 1.46 × 10−11 s s−1 for the pairs β–d/a. The black and blue horizontal lines represent the minimum and maximum luminosity of the source during the observations. The minimum and maximum value of d/a correspond to the inner Lagrangian point and the CS position with respect to the center of mass. The observed luminosities suggest that the mass transfer is highly non-conservative and that the ejected matter leaves the system close to the inner Lagrangian point.

A further constraint on the ejection point and on the fraction of matter accreting onto the NS is given by the luminosity determined in Eq. (23). The distance of the source can be estimated using the flux during the type I X-ray bursts observed using Rossi-XTE data, as done by Galloway et al. (2008). Using Eq. (10) in Iaria et al. (2015), given that the flux measured during the bursts showing photospheric radius expansion (PRE) is (2.9 ± 0.4) × 10⁻⁸ erg s⁻¹ cm⁻² and that the NS photospheric radius is r_PRE = 1.1 in units of 10 km for XB 1916-053, we find that the distance D to the source is 8.9 ± 0.6 kpc, 9.3 ± 0.6 kpc, and 9.6 ± 0.7 kpc for a NS mass of 1.4, 1.6, and 1.8 M_⊙, respectively. Adopting the 0.1–100 keV unabsorbed flux obtained from spectra A, B, C, and D, we infer that the luminosity varies between 0.41 × 10³⁷ and 1.46 × 10³⁷ erg s⁻¹ for a NS mass of 1.4 M_⊙, between 0.45 × 10³⁷ and 1.6 × 10³⁷ erg s⁻¹ for a NS mass of 1.6 M_⊙, and between 0.48 × 10³⁷ and 1.7 × 10³⁷ erg s⁻¹ for a NS mass of 1.8 M_⊙ by taking into account the error associated with the distance and considering a 10% relative error for the flux. A similar ranges of luminosities were also observed by Boirin et al. (2000), which confirmed that the source belongs to the Atoll class. The 0.1–100 keV unabsorbed luminosities with the corresponding uncertainties are shown in Table 9 for spectra A, B, C, and D.

We plotted the boundary values of luminosity for q = 0.048 in Fig. 8; since the luminosity depends on β but not on the quantity α = (d/a)² (see Eq. (23)), it is represented with horizontal lines. The observed luminosities constrain β to be lower than 15%, independently of the assumed NS mass. The variation in observed luminosity can be explained with a change of the ejection point from which the matter leaves the binary system. At low luminosity the matter leaves the system close to the inner Lagrangian point, and the fraction of mass accretion rate is only 3%, while at higher luminosity the matter leaves the system from a larger distance and β = 0.11. We show a sketch of the top-view geometry (in scale) of XB 1916-053 assuming q = 0.048 and a NS mass of 1.4 M_⊙ in Fig. 9. The left panel shows the whole binary system, while the right panel shows the companion star and the region close the inner Lagrangian point. The two black dashed segments indicate the boundaries of the ejection region measured from L₁, the segment distant 0.065a from L₁. The lower boundary corresponds to a luminosity of 0.41 × 10³⁷ erg s⁻¹, while the upper boundary (0.101a) corresponds to a luminosity of 1.46 × 10³⁷ erg s⁻¹. Similar results are obtained assuming a NS mass of 1.6 and 1.8 M_⊙. Hence, we have shown that abandoning the stringent hypothesis that the matter could leave the system only from the inner Lagrangian point, as supposed by Iaria et al. (2015), a NS mass of 1.4 M_⊙ also describes the observed scenario well.

Fig. 9. Top view in scale of XB 1916-053 for a mass ratio q = 0.048. Left and right panel: whole binary system and a zoom of it close to the inner Lagrangian point, respectively. The blue curve indicates the Roche lobes of the two stars. The dashed and red circle represent the NS Roche-lobe radius and CS Roche-lobe radius. The magenta, black, and green diamond points are the position of the center of mass (CM), the companion star, and the inner Lagrangian point L1. The dotted circle is the radius of the companion star, without thermal inflation. The green circle represents the accretion disk truncated at 0.47a (see text). The two black dashed segments indicate the distances from L1 at which the matter leaves the binary system for a NS mass of 1.4 M⊙ at a luminosity of 0.41 × 1037 erg s−1 and 1.46 × 1037 erg s−1 going from the closest to the farthest.

We note that this scenario assumes that the observed luminosity L_obs is that emitted from the source L_real; however, XB 1916-053 is a dipping source, and so we expect an inclination angle larger than 60°. Since we do not observe eclipses, we can estimate the upper limit on the inclination angle i from the relation i = 90° −arctan(R_L2/a), where R_L2 is the Roche lobe radius of the CS, finding that i <  80°. If the emission is not spherical from the innermost region of the system then L_obs = L_real cos i and we should take into account this amplification factor. By assuming the arbitrary value of i = 70°, the real luminosity varies between 1.2 × 10³⁷ erg s⁻¹ and 4.3 × 10³⁷ erg s⁻¹. These two values of luminosity combined with the necessity to reproduce an orbital period derivative of 1.46 × 10⁻¹¹ s s⁻¹ gives a value of β between 0.08 and 0.34. In this scenario the matter leaves the system at a distance from L₁ close the CS position.

We note that the mass transfer rate is close to the Eddington limit in this evolutive stage of the source. However, the secular evolution predicts that this system should have a mass transfer rate that is two orders of magnitude lower (see, e.g., Heinke et al. 2013) for a conservative mass transfer scenario.

Finally, Lasota et al. (2008) discussed the stability of the helium accretion disk in ultracompact X-ray binary systems. From their estimations, the authors suggested that XB 1916-053 should be transient; however, the source is observed as persistent. The authors showed three possible explanations of that incongruity: (i) the adopted mass accretion rates (bolometric luminosities) could be underestimated, (ii) the outer disk radius is overestimated, (iii) the CS is not a pure-helium star, but still contains some hydrogen. From our analysis all the caveats highlighted by the authors are verified. Lasota et al. (2008) assumed a mass accretion rate of 7.6 × 10⁻¹⁰ M_⊙ yr⁻¹, while we find from the estimated bolometric fluxes that the mass accretion rate is between 3 × 10⁻¹⁰ and 1.2 × 10⁻⁹ M_⊙ yr⁻¹ for a NS mass of 1.4 M_⊙ and NS radius of 10 km. We show that the outer radius of the accretion disk is not truncated at r_T = R_T/a = 0.6/(1 + q), which corresponds to R_T = 0.57 a for q = 0.048; it is smaller because of the 3:1 resonance, and we find r = 0.47 a. Finally, the fraction of hydrogen is close to 0.2 for a NS mass of 1.4 M_⊙, and increases for higher values of NS mass. We summarized the main parameters of the components of XB 1916-053 in Table 8.

In the following we arbitrarily assume a NS mass of 1.4 M_⊙. Using the mass ratio q = 0.048, the CS mass is 0.07 M_⊙, the bloating thermal factor of the CS is b_f = 1.92, and the fraction of hydrogen of the CS is close to 20%. We find that the orbital separation of the binary system is a ≃ 3.54 × 10¹⁰ cm, while the outer radius of the accretion disk is r_disk = 1.7 × 10¹⁰ cm.

Iaria et al. (2015) suggest that the sinusoidal modulation obtained from the LQS ephemeris could due to be the presence of a third body gravitationally bound to the X-ray binary system. Assuming the existence of a third body of mass M₃, the binary system orbits around the new CM of the triple system. The distance of the binary system from the new CM is given by a_x = a_bin sin i = A c, where i is the inclination of the binary plane to the plane of the sky, A is the amplitude of the sinusoidal function obtained from the ephemeris shown in Eq. (1), and c is the speed of the light. We obtain a_x = (3.9 ± 0.4) × 10¹² cm. To infer the mass of the third body we can write its mass function as

where M_bin is the mass of the binary system. Assuming arbitrarily i = 70°, we found M₃ ≃ 45M_J under the hypothesis that the orbit of the third body is coplanar to the binary system.

6.4. Ionized absorber

From the spectroscopic analysis we investigated the nature of the ionized absorber around the compact object.

Initially we estimated the equivalent column densities of the ions producing the absorption lines in spectra A, B, C, and D. To do this we used the relation shown by Spitzer (1978)

where N_j is the equivalent column density for the relevant species, f_ij is the oscillator strength, W_λ is the equivalent width of the line, λ is the wavelength in centimeters, e is the electron charge, m_e is the electron mass, and c is the speed of the light. A similar analysis was done for the dipping source XB 1254-690 (Iaria et al. 2007a) and X 1624-690 (Iaria et al. 2007b) and, recently, for XB 1916-053 using Suzaku data (Gambino et al. 2019). To infer the column densities of the ions we assumed f_ij = 0.416 for all the H-like ions and f_ij = 0.798 for the Fe XXV ion (Verner et al. 1996b), and we used the equivalent widths shown in Table 6; the inferred N_j values are shown in Table 9. We find that the equivalent ion column density is in the range between 10¹⁶ and (7.2 − 8.5) × 10¹⁷ cm⁻².

We estimated the abundances of each ion with respect to the neutral element N_ion/N_el. To this end, we used the values of the neutral hydrogen equivalent column density of the ionized absorber, N_Hξ, shown in Table 5 and we adopted the relation N_Hξ = N_Hξ/N_el N_el/N_ion N_ion, where N_Hξ/N_el is the inverse of a given element with respect to the hydrogen and N_ion is the value of the equivalent column density associated with the ions shown in Table 9. To infer the value of N_ion/N_el we adopted the solar abundance shown by Asplund et al. (2009). The values of N_ion/N_el for spectra A, B, C, and D are shown in Table 9. We find that the ratio N_ion/N_el is 0.2%, 0.4%, 1–2%, 3–4%, (24 ± 10)%, 2.5%, and 10–20% for Ne X, Mg XII, Si XIV, S XVI, Ca XX, FeXXV, and Fe XXVI, respectively. In a scenario where the NS mass is 1.4 M_⊙, the fraction of hydrogen of the CS is close to 20%; because X ≃ 0.7 for solar abundance, we should roughly multiply the adopted abundances by a factor of three. Consequently, the ion population with respect to a given neutral element should be reduced by a factor of 1/3.

To estimate the region where the ionized absorber is located, we used the definition of ξ in ZXIPCF from which we get r² = L_x/(n_Hξ), where L_x is the unabsorbed luminosity in the 0.1–100 keV energy range, r is the distance of the absorber from the central source, ξ is the ionization parameter obtained by our fits, and n_H is the hydrogen atom density of the absorber. By combining the last equation with N_Hξ = n_HΔr, where Δr is the geometrical thickness of the absorber along the line of sight and N_Hξ is the equivalent hydrogen column density associated with the absorber, and by imposing r >  Δr, we obtained an upper limit on r that is r <  L_x/(N_Hξξ). Furthermore, we expect that the optical depth τ = σ_TN_Hξ of the absorber is lower than one, where σ_T is the Thomson cross section. By rearranging r² = L_x/(n_Hξ), we obtain r² = L_xσ_TΔr/(τξ), and consequently r >  (L_xσ_TΔr/ξ)^1/2. By using the best-fit values of N_Hξ and ξ shown in Table 5 and the 0.1–100 keV unabsorbed luminosities shown in Table 9 we are able to constrain the region where the ionized absorber is placed.

We show the distance of the absorber from the NS in units of 10⁹ cm with respect to the thickness Δr in units of 10⁹ cm in Fig. 10. The condition r <  L_x/(N_Hξξ) is represented with a blue line, r >  Δr is indicated with a black line, and r >  (L_xσ_TΔr/ξ)^1/2 with a red line. The gray area gives the pairs of r and Δr satisfying the conditions. Since the outer radius of the accretion disk is 1.7 × 10¹⁰ cm and because the maximum value of r is roughly 5 × 10⁹ cm, we conclude that the ionized absorber is placed in the innermost region of the system.

Fig. 10. Distance r of the absorber, in units of 109 cm, from the NS with respect to its thickness Δr (in units of 109 cm). The gray area indicates the possible pairs of r and Δr. The horizontal dashed line indicates the best-fit value of r obtained assuming that the observed redshift is due to gravitational redshift for a NS mass of 1.4 M⊙; the green area takes into account the associated errors. The black, red, and blue lines are the constraint conditions (see text).

After analyzing the same Chandra data, Trueba et al. (2020) suggests that the observed redshift in the absorption lines is probably due to gravitational redshift, and finds the distance of the absorber from the NS is (2.7 ± 1.7) × 10⁸ cm and (2.3 ± 1.2) × 10⁸ cm for a NS mass of 1.4 M_⊙; these values refer to spectra B and Γ reported by the authors. By adopting the same scenario and using the best-fit values of the observed redshift shown in Table 5 we find that the absorber is placed at (1.1 ± 0.1) × 10⁸ cm, (1.3 ± 0.2) × 10⁸ cm, (1.2 ± 0.2) × 10⁸ cm, and (1.3 ± 0.7) × 10⁸ cm for spectra A, B, C, and D, respectively. Our estimations are almost a factor of two smaller than the value estimated by Trueba et al. (2020), even if the values are compatible with each other. We show the best-fit value of r and the associated errors with a dashed line and the green area in Fig. 10. We find that the ionized absorber can reach the NS surface only in spectrum D when the source is dimmer; in the other cases the absorber forms a torus around the compact object. From N_Hξ = n_HΔr we can infer the lower limit of n_H assuming the maximum possible thickness of the absorber; we obtain n_H >  5 × 10¹⁵ cm⁻³, n_H >  2 × 10¹⁵ cm⁻³, n_H >  1.5 × 10¹⁵ cm⁻³, and n_H >  0.4 × 10¹⁵ cm⁻³ for spectra A, B, C, and D.

We repeated the same procedure to estimate the distance of the absorber during the dip using the values in Table 7 and knowing that the 0.1–100 keV unabsorbed luminosity is (0.6 ± 0.1) × 10³⁷ erg s⁻¹. We find a weak constraint for r; it is between 2 × 10⁸ and 1.7 × 10¹⁰ cm (the outer radius of the disk). Because of the low value of the ionization parameter it is possible that the absorber is placed at the outer radius; however, we cannot exclude inner radii. We note that the circularization radius r_circ = 0.0859 a q^−0.426 (valid for q between 0.05 and 1; Hessman & Hopp 1990) is r_circ ≃ 1.1 × 10¹⁰ cm for q = 0.048, compatible with the obtained distances.

Finally, we explored the nature of the O VIII absorption edge observed in the spectra during the persistent emission. Since its optical depth is 0.16, assuming that n_e ∼ n_i, we can write τ = N_Hσ_T, where σ_T is the Thomson cross section. We inferred an equivalent hydrogen column density of ∼2.3 × 10²³ cm⁻². This value is compatible with the equivalent hydrogen column associated with the ionized absorber. We conclude that a small fraction of O VIII could be present in the absorber.

7. Conclusions

We updated the orbital ephemeris of the compact dipping source XB 1916-053 by analyzing ten new Chandra observations and one Swift/XRT observation. Three new dip arrival times were extracted, allowing us to extend the available baseline from 36 to 40 years. Our analysis definitively excludes several models of ephemeris suggested previously, leaving as a more likely candidate the quadratic orbital ephemeris including a large periodic modulation close to 25 years (LQS ephemeris in this work).

The quadratic term of the orbital ephemeris implies the presence of an orbital period derivative of 1.46(3) × 10⁻¹¹ s s⁻¹, and allows us to refine the orbital period measurement to 3000.662(1) s. We confirm that such a fast orbital expansion implies that mass transfer should be highly non-conservative. Moreover, we find that for a NS mass of 1.4 M_⊙ it is possible to obtain a compatible orbital period derivative and luminosity values, releasing the stringent hypothesis that the mass leaves the system at the inner Lagrangian point. For a NS mass of 1.4 M_⊙ we find that the ejection point is distant from the inner Lagrangian point between 3 × 10⁸ and 1.5 × 10⁹ cm.

We show that the mass ratio of the system is q = 0.048. This value is obtained from the observed infrahump period, as already discussed in the literature by Hu et al. (2008), and from the observed superhump period, which introduces the pressure term due to a spiral wave to explain the apsidal precession period of the accretion disk. Our study allows us to explain the infrahump and superhump periods (Chou et al. 2001), assuming that the accretion disk has a prograde apsidal precession of 3.9 days and a retrograde nodal precession of 4.86 days. The outer radius of the accretion disk is truncated where the 3:1 resonance occurs (i.e., 1.7 × 10¹⁰ cm for a NS mass of 1.4 M_⊙).

Assuming a NS mass of 1.4 M_⊙, we obtain that the CS mass is 0.07 M_⊙. The thermal bloating of the CS is 1.92 and its percentage of content hydrogen is 20%.

We explain the periodic modulation observed in the orbital ephemeris as the presence of a third body orbiting around the binary system. Assuming a co-planar orbit with respect to the binary system, an inclination angle of 70°, and a hierarchical triple system, we find that the mass of the third body is M₃ ≃ 45 Jovian masses.

To perform our spectroscopic analysis we also included an old Chandra observation, previously analyzed by Iaria et al. (2006). We observed a large variation in the count rate for the different observations. Using a Comptonized component for the continuum we estimated that the 0.1–100 keV unabsorbed luminosity varies from 0.41 × 10³⁷ erg⁻¹ to 1.46 × 10³⁷ erg⁻¹ during the persistent emission.

The spectra show the presence of prominent absorption lines associated with Ne X, Mg XII, Si XIV, SXVI, and Fe XXVI. To fit these lines we included in the model the component ZXIPCF that assumes the presence of an ionized absorber between the central object and the observer along the line of sight. We find that the equivalent hydrogen column density associated with the ionized absorber ranges from 9 × 10²² cm⁻² to 2 × 10²³ cm⁻², going from the dimmest to the brightest spectrum. The ionization parameter is close to log(ξ) = 4.3. From our diagnostic, under the assumption that the observed redshift of the line produced in the absorber is due to gravitational redshift, as recently discussed by Trueba et al. (2020), we find that the absorber has a hydrogen atom density n_H higher than 10¹⁵ cm⁻³ and it is located far from the central source at 1 × 10⁸ cm for a NS mass of 1.4 M_⊙.

We extracted the dip spectrum finding that the equivalent hydrogen column density associated with the ionized absorber is 6 × 10²³ cm⁻² and the ionization parameter is close to log(ξ) = 2.8. Finally, we find that during the dip the ionized absorber is placed at a greater distance, which could be compatible with the edge of the accretion disk (1.7 × 10¹⁰ cm) even if we cannot exclude that it is placed at the circularization radius (1.0 × 10¹⁰ cm).

Acknowledgments

This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. The authors acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.0, from INAF mainstream (PI: T. Belloni; PI: A. De Rosa) and from the HERMES project financed by the Italian Space Agency (ASI) Agreement n. 2016/13 U.O. RI and TDS acknowledge the research grant iPeska (PI: Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016.

References

All Tables

All Figures

	Fig. 1. First-order MEG+HEG light curves of the observations in the 0.3–10 keV energy range. The bin time is 128 s.
In the text

	Fig. 2. Top panel: delays vs time (red points). The black curve is the best fit of the LQ function. The brown curve is the best fit of the LQS function. Bottom panel: residuals in units of σ vs time, corresponding to the LQS function
In the text

	Fig. 3. First-order MEG+HEG folded light curves in the 0.3–10 keV energy range adopting the LQS ephemeris. An ad hoc phase shift of 0.5 is imposed. The label in the plot shows the correspondence between colors and observations. The blue dashed horizontal lines indicate the arbitrary count rate threshold of 3.5 c s−1 and 5.5 c s−1 adopted to select the four sets of observations (see text).
In the text

Fig. 4. Unfolded spectrum and residuals adopting Model 1 (left panels) and Model 2 (right panels). The black, red, green, and blue data correspond to the first-order MEG of spectra A, B, C, and D, respectively. The cyan, purple, yellow, and orange data correspond to first-order HEG spectrum of spectra A, B, C, and D, respectively. The prominent absorption lines in residuals of Model 1 (bottom left panel) are associated with the presence of Ne X, Mg XII, Si XIV, S XVI, and Fe XXVI ions.

In the text

	Fig. 5. Unfolded spectrum and residuals of the dip spectrum adopting Model 2. The black and red data correspond to the first-order MEG and HEG spectra, respectively.
In the text

Fig. 6. Characteristic radii of XB 1916-053 in units of orbital separation vs the mass ratio q. Show are the NS Roche-lobe radius (blue curve), the radius of the disk at which a 2:1 and 3:1 resonance could occur (yellow and black curves), the radius associated with the apsidal rate of the disk for a 3:1 resonance without and with the pressure term due to the presence of a spiral wave in the disk (green and brown curve), and the radius associated with the nodal precession period for a 3:1 resonance (red curve).

In the text

	Fig. 7. Pitch angle of the spiral wave vs the normalized sound of speed (black curve). The red horizontal lines correspond to the angles of pitch between 13° and 21°, as obtained by Montgomery (2001) from simulations.
In the text

Fig. 8. Constraint of the matter-ejection point position and β for q = 0.048. Left, middle, and right panels describe the constraints for m1 values 1.4, 1.6, and 1.8 M⊙. The red curve represents the orbital period derivative of 1.46 × 10−11 s s−1 for the pairs β–d/a. The black and blue horizontal lines represent the minimum and maximum luminosity of the source during the observations. The minimum and maximum value of d/a correspond to the inner Lagrangian point and the CS position with respect to the center of mass. The observed luminosities suggest that the mass transfer is highly non-conservative and that the ejected matter leaves the system close to the inner Lagrangian point.

In the text

Fig. 9. Top view in scale of XB 1916-053 for a mass ratio q = 0.048. Left and right panel: whole binary system and a zoom of it close to the inner Lagrangian point, respectively. The blue curve indicates the Roche lobes of the two stars. The dashed and red circle represent the NS Roche-lobe radius and CS Roche-lobe radius. The magenta, black, and green diamond points are the position of the center of mass (CM), the companion star, and the inner Lagrangian point L1. The dotted circle is the radius of the companion star, without thermal inflation. The green circle represents the accretion disk truncated at 0.47a (see text). The two black dashed segments indicate the distances from L1 at which the matter leaves the binary system for a NS mass of 1.4 M⊙ at a luminosity of 0.41 × 1037 erg s−1 and 1.46 × 1037 erg s−1 going from the closest to the farthest.

In the text

Fig. 10. Distance r of the absorber, in units of 109 cm, from the NS with respect to its thickness Δr (in units of 109 cm). The gray area indicates the possible pairs of r and Δr. The horizontal dashed line indicates the best-fit value of r obtained assuming that the observed redshift is due to gravitational redshift for a NS mass of 1.4 M⊙; the green area takes into account the associated errors. The black, red, and blue lines are the constraint conditions (see text).

In the text