[proofplan]
The identity estimator $x \mapsto x$ gives the upper bound by computing the second moment of a centred Gaussian error. For the lower bound, fix a prior scale parameter $\tau > 0$ and compare the minimax risk with the Bayes risk under the Gaussian prior $\mathcal N(0,\tau^2 I_d)$. The posterior mean is the linear shrinkage estimator
\begin{align*}
X \mapsto \frac{\tau^2}{\sigma^2+\tau^2}X,
\end{align*}
and its Bayes risk is
\begin{align*}
d\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}.
\end{align*}
Letting the prior scale parameter $\tau$ tend to infinity forces every minimax estimator to have worst-case risk at least $d\sigma^2$.
[/proofplan]
[step:Compute the risk of the identity estimator]
Define the measurable estimator $\delta_0: \mathbb R^d \to \mathbb R^d$ by
\begin{align*}
\delta_0(x)=x.
\end{align*}
For fixed $\theta \in \mathbb R^d$, write $Z := X-\theta$. Under $\mathbb P_\theta$, the random vector $Z: \Omega \to \mathbb R^d$ has distribution $\mathcal N(0,\sigma^2 I_d)$. If $Z_i$ denotes its $i$-th coordinate, then $\mathbb E_\theta[Z_i^2]=\sigma^2$ for each $i \in \{1,\dots,d\}$. Therefore $|X-\theta|^2=\sum_{i=1}^d Z_i^2$, and linearity of expectation gives
\begin{align*}
R(\theta,\delta_0)=\mathbb E_\theta[|X-\theta|^2]=\mathbb E_\theta\left[\sum_{i=1}^d Z_i^2\right]=\sum_{i=1}^d \mathbb E_\theta[Z_i^2]=d\sigma^2.
\end{align*}
Taking the supremum over $\theta \in \mathbb R^d$ and then the infimum over all measurable estimators gives
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2) \le d\sigma^2.
\end{align*}
[/step]
[step:Lower-bound the minimax risk by Gaussian Bayes risks]
Fix $\tau > 0$. Let $\Pi_\tau$ be the probability measure $\mathcal N(0,\tau^2 I_d)$ on $\mathbb R^d$. For any measurable estimator $\delta: \mathbb R^d \to \mathbb R^d$,
\begin{align*}
\sup_{\theta \in \mathbb R^d} R(\theta,\delta)
\ge \int_{\mathbb R^d} R(\theta,\delta)\,d\Pi_\tau(\theta).
\end{align*}
Taking the infimum over $\delta$ yields
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2)
\ge
\inf_{\delta:\mathbb R^d \to \mathbb R^d \text{ measurable}}
\int_{\mathbb R^d} R(\theta,\delta)\,d\Pi_\tau(\theta).
\end{align*}
Equivalently, let $\Theta: \Omega \to \mathbb R^d$ have distribution $\Pi_\tau$, and conditionally on $\Theta=\theta$, let $X: \Omega \to \mathbb R^d$ have distribution $\mathcal N(\theta,\sigma^2 I_d)$. Since the loss $(\omega,\theta) \mapsto |\delta(X(\omega))-\theta|^2$ is nonnegative and measurable, [Tonelli's theorem](/page/Tonelli%20Theorem) justifies exchanging the prior integral with the [conditional expectation](/page/Conditional%20Expectation), even when the value is infinite. Hence
\begin{align*}
\int_{\mathbb R^d} R(\theta,\delta)\,d\Pi_\tau(\theta)=\mathbb E[|\delta(X)-\Theta|^2].
\end{align*}
Thus
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2)
\ge
\inf_{\delta:\mathbb R^d \to \mathbb R^d \text{ measurable}}
\mathbb E[|\delta(X)-\Theta|^2].
\end{align*}
[/step]
[step:Identify the posterior mean under the Gaussian prior]
Define the shrinkage map $m_\tau: \mathbb R^d \to \mathbb R^d$ by
\begin{align*}
m_\tau(x)=\frac{\tau^2}{\sigma^2+\tau^2}x.
\end{align*}
We claim that $m_\tau(X)=\mathbb E[\Theta \mid X]$.
Let $\mathcal L^d$ denote $d$-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb R^d$. Indeed, the joint density of $(\Theta,X)$ with respect to $\mathcal L^d \otimes \mathcal L^d$ is proportional to
\begin{align*}
(\theta,x)
\mapsto
\exp\left(
-\frac{|\theta|^2}{2\tau^2}
-\frac{|x-\theta|^2}{2\sigma^2}
\right).
\end{align*}
For fixed $x \in \mathbb R^d$, complete the square in $\theta$:
\begin{align*}
\frac{|\theta|^2}{\tau^2}+\frac{|x-\theta|^2}{\sigma^2}=\left(\frac{1}{\tau^2}+\frac{1}{\sigma^2}\right)\left|\theta-\frac{\tau^2}{\sigma^2+\tau^2}x\right|^2+\frac{|x|^2}{\sigma^2+\tau^2}.
\end{align*}
For this fixed $x$, the integral of the joint density over $\theta \in \mathbb R^d$ is finite and positive, so the conditional density of $\Theta$ given $X=x$ is obtained by normalising the displayed expression as a function of $\theta$.
Hence the conditional law of $\Theta$ given $X=x$ is the Gaussian law with mean
\begin{align*}
m_\tau(x)=\frac{\tau^2}{\sigma^2+\tau^2}x
\end{align*}
and covariance matrix
\begin{align*}
\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}I_d.
\end{align*}
Therefore $m_\tau(X)=\mathbb E[\Theta\mid X]$.
[guided]
The reason for introducing the Gaussian prior is that it makes the conditional distribution of $\Theta$ given $X$ computable. We define the map $m_\tau: \mathbb R^d \to \mathbb R^d$ by
\begin{align*}
m_\tau(x)=\frac{\tau^2}{\sigma^2+\tau^2}x.
\end{align*}
We now verify that this map is the posterior mean.
Let $\mathcal L^d$ denote $d$-dimensional Lebesgue measure on $\mathbb R^d$. The prior density of $\Theta$ is proportional to $\theta \mapsto \exp(-|\theta|^2/(2\tau^2))$, and the conditional density of $X$ given $\Theta=\theta$ is proportional to $x \mapsto \exp(-|x-\theta|^2/(2\sigma^2))$. Therefore the joint density of $(\Theta,X)$ with respect to $\mathcal L^d \otimes \mathcal L^d$ is proportional to
\begin{align*}
(\theta,x)
\mapsto
\exp\left(
-\frac{|\theta|^2}{2\tau^2}
-\frac{|x-\theta|^2}{2\sigma^2}
\right).
\end{align*}
For fixed $x \in \mathbb R^d$, the conditional density of $\Theta$ given $X=x$ is obtained by viewing this expression as a function of $\theta$. First expand the quadratic expression:
\begin{align*}
\frac{|\theta|^2}{\tau^2}+\frac{|x-\theta|^2}{\sigma^2}=\left(\frac{1}{\tau^2}+\frac{1}{\sigma^2}\right)|\theta|^2-\frac{2x\cdot\theta}{\sigma^2}+\frac{|x|^2}{\sigma^2}.
\end{align*}
Completing the square in $\theta$ gives
\begin{align*}
\frac{|\theta|^2}{\tau^2}+\frac{|x-\theta|^2}{\sigma^2}=\left(\frac{1}{\tau^2}+\frac{1}{\sigma^2}\right)\left|\theta-\frac{\tau^2}{\sigma^2+\tau^2}x\right|^2+\frac{|x|^2}{\sigma^2+\tau^2}.
\end{align*}
The final term $\frac{|x|^2}{\sigma^2+\tau^2}$ does not depend on $\theta$, so after normalising the conditional density it only affects the normalising constant. The normalisation is legitimate because, for each fixed $x \in \mathbb R^d$, the integral of the joint density over $\theta \in \mathbb R^d$ is a finite positive [Gaussian integral](/theorems/1140); equivalently, the marginal density of $X$ at $x$ is finite and positive. Thus the conditional law of $\Theta$ given $X=x$ is Gaussian with mean
\begin{align*}
\frac{\tau^2}{\sigma^2+\tau^2}x
\end{align*}
and covariance matrix
\begin{align*}
\left(\frac{1}{\tau^2}+\frac{1}{\sigma^2}\right)^{-1}I_d
=
\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}I_d.
\end{align*}
Therefore
\begin{align*}
\mathbb E[\Theta\mid X=x]=m_\tau(x),
\end{align*}
and hence $m_\tau(X)=\mathbb E[\Theta\mid X]$.
[/guided]
[/step]
[step:Compute the Gaussian Bayes risk]
Let $\delta: \mathbb R^d \to \mathbb R^d$ be any measurable estimator with finite Bayes risk. This finiteness implies that $\delta(X)-\Theta$ is square-integrable, so the [conditional expectation orthogonality identity](/page/Conditional%20Expectation%20Projection%20Identity) applies to the projection of $\Theta$ onto the $\sigma$-algebra generated by $X$. Since $m_\tau(X)=\mathbb E[\Theta\mid X]$, it gives
\begin{align*}
\mathbb E[|\delta(X)-\Theta|^2]=\mathbb E[|\delta(X)-m_\tau(X)|^2]+\mathbb E[|m_\tau(X)-\Theta|^2].
\end{align*}
The first term on the right is nonnegative, hence
\begin{align*}
\mathbb E[|\delta(X)-\Theta|^2]\ge\mathbb E[|m_\tau(X)-\Theta|^2].
\end{align*}
Estimators with infinite Bayes risk cannot decrease the infimum, so restricting first to finite Bayes risk loses no lower bound.
The conditional covariance matrix of $\Theta$ given $X$ is
\begin{align*}
\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}I_d.
\end{align*}
Therefore, by the [tower property of conditional expectation](/page/Tower%20Property%20of%20Conditional%20Expectation),
\begin{align*}
\mathbb E[|m_\tau(X)-\Theta|^2]=\mathbb E\left[\mathbb E\left[|\Theta-\mathbb E[\Theta\mid X]|^2\mid X\right]\right].
\end{align*}
The inner conditional expectation is the trace of the conditional covariance matrix, so
\begin{align*}
\mathbb E\left[\mathbb E\left[|\Theta-\mathbb E[\Theta\mid X]|^2\mid X\right]\right]=\operatorname{tr}\left(\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}I_d\right)=d\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}.
\end{align*}
Thus
\begin{align*}
\inf_{\delta:\mathbb R^d \to \mathbb R^d \text{ measurable}}
\mathbb E[|\delta(X)-\Theta|^2]
=
d\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}.
\end{align*}
[/step]
[step:Let the prior variance diverge to recover the unrestricted lower bound]
Combining the previous steps, for every $\tau > 0$,
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2)
\ge
d\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}.
\end{align*}
Because this lower bound holds for every $\tau > 0$, the minimax risk is bounded below by the supremum of the right-hand side over $\tau > 0$. Taking the limit as $\tau \to \infty$ gives
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2)\ge\sup_{\tau>0} d\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}=\lim_{\tau\to\infty} d\frac{\sigma^2\tau^2}{\sigma^2+\tau^2}=d\sigma^2.
\end{align*}
Together with the upper bound from the identity estimator,
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2) \le d\sigma^2,
\end{align*}
we obtain
\begin{align*}
\mathfrak M(\mathbb R^d,|\cdot|^2)=d\sigma^2.
\end{align*}
[/step]
