# The clustering of galaxies in the completed SDSS-IV extended Baryon
Oscillation Spectroscopic Survey: Primordial non-Gaussianity in Fourier Space

Eva-Maria Mueller1,9, Mehdi Rezaie2, Will J. Percival3,4,5, Ashley J. Ross6,
Rossana Ruggeri7, Hee-Jong Seo2, Héctor Gil-Marín8, Julian Bautista9,10, Joel
R. Brownstein11, Kyle Dawson11, Axel de la Macorra12, Nathalie Palanque-
Delabrouille13, Graziano Rossi14, Donald P. Schneider15,16, Christophe Yèche13
1 Department of Physics, University of Oxford, Denys Wilkinson Building, Keble
Road, Oxford OX1 3RH, UK
2 Department of Physics and Astronomy, Ohio University, 251B Clippinger Labs,
Athens, OH 45701, USA
3 Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave
W, Waterloo, ON N2L 3G1, Canada
4 Department of Physics and Astronomy, University of Waterloo, 200 University
Ave W, Waterloo, ON N2L 3G1, Canada
5 Perimeter Institute for Theoretical Physics, 31 Caroline St. North,
Waterloo, ON N2L 2Y5, Canada
6 Center for Cosmology and Astro-Particle Physics, Ohio State University,
Columbus, Ohio, USA
7 School of Mathematics and Physics, University of Queensland, Brisbane, QLD
4072, Australia
8ICC, University of Barcelona, IEEC-UB, Martí i Franquès, 1, E-08028
Barcelona, Spain
9 Institute of Cosmology and Gravitation, Dennis Sciama Building, University
of Portsmouth, Portsmouth PO1 3FX, UK
10 Aix Marseille University, CNRS/IN2P3, CPPM, Marseille, France
11 Department of Physics and Astronomy, University of Utah, 115 S. 1400 E.,
Salt Lake City, UT 84112, USA
12Instituto de Física, Universidad Nacional Autónoma de México, Apdo. Postal
20-364, México
13IRFU, CEA, Centre d’Etudes Saclay, 91191 Gif-Sur-Yvette Cedex, France
14 Department of Astronomy and Space Science, Sejong University, 209,
Neungdong-ro, Gwangjin-gu, Seoul, South Korea
15Department of Astronomy and Astrophysics, The Pennsylvania State University,
University Park, PA 16802
16Institute for Gravitation and the Cosmos, The Pennsylvania State University,
University Park, PA 16802
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present measurements of the local primordial non-Gaussianity parameter
$f_{\mathrm{NL}}^{\mathrm{local}}\ $from the clustering of 343,708 quasars
with redshifts 0.8 < z < 2.2 distributed over 4808 square degrees from the
final data release (DR16) of the extended Baryon acoustic Oscillation
Spectroscopic Survey (eBOSS), the largest volume spectroscopic survey up to
date. Our analysis is performed in Fourier space, using the power spectrum
monopole at very large scales to constrain the scale dependent halo bias. We
carefully assess the impact of systematics on our measurement and test
multiple contamination removal methods. We demonstrate the robustness of our
analysis pipeline with EZ-mock catalogues that simulate the eBOSS DR16 target
selection. We find $f_{\mathrm{NL}}=-12\pm 21$ (68% confidence) for the main
clustering sample including quasars with redshifts between 0.8 and 2.2, after
exploiting a novel neural network scheme for cleaning the DR16 sample and in
particular after applying redshift weighting techniques, designed for non-
Gaussianity measurement from large scales structure, to optimize our analysis,
which improve our results by 37%.

###### keywords:

cosmology: observations - inflation - large-scale structure of Universe

††pubyear: 2021††pagerange: The clustering of galaxies in the completed SDSS-
IV extended Baryon Oscillation Spectroscopic Survey: Primordial non-
Gaussianity in Fourier Space–The clustering of galaxies in the completed SDSS-
IV extended Baryon Oscillation Spectroscopic Survey: Primordial non-
Gaussianity in Fourier Space

##  1 Introduction

The aim of this paper is to constrain primordial non-Gaussianity (PNG) using
the large scale structure (LSS) in the clustering of quasars covering the
largest volume of the Universe surveyed to date. We use the quasars within the
final data release (DR16; Ahumada et al. 2020) of the extended Baryon
Oscillation Spectroscopic Survey (eBOSS; Dawson et al. 2016), which is a part
of the Sloan Digital Sky Survey (SDSS; Blanton et al. 2017) program.

PNG arise during inflation, a phase of rapid expansion in the very early
Universe, and can be used to differentiate between different underlying
physical models. Slow roll, single field inflation generates a primordial
gravitational potential that can be well approximated by a Gaussian random
field (Bardeen et al., 1986). Alternative inflationary models (e.g., multiple
field) predict there to be significant non-Gaussian components to the
potential (Bartolo et al., 2004). In this paper, we focus on the simplest case
of PNG, assuming local scale-independent PNG parameterized by
$f_{\mathrm{NL}}$. The tightest constraints on this form of PNG are currently
from the Planck satellite (Planck Collaboration et al., 2020) measuring the
cosmic microwave background (CMB), who found $f_{\mathrm{NL}}=-0.9\pm 5.1$.
However, the LSS in the clustering of galaxies or quasars offers an
alternative pathway to detect PNG that has the potential to improve upon these
constraints in the future (Doré et al., 2014; Yamauchi et al., 2014;
Ferramacho et al., 2014; Percival et al., 2019).

One effect of PNG is to add a scale-dependence to the large scale linear bias
of galaxies with respect to the underlying dark matter (Dalal et al. 2008;
Matarrese & Verde 2008; Slosar et al. 2008, and discussed further in Section
2). The unique $f_{\mathrm{NL}}/k^{2}$ dependence of the galaxy bias
introduced by PNG means that the signal increases to large scales, and
requires a large volume survey with exquisite control of systematic errors.
For any sample of galaxies, the bias and hence the PNG signal will change with
redshift. Weights were introduced by Mueller et al. (2019) to optimally
include all of the information across redshift and are used in our paper.
Including constraints from the bispectrum and trispectrum would further
improve these predictions (Karagiannis et al., 2018; Gualdi et al., 2021), but
we do not do this because of the increased complexity of such an analysis.

Previous measurements of PNG using LSS include Slosar et al. (2008); Ross et
al. (2013a); Giannantonio et al. (2014b); Ho et al. (2015); Leistedt et al.
(2014); Castorina et al. (2019). In general these studies are systematic
limited, with variations in observed large-scale galaxy density driven by such
effects as extinction, depth, air-mass or star density. Because the clustering
only needs to be measured on very large-scales, the use of photometric
redshifts does not significantly degrade the statistical error. However, our
ability to understand and correct for problems in the data is significantly
improved using spectroscopic redshifts.

Quasars offer an ideal large-scale structure tracer for measuring PNG because
they are highly biased (Laurent et al., 2017) and their brightness means they
can be observed over large volumes. The eBOSS quasar sample we use (Ross et
al., 2020a), includes objects up to a redshift of 2.2 over 4808 deg2, with a
physical volume of 68.7 Gpc3 (assuming our baseline $\Lambda$CDM model), and
is the largest LSS sample in volume to date. Further details are provided in
Section 3. This is a subset of the full sample of quasars observed by SDSS,
which is described in Lyke et al. (2020), and is subsampled in order to
provide a catalogue for clustering measurements for which we can model the
selection function. The large volume covered by the quasar catalogue is
crucial for the detection of local PNG through its effect on the large-scale
bias. Hence, while the eBOSS quasar sample was primarily designed (Zhao et
al., 2016) for the measurements of Baryon Acoustic Oscillations (BAO) and
Redshift-Space Distortions (RSD), it is also ideal to study PNG.

Castorina et al. (2019) measured $f_{\mathrm{NL}}$ using the eBOSS DR14 quasar
sample, finding $-51<f_{\mathrm{NL}}<21$ at 95 % confidence. Their sample was
taken when eBOSS was part-complete covered $2114$ square degrees, about half
of the DR16 sample we analyse. Our analysis also differs in that we present a
novel neural network based method to remove systematic density fluctuations in
this sample (Rezaie et al., 2021), which is summarised in Section 4. Standard
BAO and RSD measurements made with this sample are presented in Hou et al.
(2021); Neveux et al. (2020), together with a blind test of their methodology
in Smith et al. (2020). At redshifts $z>2.1$, eBOSS also made Lyman-$\alpha$
forest measurements using a high redshift quasar sample (du Mas des Bourboux
et al., 2020). The BAO measurement using the same neural network catalog we
implement in this paper is presented in Merz et al. (2021).

