By gaussian I assume you mean the normal distribution, mean 0, std deviation 'a'
f(z) = N*exp(-(z/a)^2/2)
with N = 1/(sqrt(2)*sqrt(pi)*a)
y = z^2 has one sign while z can have either sign. Since the normal distribution is symmetric, to make life easier assume the initial distribution is one-sided, with z>=0 only, and double the height of the distribution to keep the normalization correct:
f(z) = 2*N*exp(-(z/a)^2/2) z>=0
Now that z >= 0 we can use z = sqrt(y). The the idea with pdfs is to find a function g(y) such that
f(z) dz = g(y) dy
with z = y^(1/2)
which means that dz = (1/2) y^(-1/2) dy
plug that in
2*N*exp(-(y/(2*a^2))) (1/2) y(-1/2) dy = g(y) dy
g(y) = N*exp(-(y/(2*a^2)))/sqrt(y)
which is a chi-squared probability distribution function (with 1 degree of freeedom) as noted by Jeff. This function has the disadvantage of being unbounded as y --> 0, although its area is still 1 as required.
If you meant y = (z/a)^2 instead of y = z^2, then the expression would be the standard chi-squared
g(y) = N*exp(-(y/2))/sqrt(y) N = 1/(sqrt(2)*sqrt(pi))