  The signals from international pulsar timing arrays have presented a hint of
gravitational stochastic background in nHz band frequency. Further confirmation
will be based on whether the signals follow the angular correlation curves
formulated by the overlap reduction functions, known as Hellings-Downs curves.
This paper investigates the non-linear corrections of overlap reduction
functions in the present of non-Gaussianity, in which the self-interaction of
gravity is first taken into considerations. Based on perturbed Einstein field
equations for the second order metric perturbations, and perturbed geodesic
equations to the second order, we obtain non-linear corrections for the timing
residuals of pulsar timing, and theoretically study corresponding overlap
reduction functions for pulsar timing arrays. There is order-one correction for
the overlap reduction functions from the three-point correlations of
gravitational waves, and thus the shapes of the overlap reduction functions
with non-linear corrections can be distinguished from the Hellings-Downs
curves.
