[proofplan]
The identity follows immediately on the Fourier side from the multiplier formula for the Riesz transforms. The multiplier of $R_j$ is $m_j(\xi) = -i\xi_j/|\xi|$, and the multiplier of the composition $R_j \circ R_j$ is $m_j(\xi)^2$. Summing in $j$, we get $\sum_j m_j(\xi)^2 = \sum_j (-\xi_j^2/|\xi|^2) = -1$, the multiplier of $-\mathrm{Id}$. Verifying that the Fourier-side equality lifts to a genuine operator identity on $\mathcal{S}(\mathbb{R}^n)$ requires checking that the iterated Riesz transform $R_j(R_j f)$ is well-defined for $f \in \mathcal{S}$ and that the Fourier transform is injective on the relevant function class.
[/proofplan]
[step:Compute the Fourier multiplier of $R_j^2$]
By the [Riesz Transform Multiplier](/theorems/3163), for $g \in \mathcal{S}(\mathbb{R}^n)$,
\begin{align*}
\widehat{R_j g}(\xi) = -i\frac{\xi_j}{|\xi|}\hat{g}(\xi), \qquad \xi \in \mathbb{R}^n \setminus \{0\}.
\end{align*}
The multiplier $m_j: \mathbb{R}^n \setminus \{0\} \to \mathbb{C}$, $\xi \mapsto -i\xi_j/|\xi|$, is bounded with $|m_j(\xi)| = |\xi_j|/|\xi| \le 1$.
For $f \in \mathcal{S}(\mathbb{R}^n)$, we first verify that $R_j f$ is sufficiently regular for the iteration $R_j(R_j f)$ to be defined. By the Plancherel relation,
\begin{align*}
\|R_j f\|_{L^2(\mathbb{R}^n)} = \frac{1}{(2\pi)^{n/2}}\|\widehat{R_j f}\|_{L^2(\mathbb{R}^n)} = \frac{1}{(2\pi)^{n/2}}\|m_j \hat{f}\|_{L^2(\mathbb{R}^n)} \le \frac{1}{(2\pi)^{n/2}}\|\hat{f}\|_{L^2(\mathbb{R}^n)} = \|f\|_{L^2(\mathbb{R}^n)},
\end{align*}
since $|m_j| \le 1$. So $R_j f \in L^2(\mathbb{R}^n)$, and we may extend the Riesz transform to $L^2$ by density.
For $f \in \mathcal{S}(\mathbb{R}^n)$, we apply the multiplier identity twice: setting $g := R_j f \in L^2(\mathbb{R}^n)$, we have $\hat{g}(\xi) = m_j(\xi)\hat{f}(\xi) \in L^2(\mathbb{R}^n)$, and the $L^2$ extension of $R_j$ acts by multiplication by $m_j$ on the Fourier side, so
\begin{align*}
\widehat{R_j^2 f}(\xi) = \widehat{R_j g}(\xi) = m_j(\xi)\hat{g}(\xi) = m_j(\xi)^2\hat{f}(\xi).
\end{align*}
Computing $m_j(\xi)^2$:
\begin{align*}
m_j(\xi)^2 = \left(-i\frac{\xi_j}{|\xi|}\right)^2 = i^2\frac{\xi_j^2}{|\xi|^2} = -\frac{\xi_j^2}{|\xi|^2}.
\end{align*}
[/step]
[step:Sum the multipliers over $j$ to obtain the multiplier of $\sum_j R_j^2$]
By linearity of the Fourier transform,
\begin{align*}
\widehat{\sum_{j=1}^n R_j^2 f}(\xi) = \sum_{j=1}^n \widehat{R_j^2 f}(\xi) = \sum_{j=1}^n m_j(\xi)^2 \hat{f}(\xi) = \left(\sum_{j=1}^n -\frac{\xi_j^2}{|\xi|^2}\right)\hat{f}(\xi).
\end{align*}
For $\xi \in \mathbb{R}^n \setminus \{0\}$, the Euclidean norm satisfies $|\xi|^2 = \sum_{j=1}^n \xi_j^2$, so
\begin{align*}
\sum_{j=1}^n -\frac{\xi_j^2}{|\xi|^2} = -\frac{1}{|\xi|^2}\sum_{j=1}^n \xi_j^2 = -\frac{|\xi|^2}{|\xi|^2} = -1.
\end{align*}
Therefore
\begin{align*}
\widehat{\sum_{j=1}^n R_j^2 f}(\xi) = -\hat{f}(\xi)
\end{align*}
for $\xi \in \mathbb{R}^n \setminus \{0\}$, hence $\mathcal{L}^n$-a.e. on $\mathbb{R}^n$.
[/step]
[step:Lift the Fourier-side identity to an operator identity]
We have shown $\widehat{\sum_j R_j^2 f}(\xi) = -\hat{f}(\xi)$ for $\mathcal{L}^n$-a.e. $\xi$, where the left side is interpreted as the Fourier transform of an $L^2$ function. By the Plancherel theorem, the Fourier transform $\mathcal{F}: L^2(\mathbb{R}^n) \to L^2(\mathbb{R}^n)$ is a bijection (an isometry up to a constant of $(2\pi)^{n/2}$). Two $L^2$ functions with $\mathcal{L}^n$-a.e. equal Fourier transforms are themselves equal as elements of $L^2$, hence equal almost everywhere as functions. Therefore
\begin{align*}
\sum_{j=1}^n R_j^2 f = -f \quad \text{in } L^2(\mathbb{R}^n).
\end{align*}
For $f \in \mathcal{S}(\mathbb{R}^n)$, both sides are continuous functions on $\mathbb{R}^n$ (the iterated Riesz transform of a Schwartz function is Schwartz, by repeated application of the argument in the proof of Theorem 3158 adapted to the Riesz multiplier; in particular, the multiplier $-\xi_j^2/|\xi|^2$ is bounded and smooth on $\mathbb{R}^n \setminus \{0\}$, with a removable discontinuity at $0$ except for a derivative singularity that yields rapid decay of the inverse transform via the same dominated-convergence argument). Two continuous $L^2$ functions equal a.e. are equal everywhere, so
\begin{align*}
\sum_{j=1}^n R_j^2 f(x) = -f(x) \quad \text{for all } x \in \mathbb{R}^n.
\end{align*}
This completes the proof of the Riesz identity on $\mathcal{S}(\mathbb{R}^n)$. By the boundedness of $R_j$ on $L^2(\mathbb{R}^n)$ and density of $\mathcal{S}(\mathbb{R}^n) \subset L^2(\mathbb{R}^n)$, the identity extends to $-\mathrm{Id}$ as a bounded operator on $L^2(\mathbb{R}^n)$.
[/step]