The eBOSS DR16 quasar sample (Lyke et al., 2020; Ross et al., 2020a) is the
largest sample of spectroscopically confirmed quasars available to date, and
is one of three target classes observed within eBOSS. Samples of Emission line
galaxies (Raichoor et al., 2020; Tamone et al., 2020; de Mattia et al., 2020)
and Luminous Red Galaxies (Bautista et al., 2020; Gil-Marín et al., 2020) were
also observed and used to make cosmological inferences (eBOSS Collaboration et
al., 2020).

As we are interested in very large scales, we have to be careful when
calculating the power spectrum and allowing for the window function and
integral constraint: our methods are described in Section 6 and tested using
mock catalogues (described in Section 5) as described in Section 7. We do not
include relativistic effects which, although they are degenerate with
$f_{\mathrm{NL}}$, are not expected to have an influence larger than
$f_{\mathrm{NL}}\pm 1$ on the scales and redshifts fitted (Wang et al., 2020).
Similarly, at the level of constraints that can be achieved with eBOSS data,
we do not allow for the non-Gaussian nature of the likelihood as discussed,
for example by Kalus et al. (2019); Wang et al. (2019). The latter assumption
is validated by our analyses of mock catalogues. We present the results from
our analysis in Section 8 and discuss them in Section 9.

##  2 Physical Model

In this Section we provide a brief summary of the mechanism by which local
non-Gaussianity causes scale dependent halo bias. Assuming a type of non-
Gaussianity that only depends on the local value of the potential, the
primordial potential can be parametrised (Gangui et al., 1994; Komatsu &
Spergel, 2001)

$\Phi=\phi+f_{\mathrm{NL}}(\phi^{2}-\langle\phi^{2}\rangle)\,,$ (1)

where $\phi$ is a Gaussian random field and $f_{\mathrm{NL}}$ describes the
amplitude of the quadratic non-Gaussian term in the potential. In Fourier
space, the potential is related to the density field by
$\delta(k)=\alpha(k)\Phi(k)$, with

$\alpha(k)=\frac{2k^{2}T(k)D(z)}{3\Omega_{m}}\frac{c^{2}}{H_{0}^{2}}\frac{g(0)}{g(\infty)}\,.$
(2)

$T(k)$ is the transfer function, $D(z)$ the linear growth factor normalised
with $D(0)=1$. $\Omega_{m}$ is the matter density and $H_{0}$ the Hubble
constant. The factor $g(\infty)/g(0)$, with $g(z)=(1+z)D(z)$, corrects for the
late-time normalisation of $D(z)$, and can be omitted if $D(z)$ is normalised
to equal the scale factor during the matter dominated era. Following Mueller
et al. (2019), we assume the CMB convention for $f_{\mathrm{NL}}$ assuming
$\Phi$ is the primordial potential. Some authors have previously adopted a
"LSS convention" that assumes $\Phi$ is extrapolated to $z=0$, with
$f_{\mathrm{NL}}^{\mathrm{LSS}}=g(\infty)/g(0)f_{\mathrm{NL}}^{\mathrm{CMB}}\approx
1.3\ f_{\mathrm{NL}}^{\mathrm{CMB}}$.

In the standard ansatz, large-scale halo bias depends on the derivative of the
halo mass function with respect to the large-scale matter overdensity field:
small changes in the overdensity can give rise to large changes in the density
of galaxies as patches of the Universe are pushed over a critical density for
collapse (Gunn & Gott, 1972). The non-Gaussian term in Eq. 1, coupled with the
$k$-dependence in Eq. 2 thus give rise to a scale dependent halo bias $\Delta
b(k)$ (Dalal et al., 2008; Slosar et al., 2008)

$\Delta b(k)=2(b-p)f_{\mathrm{NL}}\frac{\delta_{\mathrm{crit}}}{\alpha(k)}\,,$
(3)

where $\delta_{\mathrm{crit}}=1.686$ and $p$ is a correction due to galaxy
selection beyond a Poisson sampling of the halos of a given mass: quasars are
thought to result from galaxy mergers, and the merger rate is proportional to
$\sigma_{8}^{-1}$. For quasars, the standard ansatz requires us to determine
$\partial\ln n_{\mathrm{merger}}/\partial\ln\sigma_{8}$ rather than
$\partial\ln n_{\mathrm{halo}}/\partial\ln\sigma_{8}$, and this leads to an
extra component in $\Delta b$ that changes $p=1$ to
$p=1+1/\delta_{\mathrm{crit}}\simeq 1.6$ (Slosar et al., 2008). The total bias
of quasars, including local non-Gaussianity is modelled as
$b_{\mathrm{tot}}=b+\Delta b(k)$.

##  3 The eBOSS DR16 QSO dataset

Our quasar sample is derived from the final data release, (DR16; Ahumada et
al. 2020), of the SDSS-IV eBOSS survey. Over 4.5 years, eBOSS used the Sloan
Telescope (Gunn et al., 2006) and BOSS spectrographs (Smee et al., 2013) to
obtain spectra of quasar and galaxy targets.

The algorithms to robustly identify quasars and accurately measure their
redshifts from SDSS spectra were developed in Lyke et al. (2020), who produced
the DR16 quasar catalogue. This quasar catalogue was used as the input to the
main clustering sample, used for the BAO and RSD analysis and detailed in Ross
et al. (2020a). It includes quasars with $0.8<z<2.2$. The catalogues were
processed separately in the North and South Galactic Cap (NGC and SGC)
footprints. The NGC sample contains 218,209 quasars over 2,924 square degrees
and the SGC sample 125,499 quasars over 1,884 square degrees. The sample is
simulated with 1000 mock realizations, as described in Zhao et al. (2021).

The quasar sample contains strong spurious fluctuations in its projected
density on the sky (Myers et al., 2015; Ata et al., 2018; Ross et al., 2020a).
These are primarily related to imaging depth and the probabilistic nature of
the target selection that is intrinsically dependent on the precision of the
photometry; see Myers et al. (2015) and Section 4 for more details. In Ross et
al. (2020a), these fluctuations were treated by determining weights based on a
linear regression including maps of the $g$-band depth, the PSF size in the
$i$-band, the sky background in the $i$-band, and the amount of Galactic dust
extinction determined by Schlegel et al. (1998a). This level of correction was
determined sufficient for the BAO and RSD analyses (Hou et al. 2021; Neveux et
al. 2020; Smith et al. 2020), and is described in more detail in Section 4.1,
but a re-evaluation of the needs for an $f_{\rm NL}$ measurement was promised
to be presented here and in Rezaie et al. (2021). This is presented in Section
4.2.

Fortunately, such spurious fluctuations in quasar target density do not couple
strongly to changes in the redshift distribution. This can be seen in Fig. 1,
where the relative number of quasars as a function redshift is shown when
splitting the sample by various imaging properties. The normalized redshift
distributions do not appear to depend on the properties of the imaging sample.

##  4 Photometric Systematics

The density variations caused by variations in the properties of imaging data
typically result in an excess clustering signal on large scales (e.g.,
$k<0.01~{}h/{\rm Mpc}$). Unfortunately, the $k^{-2}$-dependent signature of
primordial non-Gaussianity is sensitive to the same scales, and thus a robust
constraint on the primordial non-Gaussianity parameter $f_{\rm NL}$ demands a
thorough and careful assessment of imaging systematics. This reality has been
recognized as one of the primary challenges in previous studies of PNG with
SDSS data (see, e.g., Ho et al., 2015).

Systematic effects can be classified in general into two categories:
obscuration and contamination. These effects are intermingled too. For
instance the effect of Galactic obscuration alters the color of objects and
can push some of the quasars out of the color-magnitude selection. Adversely,
it might also allow some astrophysical artifacts such as stars to be miss-
classified as quasars and enter the selection as contaminants (see, e.g.,
Crocce et al., 2015, for a discussion).

In what follows, we briefly discuss two map-based regression approaches for
the mitigation of imaging systematics. These approaches are computationally
efficient, however, their robustness and efficiency in eliminating systematics
relies on the available imaging maps and the flexibility of the model.

Figure 1: Histograms showing the quasar distribution $n(z)$ for the NGC as a
function of redshift for bin of high and low photometric properties. The bins
were chosen so that the number of quasars in each bin is equal. Nstar and
E(B-V) refer to the stellar density and extinction respectively, while
Depth-g, PSF-g and Sky-i stands for the depth of the imaging data, the seeing
and sky background in the g-band and i-band. We find no dependence of $n(z)$
on the photometric systematics, supporting our assertion that we can correct
for the systematics with a function that only depends on angular position.

###  4.1 Method I: Linear weights

