© The Authors 2025

1. Introduction

Observations of the large-scale distribution of galaxies are among the most promising probes for accurately inferring the cosmological parameters of our Universe. Past large-scale structure surveys have helped shape our current model of the Universe by providing precise measurements of the angular diameter distance and growth rate at different epochs (see e.g. Alam et al. 2021). Forthcoming surveys such as the Dark Energy Instrument (DESI, Levi et al. 2013) and Euclid (Amendola et al. 2018) will measure the sky positions and redshifts (both photometric and spectroscopic) of an unprecedented amount of galaxies. To reliably interpret these datasets, it is of crucial importance to accurately model the spatial distributions of galaxies as a function of cosmology.

The evolution of the large-scale distribution of matter is predominantly driven by gravity and can be predicted very reliably through perturbation theory (see Bernardeau et al. 2002, for a review) or by gravity-only N-body simulations (see e.g. Angulo & Hahn 2022, for a review). However, observed galaxies trace the matter distribution only in a ‘biased’ way (Kaiser 1984; Mo & White 1996). Optimally exploiting the information from galaxy surveys therefore requires not only accurately modelling gravity, but also to account for the formation and evolution of galaxies.

There are a variety of methodologies to model the formation of galaxies. The most detailed approach involves employing hydrodynamical simulations to explicitly follow gas dynamics and the formation and evolution of stars, black holes, and galaxies (see e.g. Vogelsberger et al. 2014; Dubois et al. 2014; Schaye et al. 2015; Davé et al. 2019). In principle, such simulations would be the ideal method to account for the formation of galaxies as long as the underlying physics could be modelled reliably and at an affordable cost. However, in practice, these simulations are limited in volume and number due to their substantial computational requirements. Additionally, they have required the assumption of sub-grid physics; for instance, to model the unresolved formation of stars, as well as the growth of active galactic nuclei and their feedback processes. The results of such simulations depend strongly on the associated parameters and are therefore severely limited in their predictive ability (e.g. Genel et al. 2019; Villaescusa-Navarro et al. 2021).

To overcome the limitations of hydrodynamic simulations, a more agnostic approach can be pursued. Semi-analytic models, for example, follow a number of empirical and semi-empirical relations that provide more versatility to the modelling of galaxy formation physics (see e.g. Kauffmann et al. 1999; Henriques et al. 2015; Stevens et al. 2016; Lacey et al. 2016; Croton et al. 2016; Lagos et al. 2018). An even more agnostic strategy involves techniques such as sub-halo abundance matching (see e.g. Conroy et al. 2006; Reddick et al. 2013; Chaves-Montero et al. 2016; Lehmann et al. 2017; Dragomir et al. 2018; Contreras et al. 2021; Ortega-Martinez et al. 2024) or halo occupation distribution (see e.g.Peacock & Smith 2000; Berlind & Weinberg 2002; Zheng et al. 2005; Cacciato et al. 2012; Salcedo et al. 2022a), which rely solely on gravity simulations and fundamental assumptions about populating collapsed structures (halos and sub-halos) based on their mass. These assumptions must be sufficiently flexible to accommodate various galaxy formation scenarios. Nonetheless, these techniques also have computational constraints and the principles they rely on may oversimplify reality, necessitating extensions to account for dependencies that go beyond solely the mass of the collapsed structure (a phenomenon known as assembly bias; Gao et al. 2005; Wechsler et al. 2006; Gao & White 2007; Croton et al. 2007; Dalal et al. 2008; Faltenbacher & White 2010; Montero-Dorta et al. 2017; Zehavi et al. 2018; Ferreras et al. 2019; Sato-Polito et al. 2019; Tucci et al. 2021; Salcedo et al. 2022b; Chaves-Montero et al. 2023, among others).

The most general approach currently known is a perturbative expansion of galaxy bias (see Desjacques et al. 2018 for a review). In this framework, the evolution of matter is assumed to be modelled accurately through purely gravitational effects, whereas the galaxy field selectively populates the matter field. The galaxy field is then written in a perturbative manner as a function of the properties of the underlying matter distribution. In Eulerian bias schemes, the galaxy number density is expanded in terms of the final properties of the density field; whereas in Lagrangian schemes, this is done in terms of the initial properties of the linear density field. This approach offers great flexibility, allowing for a single model to describe biased tracers with vastly different properties. Bias approaches based on effective field theory are aimed at describing the clustering behaviour of galaxies down to k ≈ 0.2h/Mpc (Baumann et al. 2012; Baldauf et al. 2016a; Vlah et al. 2016). They have been used extensively to extract robust cosmological constraints from surveys (Ivanov et al. 2020; d’Amico et al. 2020; Colas et al. 2020; Nishimichi et al. 2020; Chen et al. 2020; Philcox & Ivanov 2022).

There exist different options for defining the biasing scheme. Traditionally, perturbation theory is used to describe the underlying gravitational evolution of the density field; whereas in the recently developed ‘hybrid approaches’, the gravitational evolution is treated exactly through N-body simulations (Modi et al. 2020; Kokron et al. 2021; Zennaro et al. 2023; Pellejero Ibañez et al. 2022, 2023; DeRose et al. 2023). However, in all bias methods the galaxy number density is expanded in a perturbative series with a set of free coefficients known as ‘bias parameters’. These parameters can be interpreted as the response of the galaxy number density to perturbations of the density field. While these parameters offer some degree of physical insights (Lazeyras et al. 2016; Lazeyras & Schmidt 2019; Barreira et al. 2021), they are generally treated as nuisance parameters that can be marginalised over to extract the cosmological information of interest.

A quantitative understanding of the bias parameters is important for analysing future galaxy surveys. On the one hand, this is necessary to determine how well the bias expansion converges. The question of whether the truncation at a given order gives accurate results depends on higher order terms being small enough that they can be neglected. On the other hand, it is of significant interest to limit the prior volume that is used when fitting to cosmological surveys to maximise the extracted cosmological information. In particular, it may significantly benefit an analysis to fix bias parameters to so called ‘co-evolution relations’ that relate the parameters to each other. Therefore, significant efforts have been placed on measuring the bias parameters and to constrain coevolution relations in simulations.

One of the most accurate methods for measuring bias parameters is based on ‘separate-universe’ simulations (Li et al. 2014; Wagner et al. 2015). In these cases, we can directly test how the number of tracers responds to an increase in the large-scale density. This technique has been used to constrain the density bias parameters of halos to very high accuracy (Lazeyras et al. 2016). Furthermore, similar ideas have been used to measure the response to non-Gaussian perturbations (Barreira 2020) and to the changes in the Laplacian of the density field (Lazeyras & Schmidt 2018). Other bias parameters have been constrained through different methods. For example, the tidal bias can be been measured through Fourier space correlations (Modi et al. 2017; Lazeyras & Schmidt 2018); alternatively, all bias parameters can be constrained simultaneously through forward modeling of the power spectrum (e.g. Zennaro et al. 2022) or the field-level galaxy distribution (e.g. Lazeyras et al. 2021).

Another successful avenue of measuring and understanding the behaviour of bias parameters is through peak theory (Bardeen et al. 1986) where the number density of peaks of the initial density field is investigated as a function of a Lagrangian smoothing scale. If galaxies (and halos) correspond to peaks of the initial density field, then bias parameters can be predicted through the response of peaks to large scale perturbations. While the mapping between structures and peaks bears significant uncertainty due to the necessary inclusion of a heuristic smoothing scale to define peaks and significant difficulty with treating the cloud-in-cloud problem (Bardeen et al. 1986), peak theory has proven as a useful tool in various aspects of biasing (Desjacques et al. 2018). For example, it predicts the importance of a Laplacian bias parameter (Bardeen et al. 1986), the scale dependence of the bias in the matter correlation function (Desjacques et al. 2010; Paranjape et al. 2013b,a) and non-zero velocity bias (Baldauf et al. 2015). All of these predictions have been established through measurements in simulations. Furthermore, peak theory and excursion set approaches have motivated the possibility of measuring biases through correlations between density and halos in Lagrangian space (Musso et al. 2012; Paranjape et al. 2013b,a; Biagetti et al. 2014). In particular, Paranjape et al. (2013a) have shown that it is possible to recover precise measurements of (scale-independent) large-scale bias parameters by mapping the scale-dependence of ‘naive’ bias parameters under the assumption of Peak theory to their large scale limit.

In this paper, we propose a novel approach for accurately measuring Lagrangian bias parameters from simulations. Here, we consider the distribution of the linear density field at the initial (Lagrangian) locations of galaxies, which we call the ‘galaxy environment distribution’. We adopted a probabilistic approach to model this galaxy environment distribution and show that the large-scale bias parameters have a simple relationship to the moments of this distribution. We used this to derive estimators of the bias parameters that can take into account spatial corrections at any order to practically eliminate any scale dependence. We derived such estimators for a broad variety of bias parameters and provide corresponding measurements for the case of halos. Operationally, our method is similar to aforementioned peak theory approaches. However, it is (by design) independent of the assumption that galaxies form in smoothed Lagrangian density peaks and it can be used to correct bias measurements at arbitrary high spatial orders.

Furthermore, we propose a new set of ‘cumulant bias parameters’, which are defined as the response of the logarithm of the galaxy number density to perturbations in the linear field. We show that these parameters have significantly improved properties when compared to their canonical counterparts. In the case of halos, we find that cumulant biases at the order of 3 and higher are consistent with zero. Therefore, the canonical co-evolution relations for halos at the order of 3 and beyond primarily appear as an effect of a sub-optimal parameterisation. We suggest that rephrasing the bias expansion in terms of the cumulant bias parameters can significantly enhance its convergence.

The vanishing of cumulants for halos at the order of 3 and beyond implies that the bias function can be well approximated through a Gaussian. We will explore such a Gaussian bias model in our companion paper (Stücker et al. 2025).

This article is organised as follows: In Sect. 2, we explain the probabilistic approach to measure Lagrangian bias parameters of scalar variables such as the density and the Laplacian. We introduce the concept of the cumulant bias expansion. In Sect. 3, we present measurements of the corresponding parameters and demonstrate how cumulant bias parameters exhibit several practical advantages. In Sect. 4, we show how to generalise the concept of the probabilistic estimators for tensorial bias variables, such as the tidal bias, b_K². In Sect. 5, we provide measurements of a few select tensorial quantities. In Sect. 6, we briefly discuss the quantitative importance of different bias terms. Finally, in Sect. 7, we summarise the benefits of our novel estimators and we discuss the circumstances where it is advantageous to phrase the bias expansion in terms of the cumulant bias parameters.

2. Theory

In this section, we introduce the necessary theory and (1) relate the bias parameters to properties of the galaxy environment distribution; (2) express the moment and cumulant generating functions of the galaxy environment distribution in terms of the large-scale bias functio; (3) introduce the concept of a cumulant bias expansion; (4) and show how to derive estimators for both the canonical and cumulant bias parameters.

2.1. Definitions

Our considerations are based on the idea of the peak-background split (PBS) which states that ‘a long-wavelength density perturbation acts as a local modification of the background density for the purposes of the formation of halos and galaxies’ (Kaiser 1984; Bardeen et al. 1986; Desjacques et al. 2018). Bias parameters describe the response of the galaxy number to such long-wavelength perturbations. In this paper, we focus exclusively on Lagrangian bias parameters as the response to perturbations in the linear density field.

An exact implementation of the PBS is given by the separate-universe approach. In this approach, we consider a generic universe with background density, ρ_bg, 0, whereby galaxies form with an average number density, n_g, 0. If we were to increase the initial background density of the universe by a linear amount δ₀ (e.g. in a separate-universe simulation, Frenk et al. 1988; Li et al. 2014; Wagner et al. 2015), in this new universe, galaxies would form with a different average number density, n_g(δ₀)¹. We refer to their ratio as:

$$\begin{array}{cc}\hfill F({\delta }_0)& =\frac{n_g({\delta }_0)}{n_{g,0}}=1+b_1{\delta }_0+\frac{1}{2}{b}_2{\delta }_0^2+\cdots +\frac{1}{n!}{b}_n{\delta }_0^n+\cdots \hfill \end{array}$$(1)

Here, we refer to the bias function or the separate-universe bias function, where δ₀ refers to a contrast in linear densities so that the new background density corresponds to

$$\begin{array}{c}\hfill {\rho }_{\mathrm{bg}}\approx {\rho }_{\mathrm{bg},0}(1+{\delta }_{0}D(a))\text{for}a\to 0,\end{array}$$(2)

where D is the linear growth factor normalised to D(a = 1) = 1 and F can be directly measured with separate-universe simulations (e.g. Lazeyras et al. 2016; Baldauf et al. 2016b). We refer to the coefficients of the indicated expansion as the ‘canonical bias parameters’ or simply ‘the bias parameters’ as:

$$\begin{array}{cc}\hfill b_n& ={\frac{{\partial }^{n}F({\delta }_0)}{\partial {\delta }_0^n}|}_{{\delta }_0=0}.\hfill \end{array}$$(3)

Therefore, the bias parameters physically describe the response of the galaxy density to small perturbations at infinitely large scales. In this paper, we want to investigate galaxy bias from a probabilistic perspective.

We considered an infinitesimally small Lagrangian volume element, about which nothing is known, apart from the linear density contrast, δ, smoothed at some scale (and possibly other features of the linear field like the Laplacian, L, or the tidal field). Neglecting the primordial non-Gaussianity, the density contrast follows a Gaussian distribution:

$$\begin{array}{cc}\hfill p(\delta )& =\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{{\delta }^2}{2{\sigma }^2}).\hfill \end{array}$$(4)

For simplicity, we assume throughout this article that the smoothed density contrast is defined with a sharp k-space filter. Most considerations can be translated to cases filtered in different ways and in a simple manner, but some additional care must be taken due to the more complicated correlation between large and small scales. This is discussed in more detail in Appendix A.

When a sufficiently small volume element is considered, then it is only possible to have either ‘0’ or ‘1’ galaxy. We may therefore speak of a binary event ‘g’ that a volume contains a galaxy. We call the average probability that a galaxy forms in such a volume element ‘p(g)’ and we call the conditional probability, given the knowledge of the linear density contrast, ‘p(g|δ)’². The excess probability,

$$\begin{array}{c}\hfill f(\delta ):=\frac{p(\mathrm{g}|\delta )}{p(\mathrm{g})},\end{array}$$(5)

is parameterised through a function f(δ), which we refer to as the ‘scale-dependent bias function’ or just the ‘bias function’ throughout this article. The bias function depends in a predictable manner on the variance of δ at the considered scale, as we show later in this paper.

Since densities at different scales add up linearly, a separate-universe style modification of the large-scale density contrast from 0 to δ₀ will immediately translate to a modification of the linear density in our volume element δ → δ + δ₀. Therefore, F and f should be related through

$$\begin{array}{cc}\hfill F({\delta }_0)& =⟨f(\delta +{\delta }_0)⟩,\hfill \end{array}$$(6)

where the angled brackets indicate an expectation value taken over the Lagrangian volume (see also Desjacques et al. 2018). The relation indicates that in a separate-universe experiment, the number of galaxies would change according to the average change in probability of forming galaxies when changing the linear density contrast everywhere in space. Therefore, the canonical bias parameters are given in terms of the scale dependent bias function as

$$\begin{array}{cc}\hfill b_n& =\langle \frac{{\partial }^{n}f}{\partial {\delta }^n}\rangle .\hfill \end{array}$$(7)

Later, we show that Equation (6) holds only approximately, since at any finite smoothing scale, the density contrast, δ, is correlated with other variables (e.g. the Laplacian), so that a change in the small-scale density contrast at the location of all galaxies is not exactly equivalent to the ‘pure’ density change in separate-universe experiments. This introduces scale-dependencies that can be accounted for.

Finally, we introduce one further object which we call the ‘galaxy environment distribution’ of p(δ|g) or halo environment distribution when we are specifically talking about halos. This quantifies the probability of the linear density contrast in an infinitesimal volume element, given that there is a galaxy at the considered location. This function can easily be measured as a histogram of the linear density at the Lagrangian locations of galaxies as is illustrated in Figure 1, where we used a damping scale of k_d = 0.15 h Mpc⁻¹ here³, leading to σ = 0.56.

Through Bayes’ theorem, the relation between the galaxy environment distribution and the bias function is given by

$$\begin{array}{cc}\hfill p(\delta |g)& =\frac{p(\delta \cap g)}{p(g)}\hfill \end{array}$$(8)

$$\begin{array}{cc}& =\frac{p(g|\delta )p(\delta )}{p(g)}\hfill \end{array}$$(9)

$$\begin{array}{cc}& =f(\delta )p(\delta ).\hfill \end{array}$$(10)

This means for example, that the bias function, f(δ), can be investigated in a non-parametric way by measuring p(δ|g) and dividing by the Gaussian background distribution, p(δ), as we will investigate in the companion paper (Stücker et al. 2025).

In this section, we explain how we use probability theory to investigate the properties of the galaxy environment distribution. In particular, we show that the bias parameters are simply related to the moments of this distribution and that probability theory motivates the application of better behaved ‘cumulant bias parameters’.

2.2. Bias estimators

Following up on Equation (6), we can write

$$\begin{array}{cc}\hfill F({\delta }_0)& =⟨f(\delta +{\delta }_0)⟩\hfill \\ \hfill & ={\int }_{-\infty }^{\infty }p(\delta )f(\delta +{\delta }_0)\mathrm{d}\delta \hfill \\ \hfill & ={\int }_{-\infty }^{\infty }p(\delta -{\delta }_0)f(\delta )\mathrm{d}\delta .\hfill \end{array}$$(11)

We note that in the last line, we made the substitution δ → δ + δ₀; also, δ represents the linearly extrapolated overdensity and, hence, it can assume values that are between negative and positive infinity. The latter Equation can be interpreted as an alternative perspective onto the separate-universe experiment: when increasing the background density, the probability of environments is changed by a factor p(δ − δ₀)/p(δ), whereas the likelihood of forming a galaxy when presupposing a given environment stays constant.

Now, combining Eqs. (3) and (11), we can evaluate the bias parameters as

$$\begin{array}{cc}\hfill b_n& ={\frac{{\partial }^n}{\partial {\delta }_0^n}{\int }_{-\infty }^{\infty }p(\delta -{\delta }_0)f(\delta )\mathrm{d}\delta |}_{{\delta }_0=0}\hfill \\ \hfill & ={(-1)}^n{\int }_{-\infty }^{\infty }{p}^{(n)}(\delta )f(\delta )\mathrm{d}\delta \hfill \\ \hfill & ={(-1)}^n{\int }_{-\infty }^{\infty }\frac{p^{(n)}(\delta )}{p(\delta )}p(\delta |g)\mathrm{d}\delta \hfill \\ \hfill & ={(-1)}^n{\langle \frac{p^{(n)}(\delta )}{p(\delta )}\rangle }_{\mathrm{g}},\hfill \end{array}$$(12)

