The STAR Collaboration reports measurements of back-to-back azimuthal correlations of di-$\pi^0$s produced at forward pseudorapidities ($2.6<\eta<4.0$) in $p$+$p$, $p+$Al, and $p+$Au collisions at a center-of-mass energy of 200 GeV. We observe a clear suppression of the correlated yields of back-to-back $\pi^0$ pairs in $p+$Al and $p+$Au collisions compared to the $p$+$p$ data. The observed suppression of back-to-back pairs as a function of transverse momentum suggests nonlinear gluon dynamics arising at high parton densities. The larger suppression found in $p+$Au relative to $p+$Al collisions exhibits a dependence of the saturation scale, $Q_s^2$, on the mass number, $A$. A linear scaling of the suppression with $A^{1/3}$ is observed with a slope of $-0.09$$\pm$$0.01$.
The correlation functions (corrected for nonuniform detector efficiency in $\phi$; not corrected for the absolute detection efficiency) vs. azimuthal angle difference between forward ($2.6<\eta<4.0$) $\pi^{0}$s in $p$+$p$ collisions at $\sqrt{s_{\mathrm{_{NN}}}}=200$ GeV at low $p_{T}$ ($p^{trig}_{T}$=2-2.5 GeV/c, $p^{asso}_{T}$=1-1.5 GeV/c)
The correlation functions (corrected for nonuniform detector efficiency in $\phi$; not corrected for the absolute detection efficiency) vs. azimuthal angle difference between forward ($2.6<\eta<4.0$) $\pi^{0}$s in $p+$Al collisions at $\sqrt{s_{\mathrm{_{NN}}}}=200$ GeV at low $p_{T}$ ($p^{trig}_{T}$=2-2.5 GeV/c, $p^{asso}_{T}$=1-1.5 GeV/c)
The correlation functions (corrected for nonuniform detector efficiency in $\phi$; not corrected for the absolute detection efficiency) vs. azimuthal angle difference between forward ($2.6<\eta<4.0$) $\pi^{0}$s in $p+$Au collisions at $\sqrt{s_{\mathrm{_{NN}}}}=200$ GeV at low $p_{T}$ ($p^{trig}_{T}$=2-2.5 GeV/c, $p^{asso}_{T}$=1-1.5 GeV/c)