  The study of azimuthal correlations in particle production at forward
rapidities in proton-nucleus collisions provides direct information about high
gluon density effects, like gluon saturation, in the nuclear wavefunction. In
the kinematical conditions for proton-lead collisions at the LHC, the forward
di-hadron production is dominated by partonic processes in which a gluon from
the proton splits into a pair of gluons, while undergoing multiple scattering
off the dense gluon system in the nucleus. We compute the corresponding
cross-section using the Colour Glass Condensate effective theory, which enables
us to include the effects of multiple scattering and gluon saturation in the
leading logarithmic approximation at high energy. This opens the way towards
systematic studies of angular correlations in two-gluon production, similar to
previous studies for quark-gluon production in the literature. We consider in
more detail two special kinematical limits: the "back-to-back correlation
limit", where the transverse momenta of the produced gluons are much larger
than the nuclear saturation momentum, and the "double parton scattering limit",
where the two gluons are produced by a nearly collinear splitting occurring
prior to the collision. We argue that saturation effects remain important even
for relatively high transverse momenta (i.e. for nearly back-to-back
configurations), leading to geometric scaling in the azimuthal distribution.
