We present measurements of elliptic and triangular azimuthal anisotropy of charged particles detected at forward rapidity $1<|\eta|<3$ in Au$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV, as a function of centrality. The multiparticle cumulant technique is used to obtain the elliptic flow coefficients $v_2\{2\}$, $v_2\{4\}$, $v_2\{6\}$, and $v_2\{8\}$, and triangular flow coefficients $v_3\{2\}$ and $v_3\{4\}$. Using the small-variance limit, we estimate the mean and variance of the event-by-event $v_2$ distribution from $v_2\{2\}$ and $v_2\{4\}$. In a complementary analysis, we also use a folding procedure to study the distributions of $v_2$ and $v_3$ directly, extracting both the mean and variance. Implications for initial geometrical fluctuations and their translation into the final state momentum distributions are discussed.

Measurements of anisotropic flow Fourier coefficients ($v_n$) for inclusive charged particles and identified hadrons $\pi^{\pm}$, $K^{\pm}$, $p$, and $\bar{p}$ produced at midrapidity in Cu+Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV are presented. The data were collected in 2012 by the PHENIX experiment at the Relativistic Heavy Ion Collider (RHIC). The particle azimuthal distributions with respect to different order symmetry planes $\Psi_n$, for $n$~=~1, 2, and 3 are studied as a function of transverse momentum $p_T$ over a broad range of collisions centralities. Mass ordering, as expected from hydrodynamic flow, is observed for all three harmonics. The charged-particle results are compared to hydrodynamical and transport model calculations. We also compare these Cu$+$Au results with those in Cu$+$Cu and Au$+$Au collisions at the same $\sqrt{s_{_{NN}}}$, and find that the $v_2$ and $v_3$, as a function of transverse momentum, follow a common scaling with $1/(\varepsilon_n N_{\rm part}^{1/3})$.

The PHENIX Collaboration has measured the ratio of the yields of $\psi(2S)$ to $\psi(1S)$ mesons produced in $p$$+$$p$, $p$$+$Al, $p$$+$Au, and $^{3}$He$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV over the forward and backward rapidity intervals $1.2<|y|<2.2$. We find that the ratio in $p$$+$$p$ collisions is consistent with measurements at other collision energies. In collisions with nuclei, we find that in the forward ($p$-going or $^{3}$He-going) direction, the relative yield of $\psi(2S)$ mesons to $\psi(1S)$ mesons is consistent with the value measured in \pp collisions. However, in the backward (nucleus-going) direction, the $\psi(2S)$ is preferentially suppressed by a factor of $\sim$2. This suppression is attributed in some models to breakup of the weakly-bound $\psi(2S)$ through final state interactions with comoving particles, which have a higher density in the nucleus-going direction. These breakup effects may compete with color screening in a deconfined quark-gluon plasma to produce sequential suppression of excited quarkonia states.

We present the first measurements of long-range angular correlations and the transverse momentum dependence of elliptic flow $v_2$ in high-multiplicity $p$$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV. A comparison of these results with previous measurements in high-multiplicity $d$$+$Au and $^3{\rm He}$$+$Au collisions demonstrates a relation between $v_2$ and the initial collision eccentricity $\varepsilon_2$, suggesting that the observed momentum-space azimuthal anisotropies in these small systems have a collective origin and reflect the initial geometry. Good agreement is observed between the measured $v_2$ and hydrodynamic calculations for all systems, and an argument disfavoring theoretical explanations based on momentum-space domain correlations is presented. The set of measurements presented here allows us to leverage the distinct intrinsic geometry of each of these systems to distinguish between different theoretical descriptions of the long-range correlations observed in small collision systems.

The PHENIX experiment at the Relativistic Heavy Ion Collider has measured the differential cross section of $\phi$(1020) meson production at forward rapidity in $p$$+$$p$ collisions at $\sqrt{s}=$510 GeV via the dimuon decay channel. The integrated cross section in the rapidity and $p_T$ ranges $1.2<|y|<2.2$ and $2<p_T<7$ GeV/$c$ is $\sigma_\phi=2.79 \pm 0.20\,{\rm (stat)} \pm 0.17\,{\rm (syst)} \pm 0.34\, {\rm (norm)} \times 10^{-2}$~mb. The energy dependence of $\sigma_\phi$ ($1.2<|y|<2.2$; $2<p_T<5$ GeV/$c$) is studied using the PHENIX measurements at $\sqrt{s}=$200 and 510 GeV and the Large-Hadron-Collider measurements at $\sqrt{s}=$2.76 and 7 TeV. The experimental results are compared to various event generator predictions (pythia6, pythia8, phojet, ampt, epos3, and epos-lhc).

