R2M

    raw 2nd moment
CSI

    channel state information
UE

    user equipment
EM

    electromagnetic
UL

    uplink
BS

    base station
TDD

    time division duplex
FDD

    frequency division duplex
ECC

    error-correcting code
MLD

    maximum likelihood decoding
HDD

    hard decision decoding
IF

    intermediate frequency
RF

    radio frequency
SDD

    soft decision decoding
NND

    neural network decoding
CNN

    convolutional neural network
ML

    maximum likelihood
GPU

    graphical processing unit
BP

    belief propagation
LTE

    Long Term Evolution
BER

    bit error rate
SNR

    signal-to-noise-ratio
ReLU

    rectified linear unit
BPSK

    binary phase shift keying
QPSK

    quadrature phase shift keying
AWGN

    additive white Gaussian noise
MSE

    mean squared error
LLR

    log-likelihood ratio
MAP

    maximum a posteriori
NVE

    normalized validation error
BCE

    binary cross-entropy
CE

    cross-entropy
BLER

    block error rate
SQR

    signal-to-quantisation-noise-ratio
MIMO

    multiple-input multiple-output
mMIMO

    massive multiple-input multiple-output
OFDM

    orthogonal frequency division multiplex
RF

    radio frequency
LOS

    line of sight
NLoS

    non-line of sight
NMSE

    normalized mean squared error
CFO

    carrier frequency offset
SFO

    sampling frequency offset
IPS

    indoor positioning system
TRIPS

    time-reversal IPS
RSSI

    received signal strength indicator
MIMO

    multiple-input multiple-output
ENoB

    effective number of bits
AGC

    automatic gain control
ADC

    analog to digital converter
ADCs

    analog to digital converters
FB

    front bandpass
FPGA

    field programmable gate array
JSDM

    Joint Spatial Division and Multiplexing
NN

    neural network
IF

    intermediate frequency
LoS

    line-of-sight
NLoS

    non-line-of-sight
DSP

    digital signal processing
AFE

    analog front end
SQNR

    signal-to-quantisation-noise-ratio
SINR

    signal-to-interference-noise-ratio
ENoB

    effective number of bits
PCB

    printed circuit board
EVM

    error vector mangnitude
CDF

    cumulative distribution function
MRC

    maximum ratio combining
MRP

    maximum ratio precoding
MRT

    maximum ratio transmission
DeepL

    deep-learning
DL

    downlink
SISO

    single-input single-output
SGD

    stochastic gradient descent
CP

    cyclic prefix
MISO

    Multiple Input Single Output
LMMSE

    linear minimum mean square error
ZF

    zero forcing
USRP

    universal software radio peripheral
RNN

    recurrent neural network
GRU

    gated recurrent unit
LSTM

    long short-term memory
NTM

    neural turing machine
DNC

    differentiable neural computer
TCN

    temporal convolutional network
FCL

    fully connected layer
MANN

    memory augmented neural network
RNN

    recurrent neural network
DNN

    deep neural network
FIR

    finite impulse response
BPTT

    back-propagation through time
GAN

    generative adversarial network
ELU

    exponential linear unit
tanh

    hyperbolic tangent
BICM

    bit-interleaved coded modulation
OTA

    over-the-air
IM

    intensity modulation
DD

    direct detection
RL

    reinforcement learning
SDR

    software-defined radio
WGAN

    Wasserstein GAN
BMD

    bit-metric decoding
BMI

    bit-wise mutual information
LDPC

    low-density parity-check
IDD

    iterative demapping and decoding
JSD

    Jensen-Shannon distance
MMSE

    minimum mean square error
FFT

    fast Fourier transform
DFT

    discrete Fourier transform
IFFT

    inverse fast Fourier transform
QAM

    quadrature amplitude modulation
EMD

    earth mover’s distance
TDL

    tapped delay line
KL

    Kullback-Leibler
PRACH

    physical random access channel
URLLC

    ultra-reliable low-latency communication
ANOMA

    asynchronous non-orthogonal multiple access
FEC

    forward error correction
PAPR

    peak-to-average power ratio
APP

    a posteriori probability
COTS

    commercial off-the-shelf
PLL

    phase locked loop
STO

    sampling time offset
SFO

    sampling frequency offset
CFO

    carrier frequency offset
CPO

    carrier phase offset
CSI

    channel state information
GNSS

    global navigation satellite system
ELAA

    extremely large aperture array
UE

    user equipment
DICHASUS

     Distributed Channel Sounder by University of Stuttgart
JCaS

    Joint Communication and Sensing
CT

    continuity
TW

    trustworthiness
KS

    Kruskal’s stress
PCA

    principal component analysis
AoA

    angle of arrival
RD

    Rajski’s distance
ADP

    angle-delay profile
MPC

    multipath component
CIR

    channel impulse response
MAE

    mean absolute error
MDS

    multidimensional scaling
t-SNE

    t-distributed stochastic neighbor embedding
SM

    Sammon’s mapping
CIRA

    channel impulse response amplitude
CS

    cosine similarity
FCF

    forward charting function
CMD

    correlation matrix distance
NR

    New Radio
3GPP

    3rd Generation Partnership Project
QuaDRiGa

    Quasi Deterministic Radio Channel Generator
TR

    technical report
VAE

    variational autoencoder
GSCM

    geometry-based stochastic channel model
GP

    gradient penalty
LTI

    linear time-invariant
RMS

    root mean square
UPA

    uniform planar array
MUSIC

    MUltiple SIgnal Classification
PDP

    power delay profile
PDF

    probability density function

# GAN-based Massive MIMO Channel Model Trained on Measured Data ††thanks: This
work is supported by the German Federal Ministry of Education and Research
(BMBF) within the projects Open6GHub (grant no. 16KISK019) and KOMSENS-6G
(grant no. 16KISK113).

Florian Euchner, Janina Sanzi, Marcus Henninger, Stephan ten Brink
Institute of Telecommunications, Pfaffenwaldring 47, University of Stuttgart,
70569 Stuttgart, Germany
<EMAIL_ADDRESS>
###### Abstract

Wireless channel models are a commonly used tool for the development of
wireless telecommunication systems and standards. The currently prevailing
geometry-based stochastic channel models (GSCMs) were manually specified for
certain environments in a manual process requiring extensive domain knowledge,
on the basis of channel measurement campaigns. By taking into account the
stochastic distribution of certain channel properties like Rician k-factor,
path loss or delay spread, they model the distribution of channel
realizations. Instead of this manual process, a generative machine learning
model like a generative adversarial network (GAN) may be used to automatically
learn the distribution of channel statistics. Subsequently, the GAN’s
generator may be viewed as a channel model that can replace conventional
stochastic or raytracer-based models. We propose a GAN architecture for a
massive MIMO channel model, and train it on measurement data produced by a
distributed massive MIMO channel sounder.

###### Index Terms:

channel model, generative adversarial network, massive MIMO

##  I Introduction

Channel modeling is essential in communications engineering. Fundamental
theoretic results are often derived starting from some of the simplest channel
models (e.g., memoryless discrete-time and discrete-value channels or additive
white Gaussian noise (AWGN) channels) and proposals for practical waveforms
and modulation formats are motivated using simple stochastic channel models
(e.g., Rician Fading, or Clarke’s model [1]). More involved channel models
have been developed for researching practical wireless communication systems,
that must be able to withstand various degrees of frequency- and/or time-
selective fading, hardware impairments and other effects. One such channel
model was proposed by the 3rd Generation Partnership Project (3GPP) for the
frequency range $0.5-100\,\mathrm{GHz}$ in technical report (TR) 38.901 [2],
which is an example for a so-called geometry-based stochastic channel model
(GSCM). These channel models define a set of stochastic parameters, including
information about the distribution of delay spreads or path loss exponents,
and guarantee some level of spatial coherence, for example by ensuring the
spatial distribution follows some three-dimensional, correlated noise process.
While stochastic channel models usually only define the model (e.g.,
procedures, distributions, parameters), channel simulators like Quasi
Deterministic Radio Channel Generator (QuaDRiGa) [3] may be used to produce
channel realizations, for example in the form of channel impulse responses.
For even better spatial consistency, raytracing-based channel simulators may
be used as channel models [4]. As evidenced by lots of successful, recent
research based on stochastic or raytracer-based channel models, these existing
solutions work well for the intended purposes. However, in both cases, there
is a lot of manual work required to come up with a model: To obtain a
stochastic channel model that describes a certain environment for a set of
system parameters (e.g., carrier frequency), one first has to obtain channel
measurements for the particular environment, to then identify relevant
stochastic parameters and to finally extract these parameters from the
measurements. The inherent danger in this process is that relevant parameters
may be overlooked or that the modeled environment simply does not match the
system the model is used for. Aforementioned TR 38.901 channel model is only
valid for a limited set of manually selected scenarios of interest, which
define environments labeled, for example, “urban macrocell” or “indoor
factory”.

Raytracing-based channel models offer more flexibility to the user, but
require significantly more effort related to the definition of the 3D model
and material properties. What is more, the user of a channel model may not
care about specifics of one particular environment, which a raytracing-based
channel model offers, but may want to rather evaluate some system on a certain
type of environment. While raytracing may be used to model an environment in a
faithful, consistent way, stochastic channel models provide more variety and
randomness in their channel realizations. Finding a trade-off between
faithfulness and randomness is a rather fundamental issue in channel modeling,
and the right balance has to be determined depending on the particular
application.

ChannelDataset $\mathcal{S}_{\mathrm{train}}$Generator $G$NoiseVector
$\mathbf{n}$DesiredPosition $\mathbf{x}$Critic
$C$$\mathbf{H}^{(l)},\mathbf{x}^{(l)}$$\mathbf{\hat{H}},\mathbf{x}$Score

Figure 1: Structural overview of the Wasserstein GAN for channel modeling. The
generator $G$ learns the CSI distribution of the channel.

This work111Partial source code is publicly available at:
https://github.com/Jeija/GAN-Wireless-Channel-Model analyzes a generative
machine learning-based approach, which can automate the procedure of
parametrizing a channel model on measurement data. Instead of specifying
relevant stochastic parameters manually, the generative channel model should
automatically learn the distribution of the measured CSI, based on samples
from one or multiple datasets. In particular, we explore the use of a GAN [5]
(a Wasserstein GAN (WGAN) [6] with gradient penalty (GP) [7], to be precise)
for automated learning of stochastic channel models. The basic structure of
such a GAN-based channel model is illustrated in Fig. 1.

###  I-A Related Work and Contributions

LinearChannelChannel Input$\mathbf{s}$Channel Output$\mathbf{y}$Generator
$G$CSI $\mathbf{H}$

(a) Linear type: Generator produces instantaneous CSI

Generator $G$Channel Input$\mathbf{s}$Channel Output$\mathbf{y}$

(b) Nonlinear type: Generator conditioned on channel input

Figure 2: Different types of GAN-based channel models.

(a)

$-10$$-6$$-2$$2$$-18$$-16$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (b)

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Training setcutoutCoordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (c)

$-10$$-5$$0$Total RX Power $\lVert\mathbf{H}^{(l)}\rVert_{\mathrm{F}}^{2}$ in
dB

Figure 3: Information about the environment the dataset was measured in: The
figure shows (a) a photograph of the environment, (b) a top view map and (c) a
scatter plot with datapoint positions $\mathbf{x}^{(l)}$ in
$\mathcal{S}_{\mathrm{meas}}$, colorized with the total received power
$\lVert\mathbf{H}^{(l)}\rVert_{\mathrm{F}}^{2}$. The antenna arrays in the map
are drawn to scale as black rectangles and their directivity is indicated by
the green sectors. The received power is normalized such that the maximum
received power is $0\,\mathrm{dB}$.

The idea of a GAN-based channel model that learns the CSI distribution, first
proposed in [8], has gained some attention in the community [9, 10, 11, 12,
13]. In [8], the concept of GAN-based channel modeling is introduced on a high
level, motivating the idea by showing that a GAN can learn the noise
distribution of an AWGN channel. In [9, 10, 11], a GAN is trained on channel
realizations (impulse responses). As an alternative to GANs, [12] investigates
the use of diffusion models for channel modeling and [13] compares a GAN-based
channel model to a variational autoencoder (VAE)-based model and a diffusion
model. In all previous research that the authors are familiar with, the
generative model is trained on data that was itself produced by another
channel model, usually related to aforementioned 3GPP TR 38.901 GSCM.

