Under what conditions is $(x^TAx)(x^TBx)The question is more or less in the title, but I was surprised (as there is no intuition in the context of the physics to which I am applying this) to be able to come up with examples of real vectors $x$ such that $x^Tx=1$ for which $(x^TAx)(x^TBx)>Tr(AB)$. Here, A and B are real-valued covariance matrices, so positive semidefinite and symmetric. A simple example is $$A=\begin{pmatrix}1 & 0.9\\\ 0.9 & 1\end{pmatrix},B=\begin{pmatrix}1 & -0.9\\\ -0.9 & 1\end{pmatrix},x=\begin{pmatrix}1 \\\ 0 \end{pmatrix}$$
Change the $0.9$s to $0.5$s and we have the reverse inequality.
My question is, is there a condition on when this is or isn't true?
Most avenues I've tried have been non-starters, so any pointers on getting this off the ground would be much appreciated.
 A: Not an answer, more of an extended comment.  It makes sense to understand what would make the left hand side large or the right hand side small.  By the spectral theorem and using positive definiteness, we may orthogonally diagonalize, so $A=UD^2U^T, B=VS^2V^T$.  Then we can write
$$\begin{split}Tr(AB)&=Tr(UD^2U^TVSV^T)\\&=Tr((DU^TVS)(SV^TUD))\\&=Tr((DU^TVS)(DU^TVS)^T)\\&=\|DU^TVS\|_F^2\end{split}$$
Using the rearrangement inequality together with the fact that every doubly stochastic matrix is in the convex hull of the permutation matrices tells us that this quantity will between $\sum \lambda_i \mu_i$ and $\sum \lambda_i \mu_{n-i}$, where $\lambda_i, \mu_i$ are the eigenvalues of $A$ and $B$ in sorted order.
The left hand side, on the other hand, is a bit tricker.  If $x$ is a common eigenvector, then the left hand side will be the product of the corresponding eigenvalues.  It will be an individual term of the sum you get calculating the Frobeninus norm of the matrix above.  But in general?  I have no clue.
I suspect that coming up with a characterization of this is difficult even when $A$ and $B$ are diagonal.

By using Lagrange multipliers, we can attempt to maximize the left hand side.  If $x$ is a maximizer, then one can show that
$$\frac{Ax}{x^TAx}+\frac{Bx}{x^TBx}=2x$$
and so $x$ is an eigenvector of a linear combination of $A$ and $B$.
In the case that $A$ and $B$ are diagonal, with eigenvalues $\lambda_1, \ldots, \lambda_n$ and $\mu_1, \ldots, \mu_n$ (in order down the diagonal, not ordered by size), this yields not only all the standard basis vectors (which our above computation shows will not violate the inequality), but additional vectors as well.  In particular, if $S\subset \{1, \ldots, n\}$ and $\alpha, \beta$ are positive constants such $\alpha\lambda_i+\beta\mu_i=2$ for all $i\in S$, then $\alpha A+\beta B$ will have an eigenvalue of $2$ with an eigenvector in $\operatorname{span}\{e_i\mid i\in S\}$, and so it is necessary to search for critical points with $x^TAx=1/\alpha, x^TBx=1/\beta$.

Let us consider $2\times 2$ diagonal matrices using this framework.
Since a common eigenvector for $A$ and $B$ is also an eigenvector for $\alpha A + \beta B$, the only way we can get more eigenvectors to check is if $\alpha A + \beta B$ has a repeated eigenvalue, in which case it will be a scalar matrix.  Diagonalizing, we may assume that $A$ and $B$ are actually diagonal.
If $A=\operatorname{diag}(a_1, a_2), B=\operatorname{diag}(b_1,b_2)$ with $a_1>a_2$ and $b_1>b_2$, then no linear combination of $A$ and $B$ with positive coefficients can be scalar, so we do not need to worry about the inequality ever being violated.
Otherwise, for simplicity, let us write $A=u_1I+u_2J, B=v_1I+v_2J$, with $I$ the identity matrix and $J=\operatorname{diag}{1,-1}$.  We are attempting to solve the following system of equations:
$$\alpha A + \beta B = 2I, x^TAx=1/\alpha, x^TBx=1/\beta, x^Tx=1.$$
If we rewrite $A, B$ in terms of $I, J$, this becomes
$$\alpha u_1+\beta v_1=2, \alpha u_2 +\beta v_2=0, x_1^2+x_2^2=1, x_1^2-x_2^2=(1 - \alpha u_1)/(\alpha u_2)=(1 - \beta v_1)/(\beta v_2).$$
When the first two equalities hold, the last out automatically does, with a value of $\frac{-(u_1v2+v_1u_2)}{2u_2v_2}$.  The system is solvable as long as this quantity is between $-1$ and $1$.  In this case, the LHS has a potential max value of $1/(\alpha\beta)=\frac{(u_1v_2-u_2v_1)^2}{-4u_2v_2}$.  This must be compared to the right hand side, $(u_1+u_2)(v_1+v_2)+(u_1-u_2)(v_1-v_2)$ to see if the desired inequality can be violated.