where the angled brackets with a ‘g’ subscript indicate an expectation value evaluated over the locations of galaxies (rather than all of Lagrangian space). Furthermore, we can use the fact that p is a Gaussian distribution, for which the derivatives are given by the (probabilisist’s) Hermite polynomials:

$$\begin{array}{cc}\hfill \frac{{\partial }^n}{\partial {\delta }^n}exp(-\frac{{\delta }^2}{2{\sigma }^2})& ={(-1)}^n{\sigma }^{-n}exp(-\frac{{\delta }^2}{2{\sigma }^2})H_n(\delta /\sigma ),\hfill \end{array}$$(13)

so that the bias estimators are:

$$\begin{array}{cc}\hfill b_{n,\mathrm{so}0}& ={\langle \frac{H_n(\delta /\sigma )}{{\sigma }^n}\rangle }_{\mathrm{g}}.\hfill \end{array}$$(14)

Here, the subscript ‘so0’ indicates that these estimators are at a ‘spatial order of 0’. In other words, they do not include corrections from higher spatial derivatives such as the Laplacian (as we explain in more detail later in this work). This expression was already used by Paranjape et al. (2013a,b) to measure the bias parameters. It was motivated by its emergence in excursion set frameworks (Musso et al. 2012) and peak statistics (Paranjape & Sheth 2012), but it is clearly also valid outside such frameworks, requiring us to only assume the PBS. On a related note, it is worth mentioning that Szalay (1988) previously proposed expanding the field in terms of Hermite polynomials.

The estimators for the first four bias parameters are expressed as:

$$\begin{array}{cc}\hfill b_{1,\mathrm{so}0}& ={\langle \frac{\delta }{{\sigma }^2}\rangle }_{\mathrm{g}},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill b_{2,\mathrm{so}0}& ={\langle \frac{{\delta }^2-{\sigma }^2}{{\sigma }^4}\rangle }_{\mathrm{g}},\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill b_{3,\mathrm{so}0}& ={\langle \frac{{\delta }^3-3\delta {\sigma }^2}{{\sigma }^6}\rangle }_{\mathrm{g}},\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill b_{4,\mathrm{so}0}& ={\langle \frac{{\delta }^4-6{\delta }^2{\sigma }^2+3{\sigma }^4}{{\sigma }^8}\rangle }_{\mathrm{g}}.\hfill \end{array}$$(18)

2.3. The moment generating function

Equation (14) shows that there exists a simple relation between moments of the galaxy environment distribution and the bias parameters. We can show this in a very general manner. Expanding the Gaussian background distribution, we obtain

$$\begin{array}{cc}\hfill p(\delta -{\delta }_0)& =\frac{1}{\sqrt{2\pi }\sigma }exp(-\frac{{\delta }^2-2\delta {\delta }_0+{\delta }_0^2}{2{\sigma }^2})\hfill \\ \hfill & =p(\delta )exp\left(\frac{\delta {\delta }_0}{{\sigma }^2}\right)exp(-\frac{{\delta }_0^2}{2{\sigma }^2}).\hfill \end{array}$$

Inserting this into Equation (11) yields

$$\begin{array}{cc}\hfill F({\delta }_0)& =exp(-\frac{{\delta }_0^2}{2{\sigma }^2}){\int }_{-\infty }^{\infty }f(\delta )p(\delta )exp\left(\frac{\delta {\delta }_0}{{\sigma }^2}\right)\mathrm{d}\delta \hfill \\ \hfill & =exp(-\frac{1}{2}{t}^2{\sigma }^2){\int }_{-\infty }^{\infty }p(\delta |g)exp\left(t\delta \right)\mathrm{d}\delta \hfill \\ \hfill & =exp(-\frac{1}{2}{t}^2{\sigma }^2){\langle exp\left(t\delta \right)\rangle }_{\mathrm{g}},\hfill \end{array}$$(19)

where we have labelled t = δ₀/σ². The last term can be identified with the moment generating function,

$$\begin{array}{c}\hfill M_g(t)={\langle exp\left(t\delta \right)\rangle }_{\mathrm{g}},\end{array}$$(20)

of the galaxy environment distribution. Therefore, the moment-generating function of the galaxy environment distribution and the separate-universe bias function can be directly converted into each other as follows:

$$\begin{array}{cc}\hfill M_g(t)& =F(t{\sigma }^2)exp\left(\frac{1}{2}{t}^2{\sigma }^2\right).\hfill \end{array}$$(21)

Furthermore, the relation between the moments and the bias parameters, as given by Equation (14) can equivalently be found by taking derivatives of the moment generating function:

$$\begin{array}{cc}\hfill {\mu }_n& :={\langle {\delta }^n\rangle }_{\mathrm{g}}={\frac{\partial }{\partial t^n}{M}_g(t)|}_{t=0}.\hfill \end{array}$$(22)

It is worth noting that this result may be related to the considerations in White (1979), where in the case of discrete tracers, the moment-generating function of the galaxy count frequency distribution is expressed through the void probability function.

2.4. Cumulant bias parameters

We may further consider how the bias parameters are related to the cumulants of the galaxy environment distribution. In probability theory, the value of cumulants are generally thought to characterise a distribution more independently than its moments. For example, if a distribution has a large first moment ⟨x⟩, then we should also expect that it has large second and third moments ⟨x²⟩ and ⟨x³⟩. However, the first cumulant of a distribution ⟨x⟩ says very little about the second and third cumulants, ⟨(x−⟨x⟩)²⟩ and ⟨(x−⟨x⟩)³⟩, respectively. For example, for a Gaussian distribution the mean is independent of its variance and the third and higher order cumulants are actually zero.

The cumulant generating function is defined as

$$\begin{array}{cc}\hfill K_g(t)& =log{M}_g(t)\hfill \\ \hfill & =\frac{1}{2}{t}^2{\sigma }^2+logF(t{\sigma }^2).\hfill \end{array}$$(23)

Cumulants of the galaxy environment distribution are simply given by the derivatives of the cumulant generating function

$$\begin{array}{cc}\hfill {\kappa }_n& ={\frac{\partial K_g}{\partial t^n}|}_{t=0}.\hfill \end{array}$$(24)

Therefore, they can be evaluated as:

$$\begin{array}{cc}\hfill {\kappa }_n& =\{\begin{array}{cc}{\beta }_n{\sigma }^{2n}\hfill & ifn\ne 2,\hfill \\ {\beta }_2{\sigma }^4+{\sigma }^2\hfill & ifn=2,\hfill \end{array}\hfill \end{array}$$(25)

where we have defined

$$\begin{array}{cc}\hfill {\beta }_n& ={\frac{{\partial }^n}{\partial {\delta }_0^n}logF({\delta }_0)|}_{{\delta }_0=0}.\hfill \end{array}$$(26)

We notice that the question of whether β₂ is above or below zero is a direct indication of whether the variance of the halo environment distribution is larger or smaller than that of the background. We refer to these parameters as ‘cumulant bias parameters’ and they are directly related to the canonical bias parameters b_n. Comparing their definition to the canonical bias parameters in Equation (3) shows that they relate to each other exactly in the same way that cumulants relate to moments (compare Equations (22) and (24)). For example, at the first four orders we have:

$$\begin{array}{cc}\hfill {\beta }_1& =b_1,\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill {\beta }_2& =b_2-b_1^2,\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill {\beta }_3& =b_3-3b_1{b}_2+2b_1^3,\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill {\beta }_4& =b_4-4b_1{b}_3-3b_2^2+12b_1^2{b}_2-6b_1^4.\hfill \end{array}$$(30)

2.5. Interpretation

Given a set of bias parameters, these relations allow us to directly find the cumulants of the galaxy environment distribution, or alternatively, they allow us to infer canonical bias parameters by measuring cumulants of the galaxy environment distribution.

However, it is also possible to phrase the bias expansion directly in terms of the cumulant bias parameters.

$$\begin{array}{cc}\hfill logF& ={\beta }_1{\delta }_0+\frac{1}{2}{\beta }_2{\delta }_0^2+\frac{1}{6}{\beta }_3{\delta }_0^3+\cdots \hfill \end{array}$$(31)

There are several reasons to believe that they may form a better set of parameters than the canonical bias parameters:

• The cumulant bias parameters are the derivatives of the logarithm of the galaxy density. Therefore, they presuppose the positivity of F and may be better behaved, especially in low density regions.

• The canonical bias parameters behave similarly to moments of the galaxy environment distribution. For example, we may expect that if b₁ is large, automatically b₄ will also be large. On the other hand, we may expect β₄ to be independent of the value of β₁ – just as cumulants are relatively independent of each other.

• If the galaxy environment distribution has the form of a Gaussian then β_n = 0 for all n > 2. On the other hand, all b_n would be non-zero in this case. Therefore, if all β_n beyond a degreeof 2 are small, this motivates the usage of a Gaussian bias model – where both log F and log f are quadratic polynomials. We show here that this is indeed the case for halos.

Furthermore, we note that the cumulant bias parameters and the probabilistic considerations in this section motivate novel approaches to parameterising the bias function at finite smoothing scales. There are several possible well motivated approaches and we only mention them here in brief, leaving a more thorough investigation to future studies.

The first and most straight-forward proposition is to assume an expansion of the logarithm of the bias function, leading to an exponential of a polynomial

$$\begin{array}{cc}\hfill logf(\delta )& ={\gamma }_0+{\gamma }_1\delta +{\gamma }_2{\delta }^2+{\gamma }_3{\delta }^3+\cdots \hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill f(\delta )& =exp({\gamma }_0+{\gamma }_1\delta +{\gamma }_2{\delta }^2+{\gamma }_3{\delta }^3+\cdots ).\hfill \end{array}$$(33)

Such a bias expansion would have to be truncated at an even order (if n ≥ 2) and the highest order coefficient must be negative to guarantee a fully normalised probability distribution. By definition, this guarantees the positivity of the bias function which is arguably a desirable property – especially if our intention is to create mocks from a bias model. Furthermore, the resulting distribution would be part of the exponential family which guarantees several desirable properties. The main limitation is that at orders of n > 2, it might be difficult to properly re-normalise the expression analytically⁴.

This choice is the only general one we know of that guarantees f(δ) > 0 for all values of δ, even if the expansion is truncated. However, there are a few other noteworthy options that do not fulfill this constraint. One option is motivated by the intriguing simplicity of the relation between cumulant generating the function in Equation (23) and F. This motivates to define the bias function directly through the cumulant generating function,

$$\begin{array}{cc}\hfill K_g(t)& =K_{\mathrm{m}}(t)+\sum _n\frac{1}{n!}{\beta }_{n}t{\sigma }^{2n,}\hfill \end{array}$$(34)

where K_m(t) is the cumulant generating function of the matter distribution. This general Ansatz might also be feasible in Eulerian space, but for the Lagrangian case, it is simply $K_{\mathrm{m}}(t)=\frac{1}{2}t{\sigma }^2$.

If the expansion is truncated at an order of n = 2 then this corresponds to a well defined Gaussian bias model, but if it is truncated at a higher order, then it is not easy to find an analytical form of f(δ). However, a numerical expression can be obtained through a Fourier transformation. Unfortunately, such a truncation at higher orders does not lead to a well defined probability density, but instead may have negative function values. In fact the cumulant generating function does not yield a positive probability density function (pdf) for any finite order polynomial of degree n > 2 (Lukacs 1970). However, this does not disfavour this approach over the canonical bias expansion which also yields a negative pdf.

Furthermore, we note that it is an option to phrase the bias expansion relative to a Gaussian distribution that has the correct cumulants up to an order of 2, as follows:

$$\begin{array}{cc}\hfill p(\delta |g)& =\frac{1}{\sqrt{2\pi }{\sigma }_g}exp(-\frac{{(\delta -{\mu }_g)}^2}{2{\sigma }_g^2})(1+{\gamma }_3{\delta }^3+{\gamma }_4{\delta }^4+\cdots ),\hfill \end{array}$$(35)

where μ_g = κ₁ = β₁σ² and σ_g² = κ₂ = β₂σ⁴ + σ². This is also where the γ_n start at an order of 3, since the lowest two orders are already specified through the Gaussian. This type of expansion is similar to a truncated Edgeworth series (which is an expansion to approximate distributions with given cumulants).

Finally, we note that it is possible to use the cumulant bias parameters in a classical polynomial bias expansion as well. For example, a third-order expansion that only uses two parameters, assuming β₃ = 0 would read

$$\begin{array}{cc}\hfill f(\delta )& ={\beta }_1\delta +({\beta }_2+{\beta }_1^2)({\delta }^2-{\sigma }^2)+(3{\beta }_1{\beta }_2-{\beta }_1^3)({\delta }^3-3\delta {\sigma }^2)+\cdots \hfill \end{array}$$(36)

This can be understood as the lower order terms predicting the likely behaviour of higher order terms.

2.6. Multivariate estimators

So far, we have only considered the bias function based on the assumption that the only known aspect of the environment is the density, δ, at the considered location. However, other aspects of the environment may be known, such as the Laplacian,

$$\begin{array}{cc}\hfill L& ={\mathrm{\nabla }}^2\delta ,\hfill \end{array}$$(37)

or higher order derivatives of the density field,

$$\begin{array}{cc}\hfill P& ={\mathrm{\nabla }}^4\delta ,\hfill \end{array}$$(38)

or the tidal field. Since the tidal field is inherently a tensorial quantity, its treatment is slightly more complicated and we therefore discuss this in more detail in Sect. 4.

Assuming that we are dealing only with scalar quantities that describe the environment, we may summarise the environment through a single vector x – for example, if we consider the density and Laplacian as variables we have x = (δ, L)^T. The majority of the relations that we have derived in the previous section hold in a similar manner for the multivariate case as well. We do not derive all of them again, since their derivation follows almost completely analogously, however, we do list all the important relations.

Assuming that all the considered environment variables follow from linear operators on the Gaussian random field, their distribution is given by a multivariate Gaussian,

$$\begin{array}{cc}\hfill p(\mathit{x})& =\frac{1}{2\pi \sqrt{det(\mathsf{C})}}exp(-\frac{1}{2}{\mathit{x}}^T{\mathsf{C}}^{-1}\mathit{x}).\hfill \end{array}$$(39)

For example, for the case of x = (δ, L)^T, which is particularly well motivated by the peak model (Bardeen et al. 1986), we have the covariance matrix:

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\delta ,L}& =\left[\begin{array}{ccc}{\sigma }_0^2& -{\sigma }_1^2{\sigma }_1^2& {\sigma }_2^2\end{array}\right],\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\delta ,L}^{-1}& =\left[\begin{array}{cc}\frac{{\sigma }_2^2}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}& \frac{{\sigma }_1^2}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}\\ \frac{{\sigma }_1^2}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}& \frac{{\sigma }_0^2}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}\end{array}\right],\hfill \end{array}$$(41)

where σ₀² = ⟨δ²⟩, σ₂² = ⟨L²⟩ and σ₁² = ⟨−δL⟩.

Now, we can conveniently write bias parameters in a matrix form

$$\begin{array}{cc}\hfill {\mathit{b}}_1& ={{\mathrm{\nabla }}_{\mathit{x}}F|}_{\mathit{x}=0},\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill {\mathsf{b}}_2& ={({\mathrm{\nabla }}_{\mathit{x}}\otimes {\mathrm{\nabla }}_{\mathit{x}})F|}_{\mathit{x}=0},\hfill \end{array}$$(43)

$$\begin{array}{cc}\hfill {\mathsf{b}}_3& ={({\mathrm{\nabla }}_{\mathit{x}}\otimes {\mathrm{\nabla }}_{\mathit{x}}\otimes {\mathrm{\nabla }}_{\mathit{x}})F|}_{\mathit{x}=0},\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill {\mathsf{b}}_n& ={{\mathrm{\nabla }}_{\mathit{x}}^{\otimes n}F|}_{\mathit{x}=0}.\hfill \end{array}$$(45)

Here, ⊗ designates an outer product, the power notation in the last line designates a repeated outer product, b₁ is a vector, b₂ is a symmetric rank two matrix, b₃ is a symmetric rank three tensor, and so on. Furthermore, ∇_x denotes a gradient with respect to the chosen variables, for example, ∇_x = (∂/∂δ, ∂/∂L)^T.

It is then straightforward to show

$$\begin{array}{cc}\hfill {\mathsf{b}}_n& ={(-1)}^n{\langle \frac{{\mathrm{\nabla }}_{\mathit{x}}^{\otimes n}p}{p}\rangle }_{\mathrm{g}},\hfill \end{array}$$(46)

which at the first two orders, can be expressed as:

$$\begin{array}{cc}\hfill {\mathit{b}}_1& ={\langle {\mathsf{C}}^{-1}\mathit{x}\rangle }_{\mathrm{g}},\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill {\mathsf{b}}_2& ={\langle {\mathsf{C}}^{-1}(\mathit{x}\otimes \mathit{x}-\mathsf{C}){\mathsf{C}}^{-1}\rangle }_{\mathrm{g}}.\hfill \end{array}$$(48)

Importantly, these expressions do not lead to the same bias estimators for the density if a second variable is present that is correlated with the density. For example, if we use the density and Laplacian as variables, we find

$$\begin{array}{cc}\hfill b_{\delta ,\mathrm{so}2}& ={\langle \frac{\delta {\sigma }_2^2+L{\sigma }_1^2}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}\rangle }_{\mathrm{g}},\hfill \end{array}$$(49)

$$\begin{array}{cc}\hfill b_{L,\mathrm{so}2}& ={\langle \frac{L{\sigma }_0^2+\delta {\sigma }_1^2}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}\rangle }_{\mathrm{g}},\hfill \end{array}$$(50)

for the components of b₁ and

