Note that $(k+1)^2 - k^2 = 2k+1$. Thus, we have
$$
((m+1)^2,(m+2)^2,\dots,(m+n+1)^2) = \\
(m^2,(m+1)^2,\dots,(m+n)^2) + (2m+1,2m+3,\dots,2m+2n+1)
$$
Conclude (inductively, if you like) that the vectors
$$
(1,2,\dots,n^2),(1,3,\dots,2n+1),(1,\dots,1)
$$
span the row/column space of $A$. Or, if you prefer, you could use the basis
$$
(1,\dots,1),(0,1,\dots,n-1), (0,1^2,\dots,(n-1)^2)
$$
More generally: (I think that) as a consequence of the Newton interpolation formula, the rank of the matrix
$$
\pmatrix{f(1)& f(2) & \cdots & f(n)\\
f(2) & f(3) & \cdots & f(n+1)\\
\vdots & \vdots & \ddots & \vdots\\
f(n) & f(n+1) & \cdots & f(2n+1)}
$$
will have rank equal to at most $1 + \deg(f)$ (exactly $1 + \deg(f)$, whenever $n \geq 1 + \deg(f)$) for any polynomial $f$.
Using the Schur complement, we know that the determinant of this matrix can be expressed as
$$
\det \pmatrix{I & B\\ B^T & 0} = \det\pmatrix{I & 0\\0 & -B^TB} = \det(-B^TB).
$$
We therefore see that it is indeed the case that this matrix has non-zero determinant if (and only if) $B$ has full column rank.
Regarding the eigenvalues of this matrix, the singular value decomposition (SVD) of $B$ can be used as a helpful intermediate step. In particular, if $B = U \Sigma V^T$ is a singular value decomposition, then the matrix of interest is similar to
$$
\pmatrix{U & 0\\0 & V}^T \pmatrix{I & B\\B^T & 0}\pmatrix{U & 0\\0 & V} = \pmatrix{I & \Sigma\\ \Sigma^T & 0}.
$$
If $m \leq n$, then there exists a permtuation matrix $P$ such that
$$
P\pmatrix{I & \Sigma\\ \Sigma^T & 0}P^T = \pmatrix{A_1 \\ & \ddots \\ && A_m\\ &&& 0},
$$
where each block $A_k$ has the form
$$
A_k = \pmatrix{1 & \sigma_k\\ \sigma_k & 0 },
$$
where $\sigma_k$ denotes the $k$th singular value of $B$. It follows that the eigenvalues of the matrix will be equal to $0$ and the solutions to
$$
\lambda^2 - \lambda - \sigma_k^2 = 0, \quad k = 1,\dots,n \implies\\
\lambda = \frac{1 \pm \sqrt{1 + 4\sigma_k^2}}{2}, \quad k = 1,\dots,n.
$$
In the case that $m>n$, we have essentially the same thing except that there is now an identity matrix of size $m-n$ added to the diagonal. Consequently, the matrix will have (in addition to the eigenvalues noted in the previous case) $1$ as an eigenvalue of multiplicity at least $m-n$.
Best Answer
One can see it in a different way :
With your notation for a column vector
$$z=\begin{bmatrix} \zeta_{1} \\ \zeta_{2} \\ \vdots \\ \zeta_{n} \end{bmatrix}$$
Let M be a matrix of the form
$$M=[z|v_2|v_3|\cdots | v_n]$$
such that the columns constitute a basis of $\mathbb{R}^n$ (it is always possible to "complete" the basis of a finite dimensional linear space in this way).
Remark : $M$ is therefore invertible.
Then, we can write $$zz^T=MJM^T \ \ \text{where} \ \ J=diag(1,0,0, \cdots 0)$$
(the "1" coefficient selects $zz^T$ and the 0 zero coefficients eliminate any influence of vectors $v_k$).
We can conclude that $rank(zz^T)=rank(J)=1$ because the rank in preserved by pre- or post- multiplication by invertible matrices.
As a consequence $\det(zz^T)=\underbrace{(1 \times 0 \times 0 \cdots \times 0)}_{\det J}\det(M)^2 = 0. :)$