[proofplan]
We first show that $b$ commutes with every spectral projection of $a$. The commutation relation $ab=ba$ gives commutation with every [continuous function](/page/Continuous%20Function) of $a$ by continuous functional calculus, and then a monotone-class argument plus strong convergence in the [Borel functional calculus](/theorems/2696) extends this to the projections $E_a(A)$. Once $E_a(A)$ commutes with $b$, the same argument applied to the operator $b$ shows that $E_a(A)$ commutes with every spectral projection $E_b(B)$.
[/proofplan]
[step:Pass from commutation with $a$ to commutation with continuous functions of $a$]
Let $\sigma(a)\subset\mathbb{R}$ denote the spectrum of $a$. Since $a$ is bounded and self-adjoint, $\sigma(a)$ is compact, and the continuous functional calculus gives a unital $*$-homomorphism
\begin{align*}
C(\sigma(a)) &\to \mathcal{L}(H)
\end{align*}
sending the coordinate function $t\mapsto t$ to $a$.
Because $ab=ba$, the operator $b$ commutes with every polynomial in $a$. By the [Weierstrass approximation theorem](/theorems/480), every $f\in C(\sigma(a))$ is a uniform limit on $\sigma(a)$ of real or complex polynomials restricted to $\sigma(a)$. The continuous functional calculus is continuous for the [uniform norm](/page/Uniform%20Norm), so if $p_n\to f$ uniformly on $\sigma(a)$, then
\begin{align*}
\|p_n(a)-f(a)\|_{\mathrm{op}}\to 0.
\end{align*}
Since each $p_n(a)$ commutes with $b$, taking the operator-norm limit gives
\begin{align*}
bf(a)=f(a)b.
\end{align*}
Thus $b$ commutes with $f(a)$ for every $f\in C(\sigma(a))$.
[/step]
[step:Extend commutation from continuous functions of $a$ to all spectral projections of $a$]
Let $\mathcal{D}_a$ be the class of Borel sets $A\subseteq\mathbb{R}$ such that
\begin{align*}
bE_a(A)=E_a(A)b.
\end{align*}
We prove that $\mathcal{D}_a=\mathcal{B}(\mathbb{R})$.
First let $F\subseteq\mathbb{R}$ be closed. If $F\cap\sigma(a)=\varnothing$, then the spectral measure of $a$ is supported on $\sigma(a)$, so
\begin{align*}
E_a(F)=E_a(F\cap\sigma(a))=E_a(\varnothing)=0.
\end{align*}
Hence $bE_a(F)=0=E_a(F)b$ in this case.
It remains to treat the case $F\cap\sigma(a)\ne\varnothing$. For any nonempty subset $S\subseteq\mathbb R$, let $d_S:\mathbb R\to[0,\infty)$ denote the distance map defined by $d_S(x)=\inf_{y\in S}|x-y|$ for $x\in\mathbb R$. Define, for each $n\in\mathbb{N}$, the function
\begin{align*}
f_n:\sigma(a)&\to[0,1]
\end{align*}
by
\begin{align*}
f_n(t)=\max\{0,1-n\,d_{F\cap\sigma(a)}(t)\}.
\end{align*}
Each $f_n$ is continuous on $\sigma(a)$, and $f_n(t)\downarrow \mathbb{1}_{F\cap\sigma(a)}(t)$ for every $t\in\sigma(a)$. By the bounded monotone convergence property of the Borel functional calculus for bounded [self-adjoint operators](/page/Self-Adjoint%20Operators),
\begin{align*}
f_n(a)\to \mathbb{1}_{F\cap\sigma(a)}(a)=E_a(F)
\end{align*}
strongly. Since every $f_n(a)$ commutes with $b$ by the previous step, for every $\xi\in H$ we have
\begin{align*}
bf_n(a)\xi=f_n(a)b\xi.
\end{align*}
Passing to the strong limit and using boundedness of $b$ gives
\begin{align*}
bE_a(F)\xi=E_a(F)b\xi.
\end{align*}
Since $\xi\in H$ was arbitrary, $bE_a(F)=E_a(F)b$. Hence every [closed set](/page/Closed%20Set) belongs to $\mathcal{D}_a$.
We next prove that $\mathcal{D}_a$ is a $\sigma$-algebra. Since $E_a(\mathbb{R})=I$, the set $\mathbb{R}$ belongs to $\mathcal{D}_a$. If $A\in\mathcal{D}_a$, then the projection-valued measure identity $E_a(\mathbb{R}\setminus A)=I-E_a(A)$ gives
\begin{align*}
bE_a(\mathbb{R}\setminus A)=E_a(\mathbb{R}\setminus A)b,
\end{align*}
so $\mathbb{R}\setminus A\in\mathcal{D}_a$.
If $A,C\in\mathcal{D}_a$, then the projection-valued measure identity $E_a(A\cap C)=E_a(A)E_a(C)$ gives
\begin{align*}
bE_a(A\cap C)=bE_a(A)E_a(C)=E_a(A)bE_a(C)=E_a(A)E_a(C)b=E_a(A\cap C)b.
\end{align*}
Thus $A\cap C\in\mathcal{D}_a$, and complements then imply $A\cup C\in\mathcal{D}_a$.
Now let $(A_n)_{n\in\mathbb{N}}$ be any sequence in $\mathcal{D}_a$. For each $m\in\mathbb{N}$, define
\begin{align*}
B_m=\bigcup_{n=1}^{m}A_n.
\end{align*}
The preceding finite-union closure gives $B_m\in\mathcal{D}_a$, and $(B_m)_{m\in\mathbb{N}}$ is increasing with union $A=\bigcup_{n=1}^{\infty}A_n$. By countable additivity of the spectral measure in the strong operator topology,
\begin{align*}
E_a(B_m)\to E_a(A)
\end{align*}
strongly. Passing to the strong limit in $bE_a(B_m)=E_a(B_m)b$ gives $bE_a(A)=E_a(A)b$, so $A\in\mathcal{D}_a$. Hence $\mathcal{D}_a$ is a $\sigma$-algebra containing all closed subsets of $\mathbb{R}$. Since the Borel $\sigma$-algebra $\mathcal{B}(\mathbb{R})$ is the smallest $\sigma$-algebra containing the closed subsets of $\mathbb{R}$, we obtain
\begin{align*}
\mathcal{D}_a=\mathcal{B}(\mathbb{R}).
\end{align*}
Therefore
\begin{align*}
bE_a(A)=E_a(A)b
\end{align*}
for every Borel set $A\subseteq\mathbb{R}$.
[guided]
The goal of this step is to upgrade from continuous functions of $a$ to characteristic functions of Borel sets. This upgrade cannot be done by uniform approximation of $\mathbb{1}_A$, because characteristic functions of arbitrary Borel sets are usually discontinuous. Instead, we use the strong-operator convergence built into the Borel functional calculus.
Define $\mathcal{D}_a$ to be the collection of Borel sets $A\subseteq\mathbb{R}$ such that
\begin{align*}
bE_a(A)=E_a(A)b.
\end{align*}
We first show that closed sets belong to $\mathcal{D}_a$. Let $F\subseteq\mathbb{R}$ be closed. There is one edge case to remove before defining a distance function. If $F\cap\sigma(a)=\varnothing$, then the spectral measure of $a$ is supported on $\sigma(a)$, so
\begin{align*}
E_a(F)=E_a(F\cap\sigma(a))=E_a(\varnothing)=0.
\end{align*}
Therefore $bE_a(F)=0=E_a(F)b$, and $F\in\mathcal{D}_a$ in this case.
Assume now that $F\cap\sigma(a)\ne\varnothing$. For any nonempty subset $S\subseteq\mathbb R$, let $d_S:\mathbb R\to[0,\infty)$ denote the distance map defined by $d_S(x)=\inf_{y\in S}|x-y|$ for $x\in\mathbb R$. For each $n\in\mathbb{N}$, define
\begin{align*}
f_n:\sigma(a)&\to[0,1]
\end{align*}
by
\begin{align*}
f_n(t)=\max\{0,1-n\,d_{F\cap\sigma(a)}(t)\}.
\end{align*}
This is a continuous function on $\sigma(a)$ because the distance to the nonempty closed subset $F\cap\sigma(a)$ of the compact [metric space](/page/Metric%20Space) $\sigma(a)$ is continuous, and the map $s\mapsto \max\{0,1-ns\}$ is continuous on $[0,\infty)$. Moreover, $f_n(t)=1$ when $t\in F\cap\sigma(a)$, and for $t\notin F\cap\sigma(a)$ the positive distance from $t$ to $F\cap\sigma(a)$ eventually forces $f_n(t)=0$. Hence
\begin{align*}
f_n(t)\downarrow \mathbb{1}_{F\cap\sigma(a)}(t)
\end{align*}
for every $t\in\sigma(a)$.
The bounded monotone convergence property of the Borel functional calculus now applies: the functions $f_n$ are uniformly bounded by $1$ and decrease pointwise to $\mathbb{1}_{F\cap\sigma(a)}$. Therefore
\begin{align*}
f_n(a)\to \mathbb{1}_{F\cap\sigma(a)}(a)
\end{align*}
strongly. By the definition of the spectral projection,
\begin{align*}
\mathbb{1}_{F\cap\sigma(a)}(a)=E_a(F).
\end{align*}
From the previous step, each $f_n(a)$ commutes with $b$. Thus, for every vector $\xi\in H$,
\begin{align*}
bf_n(a)\xi=f_n(a)b\xi.
\end{align*}
Taking limits in $H$ gives
\begin{align*}
bE_a(F)\xi=E_a(F)b\xi,
\end{align*}
because $f_n(a)\xi\to E_a(F)\xi$, $f_n(a)b\xi\to E_a(F)b\xi$, and $b$ is bounded. Since this equality holds for every $\xi\in H$, we obtain
\begin{align*}
bE_a(F)=E_a(F)b.
\end{align*}
So every closed set is in $\mathcal{D}_a$.
It remains to pass from closed sets to all Borel sets. We prove the stronger fact that $\mathcal{D}_a$ is a $\sigma$-algebra. First, $\mathbb{R}\in\mathcal{D}_a$ because $E_a(\mathbb{R})=I$ and $bI=Ib$. If $A\in\mathcal{D}_a$, then
\begin{align*}
E_a(\mathbb{R}\setminus A)=I-E_a(A)
\end{align*}
by finite additivity of the projection-valued measure $E_a$. Hence
\begin{align*}
bE_a(\mathbb{R}\setminus A)=b(I-E_a(A))=(I-E_a(A))b=E_a(\mathbb{R}\setminus A)b,
\end{align*}
so $\mathbb{R}\setminus A\in\mathcal{D}_a$.
Next, if $A,C\in\mathcal{D}_a$, then the multiplicativity property of a projection-valued measure gives
\begin{align*}
E_a(A\cap C)=E_a(A)E_a(C).
\end{align*}
Using the commutation of $b$ with both factors, we compute
\begin{align*}
bE_a(A\cap C)=bE_a(A)E_a(C)=E_a(A)bE_a(C)=E_a(A)E_a(C)b=E_a(A\cap C)b.
\end{align*}
Thus $A\cap C\in\mathcal{D}_a$. Since $A\cup C=\mathbb{R}\setminus((\mathbb{R}\setminus A)\cap(\mathbb{R}\setminus C))$, closure under complements and finite intersections gives closure under finite unions.
Now let $(A_n)_{n\in\mathbb{N}}$ be any sequence in $\mathcal{D}_a$. Define the finite unions
\begin{align*}
B_m=\bigcup_{n=1}^{m}A_n
\end{align*}
for $m\in\mathbb{N}$. Each $B_m$ belongs to $\mathcal{D}_a$, and the sequence $(B_m)_{m\in\mathbb{N}}$ increases to $A=\bigcup_{n=1}^{\infty}A_n$. Countable additivity of the spectral measure in the strong operator topology gives
\begin{align*}
E_a(B_m)\to E_a(A)
\end{align*}
strongly. Passing to the strong limit in
\begin{align*}
bE_a(B_m)=E_a(B_m)b
\end{align*}
gives $bE_a(A)=E_a(A)b$, because $b$ is bounded. Therefore $A\in\mathcal{D}_a$. We have shown that $\mathcal{D}_a$ is a $\sigma$-algebra.
Since $\mathcal{D}_a$ contains all closed subsets of $\mathbb{R}$, and since $\mathcal{B}(\mathbb{R})$ is the smallest $\sigma$-algebra containing the closed subsets of $\mathbb{R}$, we conclude that
\begin{align*}
\mathcal{D}_a=\mathcal{B}(\mathbb{R}).
\end{align*}
Therefore $b$ commutes with $E_a(A)$ for every Borel set $A\subseteq\mathbb{R}$.
[/guided]
[/step]
[step:Apply the same Borel functional calculus argument to $b$]
Fix a Borel set $A\subseteq\mathbb{R}$ and define
\begin{align*}
P=E_a(A).
\end{align*}
The operator $P$ is an [orthogonal projection](/theorems/437), hence $P\in\mathcal{L}(H)$ and $P=P^*$. The preceding step gives
\begin{align*}
Pb=bP.
\end{align*}
Therefore $P$ commutes with every polynomial in $b$.
Let $\sigma(b)\subset\mathbb{R}$ denote the spectrum of $b$. Since $b$ is bounded and self-adjoint, $\sigma(b)$ is compact. If $g\in C(\sigma(b))$ and $(q_n)_{n\in\mathbb{N}}$ is a sequence of polynomials with $q_n\to g$ uniformly on $\sigma(b)$, then the continuous functional calculus gives
\begin{align*}
\|q_n(b)-g(b)\|_{\mathrm{op}}\to 0.
\end{align*}
Since $Pq_n(b)=q_n(b)P$ for every $n\in\mathbb{N}$, passing to the operator-norm limit gives
\begin{align*}
Pg(b)=g(b)P.
\end{align*}
Thus $P$ commutes with every continuous function of $b$.
Let $\mathcal{D}_b$ be the class of Borel sets $B\subseteq\mathbb{R}$ such that
\begin{align*}
PE_b(B)=E_b(B)P.
\end{align*}
We prove that $\mathcal{D}_b=\mathcal{B}(\mathbb{R})$.
First let $F\subseteq\mathbb{R}$ be closed. If $F\cap\sigma(b)=\varnothing$, then the spectral measure of $b$ is supported on $\sigma(b)$, so $E_b(F)=0$, and hence $PE_b(F)=0=E_b(F)P$. If $F\cap\sigma(b)\ne\varnothing$, let $d_{F\cap\sigma(b)}:\mathbb R\to[0,\infty)$ be the distance map defined by $d_{F\cap\sigma(b)}(x)=\inf_{y\in F\cap\sigma(b)}|x-y|$. For each $n\in\mathbb{N}$, define $h_n:\sigma(b)\to[0,1]$ by
\begin{align*}
h_n(t)=\max\{0,1-n\,d_{F\cap\sigma(b)}(t)\}.
\end{align*}
Then $h_n$ is continuous on $\sigma(b)$ and $h_n(t)\downarrow\mathbb{1}_{F\cap\sigma(b)}(t)$ for every $t\in\sigma(b)$. By the bounded monotone convergence property of the Borel functional calculus for bounded self-adjoint operators,
\begin{align*}
h_n(b)\to \mathbb{1}_{F\cap\sigma(b)}(b)=E_b(F)
\end{align*}
strongly. Since $P$ commutes with every $h_n(b)$, for every $\xi\in H$ we have $Ph_n(b)\xi=h_n(b)P\xi$. Passing to the strong limit gives $PE_b(F)\xi=E_b(F)P\xi$. Since $\xi\in H$ was arbitrary, $F\in\mathcal{D}_b$.
We next show that $\mathcal{D}_b$ is a $\sigma$-algebra. Since $E_b(\mathbb{R})=I$, one has $\mathbb{R}\in\mathcal{D}_b$. If $B\in\mathcal{D}_b$, then $E_b(\mathbb{R}\setminus B)=I-E_b(B)$, so $\mathbb{R}\setminus B\in\mathcal{D}_b$. If $B,C\in\mathcal{D}_b$, then $E_b(B\cap C)=E_b(B)E_b(C)$, and using the commutation of $P$ with both factors gives $P E_b(B\cap C)=E_b(B\cap C)P$. Thus $\mathcal{D}_b$ is closed under finite intersections and finite unions.
Finally, if $(B_n)_{n\in\mathbb{N}}$ is a sequence in $\mathcal{D}_b$ and $C_m=\bigcup_{n=1}^{m}B_n$, then $C_m\in\mathcal{D}_b$ and $C_m$ increases to $B=\bigcup_{n=1}^{\infty}B_n$. Strong countable additivity of the spectral measure gives $E_b(C_m)\to E_b(B)$ strongly. Passing to the strong limit in $PE_b(C_m)=E_b(C_m)P$ gives $PE_b(B)=E_b(B)P$, so $B\in\mathcal{D}_b$. Therefore $\mathcal{D}_b$ is a $\sigma$-algebra containing all closed subsets of $\mathbb{R}$, hence $\mathcal{D}_b=\mathcal{B}(\mathbb{R})$. Thus, for every Borel set $B\subseteq\mathbb{R}$,
\begin{align*}
PE_b(B)=E_b(B)P.
\end{align*}
Substituting $P=E_a(A)$ gives
\begin{align*}
E_a(A)E_b(B)=E_b(B)E_a(A).
\end{align*}
Since $A,B\in\mathcal{B}(\mathbb{R})$ were arbitrary, the spectral projections of $a$ and $b$ commute for all Borel sets.
[/step]