$$\begin{array}{cc}\hfill b_{\delta \delta ,\mathrm{so}2}& ={\langle \frac{{(\delta {\sigma }_2^2+L{\sigma }_1^2)}^2-{\sigma }_2^2({\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4)}{{({\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4)}^2}\rangle }_{\mathrm{g}},\hfill \end{array}$$(51)

$$\begin{array}{cc}\hfill b_{LL,\mathrm{so}2}& ={\langle \frac{{\sigma }_0^2(-{\sigma }_0^2{\sigma }_2^2+{\sigma }_1^4)+{(L{\sigma }_0^2+\delta {\sigma }_1^2)}^2}{{({\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4)}^2}\rangle }_{\mathrm{g}},\hfill \end{array}$$(52)

$$\begin{array}{cc}\hfill b_{\delta L,\mathrm{so}2}& ={\langle \frac{(L{\sigma }_0^2+\delta {\sigma }_1^2)(L{\sigma }_1^2+\delta {\sigma }_2^2)-{\sigma }_1^2({\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4)}{{({\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4)}^2}\rangle }_{\mathrm{g}},\hfill \end{array}$$(53)

for the components of b₂, where we have written, for instance, b_δδ as a symbol for what was previously referred to as b₂. Furthermore, we note that a general form for the Nth density bias parameter is given by

$$\begin{array}{cc}\hfill b_{{\delta }^N,\mathrm{so}2}& ={\langle \frac{H_N((\delta +L{\sigma }_1^2/{\sigma }_2^2)/{\sigma }_{\ast })}{{\sigma }_{\ast }^N}\rangle }_{\mathrm{g}},\hfill \end{array}$$(54)

with ${\sigma }_{\ast }^2={\sigma }_0^2-{\sigma }_1^4/{\sigma }_2^2$. We note that these relations have already been derived through the response of density peaks to large scale perturbations (Bardeen et al. 1986; Mo & White 1996; Desjacques et al. 2010) and were later been shown to apply also to more general tracers (Lazeyras et al. 2016).

The difference between Eqs. (49)–(54) and the estimators from Eqs. (15)–(18) is that they correspond to partial derivatives of the density at fixed values of the Laplacian. The previous estimators are the derivatives of the projected distribution, which is not the same due to the correlation between density and Laplacian. As we see in Sect. 3.3, the new estimators are less scale-dependent, since they are closer to pure partial derivatives with respect to the density (in the separate-universe sense).

We refer to the estimators from Equations (49)–(54) as the estimators with spatial corrections of the second order, whereas the previous estimators are without any spatial corrections (i.e. at the spatial order of 0). It is straightforward to obtain higher order spatial corrections, for instance, by considering the covariance matrix of the distribution of x = (δ, L, P)^T with P as in Equation (38). For example, the estimators of density bias parameters at a spatial order of 4 can be expressed as

$$\begin{array}{cc}\hfill b_{{\delta }^N}& ={\langle \frac{H_N({\delta }_{\ast 4}/{\sigma }_{\ast 4})}{{\sigma }_{\ast 4}^N}\rangle }_{\mathrm{g}},\hfill \\ \hfill {\delta }_{\ast 4}& =\delta +L\frac{({\sigma }_1^2{\sigma }_4^2-{\sigma }_2^2{\sigma }_3^2)}{{\sigma }_2^2{\sigma }_4^2-{\sigma }_3^4}+P\frac{({\sigma }_1^2{\sigma }_3^2-{\sigma }_2^4)}{{\sigma }_2^2{\sigma }_4^2-{\sigma }_3^4},\hfill \\ \hfill {\sigma }_{\ast 4}& =\frac{{\sigma }_0^2{\sigma }_2^2{\sigma }_4^2-{\sigma }_0^2{\sigma }_3^4-{\sigma }_1^4{\sigma }_4^2+2{\sigma }_1^2{\sigma }_2^2{\sigma }_3^2-{\sigma }_2^6}{{\sigma }_2^2{\sigma }_4^2-{\sigma }_3^4},\hfill \end{array}$$(55)

where σ₃² = −⟨δ⋅P⟩ and σ₄² = ⟨P²⟩. In practice, it is easier to evaluate such high-order estimators numerically, rather than deriving explicit expressions for them. The covariance matrix may easily be measured from a given linear density field and then used to evaluate the estimators as e.g. in Equation (48), so that in principle it is not very difficult to obtain corrections at any spatial order. Throughout this article, we mainly focus on the estimators of a spatial order of 2, since they are already sufficiently accurate to obtain good measurements of biases, but we have selectively included both lower or higher order estimators throughout this article to demonstrate the convergence.

We note that the estimators with second-order spatial corrections are quite similar in spirit to the considerations by Musso et al. (2012), Paranjape & Sheth (2012), Paranjape et al. (2013b,a) to map scale-dependent measurements of biases as in Equation (14) under the assumption of excursion set (peaks) models to scale-independent large-scale parameters. However, here we can see that such a mapping to the large-scale limit can be done without the assumption of Peak or excursion set models and for spatial corrections of any order. Consider for example that the relation between the zeroth-order spatial estimator and second-order spatial estimator of b₁ is given by

$$\begin{array}{cc}\hfill {\langle \delta \rangle }_{\mathrm{g}}=\mathsf{C}{\mathit{b}}_1& ={\sigma }_0^2{b}_{1,\mathrm{so}0}\hfill \\ \hfill & ={\sigma }_0^2{b}_{1,\mathrm{so}2}-{\sigma }_1^2{b}_{L,\mathrm{so}2}.\hfill \end{array}$$

Therefore, b_L, so2 predicts the scale dependence of b_1, so0 and b_1, so2 posits a more scale-independent estimate:

$$\begin{array}{cc}\hfill b_{1,\mathrm{so}0}& =b_{1,\mathrm{so}2}-\frac{{\sigma }_1^2}{{\sigma }_0^2}{b}_{L,\mathrm{so}2}.\hfill \end{array}$$(56)

We note that this argument can easily be pushed to higher spatial orders:

$$\begin{array}{cc}\hfill b_{1,\mathrm{so}0}& =b_{1,\mathrm{so}4}-\frac{{\sigma }_1^2}{{\sigma }_0^2}{b}_{L,\mathrm{so}4}+\frac{{\sigma }_2^2}{{\sigma }_0^2}{b}_{P,\mathrm{so}4}.\hfill \end{array}$$(57)

The additional cost of our approach (e.g. in comparison to Paranjape et al. 2013a) is that, beneath the density, also the Laplacian has to be evaluated in Lagrangian space, but the benefits are model-independence and a much reduced mathematical complexity.

2.7. Multivariate cumulants

Analogously to the monovariate case, we define multivariate cumulant bias parameters as

$$\begin{array}{c}\hfill {\mathit{\beta }}_N={{\mathrm{\nabla }}_{\mathit{x}}^{N}log(F)|}_{\mathit{x}=0},\end{array}$$(58)

which immediately leads to

$$\begin{array}{cc}\hfill {\mathit{\beta }}_1& ={\mathit{b}}_1,\hfill \end{array}$$(59)

$$\begin{array}{cc}\hfill {\mathit{\beta }}_2& ={\mathsf{b}}_2-{\mathit{b}}_1\otimes {\mathit{b}}_1,\hfill \end{array}$$(60)

$$\begin{array}{cc}\hfill {\beta }_{3,ijk}& =b_{3,ijk}-(b_{1,i}{b}_{2,jk}+b_{1,j}{b}_{2,ki}+b_{1,k}{b}_{2,ij})+2b_{1,i}{b}_{1,j}{b}_{1,k},\hfill \end{array}$$(61)

where we gave β₃ in index notation, since the central term is difficult to express in vectorial notation. We note that this leads to the same relations between the density cumulant biases and their canonical biases as in Equations (27)–(30), but it also includes additional relations, such as:

$$\begin{array}{cc}\hfill {\beta }_{\delta \delta L}& =b_{\delta \delta L}-b_{\delta \delta }{b}_L-2b_{\delta L}{b}_{\delta }+2b_{\delta }^2{b}_L.\hfill \end{array}$$(62)

Finally, the relation between the separate-universe bias function and the moment and cumulant generating functions is analagous to the expression in Equations (21) and (23), given by

$$\begin{array}{cc}\hfill M(\mathit{t})& =exp\left(\frac{1}{2}{\mathit{t}}^T\mathsf{C}\mathit{t}\right)F({\mathit{x}}_0=\mathsf{C}\mathit{t}),\hfill \end{array}$$(63)

$$\begin{array}{cc}\hfill K(\mathit{t})& =\frac{1}{2}{\mathit{t}}^T\mathsf{C}\mathit{t}+log(F(\mathsf{C}\mathit{t})),\hfill \end{array}$$(64)

which can be differentiated to show that the cumulants are given in an easy form through the cumulant bias parameters, for instance,

$$\begin{array}{cc}\hfill {\mathit{\kappa }}_1& =\mathsf{C}{\mathit{\beta }}_1,\hfill \end{array}$$(65)

$$\begin{array}{cc}\hfill {\mathit{\kappa }}_2& =\mathsf{C}{\mathit{\beta }}_2\mathsf{C}+\mathsf{C},\hfill \end{array}$$(66)

$$\begin{array}{cc}\hfill {\kappa }_{3,ijk}& =\sum _{\mathit{abc}}{C}_{\mathit{ai}}{C}_{\mathit{bj}}{C}_{\mathit{ck}}{\beta }_{\mathit{abc}}.\hfill \end{array}$$(67)

However, instead of measuring directly the cumulants of the galaxy environment distribution and inverting Equations (65)–(67), it is in practice simpler to instead define the variable,

$$\begin{array}{cc}\hfill \mathit{u}& ={\mathsf{C}}^{-1}\mathit{x}.\hfill \end{array}$$(68)

We then measure its cumulants, κ_u, which relate to the cumulant biases as

$$\begin{array}{cc}\hfill {\mathit{\beta }}_{ijk\dots }& =\{\begin{array}{cc}{\kappa }_{u,ijk\dots }\hfill & ifi+j+k+\cdots \ne 2,\hfill \\ {\kappa }_{u,ijk\dots }-C_{\mathit{ab}}^{-1}\hfill & ifi+j+k+\cdots =2,\hfill \end{array}\hfill \end{array}$$(69)

where a and b indicate the indices of the non-zero variables in the second-order case.

We show in Appendix A that it is also possible to derive similar estimators if the density field was filtered with a function different than a sharp k-space filter. The only additional complication in that case that it is necessary to account for the correlation matrix between the smoothed large scales and unsmoothed small scales. However, for simplicity, we focus only on measurements with the sharp k-space filter in this work.

3. Density bias measurements

In this section, we explain our evaluation of the Lagrangian density bias parameters for different sets of halos. We want to verify the consistency of our estimators by comparing them to the vast literature available on the subject. Furthermore, we also compare the canonical bias parameters, b_n, with the cumulant bias parameters, β_n.

3.1. Simulation

For the analysis through this paper, we used a single cosmological box simulation with high resolution. This simulation is part of the ‘BACCO simulation project’ that was first introduced in Angulo et al. (2021). It has a box size of L = 1440h⁻¹Mpc with 4320³ particles leading to a mass resolution of m_p = 3.2 × 10⁹ h⁻¹ M_⊙. The cosmological parameters are Ω_m = 0.307, Ω_Λ = 0.693, Ω_b = 0.048., n_s = 0.9611, σ₈ = 0.9, h = 0.677 which are similar to the Planck Collaboration VI (2020) cosmology except for the roughly 10% larger value of σ₈. By default we use this simulation at a scale-factor of a = 1.08, corresponding approximately to a = 1.

To identify halos, the simulation code uses a modified version of SUBFIND (Springel et al. 2001), which first identifies halos via a friends of friends (FoF) algorithm and subsequently calculates for each FoF group, the mass, M_200b, in a region that encloses 200 times the mean density of the Universe.

3.2. Bias measurements

We defined sets of biased tracers by considering halos selected by their mass. For this, we considered 21 equally log-spaced halo masses between 10¹² h⁻¹ M_⊙ and 10¹⁶ h⁻¹ M_⊙ and for each halo mass, we considered all halos with masses in the range of

$$\begin{array}{c}\hfill M_{200b}\in [M_i/1.25,M_i·1.25].\end{array}$$(70)

To evaluate our Lagrangian bias estimators, we only need to know the linear field evaluated at the Lagrangian locations of the tracers. We approximated the Lagrangian location of each halo through the Lagrangian coordinate of its most bound particle. Since the simulation started from a Lagrangian grid, the Lagrangian origin of the most bound particle can easily be inferred from its id i_mb as

$$\begin{array}{cc}\hfill {\mathit{q}}_{\mathrm{mb}}& =\frac{L}{N_{\mathrm{grid}}}\left(\begin{array}{c}{i}_x\\ i_y\\ i_z\end{array}\right),\hfill \\ \hfill i_{\mathrm{mb}}& =i_x{N}_{\mathrm{grid}}^2+i_y{N}_{\mathrm{grid}}+i_z,\hfill \end{array}$$(71)

where N_grid = 4320 is the number of particles per dimension.

We know the linear density field of the simulation through the initial condition generator. To save on the computation time, we can create a low resolution grid representation of the linear density field with N_lin³ grid points. For fields different than the linear density field, we additionally multiply by the correct operator (e.g. −k² for the Laplacian) in Fourier space. We created a smoothed version of this field by multiplying with a sharp k−filter in Fourier space

$$\begin{array}{cc}\hfill {\delta }_k& ={\delta }_{\mathrm{lin},\mathrm{k}}·\mathrm{\Theta }(k_{\mathrm{d}}-k),\hfill \end{array}$$(72)

where Θ is the heavy-side function and we tested, for each measurement, different damping scales k_d ∈ [0.1, 0.15, 0.2, 0.25, 0.3] h⁻¹ Mpc. We then deconvolved this field with a linear interpolation kernel and interpolated it to the Lagrangian locations of our tracer set. We chose an N_lin value that is sufficiently large that the resulting interpolated values are virtually independent of this discretisation; for instance, N_lin = 183 at k_d = 0.1 h⁻¹ Mpc and N_lin = 549 at k_d = 0.3 h⁻¹ Mpc.

With the linear densities at the Lagrangian locations of tracers, it is easy to evaluate the zeroth-order spatial estimators of the biases by averaging as in Equation (14). However, for higher spatial order estimators we also need to know the values of the Laplacian (Eqs. (49)–(53)) or that of the fourth derivative (Equation (55)). The Laplacian is inferred by multiplying the linear density field in Fourier space by −k² and then interpolated to the tracer locations. Other variables such as the fourth derivatives can be evaluated in a similar manner. This allows us to evaluate any of the bias estimators (from Sect. 2) directly as simple expectation values over tracers. We note that this is a major difference compared to the measurement process in Paranjape et al. (2013a), where the zeroth-order spatial estimators, as those in Equation (14), were evaluated, but then mapped to scale-independent parameters through peak theory arguments.

To estimate the covariance of a set of measured bias parameters we use a Jackknife technique. For this we divide the box in Lagrangian space into N_jk³ subboxes with N_jk = 4. We perform 64 measurements of the vector of bias parameters b_i by subsequently leaving out all tracers in one of the subboxes. Then we estimate the covariance through

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathit{b}}& =\frac{1}{N_{\mathrm{jk}}^3-1}\sum _i({\mathit{b}}_i-{\mathit{b}}_0)\otimes ({\mathit{b}}_i-{\mathit{b}}_0),\hfill \end{array}$$(73)

$$\begin{array}{cc}\hfill {\mathit{b}}_0& =\frac{1}{N_{\mathrm{jk}}^3}\sum _i{\mathit{b}}_i.\hfill \end{array}$$(74)

It is worth noting that the Jackknife estimator results in a more reliable estimate of the uncertainty of the measurement than, for instance, a simple bootstrap would yield. When comparing these estimators we found that the Jackknife gives larger (more conservative) error estimates, since it also accounts for the uncertainty induced by cosmic variance, which is sampled by leaving out spatially correlated parts of the data sets.

3.3. b₁ and the scale dependence of estimators

In Figure 2, we show the measurements of b₁ for halos as a function of M_200b using the estimator from Equation (15) without spatial corrections (top panel) versus the estimator from Equation (49) with spatial corrections at the second order (bottom panel). The dashed lines indicate the fitting function from Tinker et al. (2010) evaluated for the cosmology of our simulation.

As expected, halos of large masses M ∼ 10¹⁵ h⁻¹ M_⊙ are highly biased b₁ ∼ 3, whereas low-mass halos (M ∼ 10¹² h⁻¹ M_⊙) have an even slightly negative (Lagrangian) bias. The zeroth-order spatial estimators have a very small degree of statistical uncertainty, but exhibit a significant scale dependence and are inconsistent across different damping scales. On the other hand, the estimators with a spatial order of 2 have a larger statistical uncertainty, but seem consistent across different damping scales – except at very large k_d and large halo masses; in this case, some scale dependence becomes visible, probably indicating the effect of higher spatial-order terms.

Our second-order spatial measurements agree at all scales well with the fitting function from Tinker et al. (2010). While they lie systematically about 5% below this fit, this is within the quoted uncertainty of the Tinker et al. (2010) fit and a similar difference can be seen in Lazeyras et al. (2016). In comparison to other studies at z = 0 (e.g. Lazeyras et al. 2016), we have lower bias values at the same masses (e.g. b₁ ∼ 0.5 at 10¹⁴ h⁻¹ M_⊙ instead of b₁ ∼ 1) because of the later time (z = −0.08) and relatively high variance (σ₈ = 0.9) of our simulation.

To further highlight the difference in the scale dependence of the estimators, we show in Figure 3 the error Δb₁ in b₁ measurements as a function of damping scale for a few different halo masses. Here, we phrase the error Δb₁ = b₁ − b_1, ref relative to the second-order spatial estimator at k_d = 0.2 h Mpc⁻¹ as b_1, ref. As expected, the zeroth-order estimator seems to approach the selected reference value on large scales. However, the zeroth-order second-order spatial estimator exhibits a significant scale dependence, whereas the second-order spatial estimators are almost scale-independent, except for the largest one at k_d ≳ 0.25 h Mpc⁻¹. We also show measurements with the fourth-order spatial estimator from Equation (55), which are scale-independent to even smaller scales, showing that the measurements converge well with adding higher order spatial corrections. However, for getting a good estimate of the bias parameters it seems sufficient to use the second-order spatial estimator at scales of k_d ≲ 0.2 h Mpc⁻¹, which bears a similar accuracy in measurement, for instance, akin to evaluating the fourth-order spatial estimator at k_d ∼ 0.3 h Mpc⁻¹.

We can estimate the scale dependence of the zeroth-order spatial estimator, as shown in Equations (56) and (57). We mark the estimated scale dependence of the zeroth-order spatial estimators using the second-order spatial bias parameters measured at k_d = 0.2 h Mpc⁻¹ as dashed lines in Figure 3 and using the fourth-order spatial parameters from k_d = 0.3 h Mpc⁻¹ as dotted lines. The scale dependence of the zeroth-order spatial estimator is predicted quite well, showing that the consideration of the Laplacian introduces a correction of the order of σ₁²/σ₀², which scales approximately as k_d² and including of the fourth derivatives provides additional corrections that become relevant on even smaller scales.

Our measurements therefore confirm the results from previous studies that including the Laplacian is vital to recover scale-independent bias on large scales (Paranjape et al. 2013a; Lazeyras et al. 2016). Furthermore, if the modelling is pushed to small scales k_d ≳ 0.25h the inclusion of additional higher spatial derivative terms may prove beneficial.

3.4. b₂ versus β₂

The first cumulant bias parameter that is different from the canonical bias parameter of the same order is β₂ = b₂ − b₁² where b₁ = β₁. Since there is a one-to-one relation between (b₁, b₂) and (β₁, β₂) it might seem that there should not be any advantage by using β₂ over b₂ as a parameter when fitting datasets. However, here we show that β₂ seems to be more independent of β₁ than b₂ of b₁, especially when their covariance matrix is considered.

In Figure 4, we show the measured co-evolution relation between b₂ and b₁ (left) versus the co-evolution between β₂ and β₁ (right) for the second-order spatial estimators. In comparison, we also show the co-evolution relations measured by Lazeyras et al. (2016) as a dashed line, which seems to match our measured relation well up to k_d ≲ 0.15 − 0.2 h Mpc⁻¹, showing that our method for measuring the bias parameters is indeed able to recover the correct large-scale limit of the bias parameters.

We note that for halos we find β₂ < 0 across all masses and with high statistical significance. This means that the width of the halo environment distribution p(δ|g) is smaller than the width of the background distribution p(δ) – showing that halo formation is more selective than a random distribution. Possibly, the assumption β₂ < 0 could be used to limit the considered prior range when fitting certain galaxy surveys. However, it is difficult to anticipate whether this should hold for every possible set of galaxies.

Furthermore, we note that the co-evolution relation between β₂ and β₁ is monotonic and roughly linear. For b₁ ≳ 2, it satisfies |β₂|< |b₂|. Therefore, β₂ versus β₁ seems slightly simpler than the b₂ to b₁ relation.

The benefit of using β₂ as a parameter is even more evident when considering the covariance matrix of the measurements. In Figure 5 we show the uncertainty of the b₂ measurements (solid lines), versus the uncertainty of the β₂ measurements. Comparing different damping scales shows that the inclusion of smaller scales in the measurements decreases the statistical uncertainty (but increases systematic error). In all cases and especially at high values of b₁ the uncertainty of β₂ is significantly smaller than that of b₂. This is so, since b₂ and b₁ are more correlated than β₂ and β₁.

To highlight this we show in Figure 6 the correlation coefficient

$$\begin{array}{cc}\hfill {\alpha }_{12}& :=\frac{C_{12}}{\sqrt{C_{11}{C}_{22}}},\hfill \end{array}$$(75)

where C is the covariance matrix between b₁ and b₂ or β₁ and β₂. For b₁ ≫ 0 the correlation coefficient of (b₁, b₂) is quite large, even close to 1 in some cases. On the other hand, the correlation coefficient of (β₁, β₂) is quite small and seems on average consistent with 0. Therefore, measurements of b₁ are quite entangled with b₂, but not so much with β₂. While we do not show it here, we have found that this is also the case for higher order correlations, for instance, those of (b₁, b₃) versus (β₁, β₃).

This can be understood when considering the difference between moments and cumulants of a probability distribution. For instance, if we knew with respect to a certain distribution p(x) that it has a large first moment ⟨x⟩, then we might also expect that the second moment ⟨x²⟩ is large. On the other hand, knowledge about the mean of a distribution, is very uninformative about its width – that is its second cumulant σ_x² = ⟨(x−⟨x⟩)²⟩ = ⟨x²⟩ − ⟨x⟩². In the same sense, we would expect the cumulant bias parameters to be more independent of each other than canonical bias parameters.

3.5. Higher order biases

While the second-order cumulant bias β₂ has already slight advantages over b₂ as highlighted in the last section, the benefits are even more significant at higher orders. Here, we compare the behaviour of b₃ and b₄ with the behaviour of β₃ and β₄.

In Figure 7, we show b₃, b₄, β₃, and β₄ as a function of b₁ for the second-order spatial estimator. The uncertainty of the measurements of b₃ and b₄ are much larger than for b₂ so that we hardly find any signal at k_d ∼ 0.1 h Mpc⁻¹. However, at k_d ∼ 0.15 h Mpc⁻¹ we can find a meaningful signal and it is consistent with the b₃ measurements of Lazeyras et al. (2016).

Strikingly, the value of β₃ is extremely small in comparison to b₃ and seems approximately consistent with 0. It is noteworthy that this is not only the case for scales where b₃ is reasonably scale independent (k_d ≲ 0.15 h Mpc⁻¹), but it is also so for much smaller scales. There does not seem to be any noteworthy relation between β₃ and β₁. Furthermore, the co-evolution relation measured by Lazeyras et al. (2016) seems consistent with β₃ = 0. Therefore, we could summarise the third-order co-evolution relation through the relation we get by using β₃ = 0 in Equation (29):

$$\begin{array}{cc}\hfill b_3& =3b_1{b}_2-2b_1^3.\hfill \end{array}$$(76)

Similarly, we find for β₄ an apparent independence of β₁ and consistency with zero across different scales. Assuming β₄ = 0 in Equation (30) leads to the co-evolution relation between b₄, b₃, b₂, and b₁ given by

$$\begin{array}{cc}\hfill b_4& =4b_1{b}_3+3b_2^2-12b_1^2{b}_2-6b_1^4.\hfill \end{array}$$(77)

We mark this relation with the measurements of b₁, b₂, and b₃ from Lazeyras et al. (2016) as a dashed line in the bottom left panel of Figure 7. This seems in good agreement with the measured b₄.

We conclude that at an order of n ≥ 3, the cumulant bias parameters are very close to 0. The co-evolution relations at these orders can be summarised simply through β_n = 0. We might therefore suggest that high-order co-evolution relations do not represent any physical insights, but rather highlight that the canonical bias parameters form a suboptimal, deeply entangled basis.

Further, the fact that all the cumulant biases beyond an order of 2 are close to zero indicates that the galaxy environment distribution is very well approximated by a Gaussian distribution. We plan to investigate the possibility of a Gaussian bias model in Stücker et al. (2025), where we also explain how this fact may arise naturally from the Gaussianity of the background distribution.

3.6. Laplacian bias

In Figure 8, we show the Lagrangian Laplacian bias parameter as a function of halo mass, as inferred from the estimator in Equation (50). The Laplacian measurements have a much stronger scale dependence than our previous measurements of density bias parameters. However, it seems that at masses M ≲ 3 × 10¹³ h⁻¹ M_⊙ the scale dependence disappears and a reliable measurement is obtained.

We compare our Lagrangian measurements with the Eulerian measurements of Lazeyras & Schmidt (2019). In general, the relation between the Lagrangian and Eulerian Laplacian bias depends on velocity bias (e.g. Desjacques et al. 2018). The velocity bias b_s quantifies the difference between the displacement field of matter s and galaxies s_g and gets as first-order contributions of the form

$$\begin{array}{c}\hfill {\mathit{s}}_{\mathrm{g}}-\mathit{s}\approx b_s\mathrm{\nabla }\delta .\end{array}$$(78)

Therefore, the galaxy density can get an extra contribution proportional to the divergence of this field, that scales as L = ∇²δ. However, for our set of tracers s_g = s by definition, so that b_s = 0 and we simply assume that their Lagrangian and Eulerian Laplacian biases are identical or at least very similar.

In the reliable range of M ≲ 3 × 10¹³ h⁻¹ M_⊙, it seems that our measurements agree well with the linear and the quadratic (in halo radius) fit that was inferred by Lazeyras & Schmidt (2019) for the Eulerian Laplacian bias. Therefore, we confirm that the Laplacian bias is negative for halos and that it can plausibly have amplitudes of the order of b_L ∼ −20 Mpc² h⁻². However, we leave more detailed considerations to future studies. In principle, more accurate measurements could be obtained by consideration of higher order corrections and through measurements at higher redshifts.

3.7. The scale-split break-down scale

The core assumption of the bias formalism is the separation of scales: the formation of galaxies and halos is assumed to depend only on the properties of the Lagrangian environment on some small length scale. Larger scale perturbations are only relevant as far as they determine the distribution of small scale environments. We refer to this as the ‘scale-split’ assumption.

This assumption makes it possible to define scale-independent bias parameters. For example, b₁ describes how the likelihood of forming a galaxy responds when changing only the linear density contrast δ of the relevant smaller scale environment, while keeping other aspects (e.g. L) constant. As all larger scale density perturbations affect this aspect in the same way b₁ is independent of scale.

Physically, the scale-split assumption has to become invalid for sufficiently small smoothing scales. For example, halo formation may respond differently to perturbations that are smaller than their Lagrangian radius, than to those which are larger.

Here, we want to show that the scale-split assumption also becomes mathematically inconsistent beyond some scale for a given set of variables: Recall that we have shown in Equation (25), how κ₂, the variance of the galaxy environment distribution, changes with scale. The predicted variance is only well defined, if

$$\begin{array}{cc}\hfill {\sigma }^2\le {\sigma }_{\max ,\delta }^2& =\frac{1}{-{\beta }_2}.\hfill \end{array}$$(79)

This becomes zero if σ = σ_max, δ. Mathematically, at this scale, the environment distribution and the bias function have to become Dirac-delta functions. We may understand σ_max as the scale where (formally) all information that is relevant for galaxy formation has been accounted for and the biasing becomes deterministic in density.

Beyond this scale, the PBS predicts a galaxy environment distribution with negative variance – making any density-only bias model mathematically inconsistent. It may seem as if it was possible to set up expansion bias models for arbitrary high σ, but this is only so, since the negativity of the bias function makes it formally possible to have a negative variance. Actual galaxies will strictly obey κ₂ > 0 and therefore, the response of density-only models has to become scale-dependent latest at σ_max, δ – but likely already earlier. We call this the ‘scale-split break-down scale’.

The break-down scale is different if additional variables are considered. For the multivariate distribution of density and Laplacian (δ, L), we have to require the covariance matrix of the galaxy environment distribution to be positive and semi-definite. For this, it is at least necessary to have:

$$\begin{array}{cc}\hfill det({\mathit{\kappa }}_2)& \ge 0,\hfill \\ \hfill \iff det(\mathsf{1}+\mathsf{C}{\mathit{\beta }}_2)& \ge 0,\hfill \end{array}$$(80)

which appears slightly more complicated, but it will also be violated at some finite damping scale.

To measure the break-down scale, we infer for each halo mass the bias parameters at three different scales k_d ∈ (0.15, 0.2, 0.25) h Mpc⁻¹. Then, we evaluate the co-variance matrix of the background distribution at 500 different equally log-spaced damping scales between 10⁻³ h Mpc⁻¹ and 10¹ h Mpc⁻¹ and determine the earliest damping scale k_max where the covariance is large enough to violate Equation (80) or Equation (81). We show the corresponding results as shaded contours in Figure 9, where we additionally mark the characteristic wave-number of halos as:

$$\begin{array}{cc}\hfill k_{\mathrm{halo}}& =\frac{\pi }{R_{\mathrm{halo}}},\hfill \end{array}$$(81)

where R_halo is the Lagrangian radius that encloses the halo mass, M_200b.

Comparing the measurements at different damping scales, we find that the inferred scale is reasonably converged with the scale that we measured the bias parameters at. The density-only break-down scale is typically a factor two smaller than k_halo and seems roughly proportional to it. We therefore conclude that any Lagrangian local in matter density (LLIMD) bias model with scale-independent bias parameters has to break down at a length scale roughly a factor two larger than Lagrangian radii of the considered halos⁵.

On the other hand, the (δ, L) case shows a notably different break-down scale. It scales only weakly with halo mass and it ranges only between k_max ∼ 0.3 − 0.5 h Mpc⁻¹. Note that this corresponds roughly to the Lagrangian scale of halos with M ∼ 2 × 10¹⁴ h⁻¹ M_⊙, where b₁ ∼ 1. For halos above M ≳ 3 × 10¹³ h⁻¹ M_⊙ including the Laplacian increases k_max relative to the density-only case, but for lower masses, it decreases it.

Comparing with the right panel of Figure 4, we notice that the bias parameters already become scale-dependent at notably smaller damping scales than k_max. For example, for M_200b ∼ 10¹⁵ h⁻¹ M_⊙ with b₁ ∼ 3, the measurements of β₂ are already scale-dependent beyond k_d ≳ 0.15 h Mpc⁻¹, whereas the mathematical break-down scale is k_max ∼ 0.3 h Mpc⁻¹. We therefore suggest that considerable care should be taken when setting up bias models close to the break-down scale. This is particularly relevant for hybrid methods, which might, in principle, allow for galaxy clustering to b described at notably smaller scales than are usually considered in a perturbative schemes (e.g. Modi et al. 2020; Zennaro et al. 2022).

4. Estimators for tensorial bias parameters

In Sect. 2, we explain how we inferred the general estimators of bias parameters associated with scalar variables (e.g. the density and Laplacian) with spatial corrections of any order. In Sect. 3, we show how these can be used to obtain reliable bias measurements from a single simulation.

However, the theory in Sect. 2 does not explain how to obtain estimators for parameters such as the tidal bias, b_K², which is defined as the response to

$$\begin{array}{cc}\hfill K^2& =\mathrm{tr}\left(\mathsf{K}·\mathsf{K}\right),\hfill \end{array}$$(82)

$$\begin{array}{cc}\hfill \mathsf{K}& =(\mathbf{\nabla }\otimes \mathbf{\nabla })\varphi -\frac{\delta }{3}{\mathsf{I}}_2,\hfill \end{array}$$(83)

where ϕ is the displacement potential, I₂ is the unit matrix, and K is the traceless tidal tensor. In this section, we present a general scheme to measure the related bias parameter b_K² and any other bias parameters that follow from contractions of derivatives of the potential field.

To achieve this, it is not optimal to consider directly the distribution of such scalar contracted quantities, since these distributions may get quite complicated. For example, p(K²) is not a Gaussian distribution, but rather a χ² distribution with five degrees of freedom. Furthermore, it is not immediately obvious how to define partial derivatives with respect to such variables, since partial derivatives may depend on what other terms are kept fixed. For example, it is not clear how a term like K³ = tr(K⋅K⋅K) would be derived with respect to K². It seems therefore difficult to generalise Equation (12) in this manner.

However, a clear and general framework for measuring such ‘tensorial’ bias terms can be developed by instead considering the full (quite high dimensional) distribution of the tidal tensor and its derivatives. Since the potential field of the early Universe is a Gaussian random field, these must follow a multivariate Gaussian distribution so that it is easy to compute derivatives of the distribution function in a general manner. The resulting ‘bias tensors’ can be decomposed into isotropic tensors that each have a one to one correspondence with traditionally used bias parameters. In this section, we introduce the needed mathematical notions step by step and will provide estimators for a few selected bias terms.

4.1. Tensorial bias expansion

We write the bias expansion in tensorial form as

$$\begin{array}{cc}\hfill F& =1+{\mathsf{B}}_{\mathsf{T}}\mathsf{T}+\frac{1}{2}\mathsf{T}{\mathsf{B}}_{\mathsf{T}\mathsf{T}}\mathsf{T}+\cdots +{\mathsf{B}}_{\mathsf{S}}\mathsf{S}+\cdots +{\mathsf{B}}_{\mathsf{R}}\mathsf{R}+\cdots \hfill \end{array}$$(84)

Here, again, F = n_g/n_g, 0 (in a separate-universe sense) and

$$\begin{array}{cc}\hfill \mathsf{T}& =(\mathrm{\nabla }\otimes \mathrm{\nabla })\varphi ,\hfill \end{array}$$(85)

$$\begin{array}{cc}\hfill \mathsf{S}& =(\mathrm{\nabla }\otimes \mathrm{\nabla }\otimes \mathrm{\nabla })\varphi ,\hfill \end{array}$$(86)

$$\begin{array}{cc}\hfill \mathsf{R}& =(\mathrm{\nabla }\otimes \mathrm{\nabla }\otimes \mathrm{\nabla }\otimes \mathrm{\nabla })\varphi ,\hfill \end{array}$$(87)

where B_T is a canonical bias tensor of rank of 2, B_TT of rank of 4, B_S of the rank of 3, and so on. This is also where an omitted product sign indicates a product over the indices of the last fully symmetric part of the first tensor and the first fully symmetric part of the second tensor. For example:

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{T}}\mathsf{T}& ={\mathsf{B}}_{\mathsf{T}}\stackrel{(2)}{·}\mathsf{T}=\sum _{\mathit{ij}}{B}_{\mathsf{T},\mathsf{ij}}{T}_{\mathit{ij}},\hfill \\ \hfill \mathsf{T}{\mathsf{B}}_{\mathsf{T}\mathsf{T}}\mathsf{T}& =\mathsf{T}\stackrel{(2)}{·}{\mathsf{B}}_{\mathsf{T}\mathsf{T}}\stackrel{(2)}{·}\mathsf{T}=\sum _{\mathit{ijkl}}{T}_{\mathit{ij}}{B}_{\mathsf{T}\mathsf{T},\mathsf{ijkl}}{T}_{\mathit{kl}},\hfill \\ \hfill {\mathsf{B}}_{\mathsf{S}}\mathsf{S}& ={\mathsf{B}}_{\mathsf{S}}\stackrel{(3)}{·}\mathsf{S}=\sum _{\mathit{ijk}}{B}_{\mathsf{S},\mathsf{ijk}}{S}_{\mathit{ijk}},\hfill \end{array}$$(88)

where we can use $\stackrel{(n)}{·}$ to explicitly denote the number of dimensions that are contracted.

We note that just as before, the bias tensors correspond to derivatives of the large-scale bias function, for instance,

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{T}}& ={\frac{\partial }{\partial \mathsf{T}}F|}_{\mathsf{T}=\mathsf{0}},\hfill \end{array}$$(89)

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{T}\mathsf{T}}& ={(\frac{\partial }{\partial \mathsf{T}}\otimes \frac{\partial }{\partial \mathsf{T}})F|}_{\mathsf{T}=\mathsf{0}}.\hfill \end{array}$$(90)

4.2. Isotropic tensors

Due to the isotropy of the Universe, each of the bias tensors has to be isotropic. That means that a bias tensor should be identical when measured from a rotated frame of reference. For example, it has to be

$$\begin{array}{cc}\hfill {\mathsf{U}}^T{\mathsf{B}}_{\mathsf{T}}\mathsf{U}& ={\mathsf{B}}_{\mathsf{T}}\hfill \end{array}$$(91)

for any rotation matrix U. From this, it follows immediately that B_T has to be proportional to the unit matrix I₂ (where the proportionality constant is equal to b₁). In general, a rank n tensor A is isotropic if it holds for any rotation matrix, U:

$$\begin{array}{cc}\hfill A_{abc\dots }{U}_{\mathit{ai}}{U}_{\mathit{bj}}{U}_{\mathit{ck}}\dots & =A_{ijk\dots }\hfill \end{array}$$(92)

Here, we have used Einstein’s sum convention, where the rotation matrix is applied to each index of A individually. To express whether some tensor is isotropic, we define 𝕌_n as the space of all tensors that are isotropic and of a rank, n.

In general, all isotropic tensors can be decomposed in index notation through different combinations of the Kronecker-delta symbol δ_ij and the Levi-Civita symbol, ϵ_ijk. For example

$$\begin{array}{cc}\hfill \mathsf{A}\in {\mathbb{U}}_2& \Rightarrow A_{\mathit{ij}}=a{\delta }_{\mathit{ij}},\hfill \end{array}$$(93)

$$\begin{array}{cc}\hfill \mathsf{A}\in {\mathbb{U}}_3& \Rightarrow A_{\mathit{ijk}}=a{ϵ}_{\mathit{ijk}},\hfill \end{array}$$(94)

$$\begin{array}{cc}\hfill \mathsf{A}\in {\mathbb{U}}_4& \Rightarrow A_{\mathit{ijkl}}=a{\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}}+b{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+c{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}},\hfill \end{array}$$(95)

where a, b, c ∈ ℝ. A compact way of writing the same type of statement can be:

$$\begin{array}{c}\hfill {\mathbb{U}}_4=\mathrm{span}(\{{\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}},{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}},{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}\}),\end{array}$$(96)

