© The Authors 2024

1. Introduction

XTE J1710−281, discovered in 1998 by the Rossi X-ray Timing Explorer (RXTE), is a transient low-mass X-ray binary (LMXB) likely associated with the ROSAT source 1RXS J171012.3−280754 (Markwardt et al. 1998). It exhibits high variability, and numerous Type-I bursts have been documented (Markwardt et al. 2001; Galloway et al. 2008), indicating that the compact object is a neutron star (NS). Galloway et al. (2008) estimated, based on the Type-I bursts, that the distance to the source is 12 kpc, assuming accreting matter with cosmic abundances, or 16 kpc, for accreting matter with helium only.

XTE J1710−281 has an orbital period of 3.28 h (Markwardt et al. 2001) as inferred from analysis of the total eclipses observed in its light curve, which also shows regular dipping phenomena, potentially attributed to occultation from the outer regions of the accretion disk. The presence of dips and total eclipses suggests that the system has an inclination angle between 75° and 80° (Frank et al. 1987). The dips in XTE J1710−281 have been studied by Younes et al. (2009); analyzing XMM-Newton data of the source, they observed that the hydrogen column density is ∼4 × 10²¹ cm⁻² during the persistent emission, from 4 × 10²¹ cm⁻² to 8 × 10²³ cm⁻² during shallow dips, and ∼1.4 × 10²³ cm⁻² during deep dips. Analyzing Chandra and Suzaku data, Raman et al. (2018) observed a broad emission line at 0.72 keV, associated with partially ionized iron, and an absorption line in the Fe-K region of the spectrum at 6.6 keV, interpreted as a blend of Fe XIX to Fe XXV transitions.

An accurate study of the eclipses of XTE J1710−281 was performed by Jain & Paul (2011) and Jain et al. (2022), who derived 78 mid-eclipse arrival times by analyzing RXTE, Chandra, Suzaku, XMM-Newton, and AstroSat observations of the source spanning 18 years. The authors estimated the duration of the eclipse to be 420 s on average and, by studying the delays with respect to the orbital cycles, they detected three orbital glitches that current theoretical models cannot explain.

X-ray binaries undergo evolutionary processes driven by various physical mechanisms. One major factor influencing their evolution is mass transfer. This process leads to the redistribution of angular momentum, which impacts the orbital dynamics (van den Heuvel 1994). Moreover, a non-conservative mass transfer can occur as a consequence of different processes, such as radiative evaporation of the secondary star or the expulsion of matter in the form of accretion disk winds, outflows, or jets (Ruderman et al. 1989; Brookshaw & Tavani 1993; Ponti et al. 2012). These mechanisms alter the overall mass content and momentum distribution within the system.

Furthermore, angular momentum lost in the system occurs via two mechanisms: gravitational wave (GR) radiation and magnetic braking (MB) of the tidally locked companion star (CS; Rappaport et al. 1983; Applegate 1992; Applegate & Shaham 1994; Verbunt 1993). These processes directly impact the orbital parameters, contributing to the dynamic evolution of X-ray binaries over time. Consequently, the orbital separation in X-ray binaries can either increase or decrease.

In this study we exploit the mid-eclipse arrival times obtained from Jain & Paul (2011) and Jain et al. (2022) and propose a different interpretation of the orbital residuals of XTE J1710−281. We suggest that the eclipse arrival times are affected by the presence of a third body orbiting the binary system with a mass of 2.7 Jupiter masses and an orbital eccentricity of 0.38.

2. Data analysis

We used the 78 mid-eclipse arrivals times shown by Jain & Paul (2011) and Jain et al. (2022). We started by assuming a reference orbital period of P₀ = 3.2810632 h and a reference time of T₀ = 54112.83194 modified Julian date (MJD), which corresponds to the mid-eclipse time observed by the RXTE observatory during observation ID 91018-01-02-00 (see Table 1 in Jain & Paul 2011). We inferred the delays corresponding to each eclipse time by taking the fractional part of (T_ecl − T₀)/P₀ and multiplying it by the value of P₀ expressed in seconds. We show the eclipse times and the corresponding delays in the top panel of Fig. 1a.

Fig. 1. Mid-eclipse time delays for different orbital solutions in units of seconds. (a) Time delays relative to the reference time T0 = 54112.83194 MJD and the orbital period P0 = 3.2810632 h, with best-fit linear function (top panel) and O−C values in units of σ (bottom panel). (b) Time delays based on the reference time T0 = 54112.831977 MJD and the orbital period P0 = 3.28106339 h, with best-fit curve as per Eq. (1) and O−C values in σ units (bottom panel).

