Conditional multivariate normal pdf with inequality $f(x_1 | x_2 > a)$

You need a normalization factor for the integral. Basically the right hand side is the weighted sum of some conditional densities (Gaussians with possibly different mean and variance) then the normalization should be the sum of weighing factors. Namely it should read:

$$ f(x_1 | x_2 > \alpha) =\frac{\int_\alpha^\infty g(x_1 | x_2 = u) \; f_2 (x_2=u)\mbox{d}u}{\int_\alpha^\infty f_2 (x_2=u)\mbox{d}u} \, $$

I took an easier approach and I think I got it. Just to be a bit more general I'll be using $(\mu_1, \mu_2 ) $ as the means of $(\epsilon_1, \epsilon_2) $ even though the question stated they were zero.

Write the standard normal PDF and CDF by $\phi$ and $\Phi$.

Define $$Pr(\epsilon_2 > \alpha) = 1 - F_2 (\alpha) $$ and multiply and divide by this term:

$$ f (\epsilon_1 | \epsilon_2 > \alpha) = {1 \over 1 - F_2 (\alpha) } f(\epsilon_1 | \epsilon_2 > \alpha) Pr (\epsilon_2 > \alpha )$$

Then: $$ = {1 \over 1 - F_2 (\alpha) } \int_\alpha^\infty f(\epsilon_1, \epsilon_2) d \epsilon_2 $$

$$ = {1 \over 1 - F_2 (\alpha) } \int_\alpha^\infty f(\epsilon_2 | \epsilon_1) f_1 (\epsilon_1) d \epsilon_2 $$

$$ = { f_1 (\epsilon_1) \over 1 - F_2 (\alpha) } \int_\alpha^\infty f(\epsilon_2 | \epsilon_1 ) d \epsilon_2 $$

Focus now on the integral. $f(\epsilon_2 | \epsilon_1 ) $ follows a normal distribution too.

$$ \int_\alpha^\infty f(\epsilon_2 | \epsilon_1 ) d \epsilon_2 = 1 - F_{\epsilon_2 | \epsilon_1} (\alpha) $$

Now that we have the whole expression we just need to convert things into standard normal PDFs and CDFs.

We can write $ f_1 (\epsilon_1) = {1 \over \sigma_{1}} \phi({\epsilon_1 - \mu_1 \over \sigma_1})$ and $F_2 (\alpha) = \Phi ( {\alpha - \mu_2 \over \sigma_2 })$

Now note that $ f(\epsilon_2 | \epsilon_1 )$ is $N(\mu_{2|1}, \sigma_{2|1}) $, where $\mu_{2|1} = \mu_2 + \rho \frac{\sigma_2}{\sigma_1} (\epsilon_1 - \mu_1)$ and $ \sigma_{2|1} = (1 - \rho^2)\sigma_2^2 $, where $\rho$ is the correlation between $\epsilon_1$ and $\epsilon_2$ (check this entry in wikipedia)

Thus, the final expression is as follows:

$$ f(\epsilon_1 | \epsilon_2 > \alpha ) = { \phi\left({\epsilon_1 - \mu_1 \over \sigma_1}\right) \over \sigma_{1} \left( 1 - \Phi \left( {\alpha - \mu_2 \over \sigma_2 }\right) \right)} \left( 1 -\Phi\left( {\alpha - \mu_{2|1} \over \sigma_{2|1}} \right) \right) $$

I might add this result to the bivariate normal entry in wikipedia so that the next poor soul like me that has this doubt can easily find it. For that to happen, please help me check that the steps I followed are correct.