which signifies that 𝕌₄ is the space of tensors that can be reached through linear combinations of the indicated tensors and we say that these are basis tensors of 𝕌₄.

It is easy to see that we should be able to decompose each of the bias tensors into a small number of independent scalars that multiply the basis tensors and that correspond to traditional bias parameters. However, before performing such a decomposition we need to further consider the symmetries of the tensors.

4.3. Symmetric isotropic tensors

Since the bias tensors correspond to derivatives of a scalar function with respect to symmetric tensors, they have to obey the same symmetries. For example, the tensor B_TT has to be symmetric in the first two and last two indices⁶:

$$\begin{array}{cc}\hfill B_{\mathsf{T}\mathsf{T},\mathsf{ijkl}}& =B_{\mathsf{T}\mathsf{T},\mathsf{jikl}}=B_{\mathsf{T}\mathsf{T},\mathsf{ijlk}}=B_{\mathsf{T}\mathsf{T},\mathsf{jilk}}.\hfill \end{array}$$(97)

Therefore, B_TT cannot be every tensor from the tensorspace 𝕌₄, but only from the subspace 𝕍₂₂ ⊆ 𝕌₄ of isotropic rank four tensors that have the 2, 2 symmetry.

To formalise this, we define 𝕍_n ⊂ 𝕌_n as the space of all isotropic tensor of rank n that are additionally fully symmetric under any permutation of indices, specifically:

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{T}}& \in {\mathbb{V}}_2,\hfill \end{array}$$(98)

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{S}}& \in {\mathbb{V}}_3,\hfill \end{array}$$(99)

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{R}}& \in {\mathbb{V}}_4.\hfill \end{array}$$(100)

