Why does the Monte Carlo estimate not depend on the dimension The Monte Carlo Estimator for some event probability (e.g., for the "failure probability") is defined as follows:
$$
\hat\mu = 1/N \sum_{i=1}^N I(\boldsymbol{x}_i),
$$
where $\boldsymbol{x}_i \in \mathbb{R}^d$ denotes a $d$-dimensional sample, drawn from a joint pdf $f_X(\boldsymbol{x})$, $N$ denotes the number of samples, and $I(x_i)$ is the indicator function of some ("failure") event.
The variance of the estimator is proportional to:
$$
Var(\hat\mu) \propto \frac{1}{\sqrt{N}}.
$$
Could someone explain to me why this variance does not depend on the dimension $d$.
Does this really mean that, for a given $N$, the level of accuracy in 3 and 1000 dimensions is indifferent?

In this Monte-Carlo setup, the indicator function $I: \mathbb{R}^d \rightarrow \{ 0,1 \}$ is defined on the domain $\mathbb{R}^d$, so it already deals with the multivariate nature of the observations and reduces this to a binary output.  So, if we have IID observations $\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3, ...$ and we set $p_I \equiv \mathbb{P}(I(\mathbf{x}_i) = 1)$ then we have:
$$I(\mathbf{x}_i) \sim \text{IID Bern}(p_I).$$
This then gives the estimator variance:
$$\begin{align}
\mathbb{V}(\hat{\mu})
&= \mathbb{V} \Bigg( \frac{1}{N} \sum_{i=1}^N I(\mathbf{x}_i) \Bigg) \\[6pt]
&= \frac{1}{N^2} \sum_{i=1}^N \mathbb{V} ( I(\mathbf{x}_i) ) \\[6pt]
&= \frac{1}{N^2} \sum_{i=1}^N p_I (1-p_I) \\[6pt]
&= \frac{p_I (1-p_I)}{N}.  \\[6pt]
\end{align}$$
As you can see from the working, the reason that this variance does not depend on the dimension $d$ is that $I(\mathbf{x}_i)$ is a Bernoulli random variable irrespective of this dimension and so its variance is determined by $p_I$.  So yes, the level of accuracy of the Monte-Carlo estimator is the same whether you have 3 dimensions or 1000.  Note that the only thing that you are estimating is the probability of a particular event, which is not something that gets more complicated in higher dimensions.
