An integral involving a Gaussian and a power of a normal cumulative distribution function

Question

Being inspired by How to evaluate $\int_0^\infty\operatorname{erfc}^n x\ \mathrm dx$? we formulated the question below.

---

Let $c \in (0,1/\sqrt{2})$ and let $n \in \Bbb{N}$. Then let $\phi(x) : =\frac{\exp(-x^2/2)}{\sqrt{2\pi}}$ be the normal probability density function and let $\Phi(x) := \frac{1}{2} \operatorname{erfc}\left(\frac{-x}{\sqrt{2}}\right)$ be the standard normal cumulative density function.

In addition define:

\begin{eqnarray} {\mathfrak F}^{(A,B)}_{a,b} &:=& \int\limits_A^B \frac{\log(z+a)}{z+b} dz\\ &=& F[B,a,b] - F[A,a,b] + 1_{t^* \in (0,1)} \left( -F[A+(t^*+\epsilon)(B-A),a,b] + F[A+(t^*-\epsilon)(B-A),a,b] \right) \end{eqnarray} where \begin{eqnarray} t^*:=-\frac{\operatorname{Im}[(A+b)(b^*-a^*)]}{\operatorname{Im}[(B-A)(b^*-a^*)]} \end{eqnarray} and \begin{equation} F[z,a,b] := \log(z+a) \log\left( \frac{z+b}{b-a}\right) + \operatorname{Li}_2\left( \frac{z+a}{a-b}\right) \end{equation}

---

Now, consider the following integral:

\begin{equation} f_n(c) := \int\limits_{-\infty}^0 \phi(\xi) \cdot [\Phi(c \cdot \xi)]^n d\xi \end{equation}

By differentiating the function above with respect to the parameter $c$ , then completing the exponential-s to a square and then integrating by parts we have derived the following recurrence relation:

\begin{equation} \frac{d f_n(c)}{d c} = - \frac{n}{(1+c^2) 2 \pi} \cdot \left[ \frac{1}{2^{n-1}} - \frac{c \cdot (n-1)}{\sqrt{1+2 c^2}} \cdot f_{n-2} (\frac{c}{\sqrt{1+2 c^2}}) \right] \quad (ii) \end{equation}

The values of $f_n(c)$ for $n=1,2$ are known and given in https://en.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions for example. Here we found the result for $n=3$. By integrating the identity $(ii)$ with respect to $c$ and taking into account the initial condition $f_n(0) = 1/2^{n+1} $ we have shown that:

\begin{equation} f_3(c) = \frac{\pi -6 \tan ^{-1}(c)}{16 \pi } - \frac{3 i}{8 \pi ^2} \sum\limits_{\xi_1={-1,1}} \sum\limits_{\xi_2={-1,1}} \sum\limits_{\eta_1={-1,1}} \sum\limits_{\eta_2={-1,1}} \frac{\xi_1 \xi_2}{\imath \eta_2} {\mathfrak F}^{1,\frac{1+i \sqrt{2} c}{\sqrt{2 c^2+1}}}_{\xi_1 \sqrt{2} + \eta_1 \sqrt{3},\frac{\xi_2+\sqrt{2}}{\eta_2 \imath}} \end{equation}

Please take a glimpse at the Mathematica code below to see the consecutive steps in the derivation of that result:

In[1292]:= n = 3; c = 
 RandomReal[{0, 1/Sqrt[2]}, WorkingPrecision -> 50]; mPrec = 30;
NIntegrate[
 PDF[NormalDistribution[], xi] CDF[NormalDistribution[], c xi]^
   n, {xi, -Infinity, 0}, WorkingPrecision -> mPrec]
1/2^(n + 1) + (-n)/(2 Pi) 1/2^(n - 1) ArcTan[c] + -(1/
   32) (-1 + n) n + ((-1 + n) n ArcTan[Sqrt[-1 + 2 (1 + c^2)]])/(
 8 \[Pi]) - 
 n (n - 1)/(2 Pi)^2 NIntegrate[ 
   xi/(1 + xi^2) 1/
     Sqrt[1 + 2 xi^2] (ArcTan[xi/Sqrt[1 + 2 xi^2]]), {xi, 0, c}, 
   WorkingPrecision -> mPrec]
1/2^(n + 1) + (-n)/(2 Pi) 1/2^(n - 1) ArcTan[c] + -(1/
   32) (-1 + n) n + ((-1 + n) n ArcTan[Sqrt[-1 + 2 (1 + c^2)]])/(
 8 \[Pi]) - 
 n (n - 1)/(2 Pi)^2 NIntegrate[ (u Sqrt[(1/(1 - 2 u^2))] )/(1 - u^2)
     ArcTan[u], {u, 0, c/Sqrt[1 + 2 c^2]}, WorkingPrecision -> mPrec]