In a separate line of work, initiated by [14], various authors have
investigated using a conditional GAN to learn the input-output relationship of
the channel itself, which is a different type of GAN-based channel model. The
distinction between [8] and [14], which is also illustrated in Fig. 2, is the
following: In [8], the GAN’s generator produces channel realizations described
by CSI in some representation (e.g., time-domain or frequency-domain), as in
Fig. 2(a). On the other hand, in [14], the generator _becomes_ the channel
itself: The generator is conditioned on the transmitted signal, and it must
attempt to produce a realistic received signal, as in Fig. 2(b), which can be
useful e.g. for end-to-end training [15]. We restrict ourselves to the first,
linear type of GAN-based channel model as described in [8]. As already
outlined in [10], the advantages of generating CSI as an intermediary product
are the improved interpretability of the GAN’s output and the better
generalization to arbitrary input signals. Furthermore, it ensures that the
learned channel is indeed a linear time-invariant (LTI) system, whereas the
approach in [14] can learn nonlinear channel models, which is undesirable for
modeling physical propagation effects.

The architecture of our neural network draws inspiration from [9], being based
on a WGAN with GP. In contrast to all this previous work, however, we train
our GAN-based channel model on real-world data measured with a massive
multiple-input multiple-output (MIMO) channel sounder, as opposed to the
previous efforts which were trained on another channel model. In addition, we
investigate the use of a conditional GANs to generate channel realizations at
defined locations, while ensuring some level of spatial consistency.

###  I-B Purpose of a GAN-based channel model

The objective of a channel model depends on the application: A GSCM may
generate CSI that follows a certain distribution, a raytracer may be used to
predict realistic CSI. The purpose of GAN-based channel models has not been
well-defined in existing literature. Therefore, we consider two different
applications: Using GAN-based channel models as an alternative to GSCMs as a
kind of _random generator_ for CSI samples, or using GAN-based channel models
as a _spatial inter- and extrapolator_ , that should ideally predict measured
CSI at new locations. While the use of a GAN as an inter-/extrapolator has
been suggested to us many times in conversations, we would like to remark that
a different neural network structure, like an autoencoder, may be better-
suited in that case. In the following sections, we will analyze both types of
applications and conclude that, in their current form, GAN-based channel
models of linear type are not particularly well-suited for either.

##  II System Model and Measurement Dataset

Our objective is to model the distribution of a channel measurement dataset
containing $L$ entries, that we refer to as _datapoints_ $l=1,\ldots,L$. We
model the propagation channel between a distributed massive MIMO antenna
system, consisting of $B$ uniform planar arrays with half-wavelength spacing,
$M_{r}$ rows and $M_{c}$ columns of antennas each, and a single-antenna user
equipment (UE). The channel is described by its complex-valued, time-domain
CIR. The time-domain CSI tensors can thus be written as
$\mathbf{H}^{(l)}\in\mathbb{C}^{B\times M_{r}\times M_{c}\times
N_{\mathrm{tap}}}$, where $N_{\mathrm{tap}}$ is the number of tapped delay
line (TDL) taps. Furthermore, we assume that a UE position
$\mathbf{x}^{(l)}\in\mathbb{R}^{2}$ is known for every CSI measurement. Hence,
one datapoint is a tuple constisting of CSI tensor and UE position, and the
whole dataset may be written as the set

$\mathcal{S}_{\mathrm{meas}}=\left\\{\left(\mathbf{H}^{(l)},\mathbf{x}^{(l)}\right)\right\\}_{l=1,\ldots,L}.$

For this work, the dataset $\mathcal{S}_{\mathrm{meas}}$ was measured by _
Distributed Channel Sounder by University of Stuttgart (DICHASUS)_, our
distributed massive multiple-input multiple-output (mMIMO) channel sounder,
whose architecture is thoroughly described in [16]. In brief, DICHASUS is a
massive MIMO\- orthogonal frequency division multiplex (OFDM) channel sounder.
It achieves long-term phase-coherence, even if antennas are distributed.
DICHASUS provides large datasets containing OFDM channel coefficients for all
antennas and all subcarriers, along with side information like timestamps and
accurate information about the positions of all antennas at all times. While
measurements are performed in frequency domain, it is only a matter of
computing the Fourier transform to obtain the time-domain representation
$\mathbf{H}^{(l)}$.

The considered dataset $\mathcal{S}_{\mathrm{meas}}$, containing $L=83893$
datapoints, is a subset of _dichasus-cf0x_ [17], and was captured in an
industrial environment with $B=4$ separate antenna arrays made up of
$M=2\times 4=8$ antennas each. A total of $N_{\mathrm{sub}}=1024$ OFDM channel
coefficients were measured at a carrier frequency of $1.272\,\mathrm{GHz}$ and
with a channel bandwidth of $50\,\mathrm{MHz}$, though we restrict ourselves
to $N_{\mathrm{tap}}=48$ time-domain taps for all subsequent analyses. The
single dipole transmit antenna is mounted on top of a robot, which travels
along a set of trajectories inside a defined, L-shaped area, with an overall
size of approximately $14\,\mathrm{m}\times 14\,\mathrm{m}$. A photo and a top
view map of the environment are shown in Fig. 3. A large metal container is
located at the inner corner of the L-shape, blocking the line-of-sight (LoS).
The datapoint positions $\mathbf{x}^{(l)}$ are shown in Fig. 3.

We define training set
$\mathcal{S}_{\mathrm{train}}\subset\mathcal{S}_{\mathrm{meas}}$ and test set
$\mathcal{S}_{\mathrm{test}}\subset\mathcal{S}_{\mathrm{meas}}$ as two
disjoint subsets of the complete measurement dataset
$\mathcal{S}_{\mathrm{meas}}$. The training set contains
$|\mathcal{S}_{\mathrm{train}}|=17857$ datapoints, the test set contains
$|\mathcal{S}_{\mathrm{test}}|=20973$ datapoints. The training set is used to
train the GAN or to parametrize the baselines. The GAN (or a baseline) then
generates CSI at the test set locations, and the generated CSI is compared to
the true, measured CSI. The test set is obtained from
$\mathcal{S}_{\mathrm{meas}}$ by choosing every fourth datapoint. We do the
same for the training set, except with an offset, and then additionally cut a
”hole” with a diameter of $4\,\mathrm{m}$ into the training set (as visible in
Fig. 3(c) and Fig. 6(a)). This hole is meant to test the ability of the GAN
and baselines to interpolate CSI at previously unseen locations.

##  III GAN-Based Channel Model

Fig. 1 shows the well-known conditional WGAN training setup: For the desired
position $\mathbf{x}$ (the condition) and as a function of some noise vector
$\mathbf{n}=(n_{1},\ldots,n_{K})$, $n_{k}\sim\mathcal{N}(0,1)$, $K=128$, the
generator $G(\mathbf{x},\mathbf{n})$ produces a _generated_ channel
realization $\mathbf{\hat{H}}\in\mathbb{C}^{B\times M_{r}\times M_{c}\times
N_{\mathrm{tap}}}$. As its input, the critic $C$ either receives a _real_
(i.e., measured) channel $\mathbf{H}^{(l)}$ from
$\mathcal{S}_{\mathrm{train}}$, or a generated channel $\mathbf{\hat{H}}$,
always in conjunction with the condition. The critic’s objective is to detect
whether the input it receives is real or fake, by assigning a score reflective
of the _realness_ of the input sample. We train the generator and critic to
operate on the _complete_ CSI tensors $\mathbf{H}^{(l)}$, which contain time-
domain tap coefficients for all $B\times M_{r}\times M_{c}$ anntennas and all
taps.

Generator and critic are dense neural networks, structured as shown in Fig.
4(a) and Fig. 4(b). In the neural network, the complex-valued CSI is processed
in separate real / imaginary parts, and conversion and reshaping layers are
omitted from Fig. 4(a) and Fig. 4(b) for the sake of brevity. In addition to
the raw generated CSI $\hat{\mathbf{H}}$, the critic is also provided antenna-
specific root mean square (RMS) delay spread values derived from
$\hat{\mathbf{H}}$, the computation of which is explained in Sec. V-A2. The
behind providing RMS delay spread values is that they may aid the critic in
scoring the realness of a generated sample. The WGAN is trained with GP [7]
with GP coefficient $\lambda=10$. Delay spread values and the conditional
positions $\mathbf{x}$ are normalized such that all numbers are in range
$[-1,1]$. The CSI, except for being split into real and imaginary parts, does
not undergo any additional treatment like normalization.

Dense, 512 Neurons, ReLUDense, 512 Neurons, ReLUDense, 1024 Neurons,
ReLUDense, 2048 Neurons, ReLUDense, 3072 Neurons,
linear$\mathbf{n}\in\mathbb{R}^{K}$Position $\mathbf{x}\in\mathbb{R}^{2}$CSI
$\mathbf{\hat{H}}\in\mathbb{C}^{B\times M_{r}\times M_{c}\times
N_{\mathrm{tap}}}$

(a) Structure of generator neural network $G$, with random noise input
$\mathbf{n}$ and condition $\mathbf{x}$ (position)

Dense, 160 Neurons, ReLUDense, 100 Neurons, ReLUDense, 50 Neurons, ReLUDense,
20 Neurons, ReLUDense, 10 Neurons, ReLUDense, 1 Neuron, linearCSI
$\mathbf{\hat{H}}\in\mathbb{C}^{B\times M_{r}\times M_{c}\times
N_{\mathrm{tap}}}$$\hat{\mathbf{DS}}_{\mathrm{rms}}\in\mathbb{C}^{B\times
M_{r}\times M_{c}}$Position $\mathbf{x}\in\mathbb{R}^{2}$Critic score
$C(\mathbf{\hat{H}})\in\mathbb{R}$

(b) Structure of critic neural network $C$, where input
$\hat{\mathbf{DS}}\in\mathbb{C}^{B\times M_{r}\times M_{c}}$ denotes the
normalized delay spread computed from $\hat{\mathbf{H}}$

Figure 4: Structure of neural networks. Reshaping layers and the conversion
from complex representation to separate real / imaginary parts are omitted.

##  IV Linear Interpolation Baseline

To judge whether a GAN-based channel model is useful as an _interpolator_ , we
compare it to a simple classical interpolation baseline. While more advanced
interpolation methods like, for example, Kriging are well-known, we employ a
simpler linear interpolation baseline: We construct the interpolant
$F_{\mathrm{interp}}:\mathbb{R}^{2}\to\mathbb{C}^{B\times M_{r}\times
M_{c}\times N_{\mathrm{tap}}}$ by first computing the Delaunay triangulation
of the two-dimensional training set positions, and then performing linear
barycentric CSI interpolation on each triangle as follows: For a single
triangle, let $\mathbf{s}\in\mathbb{R}^{3}$ denote the barycentric coordinate
vector such that $\sum_{i}\mathbf{s}_{i}=1$ and let $\mathbf{H}^{(\Delta 1)}$,
$\mathbf{H}^{(\Delta 2)}$, $\mathbf{H}^{(\Delta 3)}$ denote the CSI tensors of
the corresponding vertices of the triangle. To account for unknown global
starting phases $\bm{\varphi}\in\mathbb{R}^{3}$ in the CSI tensors, we find
the interpolated CSI tensor $\hat{\mathbf{H}}$ as the solution to the
optimization problem

$(\hat{\mathbf{H}},\bm{\varphi})=\operatorname*{arg\,min}_{(\mathbf{H},\bm{\varphi})}\sum_{i=1}^{3}s_{i}\left\lVert\mathbf{H}^{(\Delta
i)}-\mathrm{e}^{\mathrm{j}\bm{\varphi}_{i}}\mathbf{H}\right\rVert_{\mathrm{2}}^{2},$

where $\lVert\cdot\rVert_{2}$ is the straightforward generalization of the
Euclidean norm to tensors and subtraction is applied elementwise. The
optimization problem is solved by separate differentiation with respect to
$\bm{\varphi}$ and $\mathbf{H}$, which yields a simple coordinate descent
interpolation algorithm.

##  V Evaluation

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{test}}$, $b=1$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{GAN,fix}}$, $b=1$, fixed
$\mathbf{n}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{test}}$, $b=2$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{GAN,fix}}$, $b=2$, fixed
$\mathbf{n}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{test}}$, $b=3$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{GAN,fix}}$, $b=3$, fixed
$\mathbf{n}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{test}}$, $b=4$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate
$\mathbf{x}_{2}~{}[\mathrm{m}]$$\mathcal{S}_{\mathrm{GAN,fix}}$, $b=4$, fixed
$\mathbf{n}$

$-25$$-20$$-15$$-10$$-5$$0$Total RX Power
$\left\lVert\mathbf{H}_{b}^{(l)}\right\rVert_{\mathrm{F}}^{2}$ in dB

