  We develop a set of coupled Boltzmann equations to describe the dynamical
evolution of heavy quarks and quarkonia inside the quark-gluon plasma. The
quarkonium dissociation and recombination terms are calculated from pNRQCD.
Their interplay drives the system to a detailed balance. The heavy quark energy
loss term is necessary for the system to reach kinematic thermalization. By
coupling the transport equations with initial particles' momenta generated by
PYTHIA and hydrodynamic medium evolutions, we can describe the $R_{AA}$ of
$\Upsilon$ family at both RHIC and LHC energies. The transverse momentum
azimuthal anisotropy of $\Upsilon$(1S) in $2.76$ TeV peripheral Pb-Pb
collisions is also studied.
