We present the theory and numerical results for least-squares reverse time migration (LSRTM) of phase-encoded supergathers, where each supergather is the superposition of phased-encoded shots. Three type of encoding functions are used in this study: random time shift, random source polarity and random source location selected from a pre-designed table. Numerical tests for the 3D SEG/EAGE Overthrust model show that multi-source LSRTM can suppress migration artifacts in the migration image and remove most of the crosstalk noise from multi-source data. Empirical results suggest that multi-source LSRTM can provide a noticeable increase in computational efficiency compared to standard RTM, when the CSGs in a supergather are modeled and migrated together with a finite-difference simulator. If the phase-encoding functions are dynamically changed after each iteration of LSRTM, the best images are obtained. The potential drawback is that the final results are very sensitive to the accuracy of the starting model.