<table><caption>**Table A2.** Functional forms to model Cov(ΔΣ, ln N<sub>gal</sub>). The radius in log-space x is transformed to x̄ ≡ (x − γ)/τ by a horizontal offset γ and a characteristic scale τ. The functions f(x̄) are normalized so that f(x̄) asymptotically goes to 1 at +∞, −1 at −∞, f(0) = 0 and f'(0) = 1. Finally, we wrap f(x̄) by the function p(f(x̄)) ≡ s(f(x̄ + g)) to include a vertical offset g and amplitude parameter s. Together θ ∈ {τ, γ, g, s} form the set of model parameters that allow us to make apple-to-apple comparisons between models.</caption><tr><td>Error function (fiducial)</td><td>s( erf(√π/2 x̄) + g )</td></tr><tr><td>Logistics function</td><td>s( 2/(1+e<sup>x̄</sup>) - 1 + g )</td></tr><tr><td>Inverse tangent</td><td>s( 2/π arctan (π/2 x̄) + g )</td></tr><tr><td>Algebraic second order</td><td>s(x̄/(1 + x̄<sup>2</sup>)<sup>1/2</sup> + g )</td></tr></table>

marginally outperforms the full error function in most cases, along with the posterior constraints of parameters shown in Fig. D1.

## APPENDIX B: PARTICLE RESOLUTION AND ITS IMPACT ON MEASUREMENT ERRORS

Using the 300 Cori Haswell node hours allocated by the NERSC to this project, we measured ΔΣ for ~5000 clusters using dark matter particles downsampled by a factor of 10, in 20 log-spaced radial bins, at a projection depth of 200 h<sup>-1</sup>Mpc.

At a downsampling rate of 10, our effective dark matter particle resolution is $M_p \approx 1.51 \times 10^{10} h^{-1} M_\odot$. The error in ΔΣ comes from three sources: (i) cosmic variance, (ii) Poisson noise, and (iii) the intrinsic diversity of haloes accounted for by secondary halo properties. In Section 6.2, we presented the contribution to ΔΣ

<figure><img src="image_16.png" /><figcaption>**Figure B1.** The standard error of ΔΣ measurements tested on a benchmark bin of $M_{200c} \in [1 \times 10^{14}, 2 \times 10^{14})$ M<sub>☉</sub> h<sup>-1</sup> at z = 0. The standard error is estimated using the bootstrap method for the N = 500 clusters with dark matter particles downsampled by a factor of 200, 100, and 10 (solid lines). At our current resolution (nth = 10, solid green line), the standard error is just above the cosmic variance at small scales and drops below the cosmic variance at large scales. The solid black line is density fluctuation estimated from the cosmic variance floor, as described in equations (B2) and (B1). In the ideal case that Poisson noise accounts for all the standard error, fully sampling all particles (nth = 1, red dotted line) will reduce the standard error by a factor of $\sqrt{10}$, rendering it just below the cosmic variance floor at small scales. In the realistic case that the standard error for ΔΣ contributes from both Poisson noise and the intrinsic diversity of halo profiles, the fully sampled standard error should be on par with the cosmic variance at small scales.</figcaption></figure>

scatter from secondary halo properties. Here, we compare the Poisson noise to the cosmic variance floor.

The cosmic variance introduces fluctuations in a 2D surface density fluctuation, given by

$$ \delta\Delta\Sigma(R) = D_p \rho_m \sigma(R), \quad (B1) $$

where $D_p = 200h^{-1}Mpc$ is the projection depth, $\rho_m(z)$ is the mean density of the universe at that redshift, and $\sigma(R)$ is the root mean squared matter density fluctuation, given by

$$ \sigma^2(R) = \int \Delta^2(k) \left( \frac{3j_1(kR)}{kR} \right)^2 d \ln k, \quad (B2) $$

which is smoothed over an area of A = $4\pi R^2$, $\Delta^2(k)$ is the matter power spectrum for a wavenumber k, and $j_1$ is the Bessel function of the first order.

Fig. B1 shows the standard error of ΔΣ at a benchmark bin $M_{200c} \in [1 \times 10^{14}, 2 \times 10^{14})$ M<sub>☉</sub> h<sup>-1</sup> at z = 0.00. At particle downsampling factors of 200, 100, and 10, the reduction in error is consistent with the Poisson term of $\sqrt{N}$, indicating that at these sampling rates, Poisson noise dominates. At our current downsampling rate of nth = 10, the standard error is just above the cosmic variance floor at small scales and drops below the cosmic variance floor at large scales. In the ideal case that Poisson noise accounts for all the standard error, fully sampling all particles (nth = 1, red dotted line) will reduce the standard error by a factor of $\sqrt{10}$, rendering it just below the cosmic variance floor at small scales. In the realistic case that the standard error for ΔΣ contributes from both Poisson noise and the intrinsic diversity of halo profiles, the fully sampled standard error should be on par with the cosmic variance at small scales. A future study with fully sampled particles should yield greater statistical constraints.

## APPENDIX C: DERIVATION OF SECOND ORDER EXPANSION AROUND THE HMF

Following the formalism from Evrard et al. (2014), we derive equation (18). The mean observable-mass scaling relation is given by the expression

$$ \langle s_i | \ln M_0 \rangle = \alpha_i \ln M_0 + \pi_i, \quad (C1) $$

for $s_i \in \{\Delta\Sigma, \ln N_{gal}\}$ and pivot mass $M_0$. We now denote the deviation from the mean relation as $\delta_i$, which from rearranging the terms in equation (C1) is given by $\delta_i = \frac{(s_i - \pi_i)}{\alpha_i} - \ln M_0$. From Evrard et al. (2014) the expression for the observable scaling relation for generic observables $\{a, b\}$ that follow a log-linear scaling relation as in equation (C1) is given by

$$ \langle \delta_b | s_a \rangle = x_a \left[ \langle \ln M | s_a \rangle + (\gamma_1 + \gamma_2 \delta_a) r_{ab} \sigma_{\ln M|a,1} \sigma_{\ln M|b,1} \right], \quad (C2) $$

where $\gamma_1$ and $\gamma_2$ are the first- and second-order coefficients of the Taylor expansion of the HMF around the pivot mass $M_0$ and $x_a = (1 + \gamma_2 \sigma_{\ln M|a,1})^{-1}$ the curvature term. The subscript 1 denotes the scatter for the HMF expanded to first order.

We now convert the left-hand side of equation (C2) from the deviation from the mean scaling relation, $\delta_b$, to the observable $s_b$ to arrive at the expression

$$ \begin{aligned} \langle s_b | s_a \rangle &= [\alpha_b x_a \langle \ln M | s_a \rangle + \ln M_0 - \pi_b] \\ & \quad + [\alpha_b x_a (\gamma_1 + \gamma_2 \delta_a) r_{ab}] \sigma_{\ln M|a,1} \sigma_{\ln M|b,1} \\ &= \langle b | a, z \rangle_{\text{fid}} \\ & \quad + \text{Cov}(a, b) \times \left[ \frac{x_a}{\alpha_a^2} (\alpha_a \gamma_1 + \gamma_2(s_a - \pi_a)) \right], \end{aligned} \quad (C3) $$

where we made use of the fact Cov(a, b|M, z) = $r_{ab}\sigma_{a|M}\sigma_{b|M}$ and that the mass scatter conditioned on the observable to first order is
<table><caption>**Table A2.** Functional forms to model Cov(ΔΣ, ln N<sub>gal</sub>). The radius in log-space x is transformed to x&#773; = (x - γ)/τ by a horizontal offset γ and a characteristic scale τ. The functions f(x&#773;) are normalized so that f(x&#773;) asymptotically goes to 1 at +∞, -1 at -∞, f(0) = 0 and f'(0) = 1. Finally, we wrap f(x&#773;) by the function p(f(x&#773;)) = s(f(x&#773; + g)) to include a vertical offset g and amplitude parameter s. Together θ ∈ {τ, γ, g, s} form the set of model parameters that allow us to make apple-to-apple comparisons between models.</caption><tr><td>Error function (fiducial)</td><td>s( erf(√π/2 x&#773;) + g )</td></tr><tr><td>Logistics function</td><td>s( 2/(1+e<sup>-x&#773;)) - 1 + g )</td></tr><tr><td>Inverse tangent</td><td>s( 2/π arctan (π/2 x&#773;) + g )</td></tr><tr><td>Algebraic second order</td><td>s(x&#773;/((1 + x&#773;)<sup>1/2</sup> + g))</td></tr></table>

marginally outperforms the full error function in most cases, along with the posterior constraints of parameters shown in Fig. D1.

## APPENDIX B: PARTICLE RESOLUTION AND ITS IMPACT ON MEASUREMENT ERRORS

Using the 300 Cori Haswell node hours allocated by the NERSC to this project, we measured ΔΣ for ~5000 clusters using dark matter particles downsampled by a factor of 10, in 20 log-spaced radial bins, at a projection depth of 200 h<sup>-1</sup>Mpc.

