  Production of large transverse momentum $\rho^{0}$ meson in high-energy
nuclear collisions is investigated for the first time at the next-leading-order
in the QCD improved parton model. The $\rho^0$ fragmentation functions (FFs) in
vacuum at any scale $Q$ are obtained, by evolving a newly developed initial
parametrization of $\rho^0$ FFs at a scale $\rm Q_{0}^2=1.5\ GeV^2$ from a
broken SU(3) model through NLO DGLAP equations. The numerical simulations of
$p_{\rm T}$ spectra of $\rho^{0}$ meson in the elementary $\rm p+p$ collisions
at NLO give a decent description of STAR $\rm p+p$ data. In $\rm A+A$ reactions
the jet quenching effect is taken into account with the higher-twist approach
by the medium-modified parton FFs due to gluon radiation in the quark-gluon
plasma, whose space-time evolution is described by a (3+1D) hydrodynamical
model. The nuclear modification factors for $\rho^{0}$ meson and its double
ratio with $\pi^\pm$ nuclear modification in central $\rm Au+Au$ collisions at
the RHIC are calculated and found to be in good agreement with STAR
measurement. Predictions of $\rho^{0}$ nuclear modification and the yield ratio
$\rho^0/\pi^0$ in central Pb+Pb at the LHC are also presented. It is shown that
the ratio $\rho^0/\pi^0$ in central Pb+Pb will approach to that in p+p
reactions when $p_{\rm T}>12$~GeV.
