Note that Petersen's Prop 15 does not directly address your question.
The assumption that there already exists an embedding from the $n$ dimensional $M$ to $(n+1)$ dimensional Euclidean space is crucial. Petersen's Proposition 15 concerns not solvability of $\mathrm{Riem} = S\star S$ but whether the solution $S$ is unique.
What you are looking for is a theorem due to Jaak Vilms. He described it in two different places:
- Local isometric imbedding of Riemannian $n$-manifolds into Euclidean $(n+1)$-space, J. Differential Geom. 12(2): 197-202 (1977). DOI: 10.4310/jdg/1214433981
- Factorization of Curvature Operators, Transactions of the American Mathematical Society Vol. 260, No. 2 (Aug., 1980), pp. 595-605. https://www.jstor.org/stable/1998025
Note that this discussion presumes that you have a given Riemannian metric (or I guess in your case a positive definite inner product on $T_pM$).
Essentially: if one looks at the Gauss-Codazzi equations, one sees that when $M$ has an isometric immersion into Euclidean space of $p$ dimensions higher, then one has a decomposition of
$$ \mathrm{Riem} = \sum_{i = 1}^p S_i \star S_i $$
And so the question you are asking for is identical the (local) isometric immersivity of $M^n$ into $E^{n+1}$.
Notation
Denote by $\mathscr{R}: \Lambda^2 T_pM \to \Lambda^2 T_pM$ the curvature operator acting on two forms. In indices its components are $R^{ij}_{kl}$ antisymmetric in $i,j$ and in $k,l$ separately.
Define
$$ \phi(\mathscr{R}) = R^{ij}_{kl} R^{kp}_{iq} R^{lq}_{jp} $$
and
$$ \psi(\mathscr{R}) = R^{ij}_{kl} R^{kl}_{pq} R^{pq}_{ij} $$
The curvature operator $\mathscr{R}$ is assumed to satisfy Bianchi identity.
The curvature operator $\mathscr{R}$ is said to "preserve decomposability" if for arbitrary 2 forms $\alpha, \beta$:
$$ \mathscr{R}\alpha \wedge \mathscr{R}\beta = 0 \iff \alpha \wedge \beta = 0 $$
Note that this condition is trivial if $n = 3$.
Given a curvature operator $\mathscr{R}$, its rank subspace is the smallest subspace $W$ of $T_pM$ such that the image of $\mathscr{R}$ lies in $\Lambda^2 W$; below when we refer to rank of $\mathscr{R}$ we mean the dimension of $W$.
Theorem
A curvature operator with rank $> 1$ decomposes as $S\star S$ if and only if
- $\mathcal{R}$ preserves decomposability
- $\phi + \frac14 \psi > 0$; in the case $n = 3 \pmod 4$ this condition can be replaced by positivity of $\det \mathscr{R}$ when restricted to the complement of its kernel.
- a technical "non singularity condition".
The non-singularity condition is a bit annoying to state; see the second paper for details. It is trivially satisfied when $\mathscr{R}$ is positive (or negative) definite, and $n$ is either 3 or $>4$. (I think it may also be satisfied when $n = 4$ with definite curvature operator, but the computations are less obvious.)