Let $X \in \mathbb{R}^{d}$ follows the standard Gaussian distribution $N(0, I\_d)$. Let $Y = \max\_{j\in[d] } X\_j$. It is not hard to see that 
\begin{align}
\mathbb{E}\left [ Y \cdot X\right] = \sum\_{j=1}^d \mathbb{P}\left( j = \arg\max\_{i \in [d]} X\_i \right) \cdot e\_j,
\end{align}
where $e\_i$ is the standard basis in $\mathbb{R}^d$. Now I was wondering how to compute $$\mathbb{E} \left[ Y\cdot (X X^\top - I\_d) \right ].$$
Is there a closed form solution?



Expectation of max of Gaussian multiplied by a functional of Gaussian
Let $n:=d$. Assume $n\ge3$. Let $F$ and $f$ denote the cdf and pdf of $N(0,1)$, and let $I\{\cdot\}$ denote the indicator function. 
Each of the diagonal entries of the matrix $EYXX^T$ equals
\begin{equation\*}
 EX\_1^2\max\_iX\_i=E\_1+E\_2,
\end{equation\*}
where 
\begin{equation\*}
 E\_1:=EX\_1^3\,I\{X\_1>\max\_2^n X\_i\},\quad E\_2:=(n-1)EX\_1^2 X\_2\,I\{X\_2>X\_1\vee\max\_3^n X\_i\}. 
