Galaxy redshift-space bispectrum: the Importance of Being Anisotropic
Davide Gualdi Licia Verde
astro-ph.CO

Detection of the strongest magnetic field in a sunspot light bridge 
    Michiel van Noort
    
	
	  NOMAD Laboratory, Fritz Haber Institute of the Max Planck Society, Berlin, Germany;    Machine Learning Group, Technische Universität Berlin, Germany;    Institute of Mathematics, Technische Universität Berlin, Germany;    Citrine Informatics, Redwood City, CA, USA    Present address: Department of Computer and Information Science, University of Konstanz, Germany.
 

===========================================================================================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION



Higher-order statistics represent the most direct tool to access the cosmological information encrypted in the non-Gaussian features of the matter density field. Three-point (3pt) statistics are the lowest order, and therefore the  simplest, tools enabling us to harvest this information. The most recent analyses for the Baryon Oscillation Spectroscopic Survey (BOSS) <cit.> spectroscopic data-set, part of the Sloan Digital Sky Survey (SDSS) <cit.>, have been performed for both the 3pt correlation function <cit.> and its Fourier transform, the bispectrum <cit.>.

Approximately four decades have passed since the first 3pt statistics studies <cit.>. However only twenty years ago a mature formalism was developed to map theoretical models to spectroscopic data-sets in redshift space <cit.>.

The window on non-linear scales was opened with the loop corrections inclusion <cit.>; higher-order terms become important at non-linear scales and are indeed necessary to study primordial non-Gaussianity using the 3pt statistics to lift the degeneracy with the other cosmological parameters <cit.>.



The general benefits of using the bispectrum in addition to the power spectrum have recently been studied in <cit.>. Originally, the bispectrum was thought to be useful in order to break the degenaracies between bias and cosmological parameters <cit.>. Lately, the bispectrum has been shown to be a key tool for reaching more specific goals, such as detecting General Relativity effects <cit.> and improving the neutrino masses  constraints <cit.>.


Nevertheless, all the analyses performed to date on spectroscopic data have been limited to the isotropic signal component, the so-called bispectrum "monopole" <cit.>. In analogy to the power spectrum, redshift space distortions (RSD) induce an  anisotropic signal with respect to the line of sight (see pioneering work by  <cit.>). This anisotropy can be expanded  in multipole moments, where the monopole is just the averaged signal in all directions around the line of sight. To date, multipoles beyond the monopole have  been used to constrain cosmological parameters only for the power spectrum <cit.>, with the state-of-the-art at the moment of writing being Refs. <cit.>. Measuring higher-order  multipoles of the redshift-space galaxy bispectrum is complicated in terms of estimators modelling <cit.>, even when simulated volumes with periodic boundary conditions are used. For galaxy surveys, these measurements are even more difficult because of the necessity of accounting for all observational effects (e.g., window function effect, varying average density number, varying line of sight, etc.). 
As in the case of the power spectrum, the anisotropic signal in the bispectrum due to RSD potentially carries key additional cosmological information.

This additional signal can be used not only to obtain tighter parameter constraints, but also as a self-consistency check for the theoretical model itself. In other words one could test whether the same set of cosmological parameters can produce an analytical template that fits well all the measured data vector components (power spectrum / bispectrum, isotropic / anisotropic).
Hence we have a strong motivation to investigate the potential of the   galaxy bispectrum multipoles.

The purpose of this work is to assess the added constraining power obtained by including higher-order multipoles of the bispectrum to the so far used large scale structure (LSS) statistics in Fourier space (power specrum multipoles and bispectrum monopole).
We study the improvements, in terms of tighter parameter constraints, in the 1- and 2-D 68% credible regions resulting from using the bispectrum multipoles. Fisher forecasts for the bispectrum multipoles, using a diagonal covariance matrix including only Gaussian terms, were realised in <cit.>. Here we go beyond this approximation, dedicating a significant effort in modelling the off-diagonal, non-Gaussian terms of the covariance matrix.


Since measurements of these new statistics (the anisotropic bispectrum modulation) either from mocks or data are at the moment extremely expensive and the necessary codes are not fully developed or publicly available yet, it is important to estimate first the potential gain these new statistics could offer. To do so however the full covariance matrix must be estimated in order to derive the cosmological parameters posterior distributions, for example through Markov Chain Monte Carlo methods (MCMC).  
We argue here that resorting to MCMC's  is indeed necessary to  explore and place constraints on parameters-space  where there are non-trivial degeneracies, such as the parameters usually constrained in a clustering analysis. Fisher forecast can not do this since are based on the first derivatives of the data-vector's theoretical model.


Note that the covariance must be estimated for the full "data vector" which comprises both the power spectrum and bispectrum multipoles;  this matrix cannot be approximated as diagonal. Traditionally, there are two approaches to estimate covariance matrices: from a large number of mocks or analytically. Adopting the first approach alone is not viable here, for the reason explained above that measurements of these statistics from a large enough set of mocks are not currently available. 
Therefore we  model the covariance matrix analytically, calibrating it, as much as possible, for the terms that are measurable, on estimates from a large number of mocks. We then use this analytical template, calibrated  and extended to terms involving the bispectrum higher multipoles, to evaluate the likelihood in the MCMC sampling of the parameters posterior distribution.

An  analytical model of the galaxy power spectrum and bispectrum multipoles covariance matrix  alternative to the one presented here, was recently derived in <cit.>, where the effect of including the anisotropic bispectrum signal was studied through the cumulative signal to noise ratio.


In order to model the galaxy bispectrum higher multipoles we use the choice of angles described in <cit.> which differs from the one recently proposed in <cit.>. The expansion is then performed in terms of the angle between the normal to the triangle's surface and the line of sight, together with the angle describing the rotation of the triangle on the plane defined by its sides (or in other words the rotation around the normal to the triangle's surface).

The structure of the paper is the following.
In Section <ref> we give an overview of the objectives and the necessary steps to achieve them.
In Section <ref> we describe the analytical expressions used to model the data vector while in Section <ref> we list the analytical terms that we included in our covariance matrix model.
Section <ref> is of critical importance since it describes the qualitative and quantitative tests that we performed in order to check the accuracy and reliability of our covariance analytical model. In Section <ref> we forecast through MCMC sampling the parameter constraints improvements obtained by adding to the data vector the galaxy bispectrum higher multipoles. We conclude in Section <ref>.












§ STATEMENT OF THE PROBLEM, OBJECTIVES AND ROADMAP


It is not at all a trivial task to develop the estimators and the associated numerical pipelines to measure the galaxy-redshift space bispectrum anisotropic signal in terms of multipoles expansion  <cit.>.
Hence it is important to ask: is it worth it?
There are two aspects to this question. One is the simple potential in reducing parameters errors or lifting parameters degeneracies. The other aspect is that in any practical application,  these extra statistics can offer a powerful consistency check of the modelling and thus  possibly of the model itself. The second aspect, while possibly of more consequence than the first, cannot be easily demonstrated and quantified before one is  able to measure these statistics on mock survey catalogs (mocks for short). However,   it is reasonable to expect that these statistics can provide a useful consistency checks if they carry significant additional information, and therefore only if they are useful to reduce parameters errors.    
This is what we set up to address in  this work.



Our goal is indeed to show whether adding the bispectrum anisotropic signal to the standard power spectrum P_g^(ℓ) (isotropic + anisotropic) and bispectrum B_g^(00) (isotropic) data vector, produces tighter 1-2D marginalised posterior distributions for the parameters usually constrained by galaxy-clustering analyses. We prove this through the path described  in the conceptual map of Figure  <ref>.

Given the absence of measurements, the first task consists in analytically modelling the bispectrum anistropic signal.  We choose a decomposition in terms of Legendre polynomials for a specific choice of angles describing the orientation with respect to the line of sight. We  then proceed to use this theoretical model as a "mock"[The meaning of "mock" here is different from the more widespread use of the term. Usually "mock" data refer to  carefully simulated data in such a way as to have properties as similar as possible to the real data. Here our generated data vector is somewhat idealized, being obtained from theoretical templates. For this reason we use the term "mock" here in quotes. In the rest of the paper we will  use the word mock in its more common meaning and refer for brevity to the template-generated  data as (fiducial) template data vector. In some applications we will use a data vector obtained from the survey mocks. We will refer to this as simulated data vector.] data vector. A similar approach was used in <cit.> to derive forecasts on the sum of the neutrino masses constraints.


In summary our "data" vector is composed by four blocks, galaxy power spectrum (isotropic plus anisotropic) and galaxy bispectrum (also isotropic plus anisotropic):
 [P^(0)_g;P^(2)_g;B^(00)_g;B^()_g], where B^()_g indicates all the possible bispectrum higher multipoles. Each power spectrum block includes all the k-bins (bandpowers) while each bispectrum block contains all the triangle configurations used.
 

As it is widely done in large-scale structure studies, we assume the likelihood of our data vector to be Gaussian, with a covariance matrix that does not depend on the cosmological parameters. This assumption has been extensively discussed, tested, and motivated, see e.g., <cit.> for an insightful discussion. 

To perform likelihood-based parameter inference, 
 we also need the covariance matrix.  The covariance matrix  is usually estimated numerically from the data vector's measurements on a large enough set of simulations or mock survey data. 
In our case we only have measurements for a reduced data vector which includes the standard galaxy clustering statistics, [P_g^(ℓ);B_g^(00)] and therefore only for a sub-matrix.  To overcome this obstacle, we derive an analytical model for the full covariance matrix to be used as a template. This template can then be calibrated using the numerically estimated covariance for the sub matrix involving [P_g^(ℓ);B_g^(00)], for which we have measurements from 1400 galaxy mock catalogues. The calibration is done by fitting parameters which are common to both the anisotropic and isotropic components of the  full data vector's covariance and therefore independent of the multipole expansion.
Two of them are amplitudes parameters (one for each of the power spectrum and bispectrum parts respectively) and the third one is an effective average galaxy number density (modelling the shot-noise). The fit is done by comparing our analytical template to the covariance matrix diagonal estimated from the mocks. For the calibration we only use the diagonal elements since  they are signal-dominated (i.e.,  their statistical noise is negligible with respect to the physical signal).
A similar approach was followed by Ref.  <cit.> to compute the three-point correlation function's covariance matrix .
In the analytic expression for the covariance matrix it is hard to include the  extra mode-coupling and  coupling between different statistics, induced by the window. This will be taken into consideration when performing the  comparison to the numerical estimate of the covariance.

