In this work, we propose a novel approach for probing the neutron skin thickness of gold (Au) by analyzing the event-by-event distribution of ${\pi }^{-}$ and ${\pi }^{+}$ yield differences. This is achieved through SMASH simulations of ultra-peripheral Au+Au collisions at $\sqrt{s_{\rm NN}}=3$ GeV. Our results demonstrate that the mean value of $\Delta n_{\pi }=n_{{\pi }^{-}}-n_{{\pi }^{+}}$, along with the Pearson correlation and mutual information between $({\pi }^{-}+{\pi }^{+})$ and $({\pi }^{-}-{\pi }^{+})$, all scale linearly with the neutron skin thickness. Moreover, the slope of the line connecting two distinct $\Delta n_{\pi }$ values in the event-by-event distribution also exhibits a linear dependence on the neutron skin thickness. The most sensitive $\Delta n_{\pi }$ pairs are identified as $(-2,-1)$, $(-1,1)$, $(0,1)$, and $(0,2)$. These findings establish a new pathway for determining the neutron skin thickness. Finally, by comparing SMASH, JAM, and UrQMD simulations under identical initial conditions, we observed a certain degree of model dependence in the outcomes, yet the overall trends remained consistent. However, by extracting slopes from multiple $\Delta n_{\pi }$ pairs in experimental event-by-event data and inferring the corresponding neutron skin thickness, one can assess which model better aligns with the true physical value.