Define a mean $m$ with $n m = \sum_{i=1}^{n}x_{i}$.
Define deviations $d_i = x_i -m$, where obviously $\sum_{i=1}^{n}d_{i} = 0$. Since $1 = \sum_{i=1}^{n}x_{i}^2 = n m^2 + \sum_{i=1}^{n}d_{i}^2$ we have $|m| < 1/\sqrt n$. Let's continue working with fixed $m$ in the following, and come back to the discussion of $m$ later.
Now write the target function
$$
F = \sum_{i=1}^{n}x_{i}\sum_{i=1}^{n}x^3_{i} = nm (nm^3 + 3 m \sum_{i=1}^{n}d_{i}^2 + \sum_{i=1}^{n}d_{i}^3) = n^2m^4 + 3 n m^2(1-nm^2) + n m \sum_{i=1}^{n}d_{i}^3
$$
For the last term, we would like to use some known result. Let $d_i = \sqrt{1 - n m^2} a_i$, then $\sum_{i=1}^{n}a_{i}^2 = 1$ and
$$
F = n^2m^4 + 3 n m^2(1-nm^2) + n m (1 - nm^2)^{3/2} \sum_{i=1}^{n}a_{i}^3
$$
So for bounding the last term, the question arises about the maximum and minimum of $\sum_{i=1}^{n}a_{i}^3$ under the constraints $\sum_{i=1}^{n}a_{i} = 0$ and $\sum_{i=1}^{n}a_{i}^2 = 1$. This has been solved for example here, showing that
$$
- \frac{n-2}{\sqrt{n(n-1)}} \le \sum_ia_i^3 \le \frac{n-2}{\sqrt{n(n-1)}}.
$$
Building this into the target function, we get that the maximum (for some $m$) is achieved for $m$ and $\sum_ia_i^3$ having the same sign, giving
$$
F_{max}(|m|) = n^2m^4 + 3 n m^2(1-nm^2) + n |m| (1 - nm^2)^{3/2} \frac{n-2}{\sqrt{n(n-1)}}
$$
Conversely, the minimum is achieved for $m$ and $\sum_ia_i^3$ having opposite sign, giving
$$
F_{min}(|m|) = n^2m^4 + 3 n m^2(1-nm^2) - n |m| (1 - nm^2)^{3/2} \frac{n-2}{\sqrt{n(n-1)}}
$$
Since $n m^2 \le 1$, one can write $nm^2 = \cos(x)^2$ and bring this into the form
$$
F(x) = \cos(x)^4 + 3 \cos(x)^2\sin(x)^2 \pm \cos(x) \sin(x)^3 \frac{n-2}{\sqrt{n-1}} \quad \tag{1}
$$
however this does not appear to provide easier solutions, so it was not followed further on. It is useful for a large-$n$ discussion, see bottom.
This leaves us with the discussion of $m$. We have that the global $
F_{max} = {\max}_{|m|}{F_{max}(|m|)}
$ and likewise the global $
F_{min} = {\min}_{|m|}{F_{min}(|m|)}
$.
For $n=2$, we have in the range $0 \le |m| \le 1/\sqrt{2}$ that $F_{max} = {\max}_{|m|}4m^4 + 6 m^2(1-2m^2) = 9/8 $ at $|m|≈0.61 < 1/\sqrt 2 $ and $F_{min} = {\min}_{|m|}4m^4 + 6 m^2(1-2m^2) = 0 $ at $|m|=0$ , as was already stated in the question.
For general $n$, closed form results are not directly obtainable. For example, for $n=3$ we have (using WolframAlpha) that
$
F_{max} = {\max}_{|m|}9 m^4 + 9 m^2 (1 - 3 m^2) + \sqrt{3/2} |m| (1 - 3 m^2)^{3/2}
≈1.22307$ at $|m|≈0.46416 < 1/\sqrt 3$
and
$
F_{min} = {\min}_{|m|}9 m^4 + 9 m^2 (1 - 3 m^2) - \sqrt{3/2} |m| (1 - 3 m^2)^{3/2}
≈-0.040397$ at $|m|≈0.065262 < 1/\sqrt 3$
For $n=4$, WolframAlpha provides $
F_{max}≈1.30667$ at $|m|≈0.38302$ and $F_{min}≈-0.102869$ at $|m|≈0.086824$.
One observes that with rising $n$, maxima get higher and minima lower.
Indeed, when observing eq. (1), for large $n$ the first two terms have no relation to $n$ and stay of order $1$, whereas the last term grows with roughly $\propto \sqrt{n}$. Since the maximum of $ \cos(x) \sin(x)^3$ is $3 \sqrt{3}/16 \simeq 0.325$, we have that $F_{max} {{n \; {\text{large}}} \atop {\to}} \simeq 0.325 \sqrt{n}$ and $F_{min} {{n \; {\text{large}}} \atop {\to}} \simeq - 0.325 \sqrt{n}$.