  A multi-time extension of a density correlation function is introduced to
reveal temporal information about dynamical heterogeneity in glass-forming
liquids. We utilize a multi-time correlation function that is analogous to the
higher-order response function analyzed in multidimensional nonlinear
spectroscopy. Here, we provide comprehensive numerical results of the
four-point, three-time density correlation function from longtime trajectories
generated by molecular dynamics simulations of glass-forming binary soft-sphere
mixtures. We confirm that the two-dimensional representations in both time and
frequency domains are sensitive to the dynamical heterogeneity and that it
reveals the couplings of correlated motions, which exist over a wide range of
time scales. The correlated motions detected by the three-time correlation
function is divided into mobile and immobile contributions that are determined
from the particle displacement during the first time interval. We show that the
peak positions of the correlations are in accord with the information on the
non-Gaussian parameters of the van-Hove self correlation function. Furthermore,
it is demonstrated that the progressive changes in the second time interval in
the three-time correlation function enable us to analyze how correlations in
dynamics evolve in time. From this analysis, we evaluated the lifetime of the
dynamical heterogeneity and its temperature dependence systematically. Our
results show that the lifetime of the dynamical heterogeneity becomes much
slower than the $\alpha$-relaxation time that is determined from the two-point
density correlation function when the system is highly supercooled.
