From this

By $\|z_1\|$, I understand you mean the maximal modulus of the $z_i$'s.

If you are interested in the process of the $\|z_i\|$'s, you have no chance for a determinantal structure since you may have two different $z_i$'s with same modulus with non zero probability.

Concerning the process of the $z_i$'s, note that even if the Ginibre Ensemble (that is the case $\mathbb{R}^2$) is indeed determinantal, its kernel is related to the polynomials orthogonal with $\exp(-\|x\|^2/2)$, that is the $(z^k)_{k\geq 0}$ ... which have trivial zeros ! My point is that except on $\mathbb{R}$ you won't get so much information concerning $\|z_1\|$ from a determinantal structure.

I don't know how prove the convergence of $\|z_1\|$, but note that from your density expression, once renormalized $z_i\rightarrow z_i/\sqrt{N}$, you still can use the Coulomb-gaz approach and to characterize the global distribution of the $z_i$'s (for example prove by proving a large deviation principle to) obtain that for the limiting distribution empirical measure) : It is given by the unique minimizer $\mu^*$ of the functional

For the distribution of $z_1$, you can compute explicitly its distribution at $N$ fixed, and it is gaussian : Integrate out $z_2, \ldots , z_N$ in $p(z_1,\ldots,z_N)$, I guess that $\|z_1\|$ should converge towards$p(z_1)= C'_N e^{-N\|z_1\|^2/2}\int\ldots\int \prod_{2\leq i < j \max \leq N}\|z_i-z_j\|^2\prod_{j=2}^N\|z_1-z_j\|^2e^{-N\|z_j\|^2/2}dz_jand do $z_j big(Supp(\mu^*)\cap \rightarrow z_j+z_1$ to obtainmathbb{R}\big)$...

From this density, once renormalized $z_i\rightarrow z_i/\sqrt{N}$, you still can use the Coulomb-gaz approach and (for example prove a large deviation principle to) obtain that the limiting distribution is given by the unique minimizer of the functional $$ \iint\log\frac{1}{\| x-y\|}d\mu(x)d\mu(y) +\frac{1}{2}\int \|x\|^2d\mu(x) $$ over probability measures $\mu$ on $\mathbb{R}^3$ (or higher).

For the distribution of $z_1$, you can compute explicitly its distribution at $N$ fixed, and it is gaussian : Integrate out $z_2, \ldots , z_N$ in $p(z_1,\ldots,z_N)$, $$ p(z_1)= C'_N e^{-N\|z_1\|^2/2}\int\ldots\int \prod_{2\leq i < j \leq N}\|z_i-z_j\|^2\prod_{j=2}^N\|z_1-z_j\|^2e^{-N\|z_j\|^2/2}dz_j $$ and do $z_j \rightarrow z_j+z_1$ to obtain$$ p(z_1)= C'_N e^{-\|z_1\|^2/2}\int\ldots\int \prod_{j=2}^N\|z_1-z_j\|^2dz_j = \frac{1}{(2\pi)^{3/2}}e^{-\|z_1\|^2/2}. $$

From this density, you still can use the Coulomb-gaz approach and (for example prove a large deviation principle to) obtain that the limiting distribution is given by the unique minimizer of the functional $$ \iint\log\frac{1}{\| x-y\|}d\mu(x)d\mu(y) +\frac{1}{2}\int \|x\|^2d\mu(x) $$ over probability measures $\mu$ on $\mathbb{R}^3$ (or higher).

For the distribution of $z_1$, you can compute explicitly its distribution at $N$ fixed, and it is gaussian : Integrate out $z_2, \ldots , z_N$ in $p(z_1,\ldots,z_N)$, and do $z_j \rightarrow z_j+z_1$ to obtain $$ p(z_1)= C'_N e^{-\|z_1\|^2/2}\int\ldots\int \prod_{j=2}^N\|z_1-z_j\|^2dz_j = \frac{1}{(2\pi)^{3/2}}e^{-\|z_1\|^2/2}. $$