A comprehensive study of orbital evolution of LMC X–4: Existence of a second derivative of the orbital period

We report here results from pulse arrival time delay analysis of the eclipsing high mass X-ray binary pulsar LMC X–4 using observations made with the Rossi X-ray Timing Explorer, XMM -Newton, NuSTAR and AstroSat . Combining the orbital parameters determined from these observations with the historical measurements dating back to 1998, we have extended the 𝑇 𝜋 / 2 epoch history of LMC X–4 by about 4600 binary orbits spanning about 18 years. We also report mid-eclipse time measurements ( 𝑇 𝑒𝑐𝑙 ) using data obtained from wide-field X-ray monitors of MAXI -GSC and Swift -BAT. Combining the new 𝑇 𝜋 / 2 and 𝑇 𝑒𝑐𝑙 estimates with all the previously reported values, we have significantly improved the orbital evolution measurement, which indicates that the orbital period is evolving at a time scale ( 𝑃 orb / (cid:164) 𝑃 orb ) of about 0.8 Myr. For the first time in an accreting X-ray pulsar system, we confirm the existence of a second derivative of the orbital period, having an evolution time scale ( (cid:164) 𝑃 𝑜𝑟𝑏 / (cid:165) 𝑃 𝑜𝑟𝑏 ) of about 55 yr. Detection of a second derivative of the orbital period in LMC X–4 makes its orbital evolution timescale more uncertain, which may also be true for other HMXBs. Independent solutions for the orbital evolution measurement using the mid-eclipse data and the pulse timing data are consistent with each other, and help us put an upper limit of 0.009 on the eccentricity of the binary system.