The standard systematics treatment, which is implemented in Ross et al.
(2020b), uses a multivariate linear model to capture the response of the
observed quasar count map to the templates of imaging properties (e.g.,
Bautista et al., 2018). The templates include Galactic extinction from
Schlegel et al. (1998b), survey depth in the $g$ band (depth-g), sky
brightness in the $i$ band (sky-i), and seeing in the $i$ band (psf-i) from
SDSS imaging; and all are pixellated in HEALPix format with nside$=512$. These
maps are identified as primary imaging properties which have induced
systematic fluctuations in the DR16 quasar sample Ross et al. (2020b). In the
standard treatment, the impact of theses templates on the target sample
density is assumed to be linear, with four parameters - the amplitude of the
effect from each template. The parameters of the linear model are trained by a
least-square minimization that reduces the sum of squared error fluctuations
in the mean quasar density against imaging properties, including the four
primary templates as well as the SDSS stellar density and airmass.

###  4.2 Method II: Neural Networks

As an alternative, we rely on an ensemble of fully connected neural networks
(NN) to obtain the systematic weights as described in Rezaie et al. (2021).
Neural networks are capable of representing complex functions, and thus this
new set of systematics weights are able to eliminate nonlinear systematic
effects. This is very crucial for cleaning the DR16 sample, especially in the
NGC, where we observe a strong nonlinear trend against Galactic extinction.

We adapt neural networks to predict the projected quasar density given the
imaging data properties imprinted on angular maps pixellated using HEALPix
(Górski et al., 2005) with nside$=512$ as our default resolution, covering the
survey. In addition to the four maps used in Ross et al. (2020b), we find that
there is a strong dependence on the local stellar density, constructed from
Gaia DR2 (Gaia Collaboration et al., 2018) with the g-band magnitude cut
$12<{\rm gmag}<17$. We find that the inclusion of the Gaia stellar density was
necessary to obtain a small $\chi^{2}$ value consistent with the systematics-
free simulations. The importance of this map was only realised after the
public released with the catalogue, so it is not included in the default
linear weights described in Section 4.1, which used a stellar template
constructed from SDSS. We also confirm that the improvements we see cannot be
reproduced with this stellar template.

A fully connected neural network is constructed from many nonlinear units,
called neurons, organized on multiple layers. The information of imaging
properties flows forward starting from the first layer through the
intermediate layers, and out to the output layer as the predicted quasar
density. Each neuron calculates a linear combination of the outputs from the
previous layer neurons, and then applies a nonlinear activation function. The
arrangement of multiple nonlinear neurons on multiple layers allows the neural
network to have the flexibility to model a wide range of functions. Our
implementation uses three hidden layers with 20 Rectified Linear Unit neurons
on each hidden layer, which is defined as ${\rm ReLU}(x)=\max(x,0)$. The
output layer has a single neuron with a Softplus activation, i.e.,
Softplus(x)=$\log(1+\exp(x))$, to ensure that the predicted quasar count is
positive. The parameters of the neural network are optimized to increase the
similarity between the neural network output and the observed projected
density field by minimizing the Poisson negative log-likelihood. We
incorporate five-fold validation for the training, validation, and testing our
neural network. Specifically, we train the parameters with $60\%$ of the data,
use $20\%$ of the data to validate and find the best model with minimum
validation error. Next, we apply the best model on the remaining $20\%$ of the
sample, which is called the test set, and then use the inverse of the
predicted quasar count as $w_{\rm systot}$. By permutation of the training,
validation, and test sets, we test the methodology on the entire footprint and
obtain the systematic weights for the entire sample. Finally, the weights are
normalized such that the total weighted number of quasars stays invariant
before and after applying the weights.

To test our implementation, we have investigated the impact of changing the
resolution of imaging maps, splitting the quasar sample into $0.8<z<1.5$ and
$1.5<z<2.2$, and using various combinations of imaging maps. For each of these
tests, we characterize the residual systematic error with the $\chi^{2}$
statistics by calculating the angular cross-correlation of the projected
quasar density and imaging maps and the fluctuations in the mean density of
quasars against imaging properties. We find no evidence for redshift-dependent
imaging systematics, and no significant impact when training the networks with
templates in nside$=256$.

Fig. 2 illustrates the mean density of quasars versus three primary imaging
properties in the NGC (top) and SGC (bottom) regions. The survey $depth-g$
causes the strongest trend in both galactic caps with chi-squared values
around $266.94$ and $828.74$, respectively, in the NGC and SGC. The NGC sample
also indicates a nonlinear trend against the extinction with
$\chi^{2}=374.95$, which remains significant ($\chi^{2}=36.25$) even after the
linear treatment (red). On the other hand, our neural network approach (black)
is able to reduce the chi-squared error to $9.9$ (see Rezaie et al. 2021 for
more details).

Figure 2: Mean density of the DR16 quasars vs imaging properties before
imaging systematics correction (blue), after linear (red), and neural network
(black) treatments in the NGC (top) and SGC (bottom) regions. The residual
fluctuations are characterized using the $\chi^{2}$ statistics and covariance
matrix constructed from the EZmocks.

##  5 Mocks

To estimate the covariance matrix we use a set of 1000 mock catalogues
described in Zhao et al. (2021), called EZmocks. These synthetic clustering
catalogues include the survey geometry, redshift evolution and systematics,
and are based on the framework of Chuang et al. (2014) using the Effective
Zel’dovich approximation. The simulation assumes a $\Lambda$CDM cosmology of
$h=0.6777$, $\Omega_{m}=0.307115$, $\Omega_{\Lambda}=0.692885$,
$\Omega_{b}=0.048206$, $\sigma_{8}=0.8225$ and $n_{s}=0.9611$ in a cube of 5
$h^{-1}$Gpc comoving sidelength. The mocks assume zero $f_{\mathrm{NL}}\ $.

The EZmocks were made from seven simulation snapshots with the clustering
signal tuned to match the DR16 quasar catalogue. The EZmocks are fast,
approximate mocks that agree well with N-body simulations at large scales, but
show some discrepancy at small scales due to non-linear gravitational
interactions. These effects only impact scales above $k>0.25\,h/{\rm Mpc}$
which are not used in this study since we only consider only linear scales.
Thus we expect no problems for our analysis resulting from the approximations
made in creating these simulations.

To quantify the impact of systematic effects present in the data, we
introduced linear fluctuations in the distribution of tracers based on the
survey depth in the $g$ band, and the SDSS stellar density map. An angular map
was generated based on these two maps and the EZmocks were down-sampled
randomly by discarding mock tracers with a probability selected according to
this map. Thus, the EZmocks are artificially contaminated with linear
systematic trends. Additionally, the effects of redshift failure and fiber
collisions were introduced into the mocks (see Zhao et al. 2021 for details),
however these effects mainly play a role at much smaller scales then used in
our analysis. Thus, the contaminated mocks include the uncertainty due to the
stochasticity of the systematic treatment. We then assigned weights to the
contaminated mocks to account for these effects using the two techniques
described Section 4, but using only the two templates used to assign
systematics to the mocks. Thus, we can test the effectiveness of the neural
network based approaches in comparison to the linear regression method, and
the effect of the different approaches, although because they differ in
details, we cannot use the mocks to measure an absolute size of the effect of
contaminants or their removal.

##  6 Analysis Technique

In this section we outline the analysis technique used in this study. The
methods are similar to previous clustering studies in Fourier space, however
the treatment of the window function and integral constraint differ slightly
to ensure a higher precision at very large scales. Additionally, we apply
redshift weights designed to optimise constraints on primordial non-
Gaussianity that are not used in the standard BAO and RSD studies. We restrict
our analysis to linear scales and can therefore model the power spectrum
without non-linear corrections.

###  6.1 Power spectrum Estimation

We estimate the power spectrum following the approach of Feldman et al.
(1994), calculating the observed over-density field using a random catalogue
with the same survey geometry, completeness and redshift distribution as our
data to approximate the expected mean density of the quasar sample. We start
with the un-normalised over-density field $F(x)$,

$F(\mathbf{x})=n_{q}(\mathbf{x})w_{q}(\mathbf{x})-\alpha
n_{r}(\mathbf{x})w_{r}(\mathbf{x})$ (4)

where $n_{q}(\mathbf{x})$ is the quasar density on a 3D cartesian grid $x$,
$n_{r}(\mathbf{x})$ the density of the random catalog, and $\alpha=\sum
w_{\mathrm{d,noFKP}}/\sum w_{\mathrm{r,noFKP}}$ matches weighted densities in
the data and random catalogues and $w_{q}$ and $w_{r}$ are the weights that
correct for various observing effects as discussed in Section 3. The objects
in the catalogues are put on a grid of size $L=6600\,{\rm Mpc}/h$ and $N=512$
using the Cloud-in-Cell interpolation scheme (see Sefusatti et al. 2016 for
instance). We then determine the multipoles of the power spectrum by taking a
series of Fast Fourier Transform (FFTs) of $F(\mathbf{x})$ and applying the
estimator of Bianchi et al. (2015) and Scoccimarro (2015). The monopole, for
instance, is then calculated as