We use the analytical model for the  power spectrum multipoles covariance only for the calibration process.  In our forecasts,  in order to minimize potential systematic errors connected with the power spectrum covariance modelling, for this part of the data vector, [P_g^(ℓ)], we will use the covariance numerically estimated from the simulations. The resulting "hybrid" covariance model will be referenced as HYB in the figures.

The following step is to verify that our covariance analytical model, when used to sample via a Markov chain Monte Carlo (MCMC)  the parameters posterior distribution, produces the same results obtained with the  numerical covariance matrix. 

We test this for the case of the anisotropic power spectrum plus isotropic bispectrum data vector (for which we have measurements from a large set of galaxy mock catalogues).
In particular we compare the posterior distributions 1D medians and 68% credible intervals  for a set of 1400 galaxy mock catalogues.

Provided that there is no significant difference between the parameters posteriors obtained using  the analytical covariance matrix template and the numerically estimated one, we  can  assume that the calibrated covariance model can  safely be extended  to the bispectrum multipoles.

The analytical model of the extended covariance matrix can  then be used to evaluate the effect  on cosmological parameter inference of including  higher  bispectrum multipoles in  the data vector.

We will show in Section <ref> two main  results. First we prove that the anisotropic bispectrum signal contains additional significant cosmological information. Second,  this improvement saturates at the level of the hexadecapole terms, meaning that measuring higher multipoles  is not expected to  further  tighten parameter constraints (i.e., it would not be worth the effort). 





§ SIGNAL MODELLING


Given the lack of measurements from simulations or "real world" surveys, in our work we simulate the bispectrum anisotropic signal using an analytical template. Therefore we employ theoretical expressions derived from standard perturbation theory (see <cit.> for a complete review) to generate our "fiducial" (but realistic) template data vector.

 The tree-level expression for the power spectrum and bispectrum are:

    P_g(𝐤)   =Z^s_1[𝐤]^2P_m( k )D_FoG^P(k)
    
    B_g(𝐤_1,𝐤_2,𝐤_3)    = [2 P_m(k_1)P_m(k_2)Z^s_1[𝐤_1]Z^s_1[𝐤_2]Z^s_2[𝐤_1,𝐤_2] .
    
       + . two cyclic  terms]×  D_FoG^B(k_1,k_2,k_3),


where the subscript "g" indicates that the statistics are computed for the galaxy field in redshift space. The standard perturbation theory kernels in redshift space used in this work (Z^s_n) are those reported in Appendix C of <cit.>.
The D_FoG terms model the Fingers-of-God effect (hereafter FoG) <cit.> introducing two effective parameters σ_P and σ_B for the power spectrum and bispectrum respectively. These parameters enter into two Lorentzian damping functions in front of the respective data vectors <cit.>:


    D_FoG^P   =11 + k^2μ^2σ_P^2/2  
    
        D_FoG^B   = 11 + (k_1^2μ_1^2+k_2^2μ_2^2+k_3^2μ_3^2)σ_B^2/2   .


Beside depending on the k-vectors coordinates, the Z_n^s kernels are also a function of the bias parameters b_1, b_2 and b_s^2 together with the growth rate f. 

As is it  widespread, we also make the approximation of the bias being local in Lagrangian space, which allows us to express the non-local bias term b_s^2 as function of the linear bias parameter b_1, b_s^2=-47(b_1-1) <cit.>.

In our modeling and analysis we include the Alcock-Paczyǹski <cit.> parameters to account for differences between the true underlying cosmology and the fiducial one used to convert redshift measurements into Euclidean distances. A discrepancy between the fiducial and true underlying cosmology induces distortions in the clustering properties. This can be modelled through the parallel and perpendicular dilation parameters  α_∥ = k_∥/p_∥ and α_⊥ = k_⊥/p_⊥ which relate the measured wave-vector 𝐤 with the corresponding wave-vector 𝐩 for a particular cosmology. Being μ and η the cosines of the angles between the line of sight and 𝐤 and 𝐩 respectively, the conversion is given by <cit.>


    p    = |𝐩| = kα_⊥[1+μ^2(F^-2-1)]^1/2
    η   = μF[1+μ^2(F^-2-1)]^-1/2 ,


where F≡α_∥ / α_⊥.










 §.§ Multipole expansion

We compute the multipole expansion of the tree-level galaxy power spectrum  with respect to the line of sight in the standard way


    P_g^(ℓ)(k) = (2ℓ + 1)2α_∥α_⊥^2∫^+1_-1dμ ℒ_ℓ(μ) P_g(p,η)


where μ is the cosine of the angle between the line of sight and the k-vector while ℒ_ℓ(μ) is the ℓ-th order Legendre polynomial.

For the bispectrum multipole expansion we adopt the choice of angles presented in <cit.>. The first angle is the one between the normal to the surface of the triangle and the line of sight, 𝐳̂, whose cosine we indicate by the letter μ. The second angle describes the rotation of the triangle around the normal to its surface, with a fantasy effort we call ν its cosine;


    μ = cosω   =   (𝐤_1×𝐤_2)·𝐳̂sinθ_12
    ν = cosϕ   =   [𝐳̂×(𝐤_1×𝐤_2)]·𝐤̂_1sinω ,


where cosθ_12=(𝐤_1·𝐤_2)/(k_1k_2). 
The origin of both angles can be for example imagined translating the line of sight and normal to the surface vectors to 𝐤_1's starting point (see Figure 1 of <cit.> for a visual representation).

For numerical reasons we still integrate over μ_1 and μ_2 (the cosine of the angles between the line of sight and the k-vectors k_1 and k_2) since the μ and ν are derived quantities from the standard bispectrum coordinates (k_1,k_2,k_3,μ_1,μ_2), conversion formula given in  Appendix <ref>, Equation <ref>. The bispectrum multipoles are then given by


    B_g^(,)(k_1,k_2,k_3) = (2+1)(2+1)4πα_∥^2α_⊥^4∫^+1_-1dμ_1∫^2π_0 dϕ ℒ_(μ)ℒ_(ν)B_g(p_1,p_2,p_3,η_1,η_2) ,


where ℒ_(μ) and ℒ_(ν) are the -th and -th order Legendre polynomials for μ and ν, respectively. The azimuthal angle ϕ around 𝐤_1 is defined in Section 3 Equation 18 of <cit.> as μ_2≡μ_1cosθ_12 - √(1-μ_1^2)√(1-cosθ_12^2)cosϕ, which was used also in <cit.>.  p_1, p_2, p_3 are the k-modules while η_1, η_2 are the cosines of the angles between 𝐤_1, 𝐤_2 and the line of sight corrected for the Alcock-Paczyǹski effect (Equation <ref>).




§ COVARIANCE MATRICES: FORMULAE


We update and extend the analytical results presented in <cit.>. The full derivation  and full expressions of the terms used in this section  are shown in Appendix <ref>.
 The covariance matrix can be divided into different blocks which, given the definition of the data vector components in Equations <ref> and <ref> we organize as follows:


    [       C^P^0_g  C^P^0_gP^2_g C^P^0_gB^00_g   C^P^0_gB^_g;  C^P^2_gP^0_g       C^P^2_g C^P^2_gB^00_g   C^P^2_gB^_g; C^P^0_gB^00_g C^P^2_gB^00_g      C^B^00_g  C^B^00_gB^_g;   C^P^0_gB^_g   C^P^2_gB^_g  C^B^_gB^00_g        C^B^_g;               ]


In the above covariance matrix sketch on the diagonal there are  blocks for both the power spectrum C^P^0_g and bispectrum multipoles C^B^_g auto-covariances. Off-diagonal blocks correspond to  the cross-covariances, between different multipoles of the same statistic, C^P^0_gP^2_g or C^B^_gB^00_g, and between different statistics, like C^P^0_gB^_g.
 The gray background indicates the covariance part that can be evaluated both analytically and numerically from the mock catalogues and therefore it can be used to calibrate the analytic expressions.

In all the auto-covariances and cross-covariances between different multipoles of the same statistic there are both Gaussian and non-Gaussian terms. The Gaussian terms are the ones obtained at lower order in perturbation theory, when the field is considered to be completely Gaussian and therefore all moments higher than the variance vanish. The Gaussian terms are non-zero only on the diagonal of each block. Non-Gaussian terms appear when higher-orders of the matter-density field are considered not only in modelling the data vector (where they are dominant for example in the case of the bispectrum), but also in the covariance matrix theoretical modelling. Without second-order terms (inducing a non-zero bispectrum) for example, there would not be any correlation between power spectrum and bispectrum data vector elements.

Differently from Ref. <cit.>, we limit our perturbation theory expansion of the non-Gaussian terms to the lowest order for each component of the covariance matrix ( ∝δ_m^6 in the cross-covariance and ∝δ_m^8 in the bispectrum auto-covariance, δ_m being the matter perturbation density variable). This implies not modelling higher-order terms such as the one proportional to the tetraspectrum (∝δ_m^8) in the cross-covariance between power spectrum and bispectrum together with the one proportional to the pentaspectrum (∝δ_m^10) in the bispectrum auto covariance.  According to  Figures 2 and 3 of Ref. <cit.> for our range of scales (k<0.09h Mpc^-1) and considering the scatter in the covariance matrix terms we will later discuss in Figure <ref>, the lack of the above mentioned higher-order terms is expected to  have negligible impact in our analysis. Moreover, since we calibrate the amplitude of the covariance  terms on a suite of mocks (see Sec. <ref>), any possible residual is likely absorbed by this normalization.
Full expressions and more extended discussion can be found in Appendix <ref>.


For the covariance between different galaxy power spectrum multipoles we  consider only the Gaussian term since on linear scales  it is dominant compared  to the non-Gaussian part  <cit.>. In particular, the covariance between the power spectrum ℓ_1-multipole for mode k_1 and the power spectrum ℓ_2-multipole for mode k_2 is given by 



    C_G^P^ℓ_1_gP^ℓ_2_g(k_1;k_2) = 
    (2ℓ_1+1)(2ℓ_2+1)
    2(2π)^3δ^K_124π k_1^2Δ kV_s  12∫^1_-1dμ_1 ℒ_ℓ_1(μ_1)ℒ_ℓ_2(μ_1)
    P_g(k_1,μ_1)^2,


