## Appendix A: Cross-correlation covariance

We generalize here the expression for the covariance of the cross-power spectrum of different tracers proposed by Smith (2009, Section 3, and Equation 43 specifically), accounting for the cross correlation of four tracers (mass bins in our case) A, B, C, and D. Similarly to Eq. (9), we define the covariance of the cross-power spectrum as

$$C_{ABCD}(\mathbf{k}_1, \mathbf{k}_2) = \left\langle \left[ \hat{P}_{AB}(\mathbf{k}_1) - \langle \hat{P}_{AB}(\mathbf{k}_1) \rangle \right] \left[ \hat{P}_{CD}(\mathbf{k}_2) - \langle \hat{P}_{CD}(\mathbf{k}_2) \rangle \right] \right\rangle. \quad (A.1)$$

Substituting the definition of the estimator for the halo cross-power spectrum,

$$\hat{P}_{AB}(\mathbf{k}_1) = V \langle \delta_A(\mathbf{k}_1) \delta_B(-\mathbf{k}_1) \rangle - \frac{1}{n}, \quad (A.2)$$

we obtain

$$C_{ABCD}(\mathbf{k}_1, \mathbf{k}_2) = V^2 \langle \delta_A(\mathbf{k}_1) \delta_B(-\mathbf{k}_1) \delta_C(\mathbf{k}_2) \delta_D(-\mathbf{k}_2) \rangle \\ - \left( P_{AB}(\mathbf{k}_1) + \frac{\delta_{AB}^K}{n_A} \right) \left( P_{CD}(\mathbf{k}_2) + \frac{\delta_{CD}^K}{n_C} \right). \quad (A.3)$$

The first term is the cross trispectrum,

$$V^2 \langle \delta_A(\mathbf{k}_1) \delta_B(-\mathbf{k}_1) \delta_C(\mathbf{k}_2) \delta_D(-\mathbf{k}_2) \rangle \\ = \frac{V^2}{N_A N_B N_C N_D} \sum_{ijkl} e^{i \mathbf{k}_1 \cdot \mathbf{r}_i - i \mathbf{k}_1 \cdot \mathbf{r}_j + i \mathbf{k}_2 \cdot \mathbf{r}_k - i \mathbf{k}_2 \cdot \mathbf{r}_l} \\ \times \left\langle \left( N_{Ai} - \langle N_{Ai} \rangle \right) \left( N_{Bj} - \langle N_{Bj} \rangle \right) \left( N_{Ck} - \langle N_{Ck} \rangle \right) \left( N_{Dl} - \langle N_{Dl} \rangle \right) \right\rangle, \quad (A.4)$$

where $N_{Ai}$ is the number of A-type halos (or halos in the Ath mass bin) at position vector $\mathbf{r}_i$.

Performing the calculations analogously to those shown in Sect. 3 of Smith (2009), we obtain the final expression for the cross-power spectrum covariance,

$$C_{ABCD}(\mathbf{k}_1, \mathbf{k}_2) = \frac{1}{V} T_{ABCD}(\mathbf{k}_1, -\mathbf{k}_1, \mathbf{k}_2, -\mathbf{k}_2) \\ + \left( P_{AC} + \frac{\delta_{AC}^K}{n_A} \right) \left( P_{BD} + \frac{\delta_{BD}^K}{n_B} \right) \delta^K(\mathbf{k}_1 + \mathbf{k}_2) \\ + \left( P_{AD} + \frac{\delta_{AD}^K}{n_A} \right) \left( P_{BC} + \frac{\delta_{BC}^K}{n_B} \right) \delta^K(\mathbf{k}_1 - \mathbf{k}_2) \\ + \frac{\delta_{AC}^K}{n_A V} B_{BDA}(\mathbf{k}_1, \mathbf{k}_2) + \frac{\delta_{AD}^K}{n_A V} B_{BCA}(\mathbf{k}_1, -\mathbf{k}_2) \\ + \frac{\delta_{BC}^K}{n_B V} B_{ADB}(\mathbf{k}_1, -\mathbf{k}_2) + \frac{\delta_{BD}^K}{n_B V} B_{ACB}(\mathbf{k}_1, \mathbf{k}_2) \\ + \frac{\delta_{CD}^K}{n_C V} B_{ABC}(\mathbf{k}_1, -\mathbf{k}_1, 0) + \frac{\delta_{AB}^K}{n_A V} B_{CDA}(\mathbf{k}_2, -\mathbf{k}_2, 0) \\ + \frac{\delta_{AC}^K \delta_{BD}^K}{n_A n_B V} P_{AB}(\mathbf{k}_1 + \mathbf{k}_2) + \frac{\delta_{AD}^K \delta_{BC}^K}{n_A n_B V} P_{AB}(\mathbf{k}_1 - \mathbf{k}_2) \\ + \frac{\delta_{AB}^K \delta_{CD}^K}{n_A n_C V} P_{AC}(0) + \frac{1}{n_A^3 V} \delta_{AB}^K \delta_{BC}^K \delta_{CD}^K \\ + \frac{1}{n_A^2 V} \left[ \delta_{AB}^K \delta_{BC}^K P_{AC}(k_2) + \delta_{AB}^K \delta_{BD}^K P_{AC}(k_2) \right. \\ \left. + \delta_{AC}^K \delta_{CD}^K P_{AB}(k_1) + \delta_{BC}^K \delta_{CD}^K P_{AB}(k_1) \right]. \quad (A.5)$$

