First, using Lemmas 1-3 of [@song:schm:1995], $$\begin{aligned}
\mathbb{E}\left[\sum\limits_{k=1}^{m}\left(\hat{\mu}_k-\hat{\mu}\right)\left(\hat{\mu}_k-\hat{\mu}\right)^{T}\right] &= \mathbb{E}\left[\sum_{k=1}^{m}\hat{\mu}_k\hat{\mu}_k^{T} -\sum_{k=1}^{m}\hat{\mu}_k\hat{\mu}^{T} - \hat{\mu}\sum_{k=1}^{m}\hat{\mu}_k^{T} + m\hat{\mu}\hat{\mu}^{T}\right]\\
&= \mathbb{E}\left[\sum_{k=1}^{m}\hat{\mu}_k\hat{\mu}_k^{T} - m\hat{\mu}\hat{\mu}^{T}\right]\\
&= m\left(\mathbb{E}\left[\bar{Y}_{1}\bar{Y}_{1}^{T}\right] - \mathbb{E}\left[\hat{\mu}\hat{\mu}^{T}\right]\right)\\
&= m\left({\text{Var}}\left(\bar{Y}_{1}\right) - {\text{Var}}\left(\hat{\mu}\right)\right)\\
&= m\left[\dfrac{\Sigma}{n} + \dfrac{\Gamma}{n^{2}} + o\left(\dfrac{1}{n^{2}}\right) - \dfrac{\Sigma}{mn} - \dfrac{\Gamma}{mn^{2}} + o\left(\dfrac{1}{mn^{2}}\right)\right]\\
&= m\left[\dfrac{\Sigma}{n}\left(1-\dfrac{1}{m}\right) + \dfrac{\Gamma}{n^{2}}\left(1-\dfrac{1}{m}\right) + o\left(\dfrac{1}{n^{2}}\right)\right]\,. {\addtocounter{equation}{1}\tag{\theequation}}\label{eq:expect_y_diff}\end{aligned}$$ By , , and the bias of the BM estimator from @vats:flegal:2018 [Theorem 1], $$\begin{aligned}
\mathbb{E}[\hat{\Sigma}_{R,b}] &= \mathbb{E}\left[\dfrac{a-1}{am-1}\sum_{k=1}^{m}\hat{\Sigma}_{k,b} + \dfrac{ab}{am-1}\sum_{k=1}^{m}\left(\hat{\mu}_k-\hat{\mu}\right)\left(\hat{\mu}_k-\hat{\mu}\right)^{T}\right]\\
&= \dfrac{m\left(a-1\right)}{am-1}\mathbb{E}\left[\hat{\Sigma}_{k,1}\right] + \dfrac{abm}{am-1}\left(\dfrac{\Sigma}{n}\left(1-\dfrac{1}{m}\right) + \dfrac{\Gamma}{n^{2}}\left(1-\dfrac{1}{m}\right) + o\left(\dfrac{1}{n^{2}}\right)\right)\\
&= \dfrac{m\left(a-1\right)}{am-1}\left[\Sigma + \dfrac{\Gamma}{b} + o\left(\dfrac{1}{b}\right)\right] + \dfrac{m\Sigma}{am-1}\left(1-\dfrac{1}{m}\right) + \dfrac{m\Gamma}{n\left(am-1\right)}\left(1-\dfrac{1}{m}\right) + \dfrac{m}{am-1}o\left(\dfrac{1}{n}\right)\\
&= \dfrac{m\Sigma}{am-1}\left[\left(a-1\right)+1-\dfrac{1}{m}\right] + \dfrac{m\Gamma}{am-1}\left(\dfrac{a-1}{b} + \dfrac{1}{n}\left(1-\dfrac{1}{m}\right)\right) + o\left(\dfrac{1}{b}\right)\\
&= \Sigma + \dfrac{m\left(a-1\right)\Gamma}{\left(am-1\right)b} + \dfrac{m\Gamma}{n\left(am-1\right)}\left(1-\dfrac{1}{m}\right) + o\left(\dfrac{1}{b}\right)\\
&= \Sigma + \dfrac{\Gamma\left(am - 1 + 1 - m\right)}{\left(am - 1\right)b} + o\left(\dfrac{1}{b}\right)\\
&= \Sigma + \dfrac{\Gamma}{b} + o\left(\dfrac{1}{b}\right) \; .\end{aligned}$$

Variance of RBM {#subsec:variance_RBM}
---------------

Recall that $\hat{\Sigma}_{R,L}$ is the lugsail RBM estimator and $\hat{\Sigma}_{A,l}$ is the lugsail averaged batch means estimator. Further, let $\hat{\Sigma}^{ij}_{R,L}$ and $\hat{\Sigma}^{ij}_{A,L}$ be the $\left(i,j\right)$th element of the matrices, respectively. We will prove that the variance of both the estimators is equivalent for large sample sizes. Due to the strong consistency proof, as $n \to \infty$, $$\label{eq:rbm_abm_consis}
 \left|\left|\hat{\Sigma}_{R,L} - \dfrac{m\left(a-1\right)}{am - 1} \hat{\Sigma}_{A,l}\right|\right| \to 0 \text{ with probability 1}\,.$$ For ease 