[proofplan]
We realize the Gaussian copula through a bivariate normal random vector and rewrite the upper tail coefficient as a quotient of joint and marginal Gaussian tail probabilities. A linear orthogonal-normal transformation diagonalizes the covariance structure and reduces the joint upper tail event to a one-dimensional standard normal tail with a larger threshold. Comparing this exponential upper bound with a matching elementary lower bound for the marginal normal tail proves that the quotient tends to zero. The lower tail follows by applying the same upper-tail estimate to $(-Z_1,-Z_2)$, which has the same correlation.
[/proofplan]
[step:Represent the copula tails as Gaussian tail probabilities]
Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space supporting a centered Gaussian random vector
\begin{align*}
Z: \Omega &\to \mathbb{R}^2\\
\omega &\mapsto (Z_1(\omega),Z_2(\omega))
\end{align*}
with covariance matrix
\begin{align*}
\Sigma_r =
\begin{pmatrix}
1 & r\\
r & 1
\end{pmatrix}.
\end{align*}
For $u \in (0,1)$, define $x_u := \Phi^{-1}(u) \in \mathbb{R}$. Since each coordinate has standard normal distribution,
\begin{align*}
1-u = \mathbb{P}(Z_1 > x_u).
\end{align*}
Also,
\begin{align*}
1 - 2u + C_r(u,u)
&= \mathbb{P}(Z_1 > x_u, Z_2 > x_u).
\end{align*}
Indeed, $C_r(u,u)=\mathbb{P}(Z_1 \le x_u,Z_2 \le x_u)$, and inclusion-exclusion gives
\begin{align*}
\mathbb{P}(Z_1 > x_u, Z_2 > x_u)
&= 1-\mathbb{P}(Z_1 \le x_u)-\mathbb{P}(Z_2 \le x_u)+\mathbb{P}(Z_1 \le x_u,Z_2 \le x_u)\\
&= 1-2u+C_r(u,u).
\end{align*}
Thus
\begin{align*}
\frac{1 - 2u + C_r(u,u)}{1-u}
=
\frac{\mathbb{P}(Z_1 > x_u,Z_2 > x_u)}{\mathbb{P}(Z_1 > x_u)}.
\end{align*}
As $u \uparrow 1$, $x_u \to \infty$, so it remains to prove
\begin{align*}
\lim_{x \to \infty}
\frac{\mathbb{P}(Z_1 > x,Z_2 > x)}{\mathbb{P}(Z_1 > x)}
=0.
\end{align*}
[guided]
The copula is defined by applying the bivariate Gaussian distribution to normal quantiles. Therefore the diagonal value $C_r(u,u)$ is exactly the probability that both Gaussian coordinates lie below the same threshold $x_u=\Phi^{-1}(u)$:
\begin{align*}
C_r(u,u)=\mathbb{P}(Z_1 \le x_u,Z_2 \le x_u).
\end{align*}
The upper tail coefficient uses the expression $1-2u+C_r(u,u)$. This expression is not accidental: by inclusion-exclusion it is the probability of the simultaneous upper tail event. Since $\mathbb{P}(Z_1 \le x_u)=\mathbb{P}(Z_2 \le x_u)=u$, we have
\begin{align*}
\mathbb{P}(Z_1 > x_u,Z_2 > x_u)
&= 1-\mathbb{P}(Z_1 \le x_u)-\mathbb{P}(Z_2 \le x_u)+\mathbb{P}(Z_1 \le x_u,Z_2 \le x_u)\\
&= 1-2u+C_r(u,u).
\end{align*}
Similarly, because $Z_1$ is standard normal,
\begin{align*}
1-u=\mathbb{P}(Z_1>x_u).
\end{align*}
Thus the upper tail coefficient is the limit of a conditional-type quotient:
\begin{align*}
\frac{1 - 2u + C_r(u,u)}{1-u}
=
\frac{\mathbb{P}(Z_1 > x_u,Z_2 > x_u)}{\mathbb{P}(Z_1 > x_u)}.
\end{align*}
The limit $u \uparrow 1$ corresponds exactly to $x_u=\Phi^{-1}(u)\to\infty$. Hence the problem becomes a comparison between a bivariate Gaussian joint tail and a one-dimensional Gaussian tail.
[/guided]
[/step]
[step:Diagonalize the correlated Gaussian vector]
Let $U: \Omega \to \mathbb{R}$ and $V: \Omega \to \mathbb{R}$ be the real-valued random variables
\begin{align*}
U &= \frac{Z_1+Z_2}{\sqrt{2(1+r)}},&
V &= \frac{Z_1-Z_2}{\sqrt{2(1-r)}}.
\end{align*}
Since $(Z_1,Z_2)$ is centered Gaussian, every linear image is Gaussian. Direct covariance computation gives
\begin{align*}
\mathbb{E}[U] &= 0,&
\mathbb{E}[V] &= 0,\\
\operatorname{Var}(U)
&= \frac{\operatorname{Var}(Z_1)+2\operatorname{Cov}(Z_1,Z_2)+\operatorname{Var}(Z_2)}{2(1+r)}
=1,\\
\operatorname{Var}(V)
&= \frac{\operatorname{Var}(Z_1)-2\operatorname{Cov}(Z_1,Z_2)+\operatorname{Var}(Z_2)}{2(1-r)}
=1,\\
\operatorname{Cov}(U,V)
&= \frac{\operatorname{Var}(Z_1)-\operatorname{Var}(Z_2)}{2\sqrt{(1+r)(1-r)}}=0.
\end{align*}
Because a centered Gaussian vector with diagonal covariance has independent coordinates, $U$ and $V$ are independent standard normal random variables.
Solving for $Z_1$ and $Z_2$ gives
\begin{align*}
Z_1 &= \frac{\sqrt{1+r}\,U+\sqrt{1-r}\,V}{\sqrt{2}},\\
Z_2 &= \frac{\sqrt{1+r}\,U-\sqrt{1-r}\,V}{\sqrt{2}}.
\end{align*}
Therefore, for $x>0$,
\begin{align*}
\{Z_1>x,\ Z_2>x\}
\subset
\left\{U>\sqrt{\frac{2}{1+r}}\,x\right\}.
\end{align*}
Indeed, adding the two inequalities $Z_1>x$ and $Z_2>x$ gives $\sqrt{2(1+r)}\,U>2x$.
[guided]
The dependence between $Z_1$ and $Z_2$ is carried by their sum and difference. Define
\begin{align*}
U &= \frac{Z_1+Z_2}{\sqrt{2(1+r)}},&
V &= \frac{Z_1-Z_2}{\sqrt{2(1-r)}}.
\end{align*}
The denominators are chosen so that both variables have variance $1$. We verify this explicitly:
\begin{align*}
\operatorname{Var}(U)
&= \frac{\operatorname{Var}(Z_1)+2\operatorname{Cov}(Z_1,Z_2)+\operatorname{Var}(Z_2)}{2(1+r)}
= \frac{1+2r+1}{2(1+r)}
=1,\\
\operatorname{Var}(V)
&= \frac{\operatorname{Var}(Z_1)-2\operatorname{Cov}(Z_1,Z_2)+\operatorname{Var}(Z_2)}{2(1-r)}
= \frac{1-2r+1}{2(1-r)}
=1.
\end{align*}
Their covariance is
\begin{align*}
\operatorname{Cov}(U,V)
&= \frac{\operatorname{Cov}(Z_1+Z_2,Z_1-Z_2)}{2\sqrt{(1+r)(1-r)}}\\
&= \frac{\operatorname{Var}(Z_1)-\operatorname{Cov}(Z_1,Z_2)+\operatorname{Cov}(Z_2,Z_1)-\operatorname{Var}(Z_2)}{2\sqrt{(1+r)(1-r)}}\\
&=0.
\end{align*}
Since $(U,V)$ is a centered Gaussian vector, zero covariance implies independence. Thus $U$ and $V$ are independent standard normal variables.
Now invert the linear transformation:
\begin{align*}
Z_1 &= \frac{\sqrt{1+r}\,U+\sqrt{1-r}\,V}{\sqrt{2}},\\
Z_2 &= \frac{\sqrt{1+r}\,U-\sqrt{1-r}\,V}{\sqrt{2}}.
\end{align*}
If both $Z_1>x$ and $Z_2>x$, then adding the two inequalities eliminates $V$:
\begin{align*}
Z_1+Z_2 > 2x.
\end{align*}
But $Z_1+Z_2=\sqrt{2(1+r)}\,U$, so
\begin{align*}
U>\sqrt{\frac{2}{1+r}}\,x.
\end{align*}
This gives the event inclusion
\begin{align*}
\{Z_1>x,\ Z_2>x\}
\subset
\left\{U>\sqrt{\frac{2}{1+r}}\,x\right\}.
\end{align*}
The important point is that $\sqrt{2/(1+r)}>1$ whenever $r<1$, so the joint event forces a strictly more extreme one-dimensional Gaussian tail.
[/guided]
[/step]
[step:Compare the joint upper tail with the marginal upper tail]
Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. Define the standard normal density $\varphi: \mathbb{R} \to (0,\infty)$ by
\begin{align*}
\varphi(t)=\frac{1}{\sqrt{2\pi}}e^{-t^2/2}.
\end{align*}
For $a:=\sqrt{2/(1+r)}>1$ and $x>0$, the previous step gives
\begin{align*}
\mathbb{P}(Z_1>x,Z_2>x)
\le \mathbb{P}(U>ax).
\end{align*}
Using the elementary upper Gaussian tail bound
\begin{align*}
\mathbb{P}(U>y)
=
\int_y^\infty \varphi(t)\,d\mathcal{L}^1(t)
\le
\frac{1}{y\sqrt{2\pi}}e^{-y^2/2},
\qquad y>0,
\end{align*}
with $y=ax$, we obtain
\begin{align*}
\mathbb{P}(Z_1>x,Z_2>x)
\le
\frac{1}{ax\sqrt{2\pi}}e^{-a^2x^2/2}.
\end{align*}
For the marginal lower bound, since $\varphi$ is decreasing on $[0,\infty)$,
\begin{align*}
\mathbb{P}(Z_1>x)
=
\int_x^\infty \varphi(t)\,d\mathcal{L}^1(t)
\ge
\int_x^{x+1/x}\varphi(t)\,d\mathcal{L}^1(t)
\ge
\frac{1}{x}\varphi\left(x+\frac{1}{x}\right)
\end{align*}
for $x>0$. Since
\begin{align*}
\left(x+\frac{1}{x}\right)^2
=
x^2+2+\frac{1}{x^2}
\le x^2+3
\end{align*}
for $x \ge 1$, we have
\begin{align*}
\mathbb{P}(Z_1>x)
\ge
\frac{e^{-3/2}}{x\sqrt{2\pi}}e^{-x^2/2},
\qquad x\ge 1.
\end{align*}
Combining the two estimates yields, for $x\ge1$,
\begin{align*}
0
\le
\frac{\mathbb{P}(Z_1>x,Z_2>x)}{\mathbb{P}(Z_1>x)}
\le
\frac{e^{3/2}}{a}
\exp\left(-\frac{a^2-1}{2}x^2\right).
\end{align*}
Because
\begin{align*}
a^2-1
=
\frac{2}{1+r}-1
=
\frac{1-r}{1+r}
>0,
\end{align*}
the right-hand side tends to $0$ as $x\to\infty$. Hence
\begin{align*}
\lim_{x\to\infty}
\frac{\mathbb{P}(Z_1>x,Z_2>x)}{\mathbb{P}(Z_1>x)}
=0.
\end{align*}
[guided]
Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. The event inclusion from the previous step reduces the joint tail to a one-dimensional Gaussian tail:
\begin{align*}
\mathbb{P}(Z_1>x,Z_2>x)
\le
\mathbb{P}\left(U>\sqrt{\frac{2}{1+r}}\,x\right).
\end{align*}
Set
\begin{align*}
a:=\sqrt{\frac{2}{1+r}}.
\end{align*}
Since $r<1$, we have $a>1$. The standard Gaussian tail estimate gives, for $y>0$,
\begin{align*}
\mathbb{P}(U>y)
=
\int_y^\infty \frac{1}{\sqrt{2\pi}}e^{-t^2/2}\,d\mathcal{L}^1(t)
\le
\frac{1}{y\sqrt{2\pi}}e^{-y^2/2}.
\end{align*}
This estimate follows by writing, for $t\ge y>0$, $1\le t/y$, and therefore
\begin{align*}
\int_y^\infty e^{-t^2/2}\,d\mathcal{L}^1(t)
\le
\frac{1}{y}\int_y^\infty t e^{-t^2/2}\,d\mathcal{L}^1(t)
=
\frac{1}{y}e^{-y^2/2}.
\end{align*}
Applying this with $y=ax$ gives
\begin{align*}
\mathbb{P}(Z_1>x,Z_2>x)
\le
\frac{1}{ax\sqrt{2\pi}}e^{-a^2x^2/2}.
\end{align*}
We also need a lower bound for the denominator $\mathbb{P}(Z_1>x)$, because the quotient could only be controlled if the marginal tail is not underestimated. Since the standard normal density
\begin{align*}
\varphi: \mathbb{R} &\to (0,\infty)\\
t &\mapsto \frac{1}{\sqrt{2\pi}}e^{-t^2/2}
\end{align*}
is decreasing on $[0,\infty)$, for $x>0$ we restrict the integration domain from $(x,\infty)$ to the smaller interval $(x,x+1/x)$ and obtain
\begin{align*}
\mathbb{P}(Z_1>x)
=
\int_x^\infty \varphi(t)\,d\mathcal{L}^1(t)
\ge
\int_x^{x+1/x}\varphi(t)\,d\mathcal{L}^1(t)
\ge
\frac{1}{x}\varphi\left(x+\frac{1}{x}\right).
\end{align*}
For $x\ge1$,
\begin{align*}
\left(x+\frac{1}{x}\right)^2
=
x^2+2+\frac{1}{x^2}
\le x^2+3,
\end{align*}
so
\begin{align*}
\mathbb{P}(Z_1>x)
\ge
\frac{e^{-3/2}}{x\sqrt{2\pi}}e^{-x^2/2}.
\end{align*}
Dividing the joint upper bound by this marginal lower bound gives
\begin{align*}
0
\le
\frac{\mathbb{P}(Z_1>x,Z_2>x)}{\mathbb{P}(Z_1>x)}
\le
\frac{e^{3/2}}{a}
\exp\left(-\frac{a^2-1}{2}x^2\right).
\end{align*}
The exponent is genuinely negative because
\begin{align*}
a^2-1
=
\frac{2}{1+r}-1
=
\frac{1-r}{1+r}
>0.
\end{align*}
This strict positivity is exactly where the assumption $r<1$ is used. Therefore the right-hand side tends to $0$, and the desired upper-tail quotient tends to $0$.
[/guided]
[/step]
[step:Conclude that the upper tail dependence coefficient is zero]
Since $x_u=\Phi^{-1}(u)\to\infty$ as $u\uparrow1$, the limit from the previous step gives
\begin{align*}
\lambda_U
&=
\lim_{u\uparrow1}
\frac{1-2u+C_r(u,u)}{1-u}\\
&=
\lim_{u\uparrow1}
\frac{\mathbb{P}(Z_1>x_u,Z_2>x_u)}{\mathbb{P}(Z_1>x_u)}
=0.
\end{align*}
Thus the Gaussian copula has no upper tail dependence for every $r\in(-1,1)$.
[/step]
[step:Apply the same upper-tail estimate to the reflected vector for the lower tail]
For $u\in(0,1)$, again write $x_u=\Phi^{-1}(u)$. Then
\begin{align*}
C_r(u,u)
=
\mathbb{P}(Z_1\le x_u,Z_2\le x_u).
\end{align*}
Define the reflected Gaussian vector $W:\Omega\to\mathbb{R}^2$ by
\begin{align*}
W(\omega)=(-Z_1(\omega),-Z_2(\omega)).
\end{align*}
The vector $W=(W_1,W_2)$ is centered Gaussian with unit variances and correlation $r$. If $u\downarrow0$, then $x_u\to-\infty$ and $y_u:=-x_u\to\infty$. Therefore
\begin{align*}
\frac{C_r(u,u)}{u}
&=
\frac{\mathbb{P}(Z_1\le x_u,Z_2\le x_u)}{\mathbb{P}(Z_1\le x_u)}\\
&=
\frac{\mathbb{P}(W_1\ge y_u,W_2\ge y_u)}{\mathbb{P}(W_1\ge y_u)}.
\end{align*}
Replacing non-strict inequalities by strict inequalities does not change these probabilities because one-dimensional and two-dimensional Gaussian distributions assign probability zero to affine hyperplanes. Applying the already proved upper-tail estimate to $W$ yields
\begin{align*}
\lim_{u\downarrow0}\frac{C_r(u,u)}{u}=0.
\end{align*}
Thus $\lambda_L=0$. Combining this with $\lambda_U=0$ proves
\begin{align*}
\lambda_L=\lambda_U=0.
\end{align*}
[/step]
