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 these reveal the couplings of correlated motions, which exist over a wide range of time scales. The correlated motions detected by the three-time correlation function are 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.