  We study the azimuthal $\cos \varphi$ and $\cos 2\varphi$ asymmetries in
heavy-quark pair leptoproduction, $lN\rightarrow l^{\prime}Q\bar{Q}X$, as
probes of linearly polarized gluons inside unpolarized proton, where the
azimuth $\varphi$ is the angle between the lepton scattering plane
$(l,l^{\prime})$ and the heavy quark production plane $(N,Q)$. First, we
determine the maximal values for the $\cos \varphi$ and $\cos 2\varphi$
asymmetries allowed by the photon-gluon fusion with unpolarized gluons; these
predictions are large, $(\sqrt{3}-1)/2$ and $1/3$, respectively. Then we
calculate the contribution of the transverse-momentum dependent gluonic
counterpart of the Boer-Mulders function, $h_{1}^{\perp g}$, describing the
linear polarization of gluons inside unpolarized proton. Our analysis shows
that the maximum values of the azimuthal distributions depend strongly on the
gluon polarization; they vary from 0 to 1 depending on $h_{1}^{\perp g}$. We
conclude that the azimuthal $\cos \varphi$ and $\cos 2\varphi$ asymmetries in
heavy-quark pair leptoproduction are predicted to be large and very sensitive
to the contribution of linearly polarized gluons. For this reason, future
measurements of the azimuthal distributions in charm and bottom production at
the proposed EIC and LHeC colliders seem to be very promising for determination
of the linear polarization of gluons inside unpolarized proton.