$P_{0}(k)=\frac{1}{A}\int\frac{d\Omega_{k}}{4\pi}\left[\hat{F}(\mathbf{k})\hat{F}^{*}(\mathbf{k})-S\right]$
(5)

where $A$ the normalisation and the short noise $S$ are given by

$\displaystyle A$ $\displaystyle=$ $\displaystyle\alpha\sum
w^{2}_{\mathrm{r,FKP}}\ n(z)\ w_{\mathrm{r,noFKP}}$ (6) $\displaystyle S$
$\displaystyle=$ $\displaystyle\sum w^{2}_{\mathrm{q}}+\alpha^{2}\sum
w^{2}_{\mathrm{r}}$ (7)

where the subscripts “q” and “r” refer to quasars and randoms, respectively.
Note that the values of $\alpha$, the normalisation $A$ and shot noise $S$
differ slightly from previous analyses due to the treatment of the
completeness in the random catalogue since the randoms are weighted instead of
down-sampled to match the data catalogue. We cut off our analysis at
$k=0.12\,h/{\rm Mpc}$, well below the Nyquist frequency of
$k_{\mathrm{Ny}}=0.24\,h/{\rm Mpc}$ and bin the power spectrum using $\Delta
k=2k_{f}=4\pi/L$. This procedure is the same as in the BAO and RSD studies in
Fourier space of the quasar Neveux et al. (2020), LRG Gil-Marín et al. (2020)
and ELG de Mattia et al. (2020) samples.

Figure 3 shows the power spectrum monopole for the NGC (left) and SGC (right)
regions after the neural network (black) and linear treatment (red). To
visualise the variance and expected mean at large scales, the power spectrum
result for the individual mocks is over-plotted in grey and the mean of the
mocks in green. We do not use the quadrupole in this analysis since it is very
noise at large scales for the DR16 QSO data.

Figure 3: Power spectrum monopole of the eBOSS DR16 QSO NGC (left) and SGC
(right) for different systematic treatments. The red line corresponds to the
standard linear treatment for systematics applied to the data, while the black
line represents the advanced neural network based systematic removal
technique. The mean of the mocks is displayed in green, while the individual
mocks are over-plotted in grey.

###  6.2 Window function and Integral constraint

The area of the survey leaves a footprint in the clustering of quasars and
needs to be taking into account when analysing two points statistics. In
Fourier space, the estimated power spectrum from the survey is not just the
quasar power spectrum but the power spectrum convolved with the survey window
function. We include this effect by convolving the model power spectrum with
the window function. We calculate the window function multipoles using random
pair-counts as outlined in Beutler et al. (2017) and Wilson et al. (2017) with
the normalisation of the window given by the same normalisation as the power
spectrum. The convolved power spectrum $\hat{P}(k)$ is then calculated from
the convolved correlation function $\hat{\xi}$ through

$\hat{P}_{l}=4\pi(-i)^{l}\int ds\ s^{2}\hat{\xi}_{l}(s)j_{l}(sk).$ (8)

Additionally, we have to account for the integral constraint. The formalism
described in Section 6.1 to estimate the over-density field $F(\mathbf{x})$
implicitly sets the mean quasar density to zero in our analysis. Taking a
Fourier transform of $F(\mathbf{x})$ then gives rise to an additional
contribution to the monopole through the convolution with the window function

$\left<F(\mathbf{k})F^{*}(\mathbf{k})\right>=\int\frac{dk^{\prime}}{(2\pi)^{3}}\left[\tilde{P}(k^{\prime})-\tilde{P}(0)\delta_{D}(k)\right]W(k)$
(9)

with the shot-noise uncorrected power spectrum $\tilde{P}$. As pointed out by
de Mattia & Ruhlmann-Kleider (2019), there is an additional contribution to
the integral constraint because the radial distribution of the randoms is
exactly matched to the data giving rise to an ‘radial’ contribution to the
integral constraint. Here we account for the radial integral constraint by
calculating the mean power spectrum of the mocks with the redshifts of the
randoms drawn for the data mocks, i.e. with the radial integral constraint
introduced into the mocks, and with the redshifts drawn from the
$n(z)$-distribution. The additional radial component is then calculated as the
fractional difference of the mean power spectra, $W_{\mathrm{RIC}}$. The total
integral constraint is then given by
$-\tilde{P}(0)W(k)-\tilde{P}(k)W_{\mathrm{RIC}}$.

###  6.3 Optimal Weights

To fully exploit all the information in the LSS of the quasar clustering we
apply redshift weights designed to minimise the uncertainty in the PNG
parameter $f_{\mathrm{NL}}\ $as derived in Mueller et al. (2018). The redshift
weighting technique was first developed for BAO measurements (Zhu et al.,
2015, 2016) and was subsequently also applied to optimised RSD measurements
(Ruggeri et al., 2017; Zhao et al., 2019). While the FKP weights balance shot
noise and cosmic variance contributions to minimise the statisical
uncertanity, the redshift weights take the redshift evolution of the
underlying physical theory into consideration to further improve the
constraints. This methods avoids a binning in redshift space and reduces the
information loss at the edges of the bins. In this study we only focus on the
power spectrum monopole and thus the redshift weight is defined as

$w_{z}^{2}(z)=(b(z)+\frac{1}{3}f(z))(b-p)D(z).$ (10)

with the fiducial bias given by $b(z)=0.53+0.29(1+z)^{2}$. We apply these
weights to each object in the data and randoms. The total weight is then

$w_{\mathrm{tot}}=w_{\mathrm{FKP}}\times w_{\mathrm{CP}}\times
w_{\mathrm{NOZ}}\times w_{z}\times w_{\mathrm{obs}}.$ (11)

The effective redshift of $z_{\mathrm{eff}}=(\sum w_{\mathrm{tot}}z)/\sum
w_{\mathrm{tot}}=1.51$ increases to $z_{\mathrm{eff}}=1.80$, reflecting the
up-weighting to higher redshift due to the redshift weights.

###  6.4 Power spectrum Model

We restrict our analysis to linear scales with the redshift power spectrum,
$P(k,\mu)$ modelled as

$P(k,\mu)=\frac{(b+\mu^{2}f)^{2}P_{m}(k)}{1+(k\mu\sigma_{v}/H_{0})^{2}}$ (12)

with the linear matter power spectrum $P_{m}(k)$ from CAMB, the galaxy bias
$b$, the growth rate $f$ and a Lorenzian damping factor given the velocity
dispersion $\sigma_{v}$. The power spectrum multipoles can then be calculated
as

$P_{\ell}(k)=\frac{2\ell+1}{2}\int_{-1}^{1}P(k,\mu)L_{\ell}(\mu)d\mu.$ (13)

Figure 4 shows the model of the power spectrum monopole (blue) using the best
fit parameters from the fit to the mean of the mocks. When the convolution
with the survey window function is included (orange) the power spectrum is
modified at large scales. If we additionally apply the integral constraint to
the model (green), we can recover the mean of the mocks (black dashed).

Figure 4: Theoretical model of the power spectrum monopole given the best fit
parameters to the mean of the mocks (blue), with the convolution of the window
function (orange) and with the integral constraints (green) applied to the
model. The black dashed line represents the mean of the mocks.

###  6.5 Covariance and Likelihood

For the main quasar clustering catalogue we estimate the covariance of the
power spectrum multipoles numerically from the 1000 EZ-mocks described in
Section 5. We invert the covariance including the Hartlap correction (Hartlap
et al., 2007), a 6.5% effect for our configuration. The constraints on
$f_{\mathrm{NL}}$ along with the nuisance parameters $b$ and $\sigma_{v}$ are
then inferred from maximizing the likelihood $\mathcal{L}\propto
e^{-\chi^{2}/2}$ using a Monte Carlo Markov Chain (MCMC) Metropolis-Hasting
algorithm with the $\chi^{2}$ given by
$\chi^{2}=(\mathbf{D}-\mathbf{T(p)})^{T}\mathbf{C}^{-1}(\mathbf{D}-\mathbf{T(p)})$
where $D$ refers to the data vector, $T(p)$ to the theory prediction and p the
set of parameters. We use flat priors on all parameters and to assure
convergence of the MCMC chains we apply the Gelma-Rubin statistics. All MCMC
chains have a converge of $R<0.01$ or lower. In our analysis we marginalised
over the shot noise S, the bias and the velocity dispersion $\sigma_{v}$ in
each cap separately as nuisance parameters. With the rest of the cosmological
parameters fixed, we vary four parameters in total in our MCMC runs.