PHENIX measurements are presented for the cross section and double-helicity asymmetry ($A_{LL}$) in inclusive $\pi^0$ production at midrapidity from $p$$+$$p$ collisions at $\sqrt{s}=510$~GeV from data taken in 2012 and 2013 at the Relativistic Heavy Ion Collider. The next-to-leading-order perturbative-quantum-chromodynamics theory calculation is in excellent agreement with the presented cross section results. The calculation utilized parton-to-pion fragmentation functions from the recent DSS14 global analysis, which prefer a smaller gluon-to-pion fragmentation function. The $\pi^{0}A_{LL}$ results follow an increasingly positive asymmetry trend with $p_T$ and $\sqrt{s}$ with respect to the predictions and are in excellent agreement with the latest global analysis results. This analysis incorporated earlier results on $\pi^0$ and jet $A_{LL}$, and suggested a positive contribution of gluon polarization to the spin of the proton $\Delta G$ for the gluon momentum fraction range $x>0.05$. The data presented here extend to a currently unexplored region, down to $x\sim0.01$, and thus provide additional constraints on the value of $\Delta G$. The results confirm the evidence for nonzero $\Delta G$ using a different production channel in a complementary kinematic region.

PHENIX reports differential cross sections of $\mu\mu$ pairs from semileptonic heavy-flavor decays and the Drell-Yan production mechanism measured in $p$$+$$p$ collisions at $\sqrt{s}=200$ GeV at forward and backward rapidity ($1.2<|\eta|<2.2$). The $\mu\mu$ pairs from $c\bar{c}$, $b\bar{b}$, and Drell-Yan are separated using a simultaneous fit to unlike- and like-sign muon pair spectra in mass and $p_T$. The azimuthal opening angle correlation between the muons from $c\bar{c}$ and $b\bar{b}$ decays and the pair-$p_T$ distributions are compared to distributions generated using {\sc pythia} and {\sc powheg} models, which both include next-to-leading order processes. The measured distributions for pairs from $c\bar{c}$ are consistent with {\sc pythia} calculations. The $c\bar{c}$ data presents narrower azimuthal correlations and softer $p_T$ distributions compared to distributions generated from {\sc powheg}. The $b\bar{b}$ data are well described by both models. The extrapolated total cross section for bottom production is $3.75{\pm}0.24({\rm stat}){\pm}^{0.35}_{0.50}({\rm syst}){\pm}0.45({\rm global})$[$\mu$b], which is consistent with previous measurements at the Relativistic Heavy Ion Collider in the same system at the same collision energy, and is approximately a factor of two higher than the central value calculated with theoretical models. The measured Drell-Yan cross section is in good agreement with next-to-leading-order quantum-chromodynamics calculations.

The invariant yields for $J/\psi$ production at forward rapidity $(1.2<|y|<2.2)$ in U$+$U collisions at $\sqrt{s_{_{NN}}}$=193 GeV have been measured as a function of collision centrality. The invariant yields and nuclear-modification factor $R_{AA}$ are presented and compared with those from Au$+$Au collisions in the same rapidity range. Additionally, the direct ratio of the invariant yields from U$+$U and Au$+$Au collisions within the same centrality class is presented, and used to investigate the role of $c\bar{c}$ coalescence. Two different parameterizations of the deformed Woods-Saxon distribution were used in Glauber calculations to determine the values of the number of nucleon-nucleon collisions in each centrality class, $N_{\rm coll}$, and these were found to give significantly different $N_{\rm coll}$ values. Results using $N_{\rm coll}$ values from both deformed Woods-Saxon distributions are presented. The measured ratios show that the $J/\psi$ suppression, relative to binary collision scaling, is similar in U$+$U and Au$+$Au for peripheral and midcentral collisions, but that $J/\psi$ show less suppression for the most central U$+$U collisions. The results are consistent with a picture in which, for central collisions, increase in the $J/\psi$ yield due to $c\bar{c}$ coalescence becomes more important than the decrease in yield due to increased energy density. For midcentral collisions, the conclusions about the balance between $c\bar{c}$ coalescence and suppression depend on which deformed Woods-Saxon distribution is used to determine $N_{\rm coll}$.