\end{equation\*}
Note that the cdf of the maximum of $k$ iid $N(0,1)$ random variables (r.v.'s) is $F^k$. So,
\begin{multline\*}
 E\_1=\int\_{-\infty}^\infty d(F(c)^{n-1})\,EX\_1^3\,I\{X\_1>c\}
 =\int\_{-\infty}^\infty d(F(c)^{n-1})\,\int\_c^\infty dx\,x^3f(x) \\ 
 =\int\_{-\infty}^\infty dc\, F(c)^{n-1}c^3f(c), 
\end{multline\*}
by integration by parts. Next, using the identities 
$$EX\_2\,I\{X\_2>u\}=f(u)$$ 
and 
$$d\_c f(u\vee c)=-cf(c)I\{c>u\}\,dc$$ 
for real $u$, and integrating by parts (twice), we have 
\begin{multline\*}
 E\_2:=\int\_{-\infty}^\infty (n-1)d(F(c)^{n-2})\,EX\_1^2 X\_2\,I\{X\_2>X\_1\vee c\} \\
 =\int\_{-\infty}^\infty (n-1)d(F(c)^{n-2})\,EX\_1^2 f(X\_1\vee c) \\ 
 =\int\_{-\infty}^\infty (n-1)dc\, F(c)^{n-2}\,cf(c) EX\_1^2\,I\{X\_1<c\} \\
 =\int\_{-\infty}^\infty (n-1)dc\, F(c)^{n-2}\,cf(c) (F(c)-cf(c)) \\ 
 =\int\_{-\infty}^\infty (n-1)dc\, F(c)^{n-1}\,cf(c) 
 -\int\_{-\infty}^\infty (n-1)dc\, F(c)^{n-2}\,f(c)c^2f(c) \\ 
 =\int\_{-\infty}^\infty (n-1)dc\, F(c)^{n-1}\,cf(c) 
 +\int\_{-\infty}^\infty dc\, F(c)^{n-1}\,(2c-c^3)f(c) \\ 
 =\int\_{-\infty}^\infty dc\, F(c)^{n-1}\,((n+1)c-c^3)f(c). 
\end{multline\*}
Thus, each of the diagonal entries of the matrix $EYXX^T$ equals
\begin{equation\*}
 EX\_1^2\max\_iX\_i=E\_1+E\_2=(n+1)\int\_{-\infty}^\infty dc\, F(c)^{n-1}\,f(c)c
\end{equation\*}. 
Since 
\begin{equation\*}
 EY=E\max\_iX\_i=n\int\_{-\infty}^\infty dc\, F(c)^{n-1}\,f(c)c, \tag{1}
\end{equation\*}
we conclude that each of the diagonal entries of the matrix $EY(XX^T-I)$ equals
\begin{equation\*}
\text{diags}= EX\_1^2\max\_iX\_i=\int\_{-\infty}^\infty dc\, F(c)^{n-1}\,f(c)c=\tfrac1n\,EY. \tag{2}
\end{equation\*}


Each of the off-diagonal entries of the matrix $EY(XX^T-I)$ equals
\begin{equation\*}
 EX\_1X\_2\max\_iX\_i=\sum\_1^n s\_j, 
\end{equation\*}
where
\begin{equation\*}
 s\_j:=EX\_1 X\_2 X\_j\,I\{X\_j\ge \max\_1^n X\_i\}.
\end{equation\*}


If $j\notin\{1,2\}$, then, similarly to the previous multi-line display, 
\begin{multline\*}
 s\_j=EX\_1 X\_2 X\_3\,I\{X\_3>X\_1\vee X\_2\vee \max\_4^n X\_i\} \\ 
 =\int\_{-\infty}^\infty d(F(c)^{n-3})\,EX\_1 X\_2 X\_3\,I\{X\_3>X\_1\vee X\_2\vee c\} \\ 
=\int\_{-\infty}^\infty d(F(c)^{n-3})\,EX\_1 X\_2 f(X\_1\vee X\_2\vee c) \\ 
=\int\_{-\infty}^\infty dc\,F(c)^{n-3}f(c)c\,EX\_1 X\_2 I\{X\_1\vee X\_2<c\} \\ 
=\int\_{-\infty}^\infty dc\,F(c)^{n-3}f(c)c\,f(c)^2 \\ 
=\frac1{n-2}\,\int\_{-\infty}^\infty dc\,F(c)^{n-2}(2c^2-1)\,f(c)^2 \\ 
=\frac2{(n-1)(n-2)}\,\int\_{-\infty}^\infty dc\,F(c)^{n-1}(2c^3-5c)\,f(c). 
\end{multline\*} 


If $j\in\{1,2\}$, then 
\begin{multline\*}
 s\_j=EX\_1 X\_2^2\,I\{X\_2>X\_1\vee \max\_3^n X\_i\} \\ 
 =\int\_{-\infty}^\infty d(F(c)^{n-2})\,EX\_1 I\{X\_1<X\_2\}X\_2^2\,I\{X\_2>c\} \\ 
 =-\int\_{-\infty}^\infty d(F(c)^{n-2})\,Ef(X\_2)X\_2^2\,I\{X\_2>c\} \\ 
 =-\int\_{-\infty}^\infty d(F(c)^{n-2})\,\int\_c^\infty dx\, f(x)^2x^2 \\ 
 =-\int\_{-\infty}^\infty dc\,F(c)^{n-2}\, f(c)c^2 \\ 
 =\frac1{n-1}\,\int\_{-\infty}^\infty dc\, F(c)^{n-1}\,(2c-c^3)f(c). 
\end{multline\*}


So, each of the off-diagonal entries of the matrix $EY(XX^T-I)$ equals
\begin{multline\*}
 \text{off-diags}=EX\_1X\_2\max\_iX\_i=2s\_1+(n-2)s\_3=
 \frac2{n-1}\int\_{-\infty}^\infty dc\,F(c)^{n-1}(c^3-3c)\,f(c) \\
 =\frac2{(n-1)n}\,E(Y^3-3Y) \tag{3} 
\end{multline\*}
-- cf. $(1)$, $(2)$. 


Mathematica claims that the value of the latter integral is $-\frac{1}{4 \sqrt{\pi }}\approx -0.141$ for $n=2,3$, 
but it cannot evaluate it in closed form for $n=4$ (an approximate value of that integral for $n=4$ is $-0.0969$). 


Yet, it is now easy to compute the diagonal and off-diagonal entries of the matrix $EY(XX^T-I)$ numerically, and it is also easy to find their asymptotics for large $n(=d)$, using standard methods of asymptotic analysis. 


Indeed, let us consider the integral 
\begin{equation\*}
 J\_n:=\int\_{-\infty}^\infty dc\,F(c)^{n-1}(c^3-3c)\,f(c)
\end{equation\*}
in $(3)$ for large $n$. Write
\begin{multline\*}
 J\_n=J\_{1n}+J\_{2n}+J\_{3n},\\
 J\_{1n}:=\int\_{-\infty}^{c\_n}\dots,\ J\_{2n}:=\int\_{c\_n}^{d\_n}\dots,\ J\_{3n}:=\int\_{d\_n}^\infty\dots, 
\end{multline\*}
where 
\begin{gather\*}
 c\_n:=F^{-1}(e^{-w\_n/n}),\quad
 d\_n:=F^{-1}(e^{-v\_n/n}) \\ 
 v\_n:=n^{-t\_n},\quad w\_n:=n^{t\_n},\quad t\_n:=1/\sqrt{\ln n}. 
\end{gather\*}
Then 
\begin{gather\*}
 0<v\_n<w\_n,\quad v\_n\to0,\quad w\_n=1/v\_n\to\infty,\quad d\_n>c\_n\to\infty,\\
 1-F(c\_n)\sim-\ln F(c\_n)=w\_n/n=n^{-1+o(1)},\quad 1-F(c\_n)=\exp\{-c\_n^2/(2+o(1))\}. 
\end{gather\*}
Hence and from similar relations for $d\_n$, we find that $c\_n\sim\sqrt{2\ln n}\sim d\_n$, and so, 
\begin{equation\*}
c\sim\sqrt{2\ln n}\text{ uniformly in $c\in[c\_n,d\_n]$}. \tag{4}
\end{equation\*}


Next, for real $c<c\_n$ we have $F(c)^{n-1}<F(c\_n)^{n-1}=e^{-(n-1)w\_n/n}\sim e^{-w\_n}
=\exp\{-e^{\sqrt{\ln n}}\}<\exp\{-\ln^2 n\}=o(1/n)$, whence 
\begin{equation\*}
 |J\_{1n}|\le F(c\_n)^{n-1}\int\_{-\infty}^{c\_n} dc\,|c^3-3c|\,f(c)=o(1/n).
\end{equation\*}


On the other hand, 
\begin{equation\*}
 |J\_{3n}|\le \int\_{d\_n}^\infty dc\,|c^3-3c|\,f(c)
 \sim d\_n^3(1-F(d\_n))\sim-d\_n^3\ln F(d\_n)\sim (\sqrt{2\ln n})^3 v\_n/n
 =o(1/n).
\end{equation\*}


Consider finally the integral $J\_{2n}$, which is the main term of the asymptotics. In that integral, let us use the substitution $v=-n \ln F(c)$, so that the condition $c\_n<c<d\_n$ can be rewritten as $v\_n<v<w\_n$. 
Note also that $v=-n \ln F(c)\iff c=C\_n(v):=F^{-1}(e^{-v/n})$. 
Also, by $(4)$, $C\_n(v)\sim\sqrt{2\ln n}$ uniformly in $v\in[v\_n,w\_n]$. So, 
\begin{equation\*}
 J\_{2n}=\frac1n\,\int\_{v\_n}^{w\_n}dv\,e^{-v}[C\_n(v)^3-3C\_n(v)]
 \sim\frac{(\sqrt{2\ln n})^3}n. 
\end{equation\*}


Collecting all the pieces, we have 
\begin{equation\*}
 \text{off-diags}=\frac2{(n-1)n}\,E(Y^3-3Y)
 \sim\frac{2^{5/2}\ln^{3/2} n}{n^2}. 
\end{equation\*}


It is similar (and slightly simpler) to derive from $(2)$ the following:  

\begin{equation\*}
\text{diags}=\tfrac1n\,EY 
\sim\frac{2^{1/2}\ln^{1/2} n}{n}. 
\end{equation\*}
So, one may also notice that 
\begin{equation\*}
 \text{off-diags}\sim\tfrac2{n^2}\,(EY)^3. 
\end{equation\*} 


Curiously, the off-diagonal entries of the matrix $EY(XX^T-I)$ seem to be negative for $n\in\{2,\dots,6\}$ but 
positive for $n\ge7$. 