In case of a single tracer ($A = B = C = D$), the above expression recovers the standard covariance given by Eq. (11).

By Fourier transforming Eq. (A.5), integrating over separation and redshift bins, and neglecting the high-order terms, we obtain the expression for the cross-2PCF covariance of different tracers (Eq. 16).

## Appendix B: Non-Poissonian shot noise

<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Halo power spectrum. Top panel: Measured (solid black lines), predicted (dashed darker lines), and fitted (dotted lighter lines) halo power spectrum for boxes at three different redshifts. Bottom panel: Percent residuals of the predictions with respect to the measured residuals.</figcaption></figure>

<table><caption>**Table B.1.** Best-fit parameters for the power spectrum parameters.</caption><thead><tr><th>Redshift</th><th>α</th><th>β</th></tr></thead><tbody><tr><td>0.5</td><td>0.012 ± 0.010</td><td>0.987 ± 0.002</td></tr><tr><td>1.0</td><td>0.114 ± 0.004</td><td>1.006 ± 0.002</td></tr><tr><td>1.5</td><td>0.104 ± 0.002</td><td>1.013 ± 0.002</td></tr><tr><td>Reference</td><td>0</td><td>1</td></tr></tbody></table>

To better assess the impact of shot-noise and its deviation from the Poissonian prediction, we studied the power spectrum, that is, the quantity that is directly affected by this correction. To avoid complications due to the redshift integration and the geometry of the survey, we measured the power spectrum from 1000 cubical boxes of size $L = 3870 h^{-1}$ Mpc with the same properties and cosmology of the light cones described in Sect. 3.1. We considered three redshifts $z = 0.5, 1.0, \text{and } 1.5$.

We computed the analytical total halo power spectrum as

$$P_{\text{h,tot}}(k) = \bar{b}^2 P_m(k) + \frac{1}{n}, \quad (B.1)$$
where the matter power spectrum wads calculated by means of the CAMB code (Lewis et al. 2000). We compared this quantity with the measured total power spectrum averaged over the 1000 boxes for the three redshift values.

In Fig. B.1 we show the comparison of the measured and predicted power spectra (to facilitate comparison, we show the halo power spectrum, i.e., the total spectrum minus the shot noise): while at low redshift the two quantities agree on almost all the scales, they clearly deviate between the observed and predicted spectra, which increase with redshift to more than 20%. We tried to correct these discrepancies by fitting the two parameters {$\alpha, \beta$}, which account for the correction to shot noise and halo bias, respectively, directly from the power spectrum; the best fits are shown in Table B.1 and the resulting power spectra are shown in Fig. B.1 (dotted lines). We obtain an agreement of the fitted power spectra within the 5% level at all the linear scales at all redshifts. Nevertheless, the values deviate by several $\sigma$ from the best-fit parameters found from the covariance fit in the corresponding redshift intervals (see Table 1). This confirms that the parameters in the covariance, in addition to correct for the incorrect prediction of bias and shot noise, also absorb the effect of the missing higher-order terms in the model.

## Appendix C: Covariance parameter fit

<figure><img src="image_14.png" /><figcaption>Fig. C.1. $\chi^2$ distribution for the numerical, analytical, and analytical with fitted parameters covariance matrices. The reference distribution is plotted in black.</figcaption></figure>