Figure 5: Spatial distribution of received powers
$\lVert\mathbf{H}^{(l)}_{b}\rVert_{\mathrm{F}}^{2}$, for measured dataset
$\mathcal{S}_{\mathrm{meas}}$ and generated dataset
$\mathcal{S}_{\mathrm{GAN,fixed}}$. The power is normalized such that the
maximum received power (over _all_ values
$\lVert\mathbf{H}^{(l)}_{b}\rVert_{\mathrm{F}}^{2}$, measured and generated)
is $0\,\mathrm{dB}$.

Once GAN training has finished, the critic is no longer needed, and the
generator can be evaluated. Depending on whether one views the generative
channel model as an _interpolator_ or a _random generator_ , one may want to
consider different evaluation criteria. In this work, we compare both the
spatial distributions of generated and measured CSI as well as stochastic CSI
distributions. In contrast to generative models in image processing, where
subjective human perception of the quality of a generated image can play in
important role in the evaluation process, humans cannot easily interpret CSI
data. Instead, we need to rely on human-interpretable quantities that can be
derived from channel tensors $\mathbf{\hat{H}}$ to judge the generated data.
For example, we can interpret and compare (to measured data) the distribution
of powers, delay spreads or angular-domain properties in the generated data.

###  V-A Evaluated Channel Parameters and Statistics

We focus our evaluation on three properties of the generated CSI, those being
the received power, the delay spread and the angle of arrival (AoA).
Furthermore, in an exemplary manner, we compare the delay spread distributions
using the Jensen-Shannon distance (JSD).

####  V-A1 Received Power

The total received power over all antennas in one antenna array $b$ will be
written as the norm

$\left\lVert\mathbf{H}_{b}\right\rVert_{\mathrm{F}}^{2}=\sum_{m_{r}=1}^{M_{r}}\sum_{m_{c}=1}^{M_{c}}\sum_{t=1}^{N_{\mathrm{tap}}}\left|\mathbf{H}_{b,m_{r},m_{c},t}\right|^{2}.$

####  V-A2 Delay Spread

For a single antenna in row $m_{r}$ and column $m_{c}$ of array $b$, the RMS
delay spread is defined as

$\begin{split}\mathrm{DS}_{\mathrm{rms},b,m_{r},m_{c}}&=\sqrt{\frac{\sum_{t}(t-\bar{t}_{b,m_{r},m_{c}})^{2}\left|\mathbf{H}_{b,m_{r},m_{c},t}\right|^{2}}{\sum_{t}\left|\mathbf{H}_{b,m_{r},m_{c},t}\right|^{2}}}\\\
\text{with}~{}~{}\bar{t}_{b,m_{r},m_{c}}&=\frac{\sum_{t}t\left|\mathbf{H}_{b,m_{r},m_{c},t}\right|^{2}}{\sum_{t}\left|\mathbf{H}_{b,m_{r},m_{c},t}\right|^{2}},\end{split}$

where $t=1,\ldots,N_{\mathrm{tap}}$ is an index for the time tap. To interpret
the results, we may also want to visualize the RMS delay spread averaged over
all antennas in one array

$\mathrm{\overline{DS}}_{\mathrm{rms},b}=\frac{1}{M_{r}M_{c}}\sum_{m_{r}=1}^{M_{r}}\sum_{m_{c}=1}^{M_{c}}\mathrm{DS}_{\mathrm{rms},b,m_{r},m_{c}}.$

By comparing the measured and generated RMS delay spread distribution, we can
test if the GAN learns meaningful channel power delay profiles.

####  V-A3 Angle of Arrival

The array correlation matrix
$\mathbf{R}^{(b)}\in\mathbb{C}^{M_{\mathrm{c}}\times M_{\mathrm{c}}}$ for the
array with index $b$ is estimated over all rows of the UPA. That is,
$\mathbf{R}^{(b)}$ may be used to estimate the azimuth AoA, but contains no
information about the elevation AoA:

$\mathbf{\hat{R}}^{(b)}_{m_{c1},m_{c2}}=\sum_{m_{\mathrm{r}}=1}^{M_{\mathrm{r}}}\sum_{t=1}^{N_{\mathrm{tap}}}\mathbf{H}_{b,m_{\mathrm{r}},m_{\mathrm{c1}},t}\mathbf{H}_{b,m_{\mathrm{r}},m_{\mathrm{c2}},t}^{*}\in\mathbb{C}$

Under the assumption of a single source, we apply the root- MUltiple SIgnal
Classification (MUSIC) algorithm to $\mathbf{\hat{A}}^{(b)}$ to find azimuth
AoA estimates $\hat{\alpha}^{(b)}$ from $\mathbf{\hat{A}}^{(b)}$ for each
datapoint and each UPA. An azimuth AoA of $\hat{\alpha}^{(b)}=0^{\circ}$
indicates that the signal arrives directly from the front. By comparing
measured and generated AoA distributions, we can test if the GAN learns a
meaningful distribution of channel coefficient phases.

####  V-A4 Jensen-Shannon Distance

Furthermore, to quantify the similarity of different stochastic distributions,
we use the JSD. For two distributions $P$ and $Q$ over the same sample space
$\mathcal{X}$, the JSD is defined as

$\mathrm{d}_{\mathrm{JS}}(P||Q)=\sqrt{\frac{d_{\mathrm{KL}}(P||M)+d_{\mathrm{KL}}(Q||M)}{2}}~{}~{}\text{with}~{}~{}M=\frac{P+Q}{2},$

where $\mathrm{d}_{\mathrm{KL}}$ is the Kullback-Leibler (KL) divergence,
given by

$\mathrm{d}_{\mathrm{KL}}(P||Q)=\sum_{x\in\mathcal{X}}P(x)\log\left(\frac{P(x)}{Q(x)}\right).$

###  V-B GAN Evaluation Procedure

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (a)
Training Set $\mathcal{S}_{\mathrm{train}}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (b)
Test Set $\mathcal{S}_{\mathrm{test}}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (c)
$\mathcal{S}_{\mathrm{GAN,fix}}$, fixed $\mathbf{n}^{(1)}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (d)
$\mathcal{S}_{\mathrm{GAN,fix}}$, fixed $\mathbf{n}^{(2)}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (e)
$\mathcal{S}_{\mathrm{GAN,var}}$, $\mathbf{n}\sim\mathcal{N}(0,1)$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (f)
Baseline $\mathcal{S}_{\mathrm{interp}}$

$30$$50$$100$$130$Mean RMS Delay Spread
$\mathrm{\overline{DS}}_{\mathrm{rms},1}$ [ns]

Figure 6: RMS delay spread $\mathrm{\overline{DS}}_{\mathrm{rms},b}$, averaged
over all antennas in array $b=1$.

To obtain a generated dataset $\mathcal{S}_{\mathrm{gen}}$ that is comparable
to the test set $\mathcal{S}_{\mathrm{test}}$, we use $G$ to produce channel
realizations for all positions $\mathbf{x}^{(l)}$ present in
$\mathcal{S}_{\mathrm{test}}$. Evaluating the generator can be approached in
two ways:

  1. (I)

The noise vector $\mathbf{n}\sim\mathcal{N}(0,1)$ may be drawn once and then
remain fixed over all positions (conditions) $\mathbf{x}$:

$\displaystyle\mathcal{S}_{\mathrm{GAN,fix}}=\left\\{\left(G\left(\mathbf{x}^{(l)},\mathbf{n}\right),\mathbf{x}^{(l)}\right)\right\\}_{(\mathbf{H}^{(l)},\mathbf{x}^{(l)})\in\mathcal{S}_{\mathrm{test}}},$
$\displaystyle~{}\text{fixed}~{}\mathbf{n}\sim\mathcal{N}(0,1)$

  2. (II)

The noise vector $\mathbf{n}^{(l)}\sim\mathcal{N}(0,1)$ may be drawn randomly
for each position (condition) $\mathbf{x}^{(l)}$:

$\mathcal{S}_{\mathrm{GAN,var}}=\left\\{\left(G\left(\mathbf{x}^{(l)},\mathbf{n}^{(l)}\right),\mathbf{x}^{(l)}\right)\right\\}_{(\mathbf{H}^{(l)},\mathbf{x}^{(l)})\in\mathcal{S}_{\mathrm{test}},~{}\mathbf{n}^{(l)}\sim\mathcal{N}(0,1)}$

The first approach should lead to a higher level of spatial consistency in the
generated CSI data, whereas the second approach can reveal what the generator
has learned about the spatial channel distribution.

###  V-C Linear Interpolator Evaluation Procedure

To generate CSI with the linear interpolation baseline, which makes use of the
entire training set $\mathcal{S}_{\mathrm{train}}$ to parametrize the
interpolant, we evaluate the interpolant $F_{\mathrm{interp}}$ at the test set
locations, yielding the interpolated CSI dataset:

$\mathcal{S}_{\mathrm{interp}}=\left\\{\left(F_{\mathrm{interp}}\left(\mathbf{x}^{(l)}\right),\mathbf{x}^{(l)}\right)\right\\}_{(\mathbf{H}^{(l)},\mathbf{x}^{(l)})\in\mathcal{S}_{\mathrm{test}}}$

###  V-D Spatial Distribution Analysis

We consider the distribution of generated CSI relative to the two-dimensional
spatial coordinates, which are passed to the GAN as condition. If the CSI
generated by the GAN is similar to the test set CSI, this could indicate that
GAN-based channel models are useful as _inter-/extrapolators_. First, we
evaluate the generator output distribution according to approach (I), i.e.,
for fixed input noise $\mathbf{n}$. The corresponding received power
distributions $\lVert\mathbf{H}^{(l)}_{b}\rVert_{\mathrm{F}}^{2}$ for each
antenna array $b=1,\ldots,4$, both measured and generated for one particular
noise realization, are shown in Fig. 5. Clearly, the WGAN has learned to
produce CSI with some spatial consistency, and the spatial distribution of
generated and measured channels is similar, but not identical. The generated
channel realizations greatly depend on the choice of receiver antenna array
$b$. However, the neural network fails to produce CSI with true spatial
consistency. For example, the concentric circles around the antenna arrays
visible in the measured data (which are likely caused by a strong floor
reflection) are not visible in the generated power distribution.

We can also analyze the RMS delay spread for approach (I). The results for
array $b=1$ are shown in Fig. 6(c) and Fig. 6(d), for two different fixed
noise vectors $\mathbf{n}^{(1)}$ and $\mathbf{n}^{(2)}$: It can be seen that
the delay spread distributions depend on the noise vector and that they are
also spatially consistent, yet distinct from the measured distribution.
However, from Fig. 6(f), it is also clear that the linear interpolation
baseline produces CSI that much more closely resembles the test set if
accurate CSI prediction is desired. Though not shown in Fig. 6, these
observations also hold for antenna arrays $b=2,3,4$. If we choose a random
noise realization for every position $\mathbf{x}$ as in approach (II), the
spatial RMS delay spread distribution more closely resembles the test set, as
visible in Fig. 6(e).

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (a)
Test Set $\mathcal{S}_{\mathrm{test}}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (b)
GAN $\mathcal{S}_{\mathrm{GAN,fix}}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (c)
GAN $\mathcal{S}_{\mathrm{GAN,var}}$

$-10$$-6$$-2$$2$$-14$$-12$$-10$$-8$$-6$$-4$$-2$Coordinate
$\mathbf{x}_{1}~{}[\mathrm{m}]$Coordinate $\mathbf{x}_{2}~{}[\mathrm{m}]$ (d)
Interpolated $\mathcal{S}_{\mathrm{interp}}$

$-\pi/2$$0$$\pi/2$AoA Estimate $\hat{\alpha}^{(2)}$ [rad]

Figure 7: Spatial azimuth AoA estimate distribution for antenna array $b=2$
for various measured / generated CSI datasets. Note that the color map is
circular.

The observed spatial distribution of received powers and delay spreads
indicates that the generator was at least able to learn meaningful PDPs of the
wireless channel. To test if the generator has also learned a meaningful
distribution of phase information, we visualize the spatial distribution of
AoA estimates generated by MUSIC for measured and generated CSI data in Fig.
7. Clearly, the GAN has captured the relationship between location and CSI
phases to some extent, as indicated by the clearly visible gradient in the AoA
estimates in Fig. 7(c). The AoA estimate for the generated CSI corresponds to
the true AoA of the LoS regions in the dataset. However, the linear
interpolation baseline also clearly outperforms the GAN if the objective is to
generate true-to-life CSI data: The AoA estimates for the interpolated CSI in
Fig. 7(d) are almost identical to the AoA estimates generated for the test set
in Fig. 7(a). Fig. 7(b), which shows the AoA estimate distribution for CSI
generated by a GAN with fixed noise realization looks rather artificial: The
individual AoA estimates for the generated CSI are all plausible for the
specific position $\mathbf{x}$, but the spatial distribution is unrealistic,
with constant AoA in some areas and rapid changes over space elsewhere. This
may indicate that the generator has failed to produce CSI with meaningful
spatial consistency.

