Well $\text{Pr}(X<a)=0.01$ for $a\approx -2.33$. Since $\text{Var}(X+\epsilon_{i})=1+b$ we need to scale by $\sqrt{1+b}$, so $\text{Pr}(X+\epsilon_{i}<a)=0.01$ for $a\approx -2.33\sqrt{1+b}$.
Then note that
\begin{equation*}
\begin{pmatrix}X_{i}\\X_{j}\end{pmatrix}=
\begin{pmatrix}1 & 1 & 0\\ 1 & 0 & 1\end{pmatrix}
\begin{pmatrix}X\\ \epsilon_{i} \\ \epsilon_{j}\end{pmatrix}
\end{equation*}
so its covariance matrix will be given by
\begin{equation*}
\begin{pmatrix}1 & 1 & 0\\ 1 & 0 & 1\end{pmatrix}
\begin{pmatrix}1 & 0 & 0\\ 0 & b & 0 \\ 0 & 0 & b\end{pmatrix}
\begin{pmatrix}1 & 1\\ 1 & 0 \\ 0 & 1\end{pmatrix}
=\begin{pmatrix}1+b & 1 \\ 1 & 1 + b\end{pmatrix}.
\end{equation*}
After that if you want numbers I think you need to move to Mathematica or something.
Edit :
Not as much as you wanted? Here's a factoid about the bivariate Normal
that's kinda relevant. If
\begin{equation*}
f(\rho,x,y)=\frac{1}{2\pi\sqrt{1-\rho^{2}}}
\exp\left\{-\frac{1}{2(1-\rho^{2})}\left(x^{2}-2\rho x y + y^{2}\right)\right\}
\end{equation*}
and
\begin{equation*}
F(\rho,s,t)=\int_{-\infty}^{s}dx
\int_{-\infty}^{t}dy\, f(\rho,x,y),
\end{equation*}
then
\begin{equation*}
\frac{\partial F}{\partial \rho}(\rho,s,t)=f(\rho,s,t).
\end{equation*}
See e.g. Sungur (1990) Dependence Information in Parameterized
Copulas, Communications in Statistics---Simulation and Computation, 19:4,
1339-1360, DOI: 10.1080/03610919008812920.
To apply this to your problem, redefine $X_{i}$ according to
\begin{equation*}
X_{i}=\frac{1}{\sqrt{1+b}}\left(X+\epsilon_{i}\right)
\end{equation*}
so that
$
\begin{pmatrix}
X_{i} \\ X_{j}
\end{pmatrix}$
has covariance matrix
$
\begin{pmatrix}
1 & \frac{1}{1+b} \\ \frac{1}{1+b} & 1
\end{pmatrix}
$.
Now $a = -2.33$ can remain independent of $b$, and
$ \rho=\frac{1}{1+b}$.
Since total correlation $\rho=1$ implies simultaneous default with probability $0.01$,
\begin{equation*}
F\left(\frac{1}{1+b},a,a\right) + \int_{\frac{1}{1+b}}^{1}\frac{\partial
F}{\partial \rho}(\rho,a,a) \, d\rho= 0.01.
\end{equation*}
and therefore
\begin{align*}
F\left(\frac{1}{1+b},a,a\right)&= 0.01 - \int_{\frac{1}{1+b}}^{1}
f(\rho,a,a) \,
d\rho \\
&= 0.01 - \int_{\frac{1}{1+b}}^{1}
\frac{1}{2\pi\sqrt{1-\rho^{2}}}\exp\left\{{-\frac{a^{2}}{1+\rho}}\right\} \, d\rho.
\end{align*}
This quantity, also known as
$\mathbb{E}\left[Y_{i}Y_{j}\right]$, is I think the one you were asking for. For small $b$,
\begin{align*}
\mathbb{E}\left[Y_{i}Y_{j}\right]&\approx 0.01 -
\frac{1}{2\pi}\exp\left\{-\frac{a^{2}}{2}\right\}\arccos\frac{1}{1+b}\\
&\approx 0.01 - 0.01063 \arccos\frac{1}{1+b}
\end{align*}
