Well, this is not actually an answer to either of the OP's questions; at most, it provides an easier way to classify the solutions for the $n=3$ case, and that might point a way towards an analysis for larger $n$, but I don't make any guarantees.
Much that I'll say to begin with works for all $n$, so I'll start with some general remarks. First, let me rephrase the problem in terms of quartic polynomials in $n$ variables, since that is what $T$ actually is. Let $x_1,\ldots,x_n$ be standard orthogonal coordinate functions on $\mathbb{R}^n$, and write $T = T_{ijkl}x_ix_jx_kx_l$ (summation convention assumed here and henceforth). Then the given condition is equivalent to a nonlinear second-order PDE that is invariant under $\mathrm{O}(n)$:
$$
R(T) := \mathrm{tr}\bigl(\mathrm{Hess}(T)^2\bigr) - 48\,T = 0.
$$
(This trivial rewrite makes it easy to do computations without having to worry about symmetrization of tensors, etc.)
Under the action of $\mathrm{O}(n)$, the homogeneous quartic polynomials break up into the sum of three irreducible representations: the harmonic quartics, the product of harmonic quadratics with $|x|^2 = {x_1}^2+\cdots+{x_n}^2$, and constant multiples of $|x|^4 = ({x_1}^2+\cdots+{x_n}^2)^2$. Thus, we can write $T = T_4 + T_2\,|x|^2 + T_0\,|x|^4$ where $\Delta(T_i)=0$. By representation theory, it is easy to see that $R(T)_4$ is a positive linear combination of the harmonic components of $\Delta^2({T_4}^2)$, $\Delta(T_2T_4)$, ${T_2}^2$, $T_0T_4$ and $-48 T_4$ and that $R(T)_0$ is a positive linear combination of $\Delta^4({T_4}^2)$, $\Delta^2({T_2}^2)$ and $T_0\left(T_0-\frac{3}{n{+}8}\right)$
As a consequence, for example, one finds that if $R(T)=0$ and $T_4=0$, then $T_2=0$ and $T_0$ is either $0$ or $\frac{3}{n{+}8}$. More generally, since the quadratic forms $\Delta^4({T_4}^2)$ and $\Delta^2({T_2}^2)$ are positive definite on their respective irreducible subspaces, $R(T)_0=0$ implies that $0\le T_0 \le \frac{3}{n{+}8}$ and that $T_0=0$ only when $T=0$ while $T_0=\frac{3}{n{+}8}$ only when $T = \frac{3}{n{+}8}\,|x|^4$ (which does satisfy $R(T)=0$). Moreover the bound on $T_0$ and the vanishing of $R(T)_0$ imply bounds on the quadratic forms $\Delta^4({T_4}^2)$ and $\Delta^2({T_2}^2)$, and, hence, the space of solutions to $R(T)_0=0$ is compact. Thus, for reasonably small $n$, a numerical search for the solutions to $R(T)=0$ has a chance of giving a good description of the space of solutions, once one has taken care of the redundancy in solutions caused by the $\mathrm{O}(n)$ invariance.
One more remark about 'trivial' solutions: If $T'$ satisfies $R(T')=0$ and only involves the variables $x_1\,\dots,x_k$ while $T''$ satisfies $R(T'')=0$ and only involves the variables $x_{k+1},\ldots,x_n$, then $R(T'+T'') = R(T')+R(T'') = 0$. I'll say that a solution $T$ to $R(T)=0$ is reducible if $T$ is $\mathrm{O}(n)$-equivalent to $T'+T''$ as above for some $k$ with $1<k<n$. (It is not entirely obvious how to test whether a given solution be reducible, but a necessary condition is that $\det(\mathrm{Hess}(T))$ either be zero or have nontrivial factors.)
Because the equation $R(T)=0$, though determined, is $\mathrm{O}(n)$-invariant, the solutions will not be finite in number when $n>1$. The most we can hope for (and even this is not clear) is that the solution set might be a finite union of $\mathrm{O}(n)$-orbits. What is needed is a way to take a cross-section of the orbits and look for solutions within that cross-section. (Another approach, which does not appear to be computationally feasible, would be to parametrize the (Hausdorff) space of $\mathrm{O}(n)$-orbits in some way and work on this reduced space.) Unfortunately, the space of quartic polynomials is not a polar representation of $\mathrm{O}(n)$ (i.e., there is no linear subspace that meets each of the $\mathrm{O}(n)$-orbits orthogonally). However, as is true for any finite-dimensional representation of a compact group $G$ for which the typical orbit has a finite stabilizer, there will be a subspace whose codimension is equal to $\dim G$ that meets each orbit, typically in a finite number of points. Usually, one wants to choose such a subspace with as large a symmetry group as possible, so that the $G$-orbits will intersect the subspace in a symmetric pattern.
To choose such a subspace, consider the fact that, when $n=1$, the only solution is $T = \frac13\,{x_1}^4$, which leads, by the above remarks, to the (reducible) 'Fermat' solution
$$
\bar T = \tfrac13\,({x_1}^4 + \cdots + {x_n}^4),
$$
whose symmetry group in $\mathrm{O}(n)$ is the signed permutations on $n$ letters, a group of order $2^n\,n!$. The tangent space at $\bar T$ to the $\mathrm{O}(n)$-orbit of $\bar T$ is spanned by the quartics $x_i{x_j}^3 - x_j{x_i}^3$ for $i<j$, and the orthogonal subspace to this tangent space is determined by the $T$ for which $T_{iiij}=T_{ijjjj}$ for all $i,j$.
Thus, every quartic $T$ is on the $\mathrm{O}(n)$-orbit of a $T'$ that satisfies these $\frac12 n(n{-}1)$ linear equations. Restricting to this subspace now yields an over-determined system of equations for the coefficients of $T$, and one can hope that there will be only a finite number of solutions.
Using this strategy for $n=2$ quickly yields that there are four $\mathrm{O}(2)$-orbits that satisfy $R(T)=0$, and they are represented by
$$
0,\quad \tfrac13\,{x_1}^4,\quad \tfrac13\,({x_1}^4+{x_2}^4),\quad \tfrac3{10}\,({x_1}^2+{x_2}^2)^2.
$$
Using this strategy for $n=3$ (and employing Maple and Groebner bases) eventually yields that there are nine $\mathrm{O}(3)$-orbits that satisfy $R(T)=0$. Six of these are reducible, so they can quickly be constructed from the solutions for $n=1$ and $n=2$. The remaining three are irreducible. One of them is the orbit of the already-mentioned $\mathrm{O}(3)$-invariant $T=\frac{3}{11}({x_1}^2+{x_2}^2+{x_3}^2)^2$. The remaining two are described as follows: One is the orbit of
$$
T = \tfrac29\bigl({x_1}^4 + {x_2}^4 + {x_3}^4
+ 3\,{x_1}^2{x_2}^2+ 3\,{x_2}^2{x_3}^2+ 3\,{x_3}^2{x_1}^2\bigr),
$$
with $\Delta^2(T) = 32$, and the other is the orbit of
$$
T = c_1\,{x_3}^4 + c_2\,{x_3}^2({x_1}^2+{x_2}^2) + c_3\,({x_1}^2+{x_2}^2)^2,
$$
where $c_1$, $c_2$, $c_3$ are real numbers that lie in the field of degree $3$ over $\mathbb{Q}$ generated by $p = \Delta^2(T)\approx 32.09$, the real root of
$$
211\,p^3 −18264\,p^2+526144\,p−5048832=0.
$$
(It's now easy to determine $c_1$, $c_2$, and $c_3$, but I won't write them out here, as the formulae are not illuminating.)
When I have time, I'll remark on the actual calculations. The key is to use Maple and Groebner basis techniques to find the polynomial equation satisfied by $s = \Delta^2(T) = 120\, T_0$. This polynomial has 9 real roots, and, using this and a few more observations, one can coax Maple into listing all of the solutions with a given $\Delta^2(T)$ that lie in the 12-dimensional cross-section. However, each of these orbits can intersect this space several times (for example, the orbit of $\tfrac1{3}({x_1}^4{+}{x_2}^4{+}{x_3}^4)$ meets this subspace in 14 distinct points), and a few simple ideas are needed to show that all of the solutions with the same value of $s = \Delta^2(T)$ lie on a single $\mathrm{O}(3)$-orbit.
What I find remarkable is that all of the orbits except for the last one mentioned above have simple rational representatives. (Also, in my initial comments, I didn't recognize that all of the 'irrational' solutions actually do lie on a single orbit. I had thought that there were two distinct orbits, but, eventually, I realized that they were equal.)
Addendum: The 9 representative solutions listed for the $\mathrm{O}(3)$-orbits all share a striking feature: They are all quadratic polynomials in ${x_1}^2,{x_2}^2,{x_3}^2$. Of course, if one knew this ahead of time, it would have been relatively easy to find all of the solution orbits. This inspired me to ask how many orbits in the $n=4$ case could be represented this way, so I tried this out. (Using the general strategy outlined above for the full $n=4$ case appears to be beyond the capabilities of my laptop.)
Now, there are 10 unknown coefficients, but Maple with Groebner bases is able to compute and factor the polynomial that must be satisfied by $\Delta^2T = 192\,T_0$. It turns out that there are 13 rational roots and 4 irrational ones: Two of the latter are the expected $p$ and $p+8$, where $p\approx32.09$ is as above in the $n=3$ case, but the remaining two, with approximate values $47.83$ and $41.16$, are the unique real roots of polynomials irreducible over $\mathbb{Q}$ of degrees 3 and 11. Based on how the answer turned out when $n=3$, I expected that there would be 17 solution orbits of this 'bi-quadratic' kind when $n=4$.
However, what turned out was quite different: The two 'new' real irrational roots do not correspond to any real $\mathrm{O}(4)$-orbit. (They do occur for some complex solutions.) Meanwhile, there are obviously two orbits with $\Delta^2(T) = 32$ and at least one each for the remaining 15 real roots. While I have not yet had time to check all of them, I suspect that each of these 15 roots belongs to a single $\mathrm{O}(4)$-orbit. If that is so, then there will be exactly 16 solution orbits that are representable by quadratic expressions in ${x_1}^2,{x_2}^2,{x_3}^2,{x_4}^2$.
This makes me wonder whether there are any known solutions with $n=4$ that are not equivalent to one of these 16 'known' ones. If this does turn out to be all of the (real) solutions, then the sequence $1,4,9,16,\ldots$ as the number of (real) solution orbits as $n$ increases is certainly intriguing, though, on the face of it, it seems improbable that it would continue to hold.