We present a detailed measurement of charged two-pion correlation functions in 0%-30% centrality $\sqrt{s_{_{NN}}}=200$ GeV Au$+$Au collisions by the PHENIX experiment at the Relativistic Heavy Ion Collider. The data are well described by Bose-Einstein correlation functions stemming from L\'evy-stable source distributions. Using a fine transverse momentum binning, we extract the correlation strength parameter $\lambda$, the L\'evy index of stability $\alpha$ and the L\'evy length scale parameter $R$ as a function of average transverse mass of the pair $m_T$. We find that the positively and the negatively charged pion pairs yield consistent results, and their correlation functions are represented, within uncertainties, by the same L\'evy-stable source functions. The $\lambda(m_T)$ measurements indicate a decrease of the strength of the correlations at low $m_T$. The L\'evy length scale parameter $R(m_T)$ decreases with increasing $m_T$, following a hydrodynamically predicted type of scaling behavior. The values of the L\'evy index of stability $\alpha$ are found to be significantly lower than the Gaussian case of $\alpha=2$, but also significantly larger than the conjectured value that may characterize the critical point of a second-order quark-hadron phase transition.

We present a systematic study of charged pion and kaon interferometry in Au$+$Au collisions at $\sqrt{s_{_{NN}}}$=200 GeV. The kaon mean source radii are found to be larger than pion radii in the outward and longitudinal directions for the same transverse mass; this difference increases for more central collisions. The azimuthal-angle dependence of the radii was measured with respect to the second-order event plane and similar oscillations of the source radii were found for pions and kaons. Hydrodynamic models qualitatively describe the similar oscillations of the mean source radii for pions and kaons, but they do not fully describe the transverse-mass dependence of the oscillations.

We have measured direct photons for $p_T<5~$GeV/$c$ in minimum bias and 0\%--40\% most central events at midrapidity for Cu$+$Cu collisions at $\sqrt{s_{_{NN}}}=200$ GeV. The $e^{+}e^{-}$ contribution from quasi-real direct virtual photons has been determined as an excess over the known hadronic contributions in the $e^{+}e^{-}$ mass distribution. A clear enhancement of photons over the binary scaled $p$$+$$p$ fit is observed for $p_T<4$ GeV/$c$ in Cu$+$Cu data. The $p_T$ spectra are consistent with the Au$+$Au data covering a similar number of participants. The inverse slopes of the exponential fits to the excess after subtraction of the $p$$+$$p$ baseline are 285$\pm$53(stat)$\pm$57(syst)~MeV/$c$ and 333$\pm$72(stat)$\pm$45(syst)~MeV/$c$ for minimum bias and 0\%--40\% most central events, respectively. The rapidity density, $dN/dy$, of photons demonstrates the same power law as a function of $dN_{\rm ch}/d\eta$ observed in Au$+$Au at the same collision energy.

Dihadron and isolated direct photon-hadron angular correlations are measured in $p$$+$$p$ collisions at $\sqrt{s}=510$ GeV. Correlations of charged hadrons of $0.7<p_T<10$ GeV/$c$ with $\pi^0$ mesons of $4<p_T<15$ GeV/$c$ or isolated direct photons of $7<p_T<15$ GeV/$c$ are used to study nonperturbative effects generated by initial-state partonic transverse momentum and final-state transverse momentum from fragmentation. The nonperturbative behavior is characterized by measuring the out-of-plane transverse momentum component $p_{\rm out}$ perpendicular to the axis of the trigger particle, which is the high-$p_T$ direct photon or $\pi^0$. Nonperturbative evolution effects are extracted from Gaussian fits to the away-side inclusive-charged-hadron yields for different trigger-particle transverse momenta ($p_T^{\rm trig}$). The Gaussian widths and root mean square of $p_{\rm out}$ are reported as a function of the interaction hard scale $p_T^{\rm trig}$ to investigate possible transverse-momentum-dependent evolution differences between the $\pi^0$-h$^\pm$ and direct photon-h$^\pm$ correlations and factorization breaking effects. The widths are found to decrease with $p_T^{\rm trig}$, which indicates that the Collins-Soper-Sterman soft factor is not driving the evolution with the hard scale in nearly back-to-back dihadron and direct photon-hadron production in $p$$+$$p$ collisions. This behavior is in contrast to Drell-Yan and semi-inclusive deep-inelastic scattering measurements.

