  Background: It has been proposed that the azimuthal distributions of heavy
flavor quark-antiquark 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\sigma/d\phi$,
including transverse momentum broadening, employing $<k_T^2>$ and fragmentation
in exclusive $Q \bar Q$ pair production. While these studies were done for
$p+p$, $p + \bar 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 $p_T$ 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_T^2>$, $p_T$ cut, and rapidity are studied at $\sqrt{s} = 7$ TeV. Finally,
calculations are compared to $Q \bar Q$ data in elementary $p+p$ and $p + \bar
p$ collisions at $\sqrt{s} = 7$ TeV and 1.96 TeV as well as to the nuclear
modification factor $R_{p {\rm Pb}}(p_T)$ in $p+$Pb collisions at
$\sqrt{s_{NN}} = 5.02$ TeV measured by ALICE. Results: The low $p_T$ ($p_T <
10$ GeV) azimuthal distributions are very sensitive to the $k_T$ broadening and
rather insensitive to the fragmentation function. The NLO contributions can
result in an enhancement at $\phi \sim 0$ absent any other effects. Agreement
with the data was found to be good. Conclusions: The NLO calculations, assuming
collinear factorization and introducing $k_T$ broadening, result in significant
modifications of the azimuthal distribution at low $p_T$ which must be taken
into account in calculations of these distributions in heavy-ion collisions.
