Integrate $q (1 - (1 - p)^6)$ over bivariate normal distribution?

Question

I have a function $F(p,q) = q(1 - (1 - p)^6)$ which tells probability of occurence based on probabilities $p, q$; and then I have estimates (a normal distribution) of logistic regression for $x_1 = \Phi^{-1}(p)$ and $x_2 = \Phi^{-1}(q)$ ($\Phi$ is the inverse link function, see below).

I want to estimate the mean of $F(p,q)$, and also its Confindence Interval.

For calculating the mean, I guess the following integral is needed? Is it possible to calculate it analytically?

$$\int^\infty_{-\infty} \int^\infty_{-\infty} \Phi(x_2) (1 - (1 - \Phi(x_1))^6) \cdot \cal{N}(x_{1,2}|\mu_{1,2}, \Sigma) dx_1 dx_2$$

Where $\Sigma$ is a 2x2 covariance matrix and $\Phi$ is the inverse link function, either:

1. Normal CDF with zero mean and sd = 1
2. Inverse logit

I don't actually care much which sigmoidal function will be used; I prefer logit, which I currently use; but if the integral is not tractable analytically for logit, I can switch to Normal CDF.