(2^(-5 - n) (16 \[Pi] - 32 n ArcTan[c] - 
    2^n (-1 + n) n (\[Pi] - 4 ArcTan[Sqrt[1 + 2 c^2]])))/\[Pi] - 
 n (n - 1)/(2 Pi)^2 NIntegrate[ (ArcTan[Sin[phi]/Sqrt[2]] Sin[phi] )/(
   2 (1 - Sin[phi]^2/2)), {phi, 0, ArcSin[Sqrt[2] c/Sqrt[1 + 2 c^2]]},
    WorkingPrecision -> mPrec]
(2^(-5 - n) (16 \[Pi] - 32 n ArcTan[c] - 
    2^n (-1 + n) n (\[Pi] - 4 ArcTan[Sqrt[1 + 2 c^2]])))/\[Pi] - 
 n (n - 1) I/(2 Pi)^2 NIntegrate[  (-1 + z^2)/(1 + 6 z^2 + z^4)
     Log[(-Sqrt[2] + 4 z + Sqrt[2] z^2)/(
     Sqrt[2] + 4 z - Sqrt[2] z^2)], {z, 1, (1 + I Sqrt[2] c)/Sqrt[
    1 + 2 c^2]}, WorkingPrecision -> mPrec]
(2^(-1 - n) (\[Pi] - 2 n ArcTan[c]))/\[Pi] - 
 n (n - 1) I/(2 Pi)^2 1/
   4 Sum[(xi2 xi1)/(eta2 I)
     FF[1, (1 + I Sqrt[2] c)/Sqrt[1 + 2 c^2], 
     xi1 Sqrt[2] + eta1 Sqrt[3], (xi2 + Sqrt[2])/(eta2 I)], {xi2, -1, 
    1, 2}, {eta2, -1, 1, 2}, {xi1, -1, 1, 2}, {eta1, -1, 1, 2}]



Out[1293]= 0.0259652081311838902701876114966

Out[1294]= 0.02596520813118389027018761149565

Out[1295]= 0.02596520813118389027018761149565

Out[1296]= 0.02596520813118389027018761149565

Out[1297]= 0.02596520813118389027018761149565 + 
 2.66030783426529888023311493014*10^-44 I

Out[1298]= 0.0259652081311838902701876114955083036534738420255 - 
 1.077372575572101842*10^-31 I

Update: Using basically the same methods as outlined above we derived the result for $n=4$. It reads:

\begin{eqnarray} &&f_4(c) -1/2^5 = -\frac{\tan ^{-1}(c)}{4 \pi } + \\ && \frac{3}{4 \pi^2} \sum\limits_{\xi_1={-1,1}} \sum\limits_{\xi_2={-1,1}} \sum\limits_{\eta_1={-1,1}} \sum\limits_{\eta_2={-1,1}} \xi_1 \xi_2 \left({\mathfrak F}^{1,\frac{1+i \sqrt{2} c}{\sqrt{2 c^2+1}}}_{\xi_2 \sqrt{2}+\eta_2 \sqrt{3},\imath(\xi_1+\eta_1 \sqrt{2})} - {\mathfrak F}^{(-1)^{1/4},\frac{1+i \sqrt{4 c^2+1}}{\sqrt{4 c^2+2}}}_{\xi_2/\sqrt{2} + \eta_2 \sqrt{3/2},\xi_1\eta_1 \imath \sqrt{2+\eta_1 \sqrt{3}}} \right) \end{eqnarray}

Again, the derivation of this result is in the Mathematica code below:

n = 4; c = 
 RandomReal[{0, 1/Sqrt[2]}, WorkingPrecision -> 50]; mPrec = 30;