At a downsampling rate of 10, our effective dark matter particle resolution is M<sub>p</sub> ≈ 1.51 × 10<sup>10</sup>h<sup>-1</sup>M<sub>☉</sub>. The error in ΔΣ comes from three sources: (i) cosmic variance, (ii) Poisson noise, and (iii) the intrinsic diversity of haloes accounted for by secondary halo properties. In Section 6.2, we presented the contribution to ΔΣ

<figure><img src="image_14.png" /><figcaption>**Figure B1.** The standard error of ΔΣ measurements tested on a benchmark bin of M<sub>200c</sub> ∈ [1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>) M<sub>☉</sub> h<sup>-1</sup> at z = 0. The standard error is estimated using the bootstrap method for the N = 500 clusters with dark matter particles downsampled by a factor of 200, 100, and 10 (solid lines). At our current resolution (nth = 10, solid green line), the standard error is just above the cosmic variance at small scales and drops below the cosmic variance at large scales. The solid black line is density fluctuation estimated from the cosmic variance floor, as described in equations (B2) and (B1). In the ideal case that Poisson noise accounts for all the standard error, fully sampling all particles (nth = 1, red dotted line) will reduce the standard error by a factor of √10, rendering it just below the cosmic variance floor at small scales. In the realistic case that the standard error for ΔΣ contributes from both Poisson noise and the intrinsic diversity of halo profiles, the fully sampled standard error should be on par with the cosmic variance at small scales.</figcaption></figure>

scatter from secondary halo properties. Here, we compare the Poisson noise to the cosmic variance floor.

The cosmic variance introduces fluctuations in a 2D surface density fluctuation, given by

$$ \delta\Delta\Sigma(R) = D_p \rho_m \sigma(R), \quad (B1) $$

where D<sub>p</sub> = 200h<sup>-1</sup>Mpc is the projection depth, ρ<sub>m</sub>(z) is the mean density of the universe at that redshift, and σ(R) is the root mean squared matter density fluctuation, given by

$$ \sigma^2(R) = \int \Delta^2(k) \left( \frac{3j_1(kR)}{kR} \right)^2 d \ln k, \quad (B2) $$

which is smoothed over an area of A = 4πR<sup>2</sup>, Δ<sup>2</sup>(k) is the matter power spectrum for a wavenumber k, and j<sub>1</sub> is the Bessel function of the first order.

Fig. B1 shows the standard error of ΔΣ at a benchmark bin M<sub>200c</sub> ∈ [1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>) M<sub>☉</sub> h<sup>-1</sup> at z = 0.00. At particle downsampling factors of 200, 100, and 10, the reduction in error is consistent with the Poisson term of √N, indicating that at these sampling rates, Poisson noise dominates. At our current downsampling rate of nth = 10, the standard error is just above the cosmic variance floor at small scales and drops below the cosmic variance floor at large scales. In the ideal case that Poisson noise accounts for all the standard error, fully sampling all particles (nth = 1, red dotted line) will reduce the standard error by a factor of √10, rendering it just below the cosmic variance floor at small scales. In the realistic case that the standard error for ΔΣ contributes from both Poisson noise and the intrinsic diversity of halo profiles, the fully sampled standard error should be on par with the cosmic variance at small scales. A future study with fully sampled particles should yield greater statistical constraints.

## APPENDIX C: DERIVATION OF SECOND ORDER EXPANSION AROUND THE HMF

Following the formalism from Evrard et al. (2014), we derive equation (18). The mean observable-mass scaling relation is given by the expression

$$ \langle s_i | \ln M_0 \rangle = \alpha_i \ln M_0 + \pi_i, \quad (C1) $$

for s<sub>i</sub> ∈ {ΔΣ, ln N<sub>gal</sub>} and pivot mass M<sub>0</sub>. We now denote the deviation from the mean relation as δ<sub>i</sub>, which from rearranging the terms in equation (C1) is given by δ<sub>i</sub> = (s<sub>i</sub> - π<sub>i</sub>) / α<sub>i</sub> - ln M<sub>0</sub>. From Evrard et al. (2014) the expression for the observable scaling relation for generic observables {a, b} that follow a log-linear scaling relation as in equation (C1) is given by

$$ \langle \delta_b | s_a \rangle = x_a \left[ \langle \ln M | s_a \rangle + (\gamma_1 + \gamma_2 \delta_a) r_{ab} \sigma_{\ln M|a,1} \sigma_{\ln M|b,1} \right], \quad (C2) $$

where γ<sub>1</sub> and γ<sub>2</sub> are the first- and second-order coefficients of the Taylor expansion of the HMF around the pivot mass M<sub>0</sub> and x<sub>a</sub> = (1 + γ<sub>2</sub>σ<sup>2</sup><sub>ln M|a,1</sub>)<sup>-1</sup> the curvature term. The subscript 1 denotes the scatter for the HMF expanded to first order.

We now convert the left-hand side of equation (C2) from the deviation from the mean scaling relation, δ<sub>b</sub>, to the observable s<sub>b</sub> to arrive at the expression

$$ \begin{aligned} \langle s_b | s_a \rangle &= [\alpha_b x_a \langle \ln M | s_a \rangle + \ln M_0 - \pi_b] \\ & \quad + [\alpha_b x_a (\gamma_1 + \gamma_2 \delta_a) r_{ab}] \sigma_{\ln M|a,1} \sigma_{\ln M|b,1} \\ &= \langle b | a, z \rangle_{\text{fid}} \\ & \quad + \text{Cov}(a, b) \times \left[ \frac{x_a}{\alpha_a^2} (\alpha_a \gamma_1 + \gamma_2(s_a - \pi_a)) \right], \end{aligned} \quad (C3) $$

where we made use of the fact Cov(a, b|M, z) = r<sub>ab</sub>σ<sub>a|M</sub>σ<sub>b|M</sub> and that the mass scatter conditioned on the observable to first order is
related to the observable scatter by $\sigma_{\ln M|a, 1} = \sigma_{a|M}/\alpha_a$, as shown in equation 4 in Evrard et al. (2014) for the multivariate case. Substituting $\ln N_{\text{gal}}$ for $a$ and $\Delta\Sigma$ for $b$ yields the expression for equation (18).

## APPENDIX D: ROBUSTNESS TESTING OF COVARIANCE MODELLING

The shape posterior is sampled by a Monte Carlo Markov Chain (MCMC) using the EMCEE package (Foreman-Mackey et al. 2013). We test for convergence by ensuring that the number of steps exceeds $t_{\text{auto}}/100$ for all parameters, where $t_{\text{auto}}$ is the integrated autocorrelation time as defined by Goodman & Weare (2010) and by ensuring that the convergence diagnostic denoted with $R$ (Gelman & Rubin 1992) across all walkers satisfy $R < 1.05$.

The posterior distribution according to Bayes theorem is given as

$$p(\theta|\{y_i\}) \propto p(\{y_i\}|\theta)p(\theta) \quad (\text{D1})$$

$$= \prod_i p(y_i|\theta)p(\theta), \quad (\text{D2})$$

where the second line assumes independent and identical distribution (i.i.d) for the data vectors. We set uniform priors $p(\theta)$ shown in Table A1 with signs and ranges motivated by the shape of the covariance (i.e. a negative $\gamma$ and positive $\tau$ to offset $f(\tilde{x})$ to the left and a positive $s$ and negative $g$ shifts the fitted curve downwards).

We measure the goodness of fit using the left-tail $p$-value for the $\chi^2$ with $N_{\text{data}} - N_{\text{dim}} = 20 - 4 = 16$ degrees of freedom. We compare between models by the Deviance Information Criterion defined as

$$\text{DIC} = 2\bar{D}(\bar{\theta}) - D(\bar{\theta}), \quad (\text{D3})$$

where $\bar{\theta}$ is the best-fitting parameters, and $D(\bar{\theta})$ is defined as

$$D(\bar{\theta}) = -2 \log (P(\{x_i\}|\theta)). \quad (\text{D4})$$

The performance between different functional forms (Table A2) is reported in Table A1.

The summary statistics for the posterior distribution of the covariance models are listed in Table D1 as plotted against measurements in Fig. 4. Among the functions, the error function has either a better or comparable fit to all other functions in all other bins, as indicated by their DIC parameters. In two bins $M_{200c} \in [5 \times 10^{14}, 1 \times 10^{15})$ at $z = 0.49$ and $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{15})$ at $z = 1.03$, the amplitude of the covariance is too small relative to their errors for shape parameters to be well-constrained. The right tail $p$-value for $\chi^2$ is $p > 0.05$ for all but one bin. For this reason, we take the full error function as the nominal functional form.

For $R \ge R_{\text{vir}}$ or $R \ge R_{200c}$, we find the covariance to be null at $p$-values >0.01. A zero covariance at large scales implies $g = -1$ which coincides with the reduced model. We compare the results of the error function of the reduced model to the full model and find their performance varies from bin to bin as indicated by the DIC (Table D2). The posteriors of the reduced model provide marginally tighter constraints than the full model (Fig. D1).

