[proofplan]
We first reduce both rank correlations to probabilities involving signs of centered normal variables. The common analytic input is the bivariate normal orthant formula: if two centered standard normal variables have correlation $\rho$, then the probability that both are positive is $\frac14+\frac{1}{2\pi}\arcsin(\rho)$. Kendall's $\tau$ follows by applying this formula to the difference of two independent copies. Spearman's $\rho_S$ follows by standardising the margins, writing $F_X(X)$ and $F_Y(Y)$ as standard normal distribution functions, and introducing independent threshold variables to convert the expectation into another bivariate normal orthant probability. Throughout the proof, $\mathcal{L}^1$ denotes Lebesgue measure on $\mathbb{R}$ and $\mathcal{L}^2$ denotes Lebesgue measure on $\mathbb{R}^2$.
[/proofplan]
[step:Establish the bivariate normal orthant formula]
We prove the auxiliary formula used twice below. Let $(U,V)$ be a centered bivariate normal random vector with $\mathbb{E}[U]=\mathbb{E}[V]=0$, $\operatorname{Var}(U)=\operatorname{Var}(V)=1$, and correlation $\rho \in [-1,1]$. Then
\begin{align*}
\mathbb{P}(U>0,\ V>0)=\frac14+\frac{1}{2\pi}\arcsin(\rho).
\end{align*}
For $|\rho|<1$, the joint density $f_\rho:\mathbb{R}^2\to[0,\infty)$ of $(U,V)$ with respect to $\mathcal{L}^2$ is
\begin{align*}
f_\rho(u,v)
=
\frac{1}{2\pi\sqrt{1-\rho^2}}
\exp\left(
-\frac{u^2-2\rho uv+v^2}{2(1-\rho^2)}
\right).
\end{align*}
Define $Q:(-1,1)\to[0,1]$ by
\begin{align*}
Q(\rho):=\int_0^\infty\int_0^\infty f_\rho(u,v)\,d\mathcal{L}^1(v)\,d\mathcal{L}^1(u).
\end{align*}
A direct differentiation of $f_\rho$ gives
\begin{align*}
\frac{\partial f_\rho}{\partial \rho}(u,v)
=
\frac{\partial^2 f_\rho}{\partial u\,\partial v}(u,v).
\end{align*}
Since $f_\rho$ and its first derivatives decay exponentially on every compact subinterval of $\rho\in(-1,1)$, differentiating under the integral sign and integrating by parts on $(0,\infty)^2$ are justified. Hence
\begin{align*}
Q'(\rho)
&=
\int_0^\infty\int_0^\infty
\frac{\partial^2 f_\rho}{\partial u\,\partial v}(u,v)
\,d\mathcal{L}^1(v)\,d\mathcal{L}^1(u)\\
&=
f_\rho(0,0)\\
&=
\frac{1}{2\pi\sqrt{1-\rho^2}}.
\end{align*}
At $\rho=0$, the variables $U$ and $V$ are independent standard normal variables, so
\begin{align*}
Q(0)=\mathbb{P}(U>0)\mathbb{P}(V>0)=\frac12\cdot\frac12=\frac14.
\end{align*}
Therefore, for $|\rho|<1$,
\begin{align*}
Q(\rho)
&=
Q(0)+\int_0^\rho Q'(s)\,d\mathcal{L}^1(s)\\
&=
\frac14+\int_0^\rho \frac{1}{2\pi\sqrt{1-s^2}}\,d\mathcal{L}^1(s)\\
&=
\frac14+\frac{1}{2\pi}\arcsin(\rho).
\end{align*}
For $\rho=1$, we have $V=U$ almost surely, hence $\mathbb{P}(U>0,V>0)=\mathbb{P}(U>0)=\frac12$. For $\rho=-1$, we have $V=-U$ almost surely, hence $\mathbb{P}(U>0,V>0)=0$. These agree with the same formula at $\rho=\pm1$.
[guided]
The orthant formula is the only analytic computation in the proof. We want to compute the probability that a centered bivariate normal vector falls in the first quadrant.
Assume first that $|\rho|<1$, so the joint distribution has a density. The density is the map
\begin{align*}
f_\rho:\mathbb{R}^2&\to[0,\infty)\\
(u,v)&\mapsto
\frac{1}{2\pi\sqrt{1-\rho^2}}
\exp\left(
-\frac{u^2-2\rho uv+v^2}{2(1-\rho^2)}
\right).
\end{align*}
Define the quadrant probability as the function
\begin{align*}
Q:(-1,1)&\to[0,1]\\
\rho&\mapsto
\int_0^\infty\int_0^\infty f_\rho(u,v)\,d\mathcal{L}^1(v)\,d\mathcal{L}^1(u).
\end{align*}
The key identity is
\begin{align*}
\frac{\partial f_\rho}{\partial \rho}(u,v)
=
\frac{\partial^2 f_\rho}{\partial u\,\partial v}(u,v).
\end{align*}
This identity is obtained by differentiating the displayed density with respect to $\rho$, $u$, and $v$ and comparing the resulting expressions. Because the Gaussian exponential gives uniform exponential decay on compact subintervals of $\rho\in(-1,1)$, we may differentiate under the integral sign. Thus
\begin{align*}
Q'(\rho)
=
\int_0^\infty\int_0^\infty
\frac{\partial f_\rho}{\partial \rho}(u,v)
\,d\mathcal{L}^1(v)\,d\mathcal{L}^1(u).
\end{align*}
Substituting the derivative identity gives
\begin{align*}
Q'(\rho)
=
\int_0^\infty\int_0^\infty
\frac{\partial^2 f_\rho}{\partial u\,\partial v}(u,v)
\,d\mathcal{L}^1(v)\,d\mathcal{L}^1(u).
\end{align*}
Integrating first in $v$ and then in $u$, the boundary terms at infinity vanish because the density and its derivatives decay exponentially. The only remaining boundary value is the value at $(0,0)$:
\begin{align*}
Q'(\rho)
=
f_\rho(0,0)
=
\frac{1}{2\pi\sqrt{1-\rho^2}}.
\end{align*}
At $\rho=0$, the covariance matrix is the identity matrix, so $U$ and $V$ are independent standard normal variables. Hence
\begin{align*}
Q(0)=\mathbb{P}(U>0)\mathbb{P}(V>0)=\frac12\cdot\frac12=\frac14.
\end{align*}
Integrating the derivative from $0$ to $\rho$ gives
\begin{align*}
Q(\rho)
&=
Q(0)+\int_0^\rho Q'(s)\,d\mathcal{L}^1(s)\\
&=
\frac14+\int_0^\rho \frac{1}{2\pi\sqrt{1-s^2}}\,d\mathcal{L}^1(s)\\
&=
\frac14+\frac{1}{2\pi}\arcsin(\rho).
\end{align*}
The endpoint cases are singular but immediate. If $\rho=1$, then $V=U$ almost surely, so the first-quadrant event is just $\{U>0\}$ and has probability $\frac12$. If $\rho=-1$, then $V=-U$ almost surely, so the event $\{U>0,V>0\}$ has probability $0$. These are exactly the values of $\frac14+\frac{1}{2\pi}\arcsin(\rho)$ at $\rho=1$ and $\rho=-1$.
[/guided]
[/step]
[step:Reduce Kendall's $\tau$ to a centered normal orthant probability]
Let $(X_1,Y_1)$ and $(X_2,Y_2)$ be independent copies of $(X,Y)$. Define the difference variables
\begin{align*}
D_X&:=X_1-X_2,\\
D_Y&:=Y_1-Y_2.
\end{align*}
Since independent linear combinations of jointly normal random variables are jointly normal, $(D_X,D_Y)$ is a centered bivariate normal random vector. Moreover,
\begin{align*}
\operatorname{Var}(D_X)&=\operatorname{Var}(X_1)+\operatorname{Var}(X_2)=2\operatorname{Var}(X),\\
\operatorname{Var}(D_Y)&=\operatorname{Var}(Y_1)+\operatorname{Var}(Y_2)=2\operatorname{Var}(Y),\\
\operatorname{Cov}(D_X,D_Y)
&=
\operatorname{Cov}(X_1,Y_1)+\operatorname{Cov}(X_2,Y_2)\\
&=
2\operatorname{Cov}(X,Y).
\end{align*}
Therefore the correlation of $D_X$ and $D_Y$ is
\begin{align*}
\frac{\operatorname{Cov}(D_X,D_Y)}
{\sqrt{\operatorname{Var}(D_X)\operatorname{Var}(D_Y)}}
=
\frac{2\operatorname{Cov}(X,Y)}
{2\sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}}
=
r.
\end{align*}
Because $(D_X,D_Y)$ has continuous marginals unless $r=\pm1$, ties have probability zero for $|r|<1$. In the endpoint cases $r=\pm1$, the conclusion follows directly from $D_Y=\pm cD_X$ almost surely for some $c>0$. Thus it remains to consider $|r|<1$.
Let
\begin{align*}
U:=\frac{D_X}{\sqrt{\operatorname{Var}(D_X)}},
\qquad
V:=\frac{D_Y}{\sqrt{\operatorname{Var}(D_Y)}}.
\end{align*}
Then $(U,V)$ is centered standard bivariate normal with correlation $r$. Since multiplication by positive constants preserves signs,
\begin{align*}
\mathbb{P}(D_X>0,D_Y>0)=\mathbb{P}(U>0,V>0).
\end{align*}
By symmetry of the centered bivariate normal distribution under $(U,V)\mapsto(-U,-V)$,
\begin{align*}
\mathbb{P}(D_X<0,D_Y<0)=\mathbb{P}(D_X>0,D_Y>0).
\end{align*}
Thus
\begin{align*}
\mathbb{P}(D_XD_Y>0)=2\mathbb{P}(U>0,V>0).
\end{align*}
Since $\mathbb{P}(D_XD_Y=0)=0$ for $|r|<1$,
\begin{align*}
\tau
&=
\mathbb{E}[\operatorname{sgn}(D_XD_Y)]\\
&=
\mathbb{P}(D_XD_Y>0)-\mathbb{P}(D_XD_Y<0)\\
&=
2\mathbb{P}(D_XD_Y>0)-1\\
&=
4\mathbb{P}(U>0,V>0)-1.
\end{align*}
Applying the orthant formula with $\rho=r$ gives
\begin{align*}
\tau
=
4\left(\frac14+\frac{1}{2\pi}\arcsin(r)\right)-1
=
\frac{2}{\pi}\arcsin(r).
\end{align*}
[guided]
Kendall's $\tau$ compares the signs of the two coordinate differences in two independent observations. We therefore introduce the random variables
\begin{align*}
D_X&:=X_1-X_2,\\
D_Y&:=Y_1-Y_2.
\end{align*}
The event $D_XD_Y>0$ means the two observations are concordant, and the event $D_XD_Y<0$ means they are discordant.
Because $(X_1,Y_1)$ and $(X_2,Y_2)$ are independent copies of a bivariate normal vector, every linear combination of their coordinates is normal. Hence $(D_X,D_Y)$ is bivariate normal. Its mean is zero, and its covariance data are
\begin{align*}
\operatorname{Var}(D_X)&=\operatorname{Var}(X_1)+\operatorname{Var}(X_2)=2\operatorname{Var}(X),\\
\operatorname{Var}(D_Y)&=\operatorname{Var}(Y_1)+\operatorname{Var}(Y_2)=2\operatorname{Var}(Y),\\
\operatorname{Cov}(D_X,D_Y)
&=
\operatorname{Cov}(X_1-X_2,Y_1-Y_2)\\
&=
\operatorname{Cov}(X_1,Y_1)+\operatorname{Cov}(X_2,Y_2)\\
&=
2\operatorname{Cov}(X,Y).
\end{align*}
The mixed covariance terms vanish because the two copies are independent. Therefore the correlation of $D_X$ and $D_Y$ is
\begin{align*}
\frac{\operatorname{Cov}(D_X,D_Y)}
{\sqrt{\operatorname{Var}(D_X)\operatorname{Var}(D_Y)}}
=
\frac{2\operatorname{Cov}(X,Y)}
{2\sqrt{\operatorname{Var}(X)\operatorname{Var}(Y)}}
=
r.
\end{align*}
Now standardise the two difference variables by defining
\begin{align*}
U:=\frac{D_X}{\sqrt{\operatorname{Var}(D_X)}},
\qquad
V:=\frac{D_Y}{\sqrt{\operatorname{Var}(D_Y)}}.
\end{align*}
The constants in the denominators are positive, so $U$ and $D_X$ have the same sign, and $V$ and $D_Y$ have the same sign. Thus
\begin{align*}
\mathbb{P}(D_X>0,D_Y>0)=\mathbb{P}(U>0,V>0).
\end{align*}
For $|r|<1$, the bivariate normal law is nondegenerate, so $\mathbb{P}(D_XD_Y=0)=0$. The centered normal law is invariant under the map $(U,V)\mapsto(-U,-V)$, so the two concordant quadrants have the same probability:
\begin{align*}
\mathbb{P}(D_X<0,D_Y<0)=\mathbb{P}(D_X>0,D_Y>0).
\end{align*}
Therefore
\begin{align*}
\mathbb{P}(D_XD_Y>0)=2\mathbb{P}(U>0,V>0).
\end{align*}
Using the definition of Kendall's $\tau$,
\begin{align*}
\tau
&=
\mathbb{E}[\operatorname{sgn}(D_XD_Y)]\\
&=
\mathbb{P}(D_XD_Y>0)-\mathbb{P}(D_XD_Y<0)\\
&=
2\mathbb{P}(D_XD_Y>0)-1\\
&=
4\mathbb{P}(U>0,V>0)-1.
\end{align*}
The orthant formula applies to $(U,V)$ because it is centered standard bivariate normal with correlation $r$. Hence
\begin{align*}
\tau
=
4\left(\frac14+\frac{1}{2\pi}\arcsin(r)\right)-1
=
\frac{2}{\pi}\arcsin(r).
\end{align*}
If $r=1$, then $Y$ is an increasing affine function of $X$ almost surely, so $D_Y=cD_X$ almost surely for some $c>0$ and $\tau=1$. If $r=-1$, then $D_Y=-cD_X$ almost surely for some $c>0$ and $\tau=-1$. These values agree with $\frac{2}{\pi}\arcsin(r)$ at $r=\pm1$.
[/guided]
[/step]
[step:Represent Spearman's $\rho_S$ using standard normal distribution functions]
Let $\mu_X:=\mathbb{E}[X]$, $\mu_Y:=\mathbb{E}[Y]$, $\sigma_X:=\sqrt{\operatorname{Var}(X)}$, and $\sigma_Y:=\sqrt{\operatorname{Var}(Y)}$. Define the standardised variables
\begin{align*}
\widetilde{X}:=\frac{X-\mu_X}{\sigma_X},
\qquad
\widetilde{Y}:=\frac{Y-\mu_Y}{\sigma_Y}.
\end{align*}
Then $(\widetilde{X},\widetilde{Y})$ is centered standard bivariate normal with correlation $r$.
Let $\Phi:\mathbb{R}\to[0,1]$ denote the standard normal distribution function. Since the marginal laws of $X$ and $Y$ are normal,
\begin{align*}
F_X(X)=\Phi(\widetilde{X}),
\qquad
F_Y(Y)=\Phi(\widetilde{Y})
\end{align*}
almost surely. The random variables $\Phi(\widetilde{X})$ and $\Phi(\widetilde{Y})$ are uniform on $[0,1]$, hence each has expectation $\frac12$ and variance $\frac{1}{12}$. Therefore
\begin{align*}
\rho_S
&=
\operatorname{Corr}(\Phi(\widetilde{X}),\Phi(\widetilde{Y}))\\
&=
\frac{\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]-\frac14}{\frac{1}{12}}\\
&=
12\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]-3.
\end{align*}
[/step]
[step:Convert the Spearman expectation into an orthant probability]
Let $Z_1$ and $Z_2$ be independent standard normal random variables, independent of $(\widetilde{X},\widetilde{Y})$. Define
\begin{align*}
A:=\widetilde{X}-Z_1,
\qquad
B:=\widetilde{Y}-Z_2.
\end{align*}
Conditioning on $(\widetilde{X},\widetilde{Y})$ gives
\begin{align*}
\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]
&=
\mathbb{P}(Z_1\leq \widetilde{X},\ Z_2\leq \widetilde{Y})\\
&=
\mathbb{P}(A\geq 0,\ B\geq 0).
\end{align*}
The vector $(A,B)$ is centered bivariate normal. Its variances and covariance are
\begin{align*}
\operatorname{Var}(A)&=\operatorname{Var}(\widetilde{X})+\operatorname{Var}(Z_1)=2,\\
\operatorname{Var}(B)&=\operatorname{Var}(\widetilde{Y})+\operatorname{Var}(Z_2)=2,\\
\operatorname{Cov}(A,B)&=\operatorname{Cov}(\widetilde{X},\widetilde{Y})=r.
\end{align*}
Hence the correlation of $A$ and $B$ is
\begin{align*}
\frac{\operatorname{Cov}(A,B)}
{\sqrt{\operatorname{Var}(A)\operatorname{Var}(B)}}
=
\frac{r}{2}.
\end{align*}
Since $(A,B)$ has continuous marginals, the events $\{A\geq0,B\geq0\}$ and $\{A>0,B>0\}$ have the same probability. Applying the orthant formula with $\rho=r/2$ gives
\begin{align*}
\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]
=
\frac14+\frac{1}{2\pi}\arcsin\left(\frac{r}{2}\right).
\end{align*}
[guided]
The expectation $\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]$ is converted into a probability by adding independent threshold variables. Let $Z_1$ and $Z_2$ be independent standard normal random variables, independent of $(\widetilde{X},\widetilde{Y})$, and define
\begin{align*}
A:=\widetilde{X}-Z_1,
\qquad
B:=\widetilde{Y}-Z_2.
\end{align*}
Condition on the values of $\widetilde{X}$ and $\widetilde{Y}$. Since $Z_1$ and $Z_2$ are independent standard normal variables,
\begin{align*}
\mathbb{P}(Z_1\leq \widetilde{X},\ Z_2\leq \widetilde{Y}\mid \widetilde{X},\widetilde{Y})
=
\Phi(\widetilde{X})\Phi(\widetilde{Y}).
\end{align*}
Taking expectations gives
\begin{align*}
\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]
=
\mathbb{P}(Z_1\leq \widetilde{X},\ Z_2\leq \widetilde{Y})
=
\mathbb{P}(A\geq0,\ B\geq0).
\end{align*}
Now compute the distribution of $(A,B)$. Since $A$ and $B$ are linear combinations of jointly normal variables, $(A,B)$ is bivariate normal. It is centered because each of $\widetilde{X},\widetilde{Y},Z_1,Z_2$ has mean zero. Its variances are
\begin{align*}
\operatorname{Var}(A)&=\operatorname{Var}(\widetilde{X})+\operatorname{Var}(Z_1)=1+1=2,\\
\operatorname{Var}(B)&=\operatorname{Var}(\widetilde{Y})+\operatorname{Var}(Z_2)=1+1=2.
\end{align*}
The cross terms vanish because $Z_1$ and $Z_2$ are independent of $(\widetilde{X},\widetilde{Y})$ and of each other. The covariance is
\begin{align*}
\operatorname{Cov}(A,B)
&=
\operatorname{Cov}(\widetilde{X}-Z_1,\widetilde{Y}-Z_2)\\
&=
\operatorname{Cov}(\widetilde{X},\widetilde{Y})\\
&=
r.
\end{align*}
Therefore the correlation coefficient of $A$ and $B$ is
\begin{align*}
\frac{\operatorname{Cov}(A,B)}
{\sqrt{\operatorname{Var}(A)\operatorname{Var}(B)}}
=
\frac{r}{2}.
\end{align*}
Because $(A,B)$ has continuous marginals, replacing $\geq0$ by $>0$ does not change the probability. The orthant formula applies to the standardised version of $(A,B)$, whose correlation is still $r/2$. Hence
\begin{align*}
\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]
=
\mathbb{P}(A\geq0,\ B\geq0)
=
\frac14+\frac{1}{2\pi}\arcsin\left(\frac{r}{2}\right).
\end{align*}
[/guided]
[/step]
[step:Substitute the orthant probability into Spearman's formula]
From the previous two steps,
\begin{align*}
\rho_S
&=
12\mathbb{E}[\Phi(\widetilde{X})\Phi(\widetilde{Y})]-3\\
&=
12\left(
\frac14+\frac{1}{2\pi}\arcsin\left(\frac{r}{2}\right)
\right)-3\\
&=
\frac{6}{\pi}\arcsin\left(\frac{r}{2}\right).
\end{align*}
Together with the Kendall computation, this proves both displayed identities.
[/step]