[proofplan]
The proof factors $\partial_{x_i}\partial_{x_j}$ through $\Delta$ on the Fourier side: since $\widehat{\partial_{x_i}\partial_{x_j} u}(\xi) = -\xi_i\xi_j\,\hat u(\xi)$ and $\widehat{\Delta u}(\xi) = -|\xi|^2\,\hat u(\xi)$, we have $\widehat{\partial_{x_i}\partial_{x_j}u}(\xi) = m_{ij}(\xi)\,\widehat{\Delta u}(\xi)$ for the multiplier $m_{ij}(\xi) := \xi_i\xi_j/|\xi|^2$. The strategy is to verify that $m_{ij}$ satisfies the Mihlin condition — bounded together with all derivatives up to order $\lfloor n/2\rfloor + 1$ scaled appropriately — and then invoke the [Mihlin Multiplier Theorem](/theorems/???) to get $L^p$ boundedness for $1 < p < \infty$. The factorisation is established first on the dense subspace $\mathcal{S}(\mathbb{R}^n)$ via the Fourier inversion formula and extended to all of $W^{2,p}(\mathbb{R}^n)$ by density.
[/proofplan]
[step:Define the multiplier $m_{ij}$ and check it is smooth and homogeneous of degree zero away from the origin]
For $i, j \in \{1, \ldots, n\}$, define
\begin{align*}
m_{ij} : \mathbb{R}^n \setminus \{0\} &\to \mathbb{R} \\
\xi &\mapsto \frac{\xi_i\,\xi_j}{|\xi|^2}.
\end{align*}
This is well-defined and smooth on $\mathbb{R}^n_0 := \mathbb{R}^n \setminus \{0\}$ since $|\xi|^2 > 0$ there. Extend $m_{ij}$ to $\mathbb{R}^n$ by setting $m_{ij}(0) := 0$; the extension is bounded but not continuous at the origin. By Cauchy--Schwarz,
\begin{align*}
|m_{ij}(\xi)| = \frac{|\xi_i|\,|\xi_j|}{|\xi|^2} \le \frac{|\xi|^2}{|\xi|^2} = 1 \qquad \text{for all } \xi \ne 0,
\end{align*}
hence $\|m_{ij}\|_{L^\infty(\mathbb{R}^n)} \le 1$.
The function $m_{ij}$ is homogeneous of degree zero: for $\lambda > 0$,
\begin{align*}
m_{ij}(\lambda \xi) = \frac{(\lambda\xi_i)(\lambda\xi_j)}{|\lambda\xi|^2} = \frac{\lambda^2\,\xi_i\xi_j}{\lambda^2\,|\xi|^2} = m_{ij}(\xi).
\end{align*}
Homogeneity of degree zero plus smoothness on $\mathbb{R}^n_0$ are precisely the conditions that yield the Mihlin estimates we now establish.
[/step]
[step:Verify the Mihlin derivative condition $|\xi|^{|\alpha|}\,|D^\alpha m_{ij}(\xi)| \le C_{\alpha,n}$]
Let $\alpha \in \mathbb{N}_0^n$ be a multi-index. Since $m_{ij}$ is homogeneous of degree zero on $\mathbb{R}^n_0$, the function $D^\alpha m_{ij}$ is homogeneous of degree $-|\alpha|$ on $\mathbb{R}^n_0$:
\begin{align*}
D^\alpha m_{ij}(\lambda\xi) = \lambda^{-|\alpha|}\,D^\alpha m_{ij}(\xi), \qquad \lambda > 0,\,\xi \ne 0.
\end{align*}
Let $\mathbb{S}^{n-1} := \{\omega \in \mathbb{R}^n : |\omega| = 1\}$ denote the unit sphere. Define
\begin{align*}
M_\alpha := \sup_{\omega \in \mathbb{S}^{n-1}} |D^\alpha m_{ij}(\omega)|.
\end{align*}
Since $m_{ij}$ is smooth on the closed set $\mathbb{S}^{n-1}$ (which lies inside $\mathbb{R}^n_0$) and $\mathbb{S}^{n-1}$ is compact, $M_\alpha < \infty$, with $M_\alpha$ depending only on $\alpha$ and $n$.
For any $\xi \ne 0$, write $\xi = |\xi|\,\omega$ with $\omega := \xi/|\xi| \in \mathbb{S}^{n-1}$. By homogeneity,
\begin{align*}
|D^\alpha m_{ij}(\xi)| = |\xi|^{-|\alpha|}\,|D^\alpha m_{ij}(\omega)| \le M_\alpha\,|\xi|^{-|\alpha|},
\end{align*}
which gives
\begin{align*}
|\xi|^{|\alpha|}\,|D^\alpha m_{ij}(\xi)| \le M_\alpha \qquad \text{for all } \xi \ne 0.
\end{align*}
This is the [Mihlin Multiplier](/page/Mihlin%20Multiplier) condition: $m_{ij}$ satisfies the Mihlin bound with constant $A_{n} := \max_{|\alpha| \le \lfloor n/2 \rfloor + 1} M_\alpha$, which depends only on $n$.
[/step]
[step:Invoke the Mihlin Multiplier Theorem to get $L^p$ boundedness of $T_{m_{ij}}$]
Define the Fourier multiplier operator
\begin{align*}
T_{m_{ij}} : \mathcal{S}(\mathbb{R}^n) &\to \mathcal{S}'(\mathbb{R}^n) \\
g &\mapsto \mathcal{F}^{-1}\bigl(m_{ij}\,\hat g\,\bigr).
\end{align*}
The [Mihlin Multiplier Theorem](/theorems/???) states: if a measurable function $m \in L^\infty(\mathbb{R}^n)$ is smooth on $\mathbb{R}^n_0$ and satisfies the bound $|\xi|^{|\alpha|}\,|D^\alpha m(\xi)| \le A$ for all multi-indices $|\alpha| \le \lfloor n/2 \rfloor + 1$, then $T_m$ extends to a bounded operator on $L^p(\mathbb{R}^n)$ for every $1 < p < \infty$, with $\|T_m\|_{\mathcal{L}(L^p(\mathbb{R}^n))} \le C_{n,p}\,A$.
We verify the hypotheses for $m = m_{ij}$: (i) $m_{ij} \in L^\infty(\mathbb{R}^n)$ with $\|m_{ij}\|_{L^\infty} \le 1$ (Step 1); (ii) $m_{ij}$ is smooth on $\mathbb{R}^n_0$ (Step 1); (iii) $|\xi|^{|\alpha|}\,|D^\alpha m_{ij}(\xi)| \le M_\alpha$ for all multi-indices $\alpha$, hence in particular for $|\alpha| \le \lfloor n/2 \rfloor + 1$ with constant $A = A_n$ (Step 2). The theorem applies and yields
\begin{align*}
\|T_{m_{ij}}\|_{\mathcal{L}(L^p(\mathbb{R}^n))} \le C_{n,p}\,A_n =: C'_{n,p},
\end{align*}
where $C'_{n,p}$ depends only on $n$ and $p$.
[/step]
[step:Establish the factorisation $\partial_{x_i}\partial_{x_j}u = T_{m_{ij}}(\Delta u)$ for Schwartz functions]
Let $u \in \mathcal{S}(\mathbb{R}^n)$. By standard properties of the Fourier transform on Schwartz functions,
\begin{align*}
\widehat{\partial_{x_i}\partial_{x_j} u}(\xi) &= (i\xi_i)(i\xi_j)\,\hat u(\xi) = -\xi_i\xi_j\,\hat u(\xi), \\
\widehat{\Delta u}(\xi) &= \sum_{k=1}^n (i\xi_k)^2\,\hat u(\xi) = -|\xi|^2\,\hat u(\xi).
\end{align*}
For $\xi \ne 0$, the second identity gives $\hat u(\xi) = -|\xi|^{-2}\,\widehat{\Delta u}(\xi)$. Substituting into the first,
\begin{align*}
\widehat{\partial_{x_i}\partial_{x_j} u}(\xi) = -\xi_i\xi_j\,\bigl(-|\xi|^{-2}\,\widehat{\Delta u}(\xi)\bigr) = \frac{\xi_i\xi_j}{|\xi|^2}\,\widehat{\Delta u}(\xi) = m_{ij}(\xi)\,\widehat{\Delta u}(\xi).
\end{align*}
The single point $\xi = 0$ contributes a $\mathcal{L}^n$-null set and so the identity holds $\mathcal{L}^n$-a.e. Applying the inverse Fourier transform (justified because both sides are tempered distributions, indeed Schwartz functions, on whose Fourier transforms the inversion formula holds),
\begin{align*}
\partial_{x_i}\partial_{x_j} u = \mathcal{F}^{-1}\bigl(m_{ij}\,\widehat{\Delta u}\bigr) = T_{m_{ij}}(\Delta u).
\end{align*}
[/step]
[step:Combine to get the inequality on Schwartz functions and extend to $W^{2,p}(\mathbb{R}^n)$ by density]
For $u \in \mathcal{S}(\mathbb{R}^n)$, combining the factorisation from Step 4 with the boundedness from Step 3,
\begin{align*}
\|\partial_{x_i}\partial_{x_j} u\|_{L^p(\mathbb{R}^n)} = \|T_{m_{ij}}(\Delta u)\|_{L^p(\mathbb{R}^n)} \le C'_{n,p}\,\|\Delta u\|_{L^p(\mathbb{R}^n)}.
\end{align*}
To extend to $u \in W^{2,p}(\mathbb{R}^n)$, recall that $\mathcal{S}(\mathbb{R}^n) \subset C^\infty_c(\mathbb{R}^n) + \mathcal{S}(\mathbb{R}^n)$ is dense in $W^{2,p}(\mathbb{R}^n)$ for $1 \le p < \infty$ (this is standard; one uses approximation by mollification combined with cutoff functions, valid because $C^\infty_c(\mathbb{R}^n)$ is dense in $W^{k,p}(\mathbb{R}^n)$ for any $k \ge 0$, $1 \le p < \infty$). Choose a sequence $(u_m)_{m \in \mathbb{N}} \subset \mathcal{S}(\mathbb{R}^n)$ with $u_m \to u$ in $W^{2,p}(\mathbb{R}^n)$, that is, $\partial_{x_i}\partial_{x_j} u_m \to \partial_{x_i}\partial_{x_j} u$ in $L^p$ and $\Delta u_m \to \Delta u$ in $L^p$.
Applying the inequality just proved on Schwartz functions to $u_m - u_\ell$ for $m, \ell \in \mathbb{N}$,
\begin{align*}
\|\partial_{x_i}\partial_{x_j} u_m - \partial_{x_i}\partial_{x_j} u_\ell\|_{L^p} \le C'_{n,p}\,\|\Delta u_m - \Delta u_\ell\|_{L^p}.
\end{align*}
The right-hand side tends to zero as $m, \ell \to \infty$, so $(\partial_{x_i}\partial_{x_j} u_m)_{m \in \mathbb{N}}$ is Cauchy in $L^p(\mathbb{R}^n)$ and hence converges; the limit must coincide with $\partial_{x_i}\partial_{x_j} u$ since $u_m \to u$ in $W^{2,p}$. Passing to the limit on both sides of the Schwartz inequality,
\begin{align*}
\|\partial_{x_i}\partial_{x_j} u\|_{L^p(\mathbb{R}^n)} = \lim_{m \to \infty}\|\partial_{x_i}\partial_{x_j} u_m\|_{L^p} \le C'_{n,p}\,\lim_{m \to \infty}\|\Delta u_m\|_{L^p} = C'_{n,p}\,\|\Delta u\|_{L^p(\mathbb{R}^n)}.
\end{align*}
Setting $C_{n,p} := C'_{n,p}$ completes the proof.
[/step]