Comparing test set distribution and generated distributions of received power,
RMS delay spread and AoA estimates, we conclude that a GAN-based channel model
is not a good spatial interpolator in our scenario, as even a simple linear
interpolation baseline outperforms the GAN: Judging by RMS delay spread and
AoA estimates, the interpolated CSI is clearly more similar to the test set
than the GAN-generated CSI.

###  V-E Stochastic Distribution Analysis

$d_{\mathrm{JS}}$ | $\mathcal{S}_{\mathrm{train}}$ | $\mathcal{S}_{\mathrm{test}}$ | $\mathcal{S}_{\mathrm{GAN,var}}$ | $\mathcal{S}_{\mathrm{GAN,fix}}$ | $\mathcal{S}_{\mathrm{interp}}$ | Gauss.
---|---|---|---|---|---|---
$\mathcal{S}_{\mathrm{train}}$ | 0.000 | 0.041 | 0.066 | 0.093 | 0.049 | 0.160
$\mathcal{S}_{\mathrm{test}}$ | 0.041 | 0.000 | 0.082 | 0.115 | 0.067 | 0.149
$\mathcal{S}_{\mathrm{GAN,var}}$ | 0.066 | 0.082 | 0.000 | 0.053 | 0.054 | 0.162
$\mathcal{S}_{\mathrm{GAN,fix}}$ | 0.093 | 0.115 | 0.053 | 0.000 | 0.070 | 0.191
$\mathcal{S}_{\mathrm{interp}}$ | 0.049 | 0.067 | 0.054 | 0.070 | 0.000 | 0.156
Gauss. | 0.160 | 0.149 | 0.162 | 0.191 | 0.156 | 0.000
TABLE I: JSD between various RMS delay spread distributions.

In the previous section, we compared the spatial distribution of generated CSI
to the test set. If one is of the opinion that the GAN-based channel model
should be used as a _random generator_ for CSI that serves to generate CSI
data with a certain stochastic distribution (and not to interpolate CSI), this
is not a fair comparison. Therefore, we now compare the stochastic
distribution of measured and generated CSI.

For the following analysis, we focus on the RMS delay spread only and
disregard all spatial aspects. Of course, these analyses could also be
performed on other CSI properties like received power or angular-domain
information, with similar conclusions. We estimate the density function of the
RMS delay spreads over all antennas in all antenna arrays for measured and
generated CSI by binning delay spreads into 150 bins. In addition, we compare
the generated RMS delay spread distribution to randomly drawing from a
Gaussian probability density function (PDF), whose first- and second-order
moments are directly estimated from the training set’s RMS delay spread
distribution.

The resulting histograms and the parametrized Gaussian distribution are shown
in Fig. 8. The high degree of similarity between the measured and generated
histograms indicates that the GAN is able to approximate the measured channel
distribution well, but so is the linear interpolator. Furthermore, the
histogram indicates that the GAN not only learned the first and second-order
moments of the delay spread distribution, but also higher-order moments. If
one were instead to sample RMS delay spreads for a generative stochastic model
from an estimated Gaussian distribution, the resulting CSI distribution would
be a bad approximation of the realistic CSI. Please note that the RMS delay
spread distribution for $\mathcal{S}_{\mathrm{GAN,fix}}$ depends on the
particular random noise input $\mathbf{n}$, so the density shown here is only
exemplary.

To quantify the similarity of the densities in Fig. 8 more objectively, we
compute the JSDs between any two density functions, resulting in the JSD
distance matrix shown in Tab. I. Unsurprisingly, the RMS delay spread
distribution that most closely resembles the test set is the training set,
followed by the interpolated CSI. The RMS delay spread distribution generated
by the GAN is less similar than the interpolated CSI, but still much better
than sampling from a Gaussian distribution.

We conclude that in our scenario, the GAN-based channel model was able to
reproduce the statistical CSI distribution, at least when judged on the basis
of the delay spread of generated CSI. However, the CSI that the GAN generates
is highly specific to the considered setup and environment.

##  VI Criticism

$0$$25$$50$$75$$100$$125$$150$$175$$200$$0$$2\cdot 10^{-3}$$4\cdot
10^{-3}$$6\cdot 10^{-3}$$8\cdot 10^{-3}$$1\cdot 10^{-2}$$1.2\cdot
10^{-2}$$\mathrm{DS}_{\mathrm{rms},b,m_{r},m_{c}}~{}[\mathrm{ns}]$Density$\mathcal{S}_{\mathrm{train}}$$\mathcal{S}_{\mathrm{test}}$Gaussian

$0$$25$$50$$75$$100$$125$$150$$175$$200$$0$$2\cdot 10^{-3}$$4\cdot
10^{-3}$$6\cdot 10^{-3}$$8\cdot 10^{-3}$$1\cdot 10^{-2}$$1.2\cdot
10^{-2}$$\mathrm{DS}_{\mathrm{rms},b,m_{r},m_{c}}~{}[\mathrm{ns}]$Density$\mathcal{S}_{\mathrm{GAN,var}}$$\mathcal{S}_{\mathrm{GAN,fix}}$$\mathcal{S}_{\mathrm{interp}}$

Figure 8: Density functions of RMS delay spread distributions for various CSI
datasets and Gaussian approximation.

Neural networks show promising results for various communications-related
applications, including channel coding [18], channel charting [19], and
channel prediction [20]. What these applications have in common is that there
is a lot of training data available, but there is no good system model due to
the unknown statistical properties of the propagation channel. Clearly, under
these circumstances, machine learning can lead to significant improvements. On
the face of it, the same logic applies to GAN-based channel modeling, where
large CSI datasets are available for training. However, upon closer
inspection, this impression is deceptive: Just because there is a lot of
training data available does not necessarily mean that the training data is
suitable for training a useful generative model. In the following, we outline
why we believe GAN-based channel modeling (and generative neural network-based
channel modeling in general) is not ready for practical use in its current
form. In Sec. V-D we already concluded that GAN-based channel models are not
good interpolators. Here, we argue that in their current form, they are also
not good GSCM replacements and that significant architectural changes are
necessary for them to become useful.

###  VI-A Lack of generalizability

A GAN-based channel model is trained on one particular dataset or a collection
of CSI datasets captured with the same setup. The current system architecture
would not allow for training the same generator on datasets captured in
different environments with different antenna setups. Therefore, the GAN-based
channel model is always specific to the environment it was trained for and, in
contrast to a GSCM, cannot generalize to a different number of antennas,
antenna deployments, carrier frequencies or other environment parameters.

###  VI-B Lack of meaningful spatial consistency

Spatial consistency is an important concept in channel modeling: In GSCMs, it
can be ensured in various ways, e.g., by defining scattering clusters [3], or
by modeling spatial correlations of stochastic distributions like, for
example, spatially correlated shadow fading [21]. Judging by the results in
Sec. V-D, it may seem like the GAN ensures at least some level of spatial
consistency. This is the case in the sense that the GAN generates a CSI sample
that is plausible for the considered position (condition) $\mathbf{x}$.
However, it does not mean that the GAN ensures any meaningful spatial
correlation between two generated CSI samples $\hat{\mathbf{H}}^{(1)}$ and
$\hat{\mathbf{H}}^{(2)}$ at different locations $\mathbf{x}^{(1)}$ and
$\mathbf{x}^{(2)}$. In fact, we found that the level of spatial correlation
depends on the range to which the condition $\mathbf{x}$ is normalized. This
indicates that spatial correlation is not inferred from the training data, but
is merely an artifact of the network architecture and training procedure.

###  VI-C Lack of interpretability

One of the strengths of channel models is the ability to see how a designed
system would perform under different conditions. For example, a system
designer may want to understand how the system performs at different carrier
frequencies, for different antenna deployments, Doppler spreads, path loss
exponents or noise distributions. Other than in GSCMs, these properties are
not controlled in GAN-based channel models. If a developed system or algorithm
shows good or poor performance on a specific GAN-based channel model, it would
be difficult to understand why that is the case.

###  VI-D Availability of a dataset

The results in Sec. V-E indicate that using a GAN-based channel model is only
applicable if the objective is to have a channel model for the particular
environment that the measured (or simulated) training dataset was captured in.
In that case, a GAN can generate CSI that follows the distribution of the
training dataset. However, if the objective is to produce CSI that is as
similar as possible to a training set, one may just as well use the training
dataset itself. If the dataset is too small, one could use various classical
data augmentation techniques (e.g., adding noise, adding random global phase
rotations) or interpolation, rendering the generative model unnecessary.

##  VII Conclusion

We have demonstrated that a GAN can be used to learn the statistical
distribution of CSI in a measured dataset. GANs can capture higher-order
channel statistics, without requiring manual extraction of parameters, and, by
conditioning, can generate plausible CSI for arbitrary positions. However, as
presented here, GAN-based channel modeling still has many downsides compared
to GSCMs, some important ones being the inferior interpretability, lack of
meaningful spatial consistency and limited flexibility of the model: The
channel model is only valid for the particular environment and antenna
deployment of the dataset. For these reasons, we have come to the conclusion
that GAN-based channel models are not good interpolators and are also not
useful as general-purpose GSCMs in their current form. Achieving that would
likely require more, and more diverse training data and a different
architectural approach. For example we could imagine using GANs in conjunction
with raytracing-based channel models, where a generative model creates a radio
environment, and the raytracer computes spatially consistent CSI.

## References

  * [1] R. H. Clarke, “A statistical theory of mobile-radio reception,” _Bell system technical journal_ , vol. 47, no. 6, pp. 957–1000, 1968.
  * [2] ETSI 3rd Generation Partnership Project (3GPP), “Study on channel model for frequencies from 0.5 to 100 GHz (3GPP TR 38.901 version 16.1.0 Release 16),” ETSI, Tech. Rep., 2020.
  * [3] F. Burkhardt, S. Jaeckel, E. Eberlein, and R. Prieto-Cerdeira, “QuaDRiGa: a MIMO Channel Model for Land Mobile Satellite,” in _The 8th European Conference on Antennas and Propagation (EuCAP 2014)_. IEEE, 2014, pp. 1274–1278.
  * [4] J. Hoydis _et al._ , “Sionna RT: Differentiable Ray Tracing for Radio Propagation Modeling,” _arXiv preprint_ , Mar. 2023.
  * [5] I. Goodfellow _et al._ , “Generative Adversarial Nets,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [6] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein Generative Adversarial Networks,” in _International conference on machine learning_. PMLR, 2017, pp. 214–223.
  * [7] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. C. Courville, “Improved Training of Wasserstein GANs,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [8] Y. Yang, Y. Li, W. Zhang, F. Qin, P. Zhu, and C.-X. Wang, “Generative-Adversarial-Network-Based Wireless Channel Modeling: Challenges and Opportunities,” _IEEE Communications Magazine_ , vol. 57, 2019.
  * [9] H. Xiao, W. Tian, W. Liu, and J. Shen, “ChannelGAN: Deep learning-based channel modeling and generating,” _IEEE Wireless Communications Letters_ , vol. 11, no. 3, pp. 650–654, 2022.
  * [10] T. Orekondy, A. Behboodi, and J. B. Soriaga, “MIMO-GAN: Generative MIMO Channel Modeling,” in _ICC 2022-IEEE International Conference on Communications_ , 2022, pp. 5322–5328.
  * [11] W. Xie, M. Xiong, Z. Yang, W. Liu, L. Fan, and J. Zou, “Real and fake channel: GAN-based wireless channel modeling and generating,” _Physical Communication_ , p. 102214, 2023.
  * [12] U. Sengupta, C. Jao, A. Bernacchia, S. Vakili, and D. Shiu, “Generative Diffusion Models for Radio Wireless Channel Modelling and Sampling,” _arXiv preprint arXiv:2308.05583_ , 2023.
  * [13] J. Juhava, “Wireless channel modeling using generative machine learning models,” Master’s thesis, Aalto University School of Engineering, 2023.
  * [14] T. J. O’Shea, T. Roy, and N. West, “Approximating the Void: Learning Stochastic Channel Models from Observation with Variational Generative Adversarial Networks,” in _2019 International Conference on Computing, Networking and Communications (ICNC)_. IEEE, 2019.
  * [15] S. Dörner, M. Henninger, S. Cammerer, and S. ten Brink, “WGAN-based Autoencoder Training Over-the-air,” in _2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_. IEEE, 2020, pp. 1–5.
  * [16] F. Euchner, M. Gauger, S. Dörner, and S. ten Brink, “A Distributed Massive MIMO Channel Sounder for ”Big CSI Data”-driven Machine Learning,” in _25th Workshop on Smart Antennas_ , 2021.
  * [17] F. Euchner and M. Gauger, “CSI Dataset dichasus-cf0x: Distributed Antenna Setup in Industrial Environment, Day 1,” 2022. [Online]. Available: https://doi.org/doi:10.18419/darus-2854
  * [18] T. Gruber, S. Cammerer, J. Hoydis, and S. ten Brink, “On Deep Learning-Based Channel Decoding,” in _2017 51st annual conference on information sciences and systems (CISS)_ , 2017, pp. 1–6.
  * [19] C. Studer, S. Medjkouh, E. Gonultaş, T. Goldstein, and O. Tirkkonen, “Channel Charting: Locating Users Within the Radio Environment Using Channel State Information,” _IEEE Access_ , vol. 6, 2018.
  * [20] W. Jiang and H. D. Schotten, “Neural Network-Based Fading Channel Prediction: A Comprehensive Overview,” _IEEE Access_ , vol. 7, 2019.
  * [21] X. Cai and G. B. Giannakis, “A Two-Dimensional Channel Simulation Model for Shadowing Processes,” _IEEE Transactions on Vehicular Technology_ , vol. 52, no. 6, pp. 1558–1567, 2003.

  *[AWGN]: additive white Gaussian noise
  *[3GPP]: 3rd Generation Partnership Project
  *[TR]: technical report
  *[GSCM]: geometry-based stochastic channel model
  *[QuaDRiGa]: Quasi Deterministic Radio Channel Generator
  *[GAN]: generative adversarial network
  *[CSI]: channel state information
  *[WGAN]: Wasserstein GAN
  *[GP]: gradient penalty
  *[GANs]: generative adversarial network
  *[VAE]: variational autoencoder
  *[LTI]: linear time-invariant
  *[MIMO]: multiple-input multiple-output
  *[GSCMs]: geometry-based stochastic channel model
  *[UE]: user equipment
  *[CIR]: channel impulse response
  *[TDL]: tapped delay line
  *[DICHASUS]: Distributed Channel Sounder by University of Stuttgart
  *[mMIMO]: massive multiple-input multiple-output
  *[OFDM]: orthogonal frequency division multiplex
  *[LoS]: line-of-sight
  *[RMS]: root mean square
  *[AoA]: angle of arrival
  *[JSD]: Jensen-Shannon distance
  *[UPA]: uniform planar array
  *[MUSIC]: MUltiple SIgnal Classification
  *[KL]: Kullback-Leibler
  *[PDPs]: power delay profile
  *[PDF]: probability density function
  *[JSDs]: Jensen-Shannon distance