(*
cc = Array[#&,100,{c-1/2,c+1/2}];
f=Interpolation[{#,NIntegrate[PDF[NormalDistribution[],xi] \
CDF[NormalDistribution[],# xi]^n,{xi,-Infinity,0},WorkingPrecision\
\[Rule]mPrec]}&/@cc];
f'[c]
-n/(2 Pi) 1/(1+c^2) (1/2^(n-1) - (n-1) c/Sqrt[1+2 c^2] 1/(2 \
Pi)(-ArcTan[c/Sqrt[1+2 c^2]] + ArcTan[Sqrt[1+2 c^2/(1+2 c^2)]]))
*)

NIntegrate[
  PDF[NormalDistribution[], xi] CDF[NormalDistribution[], c xi]^
    n, {xi, -Infinity, 0}, WorkingPrecision -> mPrec] - 1/2^(n + 1)

-n/(2 Pi) NIntegrate[
   1/(1 + xi^2) (1/
       2^(n - 1) - (n - 1) xi/
        Sqrt[1 + 2 xi^2] 1/(2 Pi) (-ArcTan[xi/Sqrt[1 + 2 xi^2]] + 
         ArcTan[Sqrt[1 + 2 xi^2/(1 + 2 xi^2)]])), {xi, 0, c}, 
   WorkingPrecision -> mPrec]
-(ArcTan[c]/(4 \[Pi])) - 
 n/(2 Pi) NIntegrate[((-1 + n) xi (ArcTan[xi/Sqrt[1 + 2 xi^2]]))/(
   2 \[Pi] (1 + xi^2) Sqrt[1 + 2 xi^2]), {xi, 0, c}, 
   WorkingPrecision -> mPrec] - 
 n/(2 Pi) NIntegrate[((-1 + 
      n) xi (-ArcTan[Sqrt[1 + (2 xi^2)/(1 + 2 xi^2)]]))/(
   2 \[Pi] (1 + xi^2) Sqrt[1 + 2 xi^2]), {xi, 0, c}, 
   WorkingPrecision -> mPrec]
-(ArcTan[c]/(4 \[Pi])) - 
 n/(2 Pi) NIntegrate[((-1 + n) u ArcTan[u])/(
   2 \[Pi] Sqrt[1 - 2 u^2] (1 - u^2)), {u, 0, c/Sqrt[1 + 2 c^2]}, 
   WorkingPrecision -> mPrec] - 
 n/(2 Pi) NIntegrate[((-1 + n) u ArcTan[u])/(
   2 \[Pi] Sqrt[2 - u^2] (-3 + u^2)), {u, 1, Sqrt[
    1 + (2 c^2)/(1 + 2 c^2)]}, WorkingPrecision -> mPrec]
-(ArcTan[c]/(4 \[Pi])) - 
 n (-1 + n) /(2 Pi) NIntegrate[(
   ArcTan[Sin[phi]/Sqrt[2]] Sin[phi])/(\[Pi] (3 + Cos[2 phi])), {phi, 
    0, ArcSin[(c Sqrt[2])/Sqrt[1 + 2 c^2]]}, 
   WorkingPrecision -> mPrec] - 
 n (-1 + n)/(2 Pi) NIntegrate[( 
   ArcTan[Sqrt[2] Sin[phi]] Sin[phi])/(\[Pi] Sqrt[
    2] (-3 + 2 Sin[phi]^2)), {phi, ArcSin[1/Sqrt[2]], 
    ArcSin[1/Sqrt[2] Sqrt[1 + (2 c^2)/(1 + 2 c^2)]]}, 
   WorkingPrecision -> mPrec]

-(ArcTan[c]/(4 \[Pi])) - 
 n (-1 + n) /(2 Pi) NIntegrate[(
   ArcTan[Sin[phi]/Sqrt[2]] Sin[
     phi])/(\[Pi] (3 + 1 - 2 Sin[phi]^2)), {phi, 0, 
    ArcSin[(c Sqrt[2])/Sqrt[1 + 2 c^2]]}, WorkingPrecision -> mPrec] -
  n (-1 + n)/(2 Pi) NIntegrate[( 
   ArcTan[Sqrt[2] Sin[phi]] Sin[phi])/(\[Pi] Sqrt[
    2] (-3 + 2 Sin[phi]^2)), {phi, ArcSin[1/Sqrt[2]], 
    ArcSin[1/Sqrt[2] Sqrt[1 + (2 c^2)/(1 + 2 c^2)]]}, 
   WorkingPrecision -> mPrec]

-((ArcTan[c] + 
   6 (ArcTan[Sqrt[1 + 2 c^2]] + 
      ArcTan[(1 + I Sqrt[1 + 4 c^2])/((1 + Sqrt[3]) Sqrt[
        1 + 2 c^2])] - 
      ArcTan[((1 + Sqrt[3]) (1 + I Sqrt[1 + 4 c^2]))/(
       2 Sqrt[1 + 2 c^2])]))/(4 \[Pi])) +
 n (-1 + n) /(16 Pi^2)  Sum[
   xi1 xi2 FF[1, (1 + I Sqrt[2] c)/Sqrt[1 + 2 c^2], 
     xi2 Sqrt[2] + eta2 Sqrt[3], I (xi1 + eta1 Sqrt[2])] , {xi2, -1, 
    1, 2}, {eta2, -1, 1, 2}, {xi1, -1, 1, 2}, {eta1, -1, 1, 
    2}] + -n (-1 + n)/(16 Pi^2)  Sum[
   xi1 xi2 FF[(-1)^(1/4), (1 + I Sqrt[1 + 4 c^2])/Sqrt[2 + 4 c^2], 
     xi2/Sqrt[2] + eta2 Sqrt[3/2], 
     xi1 eta1 I Sqrt[2 + eta1 Sqrt[3]]], {xi2, -1, 1, 2}, {eta2, -1, 
    1, 2}, {xi1, -1, 1, 2}, {eta1, -1, 1, 2}] 

Update 1: We obtained a partial answer in case $n=5$. It reads:

\begin{eqnarray} &&f_5(c) -1/2^6 = \frac{5 \left(-4 \tan ^{-1}\left(\sqrt{\frac{1}{2 c^2+1}}\right)-\tan ^{-1}(c)+\pi \right)}{32 \pi } + \\ && \frac{3(n)_{(2)}}{64 \pi^2} \sum _{\xi_1=-1,1} \sum _{\eta_1=-1,1} \sum _{\xi_2=-1,1} \sum _{\eta_2=-1,1} \eta_2 \xi_1 {\mathfrak F}^{0,\frac{\sqrt{2} c}{1+\sqrt{1+2 c^2}}}_{\frac{i \left(\sqrt{3} \eta_1+\xi_1\right)}{\sqrt{2}},\xi_2 (-1)^{\frac{\eta_2}{4}}} +\\ && -\frac{(n)_{(4)} \tan ^{-1}\left(\sqrt{2 c^2+1}\right)}{32 \pi ^3} \cdot \sum _{\xi_2=-1,1} \sum _{\eta_2=-1,1} \sum _{\xi_1=-1,1} \sum _{\eta_1=-1,1} \frac{(\xi_2 \xi_1) }{\eta_2} {\mathfrak F}^{1,\frac{i \sqrt{2} c+\sqrt{1+2 c^2}}{\sqrt{1+4 c^2}}}_{\xi_1 \sqrt{2}+\eta_1 \sqrt{3},\frac{\xi_2+\sqrt{2}}{\eta_2 i}} + \\ &&\frac{(n)_{(4)}}{8 \pi^3} \int\limits_0^c \frac{t \tan ^{-1}\left(\sqrt{2 t^2+1}\right) \left(\tan ^{-1}\left(\frac{t}{\sqrt{4 t^2+1}}\right)-\frac{\pi }{2}\right)}{\sqrt{2 t^2+1} \left(3 t^2+1\right) \sqrt{4 t^2+1}} dt \end{eqnarray}

n = 5; c = RandomReal[{0, 1/3}, WorkingPrecision -> 50]; mPrec = 40;
NIntegrate[
  PDF[NormalDistribution[], xi] CDF[NormalDistribution[], c xi]^
    n, {xi, -Infinity, 0}, WorkingPrecision -> mPrec] - 1/2^(n + 1)

(5 (\[Pi] - ArcTan[c] - 4 ArcTan[Sqrt[1/(1 + 2 c^2)]]))/(32 \[Pi]) +
 (3 (-1 + n) n)/(64 \[Pi]^2)
    Sum[eta2 xi1 FF[0, (Sqrt[2] c)/(1 + Sqrt[1 + 2 c^2]), (
     I (Sqrt[3] eta1 + xi1))/Sqrt[2], xi2 (-1)^(eta2 1/4)], {xi1, -1, 
    1, 2}, {eta1, -1, 1, 2}, {xi2, -1, 1, 2}, {eta2, -1, 1, 
    2}] + -(((-3 + n) (-2 + n) (-1 + n) n)/(32 \[Pi]^3)) ArcTan[Sqrt[
   1 + 2 c^2]] Sum[(xi2 xi1)/
    eta2  (FF[1, (I Sqrt[2] c + Sqrt[1 + 2 c^2])/Sqrt[1 + 4 c^2], 
      xi1 Sqrt[2] + eta1 Sqrt[3], (xi2 + Sqrt[2])/(
      eta2 I)]), {xi2, -1, 1, 2}, {eta2, -1, 1, 2}, {xi1, -1, 1, 
    2}, {eta1, -1, 1, 2}] + ((-3 + n) (-2 + n) (-1 + n) n)/(
  32 \[Pi]^3)
   4 (NIntegrate[
    t/(Sqrt[1 + 2 t^2] (1 + 3 t^2) Sqrt[
      1 + 4 t^2]) ((-Pi/2 + ArcTan[ t/Sqrt[1 + 4 t^2]]) ArcTan[Sqrt[
        1 + 2 t^2]]), {t, 0, c}, WorkingPrecision -> mPrec])

Unfortunately we are not able to evaluate the remaining integral in the right hand side.

Having said all this my question is the following. What is the result like for $n > 5$?