Time-domain elastic least-squares reverse time migration (LSRTM) is formulated as a linearized elastic full-waveform inversion problem. The elastic Born approximation and elastic reverse time migration (RTM) operators are derived from the time-domain continuous adjoint-state method. Our approach defines P- and S-wave impedance perturbations as unknown elastic images. Our algorithm is obtained using continuous functional analysis in which the problem is discretized at the final stage (optimize-before-discretize approach). The discretized numerical versions of the elastic Born operator and its adjoint (elastic RTM operator) pass the dot-product test. The conjugate gradient least-squares method is used to solve the least-squares migration quadratic optimization problem. In other words, the Hessian operator for elastic LSRTM is implicitly inverted via a matrix-free algorithm that only requires the action of forward and adjoint operators on vectors. The diagonal of the pseudo-Hessian operator is used to design a preconditioning operator to accelerate the convergence of the elastic LSRTM. The elastic LSRTM provides higher resolution images with fewer artifacts and a superior balance of amplitudes when compared with elastic RTM. More important, elastic LSRTM can mitigate crosstalk between the P- and S-wave impedance perturbations given that the off-diagonal elements of the Hessian are attenuated via the inversion.