# Multi-tracing the primordial Universe with future surveys

Mponeng Kopana1,a, Sheean Jolicoeur1,2,b, Roy Maartens1,3,4,c

###### Abstract

The fluctuations generated by Inflation are Gaussian in the simplest models,
but may be non-Gaussian in more complex models. Since these fluctuations seed
the large-scale structure, any primordial non-Gaussianity of local type will
leave an imprint on the tracer power spectrum on ultra-large scales. In order
to combat the problem of growing cosmic variance on these scales, we use a
multi-tracer analysis that combines different tracers of the matter
distribution to maximise any primordial signal. In order to illustrate the
advantages that can be brought by multi-tracing future surveys, we consider
two pairs of (galaxy survey, 21cm intensity mapping survey), one at high
redshift ($1\leq z\leq 2$) and one at very high redshift ($2\leq z\leq 5$).
The 21cm surveys are in interferometer mode, and are idealised versions of
HIRAX and PUMA. We implement foreground avoidance filters and use non-trivial
models of the interferometer thermal noise. The galaxy surveys are idealised
versions of Euclid and MegaMapper. Via a simple Fisher forecast we illustrate
the potential of the multi-tracer. Our results show a $\sim 20-70\%$
improvement in precision on local primordial non-Gaussianity from the multi-
tracer.

##  1 Introduction

Galaxy surveys in combination with cosmic microwave background (CMB) surveys
have laid the foundations for precision cosmology. The next generation of
surveys will deliver even higher precision by covering larger sky areas and
deeper redshifts with better sensitivity. Each of these surveys is destined to
make major advances in the precision and accuracy of the cosmological
parameters of the standard LCDM model. At the same time, future surveys will
look for signals beyond LCDM. One such signal arises if there is significant
non-Gaussianity in the primordial perturbations that seed the CMB fluctuations
and the large-scale structure.

Primordial non-Gaussianity (PNG) is a key probe of Inflation, which is
currently the best framework that we have for the generation of seed
fluctuations. The local type of PNG, described by the parameter $f_{\rm NL}$,
if found to be nonzero, will rule out the simplest Inflation models [1].
Furthermore, if $0<|f_{\rm NL}|\lesssim 1$, then many other models can also be
ruled out [2, 3]. In order to achieve this, a precision of $\sigma(f_{\rm
NL})<1$ is required.

The best current 1$\sigma$ constraint comes from the Planck survey, via the
CMB bispectrum [4]:

$f_{\rm{NL}}=-0.9\pm 5.1\,.$ (1.1)

The CMB power spectrum does not carry an imprint of local PNG, and the same is
true for the matter power spectrum. However, the power spectrum of a biased
tracer, such as galaxies, does contain a signal of local PNG on ultra-large
scales – due to the phenomenon of scale-dependent bias [5, 6]. Constraints
from future CMB surveys will not be able to reach $\sigma(f_{\rm NL})<1$ due
to cosmic variance. In the case of galaxy surveys, the current best
constraints are $\sigma(f_{\rm NL})\gtrsim 20$. Future galaxy surveys will
improve on this, potentially reaching the CMB level of precision and beyond
[7, 8]. Although galaxy power spectrum constraints on $f_{\rm NL}$ will be a
powerful independent check on the CMB results, single-tracer constraints are
unlikely to achieve $\sigma(f_{\rm NL})<1$, again because of cosmic variance,
which is highest on precisely the scales where the local PNG signal is
strongest.

Since the local PNG signature in large-scale structure depends on the
particular tracer, it is possible to improve the precision on $f_{\rm NL}$ by
combining the power spectra of different tracers. This multi-tracer approach
can effectively remove cosmic variance [9, 10, 11, 12]. As a result, the
multi-tracer of future galaxy surveys with intensity mapping of the 21cm
emission of neutral hydrogen (HI IM surveys), can deliver significant
improvements in the precision on $f_{\rm NL}$ (e.g. [13, 14, 15]) and in some
cases it can in principle reach $\sigma(f_{\rm NL})<1$ (e.g. [16, 17, 18, 19,
20]). The multi-tracer can also mitigate some of the observational systematics
in single-tracer surveys.

Here we follow [15] in considering pairs made of one spectroscopic galaxy
survey and one HI intensity mapping survey. However, we consider HI IM surveys
in interferometer mode, rather than the single-dish mode surveys considered in
[15]. In addition, we consider a pair at much higher redshift than in [15]. We
use nominal simplified surveys since we perform Fisher forecasting which does
not take into account most systematics. Our main aim is to demonstrate the
advantage of the multi-tracer at high to very high redshifts with combined
galaxy–HI probes, rather than to make forecasts for specific surveys.

The paper is structured as follows. In section 2, we provide a brief
demonstration of multi-tracer power spectrum of two tracers. In section 3, we
presented the characteristics of the galaxy and HI intensity mapping surveys
we are using. We discuss the survey’s instrumentation and brief noise
requirements in section 3.1. The avoidance of foreground contamination for HI
intensity mapping is presented in section 3.2. We also discuss the limited
range of scales needed to obtain the cosmological information from both
spectroscopic and radio surveys in section 3.3. The Fisher forecast is
discussed in Section 4, along with how we combined the data from the two
surveys. Section 5 provides a summary of our main conclusions.

##  2 Multi-tracer power spectra

At linear order in perturbations, the Fourier number density or brightness
temperature contrast of tracer $A$ is

$\delta_{A}(z,\bm{k})=\Big{[}b_{A}(z)+f(z)\mu^{2}\Big{]}\delta(z,\bm{k})\;,$
(2.1)

where $\delta$ is the matter density contrast, $b_{A}$ is the (Gaussian)
clustering bias for the tracer population $A$, $\mu=\bm{\hat{k}}\cdot\bm{{n}}$
is the projection along the line-of-sight direction $\bm{n}$, and
$f=\Omega_{m}^{\gamma}$ is the linear growth rate, where $\gamma=0.55$ for
LCDM. We use a single-parameter model for the biases,

$b_{A}(z)=b_{A0}\,\beta_{A}(z)\;,$ (2.2)

where we assume that the redshift evolution $\beta_{A}(z)$ is known. The power
spectra for 2 tracers are given by

$\big{\langle}\delta_{A}(z,\bm{k})\delta_{B}(z,\bm{k}^{\prime})\big{\rangle}=(2\pi)^{3}P_{AB}(z,\bm{k})\delta^{\mathrm{D}}(\bm{k}+\bm{k}^{\prime})\;,$
(2.3)

so that

$P_{AB}=\big{(}b_{A}+f\mu^{2}\big{)}\big{(}b_{B}+f\mu^{2}\big{)}P\;,$ (2.4)

where $P$ is the matter power spectrum (computed with CLASS [21]). When $B=A$
we use the shorthand $P_{A}\equiv P_{AA}$.

In the presence of local PNG, the bias acquires a scale-dependent correction:

$\hat{b}_{A}(z,k)=b_{A}(z)+b_{A\phi}(z)\frac{f_{\rm
NL}}{\mathcal{M}(z,k)}\,,\quad\mathcal{M}(z,k)=\frac{2}{3\Omega_{m0}H_{0}^{2}}\frac{D(z)}{g_{\rm
in}}T(k)\,k^{2}\;,$ (2.5)

where $T(k)$ is the matter transfer function (normalized to 1 on very large
scales), $D$ is the growth factor (normalized to 1 today) and $g_{\rm
in}=(1+z_{\rm in})D_{\rm in}$ is the initial growth suppression factor, with
$z_{\rm in}$ deep in the matter era. For $\Lambda$CDM, $g_{\rm in}=1.274$. The
modelling of the non-Gaussian bias factor $b_{A\phi}$ is in principle
determined from halo simulations, but currently it remains uncertain [22, 23,
24, 25, 26, 27]. A simplified model is given by

$b_{A\phi}(z)=2\delta_{\rm c}\big{[}b_{A}(z)-p_{A}\big{]}\,,$ (2.6)

where the critical collapse density is given by $\delta_{\rm c}=1.686\,$ and
$p_{A}$ is a halo-dependent constant. For the simplest case of a universal
halo model, $p_{A}=1$ for all tracers $A$ [5, 6]. We use (2.6) with the
following $p_{A}$ (derived from simulations): for galaxies chosen by stellar
mass [22],

$p_{g}\approx 0.55\;,$ (2.7)

and for HI intensity mapping [23]

$p_{\rm HI}\approx 1.25\;.$ (2.8)

##  3 Surveys

We consider two pairs of nominal spectroscopic surveys, each pair with the
same redshift range and with sky area $\Omega_{\rm sky}=10,000~{}{\rm
deg}^{2}$:

  * •

galaxy survey ‘H$\alpha$’, and HI intensity mapping survey ‘H’, $1\leq z\leq
2$ (similar to $Euclid$ [28] and HIRAX [29, 30] respectively);

  * •

galaxy survey ‘LBG’ and HI intensity mapping survey ‘P’, $2\leq z\leq 5$
(similar to MegaMapper [31] and PUMA [32] respectively).

The Gaussian clustering biases are as follows:

$\displaystyle b_{{\rm H}\alpha}(z)$ $\displaystyle=b_{{\rm H}\alpha 0}(1+z)$
$\displaystyle\text{fiducial $b_{{\rm H}\alpha 0}=0.7$
\cite[cite]{[\@@bibref{Number}{Merson:2019vfr}{}{}]}},$ (3.1) $\displaystyle
b_{\rm LBG}(z)$ $\displaystyle=b_{\rm
LBG0}\big{(}1+1.155\,z+0.155\,z^{2}\big{)}$ $\displaystyle\text{fiducial
$b_{\rm LBG0}=0.710$ \cite[cite]{[\@@bibref{Number}{Sailer:2021yzm}{}{}]}},$
(3.2) $\displaystyle b_{\rm H}(z)$ $\displaystyle=b_{\rm P}(z)=b_{\rm
HI}(z)=b_{\rm HI0}\big{(}1+0.823z-0.0546z^{2}\big{)}$
$\displaystyle\text{fiducial $b_{\rm HI0}=0.842$
\cite[cite]{[\@@bibref{Number}{Villaescusa-Navarro:2018vsg}{}{}]}}.$ (3.3)