<table><caption>**Table D1.** Summary statistics for $\text{Cov}(\Delta\Sigma, N_{\text{gal}}|M, z)$ binned by $R_{200c}$ and $M_{200c}$, with $N_{\text{gal}}$ defined inside the halo $R_{200c}$. Columns 2–5 are the best-fitting parameters for the nominal error function and their $1\sigma$ ranges. Columns 6–8 are the difference between the DIC of the logistics, algebraic, and inverse tangent models with the nominal error function, respectively. Column 9 is the right-tail $p$-value as measured by the $\chi^2$ statistic with $20-4 = 16$ degrees of freedom. Across all bins with applicable posterior constraints, the error function out-performs or is comparable to alternative functional forms as indicated by the difference in DIC, and has $p \ge 0.01$ in all but one bin. In two bins $M_{200c} \in [5 \times 10^{14}, 1 \times 10^{15})$ at $z = 0.49$ and $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{15})$ at $z = 1.03$ the size of the covariance is too small relative to the size of their errors for shape parameters to be constrained. The covariance in these two bins is consistent with null at $p = 0.01$ and $p = 0.05$ levels, respectively.</caption><thead><tr><th rowspan="2">Mass and redshift</th><th rowspan="2">τ</th><th rowspan="2">γ</th><th rowspan="2">g</th><th colspan="3">Error function (full): s( erf(π/2 x̄) + g )</th><th rowspan="2">ΔDIC<sub>alg</sub></th><th rowspan="2">ΔDIC<sub>arctan</sub></th><th rowspan="2">p-value</th></tr><tr><th>10<sup>12</sup> x</th><th>ΔDIC<sub>log</sub></th><th>ΔDIC<sub>alg</sub></th></tr></thead><tbody><tr><td>[5 x 10<sup>13</sup>, 1 x 10<sup>14</sup>), z = 0.00</td><td>0.47<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.66<sup>+0.08</sup><sub>-0.11</sub></td><td>-0.99<sup>+0.005</sup><sub>+0.004</sub></td><td>2.88<sup>+0.71</sup><sub>-0.47</sub></td><td>4.8</td><td>17.8</td><td>31.5</td><td>0.27</td></tr><tr><td>[1 x 10<sup>14</sup>, 2 x 10<sup>14</sup>), z = 0.00</td><td>0.57<sup>+0.11</sup><sub>-0.08</sub></td><td>-0.76<sup>+0.14</sup><sub>-0.22</sub></td><td>-1.008<sup>+0.0005</sup><sub>-0.0007</sub></td><td>2.73<sup>+0.13</sup><sub>-0.65</sub></td><td>-2.9</td><td>-1.7</td><td>2.8</td><td>0.88</td></tr><tr><td>[2 x 10<sup>14</sup>, 5 x 10<sup>14</sup>), z = 0.00</td><td>0.41<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.51<sup>+0.05</sup><sub>-0.06</sub></td><td>-0.997<sup>+0.004</sup><sub>-0.0004</sub></td><td>1.93<sup>+0.24</sup><sub>-0.18</sub></td><td>4.4</td><td>17.4</td><td>35.0</td><td>0.79</td></tr><tr><td>[5 x 10<sup>14</sup>, 1 x 10<sup>15</sup>), z = 0.00</td><td>0.27<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.40<sup>+0.04</sup><sub>-0.05</sub></td><td>-1.006<sup>+0.005</sup><sub>-0.0005</sub></td><td>1.03<sup>+0.13</sup><sub>-0.11</sub></td><td>-0.27</td><td>14.3</td><td>29.1</td><td>0.0009</td></tr><tr><td>[5 x 10<sup>13</sup>, 1 x 10<sup>14</sup>), z = 0.49</td><td>0.28<sup>+0.6</sup><sub>-0.05</sub></td><td>-0.36<sup>+0.05</sup><sub>-0.07</sub></td><td>-0.990<sup>+0.013</sup><sub>-0.013</sub></td><td>0.75<sup>+0.11</sup><sub>-0.09</sub></td><td>2.7</td><td>9.0</td><td>15.1</td><td>0.72</td></tr><tr><td>[1 x 10<sup>14</sup>, 2 x 10<sup>14</sup>), z = 0.49</td><td>0.35<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.43<sup>+0.07</sup><sub>-0.09</sub></td><td>-1.022<sup>+0.016</sup><sub>-0.017</sub></td><td>0.81<sup>+0.15</sup><sub>-0.10</sub></td><td>1.7</td><td>5.4</td><td>9.3</td><td>0.97</td></tr><tr><td>[2 x 10<sup>14</sup>, 5 x 10<sup>14</sup>), z = 0.49</td><td>0.24<sup>+0.08</sup><sub>-0.07</sub></td><td>-0.49<sup>+0.07</sup><sub>-0.10</sub></td><td>-0.995<sup>+0.012</sup><sub>-0.011</sub></td><td>0.51<sup>+0.14</sup><sub>-0.01</sub></td><td>0.6</td><td>2.4</td><td>4.0</td><td>0.63</td></tr><tr><td>[5 x 10<sup>14</sup>, 1 x 10<sup>15</sup>), z = 0.49</td><td></td><td></td><td></td><td>NA</td><td></td><td></td><td></td><td></td></tr><tr><td>[5 x 10<sup>13</sup>, 1 x 10<sup>14</sup>), z = 1.03</td><td>0.17<sup>+0.05</sup><sub>-0.04</sub></td><td>-0.31<sup>+0.05</sup><sub>-0.05</sub></td><td>-1.019<sup>+0.020</sup><sub>-0.022</sub></td><td>0.35<sup>+0.06</sup><sub>-0.05</sub></td><td>4.3</td><td>6.5</td><td>8.3</td><td>0.05</td></tr><tr><td>[1 x 10<sup>14</sup>, 2 x 10<sup>14</sup>), z = 1.03</td><td>0.21<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.42<sup>+0.05</sup><sub>-0.06</sub></td><td>-1.024<sup>+0.013</sup><sub>-0.0014</sub></td><td>0.51<sup>+0.09</sup><sub>-0.07</sub></td><td>-0.3</td><td>3.3</td><td>6.1</td><td>0.01</td></tr><tr><td>[2 x 10<sup>14</sup>, 5 x 10<sup>14</sup>), z = 1.03</td><td></td><td></td><td></td><td>NA</td><td></td><td></td><td></td><td></td></tr></tbody></table>
related to the observable scatter by $\sigma_{\ln M|a, 1} = \sigma_{a|M}/\alpha_a$, as shown in equation 4 in Evrard et al. (2014) for the multivariate case. Substituting $\ln N_{\text{gal}}$ for $a$ and $\Delta\Sigma$ for $b$ yields the expression for equation (18).

## APPENDIX D: ROBUSTNESS TESTING OF COVARIANCE MODELLING

The shape posterior is sampled by a Monte Carlo Markov Chain (MCMC) using the EMCEE package (Foreman-Mackey et al. 2013). We test for convergence by ensuring that the number of steps exceeds $t_{\text{auto}}/100$ for all parameters, where $t_{\text{auto}}$ is the integrated autocorrelation time as defined by Goodman & Weare (2010) and by ensuring that the convergence diagnostic denoted with $R$ (Gelman & Rubin 1992) across all walkers satisfy $R < 1.05$.

The posterior distribution according to Bayes theorem is given as

$$p(\theta|\{y_i\}) \propto p(\{y_i\}|\theta)p(\theta) \quad (\text{D1})$$

$$= \prod_i p(y_i|\theta)p(\theta), \quad (\text{D2})$$

where the second line assumes independent and identical distribution (i.i.d) for the data vectors. We set uniform priors $p(\theta)$ shown in Table A1 with signs and ranges motivated by the shape of the covariance (i.e. a negative $\gamma$ and positive $\tau$ to offset $f(\tilde{x})$ to the left and a positive $s$ and negative $g$ shifts the fitted curve downwards).

We measure the goodness of fit using the left-tail $p$-value for the $\chi^2$ with $N_{\text{data}} - N_{\text{dim}} = 20 - 4 = 16$ degrees of freedom. We compare between models by the Deviance Information Criterion defined as

$$\text{DIC} = 2\bar{D}(\bar{\theta}) - D(\bar{\theta}), \quad (\text{D3})$$

where $\bar{\theta}$ is the best-fitting parameters, and $D(\bar{\theta})$ is defined as

$$D(\bar{\theta}) = -2 \log (P(\{x_i\}|\theta)). \quad (\text{D4})$$

The performance between different functional forms (Table A2) is reported in Table A1.

The summary statistics for the posterior distribution of the covariance models are listed in Table D1 as plotted against measurements in Fig. 4. Among the functions, the error function has either a better or comparable fit to all other functions in all other bins, as indicated by their DIC parameters. In two bins $M_{200c} \in [5 \times 10^{14}, 1 \times 10^{15})$ at $z = 0.49$ and $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{15})$ at $z = 1.03$, the amplitude of the covariance is too small relative to their errors for shape parameters to be well-constrained. The right tail $p$-value for $\chi^2$ is $p > 0.05$ for all but one bin. For this reason, we take the full error function as the nominal functional form.

