Background: It has been proposed that the azimuthal distributions of heavy flavor quarkantiquark pairs may be modified in the medium of a heavy-ion collision. Purpose: This work tests this proposition through next-to-leading order (NLO) calculations of the azimuthal distribution, dσ/dφ, including transverse momentum broadening, employing k 2 T and fragmentation in exclusive QQ pair production. While these studies were done for p+p, p+p and p+Pb collisions, understanding azimuthal angle correlations between heavy quarks in these smaller, colder systems is important for their interpretation in heavy-ion collisions. Methods: First, single inclusive pT distributions calculated with the exclusive HVQMNR code are compared to those calculated in the fixed-order next-to-leading logarithm approach. Next the azimuthal distributions are calculated and sensitivities to k 2 T , pT cut, and rapidity are studied at √ s = 7 TeV. Finally, calculations are compared to QQ data in elementary p + p and p + p collisions at √ s = 7 TeV and 1.96 TeV as well as to the nuclear modification factor R pPb (pT ) in p+Pb collisions at √ sNN = 5.02 TeV measured by ALICE. Results: The low pT (pT < 10 GeV) azimuthal distributions are very sensitive to the kT broadening and rather insensitive to the fragmentation function. The NLO contributions can result in an enhancement at φ ∼ 0 absent any other effects. Agreement with the data was found to be good. Conclusions: The NLO calculations, assuming collinear factorization and introducing kT broadening, result in significant modifications of the azimuthal distribution at low pT which must be taken into account in calculations of these distributions in heavy-ion collisions.