To correct the values of P₀ and T₀, we fitted the delays with a linear function y = m (t − T_s)+q. The term m = ΔP/P gives the correction to P₀, q = ΔT₀ is the correction to the reference time T₀ and, finally, T_s is fixed to 54112.83194 MJD and corresponds to a temporal shift applied to the eclipse times.

From the fit we obtained a χ²(d.o.f.) of 1298(76), m = 0.0052(3) s days⁻¹ and q = 3.2(5) s. The errors are reported at a 68% (1σ) confidence level. We show the linear best-fit in the top panel of Fig. 1a. The observed minus calculated (O−C) values in units of σ displayed in the bottom panel of Fig. 1a show that the correction is not sufficient to describe the temporal evolution of the source as the delays deviate from the best-fit model up to 10σ. However, a sinusoidal modulation seems to be present in the residuals.

We implemented the corrections from the linear fit by obtaining P₀ = 3.281063396(11) h and T₀ = 54112.831977(6) MJD and recalculated the delays corresponding to the mid-eclipse times accordingly. We plot the delays versus time in the top panel of Fig. 1b. To fit the delays, we adopted a linear plus sinusoidal function (hereafter LS function) defined as

$$\begin{array}{c}\hfill y=q+m(t-T_{\mathrm{s}})+Asin\left[2\pi \frac{[(t-T_{\mathrm{s}})-t_0]}{P_{\mathrm{m}}}\right],\end{array}$$(1)

where the sinusoidal term takes into account the modulation observed in the O−C values (bottom panel of Fig. 1a). By fitting the data, we found a χ²(d.o.f.) of 253(73) that translates to an F-test probability of chance improvement of 7.4 × 10⁻²⁶ compared to the linear model. The associated errors of the best-fit parameters were scaled by the factor $\sqrt{{\chi }_{\mathrm{red}}^2}$ to take into account a value of the ${\chi }_{\text{red}}^2$ of the best-fit model larger than one. We found the following best-fit values: q = 0.16(59) s, m = 0.00129(34) s days⁻¹, A = 6.2(7) s, P_m = 5851(576) days (i.e., 16.0 ± 1.6 years) and t₀ = 2285(257) days. By correcting P₀ and T₀ with the best-fit values of m and q, we obtain T₀ = 54112.831979(7) MJD and P₀ = 3.28106345(13) h. We show the best-fit function in red in the top panel of Fig. 1b and the corresponding O−C values in units of sigma in the bottom panel. We report the best-fit values of the parameters in the second column of Table 1.

The ephemeris of the source obtained using the LS function is

$$\begin{array}{cc}\hfill T_{\mathrm{ecl}}=& \mathrm{MJD}(\mathrm{TDB})54112.831979(7)+\frac{3.28106345(13)}{24}N\hfill \\ \hfill & +\frac{6.2(7)}{86400}sin\left[2\pi \frac{(t-54112.83194)-2285(257)}{5851(576)}\right],\hfill \end{array}$$

where N indicates the orbital cycle and t is expressed in MJD.

We recalculated the delays using the best-fit LS values of T₀ and P₀. We then fitted the delays with the LS function, ensuring that the parameters q and m had best-fit values compatible with zero. The other fit parameters yield the same best-fit values as reported in the second column of Table 1. In light of the significantly high χ²/d.o.f. value, attributable to residuals that still exhibit clear deviations from the LS model, we considered the possible scenario in which the sinusoidal modulation may be associated with the presence of a third body orbiting around the binary system with eccentricity non-null in a hierarchical triple system.

The determination of the orbital eccentricity was achieved through the application of a numerical solution to Kepler’s equation. The utilized technique involved modeling the delay modulation with an elliptical orbit. The eccentricity, e, was extracted by employing an iterative approach to solve Kepler’s equation. The iterative solution, known as the Newton-Raphson method, facilitated the computation of the eccentric anomaly E = M − e sin E (where M is the mean anomaly), providing a quantitative measure of e. M is defined as π[(t − T₀ − T_s)−T_p]/P_m, where T_p is the periastron passage time and P_m the period of the delay modulation.

We call the function adopted to fit the delays “LSe”, which includes a constant term, q, a linear term, m, and an eccentric sinusoidal term. We fitted the data, obtaining a χ²(d.o.f.) of 171.8(71) and a value of the parameter m of 3.3(49)×10⁻² s days⁻¹ that is compatible with zero as we expect. To avoid issues related to parameter degeneracy, we fixed the value of m to zero, assessing later its influence on the χ² value.