where ℓ_1,ℓ_2=0,2. V_s is the survey effective volume (where for the analytical model we approximate V_s=(2π)^3k_f^-3), Δ k is chosen resolution in Fourier space for the power spectrum data vector. 

In order to correct for the shot-noise, in all our modelling we implicitly adopt the expressions derived in <cit.>, for example P_g(k_1,μ_1)≡P^N_g(k_1,μ_1)=P_g(k_1,μ_1) + 1/n_g. This  is appropriate because the estimators used to measure the power spectrum and bispectrum statistics from galaxy catalogues and data have the Poisson shot-noise subtraction implemented. However the shot-noise contribution must still be included in the covariance matrix.   


The Kronecker's delta  δ^K_12 is selecting only equal k's (in the sense of |k_a-k_b|<Δ k/2). The above basic model for the power spectrum multipoles covariance matrix diagonal only serves to the purpose of fitting the amplitude parameter needed for the cross-correlation term between power spectrum and bispectrum as described later in Section <ref>. Indeed in our analytical covariance matrix, for the power spectrum multipoles auto-covariance, we use the one estimated from the mocks, hence the name "hybrid" covariance matrix (referenced as HYB in the figures). In this way we avoid any possible source of systematic error linked with a non-optimal modelling of the power spectrum covariance. For a recent specific work on the accurate modelling of the power spectrum multipoles covariance matrix see <cit.>.

In order to deal with the many indexes needed for the galaxy bispectrum multipole expansion we introduce the following notation.
Let's say we want to  compute the covariance between two  different bispectrum configurations a, b, and  multipoles, B_g^(αβ) and B_g^(κλ). The indices α, β etc. correspond to the multipoles ℓ=0,2,4, where the first index refers to the angle cosine μ and the second to ν. 
Therefore the Legendre polynomials for both angle cosines μ and ν will be written as ℒ^(α)_μ_a, ℒ^(β)_ν_a for the  triangle a, and ℒ^(κ)_μ_b, ℒ^(λ)_ν_b for the triangle b. For the sake of brevity we introduce the short notation ℳ^(αβκλ)_μν,ab =ℒ^(α)_μ_aℒ^(β)_ν_aℒ^(κ)_μ_bℒ^(λ)_ν_b with the added benefit of appearing cool general relativity-like. In the same way we write the compact pre-factor 𝒞^αβ_κλ = (2α+1)(2β+1)(2κ+1)(2λ+1). 


The bispectrum covariance matrix at lowest order in terms of matter overdensity field perturbative expansion is composed of three terms: the first is proportional to the product of three power spectra (Gaussian), the second and the third are proportional to the product of two bispectra and power spectrum times trispectrum (non-Gaussian), respectively. The six-points correlator generating these contributions is reported in Equation <ref>. 
The galaxy bispectrum multipoles covariance Gaussian term (i.e., first line of Eq. <ref>) is then given by


    C_∝P^3^B^αβ_gB^κλ_g(k_1,k_2,k_3;k_4,k_5,k_6)    =
    D^123_456 (2π)^6𝒞^αβ_κλV^t_123V_s
       ×14π∫^+1_-1dμ_1∫^2π_0dϕℳ^(αβκλ)_μν,aaP_g(k_1,μ_1)P_g(k_2,μ_2)P_g(k_3,μ_3),


where D^123_456 stands for all the possible indexes permutations of the product δ^K_14δ^K_25δ^K_36 and has values 6,2,1 respectively for equilateral, isosceles and scalene triangles. Finally V^t_123=8π^2k_1k_2k_3Δ k^3 is the integration volume in Fourier space for the bispectrum. Notice that for the Gaussian term to be non-zero the triangles a and b are identical; that is why we wrote the subscript aa for  ℳ^(αβκλ)_μν,aa.

Up to order ∝δ^8 in terms of the matter over-density perturbation,  there are only two non-Gaussian terms contributing to the bispectrum covariance (see eq. <ref>) since the third one, proportional to the pentaspectrum, is of order ∝δ^10 <cit.>. The first is proportional to the product of two bispectra:


    C_∝BB^B^αβ_gB^κλ_g(k_1,k_2,k_3;k_4,k_5,k_6) = 
    (2π)^3𝒞^αβ_κλδ^K_344π k_3k_4Δ k V_s
       ×18∫_-1^+1dμ_3dμ_1dμ_5 
    ℳ^(αβκλ)_μν,abB_g^s(k_1,k_2,k_3,μ_3,μ_1)B_g^s(k_3,k_5,k_6,μ_3,μ_5)   + 8  perm.
    ,


for which in total there are nine possible permutations. The second term is proportional to the product of the power spectrum and trispectrum:



    C_∝TP^B^αβ_gB^κλ_g(k_1,k_2,k_3;k_4,k_5,k_6) = 
    δ^K_34(2π)^3𝒞^αβ_κλΔ D4π k_3k_4Δ k_3Δ k_4V_s
       ×18∫^+1_-1dμ_Ddμ_1dμ_5
    T_g(k_1,k_2,k_5,k_6,D,μ_D,μ_1,μ_5)P_g(D,-μ_D)
      + 8  perm. ,


where k_3 = k_4 = D is the diagonal of the quadrilateral configuration defining the trispectrum. Before the angular integration, eight coordinates are needed in order to describe a galaxy trispectrum configuration. Differently from what we had for the previous non-Gaussian term in Equation (<ref>), δ_D(𝐪_3-𝐪_4), in this case the Kronecker's delta condition derives from δ_D(𝐪_3+𝐪_4) and that is the reason for having -μ_D as second argument of the galaxy power spectrum in the integrand (even if it actually does not make any difference at tree-level). More detailed explanation about the trispectrum are reported in Appendix <ref>.

Finally we model at lowest order the cross-correlation between the galaxy power spectrum and bispectrum multipoles by computing (for the five-point correlator originating this term see Equation <ref>):



    C_∝PB^P^ℓ_gB^αβ_g(k_1;k_2,k_3,k_4)    =
    2 (2π)^3𝒞_αβ^ℓδ^K_124π k_2^2 Δ kV_s
       ×14π∫^+1_-1∫^2π_0dϕℳ^(ℓαβ)_μ_aμ_bν_bP_g^s(k_2,μ_2)
    B_g^s(k_2,k_3,k_4,μ_2,ϕ)   + 2  perm. .


 The Δ k appearing at the denominator corresponds to the  Fourier space  resolution adopted for the power spectrum  (i.e., giving the number of  power spectrum modes available in Fourier space). The bispectrum has its own Δ k which does not have to coincide with the one chosen for the power spectrum, but does not appear here[As shown in Equation  <ref> the bispectrum  Δ k cancels out.]. For the power spectrum, the multipole expansion is done in terms of μ_a, the cosine of the angle between the k-vector and the line of sight. The angle ϕ is the same as the one defined below Equation (<ref>).










§ NUMERICAL AND HYBRID COVARIANCE MATRIX FOR 



We compare the analytical covariance matrix for the data vector [P_g^(0,2);B_g^(00)] with a numerical one derived by measuring the above statistics from 1400 realisations of the MultiDark Patchy BOSS DR12 galaxy catalogues <cit.>. The Patchy mocks reproduce the number density, selection function, survey geometry and clustering properties of BOSS DR12 galaxy sample <cit.>; for each mock the power spectrum monopole and quadrupole are available and we have calculated the bispectrum monopole for both North and South Galactic Caps (NGC and SGC).

The same fiducial cosmology of the mocks, Ω_Λ = 0.693, Ω_m(z=0) = 0.307, Ω_b(z=0) = 0.048, σ_8(z=0)=0.829, n_s = 0.96, h_0 = 0.678, has been used the generate the necessary analytical matter power spectrum <cit.>  for our modeling.

Note that the covariance matrix obtained from the  survey mocks has a contribution due to the survey  window function. The analytic expression for the covariance matrix does not account for the  extra mode-coupling and thus coupling between different statistics, induced by the window. 
In  our analytically-generated data vector used for  parameter  forecast  we do not include the survey window effect, hence we do not worry too much about modeling the same effect in the covariance matrix.  More discussion on the effect of the window can be found below and in Appendix <ref>.
In order to model the window additional contribution to statistical error one could proceed as in <cit.> for the CMB angular power spectrum: under the assumption of the mask correction being small with respect to the signal covariance, the total covariance could be written as a sum C_tot = C_signal + C_mask.

Similarly to  <cit.> we set the bias parameters values b_1 = 2.0 and b_2=0., together with using the same effective survey volume V_s = 1.76 [h^-1 Gpc]^3. The choice for the linear bias b_1 is not too dissimilar from what obtained in the BOSS clustering analysis <cit.> while the survey volume derived in <cit.> is not too far from the one (V_s = 1.92 [h^-1 Gpc]^3) fitted from comparing the numerical covariance estimated from the mocks and its theoretical model in the 3-pt correlation function analysis done in <cit.>. 





In the power spectrum multipoles covariance matrix expression we set Δ k = 0.01 h/Mpc, which is the same value used in numerical measurement pipeline. For the galaxy bispectrum we use two different Δ k bin sizes proportional to the fundamental frequency k_f = (2π)^3/V_b where V_b = (3500 Mpc/h)^3 is the survey volume for the mock catalogues cubic box. We select the two cases: Δ k_6,5 = 6,5 × k_f which given our algorithm <cit.> to generate (k_1,k_2,k_3) triplets imposing the triangle closure condition and the range of scales considered,  corresponds to using 32 and 68 triangle configurations, respectively. In the main body of this work all the results for the Δ k_5 case have been presented. The key figures and table for the Δ k_6 case are shown in Appendix <ref>.