##  7 Robustness tests with mocks

###  7.1 Testing the model

We test the reliability of our analysis pipeline, in particular, the modelling
of the window function and integral constraints, using the EZ-mocks described
in Section 5, with no observational systematics added, which we call ‘null’
mocks.

Figure 5: Robustness test I: Best-fit $f_{\mathrm{NL}}\ $(left) and
corresponding $\chi^{2}$ (right) from performing the MCMC analysis for 200
null mocks constructed without systematic errors and including optimal weights
in the analysis. The red dashed lines correspond to the best fit results from
the data.

Figure 5 shows the minimum $\chi^{2}$-distribution (right) and the best-fit
$f_{\mathrm{NL}}\ $(left) of 200 individual mocks. We recover zero
$f_{\mathrm{NL}}\ $as expected from the mocks and find no bias within our
pipeline. The red dashed line shows the best fit $\chi^{2}$ and the best fit
$f_{\mathrm{NL}}\ $value of the data. Figure 6 shows the constraints on
$f_{\mathrm{NL}}$ from the mean power spectrum of the 1000 null mocks with the
optimal weights (black) and without (red). We can recover zero non-Gaussianity
within 1-$\sigma$ and find no significant bias in the pipeline.

Figure 6: Marginalized posterior distribution of $f_{\mathrm{NL}}\ $from the
mean power-spectrum monopole of the null mocks but assuming the data
covariance. Results for the null mocks with (black) and without (red) optimal
weights.

###  7.2 Observational systematics

Using mocks with and without linear systematics added, as described in Section
5, we now test our ability to recover the correct value of $f_{\mathrm{NL}}\
$. The biggest concern is that, unfortunately any treatment of systematic to
remove contamination also has the potential to remove part of the signal and
decrease the power especially at large scales, where there are only very few
modes. Here we study this effect by applying the neural network based
correction to the contaminated mocks and compare the results to the null
mocks. The left-hand panel of Figure 7 shows the mean power spectrum monopole
of the contaminated mocks after neural net removal of systematics (black) and
the null mocks (green) for the NGC (upper panel) and SGC (lower panel). The
power spectrum at large scales has a reduced power after the neural network
treatment compared to the null, i.e. uncontaminated, mocks. This
consequentially shifts the $f_{\mathrm{NL}}\ $maximum likelihood to negative
values resulting in a mild bias of the $f_{\mathrm{NL}}\ $constraints. Note,
that we do not use the first $k$-bin in the MCMC analysis due to its
instability to the systematic correction. The maximum likelihood value is
$f_{\mathrm{NL}}=-12$ with optimal weights and $f_{\mathrm{NL}}=-19$. This
gives us an estimate of the overall bias caused by the systematic treatment.

Rezaie et al. (2021) present a method to correct for this over-subtraction,
calibrated on the mocks. This removes this small bias. However, the systematic
errors are only applied to the mocks using a linear model with two templates,
and the effect of the NN correction is significantly smaller than for the
data. Thus, we do not trust this correction when applied to the data. Given
the small amplitude of this effect, we choose to not apply a correction,
rather than apply an incorrect correction.

Figure 7: Impact of the treatment of systematics on the mocks: Left panel:
Mean power spectrum monopole from the null mocks (green) and contaminated
mocks (black) with the neural network removal technique applied to each
individual mock. Right panel: Corresponding $f_{\mathrm{NL}}\ $constraint from
the mean of the contaminated mocks (with errors as for a single mock).

##  8 Results

In this section we present the results of measuring the local, primordial non-
Gaussianity parameter $f_{\mathrm{NL}}\ $, summarised in Table 2. We discuss
the impact of optimal weights to improve the constraints, the impact of
photometric systematics and the merging history of halos. Next, we evaluate
the robustness of our results by varying our default analysis assumptions
including the theoretical model as well as scale cuts and priors (see Table
3).

###  8.1 Constraints on the local PNG

We constrain the local PNG parameter $f_{\mathrm{NL}}\ $to
$-33<f_{\mathrm{NL}}<10$ at 68% C.L and $-51<f_{\mathrm{NL}}<28$ at 95% C.L
using the eBOSS quasar sample with a redshift range of $0.8<z<2.2$ and a
minimum scale of $k_{\mathrm{min}}=0.019\,h/{\rm Mpc}$ included in the
analysis. These are the tightest constraints on PNG from a spectroscopic
survey to date covering the largest survey volume and using information on the
largest cosmological scales. We find no evidence of PNG with the
$f_{\mathrm{NL}}\ $parameter agreeing with zero within the 1-$\sigma$
uncertainty. We improve upon the eBOSS DR14 results of
$-81<f_{\mathrm{NL}}<26$ at 95% C.L by 35%. The default settings for our
analysis are summarised in Table 1. Figure 10 shows the power spectrum
monopole in the NGC (upper panel) and SGC (lower panel) with the best fit
model (blue line) and the 1-$\sigma$ uncertainty in $f_{\mathrm{NL}}\ $(grey
contours). Overall, the best fit of our model to the data is very good with
$\chi^{2}$/dof = 94/119.

Here we have applied a set of optimal weights including redshift weights and
optimal FKP weights as discussed in Section 6.3. The redshift weighting
technique includes information of the underlying theoretical model and
increases the effective redshift of the sample from $z_{\mathrm{eff}}=1.51$ to
$z_{\mathrm{eff}}=1.82$. Since the impact of PNG is stronger at high redshift
the constraints can be improved. Additionally the FKP weight is chosen at
$P_{0}=20000\ (\mathrm{Mpc}/h)^{3}$ instead of the standard $P_{0}=6000\
(\mathrm{Mpc}/h)^{3}$ that is used in BAO and RSD analysis to balance shot
noise and cosmic variance at large scales, which are of interest in this
study. Figure 9 shows the the constraints on $f_{\mathrm{NL}}\ $with (black
line) and without (red line) the optimal weights; our results using the NN
weights and merging factor $p=1.6$ are shown in the upper-left panel. The
optimal weights lead to an improvement of 37% compared to the constraint
without optimal weights of $-28<f_{\mathrm{NL}}<31$ at 68% C.L and a mild
shift to higher $f_{\mathrm{NL}}\ $values. The redshift weighting technique is
particularly effective for the quasar sample because of its large redshift
range but also because of the strong redshift dependency of the quasar bias.

Next, we study the effect of the photometric contaminations in the sample and
different systematics removal techniques. The photometric systematics in
particular effect the largest scales of the quasar power spectrum monopole
(see Figure 3) and their correct removal is crucial for any study of PNG. In
this study we focus on two systematic treatments:

  * •

Linear treatment: standard linear weights as used in the BAO analysis

  * •

Neural network treatment: advanced machine learning based approach

In Section 4 we summarise each technique for the data. The details of the
systematic treatment used for this analysis discussed in detail in the
accompanying paper of Rezaie et al. (2021). The main differences between the
two approaches is that the neural network technique allows us to remove non-
linear effects caused by the photometric contaminations while the standard
weights only remove linear trends. Additionally, Rezaie et al. (2021)
identified that the stellar density was better matched using a new template
constructed using Gaia data, whereas the standard treatment used a template
based on SDSS data (see Section 4.2). Figure 9 shows the constraints on
$f_{\mathrm{NL}}\ $for the linear systematics treatment. With the optimal
weights, the linear $f_{\mathrm{NL}}\ $measurement is in agreement with zero
at 1-$\sigma$ and is consistent with the neural network systematics removal.
However, without the optimal weights the maximum likelihood value of
$f_{\mathrm{NL}}\ $shifts by almost 4 $\sigma$, hinting at systematic
contamination at low redshifts. For the linear systematic removal, our
constraints on $f_{\mathrm{NL}}\ $are not robust and confirm the findings of
residual systematics in the standard catalogue as found in Rezaie et al.
(2021). This potential bias in $f_{\mathrm{NL}}\ $highlights the need for
robust systematics removal techniques to account for photometric contaminants.

Third, we discuss the impact of the merger history of quasars on the
$f_{\mathrm{NL}}\ $constraints parametrised as $p$ in Eq. 3. Tracers that
mainly populate recently merged halos lead to a weaker scale dependency of the
halo bias with $p\sim 1.6$. While this is our default assumption, it is likely
that not all quasars are in halos that resulted from recent mergers. Our
default results can therefore be considered the most conservative estimate. If
we instead assume $p=1$, i.e. no recent merger history, our constraints
tighten to $-18<f_{\mathrm{NL}}<11$ with the optimal weighing technique. This
is the most optimistic estimate of the $f_{\mathrm{NL}}\ $uncertainties. Note,
that the effective redshift in this case is $z_{\mathrm{eff}}=1.74$ and the
improvement of the optimal weights compared to the standard analysis without
optimal weights is 10%. The redshift evolution of the optimal weights with
$p=1$ is weaker that for $p=1.6$ resulting in less improvement and a lower
overall effective redshift. Figure 8 shows the redshift distribution, $n(z)$
for different weights.