For $R \ge R_{\text{vir}}$ or $R \ge R_{200c}$, we find the covariance to be null at $p$-values >0.01. A zero covariance at large scales implies $g = -1$ which coincides with the reduced model. We compare the results of the error function of the reduced model to the full model and find their performance varies from bin to bin as indicated by the DIC (Table D2). The posteriors of the reduced model provide marginally tighter constraints than the full model (Fig. D1).

<table><caption>**Table D1.** Summary statistics for Cov($\Delta\Sigma$, $N_{\text{gal}}|M$, $z$) binned by $R_{200c}$ and $M_{200c}$, with $N_{\text{gal}}$ defined inside the halo $R_{200c}$. Columns 2–5 are the best-fitting parameters for the nominal error function and their $1\sigma$ ranges. Columns 6–8 are the difference between the DIC of the logistics, algebraic, and inverse tangent models with the nominal error function, respectively. Column 9 is the right-tail $p$-value as measured by the $\chi^2$ statistic with $20-4 = 16$ degrees of freedom. Across all bins with applicable posterior constraints, the error function out-performs or is comparable to alternative functional forms as indicated by the difference in DIC, and has $p \ge 0.01$ in all but one bin. In two bins $M_{200c} \in [5 \times 10^{14}, 1 \times 10^{15})$ at $z = 0.49$ and $M_{200c} \in [2 \times 10^{14}, 5 \times 10^{15})$ at $z = 1.03$ the size of the covariance is too small relative to the size of their errors for shape parameters to be constrained. The covariance in these two bins is consistent with null at $p = 0.01$ and $p = 0.05$ levels, respectively.</caption><thead><tr><th rowspan="2">Mass and redshift</th><th rowspan="2">τ</th><th rowspan="2">γ</th><th rowspan="2">g</th><th colspan="3">Error function (full): s(erf(√π/2̄x) + g)</th><th rowspan="2">ΔDIC<sub>alg</sub></th><th rowspan="2">ΔDIC<sub>arctan</sub></th><th rowspan="2">p-value</th></tr><tr><th>10<sup>12</sup> × s</th><th>ΔDIC<sub>log</sub></th><th>ΔDIC<sub>alg</sub></th></tr></thead><tbody><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>), z = 0.00</td><td>0.47<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.66<sup>+0.08</sup><sub>-0.11</sub></td><td>-0.99<sup>+0.005</sup><sub>+0.004</sub></td><td>2.88<sup>+0.71</sup><sub>-0.47</sub></td><td>4.8</td><td>17.8</td><td>31.5</td><td>0.27</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>), z = 0.00</td><td>0.57<sup>+0.11</sup><sub>-0.08</sub></td><td>-0.76<sup>+0.14</sup><sub>-0.22</sub></td><td>-1.008<sup>+0.0005</sup><sub>-0.0007</sub></td><td>2.73<sup>+0.13</sup><sub>-0.65</sub></td><td>-2.9</td><td>-1.7</td><td>2.8</td><td>0.88</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>), z = 0.00</td><td>0.41<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.51<sup>+0.05</sup><sub>-0.06</sub></td><td>-0.997<sup>+0.004</sup><sub>-0.0004</sub></td><td>1.93<sup>+0.24</sup><sub>-0.18</sub></td><td>4.4</td><td>17.4</td><td>35.0</td><td>0.79</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>), z = 0.00</td><td>0.27<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.40<sup>+0.04</sup><sub>-0.05</sub></td><td>-1.006<sup>+0.005</sup><sub>-0.0005</sub></td><td>1.03<sup>+0.13</sup><sub>-0.11</sub></td><td>-0.27</td><td>14.3</td><td>29.1</td><td>0.0009</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>), z = 0.49</td><td>0.28<sup>+0.6</sup><sub>-0.05</sub></td><td>-0.36<sup>+0.05</sup><sub>-0.07</sub></td><td>-0.990<sup>+0.013</sup><sub>-0.013</sub></td><td>0.75<sup>+0.11</sup><sub>-0.09</sub></td><td>2.7</td><td>9.0</td><td>15.1</td><td>0.72</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>), z = 0.49</td><td>0.35<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.43<sup>+0.07</sup><sub>-0.09</sub></td><td>-1.022<sup>+0.016</sup><sub>-0.017</sub></td><td>0.81<sup>+0.15</sup><sub>-0.10</sub></td><td>1.7</td><td>5.4</td><td>9.3</td><td>0.97</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>), z = 0.49</td><td>0.24<sup>+0.08</sup><sub>-0.07</sub></td><td>-0.49<sup>+0.07</sup><sub>-0.10</sub></td><td>-0.995<sup>+0.012</sup><sub>-0.011</sub></td><td>0.51<sup>+0.14</sup><sub>-0.01</sub></td><td>0.6</td><td>2.4</td><td>4.0</td><td>0.63</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>), z = 0.49</td><td></td><td></td><td></td><td>NA</td><td></td><td></td><td></td><td></td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>), z = 1.03</td><td>0.17<sup>+0.05</sup><sub>-0.04</sub></td><td>-0.31<sup>+0.05</sup><sub>-0.05</sub></td><td>-1.019<sup>+0.020</sup><sub>-0.022</sub></td><td>0.35<sup>+0.06</sup><sub>-0.05</sub></td><td>4.3</td><td>6.5</td><td>8.3</td><td>0.05</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>), z = 1.03</td><td>0.21<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.42<sup>+0.05</sup><sub>-0.06</sub></td><td>-1.024<sup>+0.013</sup><sub>-0.0014</sub></td><td>0.51<sup>+0.09</sup><sub>-0.07</sub></td><td>-0.3</td><td>3.3</td><td>6.1</td><td>0.01</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>), z = 1.03</td><td></td><td></td><td></td><td>NA</td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>**Table D2.** Summary statistics for Cov(ΔΣ, N<sub>gal</sub>|M, z) binned by R<sub>200c</sub> and M<sub>200c</sub> with N<sub>gal</sub> defined inside the halo R<sub>200c</sub> for the reduced error function model. Compared with the full error function model, the performance of the reduced model varies from bin to bin – using ΔDIC > 3 as a statistically significant result, it outperforms the full model in 5/9 overlapping bins, under-performs in 3/9 bins, and is comparable in 2 bins. The reduced model is able to yield convergent chains for M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>) at z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>) at z = 1.03 but with poor constraints on the parameters.</caption><thead><tr><th rowspan="2">Mass & redshift</th><th colspan="5">Error function (reduced): s (erf(√π/2 x̄) - 1)</th></tr><tr><th>τ</th><th>γ</th><th>10<sup>12</sup> x̄</th><th>ΔDIC<sub>erf-full</sub></th><th>p-value</th></tr></thead><tbody><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>), z = 0.00</td><td>0.44<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.63<sup>+0.07</sup><sub>-0.09</sub></td><td>2.66<sup>+0.56</sup><sub>-0.39</sub></td><td>-5.2</td><td>0.06</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>), z = 0.00</td><td>0.69<sup>+0.11</sup><sub>-0.09</sub></td><td>-0.98<sup>+0.20</sup><sub>-0.27</sub></td><td>4.10<sup>+2.31</sup><sub>-1.19</sub></td><td>4.5</td><td>0.96</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>), z = 0.00</td><td>0.40<sup>+0.04</sup><sub>-0.03</sub></td><td>-0.50<sup>+0.05</sup><sub>-0.06</sub></td><td>1.90<sup>+0.26</sup><sub>-0.20</sub></td><td>0.2</td><td>0.78</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>), z = 0.00</td><td>0.29<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.40<sup>+0.04</sup><sub>-0.04</sub></td><td>1.12<sup>+0.13</sup><sub>-0.11</sub></td><td>6.1</td><td>0.005</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>), z = 0.49</td><td>0.25<sup>+0.04</sup><sub>-0.03</sub></td><td>-0.35<sup>+0.04</sup><sub>-0.04</sub></td><td>0.79<sup>+0.08</sup><sub>-0.07</sub></td><td>6.2</td><td>0.57</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>), z = 0.49</td><td>0.38<sup>+0.14</sup><sub>-0.09</sub></td><td>-0.48<sup>+0.12</sup><sub>-0.24</sub></td><td>0.85<sup>+0.46</sup><sub>-0.19</sub></td><td>-6.9</td><td>0.81</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>), z = 0.49</td><td>0.25<sup>+0.21</sup><sub>-0.12</sub></td><td>-0.55<sup>+0.14</sup><sub>-0.42</sub></td><td>0.40<sup>+0.69</sup><sub>-0.13</sub></td><td>NA</td><td>0.48</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>), z = 0.49</td><td>0.61<sup>+0.27</sup><sub>-0.24</sub></td><td>-0.90<sup>+0.49</sup><sub>-0.76</sub></td><td>1.01<sup>+3.05</sup><sub>-0.60</sub></td><td>-5.7</td><td>0.99</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>), z = 1.03</td><td>0.20<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.35<sup>+0.06</sup><sub>-0.08</sub></td><td>0.38<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.75</td><td>0.04</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>), z = 1.03</td><td>0.22<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.42<sup>+0.06</sup><sub>-0.08</sub></td><td>0.54<sup>+0.11</sup><sub>-0.09</sub></td><td>3.9</td><td>0.3</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>), z = 1.03</td><td>5.97<sup>+2.76</sup><sub>-3.01</sub></td><td>-0.55<sup>+3.30</sup><sub>-2.99</sub></td><td>0.02<sup>+0.02</sup><sub>-0.008</sub></td><td>NA</td><td>0.15</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Figure D1.** Posterior distribution of shape parameters in a benchmark bin of M<sub>200c</sub>/M<sub>vir</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>) at z = 0.00 under different binning schemes r<sub>p</sub> and N<sub>gal</sub> models with different halo boundaries R<sub>halo</sub>. The marginalized parameter constraints for the full model closely overlap one another, and using the reduced model with g = −1 marginally improves the posterior constraints. The plot was generated using pygtc (Bocquet & Carter 2016).</figcaption></figure>