We show here the results of the covariance parameters fit, following the method proposed in Fumagalli et al. (2022). In Fig. C.1 we show the $\chi^2$ distributions computed with respect to the measurements from the 1000 light cones for the three covariance matrices. In Table C.1 we report the mean and standard deviation for each distribution, with the corresponding $1\sigma$ uncertainties computed with the bootstrap technique. By construction, we expect the numerical matrix to perfectly follow the reference distribution. While this is true for the mean value, the variance differs by ~4$\sigma$ from the expected value; this discrepancy is ascribed to the noise in the numerical matrix that tightens the distribution. Because the errors are quite small, this distribution

<table><caption>Table C.1. $\chi^2$ distribution test values.</caption><thead><tr><th></th><th>mean</th><th>variance</th></tr></thead><tbody><tr><td>Numerical</td><td>149.0 ± 0.5</td><td>255.0 ± 11.7</td></tr><tr><td>Model</td><td>195.0 ± 0.7</td><td>494.0 ± 22.8</td></tr><tr><td>Model + fit</td><td>152.0 ± 0.6</td><td>302.0 ± 13.9</td></tr><tr><td>Reference</td><td>150</td><td>300</td></tr></tbody></table>

can be considered to agree well with the expected one. Instead, the model of Eq. (12) produces a distribution that is several $\sigma$ different from the expectation, confirming that this model is not suitable for describing the covariance of the data. Finally, the fitted matrix agrees well with the reference distribution for the mean value and the variance. This proves the goodness of our fit and ensures that the resulting model is able to correctly describe the covariance as well as if not better than the numerical matrix. The fitting process provides consistent results when the fit is performed with only 100 simulations.

## Appendix D: Cosmology-dependent matrix

To further explore the comparison between the fully cosmology-dependent likelihood analysis and the iterative method, we generated 100 synthetic light cones starting from a Gaussian distribution, with the amplitude given by the covariance model at the input cosmology. In this way, we ensured that the Gaussian distribution was the true likelihood describing the data and not just an approximation. We repeated the analysis of the 100 light cones described in Sect. 4.4, finding a mean value $\langle\Delta\text{DIC}\rangle_{\text{synth}} = -8.4 \pm 1.7$, to be compared with the value from the analysis of the 100 PINOCCHIO mocks, that is, $\langle\Delta\text{DIC}\rangle_{\text{sims}} = -6.3 \pm 0.9$. Moreover, by comparing the FoM of the cosmology-dependent covariance and the iterative method, we obtain a mean variation of $\langle\Delta\text{FoM}\rangle_{\text{synth}} = 176 \pm 38\%$ for the synthetic catalogs and $\langle\Delta\text{FoM}\rangle_{\text{sims}} = 142 \pm 33\%$ for the PINOCCHIO mocks. Both the DIC and the FoM analyses indicate that the two analysis are fully consistent. Although this result still does not define which posteriors are correct when the true likelihood is unknown, it shows that at least for this particular analysis, the narrowing of the posteriors does not primarily depend on some incorrect approximation of the likelihood function. In other words, when a Gaussian likelihood is assumed, it is possible to extract information from the cosmology dependence of the covariance.

## Appendix E: Cosmology dependence of the number count covariance

We discuss here some additional considerations about the cosmology dependence of the covariance, using as an example the cluster number counts. It has been shown in Euclid Collaboration (2021) that the use of a fixed covariance in the likelihood analysis can bring an under/overestimation of the FoM by more than 40%, if the cosmology at which the covariance is evaluated deviates from the fiducial values of an amount of $2\sigma$ from Planck Collaboration VI (2020). Instead, the input and varying covariance cases do not vary significantly. We performed the same test on the cosmology dependence of the covariance described in Sect. 4.4, obtaining the posteriors shown in Fig. E.1. Unlike for the clustering case, for the number counts, the mean value is much more constraining than its covariance,
where the matter power spectrum wads calculated by means of the CAMB code (Lewis et al. 2000). We compared this quantity with the measured total power spectrum averaged over the 1000 boxes for the three redshift values.

