$\DeclareMathOperator\erf{erf}$ Let's consider the Gaussian distribution $P_X(x)= \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{x^2}{2 \sigma^2}}$. Now consider the random variable $W \equiv \max \{ X_1, \dots , X_m \}$, where $X_i \sim P_X(x)$. So $w$ indicates the maximum over $m$ i.i.d. Gaussian draws; we call $P_W(w)$ the corresponding distribution. $P_W(w)$ can be easily calculated from the formula $P_Y(y)=mP_X(y)F_X(y)^{m-1}$, where $F_X(y)$ is the c.d.f. of $P_X(x)$. We get:

$$P_Y(y)= m \frac{1}{\sqrt{2 \pi \sigma^2}} e^{- \frac{y^2}{2 \sigma^2}} \left(\frac{1}{2} \left( 1 + \erf \left( \frac{y}{\sqrt{2} \sigma}\right) \right) \right)^{m -1}$$

I'm interested in the convolution between $P_X(x)$ and $P_W(w)$, i.e.:

$$I(\Delta) = \int_{-\infty}^\infty P_X(x') P_W(x'+ \Delta) \, dx'$$

Is it possible to get an analytical closed form or even an approximated solution for $I(\Delta)$?