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)$...