We note that 𝕍₃ = {0}, namely, the only symmetric rank three isotropic tensor is 0. Furthermore, we define 𝕍_nm as the space of all isotropic tensors of rank n + m that are symmetric in the first n indices and the last m indices. For example,

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{T}\mathsf{T}}& \in {\mathbb{V}}_{22}.\hfill \end{array}$$(101)

To find a basis for some symmetric isotropic tensor space of rank n, we can consider all basis tensors of 𝕌_n, symmetrise these and then discard any tensors that are duplicate or zero. We define the symmetrisation operator S_n which symmetrises a tensor in n indices. For instance, we have:

$$\begin{array}{cc}\hfill S_2(M_{\mathit{ij}})& =\frac{1}{2}(M_{\mathit{ij}}+M_{\mathit{ji}}),\hfill \end{array}$$(102)

$$\begin{array}{cc}\hfill S_3(M_{\mathit{ijk}})& =\frac{1}{6}(M_{\mathit{ijk}}+M_{\mathit{ikj}}+M_{\mathit{jik}}+M_{\mathit{jki}}+M_{\mathit{kij}}+M_{\mathit{jki}}).\hfill \end{array}$$(103)

Furthermore, we defined a double symmetrisation operator, S_nm, that symmetrises the first n indices and the last m indices. For example, the effect of S₂₂ onto a rank four tensor M is given in the index notation by

$$\begin{array}{cc}\hfill S_{22}(M_{\mathit{ijkl}})& =\frac{1}{4}(M_{\mathit{ijkl}}+M_{\mathit{jikl}}+M_{\mathit{ijlk}}+M_{\mathit{jilk}}).\hfill \end{array}$$(104)

It acts on the basis tensors of 𝕌₄ in the following ways:

$$\begin{array}{cc}\hfill S_{22}({\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}})& ={\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}}=:I_{22,ij},\hfill \end{array}$$(105)

$$\begin{array}{cc}\hfill S_{22}({\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}})& =\frac{1}{2}({\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}+{\delta }_{\mathit{il}}{\delta }_{\mathit{jk}})=:I_{2=2,ijkl},\hfill \end{array}$$(106)

$$\begin{array}{cc}\hfill S_{22}({\delta }_{\mathit{il}}{\delta }_{\mathit{jk}})& =\frac{1}{2}({\delta }_{\mathit{il}}{\delta }_{\mathit{jk}}+{\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}})=I_{2=2,ijkl}\hfill \end{array}$$(107)

where we have assumed that δ_ij is itself symmetric. The symmetrised second and third basis tensor are identical. Therefore, 𝕍₂₂ has (unlike 𝕌₄) only two basis tensors:

$$\begin{array}{cc}\hfill {\mathbb{V}}_{22}& =\mathrm{span}(\{{\mathsf{I}}_{22},{\mathsf{I}}_{2=2}\}).\hfill \end{array}$$(108)

The main difference between δ_ijδ_kl and the two tensors δ_ikδ_jl and δ_ilδ_jk is that the delta symbols used to define δ_ijδ_kl each connect internally inside the groups of symmetric indices and have no connections between the groups. (Recall that i ↔ j and k ↔ l are to be symmetrised.) However, δ_ikδ_jl and δ_ilδ_jk each have zero group internal connections, but two connections between the symmetry groups. In fact, the symmetrisation operation identifies all terms that have the same number of intra- and inter- group connections and therefore we only need to consider how many independent possibilities exist to connect the symmetry groups to identify the basis tensors of any space, 𝕍_mn. We can therefore represent the basis tensors of each tensor space through a simple diagram as illustrated in Figure 10. We label these tensors through a symbol that shows the number of inter group connections in the index:

$$\begin{array}{cc}\hfill {\mathbb{V}}_{33}& =\mathrm{span}(\{{\mathsf{I}}_{3-3},{\mathsf{I}}_{3\equiv 3}\}),\hfill \end{array}$$(109)

$$\begin{array}{cc}\hfill {\mathbb{V}}_{44}& =\mathrm{span}(\{{\mathsf{I}}_{44},{\mathsf{I}}_{4=4},{\mathsf{I}}_{4≣4}\}),\hfill \end{array}$$(110)

$$\begin{array}{cc}\hfill {\mathbb{V}}_{24}& =\mathrm{span}(\{{\mathsf{I}}_{24},{\mathsf{I}}_{2=4}\}).\hfill \end{array}$$(111)

We explain in Appendix B.1 how to construct tensors with three or more symmetry groups. For these cases the procedure is mostly analogous, but it additionally needs to be considered whether there is symmetry with respect to permutations of the different symmetry groups.

4.4. Orthogonal basis

The basis tensors highlighted in the last section are sufficient to uniquely decompose the bias tensors. However, it is advantageous to define the basis tensors in such a way that they are orthogonal to each other.

For example, we can decompose an isotropic tensor M ∈ 𝕍₂₂ as

$$\begin{array}{cc}\hfill \mathsf{M}& =A{\mathsf{I}}_{22}+B{\mathsf{I}}_{2=2}.\hfill \end{array}$$(112)

Its full contraction with I₂₂ is

$$\begin{array}{cc}\hfill \mathsf{M}\stackrel{(4)}{·}{\mathsf{I}}_{22}& =A({\mathsf{I}}_{22}\stackrel{(4)}{·}{\mathsf{I}}_{22})+B({\mathsf{I}}_{2=2}\stackrel{(4)}{·}{\mathsf{I}}_{22}),\hfill \end{array}$$(113)

$$\begin{array}{cc}\hfill & =9A+3B,\hfill \end{array}$$(114)

which includes contributions of both A and B. We could infer A and B from M by additionally contracting with I_2 = 2 and then solving the resulting system of equations.

However, it would be desirable to have a simple way of finding the coefficients. Therefore, we define an orthogonal basis such that

$$\begin{array}{cc}\hfill {\mathsf{J}}_a\stackrel{(r)}{·}{\mathsf{J}}_b& =\{\begin{array}{c}\left||{\mathsf{J}}_a\right|{|}^2\text{if}a=b,\hfill \\ 0\text{otherwise},\hfill \end{array}\hfill \end{array}$$(115)

where the contraction goes over the full rank r of J_a and J_b (which must have equal rank). Such a basis can be found through Gram-Schmidt orthogonalisation and is for example given for 𝕍₂₂ by:

$$\begin{array}{cc}\hfill {\mathsf{J}}_{22}& ={\mathsf{I}}_{22},\hfill \end{array}$$(116)

$$\begin{array}{cc}\hfill {\mathsf{J}}_{2=2}& ={\mathsf{I}}_{2=2}-\frac{1}{3}{\mathsf{I}}_{22}.\hfill \end{array}$$(117)

Then we can easily find the coefficients as

$$\begin{array}{cc}\hfill \mathsf{M}& =\alpha {\mathsf{J}}_{22}+\beta {\mathsf{J}}_{2=2},\hfill \end{array}$$(118)

$$\begin{array}{cc}\hfill \alpha & =\frac{\mathsf{M}\stackrel{(\mathsf{4})}{·}{\mathsf{J}}_{\mathsf{22}}}{\left||{\mathsf{J}}_{22}\right|{|}^2}=\frac{1}{9}\mathsf{M}\stackrel{(4)}{·}{\mathsf{J}}_{22},\hfill \end{array}$$(119)

$$\begin{array}{cc}\hfill \beta & =\frac{\mathsf{M}\stackrel{(\mathsf{4})}{·}{\mathsf{J}}_{\mathsf{2}=\mathsf{2}}}{\left||{\mathsf{J}}_{2=2}\right|{|}^2}=\frac{1}{5}\mathsf{M}\stackrel{(4)}{·}{\mathsf{J}}_{2=2}.\hfill \end{array}$$(120)

We note that the coefficients inferred in this way are different than the A and B and the chosen basis will also affect the inferred bias parameters.

It is worth noting that a decomposition of a bias tensor into different bases,

$$\begin{array}{cc}\hfill {\mathsf{B}}_{\mathsf{T}\mathsf{T}}& =b_{{\mathsf{I}}_{22}}{\mathsf{I}}_{22}+b_{{\mathsf{I}}_{2=2}}{\mathsf{I}}_{2=2},\hfill \\ \hfill & =b_{{\mathsf{J}}_{22}}{\mathsf{J}}_{22}+b_{{\mathsf{J}}_{2=2}}{\mathsf{J}}_{2=2},\hfill \end{array}$$(121)

corresponds to different independent scalar variables appearing in the bias expansion. For example, if we track the corresponding term from Equation (85), we find

$$\begin{array}{cc}\hfill \mathsf{T}{\mathsf{B}}_{\mathsf{T}\mathsf{T}}\mathsf{T}& =b_{{\mathsf{I}}_{22}}\mathrm{tr}{\left(\mathsf{T}\right)}^2+b_{{\mathsf{I}}_{2=2}}\mathrm{tr}\left({\mathsf{T}}^2\right),\hfill \\ \hfill & =b_{{\mathsf{I}}_{22}}{\delta }^2+b_{{\mathsf{I}}_{2=2}}{T}^2,\hfill \end{array}$$(122)

$$\begin{array}{cc}\hfill \mathsf{T}{\mathsf{B}}_{\mathsf{T}\mathsf{T}}\mathsf{T}& =b_{{\mathsf{J}}_{22}}\mathrm{tr}{\left(\mathsf{T}\right)}^2+b_{{\mathsf{J}}_{2=2}}(\mathrm{tr}\left({\mathsf{T}}^2\right)-\frac{{\delta }^2}{3}),\hfill \\ \hfill & =b_{{\mathsf{J}}_{22}}{\delta }^2+b_{{\mathsf{J}}_{2=2}}{K}^2.\hfill \end{array}$$(123)

Thus, the non-orthogonal basis leads to δ and T² as independent terms, whereas the orthogonal basis leads to using δ and K² as independent terms. The corresponding bias parameters can be directly identified as b_J22 = b₂ and $\frac{1}{2}{b}_{{\mathsf{J}}_{2=2}}=b_{K^2}$.

We list an overview of the orthogonal basis tensors that we use here in Table 1. While it is possible to derive the full algebra of products between isotropic tensors from the index representations, this is rather cumbersome. Therefore, we have written a short code that creates explicit numerical representations of these tensors and with which it is easy to evaluate (and decompose) different types of products of these tensors. We use this together with the symbolic computer algebra system SYMPY (Meurer et al. 2017) to systematically compute symbolic representations of expressions in the following sections. This code is openly available⁷.

4.5. Bias estimators

Given the orthogonal basis tensors, we define the tensorial bias parameters through

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_X}& ={\frac{{\partial }^{N}F}{\partial {\mathsf{T}}_0^N}\stackrel{(2N)}{·}\frac{{\mathsf{J}}_X}{\left||{\mathsf{J}}_X\right|{|}^2}|}_{{\mathsf{T}}_0=0},\hfill \end{array}$$(124)

where the rank of J_X is 2N. While for scalar terms, as in Equation (3), the bias parameters are simply given by derivatives of the galaxy number with respect to the corresponding scalar, more generally we define bias parameters for any tensorial terms as derivatives of the galaxy number with respect to the tidal tensor (or higher spatial order tensors) that are contracted with the corresponding isotropic tensor.

In complete analogy to Equation (46) it follows that these parameters can be estimated as

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_X}& ={(-1)}^N{\langle \frac{1}{p}\frac{\partial p}{\partial {\mathsf{T}}^N}\stackrel{(2N)}{·}\frac{{\mathsf{J}}_X}{\left||{\mathsf{J}}_X\right|{|}^2}\rangle }_g,\hfill \end{array}$$(125)

where p is the full distribution of the tidal tensor, T.

Bias terms corresponding to higher spatial derivatives like R and S can be defined analogously.

4.6. Tidal bias

We show in Appendices B.2 and B.3 that the distribution of the tidal tensor is given by

$$\begin{array}{cc}\hfill p(\mathsf{T})& =Nexp(-\frac{1}{2}{\mathsf{T}}^T{\mathsf{C}}_{\mathsf{T}}^{+}\mathsf{T}),\hfill \end{array}$$(126)

where C_T⁺ (the generalised inverse of the covariance matrix, as explained in the Appendix) is given by

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}}^{+}& =\frac{1}{{\sigma }_0^2}{\mathsf{J}}_{22}+\frac{15}{2{\sigma }_0^2}{\mathsf{J}}_{2=2}.\hfill \end{array}$$(127)

Taking derivatives yields

$$\begin{array}{cc}\hfill \frac{1}{p}\frac{\partial p}{\partial \mathsf{T}}& =-{\mathsf{C}}_{\mathsf{T}}^{+}\mathsf{T},\hfill \end{array}$$(128)

$$\begin{array}{cc}\hfill \frac{1}{p}\frac{{\partial }^{2}p}{\partial {\mathsf{T}}^2}& =(\mathsf{T}{\mathsf{C}}_{\mathsf{T}}^{+})\otimes ({\mathsf{C}}_{\mathsf{T}}^{+}\mathsf{T})-{\mathsf{C}}_{\mathsf{T}}^{+}\hfill \end{array}$$(129)

and we find

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_2}=b_1& ={\langle \frac{{\mathsf{J}}_2{\mathsf{C}}^{+}\mathsf{T}}{\left||{\mathsf{J}}_2\right|{|}^2}\rangle }_{\mathrm{g}}={\langle \frac{{\mathsf{J}}_2\mathsf{T}}{{\sigma }_0^2}\rangle }_{\mathrm{g}}={\langle \frac{\delta }{{\sigma }_0^2}\rangle }_{\mathrm{g}}\hfill \end{array}$$(130)

where we used the inner product between two isotropic tensors: ${\mathsf{J}}_2\stackrel{(2)}{·}{\mathsf{J}}_{22}=3{\mathsf{J}}_2$. We evaluated the terms of this type systematically with numerical representations as described in Sect. 4.4. Naturally, this estimator of b_J2 = b₁ is consistent with the one that we have inferred previously in Equation (15).

At the second order, we find the following estimators:

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_{22}}=b_2& ={\langle \frac{{\mathsf{T}}^T{\mathsf{J}}_{22}\mathsf{T}-{\sigma }_0^2}{{\sigma }_0^4}\rangle }_{\mathrm{g}}\hfill \\ \hfill & ={\langle \frac{{\delta }^2-{\sigma }_0^2}{{\sigma }_0^4}\rangle }_{\mathrm{g}},\hfill \end{array}$$(131)

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_{2=2}}=2b_{K^2}& =\frac{15}{4{\sigma }_0^4}{\langle 3{\mathsf{T}}^T{\mathsf{J}}_{2=2}\mathsf{T}-2{\sigma }_0^2\rangle }_{\mathrm{g}},\hfill \end{array}$$(132)

$$\begin{array}{cc}& =\frac{15}{4{\sigma }_0^4}{\langle 3K^2-2{\sigma }_0^2\rangle }_{\mathrm{g}}.\hfill \end{array}$$(133)

We note that conventionally the bias parameter, b_K², is defined so that it appears with a pre-factor K² in the bias expansion, whereas our parameter b_J2 = 2 appears with a factor $\frac{1}{2}{K}^2$ in the expansion (after contracting the corresponding tensors), so that it is $b_{K^2}=\frac{1}{2}{b}_{{\mathsf{J}}_{2=2}}$. Although we think that our fore-factor convention makes more sense in principle, since b_J2 = 2 corresponds to a second derivative term, we still present our results in terms of the conventional notation (e.g. b_K²) throughout this paper.

We note that just as the large-scale value of δ/σ² at an object’s location is motivated by Equation (15) and sometimes called the ‘bias-per-object’ (e.g. Paranjape et al. 2018; Contreras et al. 2023), we could refer in a similar manner to the value of

$$\begin{array}{cc}\hfill b_{K^2,1h}& =\frac{15}{8{\sigma }_0^4}(3K^2-2{\sigma }_0^2)\hfill \end{array}$$(134)

as the value of the ‘tidal bias-per-object’.

Just as estimators with spatial corrections could be obtained for the density biases by considering the joint distribution of p(δ, L), we can obtain higher spatial order estimators for the tidal bias by considering the joint distribution of second and fourth potential derivatives p(T, R) and evaluating Equation (125) for this distribution. We show how we derived the bias estimator for this case in appendices B.5 and B.6. The resulting estimator is:

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_{2=2}}& =\frac{15}{4{\sigma }_{\ast }^8}{\langle 3K^2{\sigma }_2^4+6{\varphi }_{2=4}{\sigma }_1^2{\sigma }_2^2+3{\varphi }_{4=4}{\sigma }_1^4-2{\sigma }_2^2{\sigma }_{\ast }^4\rangle }_{\mathrm{g}},\hfill \end{array}$$(135)

where

$$\begin{array}{cc}\hfill {\varphi }_{2=4}& :=\mathsf{T}{J}_{2=4}\mathsf{R}=\sum _{\mathit{ij}}({\partial }_i{\partial }_j\varphi )({\partial }_i{\partial }_j\delta )-\frac{1}{3}\delta {\mathrm{\nabla }}^2\delta ,\hfill \end{array}$$(136)

