Let $a,b,c\in\mathbb{O}$ be octonions and consider the linear map $L:\mathbb{O}\to\mathbb{O}$ defined by
$$
L(x) = (b(cx))a = R_aL_bL_c(x).
$$
One desires a formula for the characteristic polynomial of $S$, the symmetric part of $L$, i.e.,
$$
S(x) = \tfrac12\bigl(R_aL_bL_c + {}^t(R_aL_bL_c) \bigr).
$$
(I note that the OP seems to have inadvertently omitted the factor of $\tfrac12$ in the formula for $S$; this becomes apparent when one compares the claimed formula for the characteristic polynomial of $S$ when $a = b = c = 1$.)
This is equivalent to knowing the symmetric functions of the eigenvalues of the quadratic form
$$
Q(x) = L(x)\cdot x = (b(cx))a\cdot x
$$
relative to the quadratic form $Q_0(x) = x\cdot x$.
First, note that, if any of $a$, $b$ or $c$ vanishes, then, of course,
$L$ and $Q$ vanish identically, and all of the eigenvalues of $S$ are equal to zero. Thus, we can assume that none of $a$, $b$, or $c$ vanishes. Then, dividing by
$|abc|\not=0$, we can assume that $|a| = |b| = |c| = 1$.
In this case, since $L$, being a product of elements of $\mathrm{SO}(8)$, belongs to $\mathrm{SO}(8)$, it follows that $L$ is conjugate in $\mathrm{SO}(8)$ to an element of the maximal torus ${\mathrm{SO}(2)}^4$, i.e., a blocked diagonal matrix
where the diagonal elements are the $2$-by-$2$ rotation matrices $R(\theta_i)$ for $1\le i\le 4$. The matrix $S$ (the symmetric part of $L$) is then diagonal with double eigenvalues $\cos(\theta_i)$, and hence the characteristic polynomial of $S$ is
$$
p(t) = (t-\cos\theta_1)^2(t-\cos\theta_2)^2(t-\cos\theta_3)^2(t-\cos\theta_4)^2.
$$
When $p(t)= t^8 + r_1\,t^7 + r_2\,t^6 + \cdots + r_8$, we have $p(t) = q(t)^2$ where
$$
q(t) = t^4 + \tfrac12\,r_1\,t^3 + \tfrac18(4r_2-{r_1}^2)\,t^2 + \cdots
$$
(I leave it to the interested reader to work out the formulae for the $t$ and constant coefficients of $q$ as polynomials in $r_1,r_2,r_3,r_4$).
Now, by a theorem of Dickson, $b$ and $c$ lie in a quaternion subalgebra $\mathbb{A}\subset\mathbb{O}$. Also, we can write $a = \cos\theta\,a_0 + \sin\theta\, u$ where $a_0\in\mathbb{A}$ is a unit vector, $u\in\mathbb{A}^\perp$ is a unit octonion, and $0\le \theta\le \tfrac12\pi$. [This expression for $a$ will be unique relative to $\mathbb{A}$ if $0<\theta<\tfrac12\pi$.]
In what follows, it will be useful to remember that elements of $\mathbb{O}$ satisfy $xy\cdot z = x\cdot z\bar y = y\cdot \bar x z$. Recalling that $\mathrm{Re}(x) = x\cdot\mathbf{1}$ (where $\mathbf{1}\in\mathbb{O}$ is the multiplicative unit), these identities imply that $\mathrm{Re}(xy) = \mathrm{Re}(yx)$ and that $\mathrm{Re}\bigl(a(bc)\bigr) = \mathrm{Re}\bigl((ab)c\bigr)$, so that $\mathrm{Re}(abc)$ is unambiguous, even though $\mathbb{O}$ is not associative. We also have
$$
\mathrm{Re}(abc) = \mathrm{Re}(bca) = \mathrm{Re}(cab)
= \mathrm{Re}(\bar c\,\bar b\,\bar a),
$$
but note that $\mathrm{Re}(abc)\not=\mathrm{Re}(bac)$ in general.
While $\mathbb{O} = \mathbb{A}\oplus \mathbb{A} u$
(note the orthogonal direct sum) is not associative,
we have the product formula of Cayley and Dickson:
$$
(a+b\,u)(c+d\,u) = \bigl( ac - \bar d b\bigr) + (da + b\bar c)\,u
$$
for all $a,b,c,d\in\mathbb{A}$.
Writing $x\in\mathbb{O}$ as $x = x_0 + x_1\,u$ where $x_i\in\mathbb{A}$
and using the Cayley-Dickson formula several times, we have
$$
\begin{aligned}
Q(x)
&= (b(cx))a\cdot x = b(cx)\cdot x\bar a \\
&= b(c(x_0+x_1\,u))\cdot (x_0+x_1\,u)(\cos\theta\,\overline{a_0} -\sin\theta\,u)\\
&= \bigl(bcx_0 + (x_1cb)\,u\bigr)\cdot\bigl((\cos\theta\,x_0\overline{a_0}+\sin\theta\,x_1) + (\cos\theta\,x_1a_0-\sin\theta\,x_0)\,u\bigr)\\
&=\cos\theta\,bcx_0\cdot x_0\overline{a_0}
+ \sin\theta(bcx_0\cdot x_1-x_1cb\cdot x_0) + \cos\theta\,x_1cb\cdot x_1a_0\\
&= \cos\theta\,bcx_0\cdot x_0\overline{a_0}
+ \sin\theta(bcx_0-x_0\bar b\bar c)\cdot x_1
+ \cos\theta\,x_1\cdot x_1a_0\bar b\bar c\\
&= \cos\theta\,bcx_0a_0\cdot x_0
+ \sin\theta(bcx_0-x_0\bar b\bar c)\cdot x_1
+ x_1{\cdot}x_1\,\mathrm{Re}(\cos\theta\,a_0\bar b\bar c)
\end{aligned}
$$
(Note that I have used $x_1\cdot x_1a_0\bar b\bar c = \overline{x_1}x_1\cdot a_0\bar b\bar c = |x_1|^2 \mathrm{Re}(a_0\bar b\bar c)$.
This expression can be further simplified. Since $b$ and $c$ are unit vectors in
$\mathbb{A}$, we can write $bc = w$, which implies that $b = w\bar c$ and hence that $\bar b\bar c = c\bar w\bar c$. Making the substitution $x_0 = y_0\bar c$
and $x_1 = y_1\bar c$ then yields $Q_0(x) = |x_0|^2+|x_1|^2 = |y_0|^2+|y_1|^2$ while setting $v = \bar c a_0 c$ yields
$$
\begin{aligned}
Q(x) &= \cos\theta\,wy_0{\bar c}a_0 \cdot y_0\bar c
+ \sin\theta\,(wy_0-y_0\bar w)\cdot y_1
+ |y_1|^2\,\mathrm{Re}(\cos\theta\,a_0\bar b\bar c)\\
& = \cos\theta\,wy_0v \cdot y_0
+ \sin\theta\,(wy_0-y_0\bar w)\cdot y_1
+ |y_1|^2\,\mathrm{Re}(\cos\theta\,v\bar w)
\end{aligned}
$$
Using this reduced form, it is straightforward to compute that the characteristic polynomial of $S$ is
$$
\bigl(t-\mathrm{Re}(\cos\theta\,v\bar w)\bigr)^4
\bigl(t^2-2\cos\theta\,\mathrm{Re}(v)\mathrm{Re}(w)\,t
+ \cos^2\theta\,\mathrm{Re}(v)^2w\bar w - |\mathrm{Im}(w)|^2\bigr)^2,
$$
which can also be written as
$$
\bigl(t-\mathrm{Re}(\cos\theta\,v\bar w)\bigr)^4
\bigl((t-\cos\theta\,\mathrm{Re}(v)\mathrm{Re}(w))^2
- |\mathrm{Im}(w)|^2(1-\cos^2\theta\,\mathrm{Re}(v)^2)\bigr)^2,
$$
so that its roots are $t = \mathrm{Re}(\cos\theta\,v\bar w)$ with multiplicity $4$
and
$$
t = \cos\theta\,\mathrm{Re}(v)\mathrm{Re}(w)
\pm |\mathrm{Im}(w)|\bigl(1-\cos^2\theta\,\mathrm{Re}(v)^2)^{1/2},
$$
each with multiplicity $2$.
Finally, tracing back through the definitions and normalizations, we see that the characteristic polynomial of $S$ can be written in the form
$$
\bigl(t-\mathrm{Re}(a\bar b\bar c)\bigr)^4
\bigl(\bigl(t-\mathrm{Re}(a)\,\mathrm{Re}(bc)\bigr)^2-\mathrm{Im}(a)^2\,\mathrm{Im}(bc)^2\bigr)^2.
$$
Note on a computation: Since the OP asked, here is a bit of
detail of the computation of the characteristic polynomial.
We can choose a basis of $\mathbb{A}\simeq\mathbb{H}$ as $(\bf{1},\bf{i},\bf{j},\bf{k})$
in such a way that $w = w_0\,{\bf{1}}+w_1\,\bf{i}$ while $v = v_0\,{\bf{1}}+v_1\,{\bf{i}}+v_2\,{\bf{j}}+v_3\,{\bf{k}}$. Now write out $y_0$
in the orthonormal basis $(\bf{1},\bf{i},\bf{j},\bf{k})$ and $y_1$
in the orthonormal basis $(-\bf{i},\bf{1},\bf{j},\bf{k})$.
Set $\lambda_1 = c(v_{{0}}w_{{0}}{-}v_{{1}}w_{{1}})$ and $\lambda_2=c(v_{{0}}w_{{0}}{+}v_{{1}}w_{{1}})$ where $c = \cos\theta$ and $s=\sin\theta$
(to make the matrix below more readable). Then the matrix of
the quadratic form $Q$ in this basis is
$$
\left(
\begin{array}{cccccccc}
\lambda_1&0&-cw_{{1}}v_{{3}}&cw_{{1}}v_{{2}}&sw_{{1}}&0&0&0\\
0&\lambda_1&-cw_{{1}}v_{{2}}&-cw_{{1}}v_{{3}}&0&sw_{{1}}&0&0\\
-cw_{{1}}v_{{3}}&-cw_{{1}}v_{{2}}&\lambda_2&0&0&0&0&0\\
cw_{{1}}v_{{2}}&-cw_{{1}}v_{{3}}&0&\lambda_2&0&0&0&0\\
sw_{{1}}&0&0&0&\lambda_2&0&0&0\\
0&sw_{{1}}&0&0&0&\lambda_2&0&0\\
0&0&0&0&0&0&\lambda_2&0\\
0&0&0&0&0&0&0&\lambda_2
\end{array}
\right)
$$
Now tell MAPLE to compute the characteristic polynomial of this matrix and tell it to factor the result. (I confess that it's a bit surprising that $\lambda_2$ turns out to be a root of multiplicity $4$ instead of just $2$. I don't have a 'theoretical' understanding of this. Looking at the matrix, you can see that by choosing the basis of $\mathbb{A}$ a bit more carefully, you could arrange that $v_2=0$, and that makes it clear that the characteristic polynomial will be a square, but we knew that already.)
