  The degree of Gaussianity of a field offers insights into its cosmological
nature, and its statistical properties serve as indicators of its Gaussianity.
In this work, we examine the signatures of Gaussianity in a gravitational wave
background (GWB) by analyzing the cumulants of the one- and two-point functions
of the relevant observable, using pulsar timing array (PTA) simulations as a
proof-of-principle. This appeals to the ongoing debate about the source of the
spatially-correlated common-spectrum process observed in PTAs, which is likely
associated with a nanohertz stochastic GWB. We investigate the distribution of
the sample statistics of the one-point function in the presence of a Gaussian
GWB. Our results indicate that, within PTAs, one-point statistics are
impractical for constraining the Gaussianity of the nanohertz GWB due to
dominant pulsar noises. However, our analysis of two-point statistics shows
promise, suggesting that it may be possible to constrain the Gaussianity of the
nanohertz GWB using PTA data. We also emphasize that the Gaussian signatures
identified in the one- and two-point functions in this work are expected to be
applicable to any gravitational wave background.
