OK, here is my argument (sorry for the delay).
First of all, $Z$ is essentially the maximum of the absolute value of the polynomial $P(z)=\prod_j(1-z_jz)$ on the unit circumference (up to a factor of $n$, but it is not noticeable on the scale we are talking about).
Second, the maximum of the absolute value of a (trigonometric) polynomial of degree $K$ can be read from any $AK$ uniformly distributed points on the unit circumference $\mathbb T$ (say, roots of unity of degree $AK$) with relative error of order $A^{-1}$.
Now let $\psi(z)=\log(1-z)$. We want to find the asymptotic distribution of the $\max_z n^{-1/2}Re\sum_j \psi(zz_j)$ where $z_j$ are i.i.d. random variables uniformly. Since it is the logarithm of $|P(z)|$, the maximum can be found using $10n$ points.
Decompose $\psi(z)$ into its Fourier series $-\sum_{k\ge 1}\frac 1kz^k$. Then, formally, we have $n^{-1/2}\sum_j \psi(zz_j)=-\sum_k\left(n^{-1/2}\sum_j z_j^k\right)\frac{z^k}{k}$. It is tempting to say that the random variables $\xi_n,k=n^{-1/2}\sum_j z_j^k$ converge to the uncorrelated standard complex Gaussians $\xi_k$ by the CLT in distribution and, therefore, the whole sum converges in distribution to the random function $F(z)=\sum_k\xi_k\frac{z^k}k$, so $n^{-1/2}\log Z$ converges to $\max_z\Re F(z)$ (the $-$ sign doesn't matter because the limiting distribution is symmetric). This argument would be valid literally if we had a finite sum in $k$ but, of course, it is patently false for the infinite series (just because if we replace $\max$ by $\min$, we get an obvious nonsense in the end). Still, it can be salvaged if we do it more carefully.
Let $K$ run over the powers of $2$. Choose some big $K_0$ and apply the above naiive argument to $\sum_{k=1}^{K_0}$. Then we can safely say that the first $K_0$ terms in the series give us essentially the random function $F_{K_0}(z)$ which is the $K_0$-th partial sum of $F$ when $n$ is large enough.
Our main task will be to show that the rest of the series cannot really change the maximum too much. More precisely, it contributes only a small absolute error with high probability.
To this end, we need
Lemma: Let $f(z)$ be an analytic in the unit disk function with $f(0)=0$, $|\Im f|\le \frac 12$. Then we have $\int_{\mathbb T}e^{\Re f}dm\le \exp\left(2\int_{\mathbb T}|f|^2dm\right)$ where $m$ is the Haar measure on $\mathbb T$.
Proof: By Cauchy-Schwartz,
$$
\left(\int_{\mathbb T}e^{\Re f}dm\right)^2\le \left(\int_{\mathbb T}e^{2\Re f}e^{-2|\Im f|^2}dm\right)\left(\int_{\mathbb T}e^{2|\Im f|^2}dm\right)
$$
Note that if$|\Im w|\le 1$, we have $e^{\Re w}e^{-|\Im w|^2}\le \Re e^w$. So the first integral does not exceed $\int_{\mathbb T}\Re e^{2f}dm=\Re e^{2f(0)}=1$. Next, $e^s\le 1+2s$ for $0\le s\le\frac 12$, so $\int_{\mathbb T}e^{2|\Im f|^2}dm\le 1+4\int_{\mathbb T}|\Im f|^2dm\le 1+4\int_{\mathbb T}|f|^2dm$. Taking the square root turns $4$ into $2$ and it remains to use that $1+s\le e^s$
The immediate consequence of Lemma 1 is a Bernstein type estimate for $G_K(z)=\sum_{k\in (K,2K]}\left(n^{-1/2}\sum_j z_j^k\right)\frac{z^k}{k}$
$$
P(\max|\Re G_K|\ge 2T)\le 20Ke^{-T^2K/9}
$$
if $0\le TK\le \sqrt n$, say.
Indeed, just use the Bernstein trick on the independent random shifts of $g_K(z)=\sum_{k\in (K,2K]}\frac{z^k}{k}$:
$$
E e^{\pm t\Re G_K(z)}\le \left(\int_{\mathbb T}e^{\Re tn^{-1/2}g_K}dm\right)^n\le e^{2t^2/K}
$$
for every $t\le \sqrt n/2$ (we used the Lemma to make the last estimate) and put $t=\frac{TK}{3}$. After that read the maximum from $10K$ points with small relative error and do the trivial union bound.
Choosing $T=K^{-1/3}$, we see that we can safely ignore the sum from $K=K_0$ to $K=\sqrt n$ if $K_0$ is large enough. Now we are left with
$$
G_K(z)=\sum_{k\ge \sqrt n}\left(n^{-1/2}\sum_j z_j^k)\right)\frac{z^k}{k}
$$
to deal with. Recall that all we want here is to show that it is small at $10n$ uniformly distributed points. Again, if $g(z)=\sum_{k\ge \sqrt n}\frac{z^k}{k}$, we have $|\Im g|\le 10$, say so we can use the same trick and get
$$
P(\max_{10n\text{ points}}|\Re G|\ge 2T)\le 20n e^{n^{-1/2}t^2-tT}
$$
if $0\le t\le \sqrt n/20$, say.
Here we do not need to be greedy at all: just take a fixed small $T$ and choose $t=\frac{2\log n}T$.
Now, returning to your original determinant problem, we see that the norm of the inverse matrix is essentially $Z/D$ where $D=\min_i\prod_{j:j\ne i}|z_i-z_j|$. We know the distribution of $\log Z$ and we have the trivial Hadamard bound $D\le n$. This already tells you that the typical $\lambda_1$ is at most $e^{-c\sqrt n}$. The next logical step would be to investigate the distribution of $\log D$.