In Fig. B.1 we show the comparison of the measured and predicted power spectra (to facilitate comparison, we show the halo power spectrum, i.e., the total spectrum minus the shot noise): while at low redshift the two quantities agree on almost all the scales, they clearly deviate between the observed and predicted spectra, which increase with redshift to more than 20%. We tried to correct these discrepancies by fitting the two parameters {$α, β$}, which account for the correction to shot noise and halo bias, respectively, directly from the power spectrum; the best fits are shown in Table B.1 and the resulting power spectra are shown in Fig. B.1 (dotted lines). We obtain an agreement of the fitted power spectra within the 5% level at all the linear scales at all redshifts. Nevertheless, the values deviate by several $σ$ from the best-fit parameters found from the covariance fit in the corresponding redshift intervals (see Table 1). This confirms that the parameters in the covariance, in addition to correct for the incorrect prediction of bias and shot noise, also absorb the effect of the missing higher-order terms in the model.

## Appendix C: Covariance parameter fit

<figure><img src="image_14.png" /><figcaption>**Fig. C.1.** $χ^2$ distribution for the numerical, analytical, and analytical with fitted parameters covariance matrices. The reference distribution is plotted in black.</figcaption></figure>

We show here the results of the covariance parameters fit, following the method proposed in Fumagalli et al. (2022). In Fig. C.1 we show the $χ^2$ distributions computed with respect to the measurements from the 1000 light cones for the three covariance matrices. In Table C.1 we report the mean and standard deviation for each distribution, with the corresponding $1σ$ uncertainties computed with the bootstrap technique. By construction, we expect the numerical matrix to perfectly follow the reference distribution. While this is true for the mean value, the variance differs by ~4$σ$ from the expected value; this discrepancy is ascribed to the noise in the numerical matrix that tightens the distribution. Because the errors are quite small, this distribution

<table><caption>**Table C.1.** $χ^2$ distribution test values.</caption><thead><tr><th></th><th>mean</th><th>variance</th></tr></thead><tbody><tr><td>Numerical</td><td>149.0 ± 0.5</td><td>255.0 ± 11.7</td></tr><tr><td>Model</td><td>195.0 ± 0.7</td><td>494.0 ± 22.8</td></tr><tr><td>Model + fit</td><td>152.0 ± 0.6</td><td>302.0 ± 13.9</td></tr><tr><td>Reference</td><td>150</td><td>300</td></tr></tbody></table>

can be considered to agree well with the expected one. Instead, the model of Eq. (12) produces a distribution that is several $σ$ different from the expectation, confirming that this model is not suitable for describing the covariance of the data. Finally, the fitted matrix agrees well with the reference distribution for the mean value and the variance. This proves the goodness of our fit and ensures that the resulting model is able to correctly describe the covariance as well as if not better than the numerical matrix. The fitting process provides consistent results when the fit is performed with only 100 simulations.

## Appendix D: Cosmology-dependent matrix

To further explore the comparison between the fully cosmology-dependent likelihood analysis and the iterative method, we generated 100 synthetic light cones starting from a Gaussian distribution, with the amplitude given by the covariance model at the input cosmology. In this way, we ensured that the Gaussian distribution was the true likelihood describing the data and not just an approximation. We repeated the analysis of the 100 light cones described in Sect. 4.4, finding a mean value $\langle\Delta\text{DIC}\rangle_{\text{synth}} = -8.4 \pm 1.7$, to be compared with the value from the analysis of the 100 PINOCCHIO mocks, that is, $\langle\Delta\text{DIC}\rangle_{\text{sims}} = -6.3 \pm 0.9$. Moreover, by comparing the FoM of the cosmology-dependent covariance and the iterative method, we obtain a mean variation of $\langle\Delta\text{FoM}\rangle_{\text{synth}} = 176 \pm 38\%$ for the synthetic catalogs and $\langle\Delta\text{FoM}\rangle_{\text{sims}} = 142 \pm 33\%$ for the PINOCCHIO mocks. Both the DIC and the FoM analyses indicate that the two analysis are fully consistent. Although this result still does not define which posteriors are correct when the true likelihood is unknown, it shows that at least for this particular analysis, the narrowing of the posteriors does not primarily depend on some incorrect approximation of the likelihood function. In other words, when a Gaussian likelihood is assumed, it is possible to extract information from the cosmology dependence of the covariance.

## Appendix E: Cosmology dependence of the number count covariance