We adopt a very conservative k-range for the sides of the triangle configurations used for the bispectrum multipoles, in particular we set k_min = 0.03  h/Mpc and k_max = 0.09  h/Mpc. The reason for choosing a low k_max is twofold. Firstly, this allows us to trust our standard perturbation theory calculations at tree-level, for example for what concerns the cross-correlation between power spectrum and bispectrum as noted in <cit.>. Secondly, this choice limits the number of triangle configurations  in the data vector. This is essential in order to make it feasible to analytically compute in a reasonable amount of time all the terms of the covariance matrix with a sufficient precision. Moreover, the longer the data vector, the less numerically stable is the inversion of its covariance matrix needed for the likelihood evaluation.  As already pointed  out in  <cit.> this is the case when several multipoles of the bispectrum are present simultaneously in the data vector without compression.

Nevertheless by choosing k_max = 0.09  h/Mpc our analysis uses only a fraction of the baryonic acoustic oscillations  information encoded in  the galaxy density field. To fully capture and exploit this information and in particular its anisotropic component using the bispectrum multipoles, it would be necessary to further improve the covariance analytical model to extend its reliability to higher k's. For this reason our finding might be conservative.

In order to choose the scalar amplitude parameter A_s needed to compute the linear matter power spectrum (with CLASS <cit.>) to be used in the analytical covariance matrix template, we run a preliminary MCMC sampling using the numerical covariance matrix and the mean of the mocks as data vector. We then leave  A_s free to vary together with the additional parameters (b_1,b_2,f,α_⊥, α_∥). We obtain A_s≃ 1.9 × 10^-9 as a best-fit value, we recognize that this is  lower than the value reported in the Planck analysis <cit.> (A_s=3.04×10^-9). This is just a shortcut to make our calculations numerically stable. In fact, as reported below, and explained in Fig. <ref> the effective amplitude of the power spectrum and bispectrum  covariances are later free parameters of the normalisation.

For the covariance computation we set the FoG parameters in Equation <ref> for the power spectrum and bispectrum multipoles to the best fit values obtained in the BOSS analysis <cit.>, σ_P=3.5 and σ_B=7.39. In fact, even if we limit our k-range to linear scales where the FoG effect should be negligible, in  some of the covariance matrix  terms (in particular, shot-noise correction terms)  the sums of two or three k-vectors  appear in the argument of the power spectrum or bispectrum.  This is not the case for the terms proportional to the trispectrum and hence no FoG parameter is  used for the trispectrum. In the data vector's model used for the MCMC samplings both σ_P and σ_B have been set to zero since we limit our analysis to linear scales. 

To calibrate the covariance matrix we adjust three parameters:
 the average number density of galaxies n_g, an effective amplitude for the power spectrum part of the covariance A_P and another effective amplitude A_B for the bispectrum part. These parameters act on the calibrated terms as shown in Figure <ref>. Figure <ref> clarifies the difference between fitting the two amplitude parameters A_P, A_B and the physical scalar amplitude A_s. While A_s is common to all the terms in the covariance, each is term proportional to a different power of the scalar perturbations amplitude parameter (see caption of Figure <ref> for details), A_P and A_B can be interpreted as separate corrections to the analytically-computed integration volumes in Fourier space for power spectrum (V_p = 4πΔ k k^2) and bispectrum (V^t_123 = 8π^2k_1k_2k_3Δ k_1Δ k_2Δ k_3). 
 The analytical expressions for these volumes, which are used to normalise the statistics estimators by the number of modes measured (see Appendix A of <cit.> for a detailed explanation), are obtained in the "infinitesimal" bin-width approximation. 
In particular using Δ_k=5,6× k_f and limiting our analysis to small k's the condition k≪Δ_k may not be optimal. Moreover, A_P, A_B can also account for a mismatch between the fiducial effective survey volume value and the real one. Therefore it is not surprising that a calibration of these quantities is needed when they are not measured but computed analytically; reassuringly, their value is of the order of unity as reported in Table <ref>.

 For what concerns n_g, we treat it as a free parameter in order to model shot-noise deviations from  Poissonian behaviour. This is non uncommon and often also used when analyzing galaxy surveys. Indeed the expressions that we use to account for the shot-noise terms in the covariance model were derived under the Poissonian approximation (see Ref. <cit.>). Even letting n_g free to vary, we  find a best fit value very similar to the one computed  in Ref. <cit.> for the same data-set.
 
 
 
 We fit these parameters by comparing the diagonal of the analytical covariance matrix with the one numerically estimated using a set of 1400 Patchy galaxy mock catalogue realisations after correction by the  Hartlap factor h_f = (N_mocks - N_dv -2)/(N_mocks-1) <cit.>, where N_dv is the data vector's size while N_mocks is the number of galaxy catalogues used to estimate the covariance (1400 in our case). 
In summary calibration is done following the expression below:


    C_tot^P^(ℓ),B^(,) = 
    
    A_P×[C_P^ℓ^n_g^0 + 1n_gC_P^ℓ^n_g^-1+1n_g^2C_P^ℓ^n_g^-2]
    ,    P^(ℓ)_g
     
    A_B×[C_B^^n_g^0 + 1n_gC_B^^n_g^-1
    +1n_g^2C_B^^n_g^-2
    +1n_g^3C_B^^n_g^-3
    +1n_g^4C_B^^n_g^-4]
    ,    B^()_g
    √(A_PA_B)×[C_P^(ℓ)_gB^^n_g^0 + 1n_gC_P^(ℓ)_gB^^n_g^-1
    +1n_g^2C_P^(ℓ)_gB^^n_g^-2
    +1n_g^3C_P^(ℓ)_gB^^n_g^-3]
    ,    P^(ℓ)_gB^()_g



In Table <ref> are reported the values for the fitted parameters both in the Δ k_5 and Δ k_6 cases.









 §.§ Visual comparison 

We begin by visually inspecting the differences between the analytical covariance matrix (after calibration) and the  numerical.  Figure <ref>  shows  1D sections (columns/rows) for both the bispectrum autocovariance and its cross-correlation with the power spectrum for the Δ k_5 case.
The shaded area is obtained by splitting the 1400 mocks into ten groups and then computing the numerical covariance matrix for each of them. From the ten estimated covariance matrices we compute each element's standard deviation as a proxy for the statistical noise. This gives us a possible estimate of the level of agreement between analytical and numerical covariance matrices.
The bispectrum auto-covariance results are well approximated (mostly within the shaded area) by the analytical model. 

The cross-covariance between bispectrum and power spectrum multipoles shows a larger discrepancy between numerical and analytical covariance matrices. In particular the amplitude of certain elements seems to be lower in the analytical case. For other "pixels" the analytical value is exactly equal to zero due to the Kronecker's delta condition from Equation <ref>. The absence in our model of higher order terms such as the one proportional to the tetraspectrum could explain the lower amplitude observed in the analytical lines in the right panel of Figure <ref>. In the pixels where the analytical covariance matrix is by definition equal to zero we see that the shaded area indicates that also the numerical estimate is compatible with zero and therefore the non-null values seen in the numerical covariance matrix could be interpreted simply as statistical noise. Another source of discrepancy could be the lack of the survey window selection effect in our covariance template (more below).
Finally the shot-noise correction appears to be relatively negligible with respect to the non-Gaussian terms for the off-diagonal elements.  


It also results useful to look at the whole 2D analytical and numerical reduced covariance matrices (i.e., element of row i and column j is divided by the square root of the product of element ii times element jj). In Figure <ref> are shown the reduced covariance matrices for the Δ k_5 binning case[Another reason to choose this binning case for the whole analysis is that there are enough triangle configurations to resolve the structure and features present in the covariance matrix but not too many to make them invisible due to the limited size of the image.]. The comparison between the numerical and analytical covariance reveals qualitative good agreement. The only relevant difference is the lack, in the calibrated hybrid template, of the features parallel to the covariance diagonal. 
These features are due to the additional mode-coupling induced by the survey window.  In Appendix <ref> we compare our hybrid model to the numerical covariance obtained from a set of 440 simulated galaxy fields in cubic boxes with periodic boundary conditions <cit.>. In Figure <ref> we see that in the absence of survey mask and window, such features are not present in  the covariance matrix. It will become evident in Section <ref> that these features have no significant effect on our analysis.











 §.§ Constraints comparison for the (reduced)  data vector 



We estimate the  cosmological posterior for the (reduced) data vector [ 𝐏_𝐠^(0,2);𝐁_𝐠^(00)] using the analytical/hybrid covariance matrix and the numerical one.
Whenever we compare the two covariance matrices we use for each test the same data vector. In this section we either use the mean measurement from the mocks or each individual mock measurement. 

To compare the two posteriors so obtained, we consider a seven parameters set:
   (b_1,b_2,f,σ_8,f_NL,α_⊥, α_∥), where
b_1 and b_2 are the first and second order Eulerian bias parameters, f is the linear growth-rate, σ_8 is the normalisation of the matter power spectrum, f_NL is the amplitude of local primordial non-Gaussianity perturbations <cit.> and α_⊥, α_∥ are the Alcock-Paczyǹsky parameters. 

Following what is done in the BOSS analyses <cit.> we use the local Lagrangian bias approximation to fix the value of the non-local tidal bias parameter b_s_2 appearing in the second order kernel of the bispectrum to b_s_2=-4(b_1-1)/7. We have checked through a specific MCMC sampling that allowing this parameter to vary does not significantly change our results.

This is shown in Figure  <ref> for the Δ k_5 case (for which there are 68 triangle configurations in the bispectrum data vector). 

Since we focus on the change in the parameter constraints width, we show our results with each parameter's central value (median) subtracted. The central value is given by the analysis performed using the numerical covariance matrix. A full parameters constraints analysis will be presented elsewhere.


From Figure <ref> we can see that qualitatively there is very good agreement between the posterior distributions obtained with the hybrid and numerical covariance matrices. 
With another separate MCMC run we tested what happens in both cases:  setting to zero all the off-diagonal terms or only the bispectrum auto-covariance ones. The resulting marginalised 1-2D posterior distributions are- not unexpetedly- significantly tighter then the ones obtained using the covariance matrix estimated from the mocks.



We go even further to quantify the performance of the hybrid covariance matrix approach: we perform a MCMC sampling on the data vector measurement from each of 1400 mock catalogues using both analytic and the hybrid covariance matrices.  