Figure 8: Redshift distribution of the eBOSS QSO sample in the NGC without
optimal weights (blue), with optimal weights for $p=1$ (orange) and optimal
weights with $p=1.6$ (green). The optimal weights increases the effective
redshift from 1.51 to 1.74 ($p=1$) and 1.83 $p=1.6$, shifting the redshift
distribution to peak at higher redshifts. The parameter $p$ depends on the
assumed merger history of the QSO sample and defines the scale dependency of
the PNG contribution to the halo bias.

Nevertheless, the constraints for $p=1$ are 80% tighter than for $p=1.6$.
Therefore, our 1-$\sigma$ results are between $\pm
14<\sigma(f_{\mathrm{NL}})<\pm 21$. Note, that the maximum likelihood value
for $f_{\mathrm{NL}}\ $does not change with p, only the uncertainty on
$f_{\mathrm{NL}}\ $varies with p.

survey | eBOSS
---|---
tracer | quasars
redshift range | $0.8<z<2.2$
effective redshift | $z_{\mathrm{eff}}=1.82$
scale range in $h$/Mpc | $0.0019<k<0.121$
k bins | 63 with $\Delta k=2\pi/L\approx 0.019$ $h$/Mpc
photometric treatment | neural network weights
nuisance parameters | {$b_{N}$, $b_{S}$, $\Delta P^{\mathrm{shot}}_{N}$, $\Delta P^{\mathrm{shot}}_{S}$, $\sigma_{v,N}$, $\sigma_{v,S}$}
merger rate p | 1.6
Weights | redshift and optimised FKP weights
Table 1: Default assumptions of the analysis

Figure 9: Marginalized posterior distribution of $f_{\mathrm{NL}}\ $for the Neural Network systematic treatment [left] and linear systematic treatment [right]. Lower panels for p=1 assuming no recent merger history and upper panels for p=1.6. These are our main results from the DR16 QSO data. Figure 10: Power spectrum monopole for the NGC [upper panel] and SGC [lower panel]. The best fit model is plotted in blue. The blue contours show the range of models within the 1-$\sigma$ uncertainty in $f_{\mathrm{NL}}\ $. p | Systematic treatment |  $f_{\mathrm{NL}}\ $68% C.L. |  $f_{\mathrm{NL}}\ $95% C.L | Max. Likelihood | Best fit | Marg. Mean | $\chi^{2}/\mathrm{dof}$
---|---|---|---|---|---|---|---
1.6 | Neural Network | $\mathbf{-33<f_{\mathrm{NL}}<10}$ | $\mathbf{-51<f_{\mathrm{NL}}<28}$ | -12 | -13 | -12 | 94/119
Linear | $-50<f_{\mathrm{NL}}<10$ | $-68<f_{\mathrm{NL}}<36$ | -26 | -6 | -17 | 87/119
w/o Optimal Weights |  |  |  |  |  |
Neural Network | $-28<f_{\mathrm{NL}}<31$ | $-61<f_{\mathrm{NL}}<51$ | 4 | 11 | -1 | 99/119
| Linear | $47<f_{\mathrm{NL}}<103$ | $-114<f_{\mathrm{NL}}<133$ | 78 | 80 | 62 | 89/119
1 | Neural Network | $-17<f_{\mathrm{NL}}<11$ | $-31<f_{\mathrm{NL}}<22$ | -2 | 1 | -4 | 98/119
Linear | $-39<f_{\mathrm{NL}}<25$ | $-51<f_{\mathrm{NL}}<33$ | 13 | 13 | -7 | 87/119
w/o Optimal Weights |  |  |  |  |  |
Neural Network | $-14<f_{\mathrm{NL}}<18$ | $-32<f_{\mathrm{NL}}<29$ | 3 | 10 | 0 | 99/119
| Linear | $25<f_{\mathrm{NL}}<52$ | $-65<f_{\mathrm{NL}}<65$ | 40 | 44 | 31 | 89/119
Table 2: 68% and 95% C.L. constraints on local PNG for the systematic
treatment using neural networks as well as the standard removal technique
using linear weights. We present the results with the optimal weights and
without, as well as for different assumptions on the recent merger history
parametrised by p. These results are also displayed in Figure 9. The first
row, showing the Neural Network treatment with optimal weights is our default
result.

###  8.2 Robustness of the results

Figure 11: Robustness tests: 68% C.L constraints for different analysis settings and assumptions on the theoretical model. The Neural Network systematic treatment gives consistent results, while the standard linear systematic removal techniques shows scatter. The grey dashed lines correspond to zero $f_{\mathrm{NL}}^{\mathrm{local}}\ $.  | Neural Networks Systematic treatment | Linear Systematic treatment
---|---|---
Assumption |  $f_{\mathrm{NL}}\ $68% C.L. | Max. Like. | $\chi^{2}/\mathrm{dof}$ |  $f_{\mathrm{NL}}\ $68% C.L. | Max. Like. | $\chi^{2}/\mathrm{dof}$
Default | $-33<f_{\mathrm{NL}}<11$ | -12 | 94/119 | $-50<f_{\mathrm{NL}}<10$ | -28 | 87/119
NGC only | $-35<f_{\mathrm{NL}}<10$ | -13 | 44/59 | $-34<f_{\mathrm{NL}}<7$ | -19 | 42/59
SGC only | $-58<f_{\mathrm{NL}}<33$ | -4 | 50/59 | $-42<f_{\mathrm{NL}}<23$ | -5 | 45/59
Covariance w/o sys. error |  |  |  |  |  |
$-33<f_{\mathrm{NL}}<14$ | -9 | 105/119 | $-46<f_{\mathrm{NL}}<23$ | 6 | 99/54
minimum k [$h$/Mpc] |  |  |  |  |  |
$k_{\mathrm{min}}=0.0$ | $-26<f_{\mathrm{NL}}<8$ | -14 | 95/121 | $-38<f_{\mathrm{NL}}<15$ | -27 | 87/121
$k_{\mathrm{min}}=0.0038$ | $-71<f_{\mathrm{NL}}<27$ | -9 | 92/117 | $-117<f_{\mathrm{NL}}<-4$ | -3 | 86/117
$k_{\mathrm{min}}=0.0057$ | $-101<f_{\mathrm{NL}}<17$ | -34 | 91/115 | $-108<f_{\mathrm{NL}}<2$ | -48 | 86/115
$k_{\mathrm{min}}=0.0076$ | $-107<f_{\mathrm{NL}}<8$ | -48 | 91/113 | $-109<f_{\mathrm{NL}}<4$ | -52 | 86/113
maximum k [$h$/Mpc] |  |  |  |  |  |
$k_{\mathrm{max}}=0.1$ | $-33<f_{\mathrm{NL}}<15.$ | -8 | 74/99 | $-46<f_{\mathrm{NL}}<25$ | 1 | 65/99
shot noise |  |  |  |  |  |
$\Delta P_{\mathrm{shot}}>0$ | $-33<f_{\mathrm{NL}}<11$ | -11 | 95/119 | $-47<f_{\mathrm{NL}}<16$ | -20 | 88/119
$-1000<\Delta P_{\mathrm{shot}}>1000$ | $-32<f_{\mathrm{NL}}<11$ | -9 | 94/119 | $-44<f_{\mathrm{NL}}<17$ | -8 | 87/119
Gaussian $\sigma_{v}$ | $-33<f_{\mathrm{NL}}<11$ | -9 | 94/119 | $-49<f_{\mathrm{NL}}<-13$ | -27 | 87/119
Table 3: Robustness tests using different analysis settings. These results are
also displayed in Figure 11. The ’Default’ values (third row) are the main
result of the paper.

In this section we investigate the effect of different assumptions in the
analysis, for instance, the covariance matrix without systematic contribution,
the theoretical model as well as scale cuts and assumptions on priors for both
the neural network systematics treatment as well as the linear systematics
treatment. Table 3 and Figure 11 summarise our results. Overall, the results
using the neural network systematics treatment are robust to all assumptions
and analysis settings, whereas the linear systematic treatment shows a wider
scatter in the resulting $f_{\mathrm{NL}}\ $constraints. In the following we
discuss each scenario in more detail.