We discuss here some additional considerations about the cosmology dependence of the covariance, using as an example the cluster number counts. It has been shown in Euclid Collaboration (2021) that the use of a fixed covariance in the likelihood analysis can bring an under/overestimation of the FoM by more than 40%, if the cosmology at which the covariance is evaluated deviates from the fiducial values of an amount of $2σ$ from Planck Collaboration VI (2020). Instead, the input and varying covariance cases do not vary significantly. We performed the same test on the cosmology dependence of the covariance described in Sect. 4.4, obtaining the posteriors shown in Fig. E.1. Unlike for the clustering case, for the number counts, the mean value is much more constraining than its covariance,
<figure><img src="image_15.png" /><figcaption>**Fig. E.1.** Same as Fig. 9 for cluster number counts.</figcaption></figure>

making the degeneracy direction of the latter totally irrelevant. In
this case, the covariance thus only contributes as an estimation of
the uncertainty, without adding further independent information.

**Appendix F: Mass dependence of the covariance**

In Fig. F.1 we show the different diagonal contribution to the covariance, composing the model of Eq. (12). The components of the model show a different dependence on the mass threshold, which drastically changes the impact of the various contributions. While at the high-mass threshold, the number density is low, making the shot-noise terms dominant, at low-mass cuts, this contribution is much lower than that from the halo power spectrum.

This mass dependence has an important consequence when considering the cosmology dependence of the covariance: Because the additional information extracted from the covariance is mainly given by the shot-noise term, we expect the impact of the cosmology-dependent covariance to be negligible when considering low-mass objects. In other words, the results shown in this work only hold for cluster clustering, while they would not have any substantial impact on galaxy clustering. Moreover, the results shown here consider the ideal case of a simple mass cut; this may vary in real surveys, where selection functions set a mass- and redshift-dependent richness threshold. The study of the impact of the covariance in real observable conditions is left for a future work.

<figure><img src="image_16.png" /><figcaption>**Fig. F.1.** Diagonal contributions to the covariance matrix composing the model of Eq. (12). The red, green, and blue lines represent the three terms of the Gaussian covariance, and the yellow line is the non-Gaussian term. The comparison is shown for two different mass thresholds: a low-mass threshold $M_{\text{th}} = 1 \times 10^{13} M_{\odot}$ in the top panel, and a high-mass threshold $M_{\text{th}} = 1 \times 10^{14} M_{\odot}$ in the bottom panel. The terms are obtained by evaluating the covariance model at the input cosmology in the redshift bin $z = 0.8 - 1.2$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. E.1.** Same as Fig. 9 for cluster number counts.</figcaption></figure>

making the degeneracy direction of the latter totally irrelevant. In
this case, the covariance thus only contributes as an estimation of
the uncertainty, without adding further independent information.

**Appendix F: Mass dependence of the covariance**

In Fig. F.1 we show the different diagonal contribution to the covariance, composing the model of Eq. (12). The components of the model show a different dependence on the mass threshold, which drastically changes the impact of the various contributions. While at the high-mass threshold, the number density is low, making the shot-noise terms dominant, at low-mass cuts, this contribution is much lower than that from the halo power spectrum.

This mass dependence has an important consequence when considering the cosmology dependence of the covariance: Because the additional information extracted from the covariance is mainly given by the shot-noise term, we expect the impact of the cosmology-dependent covariance to be negligible when considering low-mass objects. In other words, the results shown in this work only hold for cluster clustering, while they would not have any substantial impact on galaxy clustering. Moreover, the results shown here consider the ideal case of a simple mass cut; this may vary in real surveys, where selection functions set a mass- and redshift-dependent richness threshold. The study of the impact of the covariance in real observable conditions is left for a future work.

<figure><img src="image_16.png" /><figcaption>**Fig. F.1.** Diagonal contributions to the covariance matrix composing the model of Eq. (12). The red, green, and blue lines represent the three terms of the Gaussian covariance, and the yellow line is the non-Gaussian term. The comparison is shown for two different mass thresholds: a low-mass threshold $M_{\text{th}} = 1 \times 10^{13} M_{\odot}$ in the top panel, and a high-mass threshold $M_{\text{th}} = 1 \times 10^{14} M_{\odot}$ in the bottom panel. The terms are obtained by evaluating the covariance model at the input cosmology in the redshift bin $z = 0.8 - 1.2$.</figcaption></figure>