Let $P = N(\vec{0}, I^d)$ be a standard multivariate Gaussian distribution in $d$ dimensions. Let $Q$ be distributed the same as $P$, except that samples from $Q$ have one of their coordinates, chosen uniformly at random from $1$ to $d$, distributed as $N(\mu, \sigma^2)$ instead.

That is, $Q$ is a mixture of $d$ multivariate Gaussians, each with weight $\frac1d$. The $i$th mixture component has mean $\vec{0} + \vec{e}_i \cdot \mu$ and variance with diagonal $\vec{1} + \vec{e}_i \cdot (\sigma^2 - 1)$.

I'd like to bound the total variation distance between $P$ and $Q$, as a function of $d$ (when $\mu$ and $\sigma$ are constants that do not depend on $d$).

My intuition is that $Q$ is very close to a multivariate Gaussian $N(\frac{\vec{\mu}}{d}, I^d \cdot (1+\frac{\sigma^2}{d}+ o(\frac1d))$, which would give a Kullback–Leibler divergence of $D_{KL}(P || Q) = O(\frac1d)$ and thus a total variation distance of $O(\frac{1}{\sqrt d})$.

The case where $\sigma^2=1$ might be easier and sufficient for what I need. Any suggestions on how I could bound the TV distance between P and Q?