The background comoving number densities for H$\alpha$ [36] and LBG [37] are
taken as

$\displaystyle\bar{n}_{{\rm H}\alpha}(z)$
$\displaystyle=\big{(}3.63\,z^{-0.910}{\rm e}^{0.402\,z}-4.14\big{)}\times
10^{-3}h^{3}~{}\mathrm{Mpc^{-3}}\,,$ (3.4) $\displaystyle\bar{n}_{\rm LBG}(z)$
$\displaystyle=\big{(}2.43-2.21\,z+0.773\,z^{2}-0.122\,z^{3}+0.00733\,z^{4}\big{)}\times
10^{-2}h^{3}\,\mathrm{Mpc^{-3}}\,.$ (3.5)

For HI IM, the background brightness temperature is modelled as [38]:

$\bar{T}_{\rm HI}(z)=0.0559+0.2324\,z-0.0241\,z^{2}~{}~{}\mathrm{mK}\,.$ (3.6)

Equations (3.1)–(3.6) are shown in Figure 1.

Figure 1: For the galaxy and HI IM surveys: clustering biases (left panel);
comoving number densities (right panel, red, left-hand $y$-axis) and HI IM
temperature (blue, right-hand $y$-axis).

###  3.1 Noise

For galaxy surveys, $g={\rm H}\alpha$ or LBG, the shot noise is assumed to be
Poissonian, so that the total observed auto-power is

$\tilde{P}_{g}(z,k,\mu)=P_{g}(z,k,\mu)+\frac{1}{\bar{n}_{g}(z)}\;.$ (3.7)

In HI IM surveys, HI = H or P, the dominant noise on linear scales is the
thermal (or instrumental) noise [39, 35], so that shot noise can safely be
neglected and the observed HI IM power spectrum is

$\tilde{P}_{\rm HI}(z,k,\mu)=P_{\rm HI}(z,k,\mu)+P_{\rm
HI}^{\mathrm{therm}}(z,k,\mu)\;,$ (3.8)

where $P_{\rm HI}^{\mathrm{therm}}$ depends on the model of the sky
temperature in the radio band, the survey specifications and the survey mode
(single-dish or interferometer). For the interferometer (IF) mode that we
consider, the thermal noise power spectrum is [40, 41, 42, 43]:

$P_{\rm
HI}^{\mathrm{therm}}(z,k_{\perp})=\frac{\Omega_{\mathrm{sky}}}{2\nu_{21}t_{\mathrm{tot}}}\,\frac{(1+z)r(z)^{2}}{\mathcal{H}(z)}\,\left[\frac{T_{\mathrm{sys}}(z)}{\bar{T}_{\rm
HI}(z)}\right]^{2}\,\left[\frac{\lambda_{21}(1+z)}{A_{\mathrm{e}}}\right]^{2}\frac{1}{n^{\mathrm{phys}}_{\mathrm{b}}(z,k_{\perp})}{\frac{1}{\theta_{\rm
fov}(z)^{2}}}\;,$ (3.9)

where $k_{\perp}=k\sqrt{1-\mu^{2}}$ is the transverse mode,
$\nu_{21}=1420\,\mathrm{MHz}$ is the rest-frame frequency of the
$\lambda_{21}=21\,$cm emission, $t_{\mathrm{tot}}$ is the total observing
time, $r$ is the radial comoving distance, $\mathcal{H}$ is the conformal
Hubble rate, $n^{\mathrm{phys}}_{\mathrm{b}}$ is the physical baseline density
of the array distribution (see below). The effective beam area is

$A_{\mathrm{e}}=0.7\,\frac{\pi}{4}D_{\mathrm{d}}^{2}\;,$ (3.10)

with dish diameter $D_{\mathrm{d}}$ and the field of view is

$\theta_{\rm fov}(z)=1.22\,\frac{\lambda_{21}(1+z)}{D_{\rm d}}\;.$ (3.11)

For a given baseline of length $L$, we use the fitting formula for close-
packed arrays given in [42] (see also [44, 43]):

$n^{\mathrm{phys}}_{\mathrm{b}}(z,k_{\perp})=\left(\frac{N_{\mathrm{s}}}{D_{\mathrm{d}}}\right)^{2}\frac{a+b(L/L_{\mathrm{s}})}{1+c(L/L_{\mathrm{s}})^{d}}\exp{\big{[}-(L/L_{\mathrm{s}})^{e}\big{]}}~{}~{}\text{with}~{}~{}L=\frac{k_{\perp}r(z)}{2\pi}\lambda_{21}(1+z)\,,$
(3.12)

where $L_{\mathrm{s}}=N_{\mathrm{s}}D_{\mathrm{d}}$. H is a square close-
packed array, while P is a hexagonal close-packed array with 50% fill factor.
This means that half the dish sites are randomly empty, so the array is
equivalent in size to that of twice the number of elements but with a quarter
of the baseline density [44]. Then it follows that $N_{\rm s}^{2}=N_{\rm d}$
for H and $N_{\rm s}^{2}=2N_{\rm d}$ for P [42, 45]. The values of the
parameters are given in Table 1.

The baseline densities for the H and P surveys are shown in Figure 2. Note
that the thermal noise, which is $\propto 1/n_{\rm b}^{\rm phys}$, tends to
infinity at the maximum baseline $D_{\rm max}$.

Table 1: Specifications of HI IM surveys.  Survey |  $N_{\mathrm{d}}$ |  $N_{\mathrm{s}}$ |  $t_{\rm tot}$ |  $T_{\mathrm{d}}$ |  $D_{\mathrm{d}}$ |  $D_{\mathrm{max}}$
---|---|---|---|---|---|---
|  |  |  [$10^{3}\rm{hr}$] |  [K] |  [m] |  [m]
H256 (1024) | 256 (1024) | 16 (32) | 12 | 50 | 6 | 141 (282)
P5k (32k) | 5k (32k) | 100 (253) | 20 | 93 | 6 | 648 (1640)
Fitting |  $a$ |  $b$ |  $c$ |  $d$ |  $e$ |
parameters |  |  |  |  |  |
H256 (1024) | 0.4847 | -0.3300 | 1.3157 | 1.5974 | 6.8390 |
P5k (32k) | 0.5698 | -0.5274 | 0.8358 | 1.6635 | 7.3177 |

Figure 2:  Physical baseline density for HI IM surveys H1024 (_left_) and P32k
(_right_).

Figure 3:  Redshift evolution of galaxy shot noise (_left_) and HI IM thermal
noise times physical baseline density (_right_).

The system temperature is modelled as [42]:

$T_{\mathrm{sys}}(z)=T_{\rm d}(z)+T_{\rm sky}(z)=T_{\rm
d}+2.7+25\bigg{[}\frac{400\,\mathrm{MHz}}{\nu_{21}}(1+z)\bigg{]}^{2.75}~{}\mathrm{K},$
(3.13)

where $T_{\rm d}$ is the dish receiver temperature, given in Table 1.

Figure 4: 2D colour plot of thermal noise for HI IM surveys: H256 (_top
left_), H1024 (_top right_), P5k (_bottom left_) and P32k (_bottom right_).

The cross-noise power spectrum between HI IM and galaxy surveys depends on the
extent of overlapping halos that host the HI emitters and the galaxies.
Typically it is neglected, as motivated in [46, 47]. Then the observed cross-
power spectra are

$\tilde{P}_{g{\rm HI}}=P_{g{\rm HI}}\,.$ (3.14)

The noise auto-power for all surveys is shown in Figure 3.

###  3.2 HI IM foreground avoidance

HI IM surveys are contaminated by foregrounds that are much larger than the HI
signal. The cleaning of foregrounds from the HI signal relies on the fact that
the foregrounds are nearly spectrally smooth (see e.g. [48, 49] and references
there in for recent advances). However, the cleaning is not perfect, and in
some regions of Fourier space the signal loss is large. For simplified Fisher
forecasts, we assume that cleaning is successful, except in the known regions
of Fourier space where cleaning is expected to be least effective. Then we use
foreground filters that avoid these regions. There are two regions of
foreground avoidance for IF-mode surveys.

  1. 1.

We filter out long-wavelength radial modes with $|k_{\parallel}|\lesssim
k_{\rm fg}$, where $k_{\rm fg}$ is a suitable minimum radial mode for future
surveys. We follow [50] and use the radial damping factor

${\cal D}_{\rm fg}(k,\mu)=1-\exp\bigg{[}-\bigg{(}\frac{k_{\parallel}}{k_{\rm
fg}}\bigg{)}^{\\!\\!{2}\,}\bigg{]}~{}~{}\text{with}~{}~{}k_{\rm
fg}=0.005~{}h{\rm Mpc}^{-1}\;.$ (3.15)

  2. 2.

A physical baseline in the IF array probes different angular scales
$k_{\perp}$ at different frequencies (and hence different $z$). As a
consequence, monochromatic emission from a foreground point source can
contaminate the signal in a wedge-shaped region of Fourier space [51, 52].
This contamination occurs within $N_{\rm w}$ primary beams of a pointing,
where $N_{\rm w}=0$ is the ideal case of no contamination, which occurs in
principle if calibration of the interferometer is near-perfect [53]. Currently
this precision in calibration and stability is not achievable. In order to
avoid the wedge region we therefore impose the filter [51, 52, 54, 41, 42]

$\left|k_{\parallel}\right|>A_{\mathrm{w}}(z)\,k_{\perp}\quad\text{where}\quad
A_{\mathrm{w}}(z)=r(z)\mathcal{H}(z)\sin\big{[}0.61N_{\mathrm{w}}\,\theta_{\mathrm{fov}}(z)\big{]}\;.$
(3.16)

In our forecasts we take $N_{\mathrm{w}}=1$.

Figure 5 illustrates ${\cal D}_{\rm fg}$ and $A_{\mathrm{w}}$. Foreground
filters lead to the following modifications of the power spectra:

$\displaystyle P_{A{\rm HI}}(z,{k},\mu)$ $\displaystyle\rightarrow$
$\displaystyle{\cal D}_{\rm
fg}(k,\mu)\,\Theta\big{(}\left|k_{\parallel}\right|-A_{\mathrm{w}}(z)k_{\perp}\big{)}\,P_{A{\rm
HI}}(z,{k},\mu)\,,$ (3.17)

where $A=\,$HI or H$\alpha$ or LBG, and $\Theta$ is the Heaviside step
function.

Figure 5: Radial foreground damping function (_left_). Schematic to illustrate
the foreground wedge; the wedge angle increases with redshift (_right_).

###  3.3 Maximum and minimum scales

Since the local PNG signal in power spectra is on ultra-large scales, we only
require linear perturbation theory. This means that we can exclude scales
where linear perturbations break down, i.e. $k<k_{\mathrm{max}}$, where [55,
56, 57, 58]

$k_{\mathrm{max}}(z)=k_{\mathrm{max},0}\big{(}1+z\big{)}^{2/(2+n_{s})}\;.$
(3.18)

Then we choose a very conservative value of
$k_{\mathrm{max},0}=0.08\;h\,\mathrm{Mpc}^{-1}$.

The minimum wavenumber $k_{\mathrm{min}}$ in each redshift bin is the
fundamental mode $k_{\mathrm{f}}$, given by the comoving volume of the bin of
width $\Delta z$ centered at $z$:

$k_{\mathrm{min}}(z)=k_{\mathrm{f}}(z)=\frac{2\pi}{V(z)^{1/3}}\quad\text{with}\quad
V(z)=\frac{\Omega_{\mathrm{sky}}}{3}\Big{[}r\big{(}z+\Delta
z/2\big{)}^{3}-r\big{(}z-\Delta z/2\big{)}^{3}\Big{]}\;.$ (3.19)

The minimum and maximum scales are shown in Figure 6.

Figure 6: Minimum and maximum wavenumbers (left) and comoving volume (right)
over the redshift range of the surveys.

##  4 Fisher forecast

For a given redshift bin, the multi-tracer Fisher matrix for the combination
of two dark matter tracers $g$ and HI is

$F_{\alpha\beta}^{\bm{P}}=\sum_{\mu=-1}^{+1}\,\sum_{k=k_{\mathrm{min}}}^{k_{\mathrm{max}}}\,\partial_{\alpha}\,\bm{P}\cdot\mathrm{Cov}(\bm{P},\bm{P})^{-1}\cdot\partial_{\beta}\,\bm{P}^{\mathrm{T}}\;,$
(4.1)

where $\partial_{\alpha}=\partial\,/\,\partial\vartheta_{\alpha}$. Here
$\bm{P}$ is the data vector of the power spectra,

$\bm{P}=\big{(}P_{g}\,,\,P_{g{\rm HI}}\,,\,P_{\rm HI}\big{)}\;.$ (4.2)

