Let $\Sigma$ denote the covariance matrix of our asset universe, $\mu$ is the vector of expected returns. Further, $\mathbb{1}$ is a vector of ones. Let's identify the vector of the minimum variance portfolio's asset weights with $w_0$ and that of the tangency portfolio with $w_m$. Let's further identify
$$
\begin{align}
a&\equiv \mathbb{1}^T\Sigma^{-1}\mathbb{1}\\
b&\equiv \mathbb{1}^T\Sigma^{-1}\mathbb{\mu}\\
c&\equiv \mathbb{\mu}^T\Sigma^{-1}\mathbb{\mu}\\
d&=ac-b^2
\end{align}
$$
Canonically, the minimum variance portfolio's optimal weights are
$$
w_0=\frac{\Sigma^{-1}\mathbb{1}}{\mathbb{1}^T\Sigma^{-1}\mathbb{1}}=\frac{\Sigma^{-1}\mathbb{1}}{a}
$$
Its expected return is $\mathrm{E}(R_0)=w_0^T\mu=b/a$, its variance is $\mathrm{Var}(R_0)\equiv \sigma_0^2=w_0^T\Sigma w_0=1/a$.
The weights of the tangency portfolio are
$$
w_m=\frac{\Sigma^{-1}\mathbb{\mu}}{\mathbb{1}^T\Sigma^{-1}\mathbb{\mu}}=\frac{\Sigma^{-1}\mathbb{\mu}}{b}
$$
Its expected return is $\mathrm{E}(R_m)=w_m^T\mu=c/b$, its variance is $\mathrm{Var}(R_m)\equiv \sigma_m^2=w_m^T\Sigma w_m=c/b^2$. The covariance between the two is $\sigma_{m,0}=1/a=\sigma_0^2$
For any portfolio on the efficient frontier, $R_i$, its expected return is a combination of these two portfolios (or any other two portfolios):
$$
\begin{align}
\mathrm{E}(R_i)&=\alpha_i\mathrm{E}(R_m)+(1-\alpha_i)\mathrm{E}(R_0)\\
\Rightarrow\qquad \alpha_i&=\frac{\mathrm{E}(R_i)-\mu_0}{\mu_m-\mu_0}\\
&=\frac{\mathrm{E}(R_i)-b/a}{c/b-b/a}\\
&=\frac{ab}{d} \left(\mathrm{E}(R_i)-b/a\right)
\end{align}
$$
We can now calculate the covariance as
$$
\begin{align}
\mathrm{Cov}(R_i,R_j)&=\left[\alpha_i w_m+\left(1-\alpha_i\right)w_0\right]\Sigma\left[\alpha_j w_m+\left(1-\alpha_j\right)w_0\right]\\
&=\alpha_i\alpha_j w_m^T\Sigma w_m\\
&+\alpha_i(1-\alpha_j)w_m^T\Sigma w_0\\
&+(1-\alpha_i)\alpha_j w_0^T\Sigma w_m\\
&+(1-\alpha_i)(1-\alpha_j)w_0^T\Sigma w_0\
\end{align}
$$
From here on, we can plug in the values for $\alpha_i,\alpha_j$, substitute the (co-)variances and arrive at
$$
\begin{align}
\mathrm{Cov}(R_i,R_j)&=\sigma_0^2+\alpha_i\alpha_j (\sigma_m^2-\sigma_0^2)\\
&=\frac{1}{a}+\frac{a^2b^2}{d^2}\left(\frac{c}{b^2}-\frac{1}{a}\right)\left(\mathrm{E}(R_i)-b/a\right)\left(\mathrm{E}(R_j)-b/a\right)\\
&=\frac{1}{a}+\frac{a}{d}\left(\mathrm{E}(R_i)-\frac{b}{a}\right)\left(\mathrm{E}(R_j)-\frac{b}{a}\right)
\end{align}
$$