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$.