We then compare in Figure <ref> the probability distributions 1D 68% confidence intervals and the medians for each model parameter. The relative difference between variances and medians obtained using numerical and hybrid covariances is for the large majority of the mocks less than ∼ 25% the value of each parameter's variance.
In Figure <ref> we show the 1 and 2D posterior distributions for the parameters and the  scatter of the median[As it is common to do in Bayesian inference, we stick to the so-called Central Credible Interval which reports  errors around the median of a distribution. Posteriors are very symmetric hence the median is very close to the mean and the maximum anyway.] values for each mock realisation.  Point and lines in red correspond to the analytic/hybrid covariance, points and lines in gray correspond to the numerical covariance. 


From these quantitative and statistical tests we conclude that our analytical model, once  calibrated (i.e., fitting the three parameters A_P, A_B, n_g), represents a realistic proxy for the covariance matrix estimated from the galaxy mock catalogues for the [P_g^(0,2);B_g^(00)] data vector.

Given this result, the main assumption of this work is to consider reliable the extension of our covariance model to the galaxy bispectrum higher multipoles. This is motivated by the fact that  our choice of fitting parameters does not depend on the multipole chosen but just on the statistic (power spectrum or bispectrum).





§ HIGHER MULTIPOLES OF THE BISPECTRUM:




We now proceed to estimate the added value given by extending the clustering statistics data vector to the galaxy bispectrum higher multipoles.  We use the covariance analytical model presented in Section <ref> to sample the model parameters posterior distribution. In Figure <ref> the resulting reduced covariance matrix for the data vector [P_g^(0,2);B^(00,20,02,40,04,22)_g] is displayed. As an example it is interesting to notice that bispectrum monopole B_g^(00) is negatively correlated with B_g^(20) while positively correlated with B_g^(02). 
The difference in the correlations between different bispectrum multipoles terms suggests that the multipoles have complementary constraining power.












 §.§ Cumulative signal to noise ratio

A first preliminary study can be done by looking at the comulative signal to noise ratio (S/N) as a function of the range of scale used (varying k_max in our case). Given a data vector 𝐱 and its covariance matrix C_𝐱, the squared cumulative signal to noise ratio is given by


    (S/N)^2   =   𝐱^⊺ C_𝐱 𝐱 .


