Calculate variance on the ratio of 2 angular power spectra

• A
Summary:
I am looking for a calculation about the variance on the ratio between 2 angular power spectra.
In the context of Survey of Dark energy stage IV, I need to evaluate the error on a new observable called "O" which is equal to :

\begin{equation}
O=\left(\frac{C_{\ell, \mathrm{gal}, \mathrm{sp}}^{\prime}}{C_{\ell, \mathrm{gal}, \mathrm{ph}}^{\prime}}\right)=\left(\frac{b_{s p}}{b_{p h}}\right)^{2}
\end{equation}

with spectroscopic ##C_{\ell, \text { gal,sp }}^{\prime}## defined by the integral :
$$C_{\ell, \mathrm{gal}, \mathrm{sp}}^{\prime}=\int_{l_{m i n}}^{l_{\max }} C_{\ell, \mathrm{gal}, \mathrm{ph}}(\ell) \mathrm{d} \ell=b_{s p}^{2} \int_{l_{\min }}^{l_{\max }} C_{\ell, \mathrm{DM}}(\ell) \mathrm{d} \ell=b_{s p}^{2} C_{\ell, \mathrm{DM}}^{\prime}$$
Same for photometric $C_{\ell, \text { gal,ph }}^{\prime}$ :
$$C_{\ell, \mathrm{gal}, \mathrm{ph}}^{\prime}=\int_{l_{m i n}}^{l_{\max }} C_{\ell, \mathrm{gal}, \mathrm{ph}}(\ell) \mathrm{d} \ell=b_{p h}^{2} \int_{l_{m i n}}^{l_{\max }} C_{\ell, \mathrm{DM}} \mathrm{d} \ell=b_{p h}^{2} C_{\ell, \mathrm{DM}}^{\prime}$$
We can infer this ratio of bias from the following definitions :
$$\begin{array}{l} C_{\ell, \mathrm{gal}, \mathrm{sp}}^{\prime}(k, z)=b_{s p}^{2} C_{\ell, \mathrm{DM}}^{\prime}(k, z) \\ C_{\ell, \mathrm{gal}, \mathrm{ph}}^{\prime}(k, z)=b_{p h}^{2} C_{\ell, \mathrm{DM}}^{\prime}(k, z) \end{array}$$

I need to estimate the variance ##\sigma_o^{2}## of this new observable. How could I begin a such task ?

For convenience, I am going to begin with ony one bin, so the integrals would disappear.

The ##C_l^{GG}## is defined by (by taking ##\gamma = G##

$$C_{i j}^{\gamma \gamma}(\ell) \simeq \frac{c}{H_{0}} \int \mathrm{d} z \frac{W_{i}^{\gamma}(z) W_{j}^{\gamma}(z)}{E(z) r^{2}(z)} P_{\delta \delta}\left[\frac{\ell+1 / 2}{r(z)}, z\right]$$
where ##i## and ##j## identify pairs of redshift bins, ##E(z)## is the dimensionless Hubble parameter of Eq. (11), ##r(z)## is the comoving distance, ##P_{\delta \delta}(k, z)## is the matter power spectrum evaluated at ##k=k_{\ell}(z) \equiv(\ell+1 / 2) / r(z)## due to the Limber approximation; we define the weight function $W^{\gamma}(z)$ as
\begin{aligned} W_{i}^{\gamma}(z) &=\frac{3}{2} \frac{H_{0}}{c} \Omega_{\mathrm{m}, 0}(1+z) \tilde{r}(z) \int_{z}^{z_{\max }} \mathrm{d} z^{\prime} n_{i}\left(z^{\prime}\right)\left[1-\frac{\tilde{r}(z)}{\tilde{r}\left(z^{\prime}\right)}\right] \\ &=\frac{3}{2} \frac{H_{0}}{c} \Omega_{\mathrm{m}, 0}(1+z) \tilde{r}(z) \widetilde{W}_{i}(z) \end{aligned}
with ##z_{\max }## the maximum redshift of the source redshift distribution. With respect to the standard formalism, we replace the comoving distance ##r(z)## with its dimensionless scaled version ##\tilde{r}(z)=r(z) /\left(c / H_{0}\right)## to highlight that the dependence of ##W_{i}^{\gamma}(z)## on the cosmological.

So, as you can see, this is not going to be easy to compute the variance on the ration between ##C\ell^{GG}##.

In the paper https://arxiv.org/abs/1910.09273, they simply put the following variance without demonstration :

\begin{equation}
\Delta C_{i j}^{\epsilon \epsilon}(\ell)=\sqrt{\frac{2}{(2 \ell+1) \Delta \ell f_{\mathrm{sky}}}} C_{i j}^{\epsilon \epsilon}(\ell)
\end{equation}

That's why any suggestion to start is welcome.