How find this matrix $A=(\sqrt{i^2+j^2})$ eigenvalue

As David Speyer's numerical experiments suggest, the quadratic form associated to this symmetric matrix is negative-definite on the hyperplane $c_1 + \cdots + c_n = 0$, and thus has signature $(1,n-1)$ because it is positive on the unit vectors. This is the special case $a_i = i^2$, $s = 1/2$ of the following result:

Proposition. Let $a_1,\ldots,a_n$ be distinct positive real numbers and $s \in (0,1)$. Then the quadratic form $$ Q(c_1,\ldots,c_n) = \sum_{i=1}^n\sum_{j=1}^n (a_i+a_j)^s c_i c_j $$ is negative-definite on the hyperplane $c_1 + \cdots + c_n = 0$.

Proof: We use the integral representation $$ a^s = \frac{s}{\Gamma(1-s)} \int_{x=0}^\infty (1-e^{-ax}) \, x^{-s} \frac{dx}{x}, $$ which holds for all $a>0$, and follows from the Gamma integral $\int_0^\infty e^{-ax} x^{-s} dx = \Gamma(1-s) \, a^{s-1}$ by integration by parts. It follows that $$ Q(c_1,\ldots,c_n) = \int_{x=0}^{\infty} \frac{s}{\Gamma(1-s)} (\,f(0)^2 - f(x)^2) \, x^{-s} \frac{dx}{x}$$ where $f(x) = \sum_{i=1}^n c_i e^{-a_i x}$. If $c_1+\cdots+c_n = 0$ then $f(0)=0$, and then the integrand is $-f(x)^2 \, x^{-s} dx/x$, which is everywhere $\leq 0$, and not identically zero unless $c_i=0$ for all $i$. Therefore $Q(c_1,\ldots,c_n) \leq 0$ with equality only at zero, QED.

Report on some basic experiments. For $n \leq 30$, there is one positive eigenvalue and all the others are negative.

I checked this with the following Mathematica command:

 mm[n_] := Table[Sqrt[i^2 + j^2], {i, 1, n}, {j, 1, n}]

 Table[Count[Sign[Eigenvalues[SetPrecision[mm[n], 50]]], 1], {n, 1, 30}]

SetPrecision tells Mathematica to treat the square roots as floating point numbers with $50$ decimal digits of accuracy. If you tell it to treat them as exact quantities, the computation times out; if you use the default accuracy it won't get the signs of the smallest eigenvalues right. The smallest eigenvalues here are around $10^{-30}$, so you need to be careful.

Probably the easiest way to prove this would be to exhibit an $n-1$ dimensional subspace on which this matrix is negative definite. I took my first guess, the span of the vectors $(1, -1,0,0,0,\ldots)$, $(0,1,-1,0,0,\ldots)$, $(0,0,1,-1,0,0,\ldots)$, ...

(* Change of basis matrix to the n-1 dimensional space. *)
ss[n_] := Table[If[i == j, 1, If[i == j + 1, -1, 0]], {i, 1, n}, {j, 1, n - 1}]

(* Quadratic form in the new basis. *)
qq[n_] := Transpose[ss[n]].mm[n].ss[n]

Table[PositiveDefiniteMatrixQ[SetPrecision[-qq[n], 50]], {n, 2, 30}]

For $n \leq 30$, the qudratic form is negative definite on this $n-1$ plane.