INTRODUCTION
It has been long realized that most of the bright galactic X-ray sources occur in binary systems (Verbunt 1993;Sana et al. 2012), and binary evolution plays a key role in understanding their stellar evolution (van den Heuvel 1994).Several key astrophysical phenomena, such as the formation of double compact binary (double black holes, black hole-neutron star and double neutron star systems), followed by the merger of the two stellar components, production of short gamma-ray bursts and eventually a possible gravitational wave detection requires a comprehensive understanding of the interaction between the binary components (Belczynski et al. 2002).Measurements of the orbital period decay are often used to place limits on the mass transfer and/or rate of mass loss from the system (Deeter et al. 1981).Orbital decay of the famous Hulse-Taylor binary pulsar provided an unprecedented insight into the loss of energy due to gravitational waves, consistent with the general theory of relativity (Taylor & Weisberg 1982).
Being progenitors of double compact objects, the orbital period of X-ray binaries and their evolution has been extensively studied.Important ingredients of such studies include the effect of mass exchange between binary components and mass loss from the binary system on the orbital parameters.The orbital evolution mechanisms largely include mass transfer (conservative as well as non-conservative) from ★ E-mail: chetanajain11@gmail.com (CJ) † E-mail: rahul1607kumar@gmail.com (RS) the companion star to the compact object (van den Heuvel 1994), tidal dissipation in close binaries (Lecar et al. 1976;Zahn 1977), loss of orbital angular momentum due to stellar winds (Brookshaw & Tavani 1993), gravitational wave radiation (Verbunt 1993) and magnetic activity associated with the companion star (Wolff et al. 2009;Jain & Paul 2011;Jain et al. 2022).
The current work is an accurate and most up-to-date study of the orbital evolution of LMC X-4 which is an eclipsing high mass Xray binary (HMXB) system located about 50 kpc away in the Large Magellanic Cloud (LMC) (Giacconi et al. 1972).It was discovered in 1972 by the UHURU observatory.This system comprises of a 1.25 M ⊙ neutron star in an almost circular orbit around a 14th magnitude OB star (Pakull & Olander 1976;Pesch et al. 1976;Kelley et al. 1983a;van der Meer et al. 2007).The timing variability of LMC X-4 includes 13.5 s coherent pulsations (Kelley et al. 1983a), ∼26 mHz quasi-periodic oscillations (Rikame et al. 2022;Sharma et al. 2023), about an hour long flaring episodes (Epstein et al. 1977;Skinner et al. 1980;Levine et al. 1991;Beri & Paul 2017), about 5 hr long X-ray eclipse (White 1978), an orbital period of ∼1.4 d (Li et al. 1978;White 1978) decaying at a rate of about 10 −6 yr −1 (Falanga et al. 2015) and 30.5 d intensity variation due to a precessing tilted accretion disk (Lang et al. 1981;Paul & Kitamoto 2002;Molkov et al. 2015).
The orbital evolution in eclipsing X-ray binary pulsars is measured using two well established techniques, viz., timing of the X-ray eclipses and the measurement of the orbital epoch using pulse arrival time delays across the orbital phases.The eclipse timing technique is based on the hypothesis that since the binary components revolve around their center of mass in Keplerian orbits, they are expected to eclipse each other after regular intervals of time.But if the orbit of the binary is perturbed, then the occurrence of eclipses is delayed (increasing orbit) or is advanced (decreasing orbit) in time.Therefore, in this technique, the mid-eclipse times (  ) are measured over a sufficiently long time base and time connecting them gives estimates on the orbital period evolution of the system.This method has been used to determine the orbital evolution in several eclipsing X-ray binaries, such as, 4U 1822-37 (Jain et al. 2010), AX J1745.6-2901(Ponti et al. 2017), EXO 0748-676 (Wolff et al. 2009), XTE J1710-281 (Jain & Paul 2011;Jain et al. 2022), MXB 1658-298 (Jain et al. 2017), 4U 1700-37 (Rubin et al. 1996;Falanga et al. 2015;Islam & Paul 2016) etc.
The pulse arrival time technique (Staubert et al. 2009) is based on correcting the arrival time of the pulses for the binary motion of the compact object.In this method, the pulse profiles are produced using the already known pulse period, P spin (using  2 maximisation and epoch folding; Leahy et al. 1983).Assuming  0 is the reference time of the first pulse, and considering non-zero first derivative of the pulse period, ignoring its higher derivatives, the expected arrival time of the  th pulse (  ) as a function of the pulse number (n), is given by Equation 1, Assuming an almost circular orbit (eccentricity,  << 1), the fourth term in this equation corresponds to the pulse arrival time delay due to the orbital motion.Here,   sin  is the projected radius of the orbit and  /2 is the mean orbital longitude of the neutron star and corresponds to the maximum delay in the pulse arrival time.An important limitation of this method is the fact that this technique requires sufficiently long observations, covering at least a significant part of the binary orbit.This method has been used in X-ray binaries, such as SAX J1808.4-3658(Jain et al. 2007;Burderi et al. 2009), Her X-1 (Deeter et al. 1991;Staubert et al. 2009), 4U 1538-52 (Baykal et al. 2006;Mukherjee et al. 2006), Cen X-3 (Kelley et al. 1983b;Raichur & Paul 2010a), SAX J1748.9-2021(Sanna et al. 2016;Sharma et al. 2020) etc.The first measurements of the pulse arrival time delay in LMC X-4 were reported by Kelley et al. (1983a) using the SAS-3 observations.Later, using the pulse timing analysis, Dennerl (1991) and Levine et al. (1991) established upper limits of ∼ 10 −6 yr −1 on the orbital decay.Safi Harb et al. (1996) and Woo et al. (1996) discussed the orbital period decay in LMC X-4 in the context of conservative mass transfer and tidal evolution, superposed by mass loss from the binary system in the form of stellar winds.Using a long RXTE observation of 1998 October, Levine et al. (2000) gave the first definite estimate of the decay rate which was refined by Naik & Paul (2004), Falanga et al. (2015) and Molkov et al. (2015) using data fetched from several X-ray missions.
In this paper, we present an update on the orbital parameters of LMC X-4 using both the methods described above.Historically, the orbital evolution measurements of LMC X-4 have been done using one of these methods or by combining them; but using only a partial set of the available data.For example, Falanga et al. (2015) used only the eclipse data, Levine et al. (2000) and Naik & Paul (2004) used only pulse arrival time data and Molkov et al. (2015) used a part of available data from both techniques.
The paper is organized in the following way.Section §2 gives a description of the observation and the data reduction procedure.In Section §3, the results from the timing analysis of LMC X-4 are presented.Our findings are discussed in Section §4.