By fitting the delays with the LSe model we found a χ²(d.o.f.) of 173(72) and Δχ² of 80 with respect to the best-fit obtained using the LS function; the F-test probability of chance improvement is 1.85 × 10⁻⁷, corresponding to a significance of 5.2σ, suggests that a non-null eccentricity is required with high statistical significance. The associated errors of the best-fit parameters were scaled by the factor $\sqrt{{\chi }_{\mathrm{red}}^2}$ to take into account a value of ${\chi }_{\text{red}}^2$ of the best-fit model larger than one. We obtained q = 1.9(1.4) s, A = 6.1(5) s, P_m = 6249(564) days corresponding to 17.1 ± 1.5 years, T_p = −1261(322) days, e = 0.38(17) and ω = 2.5(4) rad, where ω is the argument of periastron. The delays, the best-fit function, and the corresponding O−C values in units of σ are shown in Fig. 2a (top and bottom panel, respectively). The best-fit values of the parameters are shown in the third column of Table 1. Setting the parameter m equal to zero has little influence on the fit; the F-test probability of chance improvement, leaving the parameter free to vary, is on the order of 0.5.

Fig. 2. Mid-eclipse time delays for different orbital solutions in units of seconds. (a) Time delays obtained by adopting T0 = 54112.831979 MJD for the reference time and P0 = 3.28106345 as the orbital period (top panel). The best-fit curve (red) is described in the text. O−C values are in units of σ (bottom panel). (b) Time delays obtained by subtracting the eccentric sinusoidal modulation (top panel). The best-fit curve is described by Eq. (2). Corresponding O−C values are in seconds (bottom panel).

We also explored the correlation of the eccentricity with other parameters but with the parameter m fixed to zero. This correlation is illustrated in Fig. 3, where contour plots for 68% (green), 90% (magenta), and 99% (red) confidence levels are shown across all six panels. Notably, the strongest correlation is observed between eccentricity and parameter q, as well as between the amplitude A of the sinusoidal function and eccentricity e, where e varies within the range of 0.25−0.55. The inclusion of eclipse arrival times distributed over a broader temporal baseline will serve to minimize these correlations, enhancing the reliability of our findings.

Fig. 3. Correlation of eccentricity (e) with other parameters of the LSe function. The correlations are modest, with the most notable observed between parameters q and e. The eccentricity is constrained to be within 0.3 and 0.5 at a 90% confidence level. Contours in green, magenta, and red represent the 68%, 90%, and 99% confidence levels.

To obtain a constraint on the derivative of the orbital period, Ṗ, we subtracted the best-fit eccentric sinusoidal function from the delays. We then fitted the modified delays using the quadratic function

$$\begin{array}{c}\hfill y=q+m(t-T_{\mathrm{s}})+c{(t-T_{\mathrm{s}})}^2,\end{array}$$(2)

where c = Ṗ/(2P₀) in s days⁻². We show the data and the best-fit function in the top panel of Fig. 2b. We obtained the following best-fit values: q = −0.07 ± 0.33 s, m = ( − 0.4 ± 2.0)×10⁻⁴ s days⁻¹ and c = (2 ± 9)×10⁻⁸ s days⁻². Adopting the best-fit value of c we found that Ṗ =(0.7 ± 2.9) × 10⁻¹³ s s⁻¹. We show the O−C residuals with respect to this model in the bottom panel of Fig. 2b.

3. Discussion

3.1. The jittered behavior in the eclipse arrival times

In our analysis, we investigated the total eclipse arrival times derived by Jain & Paul (2011) and Jain et al. (2022), who utilized a step-and-ramp function to model the eclipse shape in the light curves of XTE J1710−281. However, it is worth noting that the eclipse shape may vary from one eclipse to another, leading to changes in its ingress/egress and duration times, or solely its ingress/egress time. The ingress, egress, and eclipse durations show a jittered behavior on the order of 15 s in EXO 0748−676 (Wolff et al. 2002), close to 5 s in MXB 1659−298 (Iaria et al. 2018) and up to 10 s in AX J1745.6−2901 (Ponti et al. 2017). We show in the bottom panel of Fig. 2b that the jittered behavior is close to 5 s for XTE J1710−281.

Wolff et al. (2007) proposed that the magnetic activity of the CS can generate extended coronal loops above the CS’s photosphere, accounting for the observed jitters. Furthermore, Ponti et al. (2017) present a different scenario for AX J1745.6−2901, where jitters were observed in the ingress and egress, while the eclipse duration remained relatively constant. They hypothesized that matter ejected from the accretion disk can interact with the CS’s atmosphere, causing displacement and hence delays in the ingress and egress times.

