Solution (source code)

= Solution

Use the dimensionally correct flat <Friedmann-Lemaître-Robertson-Walker metric>
$$
ds^2=-c^2dt^2+a^2(t)\left[d\chi^2+\chi^2(d\vartheta^2+\sin^2\vartheta\,d\varphi^2)\right].
$$
The source's $-c\,dt^2$ is missing a factor $c$, and its angular term must mean the unit-sphere line element $d\Omega^2$, not the <solid angle> element $d\Omega=\sin\vartheta\,d\vartheta\,d\varphi$. For example, $c\,dt^2$ has dimensions length times time, while $ds^2$ has dimensions length squared. The corrected convention above removes both ambiguities.

A radial light ray has $ds^2=0$ and $|d\chi|=c|dt|/a$. Thus
$$
\frac{d\chi}{dz}=\frac c{H(z)},\qquad
\chi(z)=\frac c{H_0}\int_0^z\frac{du}{(1+u)^{3/2}}
=\frac{2c}{H_0}\left(1-\frac1{\sqrt{1+z}}\right).
$$
A proper shell volume at emission is $dV_{\rm proper}=a^3\chi^2d\chi\,d\Omega$. Multiplying by $n(z)=n_0a^{-3}$ cancels the expansion factors, equivalently using the comoving volume element directly. Hence
$$
\frac{d^2N}{d\Omega\,dz}=n_0\chi^2\frac{d\chi}{dz}
=\boxed{\frac{4c^3n_0(\sqrt{1+z}-1)^2}{H_0^3(1+z)^{5/2}}.}
$$
This counts every galaxy in the shell per <solid angle>, rather than intersections weighted by individual galaxy areas. Near $z=0$ it reduces to $n_0(c/H_0)^3z^2$, as expected for a Euclidean volume shell.