[proofplan]
We define $v(x) = u(Ox)$ and compute $\Delta v$ by applying the multivariable chain rule twice. The first application expresses $\partial_{x_i} v$ as a sum involving the entries $O_{ji}$ of the orthogonal matrix. The second application produces a double sum involving second-order partial derivatives of $u$ and products of entries of $O$. The orthogonality condition $O^\top O = I_n$ collapses the double sum to $\Delta u$, which vanishes by hypothesis.
[/proofplan]
[step:Define the composed function and compute its first partial derivatives via the chain rule]
Define the rotated coordinates $y = Ox$, so that $y_k = \sum_{j=1}^n O_{kj} x_j$ for each $k \in \{1, \ldots, n\}$. The composed function is $v(x) = u(y)$ where $y = Ox$. Since $O$ is a constant matrix, each $y_k$ is a linear function of $x$, and in particular
\begin{align*}
\frac{\partial y_k}{\partial x_i} = O_{ki}.
\end{align*}
Applying the [Multivariable Chain Rule](/theorems/830) to the composition $v = u \circ L_O$, where $L_O: \mathbb{R}^n \to \mathbb{R}^n$ is the linear map $x \mapsto Ox$, we obtain
\begin{align*}
\frac{\partial v}{\partial x_i} = \sum_{j=1}^n \frac{\partial u}{\partial y_j}(Ox) \cdot \frac{\partial y_j}{\partial x_i} = \sum_{j=1}^n \frac{\partial u}{\partial y_j}(Ox) \cdot O_{ji}.
\end{align*}
The chain rule applies because $u \in C^2(U)$ and $L_O$ is smooth (linear).
[guided]
We want to compute $\Delta v = \sum_{i=1}^n \frac{\partial^2 v}{\partial x_i^2}$, so we first need the first partial derivatives of $v$ with respect to the $x$-coordinates.
Define $y = Ox \in \mathbb{R}^n$, so that componentwise $y_k = \sum_{j=1}^n O_{kj} x_j$ for each $k \in \{1, \ldots, n\}$. The function $v$ is the composition $v = u \circ L_O$ where $L_O: \mathbb{R}^n \to \mathbb{R}^n$ is the linear map $x \mapsto Ox$. Since $L_O$ is linear, its Jacobian is the constant matrix $O$ itself, so
\begin{align*}
\frac{\partial y_k}{\partial x_i} = O_{ki}.
\end{align*}
Why does the [Multivariable Chain Rule](/theorems/830) apply? The chain rule requires that $L_O$ be differentiable at $x$ and that $u$ be differentiable at $L_O(x) = Ox$. The first condition holds because $L_O$ is linear (hence smooth). The second holds because $u \in C^2(U)$ and $Ox \in U$ by the hypothesis $O(U) = U$.
Applying the chain rule to $v(x) = u(Ox)$, we differentiate with respect to $x_i$. Each intermediate variable $y_j$ depends on $x_i$ via $\frac{\partial y_j}{\partial x_i} = O_{ji}$, so:
\begin{align*}
\frac{\partial v}{\partial x_i} = \sum_{j=1}^n \frac{\partial u}{\partial y_j}(Ox) \cdot \frac{\partial y_j}{\partial x_i} = \sum_{j=1}^n \frac{\partial u}{\partial y_j}(Ox) \cdot O_{ji}.
\end{align*}
Note that $O_{ji}$ is the $(j, i)$-entry of $O$, not the $(i, j)$-entry: the index $j$ runs over the output coordinates $y_j$, and $i$ is the input coordinate $x_i$ we differentiate with respect to.
[/guided]
[/step]
[step:Differentiate a second time to express $\Delta v$ as a double sum over entries of $O$]
Differentiating $\frac{\partial v}{\partial x_i} = \sum_{j=1}^n \frac{\partial u}{\partial y_j}(Ox) \cdot O_{ji}$ once more with respect to $x_i$, and noting that $O_{ji}$ is a constant, we obtain
\begin{align*}
\frac{\partial^2 v}{\partial x_i^2} = \sum_{j=1}^n O_{ji} \cdot \frac{\partial}{\partial x_i} \left( \frac{\partial u}{\partial y_j}(Ox) \right).
\end{align*}
The inner expression $\frac{\partial u}{\partial y_j}(Ox)$ is itself a composition of the function $\frac{\partial u}{\partial y_j}: U \to \mathbb{R}$ with the linear map $L_O$. Applying the [Multivariable Chain Rule](/theorems/830) a second time — valid because $u \in C^2(U)$ guarantees that $\frac{\partial u}{\partial y_j} \in C^1(U)$ — we get
\begin{align*}
\frac{\partial}{\partial x_i} \left( \frac{\partial u}{\partial y_j}(Ox) \right) = \sum_{k=1}^n \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox) \cdot O_{ki}.
\end{align*}
Substituting back:
\begin{align*}
\frac{\partial^2 v}{\partial x_i^2} = \sum_{j=1}^n \sum_{k=1}^n O_{ji} \, O_{ki} \, \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox).
\end{align*}
Summing over $i = 1, \ldots, n$ to form the Laplacian:
\begin{align*}
\Delta v = \sum_{i=1}^n \frac{\partial^2 v}{\partial x_i^2} = \sum_{j=1}^n \sum_{k=1}^n \left( \sum_{i=1}^n O_{ji} \, O_{ki} \right) \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox).
\end{align*}
[guided]
We now differentiate the expression $\frac{\partial v}{\partial x_i} = \sum_{j=1}^n \frac{\partial u}{\partial y_j}(Ox) \cdot O_{ji}$ with respect to $x_i$ a second time. Since $O_{ji}$ is a constant (the entries of the orthogonal matrix do not depend on $x$), the derivative passes through to the $u$-dependent factor:
\begin{align*}
\frac{\partial^2 v}{\partial x_i^2} = \sum_{j=1}^n O_{ji} \cdot \frac{\partial}{\partial x_i} \left( \frac{\partial u}{\partial y_j}(Ox) \right).
\end{align*}
The expression $\frac{\partial u}{\partial y_j}(Ox)$ is the composition of $\partial_{y_j} u: U \to \mathbb{R}$ with the linear map $L_O: x \mapsto Ox$. Since $u \in C^2(U)$, the first partial $\partial_{y_j} u$ belongs to $C^1(U)$, so the [Multivariable Chain Rule](/theorems/830) applies a second time. The result is:
\begin{align*}
\frac{\partial}{\partial x_i} \left( \frac{\partial u}{\partial y_j}(Ox) \right) = \sum_{k=1}^n \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox) \cdot \frac{\partial y_k}{\partial x_i} = \sum_{k=1}^n \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox) \cdot O_{ki}.
\end{align*}
Substituting into the expression for $\frac{\partial^2 v}{\partial x_i^2}$:
\begin{align*}
\frac{\partial^2 v}{\partial x_i^2} = \sum_{j=1}^n \sum_{k=1}^n O_{ji} \, O_{ki} \, \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox).
\end{align*}
To form $\Delta v$, we sum over $i = 1, \ldots, n$. The key structural observation is that we can rearrange the order of summation: the indices $j$ and $k$ are "output-side" indices (they index derivatives of $u$ with respect to $y$), while $i$ is an "input-side" index. Pulling the $i$-sum inward:
\begin{align*}
\Delta v = \sum_{i=1}^n \frac{\partial^2 v}{\partial x_i^2} = \sum_{j=1}^n \sum_{k=1}^n \underbrace{\left( \sum_{i=1}^n O_{ji} \, O_{ki} \right)}_{\text{inner product of rows } j \text{ and } k \text{ of } O} \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox).
\end{align*}
The parenthesized sum is the inner product of the $j$-th and $k$-th rows of $O$. This is where orthogonality will enter.
[/guided]
[/step]
[step:Apply the orthogonality condition $O^\top O = I_n$ to collapse the double sum to $\Delta u$]
The inner sum $\sum_{i=1}^n O_{ji} \, O_{ki}$ is the $(j, k)$-entry of the matrix product $O O^\top$. By the [Orthogonal Group and Its Properties](/theorems/805), $O$ satisfies $O^\top O = I_n$, which is equivalent to $O O^\top = I_n$ (every left inverse of a square matrix is a right inverse). Therefore
\begin{align*}
\sum_{i=1}^n O_{ji} \, O_{ki} = (O O^\top)_{jk} = (I_n)_{jk} = \delta_{jk},
\end{align*}
where $\delta_{jk}$ denotes the Kronecker delta. Substituting into the expression for $\Delta v$:
\begin{align*}
\Delta v = \sum_{j=1}^n \sum_{k=1}^n \delta_{jk} \, \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox) = \sum_{j=1}^n \frac{\partial^2 u}{\partial y_j^2}(Ox) = (\Delta u)(Ox).
\end{align*}
[guided]
We now evaluate the inner sum $\sum_{i=1}^n O_{ji} \, O_{ki}$. Writing this in matrix notation: the $(j, k)$-entry of $O O^\top$ is
\begin{align*}
(O O^\top)_{jk} = \sum_{i=1}^n O_{ji} \, (O^\top)_{ik} = \sum_{i=1}^n O_{ji} \, O_{ki}.
\end{align*}
This is precisely our inner sum.
By the [Orthogonal Group and Its Properties](/theorems/805), $O$ is orthogonal, meaning $O^\top O = I_n$. Since $O$ is a square matrix with $O^\top O = I_n$, the matrix $O^\top$ is both a left and right inverse of $O$, so $O O^\top = I_n$ as well. Therefore:
\begin{align*}
\sum_{i=1}^n O_{ji} \, O_{ki} = (O O^\top)_{jk} = (I_n)_{jk} = \delta_{jk},
\end{align*}
where $\delta_{jk}$ is the Kronecker delta ($\delta_{jk} = 1$ if $j = k$ and $\delta_{jk} = 0$ if $j \neq k$).
This is the heart of the proof: the orthogonality of $O$ forces all cross-terms (those with $j \neq k$) to vanish. Substituting:
\begin{align*}
\Delta v = \sum_{j=1}^n \sum_{k=1}^n \delta_{jk} \, \frac{\partial^2 u}{\partial y_k \partial y_j}(Ox) = \sum_{j=1}^n \frac{\partial^2 u}{\partial y_j^2}(Ox) = (\Delta u)(Ox).
\end{align*}
The double sum collapses to a single sum of pure second partials — precisely the Laplacian of $u$ evaluated at $Ox$.
[/guided]
[/step]
[step:Conclude that $\Delta v = 0$ using the hypothesis $\Delta u = 0$]
Since $u$ satisfies $\Delta u = 0$ on $U$ and $Ox \in U$ for every $x \in U$ (by the hypothesis $O(U) = U$), we have $(\Delta u)(Ox) = 0$ for every $x \in U$. Combining with the computation above:
\begin{align*}
\Delta v(x) = (\Delta u)(Ox) = 0 \quad \text{for all } x \in U.
\end{align*}
Therefore $v$ is harmonic on $U$, and [Laplace's equation](/page/Laplace's%20Equation) is invariant under orthogonal transformations.
[/step]