Both scenarios suggest that the presence of magnetic activity and/or the atmosphere of the CS may slightly alter the estimates of the eclipse arrival times. This can increase the reduced chi-square, which may be larger than 1, as obtained for our best-fit model.

3.2. Constraints on the companion star mass

XTE J1710−281 shows total eclipse and dips; therefore, its inclination angle must be between 75° and 80° (see Younes et al. 2009, and references therein). Moreover, the duration of the eclipse ΔT_ecl is on average about 420 s (Jain & Paul 2011). Using these observational results, we can provide an estimate of the mass ratio q = M₂/M₁, where M₁ and M₂ are the masses of the NS and the CS.

Knowing that the eclipse duration is ΔT_ecl = 420 s, we can estimate the size of the occulted region x (see Fig. 5 in Iaria et al. 2018) by using the expression

$$\begin{array}{c}\hfill \frac{2\pi }{P}a=\frac{2x}{\mathrm{\Delta }{T}_{\mathrm{ecl}}},\end{array}$$(3)

where P and a are the orbital period and the orbital separation of the system. The value of x depends on a, which is unknown. The angle, θ, complementary to the inclination angle i of the binary system, can be estimated using the following relationship:

$$\begin{array}{c}\hfill tan\theta ={\left[\frac{R_2^2-x^2}{a^2-(R_2^2-x^2)}\right]}^{1/2},\end{array}$$(4)

where R₂ is the CS radius. In light of the system accreting via the inner Lagrangian point, the CS fills its lobe, resulting in the radius of the CS coinciding with the Roche lobe radius of the CS, which is given by Eggleton (1983)

$$\begin{array}{c}\hfill R_{\mathrm{L}2}=a\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})}·\end{array}$$(5)

By substituting Eqs. (3) and (5) into Eq. (4), the angle, θ, depends only on P, which is the known orbital period, on ΔT_ecl, which is 420 s, and on q, our unknown variable. We show the dependence of the inclination angle on q in Fig. 4.

Fig. 4. Variation in the inclination angle with the mass ratio, q (red curve), as described by Eq. (4).

Knowing that the inclination angle of the system ranges between 75° and 80°, we deduce that the mass ratio q falls within the range of 0.1−0.3. Assuming a NS mass of 1.4 M_⊙, the CS mass ranges between 0.14 M_⊙ and 0.43 M_⊙ for an inclination angle of 80° and 75°, respectively.

We can estimate the CS mass with further precision using the mass–radius relationship for stars in thermal equilibrium derived from the study of cataclysmic variables in Knigge et al. (2011). Below, we adopt the mass–radius relationship for binary systems with an orbital period greater than 3.18 h (the orbital period of XTE J1710−281 is close to 3.28 h). The relationship is

$$\begin{array}{c}\hfill R_2=0.293\pm 0.010{\left(\frac{M_2}{M_{\mathrm{conv}}}\right)}^{0.69\pm 0.03}{R}_{\odot },\end{array}$$(6)

where M_conv has a value of 0.20 ± 0.02 M_⊙ and represents the mass of the convective region of the CS. This region fills its Roche lobe, R_L2, so we can assume that R₂ = R_L2. Since we estimated that q is between 0.1 and 0.3, we can use the expression for the Roche lobe radius proposed by Paczyński (1971). Combining the expression for the Roche lobe radius with Kepler’s third law, we find that

$$\begin{array}{c}\hfill R_{\mathrm{L}2}=0.233m_2^{1/3}{P}_{\mathrm{h}}^{2/3}{R}_{\odot },\end{array}$$(7)

where m₂ is the CS mass in units of solar masses, and P_h is the orbital period in hours. Combining Eqs. (6) and (7), we find M₂ = 0.22 ± 0.07 M_⊙.

Assuming a NS mass of 1.4 M_⊙, we infer that q = 0.16 ± 0.05 corresponds to an inclination angle of $i=78.1{}_{-1.2}^{+1.5}°$. For a NS mass of 2 M_⊙ we obtain q = 0.11 ± 0.03 corresponding to an inclination angle of $i=79.6{}_{-1.2}^{+1.6}°$. Using Kepler’s third law, we estimate that the orbital separation of XTE J1710−281 is a = 9.1 × 10¹⁰ cm for an inclination angle of 78°, M₁ = 1.4 M_⊙, and M₂ = 0.22 M_⊙.

3.3. Self-consistence test for the thermal equilibrium of the CS