Note that the noise is excluded since it is independent of the cosmological
and nuisance parameters. The noise enters the covariance [59, 60, 61, 62]:

$\mathrm{Cov}(\bm{P},\bm{P})=\frac{k_{\rm f}^{3}}{4\pi k^{2}\Delta
k}\,\frac{2}{\Delta\mu}\,\begin{pmatrix}\tilde{P}_{g}^{2}&&\tilde{P}_{g}\tilde{P}_{g{\rm
HI}}&&\tilde{P}_{g{\rm HI}}^{2}\\\ &&&&\\\ \tilde{P}_{g}\tilde{P}_{g{\rm
HI}}&&\frac{1}{2}\big{(}\tilde{P}_{g}\tilde{P}_{\rm HI}+\tilde{P}_{g{\rm
HI}}^{2}\big{)}&&\tilde{P}_{\rm HI}\tilde{P}_{g{\rm HI}}\\\ &&&&\\\
\tilde{P}_{g{\rm HI}}^{2}&&\tilde{P}_{\rm HI}\tilde{P}_{g{\rm
HI}}&&\tilde{P}_{\rm HI}^{2}\end{pmatrix}\;.$ (4.3)

Here $\Delta k$ and $\Delta\mu$ are the bin-widths for $k$ and $\mu$, which we
choose following [63, 64]:

$\Delta\mu=0.04\quad\text{and}\quad\Delta k=k_{\mathrm{f}}\;.$ (4.4)

In this paper our focus is on the improvements that the multi-tracer can
deliver for pairs (galaxy, HI IM) of high and very high redshift surveys. We
assume that most of the cosmological parameters have been measured by surveys
that focus on medium to small scales. Since $f_{\mathrm{NL}}$ affects the
amplitude and shape of the power spectrum on ultra-large scales, we include in
the Fisher analysis the cosmological parameters $A_{s}$ and $n_{s}$ that
directly affect these properties of the power spectrum. We further assume that
(a) $\bar{T}_{\rm HI}$, (b) the redshift evolution of the Guassian clustering
biases $b_{A}$, and (c) the degeneracies between $b_{A}$ and $A_{s}$ and
between $\bar{T}_{\rm HI}$ and $A_{s}$, have been broken by other surveys that
focus on medium to small scales. Then we consider the following cosmological
and nuisance parameters:

$\vartheta_{\alpha}=\big{(}f_{\mathrm{NL}},\;A_{s},\;n_{s},\;b_{g0},\;b_{{\rm
HI}0}\big{)}\,.$ (4.5)

The fiducial values for the LCDM cosmological parameters are taken from Planck
[65] and we assume $\bar{f}_{\rm NL}=0$.

Derivatives for $f_{\mathrm{NL}}$, $A_{s}$ and $b_{A0}$ are computed
analytically, while the $n_{s}$ derivative is computed using the 5-point
stencil method, with step-length 0.1 [66].

The Fisher forecast results are shown in Figure 7 and Figure 8, which display
the $1\sigma$ error contours for the multi-tracer pairs $(A,B)=({\rm
H}\alpha$, H) and (LBG, P), where we have marginalised over the bias nuisance
parameters. The numerical values of $\sigma(f_{\rm NL})$ are shown in Table 2.

Figure 7: 1$\sigma$ marginalised contours for cosmological parameters from the
H$\alpha$ and H surveys.

Figure 8: 1$\sigma$ marginalised contours for cosmological parameters from the LBG and P surveys. Table 2: Marginalised errors on $f_{\mathrm{NL}}$. Final column shows the percentage improvement in precision for each single tracer from the multi-tracer.  Surveys $A,\,A\otimes B$ | $\sigma(f_{\mathrm{NL}})$ | $\big{[}\sigma(A)-\sigma(A\otimes B)\big{]}/\sigma(A)$
---|---|---
|  | $\times 100\%$
H$\alpha$ | 3.83 | 18 (24)
LBG | 0.787 | 23 (25)
H256 (1024) | 9.42 (7.74) | 68 (62)
P5k (32k) | 1.31 (1.24) | 54 (52)
H$\alpha$ $\otimes$ H256 (1024) | 3.13 (2.92) |
LBG $\otimes$ P5k (32k) | 0.606 (0.591) |

##  5 Conclusion

Our Fisher analysis shows the effectiveness of combining tracers to enhance
the precision in future measurements of local PNG. The multi-tracer also
improves precision on $n_{s}$ and $A_{s}$. It is noticeable that the HI IM
surveys provide stronger constraints on $n_{s},A_{s}$ than the galaxy surveys,
while the reverse is true for $f_{\rm NL}$. The reason is that the
$n_{s},A_{s}$ signals are mainly dependent on medium scales, while the $f_{\rm
NL}$ signal is based on ultra-large scales. The HI IM IF-mode surveys provide
better constraints on medium scales than the galaxy surveys. However, on
ultra-large scales these surveys lose constraining power due to the foreground
filters.

Our assumption on the halo-dependent PNG bias, i.e. that $b_{A\phi}\propto
b_{A}-p_{A}$ as in (2.7) and (2.8), also affects galaxies and HI IM
differently. In the simplest universality model, $p_{A}=1$, there is no
difference between galaxies and HI IM in the form of the PNG bias. In our
case, since $p_{g}<1$, it follows that $b_{A\phi}$ is larger than in the $p=1$
model. Consequently, the constraints for galaxy surveys on $f_{\rm NL}$ are
improved compared to $p=1$. By contrast, $p_{\rm HI}>1$ means that $f_{\rm
NL}$ constraints are degraded for HI IM surveys compared to $p=1$.

Single-tracer Fisher constraints for surveys similar to our mock surveys have
been computed previously: e.g., [67] for Euclid; [34] for MegaMapper and PUMA;
[44] for HIRAX and PUMA. These constraints are qualitatively compatible with
ours – taking into account the following differences:

  * •

Previous work uses the universality model for $b_{A\phi}$, while we use (2.7)
and (2.8).

  * •

Our mock surveys have differences in the redshift range, sky area and
observing time (for HI IM) compared to the surveys considered in previous
work.

  * •

There are differences amongst the previous works, and with our paper, in the
cosmological and nuisance parameters that are used and marginalised over.

Previous works, to our knowledge, do not treat the multi-tracer pairs that we
analyse, so that comparisons of our multi-tracer results with previous results
is not possible. Table 2 shows (final column) the improvement in precision
from the multi-tracer combinations for each single tracer case: this ranges
from $\sim 20$ to $\sim 70\%$. The best single-tracer and multi-tracer
$f_{\mathrm{NL}}$ precision is delivered at the highest redshifts, by LBG and
LBG$\,\otimes\,$P (Figure 8). The LBG survey on its own already achieves
$\sigma(f_{\mathrm{NL}})<1$ (this is compatible with [34]), while P
constraints (which are compatible with [44]) are weaker because of foreground
contamination on ultra-large scales. When LBG and P are combined, we achieve
an improvement in $\sigma(f_{\mathrm{NL}})$ of $\gtrsim 50\%$ when compared to
the single-tracer P, and $\gtrsim 20\%$ when compared to the single-tracer
LBG.

Regarding foreground avoidance filters, we have used the following values for
the radial long-wavelength mode limit and the wedge primary beam number:

$\displaystyle k_{\rm fg}=0.005\,h{\rm Mpc}^{-1}\,,\quad N_{\rm w}=1\,.$ (5.1)

Although this is fairly optimistic compared to some previous treatments, it
should be noted that further advances in foreground treatments will emerge. In
particular:

  * •

density reconstruction methods, such as forward model reconstruction, are
being developed to retrieve the long wavelength modes lost to foregrounds (e.g
[68, 58, 69]);

  * •

advances in the stability and calibration of interferometers will reduce the
impact of the wedge and in principle will be able in future to reduce its
impact to negligible.

Our Fisher forecasts are of course highly simplified compared to observational
and data reality. In practice there are numerous systematics to be modelled
and there are major data constraints on the simulations necessary to build a
data pipeline. Nevertheless, the simplified Fisher analysis has shown that
multi-tracing galaxy and intensity HI IM surveys at high to very high
redshifts is worth exploring further. Furthermore, the multi-tracer will
alleviate the impact of nuisance parameters and systematics in the galaxy and
HI IM samples.

Finally, we note that our forecasts use a plane-parallel approximation in the
Fourier power spectra. Future work will relax this assumption by including
wide-angle effects (see e.g. [70, 71, 72]).

Acknowledgements
We thank Dionysios Karagiannis, Stefano Camera and Lazare Guedezounme for very
helpful comments. We are supported by the South African Radio Astronomy
Observatory (SARAO) and the National Research Foundation (Grant No. 75415).