$$\begin{array}{cc}\hfill {\varphi }_{4=4}& :=\mathsf{T}{J}_{4=4}\mathsf{R}=\sum _{\mathit{ij}}{({\partial }_i{\partial }_j\delta )}^2-\frac{1}{3}{({\mathrm{\nabla }}^2\delta )}^2,\hfill \end{array}$$(137)

and ${\sigma }_{\ast }^2={\sigma }_0^2-{\sigma }_1^4/{\sigma }_2^2$.

4.7. Estimators for third derivative terms

The distribution of third derivatives of the potential S is derived in Appendix B.4 and is given by a multivariate Gaussian with the generalised inverse covariance matrix

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{S}}^{+}& =\frac{3}{{\sigma }_1^2}{\mathsf{J}}_{3-3}+\frac{35}{2{\sigma }_1^2}{\mathsf{J}}_{3\equiv 3}.\hfill \end{array}$$(138)

We note that it is fine to consider the third spatial derivatives independently of the second derivatives, since the joint distribution factorises as:

$$\begin{array}{cc}\hfill p(\mathsf{T},\mathsf{S})& =p(\mathsf{T})p(\mathsf{S}).\hfill \end{array}$$(139)

We find bias estimators by evaluating Equation (125) as:

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_{3-3}}& =\frac{3}{{\sigma }_1^4}{\langle \mathsf{S}{\mathsf{J}}_{\mathsf{3}-\mathsf{3}}\mathsf{S}-{\sigma }_{\mathsf{1}}^{\mathsf{2}}\rangle }_{\mathrm{g}}\hfill \\ \hfill & =\frac{3}{{\sigma }_1^4}{\langle {(\mathrm{\nabla }\delta )}^2-{\sigma }_1^2\rangle }_{\mathrm{g}},\hfill \end{array}$$(140)

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_{3\equiv 3}}& =\frac{35}{{\sigma }_1^4}{\langle 5\mathsf{S}{\mathsf{J}}_{3\equiv 3}\mathsf{S}-2{\sigma }_1^2\rangle }_{\mathrm{g}}\hfill \\ \hfill & =\frac{35}{{\sigma }_1^4}{\langle 5{\varphi }_{3\equiv 3}-2{\sigma }_1^2\rangle }_{\mathrm{g}},\hfill \end{array}$$(141)

where

$$\begin{array}{cc}\hfill {\varphi }_{3\equiv 3}& =\sum _{\mathit{ijk}}{({\partial }_{\mathit{ijk}}\varphi )}^2-\frac{3}{5}\sum _i{({\partial }_i\delta )}^2,\hfill \end{array}$$(142)

and where we note again that conventions may differ by a factor 2 so that we have, for instance, b_J3 − 3 = 2b_(∇δ)² if b_(∇δ)² is the bias parameter that would be multiplied (∇δ²) without any fore-factor in a bias expansion.

4.8. Tensorial cumulant biases

As for the scalar case, it is possible to define tensorial cumulant bias parameters instead of tensorial canonical bias parameters. Unfortunately, interesting differences arise only at orders that are higher than the ones that we have considered here. However, for the benefit of future studies, we still want to briefly line out, how to define tensorial cumulant biases.

Analogously to Equation (85), we considered a tensorial expansion of log F as

$$\begin{array}{c}\hfill logF=1+{\beta }_{\mathsf{T}}\mathsf{T}+\frac{1}{2}\mathsf{T}{\beta }_{\mathsf{T}\mathsf{T}}\mathsf{T}+\cdots +{\beta }_{\mathsf{S}}\mathsf{S}+\cdots +{\beta }_{\mathsf{R}}\mathsf{R}+\cdots \end{array}$$(143)

Taking derivatives at zero and identifying coefficients yields the relations

$$\begin{array}{cc}\hfill {\beta }_{\mathsf{T}}& ={\mathsf{B}}_{\mathsf{T}},\hfill \end{array}$$(144)

$$\begin{array}{cc}\hfill {\beta }_{\mathsf{T}\mathsf{T}}& ={\mathsf{B}}_{\mathsf{T}\mathsf{T}}-{\mathsf{B}}_{\mathsf{T}}\otimes {\mathsf{B}}_{\mathsf{T}}.\hfill \end{array}$$(145)

To decompose this into isotropic tensors we may contract these equations with the isotropic tensors to find the relations

$$\begin{array}{cc}\hfill {\beta }_1={\beta }_{{\mathsf{J}}_2}& =b_{{\mathsf{J}}_2}=b_1,\hfill \end{array}$$(146)

$$\begin{array}{cc}\hfill {\beta }_2={\beta }_{{\mathsf{J}}_{22}}& =b_{{\mathsf{J}}_{22}}-{\beta }_{{\mathsf{J}}_2}^2=b_2-b_1^2,\hfill \end{array}$$(147)

$$\begin{array}{cc}\hfill 2{\beta }_{K^2}={\beta }_{{\mathsf{J}}_{2=2}}& =b_{{\mathsf{J}}_{2=2}}=2b_{K^2}.\hfill \end{array}$$(148)

This is consistent with what we derived earlier. Tensorial cumulant biases become more interesting at third order, where third derivative of the bias Function yields

$$\begin{array}{cc}\hfill {\beta }_{\mathsf{T}\mathsf{T}\mathsf{T}}& ={\mathsf{B}}_{\mathsf{T}\mathsf{T}\mathsf{T}}-3S_{222}({\mathsf{B}}_{\mathsf{T}}\otimes {\mathsf{B}}_{\mathsf{T}\mathsf{T}})+2S_{222}({\mathsf{B}}_{\mathsf{T}}\otimes {\mathsf{B}}_{\mathsf{T}}\otimes {\mathsf{B}}_{\mathsf{T}}).\hfill \end{array}$$(149)

If we contract this with the isotropic tensor S₂₂₂(J_2 = 2 ⊗ J₂), we have

$$\begin{array}{cc}\hfill {\beta }_{{\mathsf{J}}_{2=22}}& =b_{{\mathsf{J}}_{2=22}}-3b_{{\mathsf{J}}_{\mathsf{2}}}{b}_{{\mathsf{J}}_{\mathsf{2}=\mathsf{2}}}.\hfill \end{array}$$(150)

Therefore, if the full multivariate bias function is close to Gaussian, also tensorial bias parameters of this form should be expected to be close to zero. Dealing with such objects goes beyond the framework of this study and we set this aside as an interesting proposition for a future work.

5. Measurements of tensorial bias parameters

In this section, we use the estimators derived in Sect. 4 to perform measurements of tensorial bias parameters. Here, we focus on the tidal bias, b_K², and on the biases associated with the third potential derivatives, $b_{{(\mathrm{\nabla }\delta )}^2}=\frac{1}{2}{b}_{{\mathsf{J}}_{3-3}}$ and b_J3 ≡ 3. For the tidal bias, we have additionally derived the estimator with spatial corrections of the order of 2, so that we can use it to validate that the presented theory indeed behaves as expected. Furthermore, we chose to briefly check the effect of assembly bias on b_K², which has previously been reported to be quite significant. For these measurements, we used the same simulations and jackknife technique (described in Sects. 3.1 and 3.2).

5.1. Tidal bias

We measured the tidal bias, $b_{K^2}=\frac{1}{2}{b}_{{\mathsf{J}}_{2=2}}$, with the estimators from Equation (133) at a zeroth-spatial order and Equation (135) at a spatial order of 2. We show the results as a function of mass in Figure 11. The zeroth-order estimator exhibits a strong scale dependence at scales k_d ≳ 0.15 h Mpc⁻¹. However, overall it seems that the tidal bias is very small in amplitude and that a slightly negative Lagrangian tidal bias is preferred at high halo masses. On the other hand, the second-order spatial estimator has a much lower scale dependence and seems mostly consistent across different damping scales, potentially with a slight inconsistency at very high halo masses and large k_d. This comes at the price of a significantly increased uncertainty in the measured bias values. However, this demonstrates that the estimators we derived in Sect. 4 do indeed behave the way we would expect. More precise measurements of the tidal bias could easily be obtained by employing a larger simulation volume or by averaging over several realisations.

In Figure 12, we show the co-evolution relation between b_K² and b₁. Here, we show for each estimator only the smallest two scales that appear reasonably unbiased in Figure 11. These are k_d = 0.1 h Mpc⁻¹ and k_d = 0.15 h Mpc⁻¹ for the zeroth-order spatial estimator and k_d = 0.25 h Mpc⁻¹ and k_d = 0.3 h Mpc⁻¹ for the second-order spatial estimator. It is noteworthy that the estimators of different orders are in agreement with each other. For the larger k_d values, the estimators have smaller error bars, but do appear slightly biased at large values of b₁ ≳ 2. A rough approximation to the observed functional shape is given by

$$\begin{array}{c}\hfill b_{K^2}=-0.05b_1^2.\end{array}$$(151)

However, the uncertainty is certainly quite high. However, in comparison to some previous literature estimates, we find very small values of bK2. The first measurement that we compare with is from Modi et al. (2017). Their relation is stated for Eulerian bias parameters bK2E(b1E) and we convert it to a Lagrangian relation via b1 = b1E − 1 and bK2 = bK2E + 2/7b1 (Sheth et al. 2013; Desjacques et al. 2018). Similar to other studies (Lazeyras & Schmidt 2018; Zennaro et al. 2022), we find a strong disagreement between our measurements and the Modi et al. (2017) measurements. Furthermore, we show as a dashed line in Figure 12 the co-evolution relation of the tidal bias that was presented in Zennaro et al. (2022). This relation has been measured at very small scales k ∼ 0.7 h Mpc−1 and is not really expected to accurately represent the large-scale tidal bias. However, it is interesting to see that the large-scale tidal response that we measure here is so much weaker than the response that can be measured at small scales. Next, we made a comparison with the excursion set prediction from Sheth et al. (2013) parameterised as $b_1^{\mathrm{E}}=0.524-0.547b_1^{\mathrm{E}}+0.046{(b_1^{\mathrm{E}})}^2$. This prediction seems to be in better agreement with our measurements, since it generally predicts smaller amplitudes for the tidal bias and a slight negative tendency. However, the predicted upturn in the function at b1 ≳ 2 does not seem to be present in any of our measurements.

Furthermore, we compared our measurements with the Lagrangian tidal bias inferred by Abidi & Baldauf (2018). For this, we use their fitted b_K²(M) relation and map masses to b₁ values through the Tinker et al. (2010) bias relation evaluated for their simulated cosmology. Our measurements are in good qualitative and quantitative agreement with the results from Abidi & Baldauf (2018) when compared at the same value of b₁.

Finally, we compared our data with the measurements from Lazeyras & Schmidt (2018), again transformed from Eulerian to Lagrangian parameters, indicated via the error bars in Figure 12. Our measurements agree roughly with that study, in that we favor a very small tidal bias with a preferentially slightly negative amplitude. We suggest therefore, that high Lagrangian tidal biases |b_K²|> 0.5 for mass-selected halos are strongly disfavored for lowly biased objects b₁ ≪ 2.

However, we also note that there are some statistically significant differences at b₁ ∼ 0, where we find tidal bias consistent with zero, whereas the data of Lazeyras & Schmidt (2018) clearly favor a negative value of b_K² ∼ −0.2. Given our agreement with Abidi & Baldauf (2018), the fact that past erroneous measurements have rather predicted too large absolute values of the tidal bias rather than too small ones, and that the three lowest b₁ data points of Lazeyras & Schmidt (2018) do not seem really consistent with the overall behaviour of the curve, we suggest that our measurements are more reliable here and that the Lagrangian tidal bias becomes indeed very small |b_K²|≪0.1 for b₁ ≲ 0.5.

We conclude that overall the amplitude of the Lagrangian b_K² is quite low for mass selected halos. Given the challenges of the measurement and the history of errors, it would be desirable to have a clear measurement that does not suffer from the usual scale dependencies and large uncertainties. This could be achieved through anisotropic separate-universe simulations as they have been presented in Stücker et al. (2021), Masaki et al. (2020), Akitsu et al. (2021). We checked whether the simulations that were presented in Stücker et al. (2021) are sufficient to measure the tidal bias at the level that we observe here. The simulation parameter choices were not optimal for measuring it, since those simulations do not have sufficient statistics of high mass objects – which are the objects that would have a significant tidal bias according to Equation (151). However, a reliable measurement could be obtained by running with larger tidal fields that are tuned to measure the tidal bias at higher redshift, where most objects have a larger b₁ and likely a larger b_K². However, that said, we find that the simulations have significant enough statistics to constraint an upper bound |b_K²|< 0.1 (roughly 1σ level) at low masses M_200b ∼ 10¹³ h⁻¹ M_⊙, consistent with our observations here.

5.2. Assembly bias in b_K²

While the tidal bias seems to be quite low for mass selected halos, it is not at all obvious that this should also be the case for galaxies. The formation and selection of galaxies may depend on other properties that are more sensitive to the tidal field, like the spin, the formation time and the halo concentration. For example, in tidal torque theory the spin of an object should depend on the Lagrangian tidal tensor (White 1984). Therefore, it is important to investigate whether there exists assembly bias with respect to such properties. Lazeyras et al. (2021) have shown that the tidal bias bK2E and the associated bK2(b1) co-evolution relation depend strongly on secondary halo properties. Here, we briefly checked whether we can reproduce these results.

We measure the spin of each halo in our simulation following Bullock et al. (2001) as

$$\begin{array}{cc}\hfill \lambda & =\frac{|\mathit{L}|}{\sqrt{2}{v}_{200c}{r}_{200c}},\hfill \end{array}$$(152)

where

$$\begin{array}{cc}\hfill \mathit{L}& ={\langle \mathit{r}\times \mathit{v}\rangle }_{\mathrm{particles}}\hfill \end{array}$$(153)

is the specific angular momentum averaged over all particles that are within the bound component of the main subhalo according to SUBFIND and r_200c is the radius within which the density is 200 times the critical density of the Universe and v_200c the circular velocity at that radius.

We then grouped halos in each mass bin by their spin into four quartiles q₁, q₂, q₃ and q₄, q₁ corresponding to the 25% of halos with the lowest spin and q₄ to the 25% with the highest spin. We then measured the bias parameters independently in each quartile, q₁. We show the resulting b_K²(b₁) relation in Figure 13. Here, we present only results for the second-order spatial estimator at a damping scale of k_d = 0.25 h Mpc⁻¹, but we have checked whether other choices lead to the same results.

The b_K² − b₁ relation exhibits a strong degree of assembly bias with halo spin. For example, at b₁ ∼ 1, the tidal bias of mass selected halos seems very small |b_K²|≪0.1, but after the selection on halo spin, values may range between −1 and 0.5. In individual mass bins, the spin selection increases both the value of b₁ and the value of b_K².

The sign, the b1 dependence, and the overall amplitude of the spin-assembly bias are all in good agreement with the measurements of Lazeyras et al. (2021). For example, at b1 = 2, we find assembly bias variations of ΔbK2 ∼ 2 (between low and high spins) consistent with their variations at $b_1^E=3$ of around ΔbK2E ∼ 2 and the assembly bias decreases towards lower b1 and increases towards higher b1. Therefore, our measurements provide a validation of the results of that study, but it also shows that our method allows us to reliably measure the tidal bias.

While we do not present it here, we have also tested whether we can reproduce the concentration assembly bias that was measured by Lazeyras et al. (2021). We have found that if we also define concentration through v_max/v_200c, then we find an assembly bias of a similar amplitude. However, if concentrations are defined through fits of the halo density profile, the assembly bias is smaller (possibly consistent with 0) and prefers the opposite sign. The concentration assembly bias seems therefore to be a bit more uncertain than that of the spin.

We conclude that even if b_K² seems quite small for mass selected halos, the tidal bias may not so easily be neglected for tracers that follow more general selection criteria. Values of the Lagrangian b_K² of an order of unity (and slightly larger) are a plausible possibility for realistic galaxy catalogues.

5.3. Biases of third derivative terms

In Figure 14, we show the measurements that we obtain by evaluating the bias estimators from Equations (140) and (141) and where we again divide by a factor of 2 to get b_(∇δ)² = b_J3 − 3/2, in accordance with the literature convention.

The value of b_(∇δ)² appears to depend strongly on the damping scale and we cannot therefore claim a reliable measurement here. However, it seems mostly b_(∇δ)² < 0 which seems consistent with the idea that halos should preferentially form at peaks where density gradients should be small. However, at very large masses, b_(∇δ)² seems to become positive. We are not quite sure how to interpret this behaviour, but it might just be an artefact of the uncorrected scale dependencies. It is worth noting that the negative amplitude of b_(∇δ)² is quite significant b_(∇δ)² ∼ 10 h⁻² Mpc² and even grows towards larger scales. We therefore discuss whether this (usually neglected) bias term is worth including under some circumstances, as described late in this paper. We note that b_(∇δ)² was also measured by Biagetti et al. (2014). However, since we have not properly corrected the scale dependence here, we abstain from a comparison.

On the other hand, b_J3 ≡ 3 appears to be quite small, |b_J3 ≡ 3/2|≲5h⁻²Mpc², for everything except the very high mass end. Therefore, the contributions of b_J3 ≡ 3 to the bias expansion should be significantly lower than b_(∇δ)².

6. Relevance of bias parameters

We aim to estimate the relevance of different contributions to the bias expansion. If we have a term of the form $\frac{1}{n!}{b}_X·X$ of the re-normalised bias expansion, then we can estimate its relative relevance by comparing the typical difference in the predicted galaxy density field between a bias function, f_X, which includes the term versus an expansion, f₀, which neglects it, namely:

$$\begin{array}{cc}\hfill {\sigma }_{ϵ,X}^2& =\langle {(f_X-f_0)}^2\rangle \hfill \\ \hfill & ={\left(\frac{b_X}{n!}\right)}^2\langle X^2\rangle ,\hfill \end{array}$$(154)

where the value of ⟨X²⟩ depends on the adopted damping scale and where σ_ϵ, X is a dimensionless quantity that can be compared across different bias terms. For example, for the contribution of b₃, we have

$$\begin{array}{cc}\hfill {\sigma }_{ϵ,{\delta }^3}^2& ={\left(\frac{b_3}{6}\right)}^2\langle {({\delta }^3-3\delta {\sigma }^2)}^2\rangle .\hfill \end{array}$$(155)

The value of σ_ϵ, X quantifies the typical error that is made at the field level (in Lagrangian space) when neglecting the corresponding term.

We chose two halo sets M ∼ 10¹⁴ h⁻¹ M_⊙ and M ∼ 10¹⁵ h⁻¹ M_⊙ as examples and measured all their bias parameters with a damping scale of k_d = 0.2 h Mpc⁻¹ and used the second-order spatial estimators for all parameters (except b_(∇δ)² and b_3 ≡ 3). We list the corresponding parameters in Table 2. The first set represents a moderately biased population, whereas the second set represents an example of strongly biased tracers. We then kepy these parameters fixed and evaluate ⟨X²⟩ at many different damping scales to obtain σ_ϵ, X, which we show in Figure 15.