We can use the bolometric X-ray luminosity of the source to estimate whether the CS is in thermal equilibrium. To do this, we imposed that the Kelvin-Helmholtz timescale, τ_KH (corresponding to the characteristic time that a star takes to reach thermal equilibrium), is equal to or less than the mass transfer timescale, τ_Ṁ. The Kelvin-Helmholtz timescale is given by the expression

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=3.1\times {10}^7{\left(\frac{M_2}{M_{\odot }}\right)}^2\left(\frac{R_{\odot }}{R_2}\right)\left(\frac{L_{\odot }}{L_2}\right)\mathrm{year}\end{array}$$(8)

(Verbunt 1993). Adopting the mass-luminosity relation for M-type stars proposed by Neece (1984) and given by the expression L₂/L_⊙ = 0.231(M₂/M_⊙)^2.61, along with the mass–radius relation from Eq. (6), we obtain ${\tau }_{\text{KH}}\simeq 1.5\times {10}^8\hspace{0.17em}{m}_2^{-1.3}$ year. The mass-transfer timescale is given by

$$\begin{array}{c}\hfill {\tau }_{\dot{M}}=\frac{M_2}{{\dot{M}}_2}=\frac{G{M}_1{M}_2}{L_{\mathrm{X}}{R}_{\mathrm{NS}}},\end{array}$$(9)

where Ṁ₂ is the mass transfer rate, G is the gravitational constant, L_X is the X-ray bolometric luminosity and R_NS is the NS radius. In the last equality of Eq. (9), we imposed that the mass transfer rate is equal to the mass accretion rate in the scenario of conservative mass transfer. This point is further discussed later on. We find that ${\tau }_{\dot{M}}\simeq 3.7\times {10}^{53}{m}_2{L}_{\mathrm{X}}^{-1}$ s for a NS radius of 10 km and a NS mass of 1.4 M_⊙.

By imposing that τ_KH ≤ τ_Ṁ, we obtain m₂ ≥ (1.28 × 10⁻³⁸ L_X)^1/2.3. The unabsorbed luminosity in the 0.2−10 keV energy band of XTE J1710−281, during the persistent emission, was estimated to be L_X ≃ 2.4 × 10³⁶ erg s⁻¹ for a distance to the source of 16 kpc (Younes et al. 2009). For this luminosity value, we determine that the CS is in thermal equilibrium when m₂ ≥ 0.22, which is the CS mass determined by our calculations.

3.4. The orbital period derivative of XTE J1710–281

Drawing upon theoretical frameworks, the mass-transfer rate Ṁ₂ within the context of long-term orbital evolution can be expressed as

$$\begin{array}{c}\hfill {\dot{m}}_{-8}=35{(3n-1)}^{-1}{m}_2\frac{{\dot{P}}_{-10}}{P_{5\mathrm{h}}},\end{array}$$(10)

where ṁ₋₈ is the mass transfer rate in units of 10⁻⁸ M_⊙ yr⁻¹, n is the mass–radius index of the CS, Ṗ₋₁₀ is the orbital period derivative in units of 10⁻¹⁰ s s⁻¹ and P_5h is the orbital period in units of 5 h (see Burderi et al. 2010, and references therein).

Assuming a conservative mass-transfer Ṁ₁ = Ṁ₂ and considering a source luminosity of 2 × 10³⁶ erg s⁻¹, in the hypothesis that the observed luminosity is a good tracer of the mass accretion rate ṁ, we find ṁ = 1.7 × 10⁻¹⁰ M_⊙ yr⁻¹ for a NS mass of 1.4 M_⊙ and a NS radius of 10 km. Consequently, the mass transfer rate is Ṁ₂ = −1.7 × 10⁻¹⁰ M_⊙ yr⁻¹. Using n = 0.69 (adopted in Eq. (6)), m₂ = 0.22 and an orbital period of 3.28 h, we find Ṗ = −1.5 × 10⁻¹³ s s⁻¹. This value is consistent within 1σ with what we obtained from our analysis, that is, Ṗ = (0.7 ± 2.9) × 10⁻¹³ s s⁻¹; therefore, we should expect that the binary system is undergoing a contraction of the orbit.

The orbital period change due to the loss of angular momentum can be associated with the emission of GRs and MB. Below, we explore both cases assuming a conservative mass-transfer scenario. The orbital period changes due to GRs is given by the relationship

$$\begin{array}{cc}\hfill {\dot{P}}_{\mathrm{grav}}=& -1.4\times {10}^{-12}{m}_1{m}_2{m}_T^{-1/3}{P}_{2\mathrm{h}}^{-5/3}\hfill \\ \hfill & \times [(n-1/3)\times (n+5/3-2q)]\mathrm{s}{\mathrm{s}}^{-1},\hfill \end{array}$$(11)

