When can the curvature operator of a Riemannian manifold (M,g) be diagonalized by a basis of the following form
'{${E_i\wedge E_j }$}' where '{${E_i}$}' is an orthonormal basis of the tangent space? If the manifold is three dimensional then it is always possible. But what about higher dimensional cases?
A sufficient condition is if the Riemannian manifold is conformally flat, this implies that the Weyl curvature vanishes, and the Riemann curvature tensor is a linear combination of the identity operator on two forms and the operator formed by the Kulkarni-Nomizu product of the Ricci curvature and the metric. Using that the Ricci curvature is a symmetric bilinear form, you can diagonalize it relative to the metric, and explicitly show (as in the 3 dimensional case) that the Kulkarni-Nomizu product of Ricci and the metric can be diagonalized over a basis formed by $\{e_i\wedge e_j\}$.
On the other hand, there are also large classes of manifolds for which it is impossible to satisfy your requirement. For example, consider the four dimensional (anti)-self-dual Einstein manifolds with nonvanishing Weyl curvature. The Einstein equation $Ric = \lambda g$ means that the Ricci and scalar parts of the curvature are just multiplies of the identity. But the self-duality of the Weyl part means any eigen-twoform of the curvature operator must be either self-dual or anti-self-dual, which rules them out from being rank two.
Here are also some possibly relevant papers.
There are still a few interesting things to say about this question, so I thought I'd add some comments.
In one sense, the answer to the question of when a Riemannian metric has an orthonormal coframing that diagonalizes the curvature in the manner requested by the OP is an algebraic problem: As is well-known, the space of curvature tensors (when regarded as quadratic forms on $\Lambda^2$ that satisfy the first Bianchi identity) has dimension $D_n = \tfrac{1}{12}n^2(n^2{-}1)$, while the set of those diagonalizable in some orthonormal coframe is the $\mathrm{O}(n)$-orbit of a linear subspace of dimension $\tfrac{1}{2}n(n{-}1)$, so (when $n\ge 3$) it is a cone $\mathcal{R}_n$ of dimension $n(n{-}1)$. Thus, a Riemannian metric will have to satisfy a set of at least $$ R_n = \tfrac{1}{12}n^2(n^2{-}1) - n(n{-}1) = \tfrac{1}{12}n(n{-}1)(n{-}3)(n{+}4) $$ polynomial relations on its curvature in order for such a diagonalization to be possible at every point. Writing out a set of generators for these relations is not likely to be easy and probably won't be enlightening, even for $n=4$, which is when it first becomes nontrivial. Moreover, there is no guarantee that this ideal $\mathcal{I}_n$ of polynomial relations is generated by only $R_n$ polynomials or that you won't still have to impose some inequalities to make sure that the curvature is diagonalizable by an element in $\mathrm{O}(n)$ rather than by an element of $\mathrm{O}(n,\mathbb{C})$ that doesn't lie in $\mathrm{O}(n)$. However, this approach would, in theory, give the answer to the OP's literal question.
On the other hand, one might want to interpret the question as asking how one could 'generate' all of the metrics that satisfy this diagonalizability property, at least locally. This is a more interesting (and more challenging) problem. Willie and Thomas have each given examples of classes of such metrics that essentially depend on one function of $n$ variables: Willie cited the conformally flat metrics, which are locally of the form $e^u g_0$ where $g_0$ is the standard metric on $\mathbb{R}^n$, and Thomas cited the induced metrics on hypersurfaces in a space form of dimension $n{+}1$, each of which, locally, can be described as the graph of one function of $n$-variables). The interesting question is whether these are, themselves, special cases of some more general class of metrics with the desired property. Might there be a class of examples that depend on more than one arbitrary function of $n$ variables? Another interesting question is whether their examples 'reach' all the curvature tensors that satisfy the relations $\mathcal{I}_n$ and, if not, whether are there other examples that do.
This latter question is easier to answer than the former. It is easy to see, just by an algebraic count, that neither the conformally flat metrics nor those induced on hypersurfaces in space forms can actually `reach' all of the $\mathrm{O}(n)$-orbits in $\mathcal{R}_n$. (The two sets of orbits that they reach overlap, and they are distinct, proper closed subsets of $\mathcal{R}_n$.) On the other hand, examples provided by É. Cartan of nondegenerate submanifolds of dimension $n$ in $\mathbb{R}^{2n}$ that have flat normal bundle turn out to have their curvature tensors in $\mathcal{R}_n$ and, using these, one can reach an open subset of the orbits in $\mathcal{R}_n$. Now, Cartan's examples depend locally on $n^2{-}n$ arbitrary functions of two variables (not $n$ variables), and it turns out that they satisfy many more differential equations (of higher order) than just the $\mathcal{I}_n$. For example, in Cartan's examples, the diagonalizing coframing $\omega=(\omega_i)$ turns out to be integrable, i.e., $\omega_i\wedge d\omega_i = 0$ for all $i$, so that the metric itself can be diagonalized in a local coordinate chart and thus is locally of the form $$ g = e^{2f_1}\ {dx_1}^2 +e^{2f_2}\ {dx_2}^2 + \cdots + e^{2f_n}\ {dx_n}^2. $$ Meanwhile, the condition for a metric in this form to have its curvature tensor be diagonal with respect to the coframing $\omega = (\omega_i) = (e^{f_i}\ dx_i)$ and, hence, take values in $\mathcal{R}_n$ turns out to be an involutive system of second order PDE for the functions $f_i$ whose general local solution depends on $n^2{-}n$ arbitrary functions of two variables. These turn out to be slightly more general than the ones that arise as Cartan's examples, and, using solutions of this type, one can reach all of the $\mathrm{O}(n)$-orbits in $\mathcal{R}_n$.
However, the question of how to 'generate' the 'general' metric whose curvature tensor takes values in $\mathcal{R}_n$ for $n\ge 4$ seems to be a very difficult problem. It is an overdetermined system for the metric that is not involutive, and computing its first two prolongations, even in the $n=4$ case, yields a system that is extremely algebraically complicated and still not involutive. Thus, I do not know (and I believe that it is not known) whether the general local solution of this problem depends (modulo diffeomorphism) on more than one arbitrary function of $n$ variables.
This is not what you asking, but maybe you want to know.
$\mathbb C\mathrm P^2$ does not admit a metric with splitting tensor in your sense. It has a non-zero Pontryagin number and any Pontryagin number can be expressed as an integral of some function of curvature tensor which is zero for splitting tensors. (In fact such a function is zero on any tensor which has orientation reversing symmetry and clearly splinting tensors are among them.)
Any riemannian metric on $\mathbb S^n$ can be $C^0$ approximated by $C^\infty$ metrics with splitting curvature tensor. This follows from Nash--Kuiper embedding theorem and part (1) in the answer of Thomas Richard.
I am not aware of any systematic treatment of this question. However there are two standard examples.
1) If $(M,g)$ is an hyper surface in a constant curvature manifold, then its curvature operator is given by $R=A\wedge A-kI$ where $I$ is the identity and $A$ is the second fundamental form. Then if $e_i$ diagonalise $A$, $e_i\wedge e_j$ diagonalise $R$.
2) Another example is rotationnaly symmetric metrics on $\mathbb{R}\times S^n$, the proof can be found in Riemannian Geometry by Petersen.