First, we investigate the effect of using the covariance calculated from the
‘null’ mocks, i.e. the purely statistical covariance matrix that does not
include any uncertainty due to the systematic treatment, in contrast to our
default covariance that includes scatter caused by the systematic weights (see
Section 6.5). In the neural network scenario, the constraints on
$f_{\mathrm{NL}}\ $are consistent with our default results with only a mild
increase in the $\chi^{2}/$dof, whereas in the linear scenario the maximum
likelihood value for $f_{\mathrm{NL}}\ $is shifted by more than 1 $\sigma$
showing that the linear results are not robust.

Next, we check the impact of using different scale cuts by changing the
minimum and maximum scales fitted, $k_{\mathrm{min}}$ and $k_{\mathrm{max}}$.
For the neural network technique, increasing $k_{\mathrm{min}}$ to remove
large scales, leads to an increase in the uncertainty in agreement with the
$1/k^{2}$ dependency of the PNG model and only mildly shifts the maximum
likelihood to lower values. The linear treatment on the other hand displays
strong changes to the $f_{\mathrm{NL}}\ $constraints, hinting at residue
systematic contamination. Removing the small scale information by decreasing
$k_{\mathrm{max}}$ has a minimal impact on the results for both scenarios,
only marginally increasing the uncertainty.

Additionally, we show that the prior on the shot noise has no influence on the
results for both scenarios. Lastly, we also note that the choice of a Gaussian
velocity dispersion in contrast to a Lorentzian behaviour does not change our
results.

###  8.3 Potential bias caused by the systematic treatment

Any treatment of systematics is liable to bias our results towards negative
$f_{\mathrm{NL}}\ $values by removing clustering modes and thus decreasing the
amplitude of the power spectrum monopole on large scales. In Section 4 we
estimated this effect using contaminated mocks to give $\Delta
f_{\mathrm{NL}}=-12$. Un-biasing our results by shifting the likelihood by
this value implies a $f_{\mathrm{NL}}\ $constraint of $f_{\mathrm{NL}}=0\pm
21\mathrm{(stat.)}\pm 12\mathrm{(sys.)}$. This over-correction is
significantly less than the statistical error. Since our estimate of the bias
relies on mocks with linear systeamtics that do not fully represent the non-
linearities of the systematics in the actual data, we do not include this bias
in our main results.

##  9 Conclusions

The quasar sample of eBOSS is the largest volume survey up to date, covering a
redshift range of $0.8<z<2.2$. The large survey volume allows the measurement
of the power spectrum at very large scales to an unprecedented precision. We
measure the amplitude of local PNG to be $f_{\mathrm{NL}}=-12\pm 21$ at 68
$\%$ confidence level finding no evidence for PNG within the 1-$\sigma$
uncertainty. This is the tightest constraint on $f_{\mathrm{NL}}\ $from low
redshift large-scale structure to date. Figure 12 shows this result in context
of previous analysis. In orange are the main measurements on $f_{\mathrm{NL}}\
$from CMB experiments over the years (Komatsu et al., 2003; Spergel et al.,
2007; Komatsu et al., 2009; Komatsu, 2010; Ade et al., 2014, 2016;
Collaboration et al., 2019), in blue constraints from LSS surveys (Slosar et
al., 2008; Ross et al., 2013b; Castorina et al., 2019) and in green results
from cross-correlations of LSS and CMB (Giannantonio et al., 2014a). We can
improve upon the previous eBOSS DR14 results (Castorina et al., 2019) by 35%,
by a factor of 3 upon the BOSS LRG DR9 results (Ross et al., 2013b) and by 28%
upon the constrains of Slosar et al. (2008), which used a combination of data
sets with the results dominated by the SDSS photometric quasar sample.

Figure 12: Summary of PNG measurements at 95 % C.L. from various CMB surveys
(Komatsu et al., 2003; Spergel et al., 2007; Komatsu et al., 2009; Komatsu,
2010; Ade et al., 2014, 2016; Collaboration et al., 2019) and LSS data (Slosar
et al., 2008; Ross et al., 2013b; Castorina et al., 2019) including this work
with $-51<f_{\mathrm{NL}}<28$. Note, that we used the median $f_{\mathrm{NL}}\
$value in case the best-fit or maximum likelihood value was not quoted in the
reference. The grey dashed line corresponds to zero
$f_{\mathrm{NL}}^{\mathrm{local}}\ $.

We apply optimised weights including redshift weights developed in Mueller et
al. (2019) as well as FKP weights more suitable for large scales to improve
our results by 37% highlighting the potential of advanced weighting
techniques.

We show that a robust removal of the photometric contamination is crucial for
a precise PNG measurement using the power spectrum monopole. Imaging
systematics such as stellar density, extinction or depth, can have a immense
effect in particular on the largest scales and therefore have the potential to
bias the $f_{\mathrm{NL}}\ $measurement.

In addition to the linear weights of the standard catalogue as were used in
the BOSS and eBOSS BAO and RSD analysis, we present results using a neural
network treatment of systematics as detailed in the companion paper Rezaie et
al. (2021). We find robust results using the advanced neural network
technique, on the other side, while the linear regression does not give robust
results. We conclude, that there are uncorrected contaminants present in the
data when using the linear regression technique. Therefore, for a PNG
measurement the standard treatment of systematics in not sufficient.

This paper is the latest step in the chain of development required to make
$f_{\mathrm{NL}}\ $measurements from galaxy redshift surveys. We have
demonstrated an increasingly sophisticated pipeline including careful choice
of weights and analysis of potential problems caused by the observations.
These developments are particularly important given future surveys including
DESI (DESI Collaboration et al., 2016a, b), SPHEREx (Doré et al., 2014) and
Euclid (Laureijs et al., 2011), which will cover larger effective volumes than
that considered here. Dedicated experiments such as a survey of high redshift
galaxies by an instrument such as the Maunakea Spectroscopic Explorer
(Percival et al., 2019) offer a potential to achieve errors of order $\Delta
f_{\mathrm{NL}}\sim 1$, which is the level at which interesting constraints
can be set on inflation (Bartolo et al., 2004). Robsutly realising this
potential will require further developments, particularly for the mitigation
of systematics arising from the observations. Including measurements of the
bispectrum and trispectrum have the potential to significantly improve on
these constraints (Karagiannis et al., 2018) further increasing the impact of
large-scale structure observations. What is clear is that over the next
decade, large-scale structure has the potential overtake the CMB as the
primary observation used to constrain $f_{\mathrm{NL}}\ $.

## Acknowledgements

E.-M. M. would like to thank Chris Blake for helpful discussions. E.-M. M.
acknowledges support from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 693024). M.R. and H.-J.S. are supported by the U.S. Department of
Energy, Office of Science, Office of High Energy Physics under DE-SC0014329
and DE-SC0019091. This work made substantial use of open-source software and
modules, such as Nbodykit, HEALPix, Fitsio, NumPy, SciPy, IPython, Jupyter,
and GitHub.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions.

SDSS-IV acknowledges support and resources from the Center for High
Performance Computing at the University of Utah. The SDSS website is
www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

## DATA AVAILABILITY

The new catalogue used in this article will be made public in Github at
https://github.com/mehdirezaie/eBOSSDR16QSOE. The details of the neural
network catalog creation are described in Rezaie et al. (2021). The standard
LSS catalog used in this work is available at
https://www.sdss.org/dr16/spectro/lss/.

