We discuss the transverse single spin asymmetries $A_N$ for very forward
neutrons measured by the PHENIX zero-degree calorimeters (ZDCs) in high-energy
polarized proton-nucleus ($p$-A) collisions at Relativistic Heavy Ion Collider
(RHIC). The ZDCs cover the psudorapidity range of $6.8<\eta<8.8$ where the neutron
emission angle with respect to the incident proton direction is limited to less
than 2.2$\,$mrad. First-ever $p$-A data taken in the RHIC-2015 run exhibit positive and remarkably
large $A_N\sim0.18$ only in $p$-Au collisions whereas nearly zero $A_N$ in
$p$-Al collisions. In this proceeding, to explain such $A_N$ in $p$-A collisions we present an
important but rather unknown mechanism: ultra-peripheral $p$-A collisions (UPCs,
also known as Primakoff effects). UPCs lead to very large $A_N$ of about 0.35
and have cross sections proportional to $Z^2$ of the nuclei.
UPCs contribute to inclusively measured $A_N$ modestly in $p$-Al collisions and
significantly in $p$-Au collisions. We show that the Monte Carlo simulations
incorporating both a one-pion exchange (OPE) and UPCs, where virtual
photon-polarized proton interactions follow the MAID 2007 isobar model,
successfully describe the PHENIX data in $p$-Al and $p$-Au collisions
simultaneously.