We report the first measurement of the full angular distribution for inclusive $J/\psi\rightarrow\mu^{+}\mu^{-}$ decays in $p$$+$$p$ collisions at $\sqrt{s}=510$ GeV. The measurements are made for $J/\psi$ transverse momentum $2<p_{T}<10$ GeV/$c$ and rapidity $1.2<y<2.2$ in the Helicity, Collins-Soper, and Gottfried-Jackson reference frames. In all frames the polar coefficient $\lambda_{\theta}$ is strongly negative at low $p_{T}$ and becomes close to zero at high $p_{T}$, while the azimuthal coefficient $\lambda_{\phi}$ is close to zero at low $p_{T}$, and becomes slightly negative at higher $p_{T}$. The frame-independent coefficient $\tilde{\lambda}$ is strongly negative at all $p_{T}$ in all frames. The data are compared to the theoretical predictions provided by nonrelativistic quantum chromodynamics models.

Dihadron and isolated direct photon-hadron angular correlations are measured in $p$$+$$p$ collisions at $\sqrt{s}=200$ GeV. The correlations are sensitive to nonperturbative initial-state and final-state transverse momentum $k_T$ and $j_T$ in the azimuthal nearly back-to-back region $\Delta\phi\sim\pi$. In this region, transverse-momentum-dependent evolution can be studied when several different hard scales are measured. To have sensitivity to small transverse momentum scales, nonperturbative momentum widths of $p_{\rm out}$, the out-of-plane transverse momentum component perpendicular to the trigger particle, are measured. These widths are used to investigate possible effects from transverse-momentum-dependent factorization breaking. When accounting for the longitudinal momentum fraction of the away-side hadron with respect to the near-side trigger particle, the widths are found to increase with the hard scale; this is qualitatively similar to the observed behavior in Drell-Yan and semi-inclusive deep-inelastic scattering interactions. The momentum widths are also studied as a function of center-of-mass energy by comparing to previous measurements at $\sqrt{s}=510$ GeV. The nonperturbative jet widths also appear to increase with $\sqrt{s}$ at a similar $x_T$, which is qualitatively consistent to similar measurements in Drell-Yan interactions. To quantify the magnitude of any transverse-momentum-dependent factorization breaking effects, calculations will need to be performed to compare to these measurements.

The PHENIX experiment at the Relativistic Heavy Ion Collider has measured 2nd and 3rd order Fourier coefficients of the azimuthal distributions of direct photons emitted at midrapidity in Au$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV for various collision centralities. Combining two different analysis techniques, results were obtained in the transverse momentum range of $0.4<p_{T}<4.0$ GeV/$c$. At low $p_T$ the second-order coefficients, $v_2$, are similar to the ones observed in hadrons. Third order coefficients, $v_3$, are nonzero and almost independent of centrality. These new results on $v_2$ and $v_3$, combined with previously published results on yields, are compared to model calculations that provide yields and asymmetries in the same framework. Those models are challenged to explain simultaneously the observed large yield and large azimuthal anisotropies.

We present measurements of $e^+e^-$ production at midrapidity in Au$+$Au collisions at $\sqrt{s_{_{NN}}}$ = 200 GeV. The invariant yield is studied within the PHENIX detector acceptance over a wide range of mass ($m_{ee} <$ 5 GeV/$c^2$) and pair transverse momentum ($p_T$ $<$ 5 GeV/$c$), for minimum bias and for five centrality classes. The \ee yield is compared to the expectations from known sources. In the low-mass region ($m_{ee}=0.30$--0.76 GeV/$c^2$) there is an enhancement that increases with centrality and is distributed over the entire pair \pt range measured. It is significantly smaller than previously reported by the PHENIX experiment and amounts to $2.3\pm0.4({\rm stat})\pm0.4({\rm syst})\pm0.2^{\rm model}$ or to $1.7\pm0.3({\rm stat})\pm0.3({\rm syst})\pm0.2^{\rm model}$ for minimum bias collisions when the open-heavy-flavor contribution is calculated with {\sc pythia} or {\sc mc@nlo}, respectively. The inclusive mass and $p_T$ distributions as well as the centrality dependence are well reproduced by model calculations where the enhancement mainly originates from the melting of the $\rho$ meson resonance as the system approaches chiral symmetry restoration. In the intermediate-mass region ($m_{ee}$ = 1.2--2.8 GeV/$c^2$), the data hint at a significant contribution in addition to the yield from the semileptonic decays of heavy-flavor mesons.