For the moderately biased tracers, the relevance of parameters decreases reasonably between different orders. However, this statement depends strongly on the considered damping scale. For instance, at k_d = 0.4 h Mpc⁻¹ many different terms reach similar amplitudes; whereas at k_d = 0.05 h Mpc⁻¹, there is a strong order of relevance among different terms. If we aim to fit a bias expansion at the field level at k_d = 0.15 h Mpc⁻¹ (marked as a vertical dotted line), then we would find that for the moderately biased set of b₁ > b₂ > b₃ > b₄, all of them lie within a factor of 2 of each other and for the strongly biased set, we even see their order has changed, namely, to: b₂ > b₁ > b₄ > b₃ and all of these lie within a factor of 3 of each other. This is certainly very problematic, since it does not at all seem clear that higher order terms are smaller than lower order terms here, which casts doubts on the convergence of the expansion. Of course this problem can be controlled by simply going to larger smoothing scales, but at the cost of a reduced constraining power. However, another possibility is to make use of the advantages of the cumulant bias parameters. With these parameters, significant parts of higher order terms get already absorbed by lower order parameters. For example, considering the first three terms of the bias expansion in terms of cumulant parameters:

$$\begin{array}{cc}\hfill f& =1+{\beta }_1\delta +\frac{1}{2}({\beta }_2+{\beta }_1^2)({\delta }^2-{\sigma }^2)\hfill \\ \hfill & +\frac{1}{6}({\beta }_3+3{\beta }_1{\beta }_2+{\beta }_1^3)({\delta }^3-\delta {\sigma }^2)+\cdots \hfill \end{array}$$(156)

Here β₁ and β₂ contribute also to the third term and to infinitely many higher order terms. If we write out the polynomial to high orders, but set β₃ = 0, then the leading order error we would make is

$$\begin{array}{cc}\hfill {\sigma }_{ϵ,{\beta }_3{\delta }^3}^2& ={\left(\frac{{\beta }_3}{6}\right)}^2\langle {({\delta }^3-3\delta {\sigma }^2)}^2\rangle ,\hfill \end{array}$$(157)

which is identical to the error term from the canonical bias terms, but with a different fore-factor. Error estimates of these types are indicated as dashed lines in Figure 15. Strikingly, these are significantly smaller than the canonical terms and it seems now plausibly to reach good convergence at k_d = 0.15 h Mpc⁻¹. For instance, to reach σ_ϵ ≤ 0.1 for the strongly biased set, we would only need to consider the terms (β₁, β₂, β_L).

It is noteworthy that the contribution of b_(∇δ)² is not much smaller than that of b_L so that it might make sense to include this parameter more commonly in the bias expansion at scales k_d ≳ 0.1 h Mpc⁻¹. However, we remind the reader that there are some other possibly important parameters here that we have not considered, so that we cannot present a comprehensive picture of the relevance of bias parameters.

We conclude that it may have significant advantages to phrase the bias expansion in terms of the cumulant parameters, β_n, instead of the canonical bias parameters. This is especially true when considering highly biased tracers or small smoothing scales, where the rate of convergence of the canonical expansion is excruciatingly slow.

7. Conclusions

In this article, we present a method for measuring bias parameters through moments of the galaxy environment distribution, p(δ|g). We have shown that such estimators can be derived both for scalar bias parameters such as b₁ and b_L and for tensorial ones such as b_K² and b_(∇δ)². Furthermore, we demonstrate that they can be made independent of scale if we consider spatial corrections at various orders. We have verified the reliability of these estimators by recovering well established relations between known parameters. Additionally, we used our estimators to measure terms that have been unknown thus far, such as b_3 ≡ 3 and b_(∇δ)², where the latter appears to be particularly large so that it might be worth including it more commonly in the bias expansion.

The main benefits of the method are its simplicity and generality. We are only required to evaluate the linear field of a simulation at the Lagrangian locations of a set of tracers (e.g. given by the most bound id in a halo catalogue) and to take simple expectation values. This measurement uses a minimal number of discretisation steps, making it numerically robust and easy to understand. While a similar approach was previously used by Paranjape et al. (2013a,b), we have shown that it is possible to incorporate spatial corrections at any order to get accurate reconstructions of large-scale bias parameters, even at fairly large damping scales, k_d. In contrast to measurement approaches that depend on forward models, for instance, of power spectra or correlation functions, measured bias parameters are significantly more independent of each other in the method at hand. For example, the measurement of b₁ is completely independent of the assumed order of the expansion in powers of δ. This makes it generally unlikely that the parameters can numerically compensate for the absence of other terms. Finally, we note that all bias parameters can be evaluated at a low computational cost with a single simulation and that it is possible to use arbitrary filtering methods with the slight alterations explained in Appendix A.

However, to avoid confusion, we want to point out a few limitations of our method for estimating the biases. First of all, we note that our method assumes that the background distribution is known analytically. Therefore, it does not immediately translate to situations where this is not the case. For example, at third order in the bias expansion, time derivatives of the tidal tensor appear and these depend on the second Lagrangian perturbation theory potential, ϕ⁽²⁾, which is not a Gaussian random field. Therefore, the formalism can only be translated to the associated bias parameters if the distribution of ϕ⁽²⁾ and its derivatives can be written (or at least adequately approximated) analytically. Furthermore, the formalism cannot directly be applied to the Eulerian galaxy environment distribution p^E(δ|g), since the non-linear Eulerian matter distribution p^E(δ) does also not follow Gaussian statistics.

Beyond the novel bias estimators, we have also newly introduced the cumulant bias parameters

$$\begin{array}{cc}\hfill {\beta }_n& ={\frac{{\partial }^n}{\partial {\delta }_0^n}log\left(\frac{n_g({\delta }_0)}{n_{g,0}}\right)|}_{{\delta }_0=0}.\hfill \end{array}$$(158)

These cumulant bias parameters are proportional to the cumulants of the galaxy environement distribution⁸ and they are related to the canonical bias parameters exactly in the same way as cumulants relate to moments. When measuring these parameters, we found them to be more independent from each other than canonical bias parameters. Furthermore, we have shown that for halos cumulant biases of the order of n ≥ 3 are very close to zero. This has several intriguing consequences:

(1) The galaxy environment distribution and the bias function of halos are very well approximated by a Gaussian. This motivates the use of a Gaussian bias model for the order of n = 2, as well as alternative approaches, for instance, expanding around a Gaussian bias model for orders of n > 2. We will discuss the Gaussian bias model and its desirable properties in detail in a companion paper (Stücker et al. 2025).

(2) The co-evolution relations of halos of the order of n ≥ 3 are equivalent to β_n = 0. Therefore, we may predict such relations between canonical bias parameters at very high orders by combining β_n = 0 with the mapping between canonical and cumulant bias parameters. However, this also suggests that such co-evolution relations may not carry much physical significance, but instead they are the artefacts of a suboptimal basis choice. An expansion around a Gaussian model would naturally incorporate all of these relations at any order.

(3) The bias expansion may have significantly improved convergence when phrased in terms of cumulants instead of canonical bias parameters. At the field level, the convergence of the canonical expansion seems already questionable at scales of k_d ≳ 0.15 h Mpc⁻¹ since the parameters may easily rise as their order rises. However, when phrased in terms of cumulant biases, we expect neglected higher order terms to be significantly less important, as they are already well captured through lower order parameters.

Overall, we expect cumulant biases to be most beneficial whenever the canonical expansion shows poor convergence, for instance, at high masses, at late times, or at small smoothing scales.

---

Acknowledgments

The authors thank Simon White, Oliver Hahn and Fabian Schmidt for helpful discussions and comments to the draft. JS thanks Oliver Philcox for helpful discussions. We acknowledge funding from the Spanish Ministry of Science and Innovation through grant number PID2021-128338NB-I00. RV acknowledges the support of the Juan de la Cierva fellowship (FJC2021-048002-I). MPI is supported by STFC consolidated grant no. RA5496.

References

Appendix A: Estimators with different filters

In Section 2, we have shown derivations under the assumption of a sharp k-space filter. For this choice all considered bias functions and bias parameters relate trivially to the separate-universe parameters that are commonly referred to in the literature. However, in general other filtering kernels could also be chosen. For example if we define the smoothed density as

$$\begin{array}{cc}\hfill {\delta }_l(k)& =W(k)\delta (k),\hfill \end{array}$$(A.1)

where δ_l is the smoothed density and δ the unsmoothed one. If the filter W is for example chosen to be a Gaussian kernel,

$$\begin{array}{c}\hfill W_{\mathrm{g}}(k)=exp(-\frac{k^2}{2k_{\mathrm{d}}^2}),\end{array}$$(A.2)

instead of a sharp k filter, then some considerations have to be adapted. In principle, all of the formula from Section 2 could still be applied, but they would refer to a different set of bias parameters that describe the response to perturbations of a different shape. Rather than introducing a new set of bias parameters, it is more convenient to directly write the commonly used parameters in terms of the new basis. Therefore, we show here how to directly measure the separate-universe bias parameters from the moments of the galaxy environment distribution that was obtained with a density field that was filtered with any kernel.

First of all, we consider how a bias function on very large scales f_l(δ_l) should relate to the bias function on some much smaller scale f_s(δ_s). If the PBS is still valid on the smaller scale, then the probability of forming a galaxy, given knowledge of δ_s and δ_l depends only on the value of δ_s:

$$\begin{array}{cc}\hfill p(\mathrm{g}|{\delta }_s,{\delta }_l)& =p(\mathrm{g}|{\delta }_s).\hfill \end{array}$$(A.3)

Under this assumption, we can re-express the large scale bias function through the small scale function marginalised over the conditional probability of the small scale densities:

$$\begin{array}{cc}\hfill f_l({\delta }_l)=\frac{p(\mathrm{g}|{\delta }_l)}{p(\mathrm{g})}& =\int \frac{p(\mathrm{g}|{\delta }_s)}{p(g)}p({\delta }_s|{\delta }_l)d{\delta }_s\hfill \\ \hfill & =\int f_s({\delta }_s)p({\delta }_s|{\delta }_l)d{\delta }_s.\hfill \end{array}$$(A.4)

Since, the joint distribution of δ_s and δ_l is a multivariate Gaussian, the conditional distribution corresponds to a multivariate Gaussian that is conditioned on the value of one variable, which itself is a Gaussian with modified mean and covariance:

$$\begin{array}{cc}\hfill p({\delta }_s|{\delta }_l)& =N({\delta }_s,{\mu }_{\ast },{\sigma }_{\ast }),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\mu }_{\ast }& =\alpha {\delta }_l,\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill {\sigma }_{\ast }^2& ={\sigma }_{\mathit{ss}}^2-{\alpha }^2{\sigma }_{\mathit{ll}}^2,\hfill \end{array}$$(A.7)

where α = σ_ls²/σ_ll² corresponds to the correlation between the smoothed large-scale and small-scale density fields, σ_ll² = ⟨δ_lδ_l⟩, σ_ls² = ⟨δ_lδ_s⟩ and σ_ss = ⟨δ_sδ_s⟩.

If we define the re-normalised bias function, F_α, through the bias function f_l as in Section 2, then we find

$$\begin{array}{cc}\hfill F_{\alpha }({\delta }_0)& =\int f_l({\delta }_l)p({\delta }_l-{\delta }_0)\mathrm{d}{\delta }_l\hfill \\ \hfill & =\int f_s({\delta }_s)p_{\mathit{ss}}({\delta }_{\mathit{ss}}-\alpha {\delta }_0)d{\delta }_s,\hfill \end{array}$$(A.8)

where p_ss is simply the background distribution of the small scale densities. This can be shown either by evaluating the integrals (corresponding to convolutions) explicitly or by considering the large scale limit of f_l where σ_* → σ_ss.

We can now relate the re-normalised bias functions obtained with different filtering methods, since we may assume that f_s is independent of the filter employed on large scales. Therefore,

$$\begin{array}{c}\hfill F_{{\alpha }_1}\left(\frac{{\delta }_0}{{\alpha }_1}\right)=F_{{\alpha }_2}\left(\frac{{\delta }_0}{{\alpha }_2}\right),\end{array}$$(A.9)

where we have indicated the usage of different filters through different α in the index. Therefore, different filters simply lead to a rescaling of the re-normalised bias function. The value of α depends on the filter and the smoothing scale. If we assume that δ_s was obtained with some small scale filter W_s and δ_l with W_l, we have

$$\begin{array}{cc}\hfill {\sigma }_{\mathit{ls}}^2& =\langle \int \delta ({\mathit{k}}_1)W_s(k_1){\mathrm{d}}^3{\mathit{k}}_1\int \delta ({\mathit{k}}_2)W_l(k_2){\mathrm{d}}^3{\mathit{k}}_2\rangle \hfill \\ \hfill & =\int \int \langle \delta ({\mathit{k}}_1)\delta ({\mathit{k}}_2)\rangle W_s(k_1)W_l(k_2){\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2\hfill \\ \hfill & =\int \int P(k_1){\delta }_D({\mathit{k}}_1-{\mathit{k}}_2)W_s(k_1)W_l(k_2){\mathrm{d}}^3{\mathit{k}}_1{\mathrm{d}}^3{\mathit{k}}_2\hfill \\ \hfill & =\int P(k)W_s(k)W_l(k){\mathrm{d}}^3\mathit{k},\hfill \end{array}$$(A.10)

$$\begin{array}{cc}& \approx \int P(k)W_l(k){\mathrm{d}}^3\mathit{k}\mathbf{,}\hfill \end{array}$$(A.11)

where the last approximation assumes that W_s is very close to 1 over the support of W_l, because W_l corresponds to a much larger scale smoothing than W_s. Furthermore, we have

$$\begin{array}{cc}\hfill {\sigma }_{\mathit{ll}}^2& =\int P(k)W_l^2(k){\mathrm{d}}^3\mathit{k}\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill \alpha & =\frac{\int P(k)W_l(k){\mathrm{d}}^3\mathit{k}}{\int P(k)W_l^2(k){\mathrm{d}}^3\mathit{k}}.\hfill \end{array}$$(A.13)

It is easy to see that α = 1 for any sharp Fourier space filter, independent of the filtering scale. However, for other filters it is important to note that α may be generally larger than 1, for example α ∼ 2 for a Gaussian filter at scales of k_d ∼ 0.1 h Mpc⁻¹. Furthermore, α depends generally on the filtering-scale. The re-normalised bias function does therefore only correspond to a well-defined limit at arbitrary large scales if α is constant (as for the sharp k filter).

Since the case α = 1 seems clearly to be the simplest and leads to a sensibly scale-independent re-normalisation, we consider it as the reference case and simply label its re-normalised bias function F, and we have focused in the main-text on the discussion of this case. However, we can still use the relations derived in the main-text to measure the properties of F directly through measuring aspects of some filter-dependent function F_α.

We directly show this for the general multivariate case here by adapting the considerations from Section 2.7. For the multivariate case the re-normalised bias functions are relate as

$$\begin{array}{cc}\hfill F_{\mathsf{A}}({\mathit{x}}_0)& =F(\mathsf{A}{\mathit{x}}_0),\hfill \end{array}$$(A.14)

where

$$\begin{array}{cc}\hfill \mathsf{A}& ={\mathsf{C}}_{\mathit{ls}}{\mathsf{C}}_{\mathit{ll}}^{-1},\hfill \end{array}$$(A.15)

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathit{ls}}& =\langle {\mathit{x}}_l\otimes {\mathit{x}}_s\rangle ,\hfill \end{array}$$(A.16)

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathit{ll}}& =\langle {\mathit{x}}_l\otimes {\mathit{x}}_l\rangle ,\hfill \end{array}$$(A.17)

where C_ls may more conveniently be measured as ⟨x_m⊗x_m⟩, where x_m is evaluated in the linear fields filtered with $\sqrt{W_l}$.

The cumulant-generating function of the galaxy environment distribution with δ_l is again given through Equation (64), just with F_A instead of F:

$$\begin{array}{cc}\hfill K_{\mathsf{A}}(\mathit{t})& =\frac{1}{2}{\mathit{t}}^T{\mathsf{C}}_{\mathit{ll}}\mathit{t}+log\left(F_{\mathsf{A}}({\mathsf{C}}_{\mathit{ll}}\mathit{t})\right)\hfill \\ \hfill & =\frac{1}{2}{\mathit{t}}^T{\mathsf{C}}_{\mathit{ll}}\mathit{t}+log(F({\mathsf{C}}_{\mathit{ls}}\mathit{t})),\hfill \end{array}$$(A.18)

For measuring the cumulant bias parameters, we may then define a more convenient variable,

$$\begin{array}{cc}\hfill {\mathit{u}}_{\ast }& ={\mathsf{C}}_{\mathit{ls}}{\mathit{x}}_0,\hfill \end{array}$$(A.19)

which has the cumulant generating function,

$$\begin{array}{cc}\hfill K_{{\mathit{u}}_{\ast }}(\mathit{t})& =\frac{1}{2}{\mathit{t}}^T{\mathsf{C}}_{\mathit{ls}}^{-1}{\mathsf{C}}_{\mathit{ll}}{\mathsf{C}}_{\mathit{ls}}^{-1}\mathit{t}+log(F(\mathit{t})).\hfill \end{array}$$(A.20)

Thus, the cumulant bias parameters can be inferred by differentiating, leading to:

$$\begin{array}{cc}\hfill {\mathit{\beta }}_{ijk...}& =\{\begin{array}{cc}{\kappa }_{u_{\ast },ijk...}\hfill & ifi+j+k+...\ne 2,\hfill \\ {\kappa }_{u_{\ast },ijk...}-C_{\ast ,ab}^{-1}\hfill & ifi+j+k+...=2,\hfill \end{array}\hfill \end{array}$$(A.21)

where ${\mathsf{C}}_{\ast }^{-1}={\mathsf{C}}_{\mathit{ls}}^{-1}{\mathsf{C}}_{\mathit{ll}}{\mathsf{C}}_{\mathit{ls}}^{-1}$. We can therefore measure the separate-universe bias parameters with an arbitrarily filtered distribution in a very similar manner to the one using a sharp k-space filter. It is worth noting that the inferred relations show also an important property: If β_n vanish for n > 2, which automatically implies that the corresponding cumulants vanish for any filtered galaxy environment distribution. Therefore, if the bias function is Gaussian for one filter, then it will also appear Gaussian for any other filter.

Appendix B: Isotropic tensors

B.1. Tensors with three or more symmetry groups

In the main text, we limit our discussion to an explanation of how to obtain basis tensors for tensor spaces with at most two symmetry groups. However, for higher order biases, for instance, the equivalent of b₃ the isotropic tensor decomposition requires consideration of additional symmetries. We define 𝕍₂₂₂ as the space of tensors that are symmetric under permutations of the first two, second two and third two indices. Furthermore, we define 𝕍₍₂₂₂₎ as the space of tensors that are additionally invariant to exchanges of these symmetry groups. These spaces have different number of basis tensors,

$$\begin{array}{cc}\hfill {\mathbb{V}}_{222}& =\mathrm{span}(\{{\mathsf{I}}_{222},{\mathsf{I}}_{2-2-2-},{\mathsf{I}}_{2=22},{\mathsf{I}}_{22=2},{\mathsf{I}}_{222=}\}),\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill {\mathbb{V}}_{(222)}& =\mathrm{span}(\{{\mathsf{I}}_{222},{\mathsf{I}}_{2-2-2-},{\mathsf{I}}_{(2=22)}\}),\hfill \end{array}$$(B.2)

where

$$\begin{array}{cc}\hfill {\mathsf{I}}_{(2=22)}& =\frac{1}{3}({\mathsf{I}}_{2=22}+{\mathsf{I}}_{22=2}+{\mathsf{I}}_{222=})\hfill \end{array}$$(B.3)

is additionally symmetrised between exchanges between the symmetry groups and the other tensors are given by

$$\begin{array}{cc}\hfill I_{222,ijklmn}& =S_{222}({\delta }_{\mathit{ij}}{\delta }_{\mathit{kl}}{\delta }_{\mathit{mn}}),\hfill \end{array}$$(B.4)

$$\begin{array}{cc}\hfill I_{2-2-2-,ijklmn}& =S_{222}({\delta }_{\mathit{ik}}{\delta }_{\mathit{lm}}{\delta }_{\mathit{nj}}),\hfill \end{array}$$(B.5)

$$\begin{array}{cc}\hfill I_{2=22,ijklmn}& =S_{222}({\delta }_{\mathit{ik}}{\delta }_{\mathit{jl}}{\delta }_{\mathit{mn}}),\hfill \end{array}$$(B.6)

$$\begin{array}{cc}\hfill I_{22=2,ijklmn}& =S_{222}({\delta }_{\mathit{ij}}{\delta }_{\mathit{km}}{\delta }_{\mathit{ln}}),\hfill \end{array}$$(B.7)

$$\begin{array}{cc}\hfill I_{222=,ijklmn}& =S_{222}({\delta }_{\mathit{im}}{\delta }_{\mathit{kl}}{\delta }_{\mathit{jn}}).\hfill \end{array}$$(B.8)

The basis tensors of 𝕍₍₂₂₂₎ are visualised through a diagram in Figure B.1. The third derivative tensor of that would appear in the bias expansion can be decomposed in this space

$$\begin{array}{c}\hfill {\frac{{\partial }^{3}F}{\partial \mathsf{T}\partial \mathsf{T}\partial \mathsf{T}}|}_{\mathsf{T}=\mathsf{0}}\in {\mathbb{V}}_{(222)}.\end{array}$$(B.9)

B.2. Covariance of potential derivatives

To find an explicit form for the distribution of the derivatives of the tidal tensor and other higher order derivatives of the potential, we need to find their covariance matrix and invert it.

For this we first consider the covariance between any derivatives of the potential. This is given by

$$\begin{array}{cc}\hfill C_{ij...,ab...}& =⟨{\varphi }_{ij...}{\varphi }_{ab...}⟩\hfill \\ \hfill & =\frac{1}{{(2\pi )}^6}\langle \int \int \frac{(i{k}_i)(i{k}_j)...(i{k}_a)(i{k}_b)...}{k^4}\hfill \\ \hfill & e^{-i({\mathit{k}}_{\mathbf{1}}-{\mathit{k}}_{\mathbf{2}})\mathit{x}}{\delta }_k^{\ast }({\mathit{k}}_1){\delta }_k({\mathit{k}}_2)\mathrm{d}{\mathit{k}}_1\mathrm{d}{\mathit{k}}_2\rangle ,\hfill \end{array}$$(B.10)

where δ_k is the Fourier transform of density field δ and we have used that δ is a real field so that δ_k(k) = δ_k^*(−k). Using the definition of the power spectrum,

$$\begin{array}{cc}\hfill ⟨{\delta }_k^{\ast }({\mathit{k}}_1){\delta }_k({\mathit{k}}_2)⟩& =(2{\pi }^3)P(\mathit{k}){\delta }_{\mathrm{D}}^3({\mathit{k}}_1-{\mathit{k}}_2),\hfill \end{array}$$(B.11)

where δ_D is the Dirac delta function, we can simplify

$$\begin{array}{cc}\hfill C_{ij...,ab...}& =\frac{1}{{(2\pi )}^3}\int \frac{(i{k}_i)(i{k}_j)...(i{k}_a)(i{k}_b)...}{k^4}P(k){\mathrm{d}}^3\mathit{k}\mathbf{.}\hfill \end{array}$$(B.12)

We can see that this expression does not depend on the order of the indices, nor whether they stem from the first or the second potential term. Therefore, we just count the total number n₁ of occurrences of k_x, n₂ of k_y and n₃ of k_z. Without loss of generality we can order the coordinates, such that n₁ ≥ n₂ ≥ n₃. Furthermore, we can switch to spherical coordinates and replace the cosine of the azimutal variable:

$$\begin{array}{cc}\hfill \mathit{k}& =\left(\begin{array}{c}k\mu \\ k{(1-\mu )}^{2}cos\varphi \\ k{(1-\mu )}^2{sin\varphi )}^T\end{array}\right),\hfill \end{array}$$(B.13)

so that we have

$$\begin{array}{cc}\hfill \mathrm{\Sigma }(ij...ab)& =\frac{{\sigma }_M^2}{4\pi }{\int }_{-1}^1{\mu }^{n_1}\mathrm{d}\mu {\int }_0^{2\pi }{cos}^{n_2}(\varphi ){sin}^{n_3}(\varphi )\mathrm{d}\varphi ,\hfill \end{array}$$(B.14)

$$\begin{array}{cc}\hfill {\sigma }_M^2& =\frac{1}{2{\pi }^2}\int P(k)k^{2M}\mathrm{d}k,\hfill \end{array}$$(B.15)

$$\begin{array}{cc}\hfill M& =\frac{n_1+n_2+n_3-4}{2}.\hfill \end{array}$$(B.16)

These integrals can easily be evaluated by hand or with a computer algebra system. We list a few examples of the covariances in Table B.1 We note that only terms where each n₁, n₂, and n₃ are even are non-zero. We note that all uneven terms are zero, so that, for instance, the two-second derivatives of the potential have zero covariance with the third derivatives, but a non-zero covariance with the fourth derivatives.

B.3. The distribution of the tidal tensor

From Equation (B.12) it can be seen that the covariance tensor of any derivatives of the potential has to be an isotropic tensor with full symmetry in all indices. Therefore, the covariance matrix of the tidal tensor is C_T ∈ 𝕍₄ so that it has to be proportional to I₄ with a coefficient that we can read off from table B.1:

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}}& =\langle \mathsf{T}\otimes \mathsf{T}\rangle \hfill \\ \hfill & =\frac{{\sigma }^2}{5}{\mathsf{J}}_4.\hfill \end{array}$$(B.17)

This tensor is not invertible, since, for instance, the T₁₂ component is perfectly degenerate with the T₂₁ component. However, in such cases where the distribution of T can still be inferred by constructing a generalised inverse ${\mathsf{C}}_{\mathsf{T}}^{+}$, which has the following property:

$$\begin{array}{cc}\hfill \mathsf{C}{\mathsf{C}}_{\mathsf{T}}^{+}\mathsf{C}& ={\mathsf{C}}_{\mathsf{T}}.\hfill \end{array}$$(B.18)

It then holds that

$$\begin{array}{cc}\hfill p(\mathsf{T})& =Nexp(-\frac{1}{2}{\mathsf{T}}^T{\mathsf{C}}_{\mathsf{T}}^{+}\mathsf{T}).\hfill \end{array}$$(B.19)

To find the generalised inverse we may assume that C_T⁺ ∈ 𝕍₂₂, so that we can make the Ansatz

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}}^{+}& =A_{22}{\mathsf{J}}_{22}+A_{2=2}{\mathsf{J}}_{2=2}.\hfill \end{array}$$(B.20)

We can use the decomposition ${\mathsf{J}}_4=\frac{5}{9}{\mathsf{J}}_{22}+\frac{2}{3}{\mathsf{J}}_{2=2}$ (determined computationally) and evaluate Equation (B.18) term by term:

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}}{\mathsf{C}}_{\mathsf{T}}^{+}& =A_{22}\frac{{\sigma }^2}{9}({\mathsf{J}}_{22}\stackrel{(2)}{·}{\mathsf{J}}_{22})+A_{2=2}\frac{2{\sigma }^2}{15}({\mathsf{J}}_{2=2}\stackrel{(2)}{·}{\mathsf{J}}_{2=2}),\hfill \end{array}$$(B.21)

$$\begin{array}{cc}& =3A_{22}\frac{{\sigma }^2}{9}{\mathsf{J}}_{22}+A_{2=2}\frac{2{\sigma }^2}{15}{\mathsf{J}}_{2=2},\hfill \end{array}$$(B.22)

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}}{\mathsf{C}}_{\mathsf{T}}^{+}{\mathsf{C}}_{\mathsf{T}}& =9A_{22}{\left(\frac{{\sigma }^2}{9}\right)}^2{\mathsf{J}}_{22}+A_{2=2}{\left(\frac{2{\sigma }^2}{15}\right)}^2{\mathsf{J}}_{2=2},\hfill \end{array}$$(B.23)

$$\begin{array}{cc}& \equiv {\mathsf{C}}_{\mathsf{T}}=\frac{{\sigma }^2}{9}{\mathsf{J}}_{22}+\frac{2{\sigma }^2}{15}{\mathsf{J}}_{2=2},\hfill \end{array}$$(B.24)

where we have evaluated the appearing tensor products computationally. By identifying coefficients, we find $A_{22}=\frac{1}{{\sigma }^2}$ and $A_{2=2}=\frac{15}{2{\sigma }^2}$ and therefore we have the generalised inverse of the covariance matrix:

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}}^{+}& =\frac{1}{{\sigma }^2}{\mathsf{J}}_{22}+\frac{15}{2{\sigma }^2}{\mathsf{J}}_{2=2}\hfill \end{array}$$(B.25)

and the distribution of the tidal tensor is given by

$$\begin{array}{cc}\hfill p(\mathsf{T})& =Nexp(-\frac{2{\mathsf{T}}^T{\mathsf{J}}_{22}\mathsf{T}+15{\mathsf{T}}^T{\mathsf{J}}_{2=2}\mathsf{T}}{4{\sigma }^2})\hfill \end{array}$$(B.26)

$$\begin{array}{cc}& =Nexp(-\frac{2{\delta }^2+15K^2}{4{\sigma }^2}).\hfill \end{array}$$(B.27)

B.4. The distribution of third derivatives

The third derivative tensor of the potential S has the covariance matrix

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{S}}& =\langle \mathsf{S}\otimes \mathsf{S}\rangle ,\hfill \end{array}$$(B.28)

$$\begin{array}{cc}& =\frac{{\sigma }_1^2}{7}{\mathsf{J}}_6,\hfill \end{array}$$(B.29)

$$\begin{array}{cc}& =\frac{3{\sigma }_1^2}{25}{\mathsf{J}}_{3-3}+\frac{2{\sigma }_1^2}{35}{\mathsf{J}}_{3\equiv 3}.\hfill \end{array}$$(B.30)

Analogously to the calculation in Section B.3, we find the generalised inverse

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{S}}^{+}& =\frac{3}{{\sigma }_1^2}{\mathsf{J}}_{3-3}+\frac{35}{2{\sigma }_1^2}{\mathsf{J}}_{3\equiv 3},\hfill \end{array}$$(B.31)

where we have again used numerical representations to evaluate products between isotropic tensors when evaluating Equation (B.18). Therefore, the distribution of S is given by

$$\begin{array}{cc}\hfill p(\mathsf{S})& =Nexp(-\frac{35\mathsf{S}{\mathsf{J}}_{3\equiv 3}\mathsf{S}-15\mathsf{S}{\mathsf{J}}_{3-3}\mathsf{S}}{4{\sigma }_1^2})\hfill \\ \hfill & =:Nexp(-\frac{35S_{3\equiv 3}-15S_{3-3}}{4{\sigma }_1^2}).\hfill \end{array}$$(B.32)

where in the last line we have defined some new scalar variables through the given contractions, where S_3 − 3 = (∇δ)².

B.5. Joint distribution of second and fourth derivatives

We can write the covariance matrix of fourth and second derivatives as a block matrix

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}\mathsf{R}}& =\left(\begin{array}{cc}\langle \mathsf{T}\otimes \mathsf{T}\rangle & \langle \mathsf{T}\otimes \mathsf{R}\rangle \\ \langle \mathsf{R}\otimes \mathsf{T}\rangle & \langle \mathsf{R}\otimes \mathsf{R}\rangle \end{array}\right),\hfill \end{array}$$(B.33)

$$\begin{array}{cc}\hfill & =\left(\begin{array}{cc}\frac{{\sigma }^2}{5}{\mathsf{J}}_4& \frac{{\sigma }_1^2}{7}{\mathsf{J}}_6\\ \frac{{\sigma }_1^2}{7}{\mathsf{J}}_6& \frac{{\sigma }_2^2}{9}{\mathsf{J}}_8\end{array}\right),\hfill \end{array}$$(B.34)

$$\begin{array}{cc}\hfill & =\left(\begin{array}{cc}\frac{{\sigma }_0^2}{9}{\mathsf{J}}_{22}+\frac{2{\sigma }_0^2}{15}{\mathsf{J}}_{2=2}& \frac{{\sigma }_1^2}{15}{\mathsf{J}}_{24}+\frac{4{\sigma }_1^2}{35}{\mathsf{J}}_{2=4}\\ \frac{{\sigma }_1^2}{15}{\mathsf{J}}_{42}+\frac{4{\sigma }_1^2}{35}{\mathsf{J}}_{4=2}& \frac{{\sigma }_2^2}{25}{\mathsf{J}}_{44}+\frac{24{\sigma }_2^2}{245}{\mathsf{J}}_{4=4}+\frac{8{\sigma }_2^2}{315}{\mathsf{J}}_{4==4}\end{array}\right).\hfill \end{array}$$(B.35)

For its inverse, we make the Ansatz

$$\begin{array}{cc}\hfill {\mathsf{C}}_{\mathsf{T}\mathsf{R}}^{+}& =\left(\begin{array}{cc}A{\mathsf{J}}_{22}+B{\mathsf{J}}_{2=2}& C{\mathsf{J}}_{24}+D{\mathsf{J}}_{2=4}\\ C{\mathsf{J}}_{42}+D{\mathsf{J}}_{4=2}& E{\mathsf{J}}_{44}+F{\mathsf{J}}_{4=4}+G{\mathsf{J}}_{4==4}\end{array}\right)\hfill \end{array}$$(B.36)

and we find the solution to its generalised inverse by solving

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathsf{T}\mathsf{R}}={\mathsf{C}}_{\mathsf{T}\mathsf{R}}{\mathsf{C}}_{\mathsf{T}\mathsf{R}}^{+}{\mathsf{C}}_{\mathsf{T}\mathsf{R}},\end{array}$$(B.37)

with a computer algebra system and find

$$\begin{array}{c}\hfill {\mathsf{C}}_{\mathsf{T}\mathsf{R}}^{+}=\left[\begin{array}{cc}\frac{15{\sigma }_2^2{\mathsf{J}}_{2=2}}{2{\sigma }_0^2{\sigma }_2^2-2{\sigma }_1^4}+\frac{{\sigma }_2^2{\mathsf{J}}_{22}}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}& -\frac{15{\sigma }_1^2{\mathsf{J}}_{2=4}}{2{\sigma }_0^2{\sigma }_2^2-2{\sigma }_1^4}-\frac{{\sigma }_1^2{\mathsf{J}}_{24}}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}\\ \frac{15{\sigma }_1^2{\mathsf{J}}_{4=2}}{2{\sigma }_0^2{\sigma }_2^2-2{\sigma }_1^4}-\frac{{\sigma }_1^2{\mathsf{J}}_{42}}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}& \frac{15{\sigma }_0^2{\mathsf{J}}_{4=4}}{2{\sigma }_0^2{\sigma }_2^2-2{\sigma }_1^4}+\frac{{\sigma }_0^2{\mathsf{J}}_{44}}{{\sigma }_0^2{\sigma }_2^2-{\sigma }_1^4}+\frac{315{\mathsf{J}}_{4==4}}{8{\sigma }_2^2}\end{array}\right].\end{array}$$(B.38)

The joint distribution of second and fourth derivatives of the potential can be written as

$$\begin{array}{cc}\hfill p(\mathsf{T},\mathsf{R})& =Nexp(-\frac{1}{2}\left(\begin{array}{cc}\mathsf{T}& \mathsf{R}\end{array}\right){\mathsf{C}}_{\mathsf{T}\mathsf{R}}^{+}\left(\begin{array}{c}\mathsf{T}\\ \mathsf{R}\end{array}\right)).\hfill \end{array}$$(B.39)

B.6. Tidal estimator with spatial correction

We find the tidal estimator with spatial corrections of the order of 2 as:

$$\begin{array}{cc}\hfill b_{{\mathsf{J}}_{\mathsf{2}=\mathsf{2}}}& ={\langle \frac{1}{p(\mathsf{T},\mathsf{R})}\frac{{\partial }^{2}p(\mathsf{T},\mathsf{R})}{\partial \mathsf{T}\partial \mathsf{T}}\stackrel{(4)}{·}\frac{{\mathsf{J}}_{\mathsf{2}=\mathsf{2}}}{\left||{\mathsf{J}}_{2=2}\right|{|}^2}\rangle }_{\mathrm{g}},\hfill \end{array}$$(B.40)

$$\begin{array}{cc}& =\frac{15}{4{\sigma }_{\ast }^8}{\langle 3K^2{\sigma }_2^4+6{\varphi }_{2=4}{\sigma }_1^2{\sigma }_2^2+3{\varphi }_{4=4}{\sigma }_1^4-2{\sigma }_2^2{\sigma }_{\ast }^4\rangle }_{\mathrm{g}},\hfill \end{array}$$(B.41)

where

$$\begin{array}{cc}\hfill {\varphi }_{2=4}& :=\mathsf{T}{J}_{2=4}\mathsf{R}=\sum _{\mathit{ij}}{\partial }_i{\partial }_j\varphi {\partial }_i{\partial }_j\delta -\frac{1}{3}\delta {\mathrm{\nabla }}^2\delta ,\hfill \end{array}$$(B.42)

$$\begin{array}{cc}\hfill {\varphi }_{4=4}& :=\mathsf{T}{J}_{4=4}\mathsf{R}=\sum _{\mathit{ij}}{({\partial }_i{\partial }_j\delta )}^2-\frac{1}{3}{({\mathrm{\nabla }}^2\delta )}^2.\hfill \end{array}$$(B.43)

While we do not list them here, we have verified that if we analogously evaluate the expressions for b_J2, b_J4, b_J22 and b_J24 that they are identical to the estimators for b₁, b_L, b₂, and b_δ, L, respectively, from Equations (49) - (53).

All Tables

All Figures