OBSERVATIONS
The observation details of the narrow field instruments from which data has been used for the current work are given in Table 1.For all the data analysed in this work, the source position (R.A. (J2000) = 05 ℎ 32  49.555  and dec.(J2000) = −66 • 22 ′ 13.202 ′′ ) was adopted from (Gaia Collaboration 2020) to convert the photon arrival times to the solar system barycenter.As described in the previous section, short observations of LMC X-4 have not been used in the current work for pulse timing analysis.All the light curves were extracted with a bin time of 0.1 s.LMC X-4 was observed with Rossi X-Ray Timing Explorer (RXTE: Bradt et al. 1993) several times between 1996 to 1999, amongst which Levine et al. (2000) analysed the longest (∼ 16 day long) observation of 1998 October to obtain a definite measurement of the orbital decay.For the current work, we have used data from the observation made with the Proportional Counter Array (PCA: Jahoda et al. (1996) in 1996 August (observation ID: 10135-01) and 1999 December (Observation ID: 40064-01).The total span for both of these observations was more than the orbital period of LMC X-4 and for each observation, a useful exposure of ∼60 ks was obtained.We have analysed data collected in the Good Xenon mode having a time resolution of 1 s.The combined 2-20 keV light curve was extracted using the seextract tool of XRONOS sub-package of FTOOLS (Blackburn et al. 1999).The background light curve was extracted from the Standard-2 mode data by using pcabackest with a bright source background model (pca_bkgd_cmbrightvle_eMv20051128).Since a different number of PCUs were operational during both the RXTE observations, therefore, we used the correctlc tool to calculate the equivalent count rate for the simultaneous operation of all five PCUs.The time series was barycenter corrected by using faxbary.
LMC X-4 was observed six times with XMM-Newton (Jansen et al. 2001), out of which, Molkov et al. (2015) have reported   measurements from the first two observations.For the current work, we have analysed the observation of 2003 September which had the longest exposure time (∼ 55 ks) covering about 45% of the binary orbit.XMM-Newton carries three focal plane European Photon Imaging Cameras (EPIC) for three X-ray telescopes, EPIC-MOS1, -MOS2 and -pn (Strüder et al. 2001;Turner et al. 2001).The raw data files were processed using version 20.0.0 of the XMM Science Analysis System (SAS).For the present analysis, we have used 0.5-10 keV data taken with the EPIC-pn detector.The source events were extracted from a circular region of radius 40 ′′ centred on the source position.The background events were extracted from a similar region centred away from the source location.The event arrival times in the background subtracted light curve were corrected to the solar system barycenter using the SAS tool barycen.
The Nuclear Spectroscopic Telescope ARray (NuSTAR: Harrison et al. 2013) is a focusing high-energy X-ray telescope which operates in the 3-79 keV energy band.It comprises of two identical focal plane modules (focal plane modules A (FPMA) and B (FPMB)).For the current work, NuSTAR observation of July 2012 was used with a duration of 62 ks covering ∼50% of the binary orbit.We have used the standard NuSTAR analysis software NuSTARDAS and the latest calibration files (version 20230124) for data reduction and analysis.The clean event files were generated using the nupipeline.Source events were extracted from a circular region of radius 100 ′′ centred at the source position.The background events were extracted from a similar region away from the source location.Barycenter correction was done using barycorr.The light curves from both detectors were added for further analysis.
AstroSat is India's first multi-wavelength (from optical to hard X-rays) astronomical mission.It was launched by Indian Space Research Organization in September 2015 (Agrawal 2006).Large Area X-ray Proportional Counters (LAXPC; Agrawal et al. 2017) is one of the primary instruments onboard the AstroSat.It has a high time resolution of 10 s and covers a broad X-ray spectral band in 3-80 keV.It consists of three co-aligned proportional counters (LAXPC10, LAXPC20 and LAXPC30), with a total effective area of 6000 cm −2 at 15 keV.The AstroSat observation of LMC X-4 (Observation Id G05_115T01_9000000634) made during 2016 August was analyzed for this work.Detailed analysis of this observation has been presented in Sharma et al. (2023).For the current analysis, we used data from LAXPC10 and LAXPC20.LAXPC30 was not used due to gain variability with the instruments (Agrawal et al. 2017;Antia et al. 2017).During this observation, LMC X-4 was observed for a duration of ∼90 ks, covering about 75% of the binary orbit.The Event Analysis (EA) mode data from LAXPC10 and LAXPC20 was processed by using the standard LAXPC software1 (LaxpcSoft: version 3.4.2).The light curves for the source and background from both LAXPC units were extracted from level 1 files by using the tool laxpcl1 and added using lcmath.

TIMING ANALYSIS
The background subtracted and barycenter corrected light curves of observations listed in Table 1 were filtered for the flaring and eclipse phases.Figure 1 shows all these five light curves to highlight the segments used and those that were filtered off for the pulse timing measurements.
3.1  /2 measurements LMC X-4 has an orbital period of ∼1.4 d and semi-amplitude of the arrival time delay due to orbital motion (  sin ) of ∼26.3 s.As a result, the pulse frequency gets modulated by the Doppler effect associated with the orbital motion and the pulses are expected to lose coherence within a few thousand seconds.Assuming a nearly circular orbit, we corrected the photon arrival time for the binary motion.The value of   sin  was fixed to the value taken from Levine et al. (2000) and for each observation,  orb was calculated by extrapolating the orbital solution of Molkov et al. (2015).
For all the pointed observations of RXTE, XMM-Newton, NuSTAR and AstroSat tabulated in Table 1, we searched for the correct orbital ephemeris (T /2 ) around the value extrapolated from Molkov et al. (2015).For this,  /2 was searched over a wide range of 0.2 d on either side of the extrapolated value in fine steps of ∼ 10 −4 d.We used epoch folding and  2 maximisation technique (Leahy et al. 1983) for each trial ephemeris.As an example, Figure 2 shows the variation in  2 over the trial range of ephemeris for the AstroSat observation.This variation in  2 was fit with a Gaussian profile over a narrow range of ±0.006 d (inset of Figure 2).We obtained T /2 = 57630.20621(12)MJD.The T /2 determined from the other observations are listed in Table 2, along with the previously reported measurements.We used the bootstrap method of Boldin et al. (2013) to estimate the error in the measurement of T /2 .Following Sharma et al. (2023), we simulated 1000 light curves and determined the value of T /2 in each of them by the epoch folding technique.The standard deviation of T /2 measured from all the simulated light curves was taken as 1  error on T /2 .
To be doubly sure and to avoid any model dependence, we varied these extrapolated values of  orb by ±2 × 10 −5 d.We did not find any significant effect of this variation on our measurement of  /2 .

𝑇 ecl measurements
LMC X-4 has been continuously monitored with Burst Alert Telescope (BAT) (Barthelmy et al. 2005;Krimm et al. 2013) on-board the Swift observatory (Gehrels et al. 2004) since 2005 and Gas Slit Camera (GSC) on-board the Monitor of All-sky X-ray Image (MAXI: Matsuoka et al. 2009;Mihara et al. 2011) since 2009.Therefore, in order to determine a few more mid-eclipse times, we have used the 18 year long, publicly available 15-50 keV Swift-BAT orbital light curve 2 and 14 year long 2-20 keV MAXI-GSC orbital light curve3 .The photon arrival times were corrected to the solar barycenter reference time by using earth2sun task.
The Swift/BAT light curve was divided into three segments spanning about 6 years, viz.MJD 53416-55680, 55680-57943 and 57943-60205.The MAXI-GSC light curve was divided into two segments, viz.MJD 55054-57032 and 57032-60212.Each of these light curves was folded at local  orb estimated using epoch folding and the mid-eclipse time was estimated using the ramp function.Since all these time segments are relatively long (∼ 2000 d), that the orbital period changed during these intervals.Therefore, we also folded the profiles at local  orb and  orb obtained from Molkov et al. (2015).There was no significant change in the eclipse time measurement.In fact, use of extrapolated orbital period from Molkov et al. (2015) solution resulted in a similar estimate of the mid-eclipse time, clearly indicating that the period evolution is too slow to have any considerable effect during these time intervals.
Figure 3 shows the eclipse profile of the first segment of the MAXI-GSC and Swift-BAT light curve fitted with a ramp function.The mideclipse times for these segments were found to be MJD 56042.990(3)and 54547.2978(11),respectively.Table 3 lists the previously re-   Falanga et al. (2015).

Orbital evolution
In order to determine the orbital evolution in LMC X-4, we referred to a constant orbital period  0 = 1.40837607 d at a reference epoch  0 = 53013.5878MJD (Molkov et al. 2015).Using these numbers, we calculated the orbital cycle () for every measurement listed in Tables 2 and 3.
Initially, the orbital epochs were fit with a linear+quadratic (LQ) model to all the 58 measurements of   and 14 measurements of  /2 , using the relation, The first two terms in this equation, correspond to the case of a constant orbital period.The third term describes the quadratic nature of the orbital decay in LMC X-4.The best fit parameters are given in Table 4 (columns 4 and 5).The best-fit LQ model had reduced  2 of 1.9 (in  ecl measurements) and 28 (in  /2 measurements).To be more conservative and to account for much larger than 1 reduced  2 , we scaled the errors associated with  0 ,  orb and  orb by respective factor of √︃  2 red (Elsner et al. 1980;van der Klis & Bonnet-Bidaud 1989;Iaria et al. 2015).Hereafter, this scaling of errors in  0 ,  orb and  orb has been done for all the models.
We also fitted the LQ model to the combined 72 measurements of  ecl and  /2 .The best-fit model had a  2 of 443 for 69 degrees of freedom (d.o.f).Table 4 (column 6) gives the best-fit parameters of the LQ model applied to these 72 measurements, along with their 1 errors (after multiplying by

√︃
2 red , i.e. ∼2.5).For these combined measurements, we also fitted the linear+quadratic+cubic (LQC) model (Eqn.2) to the orbital epochs, where  orb is the second derivative of the orbital period.The best fit had a  2 of 296 for 68 d.o.f.The inclusion of the cubic term improved the fit statistics by Δ 2 ∼147 for one additional d.o.f with an -test probability of 1.7 × 10 −7 which corresponds to more than 3 significance.We found the second-period derivative to be 2.5(4) × 10 −13 d d −2 and the evolution time scale of the period derivative (  orb /  orb ) to be 0.018(3) yr −1 .Figure 4 shows the Observed-Calculated (O-C) residual of orbital epochs of LMC X-4 relative to a constant orbital period, for both LQ and LQC models.For completeness, in order to observe the individual contribution of the quadratic and cubic components, to the net  2 , Figure 5 shows the residuals (in units of ) w.r.t.LQ model (upper panel) and the LQC model (lower panel).
Amongst the observations with the narrow field instruments analysed in this work, only the observation from XMM-Newton covered a complete X-ray eclipse of LMC X-4.The results from the eclipse timing are already reported in Molkov et al. (2015).However, the errors quoted in Molkov et al. (2015) are as large as 1300 s.This is perhaps due to the post-eclipse dips, which can affect the determination of the eclipse egress profile.The difference between the two  measurement techniques and the two orbital evolution models, with and without a second derivative of the orbital period is shown in Figure 6.The vertical dashed lines in this figure correspond to mideclipse time measurement obtained from the eclipse timing technique (magenta color) by Molkov et al. (2015),  /2 measured from the pulse timing (reported in Table 2) in red color,   calculated from the best-fit LQ model (blue color) and LQC model (green color).The  /2 measurement is closest to the estimates from the LQC model.This inference also supports the fact the pulse arrival time analysis is much more accurate than the mid-eclipse measurements.

DISCUSSION
In this paper, we have revisited a very well-studied HMXB system, LMC X-4 and have made robust estimates of its orbital evolution, by adding new measurements from light curves obtained with the pointed observations made with RXTE, XMM-Newton, NuS-TAR, AstroSat, and long term light curves obtained with Swift-BAT and MAXI-GSC data and thereby spanning the ephemerides history over almost 45 years.We have determined an orbital period of 1.40837655(15) d at MJD 53013.5894(2),decaying at a rate of -5.13(4)×10 −9 d d −1 , with an orbital evolution timescale of about 0.8 Myr.We have also estimated a second derivative of orbital period 2.5(4) × 10 −13 d d −2 , implying that the orbital decay rate is evolving at a time scale of just about 55 years.We discuss these results in the light of, (i) Orbital evolution mechanisms at play in LMC X-4.(ii) Importance of existence of    .(iii) Importance of a non-zero  cos  term.
Based on the orbital modulation, several authors have studied the changes in the orbital period of Cyg X-3 (van der Klis & Bonnet-Bidaud 1989;Kitamoto et al. 1992Kitamoto et al. , 1995;;Singh et al. 2002;Bhargava et al. 2017).Some earlier works based on data before 1990 (van der Klis & Bonnet-Bidaud 1989;Kitamoto et al. 1992), reported a second derivative of the orbital period but detection was marginal and results were inconclusive.The detection significance of  orb crucially depended on the intrinsic scatter in the data.The inclusion of more data (up to 1993) gave a smaller value  orb (Kitamoto et al. 1995).It was even smaller from data up to 2001 (Singh et al. 2002).The cubic fit to the O-C curve was only marginally better than quadratic fit.The latest works by Bhargava et al. (2017) over 45 years of time base are biased towards a secular variation in the orbital period without any requirement of a second derivative.They have however hinted towards short-term local period variations linked with jet emission.
Assuming a conservative mass transfer in LMC X-4, where the neutron star accretes all the matter lost by its companion, the rate of change of the orbital period is given by (van den Heuvel & de Loore 1973), This gives   ∼ −7.6 × 10 −7  ⊙ yr −1 , assuming orbital decay rate of  orb / orb = −1.31× 10 −6 yr −1 , neutron star mass of    = 1.57 ⊙ , and a companion mass of    = 18 ⊙ (Falanga et al. 2015).This estimate is nearly 3 times more than the theoretical mass loss estimate of ∼ −2.4 × 10 −7  ⊙ yr −1 (Falanga et al. 2015); and it exceeds the Eddington mass accretion limit by a large factor.Clearly, other mass loss mechanisms are at work in addition to conservative mass transfer and tidal decay.Assuming a non-conservative mass transfer, where only a fraction of the ejected matter is accreted by the neutron star; following van den Heuvel (1994) and Jenke et al. (2012), and referring to orbital parameters of LMC X-4, we have determined a lower limit on the angular momentum transferred through stellar winds from the companion by using, where  is the ratio of escaping angular momentum per unit mass to the total angular momentum per unit mass and given by  = where  is the semi-major axis of the orbit and   is the radius beyond the L 2 point of the system of escaping material.Assuming   > 1.2, and neutron star mass of    = 1.57 ⊙ , a companion mass of    = 18 ⊙ and mass loss rate of   ∼ −2.4 × 10 −7  ⊙ yr −1 (Falanga et al. 2015), we get a decay rate of > -1.4×10 −6 yr −1 .This is well consistent with the observed    /   value of -1.3×10 −6 yr −1 .
For most binary systems, one can measure the period derivative and a first order estimate for the orbital evolution timescale.However, if the period derivative itself changes over time, then this estimation is inaccurate.We, for the first time, have determined a second derivative of the orbital period, which turns out to be quite short.So the period evolution timescale that is measured today is perhaps not valid in a few decades.If the same is true for the other HMXBs, then the long-term evolution of HMXBs can not be estimated from the current measurements of  orb and  orb .However, LMC X-4 should be observed extensively over the next decade to have a more secure determination of the second period derivative.
For eccentric orbits,  ecl (measured from the eclipse timing technique) generally does not coincide with  /2 (determined from the pulse arrival time technique).This time delay, to the first order in eccentricity, is given by where  is the argument of periapsis (Falanga et al. 2015).In case of LMC X-4, the difference in the best-fit value of T 0 from T  and T /2 measurements is 0.0041 (14) d, which implies a value of ∼0.009 (3) for  cos .
The tidal interaction between the binary components of HMXBs is known to circularize their orbit.As a result, the systems amongst these with small orbital periods (of up to 4 days) have very low eccentricity.Updated orbital solution for HMXBs is important because a good estimate of eccentricity of the orbit has a potential to constrain the age of the binary system and also serve as indirect detection of gravitational wave emission from the system.Likewise, a good estimate of the periastron angle and advancement of the periastron angle with time (e.g., 4U 0115+63: Raichur & Paul 2010b) is a clue to understanding the stellar structure models of these massive binary systems.Our estimate of  cos  in the case of LMC X-4 can serve as an important ingredient in determining these numbers.

Figure 1 .Figure 2 .
Figure1.The light curve of LMC X-4 obtained from all the observations listed in Table1.The respective X-ray mission, observation ID and energy range are mentioned in each panel.The blue colored segments were used for the pulse timing measurements.The segments shown in red color were not used.

Figure 3 .
Figure 3.The eclipse profile of the light curve of LMC X-4 obtained from the first segment of the long term MAXI-GSC (top) and Swift-BAT (below) light curves.The red solid line represents the best-fit ramp model.

Figure 5 .
Figure 5.The residual in units of sigma (=(data-model)/error) of the O-C history of LMC X-4 for LQ and LQC models, depicting the individual contribution to the net  2 .The labelling with black and magenta points is similar to that in Figure 4.

Figure 6 .
Figure 6.Eclipse profile of the XMM-Newton observation (Obs.Id. 0142800101) used in the current work.The vertical lines correspond to the expected mid-eclipse time from various techniques.Red:  /2 , Blue: T  from LQ model, Green: T  from LQC, and Magenta: T  from Molkov et al. (2015).The inset figure shows the zoomed-in section near the centre of the eclipse.

Table 1 .
Log of X-ray observations of LMC X-4 used in this work. Net exposure after removing flaring, eclipse and orbital gaps, if any.
Also reported by Figure 4. O-C history of LMC X-4 as a function of the orbital cycle, relative to MJD 53013.5878.The data points related to   and  /2 have been marked with black and magenta points, respectively, in the top panel.The best-fit LQ and LQC models are plotted in red and blue colours, respectively.Below two panels show the residuals with LQ and LQC models, respectively.

Table 4 .
Updated orbital parameters of LMC X-4.  orb =  orb /  orb is evolution time scale of orbital period.   orb =  orb /  orb is evolution time scale of orbital period derivative. Error in T 0 , P orb and  orb have been artificially increased such that respective  2