The cross section and transverse single-spin asymmetries of $\mu^{-}$ and $\mu^{+}$ from open heavy-flavor decays in polarized $p$+$p$ collisions at $\sqrt{s}=200$ GeV were measured by the PHENIX experiment during 2012 at the Relativistic Heavy Ion Collider. Because heavy-flavor production is dominated by gluon-gluon interactions at $\sqrt{s}=200$ GeV, these measurements offer a unique opportunity to obtain information on the trigluon correlation functions. The measurements are performed at forward and backward rapidity ($1.4<|y|<2.0$) over the transverse momentum range of $1.25<p_T<7$ GeV/$c$ for the cross section and $1.25<p_T<5$ GeV/$c$ for the asymmetry measurements. The obtained cross section is compared to a fixed-order-plus-next-to-leading-log perturbative-quantum-chromodynamics calculation. The asymmetry results are consistent with zero within uncertainties, and a model calculation based on twist-3 three-gluon correlations agrees with the data.

Measurements of the fractional momentum loss ($S_{\rm loss}\equiv{\delta}p_T/p_T$) of high-transverse-momentum-identified hadrons in heavy ion collisions are presented. Using $\pi^0$ in Au$+$Au and Cu$+$Cu collisions at $\sqrt{s_{_{NN}}}=62.4$ and 200 GeV measured by the PHENIX experiment at the Relativistic Heavy Ion Collider and and charged hadrons in Pb$+$Pb collisions measured by the ALICE experiment at the Large Hadron Collider, we studied the scaling properties of $S_{\rm loss}$ as a function of a number of variables: the number of participants, $N_{\rm part}$, the number of quark participants, $N_{\rm qp}$, the charged-particle density, $dN_{\rm ch}/d\eta$, and the Bjorken energy density times the equilibration time, $\varepsilon_{\rm Bj}\tau_{0}$. We find that the $p_T$ where $S_{\rm loss}$ has its maximum, varies both with centrality and collision energy. Above the maximum, $S_{\rm loss}$ tends to follow a power-law function with all four scaling variables. The data at $\sqrt{s_{_{NN}}}$=200 GeV and 2.76 TeV, for sufficiently high particle densities, have a common scaling of $S_{\rm loss}$ with $dN_{\rm ch}/d\eta$ and $\varepsilon_{\rm Bj}\tau_{0}$, lending insight on the physics of parton energy loss.

The PHENIX experiment at the Relativistic Heavy Ion Collider (RHIC) has measured $\phi$ meson production and its nuclear modification in asymmetric Cu$+$Au heavy-ion collisions at $\sqrt{s_{NN}}=200$ GeV at both forward Cu-going direction ($1.2<y<2.2$) and backward Au-going direction ($-2.2<y<-1.2$), rapidities. The measurements are performed via the dimuon decay channel and reported as a function of the number of participating nucleons, rapidity, and transverse momentum. In the most central events, 0\%--20\% centrality, the $\phi$ meson yield integrated over $1<p_T<5$ GeV/$c$ prefers a smaller value, which means a larger nuclear modification, in the Cu-going direction compared to the Au-going direction. Additionally, the nuclear-modification factor in Cu$+$Au collisions averaged over all centrality is measured to be similar to the previous PHENIX result in $d$$+$Au collisions for these rapidities.

We report the transverse single-spin asymmetries of $J/\psi$ production at forward and backward rapidity, $1.2<|y|<2.2$, as a function of $J/\psi$ transverse momentum ($p_T$) and Feynman-$x$ ($x_F$). The data analyzed were recorded by the PHENIX experiment at the Relativistic Heavy Ion Collider in 2015 from $p$$+$$p$, $p$$+$Al, and $p$$+$Au collisions with transversely polarized proton beams at $\sqrt{s_{_{NN}}}=200$ GeV. At this collision energy, single-spin asymmetries for heavy-flavor particle production of $p$$+$$p$ collisions provide access to the spin-dependent gluon distribution and higher-twist correlation functions inside the nucleon, such as the gluon Qiu-Sterman and trigluon correlation functions. Proton+nucleus collisions offer an excellent opportunity to study nuclear effects on the correlation functions. The data indicate negative asymmetries at the two-standard-deviation level in the $p$$+$Au data for $p_T<2$ GeV/$c$ at both forward and backward rapidity, while in $p$$+$$p$ and $p$$+$Al collisions the asymmetries are consistent with zero within the range of experimental uncertainties.