By removing from the bispectrum data vector those triangles having at least one side larger that the chosen k_max we show in Figure <ref> the cumulative signal to noise ratio as a function of the considered scales range.

 In Figure <ref> we look into the signal to noise ratio for different combinations of the data vector's terms (by rows) and different versions of the covariance matrices (by columns). The first message that this plot conveys is that the comulative signal to noise ratio seems to saturate when all three bispectrum hexadecapoles terms are employed in the analysis. Moreover, moving from left to right from the first to second column, it results evident that the overall decrease in signal to noise is caused by the addition of the off-diagonal terms to the covariance matrix. In comparison, adding the cross-correlation between power spectrum and bispectrum has a smaller impact (from second to third column).

 On the third row on the right, when the whole covariance is used (off-diagonal terms and cross-covariances between power spectrum and bispectrum multipoles) there is a  small discrepancy between the (S/N) ratios computed using the numerical and the analytical covariances. The reason for this could be the fact that, as we previously saw in Figure <ref>, in both the analytically computed bispectrum covariance and cross-covariance terms certain elements are by definition (Kronecker's deltas relations) set to zero, while, because of noise, they are not zero in the numerical covariance. 

Comparing this to the cumulative signal to noise ratio values displayed in Figure 7 of <cit.> we find larger values for the same k-range because we used an analytical model as data vector 𝐱 while in <cit.> the authors used the measurement from the mean of the mocks which contains noise. An example of this can be seen by comparing the contours for the same data vector   [P_g^(0,2);B_g^(00)] (gray and red) in Figures <ref> and <ref> where in the first case the mean of the mocks has been used as data vector measurement while in the second it has been used the analytical model. The 1- and 2-D marginalised contours are much wider in the first case, exactly because of the noise present in a measurement from simulations/data.











 §.§ Improvement on  parameter constraints



Because of the lack of measurements for the multipoles B_g^(20), B_g^(02), B_g^(40) , B_g^(04) and B_g^(22)  we use an analytical  data vector  computed for our fiducial cosmology. 
Figure <ref> shows the 1 and 2D posterior distributions of the cosmological parameters for data vectors of different lengths, starting from the basic power spectrum monopole and quadruple and bispectrum monopole, higher bispectrum multiplies are gradually added.  Since in the previous tests  comparing analytic/hybrid and numerical  covariance matrix the data vector was taken from the mocks, in Figure <ref> (gray and red) we  repeat the comparison of sec. <ref> using analytic the data vector. 
Also in the case of the analytical template for the data vector, the results obtained with numerical and hybrid covariance matrices are in very good agreement. 

The 1-2D marginalised contours shown in Figure <ref> confirm what deduced from the signal to noise study (Figure <ref>): saturation is achieved  at the bispectrum hexadecapole level.

We observe a  shrinkage of the 1- and 2-D marginalised posterior distributions as galaxy bispectrum higher multipoles are employed with respect to the monopole term alone. This aspect, together with the saturation, is quantitatively described in Table <ref>.
This analysis suggests that including all the terms up to the galaxy bispectrum hexadecapoles (together with the power spectrum monopole and quadrupole), results in  improving the 1D 68% credible regions by a factor of ∼ 30% on average – equivalent to observing a volume in the sky ∼ 2 times larger (V_()∝ 1/σ^2_≃ 1/(0.7 σ_(00))^2 = 2 × V_(00)).


Recall that our analysis is very conservative in terms of scales, by setting k_max = 0.09 h/Mpc we restrict our investigation to the linear regime. Extension to the mildly non-linear regime may offer additional gains. The  cumulative signal to noise shown in Figure <ref> seems to hint that, as k_max increases, the difference between the isotropic (monopole) and isotropic plus anistropic (monopole plus higher multipoles) bispectrum signal may also increase. 


This trend will eventually stop, as Figure 7 of <cit.> seems to indicate around k_max≃ 0.20   h/Mpc, however it would leave a non negligible gap between the isotropic and anisotropic signal to noise ratios. Therefore as smaller and smaller scales are considered, the signal to noise difference would likely translate into a stronger impact of the anisotropic bispectrum signal in terms of parameter constraints improvements. Full investigation in this direction  is left for future work.













 §.§ Local Primordial non-Gaussianity


The bispectrum has been widely studied in the literature <cit.> as a powerful tool to constrain potential deviations from Gaussianity of the primordial gravitational potential <cit.>.
By reducing the degeneracies between cosmological parameters, the galaxy bispectrum higher multipoles can play an important role in obtaining late-times constraints on the primordial non-Gaussianity amplitude parameter f_NL <cit.> competitive with the early times ones produced by cosmic microwave background (CMB) experiment such as Planck <cit.>.
Moreover constraining f_NL using the bispectrum multipoles does not rely on extremely large scales and on highly biased tracers as the current state-of-the-art analysis of the power spectrum scale dependent bias <cit.>.

In our estimate of the signature of primordial non-Gaussianity in the anisotropic bispectrum we do not include the modelling of the terms associated with the scale dependent bias induced by a local f_NL <cit.> and arising from the multivariate bias expansion <cit.> applied to the bispectrum. This effect is particularly strong at very large scales and an analysis aimed at deriving constraints on the local non-Gaussianity amplitude parameter f_NL using the bispectrum should implement the modelling developed in <cit.> which also includes the multivariate biasing expansion in redshift-space.
Our analysis can be interpreted as a conservative estimate using only the primordial signal appearing from a non vanishing initial bispectrum, underestimating the actual signature given by f_NL≠0. This would be larger and likely less correlated with the other cosmological parameters, especially on large scales, if we would include the non-local bias effect.
This is possibly one of the reasons why we obtain larger confidence intervals for f_NL than the forecasts made through a Fisher information analysis in <cit.>.


In particular redshift space distortions add an additional difference in the standard perturbation theory kernels dependence between the gravitational collapse (B^()_g,112∝Z^s_1Z^s_1Z^s_2)  and primordial (B^()_g,111∝Z^s_1Z^s_1Z^s_1) components of the bispectrum which becomes more explicit when the anisotropic signal is studied.

In Figure <ref> we show the local primordial non-Gaussianity contributions to the bispectrum <cit.> for the different multipoles, setting f_NL=5. These terms have been computed using the expression reported in the Appendix of <cit.>. The first row shows, for the different multipoles, the ratio between the primordial signal[It is interesting to notice that the multipoles that appear to be more sensitive to (richer in signal for) local primordial non-Gaussianity, are the ones relative to the angle describing the triangle rotation around the normal to the surface. An hint for a potential explanation to this effect could be searched in the angle conversion Equation <ref>. The cosine ν of the angle describing the rotation of the triangle around the normal to the surface depends explicitly on cosα_ab. This is the cosine of the internal angle of the triangle projection on the plane perpendicular to the line of sight. Therefore ν explicitly depends on how much squeezed is the triangle, which is a factor known to influence the sensitivity to primordial non-Gaussianity <cit.>.] 
B^()_g,111 and the gravitational collapse part of the bispectrum B^()_g,112.
In the second panel from the top, the ratio is between B^()_g,111 and the variance computed using the analytical covariance matrix model used in the paper's analysis. On the third row, the comulative signal-to-noise ratio is shown as a function of the number of triangle configurations included in the data vector. Finally in the bottom panel it is possible to see the  total signal (primordial plus gravitational) for the different multipoles as a function of the triangle configurations considered.




§ CONCLUSIONS



At the power spectrum level, it is well known that  the anisotropic modulation of the signal with respect to the line-of-sight induced by redshift space distortions  carries valuable cosmological information. This anisotropy is usually expressed in multipole moments: the monopole  –being the averaged signal in all directions around the line of sight– and the quadruple. 
The bispectrum has been proved useful for breaking degeneracies among parameters and detecting specific physical signatures (e.g., the so-called GR effects or neutrino masses). To date, only the bispectrum monopole has been used to constrain cosmological parameters from galaxy spectroscopic surveys. Measuring higher-order multipoles of the redshift-space galaxy bispectrum is indeed challenging and computationally expensive <cit.>. In this paper we have addressed the question: given these challenges, is it worth it?
We consider that  the answer is "yes", if adding the bispectrum anisotropic signal to the standard power spectrum P_g^(ℓ) (isotropic + anisotropic) and bispectrum B_g^(00) (isotropic) data vector, produces tighter 1- and 2-D marginalised posterior distributions for the parameters usually constrained by galaxy-clustering analyses.

To do so, we developed  an analytical model for the  full covariance matrix including power spectrum monopole and quadruple and bispectrum monopole and higher-orders multipoles up to the hexadecapole, for a conservative range of scales (k_max=0.09 h/Mpc). We concentrated on the set of cosmological parameters  (b_1,b_2,f,σ_8,f_NL,α_⊥, α_∥).

We first calibrate the analytical covariance matrix on the one estimated  numerically from a suite of 1400 mock survey catalogs. The calibration is done on the sub-matrix corresponding to the reduced data vector [P_g^(0,2);B_g^(00)]. We have shown that  the parameter's 1- and 2-D marginalised  posteriors distributions obtained with our analytical covariance and the numerical one are very similar and basically indistinguishable for  practical applications of reporting error-bars on parameters values (Figures  <ref>, <ref>, <ref>).

Forecasts of the improvements in terms of parameter constraints induced by adding to the data vector the galaxy bispectrum anisotropic signal were presented in Section <ref>. These have been obtained by running standard MCMC samplings together with using the analytical covariance matrix to evaluate the likelihood.
For a BOSS-like survey, we found on average ∼ 22 % tighter 1D 68% credible regions by just considering the two bispectrum quadrupoles in addition to the monopole (Figure <ref> and Table <ref>).

Including the last galaxy bispectrum hexadecapole term B_g^(22) only returns a further smaller improvement, suggesting that the constraining power of the bispectrum saturates at this level when expanded in terms of Legendre polynomials. This hypothesis is supported by the behaviour  of the  signal to noise ratio in Figure <ref> and 1-2D marginalised posterior distributions (again Figure <ref> and Table <ref>), where adding the third hexadecapole term, B_g^(22), does not provide any significant information gain.
Improving by ∼ 30% on average the parameter constraints by employing the bispectrum anisotropic signal approximately corresponds to observing a  survey volume two times larger for a BOSS-like survey.

We stress that these results are conservative since we restrict our analysis to linear scales k_max =0.09 h/Mpc. Cumulative signal to noise studies in the literature indicate that the benefits of both employing the anisotropic bispectrum <cit.> and adding the bispectrum to the power spectrum analyses <cit.>, increase for larger k_max.

In Appendix <ref> we report the key Figures displaying the results obtained when the same analysis done in the main body of the paper is also performed on the bispectrum data vector obtained for a larger k-bin width (Δ k_6), hence with less triangle configurations.
Even if the overall improvement on the parameters constraints appears to be larger, information gain saturation is again reached at the bispectrum hexadecapole level. Moreover, comparing the 1- and 2-D marginalised posterior distributions for the two k-bin width cases, Figure <ref> for Δ k_5 and Figure <ref> for Δ k_6, we appreciate that  the contours for the full data vector including all the multipoles up to the hexadecapole level are qualitatively very similar.

We  therefore conclude that the (significant) effort necessary to include the bispectrum quadruple in clustering analysis of galaxy redshift surveys is very likely to pay off.  We envision that this extra information could be used not only to reduce parameter errors; these additional statistics can provide a useful consistency checks of the modelling itself but possibly also of the model.





§ COVARIANCE MATRIX TERMS DERIVATION

In this Appendix are reported the full derivations of the covariance terms used in this work. We build on and improve what already presented in a previous paper <cit.>. In particular we are interested in computing the covariance expressions for both power spectrum and bispectrum expansion in terms of Legendre polynomials encoding the relation between the signal and the k-vectors orientation with respect to the line of sight. 





 §.§ Power spectrum multipoles auto-covariance

In order to compute the covariance of the power spectrum multipoles it is necessary to evaluate the four-point correlator in Fourier space:


    ⟨δ_g(𝐤_1)δ_g(-𝐤_1)δ_g(𝐤_2)δ_g(-𝐤_2) ⟩   = ⟨δ_g(𝐤_1)δ_g(-𝐤_2)⟩_c⟨δ_g(𝐤_2)δ_g(-𝐤_1)⟩_c + 1 perm.
       +⟨δ_g(𝐤_1)δ_g(𝐤_2)δ_g(-𝐤_1)δ_g(-𝐤_2)⟩_c
    ,

where the term in the  second line is proportional to the trispectrum for parallelogram configurations. Since we limit our analysis to linear scales, we only use the first Gaussian term, which dominates over the non-Gaussian one. For a recent work regarding the power spectrum covariance matrix see <cit.>. 
It is possible to compute the Gaussian term of the covariance for the power spectrum monopole and quadrupole as follows:


    C_G^P^(ℓ)_gP^(ℓ)_g(k_1;k_2) = 
    (2ℓ+1)^2(2π)^-6N_p(k_1)N_p(k_2)∫_V_𝐪_1∫_V_𝐪_2∫_V_𝐩_1∫_V_𝐩_2
     d^3𝐪_1 d^3𝐪_2 d^3𝐩_1 d^3𝐩_2
     
       ×ℒ_ℓ(μ_1)ℒ_ℓ(μ_2)
    δ_D(𝐪_1+𝐩_1)δ_D(𝐪_2+𝐩_2) 
    2(2π)^6δ_D(𝐪_1+𝐪_2)δ_D(𝐩_1+𝐩_2)
    P_g^s(𝐪_1)P_g^s(𝐩_2)
    
       =
    (2ℓ+1)^22N_p(k_1)N_p(k_2)∫_V_𝐪_1∫_V_𝐪_2
     d^3𝐪_1 d^3𝐪_2
    ℒ_ℓ(μ_1)ℒ_ℓ(μ_2)
    δ_D(𝐪_1+𝐪_2)^2
    P_g^s(𝐪_1)P_g^s(𝐪_2)
    
       =
    (2ℓ+1)^22N_p(k_1)V_p(k_2)∫_V_𝐪_1∫_V_𝐪_2
     d^3𝐪_1 d^3𝐪_2
    ℒ_ℓ(μ_1)ℒ_ℓ(μ_2)
    δ_D(𝐪_1+𝐪_2)
    P_g^s(𝐪_1)P_g^s(𝐪_2)
    
       ≈(2ℓ+1)^22k_f^3V_p(k_1)  12∫^1_-1dμ_1 ℒ_ℓ(μ_1)^2
    P_g(k_1,μ_1)^2
    
       =
    (2ℓ+1)^2
    2k_f^3δ^K_124π k_1^2Δ k  12∫^1_-1dμ_1 ℒ_ℓ(μ_1)^2
    P_g(k_1,μ_1)^2,


where again we used the approximation made in <cit.> that δ_D^2≈V_s(2π)^3δ_D = k_f^-3δ_D. V_s is the survey volume k_f the fundamental frequency.  δ^K_12 is the Kronecker delta indicating that the vector 𝐪_1 and 𝐪_2 are identical (in the second step trivial δ_K have been omitted in order to avoid making the notation heavier by adding also the wave-vector letter).  The Fourier integration volume for the power spectrum is defined as V_p = 4πΔ k k^2 and therefore the number of modes ("pairs") as N_p = V_p / k_f^3. L denotes the Lagrange polynomials. In the last steps we made the approximation that the power spectrum monopole and quadrupoles do not significantly vary when integrated over the bin in Fourier space.



 §.§ Bispectrum multipoles auto-covariance


We write these expression as functions of the multipoles of two angles defined in <cit.>. The first is the angle ω between the normal to the surface of the triangle and the line of sight. The second is the angle ϕ between the first k-vector and one of the two axis defining the plane on which the triangle lies. We will then use the cosines for the multipoles μ = cosω and ν = cosϕ. The multipoles indexes will be given then by ℒ^(α)_μ_a, ℒ^(β)_ν_a for the first triangle and ℒ^(κ)_μ_b, ℒ^(λ)_ν_b for the second one. We introduce the short notation ℳ^(αβκλ)_μν,ab =ℒ^(α)_μ_aℒ^(β)_ν_aℒ^(κ)_μ_bℒ^(λ)_ν_b. In the same way we  will write the compact pre-factor 𝒞^αβ_κλ = (2α+1)(2β+1)(2κ+1)(2λ+1).

Following the angles definition in <cit.>, the conversion from the standard bispectrum coordinates (k_1,k_2,k_3,μ_1,μ_2) to the set used in this work  (k_1,k_2,k_3,μ,ν) is given by


    μ = cosω   = √(2(a_⊥^2b_⊥^2 + a_⊥^2c_⊥^2 + b_⊥^2c_⊥^2) - a_⊥^4 - b_⊥^4 - c_⊥^42(a^2b^2 + a^2c^2 + b^2c^2)-a^4-b^4-c^4) = √(Δ_⊥Δ)
    ν = cosϕ    = (1-μ_a^2)μ_b + cosα_abμ_a√(1-μ_a^2)√(1-μ_b^2)√(1-μ^2) ,


where for notation simplicity we used (a,b,c) instead of (k_1,k_2,k_3). The '⊥' indicates the 2D vector obtained by projecting the 3D k-vector onto to the plane perpendicular to the line of sight. α_ab is the internal angle between the sides a_⊥ and b_⊥ of the a_⊥ b_⊥ c_⊥ triangle. When |μ| = 1 the triangle normal to the surface is either parallel or anti-parallel to the line of sight. Therefore the anisotropic bispectrum signal is the same for all the potential values of ν since the triangle plane coincides with the one perpendicular to the line of sight.


In order to compute the covariance between two galaxy bispectrum multipoles, the six-points correlator in Fourier space given by the respective galaxy density variables needs to be evaluated <cit.>:



    ⟨   δ_g(𝐤_1)δ_g(𝐤_2)δ_g(𝐤_3)δ_g(𝐤_4)δ_g(𝐤_5)δ_g(𝐤_6) ⟩
       = ⟨δ_g(𝐤_1)δ_g(𝐤_4)⟩_c⟨δ_g(𝐤_2)δ_g(𝐤_5)⟩_c⟨δ_g(𝐤_3)δ_g(𝐤_6)⟩_c + 8 perms.
       +⟨δ_g(𝐤_1)δ_g(𝐤_2)δ_g(𝐤_4)⟩_c⟨δ_g(𝐤_3)δ_g(𝐤_5)δ_g(𝐤_6)⟩_c + 8 perms.
       +⟨δ_g(𝐤_1)δ_g(𝐤_2)δ_g(𝐤_5)δ_g(𝐤_6)⟩_c⟨δ_g(𝐤_3)δ_g(𝐤_4)⟩_c+ 8 perms.
       +⟨δ_g(𝐤_1)δ_g(𝐤_2)δ_g(𝐤_3)δ_g(𝐤_4)δ_g(𝐤_5)δ_g(𝐤_6)⟩_c.


From the above expression the bispectrum covariance results composed by four terms. The first is the so called Gaussian term which being proportional to the cube of the power spectrum represents the lowest order in the perturbation expansion (∝δ_m^6 where "m" indicates the matter density field). The second and third terms are the lowest order non-Gaussian contributions (∝δ_m^8), proportional to the product of two bispectra and power spectrum times trispectrum, respectively.
The last is proportional to the pentaspectrum and since it is an higher order term than the previous ones (∝δ_m^10), it will not be computed explicitly. 





  
Gaussian term:
Here we update and extend  the work presented in <cit.>.
Analogously to what done  above, we now compute the diagonal term of the bispectrum multipoles covariance matrix:


    C_∝P^3^B^αβ_gB^κλ_g(k_1,k_2,k_3;k_4,k_5,k_6) = 
    
       =
    (2π k_f)^-6𝒞^αβ_κλN^t_123N^t_456∏^6_i=1∫_V_𝐪_i d^3𝐪_i
    ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_3)δ_D(𝐪_4+𝐪_5+𝐪_6)
    
       ×
    (2π)^9δ_D(𝐪_1+𝐪_4)δ_D(𝐪_2+𝐪_5)δ_D(𝐪_3+𝐪_6)
    P_g^s(𝐪_1)P_g^s(𝐪_2)P_g^s(𝐪_3)   + 5  perm.
       =
    D^123_456(2π)^3 k_f^-6𝒞^αβ_κλN^t,2_123∏^3_i=1∫_V_𝐪_i d^3𝐪_i
    ℳ^(αβκλ)_μν,aaδ_D(𝐪_1+𝐪_2+𝐪_3)^2
    P_g^s(𝐪_1)P_g^s(𝐪_2)P_g^s(𝐪_3)
    
       =
    D^123_456(2π)^3 k_f^-9𝒞^αβ_κλN^t,2_123∏^3_i=1∫_V_𝐪_i d^3𝐪_i
    ℳ^(αβκλ)_μν,aaδ_D(𝐪_1+𝐪_2+𝐪_3)
    P_g^s(𝐪_1)P_g^s(𝐪_2)P_g^s(𝐪_3)
    
       =
    D^123_456(2π)^3 k_f^3𝒞^αβ_κλV^t,2_123∏^3_i=1∫_V_𝐪_i d^3𝐪_i
    ℳ^(αβκλ)_μν,aaδ_D(𝐪_1+𝐪_2+𝐪_3)
    P_g^s(𝐪_1)P_g^s(𝐪_2)P_g^s(𝐪_3)
    
       ≈D^123_456 (2π)^3k_f^3𝒞^αβ_κλV^t_12314π∫^+1_-1dμ_1∫^2π_0dϕℳ^(αβκλ)_μν,aaP_g(k_1,μ_1)P_g(k_2,μ_2)P_g(k_3,μ_3),