## References

  * [1] A. Achúcarro et al., Inflation: Theory and Observations, arXiv:2203.08128.
  * [2] R. de Putter, J. Gleyzes, and O. Doré, Next non-Gaussianity frontier: What can a measurement with $\sigma$(fNL)$\lesssim$1 tell us about multifield inflation?, Phys. Rev. D 95 (2017), no. 12 123507, [arXiv:1612.05248].
  * [3] P. D. Meerburg et al., Primordial Non-Gaussianity, Bull. Am. Astron. Soc. 51 (2019), no. 3 107, [arXiv:1903.04409].
  * [4] Planck Collaboration, Y. Akrami et al., Planck 2018 results. IX. Constraints on primordial non-Gaussianity, Astron. Astrophys. 641 (2020) A9, [arXiv:1905.05697].
  * [5] N. Dalal, O. Dore, D. Huterer, and A. Shirokov, The imprints of primordial non-gaussianities on large-scale structure: scale dependent bias and abundance of virialized objects, Phys. Rev. D77 (2008) 123514, [arXiv:0710.4560].
  * [6] S. Matarrese and L. Verde, The effect of primordial non-Gaussianity on halo bias, Astrophys. J. Lett. 677 (2008) L77–L80, [arXiv:0801.4826].
  * [7] A. Raccanelli, F. Montanari, D. Bertacca, O. Doré, and R. Durrer, Cosmological Measurements with General Relativistic Galaxy Correlations, JCAP 05 (2016) 009, [arXiv:1505.06179].
  * [8] D. Alonso, P. Bull, P. G. Ferreira, R. Maartens, and M. Santos, Ultra large-scale cosmology in next-generation experiments with single tracers, Astrophys. J. 814 (2015) 145, [arXiv:1505.07596].
  * [9] U. Seljak, Extracting primordial non-gaussianity without cosmic variance, Phys. Rev. Lett. 102 (2009) 021302, [arXiv:0807.1770].
  * [10] P. McDonald and U. Seljak, How to measure redshift-space distortions without sample variance, JCAP 10 (2009) 007, [arXiv:0810.0323].
  * [11] N. Hamaus, U. Seljak, and V. Desjacques, Optimal Constraints on Local Primordial Non-Gaussianity from the Two-Point Statistics of Large-Scale Structure, Phys. Rev. D 84 (2011) 083509, [arXiv:1104.2321].
  * [12] L. R. Abramo and K. E. Leonard, Why multi-tracer surveys beat cosmic variance, Mon. Not. Roy. Astron. Soc. 432 (2013) 318, [arXiv:1302.5444].
  * [13] SKA Collaboration, D. J. Bacon et al., Cosmology with Phase 1 of the Square Kilometre Array: Red Book 2018: Technical specifications and performance forecasts, Publ. Astron. Soc. Austral. 37 (2020) e007, [arXiv:1811.02743].
  * [14] Z. Gomes, S. Camera, M. J. Jarvis, C. Hale, and J. Fonseca, Non-Gaussianity constraints using future radio continuum surveys and the multitracer technique, Mon. Not. Roy. Astron. Soc. 492 (2020), no. 1 1513–1522, [arXiv:1912.08362].
  * [15] S. Jolicoeur, R. Maartens, and S. Dlamini, Constraining primordial non-Gaussianity by combining next-generation galaxy and 21 cm intensity mapping surveys, Eur. Phys. J. C 83 (2023), no. 4 320, [arXiv:2301.02406].
  * [16] D. Alonso and P. G. Ferreira, Constraining ultralarge-scale cosmology with multiple tracers in optical and radio surveys, Phys. Rev. D 92 (2015), no. 6 063525, [arXiv:1507.03550].
  * [17] J. Fonseca, S. Camera, M. Santos, and R. Maartens, Hunting down horizon-scale effects with multi-wavelength surveys, Astrophys. J. 812 (2015), no. 2 L22, [arXiv:1507.04605].
  * [18] M. Ballardini, W. L. Matthewson, and R. Maartens, Constraining primordial non-Gaussianity using two galaxy surveys and CMB lensing, Mon. Not. Roy. Astron. Soc. 489 (2019), no. 2 1950–1956, [arXiv:1906.04730].
  * [19] W. d’Assignies, C. Zhao, J. Yu, and J.-P. Kneib, Cosmological Fisher forecasts for next-generation spectroscopic surveys, Mon. Not. Roy. Astron. Soc. 521 (2023), no. 3 3648, [arXiv:2301.02289].
  * [20] M. B. Squarotti, S. Camera, and R. Maartens, Radio-optical synergies at high redshift to constrain primordial non-Gaussianity, arXiv:2307.00058.
  * [21] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, JCAP 1107 (2011) 034, [arXiv:1104.2933].
  * [22] A. Barreira, G. Cabass, F. Schmidt, A. Pillepich, and D. Nelson, Galaxy bias and primordial non-Gaussianity: insights from galaxy formation simulations with IllustrisTNG, JCAP 12 (2020) 013, [arXiv:2006.09368].
  * [23] A. Barreira, The local PNG bias of neutral Hydrogen, H${I}$, JCAP 04 (2022), no. 04 057, [arXiv:2112.03253].
  * [24] A. Barreira, Can we actually constrain fNL using the scale-dependent bias effect? An illustration of the impact of galaxy bias uncertainties using the BOSS DR12 galaxy power spectrum, JCAP 11 (2022) 013, [arXiv:2205.05673].
  * [25] A. Barreira and E. Krause, Towards optimal and robust f_nl constraints with multi-tracer analyses, JCAP 10 (2023) 044, [arXiv:2302.09066].
  * [26] E. Fondi, L. Verde, F. Villaescusa-Navarro, M. Baldi, W. R. Coulton, G. Jung, D. Karagiannis, M. Liguori, A. Ravenni, and B. D. Wandelt, Taming assembly bias for primordial non-Gaussianity, arXiv:2311.10088.
  * [27] A. G. Adame, S. Avila, V. Gonzalez-Perez, G. Yepes, M. Pellejero, C.-H. Chuang, Y. Feng, J. Garcia-Bellido, A. Knebe, and M. S. Wang, PNG-UNITsims: the response of halo clustering to Primodial Non-Gaussianities as a function of mass, arXiv:2312.12405.
  * [28] Euclid Collaboration, A. Blanchard et al., Euclid preparation: VII. Forecast validation for Euclid cosmological probes, Astron. Astrophys. 642 (2020) A191, [arXiv:1910.09273].
  * [29] L. B. Newburgh et al., HIRAX: A Probe of Dark Energy and Radio Transients, Proc. SPIE Int. Soc. Opt. Eng. 9906 (2016) 99065X, [arXiv:1607.02059].
  * [30] D. Crichton et al., Hydrogen Intensity and Real-Time Analysis Experiment: 256-element array status and overview, J. Astron. Telesc. Instrum. Syst. 8 (2022) 011019, [arXiv:2109.13755].
  * [31] D. J. Schlegel et al., The MegaMapper: A Stage-5 Spectroscopic Instrument Concept for the Study of Inflation and Dark Energy, arXiv:2209.04322.
  * [32] PUMA Collaboration, A. Slosar et al., Packed Ultra-wideband Mapping Array (PUMA): A Radio Telescope for Cosmology and Transients, Bull. Am. Astron. Soc. 51 (2019) 53, [arXiv:1907.12559].
  * [33] A. Merson, A. Smith, A. Benson, Y. Wang, and C. M. Baugh, Linear bias forecasts for emission line cosmological surveys, Mon. Not. Roy. Astron. Soc. 486 (2019), no. 4 5737–5765, [arXiv:1903.02030].
  * [34] N. Sailer, E. Castorina, S. Ferraro, and M. White, Cosmology at high redshift – a probe of fundamental physics, JCAP 12 (2021), no. 12 049, [arXiv:2106.09713].
  * [35] F. Villaescusa-Navarro et al., Ingredients for 21 cm Intensity Mapping, Astrophys. J. 866 (2018), no. 2 135, [arXiv:1804.09180].
  * [36] R. Maartens, J. Fonseca, S. Camera, S. Jolicoeur, J.-A. Viljoen, and C. Clarkson, Magnification and evolution biases in large-scale structure surveys, JCAP 12 (2021), no. 12 009, [arXiv:2107.13401].
  * [37] S. Ferraro et al., Inflation and Dark Energy from Spectroscopy at $z>2$, Bull. Am. Astron. Soc. 51 (2019), no. 3 72, [arXiv:1903.09208].
  * [38] MeerKLASS Collaboration, M. G. Santos et al., MeerKLASS: MeerKAT Large Area Synoptic Survey, in MeerKAT Science: On the Pathway to the SKA, 9, 2017.  arXiv:1709.06099.
  * [39] E. Castorina and F. Villaescusa-Navarro, On the spatial distribution of neutral hydrogen in the Universe: bias and shot-noise of the HI power spectrum, Mon. Not. Roy. Astron. Soc. 471 (2017), no. 2 1788–1796, [arXiv:1609.05157].
  * [40] P. Bull, P. G. Ferreira, P. Patel, and M. G. Santos, Late-time cosmology with 21cm intensity mapping experiments, Astrophys. J. 803 (2015), no. 1 21, [arXiv:1405.1452].
  * [41] D. Alonso, P. G. Ferreira, M. J. Jarvis, and K. Moodley, Calibrating photometric redshifts with intensity mapping observations, Phys. Rev. D 96 (2017), no. 4 043515, [arXiv:1704.01941].
  * [42] Cosmic Visions 21 cm Collaboration, R. Ansari et al., Inflation and Early Dark Energy with a Stage II Hydrogen Intensity Mapping experiment, arXiv:1810.09572.
  * [43] S. Jolicoeur, R. Maartens, E. M. De Weerd, O. Umeh, C. Clarkson, and S. Camera, Detecting the relativistic bispectrum in 21cm intensity maps, JCAP 06 (2021) 039, [arXiv:2009.06197].
  * [44] D. Karagiannis, A. Slosar, and M. Liguori, Forecasts on Primordial non-Gaussianity from 21 cm Intensity Mapping experiments, JCAP 11 (2020) 052, [arXiv:1911.03964].
  * [45] E. Castorina et al., Packed Ultra-wideband Mapping Array (PUMA): Astro2020 RFI Response, arXiv:2002.05072.
  * [46] J.-A. Viljoen, J. Fonseca, and R. Maartens, Constraining the growth rate by combining multiple future surveys, JCAP 09 (2020) 054, [arXiv:2007.04656].
  * [47] S. Casas, I. P. Carucci, V. Pettorino, S. Camera, and M. Martinelli, Constraining gravity with synergies between radio and optical cosmological surveys, Phys. Dark Univ. 39 (2023) 101151, [arXiv:2210.05705].
  * [48] S. Cunnington, S. Camera, and A. Pourtsidou, The degeneracy between primordial non-Gaussianity and foregrounds in 21 cm intensity mapping experiments, Mon. Not. Roy. Astron. Soc. 499 (2020), no. 3 4054–4067, [arXiv:2007.12126].
  * [49] M. Spinelli, I. P. Carucci, S. Cunnington, S. E. Harper, M. O. Irfan, J. Fonseca, A. Pourtsidou, and L. Wolz, SKAO HI intensity mapping: blind foreground subtraction challenge, Mon. Not. Roy. Astron. Soc. 509 (2021), no. 2 2048–2074, [arXiv:2107.10814].
  * [50] S. Cunnington, C. Watkinson, and A. Pourtsidou, The H i intensity mapping bispectrum including observational effects, Mon. Not. Roy. Astron. Soc. 507 (2021), no. 2 1623–1639, [arXiv:2102.11153].
  * [51] J. C. Pober et al., What Next-Generation 21 cm Power Spectrum Measurements Can Teach Us About the Epoch of Reionization, Astrophys. J. 782 (2014) 66, [arXiv:1310.7031].
  * [52] J. C. Pober, The Impact of Foregrounds on Redshift Space Distortion Measurements With the Highly-Redshifted 21 cm Line, Mon. Not. Roy. Astron. Soc. 447 (2015), no. 2 1705–1712, [arXiv:1411.2050].
  * [53] A. Ghosh, F. Mertens, and L. V. E. Koopmans, Deconvolving the wedge: maximum-likelihood power spectra via spherical-wave visibility modelling, Mon. Not. Roy. Astron. Soc. 474 (2018), no. 4 4552–4563, [arXiv:1709.06752].
  * [54] A. Obuljen, E. Castorina, F. Villaescusa-Navarro, and M. Viel, High-redshift post-reionization cosmology with 21cm intensity mapping, JCAP 05 (2018) 004, [arXiv:1709.07893].
  * [55] VIRGO Consortium Collaboration, R. E. Smith, J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchmann, Stable clustering, the halo model and nonlinear cosmological power spectra, Mon. Not. Roy. Astron. Soc. 341 (2003) 1311, [astro-ph/0207664].
  * [56] R. Maartens, S. Jolicoeur, O. Umeh, E. M. De Weerd, C. Clarkson, and S. Camera, Detecting the relativistic galaxy bispectrum, JCAP 03 (2020), no. 03 065, [arXiv:1911.02398].
  * [57] J. Fonseca, J.-A. Viljoen, and R. Maartens, Constraints on the growth rate using the observed galaxy power spectrum, JCAP 1912 (2019), no. 12 028, [arXiv:1907.02975].
  * [58] C. Modi, M. White, A. Slosar, and E. Castorina, Reconstructing large-scale structure with neutral hydrogen surveys, JCAP 11 (2019) 023, [arXiv:1907.02330].
  * [59] M. White, Y.-S. Song, and W. J. Percival, Forecasting Cosmological Constraints from Redshift Surveys, Mon. Not. Roy. Astron. Soc. 397 (2008) 1348–1354, [arXiv:0810.1518].
  * [60] G.-B. Zhao et al., The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: a multitracer analysis in Fourier space for measuring the cosmic structure growth and expansion rate, Mon. Not. Roy. Astron. Soc. 504 (2021), no. 1 33–52, [arXiv:2007.09011].
  * [61] A. Barreira, On the impact of galaxy bias uncertainties on primordial non-Gaussianity constraints, JCAP 12 (2020) 031, [arXiv:2009.06622].
  * [62] D. Karagiannis, R. Maartens, J. Fonseca, S. Camera, and C. Clarkson, Multi-tracer power spectra and bispectra: Formalism, arXiv:2305.04028.
  * [63] D. Karagiannis, A. Lazanu, M. Liguori, A. Raccanelli, N. Bartolo, and L. Verde, Constraining primordial non-Gaussianity with bispectrum and power spectrum from upcoming optical and radio surveys, Mon. Not. Roy. Astron. Soc. 478 (2018), no. 1 1341–1376, [arXiv:1801.09280].
  * [64] V. Yankelevich and C. Porciani, Cosmological information in the redshift-space bispectrum, Mon. Not. Roy. Astron. Soc. 483 (2019), no. 2 2078–2099, [arXiv:1807.07076].
  * [65] Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [66] S. Yahia-Cherif, A. Blanchard, S. Camera, S. Casas, S. Ilić, K. Markovic, A. Pourtsidou, Z. Sakr, D. Sapone, and I. Tutusaus, Validating the Fisher approach for stage IV spectroscopic surveys, Astron. Astrophys. 649 (2021) A52, [arXiv:2007.01812].
  * [67] J.-A. Viljoen, J. Fonseca, and R. Maartens, Multi-wavelength spectroscopic probes: prospects for primordial non-Gaussianity and relativistic effects, JCAP 11 (2021) 010, [arXiv:2107.14057].
  * [68] H.-M. Zhu, U.-L. Pen, Y. Yu, and X. Chen, Recovering lost 21 cm radial modes via cosmic tidal reconstruction, Phys. Rev. D 98 (2018), no. 4 043511, [arXiv:1610.07062].
  * [69] S. Cunnington et al., The foreground transfer function for H i intensity mapping signal reconstruction: MeerKLASS and precision cosmology applications, Mon. Not. Roy. Astron. Soc. 523 (2023), no. 2 2453–2477, [arXiv:2302.07034].
  * [70] J.-A. Viljoen, J. Fonseca, and R. Maartens, Multi-wavelength spectroscopic probes: biases from neglecting light-cone effects, JCAP 12 (2021), no. 12 004, [arXiv:2108.05746].
  * [71] M. Noorikuhani and R. Scoccimarro, Wide-angle and relativistic effects in Fourier-space clustering statistics, Phys. Rev. D 107 (2023), no. 8 083528, [arXiv:2207.12383].
  * [72] P. Paul, C. Clarkson, and R. Maartens, Wide-angle effects in multi-tracer power spectra with Doppler corrections, JCAP 04 (2023) 067, [arXiv:2208.04819].