  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 $T_{\pi/2}$ epoch
history of LMC X-4 by about 4600 binary orbits spanning about 18 years. We also
report mid-eclipse time measurements ($T_{ecl}$) using data obtained from
wide-field X-ray monitors of MAXI-GSC and Swift-BAT. Combining the new
$T_{\pi/2}$ and $T_{ecl}$ 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 ($P_{\rm orb}/\dot{P}_{\rm
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 ($\dot{P}_{orb}/\ddot{P}_{orb}$) 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.