## APPENDIX E: MODELLING SECONDARY PROPERTIES

We describe the linear regression model for richness used in Section 6.1. The same methodology is applied to ΔΣ in Section 6.2.

To model the expected natural logarithm of galaxy count (ln N<sub>gal</sub>), we decompose it linearly using secondary halo parameters listed in Table 3, as shown in equation (20). We employ the least squares method for linear regression and examine parameter redundancies. For over half of the bins, the parameters Γ<sub>inst</sub>, Γ<sub>100Myr</sub>, Γ<sub>2dyn</sub>, and Γ<sub>peak</sub> exhibit collinearity, with VIF exceeding 5. This outcome is expected,

as these quantities represent the same physical quantities smoothed over different time-scales. As for the reduced set of non-collinear parameters, their correlation coefficient are quantified in Shin & Diemer (2023) using the Erebos simulation suite. To determine which parameters to retain, we utilize the partial F-statistic.

Table E1 demonstrates the diminishing explanatory power of Π on the richness as seen by the diminishing R<sup>2</sup> and F<sub>partial</sub>. We consider the partial F-statistic serves as a heuristic measure for the explanatory power of a variable, defined as

$$F_{\text{partial}} = \frac{(\text{RSSE}_{\text{reduced}} - \text{RSSE}_{\text{full}})/p}{\text{RSSE}_{\text{full}}/(n-k)}, \quad \text{(E1)}$$

where RSSE is the residual sum of squared errors for the reduced model after removing the parameter in question and the full model containing Π ⊂ {a<sub>1/2</sub>, c<sub>vir</sub>, T/U, Γ<sub>2dyn</sub>, X<sub>off</sub>}, p is the number of parameters removed from the full model which in our case is by construction set to p = 1, n the number of data points, and k is the number of parameters in the full backward model. This statistic can be shown to be proportional to the contribution to the total R<sup>2</sup> uniquely explained by this parameter alone.

A partial F-statistic test reveals that Γ<sub>2dyn</sub> exhibits the highest partial F-statistic of all accretion rate parameters. Therefore, we retain this parameter in the reduced dimensional linear regression. The final model includes the following parameters Π ⊂ {a<sub>1/2</sub>, c<sub>vir</sub>, T/U, Γ<sub>2dyn</sub>, X<sub>off</sub>}. To ensure the robustness of the linear model across all bins, we perform the following tests:

(i) VIF test with a cut-off of 5 to detect multicollinearity.

(ii) Global F-statistic for the entire model with a significance level of 0.05 to examine the correlation between the dependent variable and all parameters.

(iii) Partial F-statistic for the entire parameter set to compare the relative importance of each parameter. The Partial F-statistic measures the additional contribution of each parameter to the multilinear fit by estimating its corresponding R<sup>2</sup> value.

(iv) T-statistic for each parameter to verify that the coefficients significantly deviate from zero at a significance level of 0.05.

(v) Breusch-Pagan Lagrange Multiplier test (Breusch & Pagan 1979) at a significance level of 0.05 to assess heteroscedasticity.

(vi) Shapiro-Wilk test (Shapiro & Wilk 1965) at a significance level of 0.05 to evaluate the Gaussianity of residuals.
<table><caption>**Table D2.** Summary statistics for Cov(ΔΣ, N<sub>gal</sub>|M, z) binned by R<sub>200c</sub> and M<sub>200c</sub> with N<sub>gal</sub> defined inside the halo R<sub>200c</sub> for the reduced error function model. Compared with the full error function model, the performance of the reduced model varies from bin to bin – using ΔDIC > 3 as a statistically significant result, it outperforms the full model in 5/9 overlapping bins, under-performs in 3/9 bins, and is comparable in 2 bins. The reduced model is able to yield convergent chains for M<sub>200c</sub> ∈ [5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>) at z = 0.49 and M<sub>200c</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>15</sup>) at z = 1.03 but with poor constraints on the parameters.</caption><thead><tr><th rowspan="2">Mass & redshift</th><th colspan="5">Error function (reduced): s (erf(√π/2 x̄) - 1)</th></tr><tr><th>τ</th><th>γ</th><th>10<sup>12</sup> × s</th><th>ΔDIC<sub>erf-full</sub></th><th>p-value</th></tr></thead><tbody><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>], z = 0.00</td><td>0.44<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.63<sup>+0.07</sup><sub>-0.09</sub></td><td>2.66<sup>+0.56</sup><sub>-0.39</sub></td><td>-5.2</td><td>0.06</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>], z = 0.00</td><td>0.69<sup>+0.11</sup><sub>-0.09</sub></td><td>-0.98<sup>+0.20</sup><sub>-0.27</sub></td><td>4.10<sup>+2.31</sup><sub>-1.19</sub></td><td>4.5</td><td>0.96</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>], z = 0.00</td><td>0.40<sup>+0.04</sup><sub>-0.03</sub></td><td>-0.50<sup>+0.05</sup><sub>-0.06</sub></td><td>1.90<sup>+0.26</sup><sub>-0.20</sub></td><td>0.2</td><td>0.78</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>], z = 0.00</td><td>0.29<sup>+0.04</sup><sub>-0.04</sub></td><td>-0.40<sup>+0.04</sup><sub>-0.04</sub></td><td>1.12<sup>+0.13</sup><sub>-0.11</sub></td><td>6.1</td><td>0.005</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>], z = 0.49</td><td>0.25<sup>+0.04</sup><sub>-0.03</sub></td><td>-0.35<sup>+0.04</sup><sub>-0.04</sub></td><td>0.79<sup>+0.08</sup><sub>-0.07</sub></td><td>6.2</td><td>0.57</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>], z = 0.49</td><td>0.38<sup>+0.14</sup><sub>-0.09</sub></td><td>-0.48<sup>+0.12</sup><sub>-0.24</sub></td><td>0.85<sup>+0.46</sup><sub>-0.19</sub></td><td>-6.9</td><td>0.81</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>], z = 0.49</td><td>0.25<sup>+0.21</sup><sub>-0.12</sub></td><td>-0.55<sup>+0.14</sup><sub>-0.42</sub></td><td>0.40<sup>+0.69</sup><sub>-0.13</sub></td><td>NA</td><td>0.48</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>], z = 0.49</td><td>0.61<sup>+0.27</sup><sub>-0.24</sub></td><td>-0.90<sup>+0.49</sup><sub>-0.76</sub></td><td>1.01<sup>+3.05</sup><sub>-0.60</sub></td><td>-5.7</td><td>0.99</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>], z = 1.03</td><td>0.20<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.35<sup>+0.06</sup><sub>-0.08</sub></td><td>0.38<sup>+0.07</sup><sub>-0.06</sub></td><td>-0.75</td><td>0.04</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>], z = 1.03</td><td>0.22<sup>+0.06</sup><sub>-0.05</sub></td><td>-0.42<sup>+0.06</sup><sub>-0.08</sub></td><td>0.54<sup>+0.11</sup><sub>-0.09</sub></td><td>3.9</td><td>0.3</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>], z = 1.03</td><td>5.97<sup>+2.76</sup><sub>-3.01</sub></td><td>-0.55<sup>+3.30</sup><sub>-2.99</sub></td><td>0.02<sup>+0.02</sup><sub>-0.008</sub></td><td>NA</td><td>0.15</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure D1.** Posterior distribution of shape parameters in a benchmark bin of M<sub>200c</sub>/M<sub>vir</sub> ∈ [2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>) at z = 0.00 under different binning schemes r<sub>p</sub> and N<sub>gal</sub> models with different halo boundaries R<sub>halo</sub>. The marginalized parameter constraints for the full model closely overlap one another, and using the reduced model with g = −1 marginally improves the posterior constraints. The plot was generated using pygtc (Bocquet & Carter 2016).</figcaption></figure>

## APPENDIX E: MODELLING SECONDARY PROPERTIES

We describe the linear regression model for richness used in Section 6.1. The same methodology is applied to ΔΣ in Section 6.2.

To model the expected natural logarithm of galaxy count (ln N<sub>gal</sub>), we decompose it linearly using secondary halo parameters listed in Table 3, as shown in equation (20). We employ the least squares method for linear regression and examine parameter redundancies. For over half of the bins, the parameters Γ<sub>inst</sub>, Γ<sub>100Myr</sub>, Γ<sub>2dyn</sub>, and Γ<sub>peak</sub> exhibit collinearity, with VIF exceeding 5. This outcome is expected,