We report a measurement of $e^+e^-$ pairs from semileptonic heavy-flavor decays in $p$+$p$ collisions at $\sqrt{s_{NN}}=200$~GeV. The $e^+e^-$ pair yield from $b\bar{b}$ and $c\bar{c}$ is separated by exploiting a double differential fit done simultaneously in dielectron invariant mass and $p_T$. We used three different event generators, {\sc pythia}, {\sc mc@nlo}, and {\sc powheg}, to simulate the $e^+e^-$ spectra from $c\bar{c}$ and $b\bar{b}$ production. The data can be well described by all three generators within the detector acceptance. However, when using the generators to extrapolate to $4\pi$, significant differences are observed for the total cross section. These difference are less pronounced for $b\bar{b}$ than for $c\bar{c}$. The same model dependence was observed in already published $d$+$A$ data. The $p$+$p$ data are also directly compared with $d$+$A$ data in mass and $p_T$, and within the statistical accuracy no nuclear modification is seen.

The PHENIX experiment has measured $\phi$ meson production in $d$$+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV using the dimuon and dielectron decay channels. The $\phi$ meson is measured in the forward (backward) $d$-going (Au-going) direction, $1.2<y<2.2$ ($-2.2<y<-1.2$) in the transverse-momentum ($p_T$) range from 1--7 GeV/$c$, and at midrapidity $|y|<0.35$ in the $p_T$ range below 7 GeV/$c$. The $\phi$ meson invariant yields and nuclear-modification factors as a function of $p_T$, rapidity, and centrality are reported. An enhancement of $\phi$ meson production is observed in the Au-going direction, while suppression is seen in the $d$-going direction, and no modification is observed at midrapidity relative to the yield in $p$$+$$p$ collisions scaled by the number of binary collisions. Similar behavior was previously observed for inclusive charged hadrons and open heavy flavor indicating similar cold-nuclear-matter effects.

Measurements of the differential production of electrons from open-heavy-flavor hadrons with charm- and bottom-quark content in $p$$+$$p$ collisions at $\sqrt{s}=200$ GeV are presented. The measurements proceed through displaced-vertex analyses of electron tracks from the semileptonic decay of charm and bottom hadrons using the PHENIX silicon-vertex detector. The relative contribution of electrons from bottom decays to inclusive heavy-flavor-electron production is found to be consistent with fixed-order-plus-next-to-leading-log perturbative-QCD calculations within experimental and theoretical uncertainties. These new measurements in $p$$+$$p$ collisions provide a precision baseline for comparable forthcoming measurements in A$+$A collisions.

We report on the nuclear dependence of transverse single-spin asymmetries (TSSAs) in the production of positively-charged hadrons in polarized $p^{\uparrow}+p$, $p^{\uparrow}+$Al and $p^{\uparrow}+$Au collisions at $\sqrt{s_{_{NN}}}=200$ GeV. The measurements have been performed at forward rapidity ($1.4<\eta<2.4$) over the range of $1.8<p_{T}<7.0$ GeV$/c$ and $0.1<x_{F}<0.2$. We observed a positive asymmetry $A_{N}$ for positively-charged hadrons in \polpp collisions, and a significantly reduced asymmetry in $p^{\uparrow}$+$A$ collisions. These results reveal a nuclear dependence of charged hadron $A_N$ in a regime where perturbative techniques are relevant. These results provide new opportunities to use \polpA collisions as a tool to investigate the rich phenomena behind TSSAs in hadronic collisions and to use TSSA as a new handle in studying small-system collisions.

During 2015 the Relativistic Heavy Ion Collider (RHIC) provided collisions of transversely polarized protons with Au and Al nuclei for the first time, enabling the exploration of transverse-single-spin asymmetries with heavy nuclei. Large single-spin asymmetries in very forward neutron production have been previously observed in transversely polarized $p$$+$$p$ collisions at RHIC, and the existing theoretical framework that was successful in describing the single-spin asymmetry in $p$$+$$p$ collisions predicts only a moderate atomic-mass-number ($A$) dependence. In contrast, the asymmetries observed at RHIC in $p$$+$$A$ collisions showed a surprisingly strong $A$ dependence in inclusive forward neutron production. The observed asymmetry in $p$$+$Al collisions is much smaller, while the asymmetry in $p$$+$Au collisions is a factor of three larger in absolute value and of opposite sign. The interplay of different neutron production mechanisms is discussed as a possible explanation of the observed $A$ dependence.