where D^123_456 stands for all the possible permutations and has values 6,2,1 respectively for equilateral, isosceles and scalene triangles.
 V^t_123 = 8π^2k_1k_2k_3Δ k_1Δ k_2Δ k_3 is the integration volume in Fourier space for the bispectrum. The number of modes is then defined as N^t_123=V^t_123/k_f^6.
Again,  we have  assumed that the power spectrum monopole does not vary significantly inside the integration volume. The angle ϕ is defined such as μ_2≡μ_1cosψ_12 - √(1-μ_1^2)√(1-cosψ_12^2)cosϕ where cosψ_12≡(𝐤_1·𝐤_2/(k_1k_2). In order to be sure of the symmetry of the numerical result we integrate for each possible combination of k-vectors (k_1,k_2), (k_1,k_3) and (k_2,k_3).







  
non-Gaussian terms:
the following term (Bispectrum cross covariance ∝ BB) was also presented in <cit.> and allows one to better capture the correlation between different triangle bins. It is one of the two non-Gaussian terms that appear in the bispectrum self-covariance matrix up to order ∝δ_m^8.


    C_∝BB^B^αβ_gB^κλ_g(k_1,k_2,k_3;k_4,k_5,k_6) = 
    
       =
    (2π k_f)^-6𝒞^αβ_κλN^t_123N^t_456∏^6_i=1∫_V_𝐪_i d^3𝐪_i
    ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_3)δ_D(𝐪_4+𝐪_5+𝐪_6)
    
       ×
    (2π)^6δ_D(𝐪_1+𝐪_2+𝐪_4)δ_D(𝐪_3+𝐪_5+𝐪_6)
    B_g^s(𝐪_1,𝐪_2,𝐪_4)B_g^s(𝐪_3,𝐪_5,𝐪_6)   + 8  perm.
       =
    k_f^-6𝒞^αβ_κλδ^K_34N^t_123N^t_456∫_V_𝐪_1∫_V_𝐪_2∫_V_𝐪_3∫_V_𝐪_5∫_V_𝐪_6
    d^3𝐪_1d^3𝐪_2d^3𝐪_3d^3𝐪_5d^3𝐪_6
     ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_3)
    
       ×δ_D(𝐪_3+𝐪_5+𝐪_6)^2
    B_g^s(𝐪_1,𝐪_2,𝐪_3)B_g^s(𝐪_3,𝐪_5,𝐪_6)   + 8  perm.
       =
    k_f^-9𝒞^αβ_κλδ^K_34N^t_123N^t_456∫_V_𝐪_1∫_V_𝐪_2∫_V_𝐪_3∫_V_𝐪_5∫_V_𝐪_6ℳ^(αβκλ)_μν,ab
    d^3𝐪_1d^3𝐪_2d^3𝐪_3d^3𝐪_5d^3𝐪_6  δ_D(𝐪_1+𝐪_2+𝐪_3)
    
       ×δ_D(𝐪_3+𝐪_5+𝐪_6)
    B_g^s(𝐪_1,𝐪_2,𝐪_3)B_g^s(𝐪_3,𝐪_5,𝐪_6)   + 8  perm.
       =
    k_f^3𝒞^αβ_κλδ^K_34V^t_123V^t_456∫_V_𝐪_1∫_V_𝐪_2∫_V_𝐪_5∫_V_𝐪_6ℳ^(αβκλ)_μν,ab
    d^3𝐪_1d^3𝐪_2d^3𝐪_5d^3𝐪_6  δ_D(𝐪_1+𝐪_2-𝐪_5-𝐪_6)
    
       ×B_g^s(𝐪_1,𝐪_2,-𝐪_1-𝐪_2)B_g^s(-𝐪_5-𝐪_6,𝐪_5,𝐪_6)   + 8  perm.
       ≈k_f^3𝒞^αβ_κλδ^K_34V^t_123V^t_456×V^q_12568∫_-1^+1dμ_3dμ_1dμ_5 
    ℳ^(αβκλ)_μν,abB_g^s(k_1,k_2,k_3,μ_3,μ_1)B_g^s(k_3,k_5,k_6,μ_3,μ_5)   + 8  perm.
       =
    k_f^3𝒞^αβ_κλδ^K_344π k_3k_4Δ k×18∫_-1^+1dμ_3dμ_1dμ_5 
    ℳ^(αβκλ)_μν,abB_g^s(k_1,k_2,k_3,μ_3,μ_1)B_g^s(k_3,k_5,k_6,μ_3,μ_5)   + 8  perm.
    ,


where the usual approximations have been used together with Equation (<ref>); in the last step these approximations have  enabled us to simplify the quadrilateral and triangles volumes ratio. Notice that in this case the k-vectors 𝐤_3 and 𝐤_4 coincide, given the implicit δ_D(𝐪_3-𝐪_4).

The non-Gaussian term proportional to the trispectrum monopole (∝ PT) presented below is original of this work.


    C_∝TP^B^αβ_gB^κλ_g(k_1,k_2,k_3;k_a,k_b,k_c) = 
    
       =
    (2π k_f)^-6𝒞^αβ_κλN^t_123N^t_abc∏^3,c_i=1,a∫_V_𝐪_i d^3𝐪_i
    ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_3)δ_D(𝐪_a+𝐪_b+𝐪_c)
    
       ×
    (2π)^6δ_D(𝐪_1+𝐪_2+𝐪_a+𝐪_b)δ_D(𝐪_3+𝐪_c)
    T_g^s(𝐪_1,𝐪_2,𝐪_a,𝐪_b)P_g^s(𝐪_3)   + 8  perm.
       =
    δ^K_3ck_f^-6𝒞^αβ_κλN^t_123N^t_abc∫ d^3𝐪_1d^3𝐪_2d^3𝐪_3d^3𝐪_ad^3𝐪_b
    ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_a+𝐪_b)δ_D(𝐪_1+𝐪_2+𝐪_3)
    
       ×δ_D(𝐪_a+𝐪_b-𝐪_3)
    T_g^s(𝐪_1,𝐪_2,𝐪_a,𝐪_b)P_g^s(𝐪_3)   + 8  perm.
       =
    δ^K_3cδ^K_|1+2|3k_f^-6𝒞^αβ_κλN^t_123N^t_abc∫ d^3𝐪_1d^3𝐪_2d^3𝐪_ad^3𝐪_b
    ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_a+𝐪_b)^2
    
       ×T_g^s(𝐪_1,𝐪_2,𝐪_a,𝐪_b)P_g^s(-𝐪_1-𝐪_2)   + 8  perm.
       =
    δ^K_3ck_f^-9𝒞^αβ_κλN^t_123N^t_abc∫ d^3𝐪_1d^3𝐪_2d^3𝐪_ad^3𝐪_b
    ℳ^(αβκλ)_μν,abδ_D(𝐪_1+𝐪_2+𝐪_a+𝐪_b)
    
       ×T_g^s(𝐪_1,𝐪_2,𝐪_a,𝐪_b)P_g^s(-𝐪_1-𝐪_2)   + 8  perm.
       ≈δ^K_3ck_f^3𝒞^αβ_κλV^t_123V^t_abc×V^q_12ab8∫^+1_-1dμ_Ddμ_1dμ_b
    T_g(k_1,k_2,k_a,k_b,D,μ_D,μ_1,μ_b)P_g(D,-μ_D)
      + 8  perm.
       =
    δ^K_3ck_f^3𝒞^αβ_κλΔ D4π k_3k_cΔ k_3Δ k_c×18∫^+1_-1dμ_Ddμ_1dμ_b
    T_g(k_1,k_2,k_a,k_b,D,μ_D,μ_1,μ_b)P_g(D,-μ_D)
      + 8  perm. ,