as these quantities represent the same physical quantities smoothed over different time-scales. As for the reduced set of non-collinear parameters, their correlation coefficient are quantified in Shin & Diemer (2023) using the Erebos simulation suite. To determine which parameters to retain, we utilize the partial F-statistic.

Table E1 demonstrates the diminishing explanatory power of Π on the richness as seen by the diminishing R<sup>2</sup> and F<sub>partial</sub>. We consider the partial F-statistic serves as a heuristic measure for the explanatory power of a variable, defined as

$$F_{\text{partial}} = \frac{(\text{RSSE}_{\text{reduced}} - \text{RSSE}_{\text{full}})/p}{\text{RSSE}_{\text{full}}/(n-k)}, \quad \text{(E1)}$$

where RSSE is the residual sum of squared errors for the reduced model after removing the parameter in question and the full model containing Π ⊂ {a<sub>1/2</sub>, c<sub>vir</sub>, T/U, Γ<sub>2dyn</sub>, X<sub>off</sub>}, p is the number of parameters removed from the full model which in our case is by construction set to p = 1, n the number of data points, and k is the number of parameters in the full backward model. This statistic can be shown to be proportional to the contribution to the total R<sup>2</sup> uniquely explained by this parameter alone.

A partial F-statistic test reveals that Γ<sub>2dyn</sub> exhibits the highest partial F-statistic of all accretion rate parameters. Therefore, we retain this parameter in the reduced dimensional linear regression. The final model includes the following parameters Π ⊂ {a<sub>1/2</sub>, c<sub>vir</sub>, T/U, Γ<sub>2dyn</sub>, X<sub>off</sub>}. To ensure the robustness of the linear model across all bins, we perform the following tests:

(i) VIF test with a cut-off of 5 to detect multicollinearity.

(ii) Global F-statistic for the entire model with a significance level of 0.05 to examine the correlation between the dependent variable and all parameters.

(iii) Partial F-statistic for the entire parameter set to compare the relative importance of each parameter. The Partial F-statistic measures the additional contribution of each parameter to the multilinear fit by estimating its corresponding R<sup>2</sup> value.

(iv) T-statistic for each parameter to verify that the coefficients significantly deviate from zero at a significance level of 0.05.

(v) Breusch-Pagan Lagrange Multiplier test (Breusch & Pagan 1979) at a significance level of 0.05 to assess heteroscedasticity.

(vi) Shapiro-Wilk test (Shapiro & Wilk 1965) at a significance level of 0.05 to evaluate the Gaussianity of residuals.
<table><caption>**Table E1.** Best-fitting parameters, global $R^2$, and explanatory power indicators for log-richness modelled in equation (20). Values in parentheses represent $1\sigma$ confidence intervals to the partial slopes $\beta$. The partial F-statistic, defined in equation (E1), is used to quantify the explanatory power of each variable. A higher partial F-statistic indicates a greater amount of predictive power uniquely attributed to that variable. Statistical significance is determined by an F-statistic of $F > 10$.</caption><thead><tr><th rowspan="2">Redshift<br>& M<sub>200c</sub> (M<sub>☉</sub>h<sup>-1</sup>)<br>z = 0.00</th><th rowspan="2">R<sup>2</sup></th><th rowspan="2">Const</th><th colspan="2">c<sub>vir</sub></th><th colspan="2">T/U</th><th colspan="2">a<sub>1/2</sub></th><th colspan="2">Γ<sup>*</sup><sub>2dyn</sub></th><th colspan="2">X<sub>off</sub></th></tr><tr><th>β<sub>Ngal</sub></th><th>F</th><th>β<sub>Ngal</sub></th><th>F</th><th>β<sub>Ngal</sub></th><th>F</th><th>10<sup>6</sup> × β<sub>Ngal</sub></th><th>F</th><th>10<sup>3</sup> × β<sub>Ngal</sub></th><th>F</th></tr></thead><tbody><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>)</td><td>0.45</td><td>-0.69</td><td>-0.035(0.004)</td><td>171</td><td>1.448(0.135)</td><td>185</td><td>0.103(0.084)</td><td>166</td><td>2.743(1.68)</td><td>134</td><td>-1.62(0.23)</td><td>69</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>)</td><td>0.48</td><td>-0.46</td><td>-0.033(0.003)</td><td>217</td><td>1.086(0.120)</td><td>179</td><td>-0.198(0.071)</td><td>170</td><td>2.668(0.628)</td><td>206</td><td>-0.76(0.15)</td><td>86</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>)</td><td>0.49</td><td>-0.41</td><td>-0.029(0.003)</td><td>196</td><td>-0.794(0.086)</td><td>145</td><td>0.292(0.007)</td><td>213</td><td>-0.722(0.223)</td><td>167</td><td>0.57(0.11)</td><td>67</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>)</td><td>0.46</td><td>-0.32</td><td>-0.019(0.002)</td><td>119</td><td>1.503(0.068)</td><td>111</td><td>0.148(0.054)</td><td>201</td><td>0.483(0.096)</td><td>226</td><td>-0.25(0.06)</td><td>51</td></tr><tr><td colspan="14">z = 0.49</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>)</td><td>0.25</td><td>-0.71</td><td>0.006(0.005)</td><td>24</td><td>0.325(0.125)</td><td>69</td><td>0.188(0.101)</td><td>142</td><td>3.90(1.64)</td><td>102</td><td>-0.3(0.19)</td><td>35</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>)</td><td>0.21</td><td>-0.45</td><td>-0.450(0.082)</td><td>27</td><td>0.314(0.135)</td><td>46</td><td>-0.353(0.158)</td><td>107</td><td>1.80(0.71)</td><td>89</td><td>-0.37(0.13)</td><td>14</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>)</td><td>0.13</td><td>-0.56</td><td>-0.003(0.004)</td><td>0</td><td>0.323(0.110)</td><td>26</td><td>0.670(0.155)</td><td>51</td><td>0.388(0.108)</td><td>26</td><td>-0.27(0.11)</td><td>2</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>)</td><td>0.13</td><td>-0.32</td><td>-0.010(0.008)</td><td>0</td><td>0.007(0.286)</td><td>0</td><td>0.697(.309)</td><td>5</td><td>-0.538(0.331)</td><td>3</td><td>-0.41(.14)</td><td>3</td></tr><tr><td colspan="14">z = 1.03</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>)</td><td>0.27</td><td>-1.12</td><td>0.015(0.005)</td><td>6</td><td>0.471(0.136)</td><td>21</td><td>0.569(0.269)</td><td>116</td><td>2.7(1.2)</td><td>63</td><td>-0.81(0.18)</td><td>0</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>)</td><td>0.22</td><td>-0.84</td><td>0.016(0.004)</td><td>10</td><td>0.374(0.135)</td><td>36</td><td>-0.233(0.234)</td><td>83</td><td>1.076(0.478)</td><td>36</td><td>-0.35(0.12)</td><td>2</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>)</td><td>0.20</td><td>-0.6</td><td>0.0015(0.005)</td><td>10</td><td>0.284(0.185)</td><td>1</td><td>-0.157(0.266)</td><td>15</td><td>0.160(0.213)</td><td>6</td><td>-0.55(0.11)</td><td>13</td></tr></tbody></table>

Across all bins, the reduced model successfully satisfies the first four tests. However, some bins fail the Shapiro-Wilk test due to a negative skew and positive kurtosis. None the less, a visual examination of Q-Q plots indicates that the residuals predominantly follow a Gaussian distribution, except for deviations at the tail ends. Q-Q plots, quantile-quantile plots, are visualization tools used to compare the quantiles of a data set to the quantiles of a theoretical distribution, typically a normal distribution. They provide a visual assessment of how well the data aligns with the assumed distribution.

## APPENDIX F: DERIVATION OF $P(\Delta\Sigma|N_{\text{GAL}}, m, z)$

We demonstrate that $P(\Delta\Sigma|N_{\text{gal}}, M, z)$ can be modelled as a multilinear relation of secondary halo parameters of mean and variance given by equations (29) and (30).

From our assertion that $P(\Delta\Sigma, \bar{\Pi} | M, z)$ is a bivariate normal, the conditional probability $P(\Pi_i|\Delta\Sigma, M, z)$ in each radial bin can be expressed as a normal distribution with mean

$$ \langle \Pi_i | \Delta\Sigma, M, z \rangle = \langle \Pi_i | M, z \rangle + \rho_{\Pi_i - \Delta\Sigma} \frac{\sigma_{\Pi_i}}{\sigma_{\Delta\Sigma}} (\Delta\Sigma - \langle \Delta\Sigma | M, z \rangle), \quad (\text{F1}) $$

and variance

$$ \sigma_{\Pi_i | \Delta\Sigma}^2 = \sigma_{\Pi_i}^2 (1 - \rho_{\Pi_i - \Delta\Sigma}^2). \quad (\text{F2}) $$

Here, we omit the radial dependence $R/R_{200c}$ for all variables and the conditional dependence on $(M, z)$ in the subscripts for $\rho$ and $\sigma$ (i.e. $\sigma_{\Delta\Sigma}$ should be explicitly written as $\sigma_{\Delta\Sigma|M, z}(R)$).