where m₁ is the NS mass in units of M_⊙, m_T is m₁ + m₂ in units of M_⊙, P_2h is the orbital period in units of two hours, and q is the mass ratio of the binary system (see di Salvo et al. 2008, and references therein). Using the values of m₁, m₂, P, and n shown above, we estimate Ṗ_grav ≃ −2.8 × 10⁻¹⁴ s s⁻¹.

The estimated mass of the CS is less than 0.3 M_⊙. Rappaport et al. (1983) proposed that the MB mechanism might not be active, attributing this to the star being fully convective. The lack of a radiative zone implies losing its magnetic field dipolar structure. If this mechanism is inhibited, the only contribution to the orbital period derivative is related to Ṗ_grav. However, both Chen (2017) and Tailo et al. (2018) indicated that the MB mechanism could remain active even in stars with masses below 0.3 M_⊙. In such a scenario, we expect the contribution of MB to the orbital period variation to be significant.

The contribution to the orbital period variation related to MB is given by the relation

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{MB}}={\dot{P}}_{\mathrm{grav}}{T}_{\mathrm{MB}},\end{array}$$(12)

where $T_{\mathrm{MB}}=41.6{(f/k_{0.277})}^{-2}{m}_2^{1/3}{m}_1^{-4/3}{P}_{5\mathrm{h}}^2$ (see Eqs. (2) and (3) in Burderi et al. 2010, and references therein). The gyration radius k is in units of 0.277, and since the CS has a mass of 0.22 M_⊙, we adopt a value k = 0.45 as reported by Wadhwa et al. (2024). The parameter f can range between 0.78 and 1.73 because it is model dependent (see the discussion in Burderi et al. 2010), and we estimate T_MB for the limit values of f finding that T_MB ranges between 5.8 and 29. The corresponding orbital period derivative associated with the MB, Ṗ_MB, is therefore between −8.2 × 10⁻¹³ s s⁻¹ and −1.6 × 10⁻¹³ s s⁻¹. The contribution of the angular momentum loss via GR and MB to the orbital period derivative yields Ṗ between −8.5 × 10⁻¹³ s s⁻¹ and −1.9 × 10⁻¹³ s s⁻¹, still compatible with what is obtained from our analysis adopting T_MB = 5.8.

3.5. The hierarchical triple system scenario

From the delays fitting, we obtain that the statistically most significant solution involves a sinusoidal modulation with an eccentricity e of 0.38 ± 0.17 and a modulation period of 17.1 ± 1.5 years. The delay modulation could be caused by a third body orbiting around XTE J1710−281, where the modulation period corresponds to the revolution period of the third body around the binary system. We can estimate the mass of the third body under the assumption that its orbit is co-planar with that of the binary system. We estimate the separation between the center of mass (CM) of the binary system and that of the triple hierarchical system to be a_x sin i = Ac, where c represents the speed of light, and A is the amplitude of the modulation; a_x sin i = 1.8 × 10¹¹ cm in our case. The mass M₃ of the third body is given by

$$\begin{array}{c}\hfill \frac{M_{3}sini}{{(M_3+M_{\mathrm{bin}})}^{2/3}}={\left(\frac{4{\pi }^2}{G}\right)}^{1/3}\frac{a_{x}sini}{P_{\mathrm{m}}^{2/3}},\end{array}$$(13)

where M_bin is the mass of the binary system, G the gravitational constant and P_m the periodic modulation inferred from the delays fit. For an inclination angle of 78° (we assume that the third body orbit lies on the binary system orbital plane), M₁ = 1.4 M_⊙ and M₂ = 0.22 M_⊙, we infer that M₃ ≃ 2.7 M_J, where M_J indicates the Jupiter mass. The separation between the third body and the CM of the hierarchical triple system is a₃ = a_x sin i M_bin/M₃ ≃ 1.1 × 10¹⁴ cm, corresponding to 7.5 AU, approximately equivalent to a distance similar to that between the Sun and halfway between Jupiter and Saturn. Finally, at the periastron passage, the separation between the CM and the third body is d = a_x sin i (1 − e) M_bin/M₃ = 7.1 × 10¹³ cm, a distance significantly exceeding the orbital separation of the binary system, represented by a = 9.1 × 10¹⁰ cm.