where δ^K_|1+2|3 is automatically satisfied and the quadrilateral diagonal is D≡ k_3,k_c. In the last passage we used the expression for the trispectrum integration volume in Fourier space given in Equation (<ref>). Once more we made the approximation that the trispectrum and power spectrum monopoles do not significantly vary inside the integration volume. We choose to integrate out the angular dependence using the cosine of the angles made by the quadrilateral diagonal (shared triangle side) and two of the quadrilateral sides with respect to the line of sight. 
Notice that in this case the k-vectors 𝐤_3 and 𝐤_c are opposite given the implicit δ_D(𝐪_3+𝐪_c).



 §.§ Power spectrum - bispectrum cross-covariance


Analogously to  Equation <ref> for the bispectrum  covariance, the cross-covariance between power spectrum and bispectrum multipoles is obtained expanding the five-points correlator in Fourier space:


    ⟨δ_g(𝐤_1)δ_g(𝐤_2)δ_g(𝐤_3)δ_g(𝐤_4)δ_g(𝐤_5) ⟩   =
    ⟨δ_g(𝐤_1)δ_g(𝐤_3)⟩_c⟨δ_g(𝐤_2)δ_g(𝐤_4)δ_g(𝐤_5)⟩_c   +    5 perm.
       +⟨δ_g(𝐤_1)δ_g(𝐤_2)δ_g(𝐤_3)δ_g(𝐤_4)δ_g(𝐤_5)⟩_c


The above expansion returns a  term proportional to the product between power spectrum and bispectrum (∝δ^6_m) and a higher order term one proportional to the tetraspectrum (∝δ^8_m).
Therefore, at lowest order, we compute  the covariance between the power spectrum and bispectrum multipoles in the following way:


    C_∝PB^P^ℓ_gB^αβ_g(k_1;k_2,k_3,k_4) = 
    
       =(2π)^-6 k_f^-3𝒞_αβ^ℓN^p_1N^t_234∏^4_i=1∫_V_𝐪_i d^3𝐪_i∫_V_𝐩_1 d^3𝐩_1
    ℳ^(ℓαβ)_μ_aμ_bν_bδ_D(𝐩_1+𝐪_1)
    δ_D(𝐪_2+𝐪_3+𝐪_4)
    
       ×
    2(2π)^6
    δ_D(𝐪_1+𝐪_2)δ_D(𝐩_1+𝐪_3+𝐪_4)
    P_g^s(𝐪_1)
    B_g^s(𝐩_1,𝐪_3,𝐪_4)   + 2  perm.
       =2 k_f^-3𝒞_αβ^ℓδ^K_12N^p_1N^t_234∏^4_i=2∫_V_𝐪_i d^3𝐪_i
    ℳ^(ℓαβ)_μ_aμ_bν_bδ_D(𝐪_2+𝐪_3+𝐪_4)^2
    P_g^s(𝐪_2)
    B_g^s(𝐪_2,𝐪_3,𝐪_4)   + 2  perm.
       =2 k_f^-6𝒞_αβ^ℓδ^K_12N^p_1N^t_234∏^4_i=2∫_V_𝐪_i d^3𝐪_i
    ℳ^(ℓαβ)_μ_aμ_bν_bδ_D(𝐪_2+𝐪_3+𝐪_4)
    P_g^s(𝐪_2)
    B_g^s(𝐪_2,𝐪_3,𝐪_4)   + 2  perm.
       ≈2𝒞_αβ^ℓδ^K_12N^p_114π∫^+1_-1∫^2π_0dϕℳ^(ℓαβ)_μ_aμ_bν_bP_g^s(k_2,μ_2)
    B_g^s(k_2,k_3,k_4,μ_2,ϕ)   + 2  perm.
       =
    2 k_f^3𝒞_αβ^ℓδ^K_124π k_2^2 Δ k14π∫^+1_-1∫^2π_0dϕℳ^(ℓαβ)_μ_aμ_bν_bP_g^s(k_2,μ_2)
    B_g^s(k_2,k_3,k_4,μ_2,ϕ)   + 2  perm.


Notice that the Δ k appearing at the denominator comes from the number of modes available in Fourier space for the power spectrum, which can be defined differently from the one relative to the bispectrum integration volume. For the power, spectrum the multipole expansion is done in terms of μ_a, the cosine of the angle between the k-vector and the line of sight. The angle ϕ is the same as the one defined below Equation (<ref>).




§ INTEGRATION VOLUMES AND NUMBER OF MODES


In this Appendix we expand the work presented in Appendix A of <cit.>.








 §.§ Trispectrum integration volume and grid points quadruplets


The key to compute the integration volume for the trispectrum consists in considering one of the two diagonals. By doing so each trispectrum configuration can be described by 8 variables (12 degrees of freedom minus 3 for the Dirac's delta closure condition minus 1 for the invariance with respect rotations around the line of sight). A good choice for these are the four sides modules (k_1,k_2,k_3,k_4), one of the two diagonals D (let's choose the one defining the triangles k_1k_2D and k_3k_4D with closure conditions δ_D(𝐤_1+𝐤_2-𝐃) and δ_D(𝐤_3+𝐤_4+𝐃) respectively) and the cosine of three angles with respect to the line of sight (μ_D,μ_1,μ_4) for the diagonal 𝐃 and the two sides 𝐤_1, 𝐤_4.  For a visual representation of the choice of coordinates see Figure <ref>.

All the other details can be obtained using as coordinates the parallel/perpendicular components of the k-vectors with respect to  the line of sight. In this way, to fully define the vector, we must specify also the azimuthal angle on the plane perpendicular to the line of sight. This allows  us to fully exploit the invariance symmetry related to the rotation of the quadrilateral around the line of sight. This is done by setting the origin of the azimuthal angle  to be in the diagonal  (γ_D=0) and by computing the other four azimuthal angles with respect to  it, using the condition that the projection on the perpendicular plane must also be a closed quadrilateral composed by two closed triangles.

Considering the two triangles mentioned above, we start from  the diagonal vector 𝐃 yielding the initial volume V_D≃ 4π D^2Δ D.
Now we can proceed for each triangle as in the case of the bispectrum integration volume (see Appendix of <cit.>) obtaining the symmetric final result:


    V_q = V_D × V_k_12× V_k_34 =16π^3 k_1k_2k_3k_4Δ k_1Δ k_2Δ k_3Δ k_4 Δ D  .


The number of grid points quadruplets is then defined as N_q=V_qk_f^9.



 §.§ Tetraspectrum integration volume and grid points quintuplets


For future reference we also include the integration volume for the tetraspectrum. This can be derived using the same trick of splitting the figure into triangles made by the sides and the diagonals. A particular pentagon configuration is described by 11 variables (15 degrees of freedom minus 3 for the Dirac's delta closure condition minus 1 for the invariance with respect rotations around the line of sight). 

Once more we choose the 5 sides modules (k_1,k_2,k_3,k_4,k_5) and two diagonals D_A and D_B such that the resulting triangles are D_Ak_2k_3, D_Ak_1D_B and D_Bk_4k_5. 
Finally we use the four cosine of the angles with respect to the line of sight (μ_D_A,μ_D_B,μ_2,μ_4).

The volume can be computed starting from the one for the quadrilateral D_Bk_1k_2k_3 which is:


    V^B123_q = 16π^3D_Bk_1k_2k_3 Δ D_BΔ k_1Δ k_2Δ k_3 .


Again applying the bispectrum formula for the triangle D_Bk_4k_5 we have that:


    V_k_45 = 2πk_4k_5D_BΔ k_4Δ k_5 ,


and therefore the tetraspectrum integration volume V_c ('c' stands for cinco/cinque) can be obtained by simply multiplying


    V_c = V^B123_q× V_k_45 = 32π^4k_1k_2k_3k_4k_5Δ k_1Δ k_2Δ k_3Δ k_4Δ k_5Δ D_A Δ D_B .


The number of grid points quintuplets is then defined as N_c=V_ck_f^12.



§ SURVEY WINDOW EFFECT ON THE COVARIANCE MATRIX



In this Appendix we investigate the origin of the features parallel to the main diagonal of the numerically estimated covariance matrix in Figure <ref>. For this purpose we use the measurements of the same data-vector (power spectrum monopole and quadrupole plus bispectrum monopole) made on simulated cubic boxes with periodic boundary conditions <cit.>. Because of the fewer number of measurements available (440 with respect to 1400), the bispectrum data-vector contains 47 triangle configurations, less than the 68 employed in Figure <ref>.

While the overall discrepancy between analytical and numerical estimates is approximately the same as in Figure <ref>, in the top-right panel of Figure <ref> there are no features parallel to the main diagonal. This difference between Figure <ref> and Figure <ref> supports our hypothesis that those features are not physical but related to the mode-coupling induced by the survey window selection function.







§  BINNING CASE FIGURES


This last Appendix contains the two key figures and improvements table for the Δ k_6 binning case.


















D.G. thanks Héctor Gil Marín for providing the data vector's measurements from the galaxy mocks and for useful discussions. D.G. and L.V. are very grateful to Cheng Zhao and Chia-Hsun Chuang for the periodic simulations used to test the survey window effect on the covariance. L.V.  and D.G. acknowledge support of European Unions Horizon 2020 research and innovation programme ERC (BePreSySe, grant agreement 725327). Funding for this work was partially provided by the Spanish Ministerio de Ciencia y Innovation y Universidades under project PGC2018-098866-B-I00.
















mnras


 