For an independent random variable $Z = X + Y$ with $X$ and $Y$ uncorrelated independent random variables with distributions of the form $\sim N(u, \sigma^2)$ and $\sim N(v, \tau^2)$, respectively, $Z$ is another Gaussian of $N(u + v, \sigma^2 + 2\rho_{X-Y}\sigma\tau + \tau^2)$. In the case that $X$ and $Y$ are correlated, we must introduce a cross term in the variance of probability distribution $Z$, namely $P(Z) \sim N(u + v, \sigma^2 + 2\rho_{X-Y}\sigma\tau + \tau^2)$.

In our specific case, we want to model the distribution $P(\sum_i \beta_{s,i} \Pi_i | \Delta\Sigma, M, z)$ from $P(\Pi_i | \Delta\Sigma, M, z)$ and the correlation here refers to the correction between secondary halo parameter $\rho_{\Pi_i, \Pi_j}$.<fnref n="3" />

From the convolution theorem and the expressions for $P(\Pi_i | \Delta\Sigma, M, z)$ in equations (F1) and (F2), we obtain the expression for $P(\ln N_{\text{gal}} | \Delta\Sigma, M, z)$ as a normal distribution with mean:

$$ \langle \ln N_{\text{gal}} | \Delta\Sigma, M, z \rangle = \langle \ln N_{\text{gal},0} | \Delta\Sigma, M, z \rangle + \sigma_{\Delta\Sigma} \left( \sum_i \frac{\beta_{N_{\text{gal},i}}}{\sigma_{\Pi_i}} \rho_{\Delta\Sigma - \Pi_i} \times (\Pi_i - \langle \Pi_i | M, z \rangle) \right), \quad (\text{F3}) $$

and variance

$$ \sigma_{\ln N_{\text{gal}} | \Delta\Sigma}^2 = \sigma_{N_{\text{gal},0}}^2 + \sum_i \beta_{N_{\text{gal},i}}^2 \sigma_{\Pi_i}^2 (1 - \rho_{\Delta\Sigma - \Pi_i}^2) + \sum_{i,j} \rho_{\Pi_i - \Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j}. \quad (\text{F4}) $$

We now want to derive the scaling relations for $P(\Delta\Sigma|N_{\text{gal}}, M, z)$. From the Bayes theorem,

$$ P(\Delta\Sigma | N_{\text{gal}}, M, z) = P(\ln N_{\text{gal}} | \Delta\Sigma, M, z) \frac{P(\Delta\Sigma | M, z)}{P(\ln N_{\text{gal}} | M, z)}, \quad (\text{F5}) $$

wherein we assume that $P(\ln N_{\text{gal}}|M, z)$ and $P(\Delta\Sigma|M, z)$ can be modelled as normal distributions, then $P(\Delta\Sigma|N_{\text{gal}}, M, z)$ is another normal distribution with mean

$$ \langle \Delta\Sigma | N_{\text{gal}}, M, z \rangle = \langle \Delta\Sigma | N_{\text{gal},0}, M, z \rangle + C_1 \sigma_{\Delta\Sigma} \left( \sum_i \frac{\beta_{N_{\text{gal},i}}}{\sigma_{\Pi_i}} \rho_{\Delta\Sigma - \Pi_i} \times (\Pi_i - \langle \Pi_i | M, z \rangle) \right) \quad (\text{F6}) $$

and variance

$$ \sigma_{\Delta\Sigma | \ln N_{\text{gal}}}^2 = \sigma_0^2 + C_2 \sum_i \beta_{N_{\text{gal},i}}^2 \sigma_{\Pi_i}^2 (1 - \rho_{\Delta\Sigma - \Pi_i}^2) + C_3 \sum_{i,j} \rho_{\Pi_i - \Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j}. \quad (\text{F7}) $$

The parameter $C_1$ for the mean relation can be explicitly derived if we know the posterior distribution of $P(\ln N_{\text{gal}}|M, z)$ and $P(\Delta\Sigma|M, z)$ by the exercise of completing the squares inside the exponents, i.e. by matching the quadratic, linear and constant terms inside the

<footnote n="3">³In Section 6.2, we show through the variance inflation test that set of reduced model parameters $\Pi \in \{a_{1/2}, c_{\text{vir}}, T/U, \Gamma_{2\text{dyn}}^*, X_{\text{off}}\}$ are not multicollinear.</footnote>
<table><caption>**Table E1.** Best-fitting parameters, global $R^2$, and explanatory power indicators for log-richness modelled in equation (20). Values in parentheses represent $1\sigma$ confidence intervals to the partial slopes $\beta$. The partial F-statistic, defined in equation (E1), is used to quantify the explanatory power of each variable. A higher partial F-statistic indicates a greater amount of predictive power uniquely attributed to that variable. Statistical significance is determined by an F-statistic of $F > 10$.</caption><thead><tr><th rowspan="2">Redshift<br>& M<sub>200c</sub> (M<sub>☉</sub>h<sup>-1</sup>)<br>z = 0.00</th><th rowspan="2">R<sup>2</sup></th><th rowspan="2">Const</th><th colspan="2">c<sub>vir</sub></th><th colspan="2">T/U</th><th colspan="2">a<sub>1/2</sub></th><th colspan="2">Γ<sup>*</sup><sub>2dyn</sub></th><th colspan="2">X<sub>off</sub></th></tr><tr><th>β<sub>Ngal</sub></th><th>F</th><th>β<sub>Ngal</sub></th><th>F</th><th>β<sub>Ngal</sub></th><th>F</th><th>10<sup>6</sup> × β<sub>Ngal</sub></th><th>F</th><th>10<sup>3</sup> × β<sub>Ngal</sub></th><th>F</th></tr></thead><tbody><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>)</td><td>0.45</td><td>-0.69</td><td>-0.035(0.004)</td><td>171</td><td>1.448(0.135)</td><td>185</td><td>0.103(0.084)</td><td>166</td><td>2.743(1.68)</td><td>134</td><td>-1.62(0.23)</td><td>69</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>)</td><td>0.48</td><td>-0.46</td><td>-0.033(0.003)</td><td>217</td><td>1.086(0.120)</td><td>179</td><td>-0.198(0.071)</td><td>170</td><td>2.668(0.628)</td><td>206</td><td>-0.76(0.15)</td><td>86</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>)</td><td>0.49</td><td>-0.41</td><td>-0.029(0.003)</td><td>196</td><td>-0.794(0.086)</td><td>145</td><td>0.292(0.007)</td><td>213</td><td>-0.722(0.223)</td><td>167</td><td>0.57(0.11)</td><td>67</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>)</td><td>0.46</td><td>-0.32</td><td>-0.019(0.002)</td><td>119</td><td>1.503(0.068)</td><td>111</td><td>0.148(0.054)</td><td>201</td><td>0.483(0.096)</td><td>226</td><td>-0.25(0.06)</td><td>51</td></tr><tr><td colspan="14">z = 0.49</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>)</td><td>0.25</td><td>-0.71</td><td>0.006(0.005)</td><td>24</td><td>0.325(0.125)</td><td>69</td><td>0.188(0.101)</td><td>142</td><td>3.90(1.64)</td><td>102</td><td>-0.3(0.19)</td><td>35</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>)</td><td>0.21</td><td>-0.45</td><td>-0.450(0.082)</td><td>27</td><td>0.314(0.135)</td><td>46</td><td>-0.353(0.158)</td><td>107</td><td>1.80(0.71)</td><td>89</td><td>-0.37(0.13)</td><td>14</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>)</td><td>0.13</td><td>-0.56</td><td>-0.003(0.004)</td><td>0</td><td>0.323(0.110)</td><td>26</td><td>0.670(0.155)</td><td>51</td><td>0.388(0.108)</td><td>26</td><td>-0.27(0.11)</td><td>2</td></tr><tr><td>[5 × 10<sup>14</sup>, 1 × 10<sup>15</sup>)</td><td>0.13</td><td>-0.32</td><td>-0.010(0.008)</td><td>0</td><td>0.007(0.286)</td><td>0</td><td>0.697(.309)</td><td>5</td><td>-0.538(0.331)</td><td>3</td><td>-0.41(.14)</td><td>3</td></tr><tr><td colspan="14">z = 1.03</td></tr><tr><td>[5 × 10<sup>13</sup>, 1 × 10<sup>14</sup>)</td><td>0.27</td><td>-1.12</td><td>0.015(0.005)</td><td>6</td><td>0.471(0.136)</td><td>21</td><td>0.569(0.269)</td><td>116</td><td>2.7(1.2)</td><td>63</td><td>-0.81(0.18)</td><td>0</td></tr><tr><td>[1 × 10<sup>14</sup>, 2 × 10<sup>14</sup>)</td><td>0.22</td><td>-0.84</td><td>0.016(0.004)</td><td>10</td><td>0.374(0.135)</td><td>36</td><td>-0.233(0.234)</td><td>83</td><td>1.076(0.478)</td><td>36</td><td>-0.35(0.12)</td><td>2</td></tr><tr><td>[2 × 10<sup>14</sup>, 5 × 10<sup>14</sup>)</td><td>0.20</td><td>-0.6</td><td>0.0015(0.005)</td><td>10</td><td>0.284(0.185)</td><td>1</td><td>-0.157(0.266)</td><td>15</td><td>0.160(0.213)</td><td>6</td><td>-0.55(0.11)</td><td>13</td></tr></tbody></table>

