  Revealing the large-scale structure from the 21cm intensity mapping surveys
is only possible after the foreground cleaning. However, most current cleaning
techniques relying on the smoothness of the foreground spectrum lead to a
severe side effect of removing the large-scale structure signal along the line
of sight. On the other hand, the clustering fossil, a coherent variation of the
small-scale clustering over large scales, allows us to recover the
long-wavelength density modes from the off-diagonal correlation between
short-wavelength modes. In this paper, we study the requirements for an
unbiased and optimal clustering-fossil estimator and show that (A) the
estimator is unbiased only when using an accurate bispectrum model for the
long-short-short mode coupling and (B) including the connected four-point
correlation functions is essential for characterizing the noise power spectrum
of the estimated long mode. The clustering fossil estimator based upon the
leading-order bispectrum yields an unbiased estimation of the long-wavelength
($k\lesssim 0.01~[h/{\rm Mpc}]$) modes with the cross-correlation coefficient
of $0.7$ at redshifts $z=0$ to $3$.
