I have an other solution, which not as slick Markus Sprecher's one, but I think the method, which is quite general, is interesting in itself.
Fix $c \in \left]-1,1\right[$ and $n \geq 2$, and consider the compact manifold
$$
\Sigma = \{ (a,b) \in S^{n-1} \times S^{n-1} \ | \ \langle a , b \rangle = c \}.
$$
It has dimension $2n - 3$, and the tangent space at $(a,b)$ is the orthogonal in $\mathbb{R}^{2n}$ to the 3-space spanned by $(a,0),(0,b),(b,a)$.
Consider the functions
$$
r_1 = \sqrt{\sum_i a_i^4}, \ \ r_2 = \sqrt{\sum_i b_i^4}, \ \ r_3 = \sum_{i} a_i^2 b_i^2
$$
as continuous functions on $\Sigma$. Let $(a,b)$ be a point of $\Sigma$ where $f = r_1 r_2 + r_3$ is maximal. We would like to show that $f(a,b) \leq 1 + c^2$.
The gradient of $f$ at $(a,b)$ must be of the form $\alpha(a,0) + \beta(0,b) + \gamma(b,a)$ for some real numbers $\alpha,\beta,\gamma$. This gives the following $2n$ equations:
$$
(1) \ \ a_i b_i^2 + a_i^3 \frac{r_2}{r_1} = \alpha a_i + \gamma b_i \\
(2) \ \ b_i a_i^2 + b_i^3 \frac{r_1}{r_2} = \beta b_i + \gamma a_i \\
$$
Multiplying $(1)$ by $a_i$ and summing ober $i$, we get $f = r_3 +r_1r_2 = \alpha + c \gamma$.
Multiplying $(2)$ by $b_i$ and summing ober $i$, we get $f = r_3 +r_1r_2 = \beta + c \gamma$. In particular $\alpha = \beta$.
Mutliply equations $(1)$ and $(2)$ by $r_1 b_i$ and $r_2 a_i$ respectively. This yields
$$
(1)' \ \ r_1 a_i b_i^3 + r_2 a_i^3 b_i = r_1 \alpha a_i b_i + r_1 \gamma b_i^2 \\
(2)' \ \ r_2 a_i^3 b_i + r_1 a_i b_i^3 = r_2 \alpha a_ib_i + r_2 \gamma a_i^2 \\
$$
The two LHSs are equal, hence so are the RHSs. Thus $(a_i,b_i)$ satisifies the quadratic equation
$$
(3) r_1 \alpha a_i b_i + r_1 \gamma b_i^2 = r_2 \alpha a_ib_i + r_2 \gamma a_i^2.
$$
There are two cases:
- if $\gamma =0$, then $\alpha = f >0$, and thus $r_1 = r_2$ follows from $(3)$. Then $(1)-(2)$ become the alternative $(a_i,b_i) =0$ or $a_i^2+ b_i^2 = \alpha$. Since $a$ and $b$ are not proportional, the latter cas happens at least twice, and thus $2 = ||a||_2^2 + ||b||_2^2 \geq 2 \alpha$. Thus $f = \alpha \leq 1$, and in particular $f \leq 1 + c^2$.
- if $\gamma \neq 0$ then there are at most $2$ projective solutions to $(3)$. Since the system $(1)-(2)$ is inhomogeneous, it has at most $2$ non-zero solutions up to sign. Thus there are non zero vectors $(a_1'',b_1'')$ and $(a_2'',b_2'')$ such that for eachi $i$ one has either $(a_i,b_i) = \pm (a_j'',b_j'')$ for some $j=1,2$, or $(a_i,b_i) = 0$. Let $[|1,n|] = I_0 \sqcup I_1 \sqcup I_2$ such that $(a_i,b_i) = 0$ for $i$ in $I_0$, and $(a_i,b_i) = \pm (a_j'',b_j'')$ if $i$ is in $I_j$ for some $j=1,2$. We then have the identities
\begin{align*}
|I_1| (a_1'')^2 + |I_2| (a_2'')^2 &= \sum_i a_i^2 = 1 \\
|I_1| a_1''b_1'' + |I_2| a_2'' b_2'' &= \sum_i a_i b_i = c \\
|I_1| (b_1'')^2 + |I_2| (b_2'')^2 &= \sum_i b_i^2 = 1
\end{align*}
Thus the vectors $a' = (|I_1|^{\frac{1}{2}} a_1'',|I_2|^{\frac{1}{2}}a_2'')$ and $b' = (|I_1|^{\frac{1}{2}} b_1'',|I_2|^{\frac{1}{2}}b_2'')$ in $\mathbb{R}^2$ satisfy
$$
||a'||_2 =1, \ \ ||b'|| = 1, \ \ \langle a',b' \rangle = c.
$$
Moreover, one has
\begin{align*}
\sum_{i=1}^2 (a_i')^4 &= |I_1|^2 (a_1'')^4 + |I_2|^2 (a_2'')^4 \geq |I_1| (a_1'')^4 + |I_2| (a_2'')^4 = r_1^2 \\
\sum_{i=1}^2 (b_i')^4 &= |I_1|^2 (b_1'')^4 + |I_2|^2 (b_2'')^4 \geq |I_1| (b_1'')^4 + |I_2| (b_2'')^4 = r_2^2 \\
\sum_{i=1}^2 (a_i')^2(b_i')^2 &= |I_1|^2 (a_1'')^2(b_1'')^2 + |I_2|^2 (a_2'')^2(b_2'')^2 \geq |I_1| (a_1'')^2(b_1'')^2 + |I_2| (a_2'')^2 (b_2'')^2 = r_3 \\
\end{align*}
Consequently, applying the case $n=2$ of the inequality to the vectors $a',b'$ yields
$$
1 + c^2 \geq r_1 r_2 + r_3 = f(a,b).
$$
Thus $f \leq 1+c^2$ on all of $\Sigma$.