## References

  * Ade et al. (2014) Ade P. A. R., et al., 2014, Astronomy & Astrophysics, 571, A24
  * Ade et al. (2016) Ade P. A. R., et al., 2016, Astronomy & Astrophysics, 594, A17
  * Ahumada et al. (2020) Ahumada R., et al., 2020, ApJS, 249, 3
  * Ata et al. (2018) Ata M., et al., 2018, Mon. Not. Roy. Astron. Soc., 473, 4773
  * Bardeen et al. (1986) Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, The Astrophysical Journal, 304, 15
  * Bartolo et al. (2004) Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Phys. Rep., 402, 103
  * Bautista et al. (2018) Bautista J. E., et al., 2018, The Astrophysical Journal, 863, 110
  * Bautista et al. (2020) Bautista J. E., et al., 2020, arXiv e-prints, p. arXiv:2007.08993
  * Beutler et al. (2017) Beutler F., et al., 2017, MNRAS, 466, 2242
  * Bianchi et al. (2015) Bianchi D., Gil-Marín H., Ruggeri R., Percival W. J., 2015, MNRAS, 453, L11
  * Blanton et al. (2017) Blanton M. R., et al., 2017, The Astronomical Journal, 154, 28
  * Castorina et al. (2019) Castorina E., et al., 2019, J. Cosmology Astropart. Phys., 2019, 010
  * Chuang et al. (2014) Chuang C.-H., Kitaura F.-S., Prada F., Zhao C., Yepes G., 2014, Monthly Notices of the Royal Astronomical Society, 446, 2621–2628
  * Collaboration et al. (2019) Collaboration P., et al., 2019, Planck 2018 results. IX. Constraints on primordial non-Gaussianity (arXiv:1905.05697)
  * Crocce et al. (2015) Crocce M., et al., 2015, Monthly Notices of the Royal Astronomical Society, 455, 4301
  * DESI Collaboration et al. (2016a) DESI Collaboration et al., 2016a, arXiv e-prints, p. arXiv:1611.00036
  * DESI Collaboration et al. (2016b) DESI Collaboration et al., 2016b, arXiv e-prints, p. arXiv:1611.00037
  * Dalal et al. (2008) Dalal N., Doré O., Huterer D., Shirokov A., 2008, Physical Review D, 77
  * Dawson et al. (2016) Dawson K. S., et al., 2016, The Astronomical Journal, 151, 44
  * Doré et al. (2014) Doré O., et al., 2014, arXiv e-prints, p. arXiv:1412.4872
  * Feldman et al. (1994) Feldman H. A., Kaiser N., Peacock J. A., 1994, The Astrophysical Journal, 426, 23
  * Ferramacho et al. (2014) Ferramacho L. D., Santos M. G., Jarvis M. J., Camera S., 2014, MNRAS, 442, 2511
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gangui et al. (1994) Gangui A., Lucchin F., Matarrese S., Mollerach S., 1994, ApJ, 430, 447
  * Giannantonio et al. (2014a) Giannantonio T., Ross A. J., Percival W. J., Crittenden R., Bacher D., Kilbinger M., Nichol R., Weller J., 2014a, Physical Review D, 89
  * Giannantonio et al. (2014b) Giannantonio T., Ross A. J., Percival W. J., Crittenden R., Bacher D., Kilbinger M., Nichol R., Weller J., 2014b, Phys. Rev. D, 89, 023511
  * Gil-Marín et al. (2020) Gil-Marín H., et al., 2020, arXiv e-prints, p. arXiv:2007.08994
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Gualdi et al. (2021) Gualdi D., Gil-Marin H., Verde L., 2021, Joint analysis of anisotropic power spectrum, bispectrum and trispectrum: application to N-body simulations (arXiv:2104.03976)
  * Gunn & Gott (1972) Gunn J. E., Gott J. Richard I., 1972, ApJ, 176, 1
  * Gunn et al. (2006) Gunn J. E., et al., 2006, Astron. J., 131, 2332
  * Hartlap et al. (2007) Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399
  * Ho et al. (2015) Ho S., et al., 2015, J. Cosmology Astropart. Phys., 2015, 040
  * Hou et al. (2021) Hou J., et al., 2021, MNRAS, 500, 1201
  * Kalus et al. (2019) Kalus B., Percival W. J., Bacon D. J., Mueller E. M., Samushia L., Verde L., Ross A. J., Bernal J. L., 2019, MNRAS, 482, 453
  * Karagiannis et al. (2018) Karagiannis D., Lazanu A., Liguori M., Raccanelli A., Bartolo N., Verde L., 2018, MNRAS, 478, 1341
  * Komatsu (2010) Komatsu E., 2010, Classical and Quantum Gravity, 27, 124010
  * Komatsu & Spergel (2001) Komatsu E., Spergel D. N., 2001, Phys. Rev. D, 63, 063002
  * Komatsu et al. (2003) Komatsu E., et al., 2003, Astrophys. J. Suppl., 148, 119
  * Komatsu et al. (2009) Komatsu E., et al., 2009, The Astrophysical Journal Supplement Series, 180, 330–376
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Laurent et al. (2017) Laurent P., et al., 2017, J. Cosmology Astropart. Phys., 2017, 017
  * Leistedt et al. (2014) Leistedt B., Peiris H. V., Roth N., 2014, Phys. Rev. Lett., 113, 221301
  * Lyke et al. (2020) Lyke B. W., et al., 2020, ApJS, 250, 8
  * Matarrese & Verde (2008) Matarrese S., Verde L., 2008, ApJ, 677, L77
  * Merz et al. (2021) Merz G., et al., 2021, arXiv e-prints, p. arXiv:2105.10463
  * Mueller et al. (2018) Mueller E.-M., Percival W., Linder E., Alam S., Zhao G.-B., Sánchez A. G., Beutler F., Brinkmann J., 2018, Monthly Notices of the Royal Astronomical Society, 475, 2122
  * Mueller et al. (2019) Mueller E.-M., Percival W. J., Ruggeri R., 2019, MNRAS, 485, 4160
  * Myers et al. (2015) Myers A. D., et al., 2015, Astrophys. J. Suppl., 221, 27
  * Neveux et al. (2020) Neveux R., et al., 2020, MNRAS, 499, 210
  * Percival et al. (2019) Percival W. J., et al., 2019, arXiv e-prints, p. arXiv:1903.03158
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A9
  * Raichoor et al. (2020) Raichoor A., et al., 2020, arXiv e-prints, p. arXiv:2007.09007
  * Rezaie et al. (2021) Rezaie M., et al., 2021, Monthly Notices of the Royal Astronomical Society
  * Ross et al. (2013b) Ross A. J., et al., 2013b, Monthly Notices of the Royal Astronomical Society, 428, 1116
  * Ross et al. (2013a) Ross A. J., et al., 2013a, MNRAS, 428, 1116
  * Ross et al. (2020a) Ross A. J., Bautista J., Tojeiro R., Brownstein J. R., Burtin E., Dawson K. S., et al. 2020a, MNRAS, 000, 000
  * Ross et al. (2020b) Ross A. J., et al., 2020b, MNRAS, 498, 2354
  * Ruggeri et al. (2017) Ruggeri R., Percival W. J., Gil-Marín H., Zhu F., Zhao G.-B., Wang Y., 2017, Monthly Notices of the Royal Astronomical Society, 464, 2698
  * Schlegel et al. (1998a) Schlegel D. J., Finkbeiner D. P., Davis M., 1998a, Astrophys. J., 500, 525
  * Schlegel et al. (1998b) Schlegel D. J., Finkbeiner D. P., Davis M., 1998b, The Astrophysical Journal, 500, 525
  * Scoccimarro (2015) Scoccimarro R., 2015, Phys. Rev. D, 92, 083532
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, MNRAS, 460, 3624
  * Slosar et al. (2008) Slosar A., Hirata C., Seljak U., Ho S., Padmanabhan N., 2008, J. Cosmology Astropart. Phys., 2008, 031
  * Smee et al. (2013) Smee S., et al., 2013, Astron. J., 146, 32
  * Smith et al. (2020) Smith A., et al., 2020, MNRAS, 499, 269
  * Spergel et al. (2007) Spergel D. N., et al., 2007, The Astrophysical Journal Supplement Series, 170, 377–408
  * Tamone et al. (2020) Tamone A., et al., 2020, arXiv e-prints, p. arXiv:2007.09009
  * Wang et al. (2019) Wang M. S., Percival W. J., Avila S., Crittenden R., Bianchi D., 2019, MNRAS, 486, 951
  * Wang et al. (2020) Wang M. S., Beutler F., Bacon D., 2020, MNRAS, 499, 2598
  * Wilson et al. (2017) Wilson M. J., Peacock J. A., Taylor A. N., de la Torre S., 2017, MNRAS, 464, 3121
  * Yamauchi et al. (2014) Yamauchi D., Takahashi K., Oguri M., 2014, Phys. Rev. D, 90, 083520
  * Zhao et al. (2016) Zhao G.-B., et al., 2016, Monthly Notices of the Royal Astronomical Society, 457, 2377
  * Zhao et al. (2019) Zhao G.-B., et al., 2019, MNRAS, 482, 3497
  * Zhao et al. (2021) Zhao C., et al., 2021, Monthly Notices of the Royal Astronomical Society, 503, 1149–1173
  * Zhu et al. (2015) Zhu F., Padmanabhan N., White M., 2015, Monthly Notices of the Royal Astronomical Society, 451, 236
  * Zhu et al. (2016) Zhu F., Padmanabhan N., White M., Ross A. J., Zhao G., 2016, Monthly Notices of the Royal Astronomical Society, 461, 2867
  * de Mattia & Ruhlmann-Kleider (2019) de Mattia A., Ruhlmann-Kleider V., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 036–036
  * de Mattia et al. (2020) de Mattia A., et al., 2020, arXiv e-prints, p. arXiv:2007.09008
  * du Mas des Bourboux et al. (2020) du Mas des Bourboux H., et al., 2020, ApJ, 901, 153
  * eBOSS Collaboration et al. (2020) eBOSS Collaboration et al., 2020, arXiv e-prints, p. arXiv:2007.08991