I will focus on the case where $|\rho| < 1$.
We are given that $\mathbb{E}[X_1] = \mathbb{E}[X_2]=0$, $\operatorname{Var}[X_1]=\operatorname{Var}[X_2]=1, $ and the correlation is $\rho$.
From there we can deduce that $\mathbb{E}[X_1^2]= \mathbb{E}[X_2^2]=1$ and $\mathbb{E}[X_1X_2]=\rho$.
Let $g(x_1, x_2; t)=x_1^2+2tx_1x_2+x_2^2$ where $|t| < 1$. We shall study some properties of $g$.
$$g(x_1, x_2;t) = x_1^2+2tx_1x_2+x_2^2 \geq |x_1|^2-2|x_1||x_2|+|x_2|^2 = (|x_1|-|x_2|)^2\ge 0$$
$\mathbb{E}(g(X_1, X_2; t)]=\mathbb{E}[X_1^2+2tX_2X_2+X_2^2]=2+2t\rho$
$g$ is convex and it has a unique stationary point which is the origin.
$\nabla g=\begin{bmatrix} 2x_1+2tx_2 \\ 2tx_1+2x_2\end{bmatrix}$ and $\nabla^2g = \begin{bmatrix} 2 & 2t \\ 2t & 2\end{bmatrix}$, the hessian is positive semidefinite since every principal determinant which are $2$ and $4(1-t^2)$ are nonnegative and by setting the gradient to zero, the stationary point is the origin.
- $\min g(x_1, x_2; t)$ subject to $x_1 \ge \epsilon$ attains minimal value of $\epsilon^2(1-t^2)$.
We have proven earlier that $g$ has a unique stationary point, hence the optimal solution for the constrained optimization problem is when the constraint is active, that is we can set $x_1 = \epsilon$ and it is equivalent to minimizing $ \epsilon^2+2t\epsilon x_2+x_2^2=(x_2+t\epsilon)^2+\epsilon^2(1-t^2)$
- $\min g(x_1, x_2; t)$ subject to $x_1 \le -\epsilon$ attains minimal value of $\epsilon^2(1-t^2)$.
The reasoning is the same as the previous paragraph. Now we set $x_1=-\epsilon$ and we have to minimize $\epsilon^2-2t\epsilon x_2+x_2^2=(x_2-t\epsilon)^2+\epsilon^2(1-t^2)$.
- Let $F= \{x \in \mathbb{R}^2: \|x\|_{\infty} < \epsilon \}$, $\min g(x_1, x_2; t)$ subject to $x \not \in F$ attains value $\epsilon^2(1-t^2)$.
This can be seen from the two optimization problems above. The cases imposing constraint on $x_2$ is similar by symmetry.
Hence we can define $$h(x_1, x_2; t) = \frac{g(x_1, x_2; t)}{\epsilon^2(1-t^2)}$$
- $h$ is nonnegative and $h(x; t) \ge 1$ if $x \not \in F$.
Let $I_A$ denote the indicator variable.
Hence we have $\forall t \in (-1, 1)$, $$\mathbb{E}[I_{F}(X_1, X_2)] + \mathbb{E}[h(X_1, X_2;t)] \ge \mathbb{E}[I_{F}(X_1, X_2)] + \mathbb{E}[I_{F^c}(X_1, X_2)] = 1$$
where the expectation is taken over the joint distribution.
$$P(|X_1| \le \epsilon , |X_2| \le \epsilon) + \frac{2+2\rho t}{\epsilon^2(1-t^2)} \ge 1$$
$$P(|X_1| \le \epsilon , |X_2| \le \epsilon) \ge 1 - \epsilon^{-2} \left[\frac{2+2\rho t}{(1-t^2)} \right]$$
If we can find a $t \in (-1, 1)$ that satisfies
$$\left[\frac{2+2\rho t}{(1-t^2)} \right] = 1+ \sqrt{1-\rho^2}$$
then we have proven the desired result.
$$(1+ \sqrt{1-\rho^2})t^2 + 2 \rho t + 1-\sqrt{1-\rho^2}=0$$
\begin{align}t &=\frac{-2\rho \pm \sqrt{4\rho^2 - 4(1+ \sqrt{1-\rho^2})(1- \sqrt{1-\rho^2})}}{2(1+ \sqrt{1-\rho^2})} \\
&= -\frac{\rho}{1+\sqrt{1-\rho^2}}\end{align} which satisfy the requirement since $|\rho| < 1.$ Hence, we have proven the result.