For $J/psi$ pair production at hadron colliders, we present the full next-to-leading order (NLO) calculations with the color-singlet channel in nonrelativistic QCD. We find that the NLO result can reasonably well describe the LHCb measured cross sect
ion, but exhibits very different behaviors from the CMS data in the transverse momentum distribution and mass distribution of $J/psi$ pair. Moreover, by adding contributions of gluon fragmentation and quark fragmentation, which occur at even higher order in $alpha_s$, it is still unable to reduce the big differences. In particular, the observed flat distribution in the large invariant mass region is hard to explain. New processes or mechanisms are needed to understand the CMS data for $J/psi$ pair production.
