I can't see why the integral $$ J=P(Z_1>0,Z_2>0,Z_1+Z_2<1) $$ can be computed as $$ J = \left(F(\frac{1}{\sqrt{2}})-\frac{1}{2}\right)^2 \qquad (*) $$ where $F$ is the cumulative probability distribution (gaussian) and $Z_1,Z_2$ are independent standard normal random variables.
What I observe : the events are not independent. $Z_1+Z_2 \sim N(0,2)$ but I can't see why $(*)$ is true.
Here is the source : Suess et al., Introduction to Probability Simulation and Gibbs Sampling with R, 2010 (page 69)
One can argue as follow by symmetry, using the transformation $(Z_1,Z_2)\to (Z_1+Z_2,Z_1-Z_2)$ (which preserves independence) : \begin{align*} \mathbb P[-1<Z_1+Z_2 < 1,-1 < Z_1-Z_2 < 1] &= \mathbb P[Z_1> 0,Z_2> 0, Z_1+Z_2 <1]\\&+\mathbb P[Z_1\leq 0,Z_2> 0, Z_1-Z_2 >-1]\\&+\mathbb P[Z_1> 0,Z_2\leq 0, Z_1-Z_2 <1]\\&+\mathbb P[Z_1\leq 0,Z_2\leq 0, Z_1+Z_2 >-1]\\ &=4\cdot \mathbb P[Z_1> 0,Z_2> 0, Z_1+Z_2 <1]\\ \end{align*} But since $Z_1+Z_2$ and $Z_1-Z_2$ are independent we get \begin{align*} \mathbb P[Z_1> 0,Z_2> 0, Z_1+Z_2 <1] &= \frac{1}{4}\cdot \mathbb P[-1<Z_1+Z_2 < 1,-1 < Z_1-Z_2 < 1]\\ &= \frac{1}{4} \cdot \mathbb P[-1<Z_1+Z_2 < 1]\mathbb P[-1 < Z_1-Z_2 < 1]\\ &=\frac{1}{4} \left( F\left(\frac{1}{\sqrt 2}\right)-F\left(-\frac{1}{\sqrt 2}\right) \right)^2\\ &=\frac{1}{4} \left( F\left(\frac{1}{\sqrt 2}\right)-1+F\left(\frac{1}{\sqrt 2}\right) \right)^2\\ &=\left( F\left(\frac{1}{\sqrt 2}\right)-\frac{1}{2} \right)^2 \end{align*}