The existence of third celestial bodies in orbit around binary systems, particularly those containing a compact object, is being increasingly confirmed by recent discoveries. Notable among these findings is the proposed detection of a third body with a mass of approximately 45 M_J (0.043 M_⊙) by Iaria et al. (2015, 2021) orbiting the dipping source XB 1916−053. Concurrently, Iaria et al. (2018) have revealed a 22 M_J mass third body orbiting the eclipsing source MXB 1659−298. The precedent for such detections was arguably set by Sigurdsson (1993), who proposed the existence of a sub-Jovian mass planet in orbit around the binary system of the millisecond radio pulsar PSR 1620−26 within the globular cluster M 4. Subsequent Hubble Space Telescope observations allowed this model to be refined, positing the third body as a planet with a mass of 2.5 ± 1.0 M_J in orbit around a binary system of a millisecond pulsar and a white dwarf companion (Sigurdsson et al. 2003).

Bailes et al. (2011) discovered a Jupiter-sized chthonian body composed primarily of carbon and oxygen in a close orbit around the millisecond pulsar PSR J1719−1438, at a distance of 0.004 AU. Such observations are in line with the seminal discovery of planets orbiting PSR 1257+12 Wolszczan & Frail (1992), which has spawned a multitude of hypotheses regarding planetary formation in the vicinity of pulsars.

Podsiadlowski (1993) discussed that planetary bodies associated with NSs are supposed to form during one of three distinct periods: (i) within a protoplanetary disk during the initial star formation (first-generation); (ii) from a fallback disk composed of debris of the supernova explosion (second-generation); (iii) or within an accretion disk that forms from matter transferred to the NS from the CS (third-generation). Planets from the first generation, if they exist, are presumed either to have been destroyed or to have had their trajectories severely altered in the wake of a supernova event. The existence of planets orbiting millisecond pulsars challenges the hypothesis of second-generation formation, as the spin acceleration implies a proximal CS, within an orbital distance of roughly 1 AU, and this would likely perturb planets orbits. This is corroborated by the instance of PSR B1257+12, which is a fully recycled pulsar, implying its companion is a low-mass star and had experienced Roche lobe overflow during its main-sequence life stage, necessitating an orbital proximity not exceeding 1 AU.

Therefore, the prevailing hypothesis for the formation of planets orbiting isolated millisecond pulsars supports the idea of third-generation formation, where planets form from the residual matter of the CS. However, this evolutionary channel clashes with the findings of this study, where we have an X-ray binary system undergoing mass transfer from a main-sequence CS. Similarly, this evolutionary path cannot easily explain the three-body system of the millisecond pulsar PSR 1620−26, where the main CS is a white dwarf. Therefore, one possibility is that these planets are solitary planets that were gravitationally captured by the binary system (Podsiadlowski et al. 1991), which is very unlikely for Galactic field systems. Alternatively, it is possible that these planets emerge from a circumbinary disk, formed from material expelled by the CS during a period of non-conservative mass transfer of the LMXB (Tavani & Brookshaw 1992).

3.6. Limitations of the Applegate mechanism in explaining orbital period modulation in XTE J1710–281

Given that the CS is a low-mass star (0.22 M_⊙), we expect its magnetic activity to be intense. However, XTE J1710−281 exhibits persistent emission with mass transfer via Roche lobe overflow, so the luminosity from the outer regions of the accretion disk may overshadow that of the CS in the visible band, hindering a detailed study of the CS Ratti et al. (2010). Consequently, the magnetic properties of the CS in XTE J1710−281 are currently unknown.

Nevertheless, an alternative explanation of the modulation of approximately 17 years observed in the orbital delays could be given by a gravitational coupling of the orbit with variations in the shape of the magnetically active CS. These variations are believed to result from the torque exerted by the magnetic activity associated with a subsurface magnetic field in the CS, interacting with its convective envelope. The convective envelope initiates a periodic exchange of angular momentum between the inner and outer regions of the CS, leading to alterations in its gravitational quadrupole moment (Applegate 1992; Applegate & Shaham 1994).

In this instance, the deduced periodicity of 6249 days and amplitude of 6.1 s translates to an orbital period variation of ΔP/P ≃ 2.2 × 10⁻⁵. Under the reasonable assumption that the CS fills its Roche lobe, and consequently, its radius coincides with the Roche-lobe radius, we can calculate the angular momentum transfer required to induce the observed orbital period change to be approximately ΔJ ≃ 4.8 × 10⁴⁴ g cm² s⁻¹.

