There is a proof of this using the theory of $\ell$-adic sheaves. There should also be an elementary proof.

This family of exponential sum is associated to a rank 2 sheaf (a slight variant of a hypergeometric sheaf). This means the exponential sum corresponds to the trace of Frobenius, a unitary matrix, on a rank two representation. The argument of the trace of a unitary matrix is half the argument of the determinant.

So your statement is equivalent to the statement that the determinant of the family of exponential sums is constant. The determinant of an individual exponential sum can be calculated using local epsilon factors. For a one-parameter exponential sum like this, if we only want to calculate the ratios between the different determinants, it can be easier. Because your family of sums can be expressed as a Fourier transform, it is easy to apply Laumon's theory of the $\ell$-adic Fourier transform to calculate the local monodromy of the sheaf, take its determinant, and verify that it is trivial. Because the determinant sheaf has trivial local monodromy it has trivial global monodromy and must be constant.

But this should all be overkill by your problem. If I have done these local calculations correctly then your sum is some scalar, probably a ratio of Gauss sums, times $$\sum_{x,y \in \mathbb F_p, xy = c(gn)^2} \chi(x/y) \zeta_p^{x+y}$$ for a constant $c \in \mathbb F_p$, and this renormalized sum is manifestly real. Multiplying by a scalar, we obtain a family of sums with fixed argument. This sort of identity, even if guessed using the $\ell$-adic methods, usually has a direct, elementary proof by some clever manipulation of terms - too clever for me to see right now, unfortunately.

EDIT - this is not actually the right formula, don't use it, I'll fix it shortly-  still, a formula like this can be computed.