Across all bins, the reduced model successfully satisfies the first four tests. However, some bins fail the Shapiro-Wilk test due to a negative skew and positive kurtosis. None the less, a visual examination of Q-Q plots indicates that the residuals predominantly follow a Gaussian distribution, except for deviations at the tail ends. Q-Q plots, quantile-quantile plots, are visualization tools used to compare the quantiles of a data set to the quantiles of a theoretical distribution, typically a normal distribution. They provide a visual assessment of how well the data aligns with the assumed distribution.

## APPENDIX F: DERIVATION OF $P(\Delta\Sigma|N_{\text{GAL}}, m, z)$

We demonstrate that $P(\Delta\Sigma|N_{\text{gal}}, M, z)$ can be modelled as a multilinear relation of secondary halo parameters of mean and variance given by equations (29) and (30).

From our assertion that $P(\Delta\Sigma, \bar{\Pi} | M, z)$ is a bivariate normal, the conditional probability $P(\Pi_i|\Delta\Sigma, M, z)$ in each radial bin can be expressed as a normal distribution with mean

$$ \langle \Pi_i | \Delta\Sigma, M, z \rangle = \langle \Pi_i | M, z \rangle + \rho_{\Pi_i - \Delta\Sigma} \frac{\sigma_{\Pi_i}}{\sigma_{\Delta\Sigma}} (\Delta\Sigma - \langle \Delta\Sigma | M, z \rangle), \quad (\text{F1}) $$

and variance

$$ \sigma_{\Pi_i | \Delta\Sigma}^2 = \sigma_{\Pi_i}^2 (1 - \rho_{\Pi_i - \Delta\Sigma}^2). \quad (\text{F2}) $$

Here, we omit the radial dependence $R/R_{200c}$ for all variables and the conditional dependence on $(M, z)$ in the subscripts for $\rho$ and $\sigma$ (i.e. $\sigma_{\Delta\Sigma}$ should be explicitly written as $\sigma_{\Delta\Sigma|M, z}(R)$).

For an independent random variable $Z = X + Y$ with $X$ and $Y$ uncorrelated independent random variables with distributions of the form $\sim N(u, \sigma^2)$ and $\sim N(v, \tau^2)$, respectively, $Z$ is another Gaussian of $N(u + v, \sigma^2 + 2\rho_{X-Y}\sigma\tau + \tau^2)$. In the case that $X$ and $Y$ are correlated, we must introduce a cross term in the variance of probability distribution $Z$, namely $P(Z) \sim N(u + v, \sigma^2 + 2\rho_{X-Y}\sigma\tau + \tau^2)$.

In our specific case, we want to model the distribution $P(\sum_i \beta_{s,i} \Pi_i | \Delta\Sigma, M, z)$ from $P(\Pi_i | \Delta\Sigma, M, z)$ and the correlation here refers to the correction between secondary halo parameter $\rho_{\Pi_i, \Pi_j}$.<fnref n="3" />

From the convolution theorem and the expressions for $P(\Pi_i | \Delta\Sigma, M, z)$ in equations (F1) and (F2), we obtain the expression for $P(\ln N_{\text{gal}} | \Delta\Sigma, M, z)$ as a normal distribution with mean:

$$ \langle \ln N_{\text{gal}} | \Delta\Sigma, M, z \rangle = \langle \ln N_{\text{gal},0} | \Delta\Sigma, M, z \rangle + \sigma_{\Delta\Sigma} \left( \sum_i \frac{\beta_{N_{\text{gal},i}}}{\sigma_{\Pi_i}} \rho_{\Delta\Sigma - \Pi_i} \times (\Pi_i - \langle \Pi_i | M, z \rangle) \right), \quad (\text{F3}) $$

and variance

$$ \sigma_{\ln N_{\text{gal}} | \Delta\Sigma}^2 = \sigma_{N_{\text{gal},0}}^2 + \sum_i \beta_{N_{\text{gal},i}}^2 \sigma_{\Pi_i}^2 (1 - \rho_{\Delta\Sigma - \Pi_i}^2) + \sum_{i,j} \rho_{\Pi_i - \Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j}. \quad (\text{F4}) $$

We now want to derive the scaling relations for $P(\Delta\Sigma|N_{\text{gal}}, M, z)$. From the Bayes theorem,

$$ P(\Delta\Sigma | N_{\text{gal}}, M, z) = P(\ln N_{\text{gal}} | \Delta\Sigma, M, z) \frac{P(\Delta\Sigma | M, z)}{P(\ln N_{\text{gal}} | M, z)}, \quad (\text{F5}) $$

wherein we assume that $P(\ln N_{\text{gal}}|M, z)$ and $P(\Delta\Sigma|M, z)$ can be modelled as normal distributions, then $P(\Delta\Sigma|N_{\text{gal}}, M, z)$ is another normal distribution with mean

$$ \langle \Delta\Sigma | N_{\text{gal}}, M, z \rangle = \langle \Delta\Sigma | N_{\text{gal},0}, M, z \rangle + C_1 \sigma_{\Delta\Sigma} \left( \sum_i \frac{\beta_{N_{\text{gal},i}}}{\sigma_{\Pi_i}} \rho_{\Delta\Sigma - \Pi_i} \times (\Pi_i - \langle \Pi_i | M, z \rangle) \right) \quad (\text{F6}) $$

and variance

$$ \sigma_{\Delta\Sigma | \ln N_{\text{gal}}}^2 = \sigma_0^2 + C_2 \sum_i \beta_{N_{\text{gal},i}}^2 \sigma_{\Pi_i}^2 (1 - \rho_{\Delta\Sigma - \Pi_i}^2) + C_3 \sum_{i,j} \rho_{\Pi_i - \Pi_j} \sigma_{\Pi_i} \sigma_{\Pi_j}. \quad (\text{F7}) $$

The parameter $C_1$ for the mean relation can be explicitly derived if we know the posterior distribution of $P(\ln N_{\text{gal}}|M, z)$ and $P(\Delta\Sigma|M, z)$ by the exercise of completing the squares inside the exponents, i.e. by matching the quadratic, linear and constant terms inside the

<footnote n="3">³In Section 6.2, we show through the variance inflation test that set of reduced model parameters $\Pi \in \{a_{1/2}, c_{\text{vir}}, T/U, \Gamma_{2\text{dyn}}^*, X_{\text{off}}\}$ are not multicollinear.</footnote>
exponents of normal distributions on the left and right hand sides of equation (F5).

The parameters C₂, C₃, and σ₀ can be explicitly derived first by using the variance of product law of correlated Gaussians, i.e. Var(Z ≡ XY) = 1 + ρ² after transforming X and Y into unit variance, zero mean Gaussians with correlation efficient ρ, and then by using the variance of quotient approximation

where μₓ, σₓ² are the mean and variances of P(ΔΣ|M, z) and μᵧ, σᵧ² the mean and variances of P(ln Ngal|M, z) in our specific case.

The exact values of C₁, C₂, C₃, and σ₀ are not essential for this paper as we aim to derive a general expression for P(ΔΣ|Ngal, M, z) as a function of secondary halo parameters.

$$ \mathrm{Var}(X/Y) = \left(\frac{\mu_X}{\mu_Y}\right)^2 \left[ \frac{\sigma_R^2}{\mu_R^2} - 2 \frac{\mathrm{Cov}(X, Y)}{\mu_X \mu_Y} + \frac{\sigma_X^2}{\mu_Y^2} \right], \quad (\mathrm{F8}) $$

This paper has been typeset from a TeX/Latex file prepared by the author.
exponents of normal distributions on the left and right hand sides of equation (F5).

The parameters C₂, C₃, and σ₀ can be explicitly derived first by using the variance of product law of correlated Gaussians, i.e. Var(Z ≡ XY) = 1 + ρ² after transforming X and Y into unit variance, zero mean Gaussians with correlation efficient ρ, and then by using the variance of quotient approximation

where μₓ, σₓ² are the mean and variances of P(ΔΣ|M, z) and μᵧ, σᵧ² the mean and variances of P(ln Ngal|M, z) in our specific case.

The exact values of C₁, C₂, C₃, and σ₀ are not essential for this paper as we aim to derive a general expression for P(ΔΣ|Ngal, M, z) as a function of secondary halo parameters.

$$ \mathrm{Var}(X/Y) = \left(\frac{\mu_X}{\mu_Y}\right)^2 \left[\frac{\sigma_R^2}{\mu_R^2} - 2 \frac{\mathrm{Cov}(X, Y)}{\mu_X \mu_Y} + \frac{\sigma_X^2}{\mu_Y^2}\right], \quad (F8) $$

This paper has been typeset from a TeX/Latex file prepared by the author.