The non-synchronicity of the companion, expressed by the ratio ΔΩ/Ω, is 6.7 × 10⁻⁵, where Ω is the orbital angular velocity of the binary system and ΔΩ represents the variation in the orbital angular velocity required to induce the change in the orbital period, ΔP. The value of ΔΩ/Ω is obtained assuming that a thin shell of mass, M_s is 10% of the CS mass, where the shell represents the limited mass in the outer part of the star that governs the quadrupole moment. For M_s mass larger than 10%, the angular momentum transfer mechanism ceases to function (Applegate 1992).

The variable part of the CS luminosity necessary to fuel the variations in the gravitational quadrupole is ΔL = 2.9 × 10⁻⁵ L_⊙. To estimate the required fraction of the CS luminosity to achieve the observed modulation in the delays, we deduce the CS luminosity using the relationship L/L_⊙ = 0.231 (M₂/M_⊙)^2.61, applicable to low-mass stars. We find that L ≃ 4.4 × 10⁻³ L_⊙ and ΔL ≃ 0.007. To achieve these brightness variations, we would need a shell mass of 0.7% of the CS mass, while the available budget should be around 10% of the CS’s luminosity. This result appears odd with the model, so we conclude that this is unlikely to be the correct explanation.

4. Conclusions

We used the mid-eclipse times reported by Jain & Paul (2011) and Jain et al. (2022) to offer a different interpretation of the orbital residuals of XTE J1710−281. By fitting the delays, we find a periodic modulation of close to 17 years described by an eccentric sinusoidal modulation with an eccentricity of about 0.38 and an amplitude of 6.1 s.

The likely scenario describing the sinusoidal modulation involves a third body orbiting the binary system with a revolution period of 17 years and an orbital eccentricity of 0.38. From these parameters, we deduce that the mass of the third body is 2.7 M_J. The presence of a third body around a binary system has also been discussed for the eclipsing LMXB MXB 1659−298 (Iaria et al. 2018); this third body has a mass exceeding 21 M_J. Moreover, the presence of a third body has also been considered by Iaria et al. (2015) for the ultra-compact LMXB XB 1916−053; in this case, the mass of the third body would exceed 45 M_J.

It is important to emphasize that the presence of a third orbiting body does not contract or expand the orbital period of a binary system. Indeed, the motion of the third body influences the position of the CM of the binary system, thereby affecting the arrival times of eclipses. The derivative of the orbital period obtained in this work is not affected at zero order by the presence of the third body; it reflects the expected evolution of an LMXB system given by the mass transfer and the loss of angular momentum via MB and gravitational radiation in the case of a conservative mass transfer.

Our interpretation of the results, statistically equivalent to that obtained by Jain et al. (2022), rules out the presence of discontinuities in the derivative of the orbital period, for which, to date, there is no theoretical model for their interpretation. Instead, our result aligns with the evolution of a binary system in which the CS of 0.22 M_⊙ is in thermal equilibrium, and conservative mass transfer tends to contract the binary system.

Acknowledgments

The authors acknowledge financial support from PRIN-INAF 2019 with the project “Probing the geometry of accretion: from theory to observations” (PI: Belloni). W.L. conducted this research during, and with the support of, the Italian national inter-university PhD program in Space Science and Technology.

References

All Tables

All Figures

Fig. 1. Mid-eclipse time delays for different orbital solutions in units of seconds. (a) Time delays relative to the reference time T0 = 54112.83194 MJD and the orbital period P0 = 3.2810632 h, with best-fit linear function (top panel) and O−C values in units of σ (bottom panel). (b) Time delays based on the reference time T0 = 54112.831977 MJD and the orbital period P0 = 3.28106339 h, with best-fit curve as per Eq. (1) and O−C values in σ units (bottom panel).

In the text

Fig. 2. Mid-eclipse time delays for different orbital solutions in units of seconds. (a) Time delays obtained by adopting T0 = 54112.831979 MJD for the reference time and P0 = 3.28106345 as the orbital period (top panel). The best-fit curve (red) is described in the text. O−C values are in units of σ (bottom panel). (b) Time delays obtained by subtracting the eccentric sinusoidal modulation (top panel). The best-fit curve is described by Eq. (2). Corresponding O−C values are in seconds (bottom panel).

In the text

	Fig. 3. Correlation of eccentricity (e) with other parameters of the LSe function. The correlations are modest, with the most notable observed between parameters q and e. The eccentricity is constrained to be within 0.3 and 0.5 at a 90% confidence level. Contours in green, magenta, and red represent the 68%, 90%, and 99% confidence levels.
In the text

	Fig. 4. Variation in the inclination angle with the mass ratio, q (red curve), as described by Eq. (4).
In the text