# Cerebral Palsy Prediction with Frequency Attention Informed Graph
Convolutional Networks

Haozheng Zhang1, Hubert P. H. Shum2 and Edmond S. L. Ho3 *This work was
supported in part by the Royal Society (Ref: IES\R2\181024 and
IES\R1\191147)1Haozheng Zhang is with Durham University, the United Kingdom
<EMAIL_ADDRESS>P. H. Shum is with Durham University, the
United Kingdom. Corresponding author<EMAIL_ADDRESS>S. L. Ho
is with Northumbria University, the United Kingdom<EMAIL_ADDRESS>
###### Abstract

Early diagnosis and intervention are clinically considered the paramount part
of treating cerebral palsy (CP), so it is essential to design an efficient and
interpretable automatic prediction system for CP. We highlight a significant
difference between CP infants’ frequency of human movement and that of the
healthy group, which improves prediction performance. However, the existing
deep learning-based methods did not use the frequency information of infants’
movement for CP prediction. This paper proposes a frequency attention informed
graph convolutional network and validates it on two consumer-grade RGB video
datasets, namely MINI-RGBD and RVI-38 datasets. Our proposed frequency
attention module aids in improving both classification performance and system
interpretability. In addition, we design a frequency-binning method that
retains the critical frequency of the human joint position data while
filtering the noise. Our prediction performance achieves state-of-the-art
research on both datasets. Our work demonstrates the effectiveness of
frequency information in supporting the prediction of CP non-intrusively and
provides a way for supporting the early diagnosis of CP in the resource-
limited regions where the clinical resources are not abundant.

##  I INTRODUCTION

General Movement Assessment (GMA) [1] is being widely used clinically for the
early prediction of cerebral palsy (CP). However, targeted GMA training for
clinicians is a time-consuming and resource-consuming task. As a result, only
a small but increasing number of clinicians have received this training in the
UK and Australia [2]. Furthermore, the process also requires manual inspection
of the infant movement and is prone to subjective assessment. Early studies
applied machine learning techniques (e.g. support vector machine, random
forest) and the optical flow-based video analysis method to propose the
automated GMA systems [3, 4]. But these works still require manual labelling
of infant joint positions. Some later studies focus on the analysis of
frequency domain data. Stahl et al. [5] used an optical flow-based approach to
assess infant movements and then applied wavelet frequency analysis to
evaluate the time-dependent trajectory signals in optical flow data. Rahmati
et al. [6] applied a motion segmentation algorithm to extract motion data from
each limb in the infant video and then classified the infants’ movements with
features obtained by frequency analysis.

Recent deep learning-based systems achieved impressive performance in CP
infants movement prediction. McCay et al. [7] proposed a fully connected deep
learning network and four Convolutional Neural Network (CNN)-based deep
learning architectures to classify the abnormal movements of CP infants by
using the histogram of joint orientation 2D and joint displacement 2D
features, achieved the highest prediction accuracy of $91.67\%$ on the MINI-
RGBD dataset [8]. Zhu [9] further applied the channel attention mechanism on
the 2D-CNN model to interpret the CP prediction outcome on the same dataset.
However, the robustness and generality of their proposed method have not been
fully evaluated since the results are obtained from a single small dataset.

Aiming at the significant difference in joints movement frequency between the
cerebral palsy infants and the healthy group, in this article, we demonstrate
a frequency-based binning mechanism and a graph convolution network to improve
the performance of CP prediction with better interpretability. Firstly, we
employ a pose estimation algorithm, namely Openpose [10] to extract the human
joint position data from the R-GBD video sequences as the input to our system.
Then we propose an automatic frequency-binning module suitable for videos with
different frame rates to reduce data noise and the percentages of high-
frequency movements information in the whole video sequence for CP prediction.
The idea is inspired by both the frequency analysis-based infants CP
prediction methods [6] and our observation. Rahmati et al. [6] provided a
result that comparing with very low or high-frequency ranges, the middle-to-
low frequency range data showed more differences between the healthy group and
the CP group. In addition, we found that the infants’ joint position data in
the high-frequency domain is mainly caused by data noise, such as the
misdetected joint position by Openpose.

We validate our system on the MINI-RGBD dataset [11] and the RVI-38 dataset
[12]. The MINI-RGBD dataset has been widely used for CP classification
performance comparison in the previous work [7, 13, 14, 15], including
synthetic video sequences of 12 normal and CP infants. The RVI-38 is a
recently collected dataset for a more challenging CP prediction task, with a
larger size of data captured during routine clinical care. Experimental
results show that our system achieves state-of-the-art CP prediction
performance on both of the dataset and allows users to interpret the weights
of movement frequencies of different joints in our prediction system.

Our contributions are as follows:

  * •

We interpret the Cerebral Palsy prediction in the joint movement frequency
domain by the attention module. In addition, we designed a new frequency-
binning module that can be applied to both deep learning and machine learning
networks for videos with different frame rates to improve the CP prediction
performance.

  * •

We propose a novel frequency attention informed graph convolutional network
(FAIGCN) for CP prediction from consumer-grade RGB-D videos. Our system
achieves state-of-the-art research on two datasets with strong robustness.

  * •

We open our source code for validation and further development:
https://github.com/zhz95/FAIGCN

##  II Dataset Processing

We verify our models on the Moving INfants In RGB-D synthetic dataset (MINI-
RGBD) [11] and RVI-38 dataset.

###  II-A The MINI-RGBD Dataset

MINI-RGBD was generated by registering and rendering the synthetic Skinned
Multi-Infant (SMIL) model [8] to the RGB-D sequences of real-world moving
infants recorded in the hospital. All 12 RGB-D video sequences were captured
when the infants were half-year-old. The MINI-RGBD dataset is a popular open
resource relating to infants CP as it consists of realistic shape, texture and
movement. It also provides precise ground truth while anonymizing the data by
replacing the raw video frames with computer graphics rendered frames. We
further obtained the annotation of each video sequence shared by [13], which
indicates the presence (i.e. labelled as “normal”) or absence (i.e. labelled
as “abnormal”) of fidgety movements in the video by an independent medical
expert using the GMA method [1].

###  II-B The RVI-38 Dataset

The RVI-38 dataset was collected from a part of routine clinical care at the
Royal Victoria Infirmary (RVI) in Newcastle upon Tyne, UK. There are 38 RGB-D
video sequences of different infants between 12-21 weeks in the RVI-38
dataset. All videos were captured by a consumer-grade handheld camera (Sony
DSC-RX100 with a resolution of 1980x1080 and the 25FPS frame rate). The length
of videos ranges between 40 seconds and 5 minutes, with an average length of 3
minutes and 36 seconds. The camera was set above the baby, and the infant’s
movement was photographed from top to bottom. All videos were annotated using
the GMA method by two experienced assessors. The annotations indicate the
presence (i.e. labelled as “normal”) or absence (i.e. labelled as “abnormal”)
of fidgety movements in the video.

###  II-C Data Preprocessing

For more effective CP predictions, we extract 2D skeleton features from the
video sequences. We apply OpenPose [10] for pose estimation due to its high
accuracy in detecting the posture of the infants, and it is less sensitive to
variations on the appearance. OpenPose returns the 2D coordinates $(x,y)$ for
18 human joint landmarks and a confidence score $C$ for each joint estimation.
However, for joints that are self-occluded or without clear visual features,
OpenPose would not be able to deduce their position, and zero values would be
returned as the joint positions. As this may impact the performance of the
prediction system, we propose to preprocess the data by replacing the zero
values in frame $f$ with the linear interpolation of neighbouring non-zero
frames.

In order to overcome the overfitting in small size dataset, we implement
several processes. 1) We calculate the global normalization of joint positions
frame by frame to reduces the infant’s global translation. To achieve this, we
set the center of the triangle of the neck and two hip joints as the global
origin, then relocate each joints by the relative distance between joints. 2)
To normalize the x-direction and y-direction pose features, we align the line
between the global orgin and the neck joint with the y-axis and keep the neck
joint above the global origin.

Figure 1: The overview of our proposed framework. Part I is the overall
network architecture, Part II is the design of each FAIGCN layer. Figure 2: An
example frame of 18-joints Openpose [10] posture layout for infant kicking in
the MINI-RGBD dataset [8]. The size of the light point represents the size of
the attention value, that is, the importance of the movement frequency of the
joint to our network in CP prediction at that frame.

##  III Frequency Attention Informed GCN System

The proposed system consists of two parts (seen in Fig.1): (1) The frequency-
binning module transforms the input joint movement features into the frequency
domain, then filters high-frequency information to make our prediction network
focus on low-to-mid infants movement frequencies. (2) The proposed Frequency
Attention Informed GCN for CP prediction and interpretation. Fig. 2 shows an
example of attention visualisation.

###  III-A The Frequency-binning Module

Given the infants’ joint movement sequences as input, we propose using
frequency operations on joint position data for CP prediction. It is motivated
by two observations. Firstly, the body movement frequencies of healthy infants
are different from infants who suffered from CP [6], and the low-frequency
range information of body movement is more critical for fidgety movements
(FMs). FMs are moderate speed movement of the neck, trunk and limbs with
different accelerations in various directions [16]. Previous work [16, 17] has
shown that the absence of FMs is an essential distinguishing feature of CP
infants from healthy infants. Frequency-binning can filter the high-frequency
(e.g., above 6 HZ) information of joint position data after FFT, thus making
the classification network focused on low-to-mid (e.g., 0-5 HZ) frequency
infant movements without eliminating raw data. Secondly, the infant movements
frequencies are generally low, and the high-frequency range from the joint
position data is mainly due to data noise, such as the misdetected joint
position and the video capture error from the datasets.

As a solution, we design a frequency-binning module that retains the critical
frequency of the joint position data while filtering the noise. The module
employs Fast Fourier Transform (FFT) to convert the time series of joint
positions into the frequency domain, then applies frequency binning to obtain
the motion frequency information mainly distributed in the low-to-mid band.
This module is adaptable for videos with a frame rate between 24 FPS to 60 FPS
and is suitable for both DNN or machine learning-based classification models.
The core of the module is the binning strategy, in which we design a formula
to use finer bins for the more crucial low-to-mid frequency and coaster bins
for higher frequency.

####  III-A1 Fast Fourier Transform (FFT)

We apply Bluestein’s FFT algorithm [18], a discrete Fourier transform
algorithm, on all 2D joints movements time series to transform original joints
position features into the frequency domain and obtain the frequency
components:

$X_{k}=\sum_{n=0}^{N-1}x_{n}e^{\frac{-i2\pi kn}{N}},\quad k=0,\dots,N-1$ (1)

where $x_{n}$ is a time series, $e^{\frac{i2\pi}{N}}$ is a primitive $N^{th}$
root of 1.

####  III-A2 The Binning Strategy

We propose a data binning strategy to emphasize the importance of low-
frequency information of the joint position data. Under the strategy, the
width of the bins are different - smaller width bins are used for low-
frequency range and increasingly larger-width bins for higher frequency range:

$b_{n}=\begin{cases}$Round$(b_{0}\cdot c^{n}),&\text{if $b_{n}\cdot
c^{n}<3$,}\\\ $Ceiling$(b_{0}\cdot c^{n}),&\text{if $b_{n}\cdot c^{n}\geq
3$,}\end{cases}$ (2)

where $b_{n}$ is the width of the $n^{th}$ bin, $b_{0}=1$, and $c$ is a
controllable parameter. Note that the width of each bin needs to be an integer
as FFT is a discrete (i.e. integer-based) system. This equation takes the
round of the bin width when the width is less than three units to increase the
density of the bins in low-to-mid frequency. According to the characteristics
of rapid exponential growth, this function distinguishes the density of the
middle frequency band and the high-frequency band for bins with a width
greater than two units by rounding up the value greater than three units.
Empirically, as shown in Sec IV, we achieve the best prediction accuracy when
$c=1.00264$ for the 25 FPS videos. The parameter $c$ could be automatically
generated for the best binning results by achieving the highest CP prediction
performance for different datasets.

As a result, after being processed by the frequency-binning module, input
joint position data are transformed into the frequency domain and endowed with
an important characteristic: low-to-mid frequency information occupies a
significantly more prominent emphasis.

###  III-B The CP Prediction Network

As shown in Fig.1 , we propose a Frequency Attention Informed Graph
Convolutional Network (FAIGCN) for CP prediction by classifying low-to-mid
frequency band infant movement frequency features with the attention
mechanism.

####  III-B1 Frequency Attention Informed Network

Most of the previous DNN-based studies on infants CP prediction are based on
traditional Convolutional Neural Networks (CNN). However, traditional discrete
convolution from CNN can only maintain translational invariance on Euclidean
data, which is not suitable for graph structure data such as the human
skeletal graph generated from OpenPose [10].

Therefore, we employ a GCN [19] to learn the infant joints dependencies from
the pose graph. Inspired by [20], we apply the pose graph which align with the
human skeletal graph $G=(V,E)$ for interpreting which joint’s movement
frequency features are considered to be important in CP prediction task. In
this graph, $\\{V={v_{bi}|b=1,...,B;i=1,...,N}\\}$ denotes the frequencies of
all joints, where $v_{bi}$ represents the b-th frequency bin of $i$-th joint.
The edge set E includes: (1) the intra-skeleton connection at each frequency
band, $\\{v_{bi}v_{bj}|(i,j)\in K\\}$, where $K$ is designed by the natural
connections of human joints. (2) the inter-frequency edges which connect the
frequency bins of a joint in the low-to-high frequency order,
$\\{v_{bi}v_{(b+1)i}\\}$.

The graph convolutional operation of FAIGCN is followed by [19], where the
propagation rule between layers can be represented by Eq. 3.

$\mathbf{H}^{(l+1)}=\sigma\left(\tilde{\mathbf{D}}^{-\frac{1}{2}}\tilde{A}\tilde{\mathbf{D}}^{-\frac{1}{2}}\mathbf{H}^{(l)}\mathbf{W}^{(l)}\right)$
(3)

where $\tilde{\mathbf{A}}=\mathbf{A}+\mathbf{I}_{L}$ is known as the adjacency
matrix of an undirected graph. $\mathbf{I}_{L}$ is an L dimensions identity
matrix. $\tilde{\mathbf{D}}_{ii}=\sum_{j}\tilde{\mathbf{A}}_{ij}$ and
$\mathbf{W}^{(l)}$ is a learnable weight matrix specified to the layer. The
nonlinear activation function $\sigma(\cdot)$ is set as ReLU in our network.

We propose the following frequency attention-informed mechanism to learn the
weight of frequency features. We aggregate the frequency features obtained
from the frequency-binning module
$\\{\mathbf{h}_{1,i},\mathbf{h}_{2,i},\dots,\mathbf{h}_{B,i}\\}$ with
attentions $\alpha_{b,i}$ by Eq. 4

$\mathbf{v}_{k}=\sum_{b=1}^{B}\alpha_{b,i}\mathbf{h}_{b,i}$ (4)

in which the frequency attention weight $\alpha_{b,i}$ is defined as:

$\alpha_{b,i}=\frac{\text{exp}\left(\sigma^{{}^{\prime}}_{n}\left(\mathbf{w}_{\alpha}^{\mathsf{T}},\mathbf{z}_{b,i}\right)\right)}{\sum_{b}\text{exp}\left(\sigma^{{}^{\prime}}\left(\mathbf{w}_{\alpha}^{\mathsf{T}},\mathbf{z}_{b,i}\right)\right)}$
(5)

$\mathbf{z}_{b,i}=tanh\left(\mathbf{W}_{z}\mathbf{h}_{b,i}\right)$ (6)

where $\sigma^{{}^{\prime}}_{n}$ is an adjustable activation function as
follows:

$\sigma^{{}^{\prime}}_{n}=\
\begin{cases}1+\left(\frac{\mathbf{w}_{\alpha}}{\left\|\mathbf{w}_{\alpha}\right\|}\right)^{\mathsf{T}}\left(\frac{\mathbf{z}_{\alpha}}{\left\|\mathbf{z}_{\alpha}\right\|}\right)&,n=1\\\
\mathbf{w}_{\alpha}^{\mathsf{T}}\mathbf{z}_{\alpha}&,n=2\end{cases}$ (7)

where $\mathbf{w}_{\alpha}$ and $\mathbf{W}_{z}$ are learnable parameters.

####  III-B2 Network Adaptation

As can be seen from Fig.1, the input layer transforms the tensor format of
input frequency data to fit in the network. Then, we use two FAIGCN layers
with 32, 64 output channels respectively. Each FAIGCN, in turn, consists of a
GCN layer, a batch normalization layer, a dropout and a ReLU layer. The kernel
sizes of FAIGCN layers $K=3,3$, and $stride=1,2$, respectively. We put a
global pooling layer after two FAIGCN layers. We applied the average pooling
as it provides the highest robustness. At the last, we put a fully connected
layer to classify features for CP prediction. The optimizer is chosen as
Adams, and we train the model with batch size $=1$, learning rate $=0.0001$
with $0.1$ decay every 100 epoches, Max Epoch $=500$ on the MINI-RGBD dataset;
batch size $=4$, learning rate $=0.001$ with $0.1$ decay every 100 epoches,
Max Epoch $=500$ on the RVI-38 dataset.

##  IV Experiments

Our experiments were run on a PC with Ubuntu 18.04 and an NVIDIA GeForce RTX
3080. The total model training time on MINI-RGBD with 12 sequences is about an
hour, including estimation of the joints position from RGB videos. But it only
takes about 50s for the CP prediction of 1000 frames ($\sim 33$s) video
sequence, which can be employed in interactive-time diagnosis.

###  IV-A Experimental Settings

In this paper, we conduct the leave-one-out cross-validation among two
datasets to evaluate our proposed system. This setting utilises all data and
ensures that the prediction system is evaluated against unseen data. Our
evaluation metrics are introduced in Sec IV-C. We report the best result for
each method to be consistent with several related works in literature [5, 17,
6, 15, 12].

###  IV-B Comparing with State-of-the-art Methods

In order to evaluate the effectiveness of our system, we compare FAIGCN with
the following methods:

  * •

FCNet [7]: This method uses fully connected deep network architectures to the
Histogram of Joint Displacement 2D (HOJD2D) and Histogram of Joint Orientation
2D (HOJO2D) calculated from human joint positions. For the HOJD2D feature, the
displacements of each joint are extracted every five frames and segmented into
16 bins. The feature of HOJO2D represents the joint orientation in 2D space,
and the feature is also segmented into 16 bins.

  * •

Conv1D-1, Conv1D-2 [7]: They are two 1D convolutional neural networks, each of
them consists of two 1D convolutional layers with differences in the output
channel sizes. They are proposed to classify the abnormal infant movements by
feature HOJO2D or HOJD2D (HOJO/D2D).

  * •

Conv2D-1, Conv2D-2 [7]: They are two 2D convolutional neural networks, each of
them consists of two 2D convolutional layers with differences in the output
channel sizes.

  * •

CANet [9]: This method proposes a 2D convolutional neural network with the
squeeze-and-excitation channel attention module. The whole system is proposed
for the CP classification task on the MINI-RGBD dataset.

  * •

ST-GCN (Spatial Temporal Graph Convolutional Network) [20]: This is a graph
convolutional neural network for human skeleton data (e.g. joint position).

  * •

STAM [21]: This is a spatial-temporal graph convolutional neural network with
the attention mechanism.

  * •

Ens-1 [13]: This method uses an ensemble classifier on the fused feature of
HOJO2D + HOJD2D (HOJO+D2D) with eight bins.

  * •

Ens-2 [12]: This method extends [13] by fusing four pose-related features and
three velocity-related features.

  * •

Ens-3 [22]: This method extends [13] by extracting features at limb-level from
small video segments to locate abnormal movements spatiotemporally.

  * •

MCI [14]: This method uses a threshold model to classify the infant CP via
Movement Complexity Index (MCI), where MCI is computed by extracting the
infant’s limb angle features.

###  IV-C Evaluation Metrics

We evaluate our system and other state-of-the-art methods by following five
metrics: the prediction accuracy (AC) shows the percentage of correctly
predicted individuals in the dataset; the sensitivity (SE) shows the
percentage of correctly predicted positive individuals among the total number
of positive individuals in the dataset; the specificity (SP) shows the
percentage of correctly predicted negative individuals among the total number
of negative individuals in the dataset; F1-Score evaluates the binary
classification performance by calculating the harmonic mean of the precision
and recall; Matthews Correlation Coefficient (MCC) [23] provides a reliable
performance metric for imbalanced dataset [24].

Table I: The comparison with state-of-the-arts on the MINI-RGBD

Method | Feature | AC | SE | SP | F1 | MCC
---|---|---|---|---|---|---
FCNet [7] | HOJD2D | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
FCNet [7] | HOJO2D | 83.33 | 75.00 | 87.50 | 75.00 | 62.50
Conv1D [7] | HOJO/D2D | 83.33 | 75.00 | 87.50 | 75.00 | 62.50
Conv2D [7] | HOJO/D2D | 83.33 | 75.00 | 87.50 | 75.00 | 62.50
Conv2D [7] | HOJO+D2D | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
Conv2D [7] | Pose | 83.33 | 75.00 | 87.50 | 75.00 | 62.50
CANet [9] | Pose | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
ST-GCN[20] | Pose | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
STAM [21] | Pose | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
Ens-1 [13] | HOJO+D2D | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
Ens-2 [12] | Velocity | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
Ens-2 [12] |  Pose* | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
Ens-2 [12] |  Vel.+Pose* | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
Ens-3 [22] | HOJO+D2D | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
MCI [14] | Limb angle | 91.67 | 100.00 | 87.50 | 88.989 | 83.67
FAIGCN | Motion Freq. | 100.00 | 100.00 | 100.00 | 100.00 | 100.00

  * *

The Pose and velocity features here fuses several hand-crafted
features including HOJOD2D

###  IV-D Result Comparison

We report the prediction results on MINI-RGBD and RVI-38 datasets on Table I
and Table II respectively. From our evaluation, we propose the following
observations: (1) Our FAIGCN system outperforms the state-of-the-art DNN based
methods in both datasets. Comparing with other non-DNN based methods, our
system also achieves state-of-the-art performance in two datasets. (2) We can
observe the advantage of the attention mechanism in DNN as CANet outperforms
Conv2D-Pose from Table I and Table II, and CANet outperforms Conv1D/Conv2D,
STAM outperforms ST-GCN from Table II. (3) From Table II, it can be seen that
ST-GCN, STAM and FAIGCN outperform all CNN-based methods (i.e. Conv1D, Conv2D
and CANet), which confirms the advantage of using graph structure to analyse
human pose data. (4) We notice that the methods that use early fusion on
features outperform those using only a single kind of feature. It can be seen
by comparing Conv2D-HOJO/D2D with Conv2D-HOJO+D2D, or comparing
Ens-2-Velocity/Pose with Ens-2-Velocity+Pose in both tables. Therefore, we
consider fusing our movement frequency features with other features in future
work. (5) An interesting finding is that the Machine learning-based methods
such as Ens-1, Ens-2 and Ens-3 outperform DNN-based methods except for FAIGCN.
On the one hand, it shows the superiority of hand-craft features in the
classification tasks; On the other hand, it inspires us to compare FAIGCN with
machine learning-based methods with the same features, seen in the Sec. IV-E
below.

Table II: The comparison with state-of-the-arts on the RVI-38

Method | Feature | AC | SE | SP | F1 | MCC
---|---|---|---|---|---|---
Conv2D [7] | Pose | 81.58 | 33.33 | 90.63 | 36.36 | 25.85
CANet [9] | Pose | 86.84 | 66.67 | 90.63 | 61.54 | 53.89
ST-GCN [20] | Pose | 89.47 | 66.67 | 93.75 | 62.50 | 60.42
STAM [21] | Pose | 92.11 | 83.33 | 93.75 | 76.92 | 72.51
Ens-1 [13] | HOJO+D2D | 94.74 | 83.33 | 96.88 | 83.33 | 80.21
Ens-2 [12] | Velocity | 94.74 | 83.33 | 96.88 | 83.33 | 80.21
Ens-2 [12] |  Pose* | 94.74 | 83.33 | 96.88 | 83.33 | 80.21
Ens-2 [12] |  Vel.+Pose* | 97.37 | 83.33 | 100.00 | 90.91 | 89.89
FAIGCN | Motion Freq. | 97.37 | 83.33 | 100.00 | 90.91 | 89.89

  * *

The Pose and velocity features here fuses several hand-crafted
features including HOJOD2D

###  IV-E Ablation study

The ablation study is divided into two parts:

In the first part, we evaluate whether there is any adverse effect on
prediction performance caused by the frequency-binning module (B.) or the
attention module (A.). The corresponding result is displayed in Table III. Due
to the limitation of the MINI-RGND dataset’s size, the contribution of the
attention mechanism is not reflected significantly. But from the larger RVI-38
dataset, it can be seen that there is a significant improvement by applying
the attention module or frequency binning module.

In the second part, we implement four machine learning-based methods with a
proposed frequency-binning module to predict the CP using the movement
frequency features from both datasets. The methods are Support Vector Machine
(SVM), Decision Tree (Tree), Logistic Regression (LR) and Linear Discriminant
Analysis (LDA). The ensemble of classification models Ens-1 [13], Ens-2 [12]
and Ens-3 [22] are not included since the types of the ensemble classifier in
Matlab was used which consists of a wide range of classifiers and handles the
late-fusion internally. Besides, we validate the effectiveness and robustness
of the frequency-binning module by eliminating it from each method. The
results are reported in Table IV. We observe the advantage of using the
proposed frequency-binning module as each method outperforms its variant of no
frequency-binning module, except in the case of SVM in the MINI-RGBD dataset.
In addition, we note that our system outperforms the implemented machine
learning-based methods. It demonstrates the effectiveness of graph
convolutional neural network in dealing with the same features.

Table III: The performance of FAIGCN and its simplified variants The MINI-RGBD
dataset
---
Method | AC | SE | SP | F1 | MCC
FAIGCN-full | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
w/o A. | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
w/o B. | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
w/o A. B. | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
|  |  |  |  |
The RVI-38 dataset
Method | AC | SE | SP | F1 | MCC
FAIGCN-full | 97.37 | 83.33 | 100.00 | 90.91 | 89.89
w/o A. | 92.11 | 83.33 | 93.75 | 76.92 | 72.51
w/o B. | 89.47 | 66.67 | 93.75 | 62.50 | 60.42
w/o A. B. | 86.84 | 66.67 | 90.63 | 61.54 | 53.89
Table IV: The comparison with machine learning based methods and their variant
without frequency-binning module The MINI-RGBD dataset
---
Methods | AC | SE | SP | F1 | MCC
SVM | 66.77 | 75.00 | 62.50 | 66.67 | 35.36
SVM w/o B. | 66.77 | 75.00 | 62.50 | 66.67 | 35.36
Tree | 75.00 | 75.00 | 75.00 | 66.67 | 47.81
Tree w/o B. | 66.77 | 75.00 | 62.50 | 66.67 | 35.36
LDA | 83.33 | 75.00 | 87.50 | 75.00 | 62.50
LDA w/o B. | 75.00 | 75.00 | 75.00 | 66.67 | 47.81
LR | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
LR w/o B. | 75.00 | 75.00 | 75.00 | 66.67 | 47.81
FAIGCN-full | 100.00 | 100.00 | 100.00 | 100.00 | 100.00
FAIGCN w/o B. | 91.67 | 100.00 | 87.50 | 88.89 | 83.67
|  |  |  |  |
The RVI-38 dataset
Methods | AC | SE | SP | F1 | MCC
SVM | 63.16 | 66.67 | 62.50 | 66.67 | 35.36
SVM w/o B. | 55.26 | 50.00 | 56.25 | 26.09 | 4.58
Tree | 81.58 | 66.67 | 84.38 | 53.33 | 43.78
Tree w/o B. | 68.42 | 50.00 | 71.89 | 33.33 | 17.16
LDA | 83.33 | 75.00 | 87.50 | 75.00 | 62.50
LDA w/o B. | 65.79 | 66.67 | 65.63 | 38.10 | 24.09
LR | 78.95 | 66.67 | 81.25 | 50.00 | 39.68
LR w/o B. | 57.89 | 66.67 | 56.25 | 33.33 | 16.74
FAIGCN-full | 97.37 | 83.33 | 100.00 | 90.91 | 89.89
FAIGCN w/o B. | 89.47 | 66.67 | 93.75 | 62.50 | 60.42
Figure 3: The robustness test compared with the state-of-the-art DNN methods.
The short vertical bar of each method in different noise-level denotes the
accuracy range between the first quartile and third quartile among all cross-
validations. The line between each bar is linked by the mean accuracy value.

###  IV-F Robustness Test

In order to evaluate the robustness of our system and other state-of-the-art
DNN-based methods [9, 20, 21], we simulate different datasets by adding
different levels of Gaussian noise to the infant joint pose data. The noise
level is divided into four levels from 15% standard deviation to 120% standard
deviation of each infant’s joint pose data. The tests results are displayed in
Fig. 3. The accuracy in the y-axis is the average accuracy among ten leave-
one-out cross-validations with ten different training seeds.

From Fig. 3, we observe that as the noise level increases, each method
decreases more slowly on the RVI-38 dataset compared to the MINI-RGBD dataset,
reflecting the stronger robustness brought by training the model on a larger
dataset. In addition, we are seeing that the accuracy of our system shows a
slower decreasing trend under different noise levels, which represents the
stronger stability and robustness of our system.

###  IV-G Qualitative Analysis

Figure 4: The visualization of attention weights of different joints among all
cross-validations on each dataset.

Fig. 4 visualises the interpretability of proposed attention module by
presenting the attention value of each joint among all leave-one-out cross-
validations on each dataset. We observe that the attention value of ‘Right
Knee’,‘Left Knee’,‘Right Wrist’ and ‘Left Wrist’ is significantly higher than
other joints on both datasets. It indicates our system pays more attention to
the movement frequencies of infants’ knees and wrists, which is convincing
since the movements of those joints have the most significant frequency change
in the video recordings. In addition, the frequency range of ‘Right Eye’,
‘Left Eye’, ‘Right Ear’ and ‘Left Ear’ is lower than other joints
significantly. One possible reason is that the self-occlusion (e.g. infant
turns head) brings the noise to Openpose estimation of these joints, so that
the attention module of our system lowers the weights to filter the noisy
data. Besides, we notice that the attention weight range of most of the joints
in the RVI-38 dataset is larger than those in the MINI-RGBD dataset. It could
be caused by more information from the larger dataset.

##  V Conclusion

In this work, we propose a novel interpretable frequency attention informed
graph convolutional network to predict cerebral palsy infants. We design a
binning module for CP data to increase the weight of the low-to-mid frequency
data to improve the CP prediction performance, which is adaptable for the
videos with a frame rate between 24 FPS to 60 FPS and suitable for both DNN or
machine learning-based classification model. Furthermore, we propose a
frequency attention module to further improve the prediction performance and
visualise the important joints that the network considers in CP prediction.
Experimental results show the importance of low-to-mid frequency data and the
effectiveness and robustness of our system in supporting the prediction of CP
non-intrusively, and provides a way for supporting the early diagnosis of CP
in the resource-limited regions where the clinical resources are not abundant.
Future work is scheduled to apply our methods to a more extensive clinical
dataset. In addition, an interesting future work could be using the
transformer to interpret the frequency features better.

##  VI Compliance with Ethical Standards

The collection of the RVI-38 dataset has been ethically approved by the host
organisation (Ref: 9865), the Research Ethics Committee (REC), the Health
Research Authority (HRA), and Health and Care Research Wales (HCRW) (Ref:
19/LO/0606, IRAS project ID: 252317). The MINI-RGBD dataset used in this study
is made open access by Fraunhofer IOSB [8], which had ethical approval.

## References

  * [1] C. Einspieler and H. F. R. Prechtl, “Prechtl’s assessment of general movements: a diagnostic tool for the functional assessment of the young nervous system,” _Mental retardation and developmental disabilities research reviews_ , vol. 11, no. 1, p. 61—67, 2005.
  * [2] D. Graham, S. P. Paget, and N. Wimalasundera., “Current thinking in the health care management of children with cerebral palsy,” _Medical Journal of Australia_ , vol. 210, no. 3, pp. 129–135, 2019.
  * [3] S. Orlandi, K. Raghuram, and C. R. S. et al., “Detection of atypical and typical infant movements using computer-based video analysis,” in _IEEE EMBC_ , 2018, pp. 3598–3601.
  * [4] E. A. F. Ihlen, R. Støen, and L. Boswell.et al., “Machine learning of infant spontaneous movements for the early prediction of cerebral palsy: A multi-site cohort study,” _Journal of Clinical Medicine_ , vol. 9, no. 1, 2020.
  * [5] A. Stahl, C. Schellewald, and O. S. et al., “An optical flow-based method to predict infantile cerebral palsy,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 20, no. 4, pp. 605–614, 2012\.
  * [6] H. Rahmati, H. Martens, O. M. Aamo, O. Stavdahl, R. Stoen, and L. Adde, “Frequency analysis and feature reduction method for prediction of cerebral palsy in young infants.” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 24, no. 11, pp. 1225–1234, 2016.
  * [7] K. D. McCay, E. S. L. Ho, H. P. H. Shum, G. Fehringer, C. Marcroft, and N. D. Embleton, “Abnormal infant movements classification with deep learning on pose-based features,” _IEEE Access_ , vol. 8, pp. 51 582–51 592, 2020\.
  * [8] N. Hesse, S. Pujades, and J. t. Romero, “Learning an infant body model from RGB-D data for accurate full body motion analysis,” in _Int. Conf. on Medical Image Computing and Computer Assisted Intervention (MICCAI)_ , 2018.
  * [9] M. Zhu, Q. Men, E. S. L. Ho, H. Leung, and H. P. H. Shum, “Interpreting deep learning based cerebral palsy prediction with channel attention,” in _IEEE EMBS International Conference on Biomedical and Health Informatics (BHI)_ , 2021, pp. 1–4.
  * [10] Z. Cao, G. Hidalgo, T. Simon, S.-E. Wei, and Y. Sheikh, “Openpose: Realtime multi-person 2d pose estimation using part affinity fields,” _arXiv e-prints_ , p. arXiv:1812.08008, Dec. 2018.
  * [11] N. Hesse, C. Bodensteiner, M. Arens, U. G. Hofmann, R. Weinberger, and A. S. Schroeder, “Computer vision for medical infant motion analysis: State of the art and RGB-D data set,” in _ECCV 2018 Workshops_ , 2018.
  * [12] K. D. McCay, P. Hu, H. P. H. Shum, W. L. Woo, C. Marcroft, N. D. Embleton, A. Munteanu, and E. S. L. Ho, “A pose-based feature fusion and classification framework for the early prediction of cerebral palsy in infants,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , pp. 1–1, 2021.
  * [13] K. D. McCay, E. S. L. Ho, C. Marcroft, and N. D. Embleton, “Establishing pose based features using histograms for the detection of abnormal infant movements,” in _IEEE EMBC_ , 2019, pp. 5469–5472.
  * [14] Q. Wu, G. Xu, F. Wei, L. Chen, and S. Zhang, “RGB-D videos-based early prediction of infant cerebral palsy via general movements complexity,” _IEEE Access_ , vol. 9, pp. 42 314–42 324, 2021.
  * [15] D. Sakkos, K. D. Mccay, C. Marcroft, N. D. Embleton, S. Chattopadhyay, and E. S. L. Ho, “Identification of abnormal movements in infants: A deep neural network for body part-based prediction of cerebral palsy,” _IEEE Access_ , vol. 9, pp. 94 281–94 292, 2021.
  * [16] F. Ferrari, G. Cioni, and C. t. Einspieler, “Cramped synchronized general movements in preterm infants as an early marker for cerebral palsy,” _Archives of Pediatrics and Adolescent Medicine_ , vol. 156, no. 5, pp. 460–467, 2002.
  * [17] C. Einspieler, R. Peharz, and P. B. Marschik, “Fidgety movements – tiny in appearance, but huge in impact,” _Jornal de Pediatria_ , vol. 92, no. 3, Supplement 1, pp. S64–S70, 2016.
  * [18] L. Bluestein, “A linear filtering approach to the computation of discrete fourier transform,” _IEEE Transactions on Audio and Electroacoustics_ , vol. 18, no. 4, pp. 451–455, 1970.
  * [19] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” _ICLR_ , 2017.
  * [20] S. Yan, Y. Xiong, and D. Lin, “Spatial temporal graph convolutional networks for skeleton-based action recognition,” _AAAI Conference on Artificial Intelligence_ , 2018.
  * [21] B. Nguyen-Thai, V. Le, C. Morgan, N. Badawi, T. Tran, and S. Venkatesh, “A spatio-temporal attention-based model for infant movement assessment from videos,” _IEEE Journal of Biomedical and Health Informatics_ , pp. 1–1, 2021\.
  * [22] K. D. McCay, E. S. L. Ho, D. Sakkos, W. L. Woo, C. Marcroft, P. Dulson, and N. D. Embleton, “Towards explainable abnormal infant movements identification: A body-part based prediction and visualisation framework,” in _2021 IEEE EMBS International Conference on Biomedical and Health Informatics (BHI)_ , 2021, pp. 1–4.
  * [23] B. Matthews, “Comparison of the predicted and observed secondary structure of t4 phage lysozyme,” _Biochimica et Biophysica Acta (BBA) - Protein Structure_ , vol. 405, no. 2, pp. 442–451, 1975.
  * [24] D. Chicco and G. Jurman, “The advantages of the matthews correlation coefficient (mcc) over f1 score and accuracy in binary classification evaluation,” _BMC Genomics_ , vol. 21, no. 1, pp. 1–13, 2020.


Incompatibility of Standard Galaxy Bias Models in General Relativity

Jaiyul Yooa,b

a Center for Theoretical Astrophysics and Cosmology, Institute for
Computational Science
University of Zürich, Winterthurerstrasse 190, CH-8057, Zürich, Switzerland

bPhysics Institute, University of Zürich, Winterthurerstrasse 190, CH-8057,
Zürich, Switzerland

Abstract
The standard model for galaxy bias is built in a Newtonian framework, and
several attempts have been made in the past to put it in a relativistic
framework. The focus of past works was, however, to use the same Newtonian
formulation, but to provide its interpretation in a relativistic framework by
either fixing a gauge condition or transforming to a local coordinate system.
Here we demonstrate that these reverse-engineered approaches do not respect
the diffeomorphism symmetry in general relativity, and we need to develop a
covariant model of galaxy bias that is diffeomorphism compatible. We consider
a simple toy model for galaxy bias and discuss the impact for measuring the
primordial non-Gaussianity.

11footnotetext<EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
  2. 2 Proper-time hypersurface for modeling galaxy bias
    1. 2.1 Gauge choice for galaxy bias: temporal and spatial
    2. 2.2 Second-order bias and primordial non-Gaussianity
    3. 2.3 Spatial gauge choices for galaxy bias in literature
  3. 3 Gauge-invariance of cosmological observables
    1. 3.1 Cosmological observables
    2. 3.2 A covariant toy model for galaxy bias
  4. 4 Conclusions

##  1 Introduction

Luminous galaxies provide a powerful way to map the large-scale structure of
the Universe, and the relation between the galaxy and the matter distributions
is known as the galaxy bias ([1, 2, 3, 4, 5], see also [6] for a recent
review). In the simplest model [1, 2, 3, 4, 5], the fluctuation
$\delta_{g}:=\delta n_{g}/\bar{n}_{g}$ in the galaxy number density $n_{g}$ is
related to the matter density fluctuation
$\delta_{m}:=\delta\rho_{m}/\bar{\rho}_{m}$ as

$\delta_{g}=b~{}\delta_{m}~{},$ (1.1)

where the galaxy bias factor $b$ is a constant. The linear bias relation in
eq. (1.1) has been widely used to extract cosmological information on large
scales (e.g., [7, 8, 9, 10]), and extensive efforts in the community have been
devoted to building a more precise model for galaxy bias on nonlinear scales
[11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30]. So far, most of the approaches to modeling galaxy bias are built in a
Newtonian framework. Since measurement precision in large-scale surveys is
highest on small scales, where the Newtonian dynamics is accurate, it is
natural and practical that most resources in terms of observational and
theoretical investigations are devoted in this direction.

However, viewed from a relativistic perspective, the simple linear bias
relation in eq. (1.1) is ambiguous, because the matter density fluctuation
gauge-transforms as

$\tilde{\delta}_{m}=\delta_{m}+3\mathcal{H}~{}T+\mathcal{O}(2)~{},$ (1.2)

under a general coordinate transformation
$\tilde{x}^{\mu}=x^{\mu}+\zeta^{\mu}$, where the infinitesimal transformation
is $\zeta^{\mu}=:(T,\mathcal{L}^{i})$, the conformal Hubble parameter is
$\mathcal{H}=a^{\prime}/a$, and the prime denotes the derivative with respect
to the conformal time $\eta$. Similarly, the galaxy number density fluctuation
also gauge-transforms as

$\tilde{\delta}_{g}=\delta_{g}-{\bar{n}_{g}^{\prime}\over\bar{n}_{g}}~{}T+\mathcal{O}(2)~{}.$
(1.3)

Therefore, if we want to maintain the linear bias relation in eq. (1.1) in all
coordinates, the evolution of the background galaxy number density
$\bar{n}_{g}$ should be related to the bias parameter as

${d\ln\bar{n}_{g}\over d\eta}=-3\mathcal{H}~{}b~{},\qquad\qquad{\rm
or}\qquad\qquad{d\ln\bar{n}_{g}\over d\ln(1+z)}=3~{}b~{},$ (1.4)

which we already know is not true in reality (see, e.g., [31, 32] for recent
observations). In contrast, note that eq. (1.4) is perfectly valid for the
case of the matter density as a model for galaxy, i.e., $n_{g}\equiv\rho_{m}$
with $b=1$. Furthermore, the relation $\rho_{m}=\bar{\rho}_{m}(1+\delta_{m})$
is valid in all coordinates or covariant. This illustrates the issues that the
standard model for galaxy bias is incompatible with the diffeomorphism
symmetry in general relativity.

Since galaxy clustering is essentially the dynamics of non-relativistic matter
components, the contributions of the relativistic effects in galaxy clustering
are naturally small compared to the matter density fluctuations, unless the
scales of interest are close to the horizon scales [33, 34]. In this respect,
the nature of the current investigation is more theoretical than of practical
applications today, and the standard approaches to modeling galaxy bias work
just fine in most cases. However, the primordial non-Gaussian signal $f_{\rm
NL}$ that might have been imprinted in the early Universe exhibits the
strongest signals in galaxy clustering on scales close to the horizon scales
[35, 36, 37], and its measurements in galaxy clustering are one of the main
goals in the upcoming large-scale surveys such as [38, 39, 40, 41, 42, 43, 44]
with the target uncertainty $\Delta f_{\rm NL}\approx$ few. In particular, the
relativistic effects on such large scales are shown to affect the inference of
the primordial non-Gaussian signal from the measurements [45, 46, 47, 48, 49,
50, 51], and hence a correct modeling of galaxy bias in a relativistic
framework is the last missing piece toward a fully relativistic description of
galaxy clustering [52, 53].

Here we critically investigate the issues in the standard approaches to
modeling galaxy bias from a relativistic perspective and discuss ways to build
galaxy bias models that are compatible with the diffeomorphism symmetry in
general relativity.

##  2 Proper-time hypersurface for modeling galaxy bias

###  2.1 Gauge choice for galaxy bias: temporal and spatial

Since observers can provide unique numbers for $\delta_{g}^{\rm obs}$ in any
given surveys, theorists need to predict unique numbers for $\delta_{g}$ in
any given models to compare to $\delta_{g}^{\rm obs}$. Given the freedom to
choose any coordinate system in general relativity, a simple way to remove the
ambiguities shown in eq. (1.3) due to gauge choice and at the same time to
maintain the linear bias relation in eq. (1.1) as in the Newtonian dynamics is
to fix a gauge for the matter density fluctuation $\delta_{m}$ in eq. (1.1) by
hand. The immediate question then arises “which gauge choice for the
theoretical prediction of $\delta_{m}$ in modeling galaxy bias?” Figure 1
illustrates that a different choice of gauge condition leads to a different
prediction for $\delta_{m}$ (hence $\delta_{g}$), which demands that a choice
of gauge condition needs to be made and physically justified for the linear
bias relation in eq. (1.1) to be valid.

The synchronous gauge was chosen [54] for the linear bias relation in the
relativistic dynamics, and this choice of gauge condition was supported [55,
56] by physical arguments that the only clock available for local observers in
galaxies is the proper time along the world line of the local observers in the
galaxies. The proper-time hypersurface of those observers moving along a time-
like geodesic coincides with those in the synchronous gauge and in the matter-
comoving gauge at the linear order in perturbations [57], though two gauge
conditions are different. The physical arguments for the proper-time
hypersurface in modeling galaxy bias are sound and indeed compatible with the
diffeomorphism symmetry in general relativity, as the four-velocity $u^{\mu}$
of the galaxies defines the time-direction for the local observers in the rest
frame of the galaxies and a hypersurface orthogonal to the four-velocity
defines the proper-time hypersurface in a covariant manner. In this way, the
temporal gauge choice for modeling galaxy bias can be made, not by a specific
hand-picked gauge choice (e.g., the time component of the metric fluctuation
is set zero in the synchronous gauge), but by a physical condition for the
proper-time hypersurface that can be expressed in any gauge choice.

Figure 1: Matter density fluctuations in various gauges. There exist an
infinite number of possibilities for a gauge choice, though the conformal
Newtonian gauge (dashed: $\chi\equiv 0$) and the comoving-synchronous gauge
(solid: $v\equiv 0$) are most popular in literature. Shown in the plot are two
additional gauge choices: uniform curvature (dotted: $\varphi\equiv 0$) and
uniform expansion (dot dashed: $\kappa\equiv 0$). While the differences in the
matter density fluctuations lie mostly on large scales, one can choose a gauge
condition to make it vanish completely, i.e., uniform density gauge
($\delta_{m}\equiv 0$). The notation for metric in the plot follows the
convention described in Appendix [58], and for the plot we assumed the best-
fit cosmological parameter set for a $\Lambda$CDM model from the Planck
measurements [59].

Beyond the linear order in perturbations, a spatial gauge choice also plays an
important role. The diffeomorphism symmetry allows four degrees of freedom in
coordinates that can be chosen by a temporal gauge and a spatial gauge. A
choice of temporal gauge amounts to a choice of a hypersurface of simultaneity
or a time-slicing, and a choice of spatial gauge amounts to a choice of a
spatial gridding in a given hypersurface. In the past, almost exclusively
spatial C-gauge has been used for a spatial gauge choice, in which there is no
off-diagonal component of scalar and vector fluctuations in the spatial
component of the metric tensor (see [57] for the terminology). In certainty
this is not the only choice one can make, and indeed various spatial gauge
choices have been explored in recent years in implementing $N$-body
simulations in a relativistic framework [60, 61, 62, 63] or developing a
relativistic version of the Lagrangian perturbation theory [64, 65]. In
addition to a spatial C-gauge choice, galaxy bias has been also modeled in
spatial B-gauge [45, 46, 66, 67, 68], in which there is no fluctuation in the
space-time component of the metric tensor (see [57] for the terminology).
Other choices of spatial gauge are also allowed in general relativity.

Now consider a general coordinate transformation

$\tilde{x}^{\mu}=e^{\zeta^{\nu}\partial_{\nu}}x^{\mu}=x^{\mu}+\zeta^{\mu}+\frac{1}{2}\zeta^{\mu}{}_{,\nu}\zeta^{\nu}+\mathcal{O}(3)~{},$
(2.1)

parametrized by an infinitesimal vector field $\zeta^{\mu}$. Any perturbation
quantities $\delta{\bf T}$ gauge transform [69, 70, 71] as

$\widetilde{\delta{\bf T}}=\delta{\bf T}-\pounds_{\zeta}\bar{{\bf
T}}+\frac{1}{2}\pounds_{\zeta}\pounds_{\zeta}\bar{{\bf
T}}-\pounds_{\zeta}\delta{\bf T}+\mathcal{O}(3)~{},$ (2.2)

where we again split a field ${\bf T}:=\bar{{\bf T}}+\delta{\bf T}$ into a
background $\bar{{\bf T}}$ and a perturbation $\delta{\bf T}$ and we
suppressed its tensorial indicies. At the linear order in perturbations, only
the Lie derivative $\pounds_{\zeta}\bar{{\bf T}}$ contributes in the gauge
transformation, and it vanishes for a spatial component of $\zeta^{\mu}$, if
$\bar{{\bf T}}$ is a scalar quantity. In other words, due to the homogeneity
and isotropy of the background universe, any background quantities $\bar{{\bf
T}}$ are just a function of time only, and hence a choice of spatial gauge
condition is irrelevant at the linear-order calculations of scalar quantities.
Beyond the linear order in perturbations, a spatial gauge choice makes a
difference in gauge transformation (this is true even at the linear order, if
${\bf T}$ is a tensor).

To be specific, we show that the ambiguities in spatial gauge choice beyond
the linear order in perturbations have a direct impact on the three-point
correlation function or the bispectrum, whose leading order contribution
arises from the second-order contribution in the field. Consider a spatial
gauge transformation $\mathcal{L}^{i}$ in a hypersurface fixed by a choice of
temporal gauge ($T\equiv 0$), and the matter density fluctuation (or any
scalar quantities) then gauge transform as

$\tilde{\delta}_{m}(x^{\mu})=\delta_{m}(x^{\mu})-\mathcal{L}^{i}~{}\nabla_{i}\delta_{m}+\mathcal{O}(3)~{},$
(2.3)

where the time coordinates in $x^{\mu}$ and $\tilde{x}^{\mu}$ are identical.
It is readily evident that the three-point correlation function (or the
bispectrum) also gauge transforms with a spatial gauge choice, even when a
temporal gauge condition is completely fixed.

For example, a single-field consistency relation in inflationary models is
obtained by computing the bispectrum of a comoving-gauge curvature
perturbation $\cal R$ in the squeezed triangular configuration [72, 73, 74,
75]. In this limit, the bispectrum is in proportion to $n_{s}-1$, as long as
there exists only one scalar degree of freedom, regardless of its dynamics,
and hence an observational test of this consistency relation would provide a
powerful way of understanding the early Universe, where $n_{s}$ is the
spectral index of the curvature perturbation. However, note that while $\cal
R$ is independent of a spatial gauge choice at the linear order, it depends on
a spatial gauge choice beyond the linear order. In computing the consistency
relation, a choice of spatial C-gauge was exclusively adopted, though other
choices of spatial gauge are possible. It was shown in [76] that this
consistency relation in the squeezed bispectrum changes according to a spatial
gauge choice and hence can be made in proportion to the running $\alpha_{s}$
of the spectral index, rather than $n_{s}$ itself, by an appropriate gauge
choice. A key point in [76] is that the consistency relation depends on a
spatial gauge choice in addition to the temporal gauge choice, because it is
not a direct observable, on which we further elaborate below.

In modeling galaxy bias, the question was “if we were to fix a gauge choice by
hand, which gauge condition do we have to choose and why?” For a choice of
temporal gauge, the physical argument for a proper-time hypersurface uniquely
fixes a time-slicing that can be expressed in any gauge choice. For a choice
of spatial gauge, however, the physical argument for a proper-time
hypersurface does not fix a spatial gridding. It was shown [57] that while
both the synchronous and the comoving gauges correspond to a proper-time
hypersurface, they differ by a spatial gauge choice, and their matter density
fluctuations are different at the second order in perturbations, leaving
ambiguities in the linear bias relation in eq. (1.1). In summary, there is
currently no physical argument for modeling galaxy bias in literature that
uniquely fixes a spatial gridding and can be expressed in any choice of gauge.

###  2.2 Second-order bias and primordial non-Gaussianity

Beyond the linear order in perturbations, there exist additional contributions
in the standard (Newtonian) galaxy bias model (see, e.g., [77, 78, 79, 21, 22,
26, 27, 28, 6]). Up to the second order, two extra bias parameters are needed

$\delta_{g}=b_{1}\delta_{m}+\frac{1}{2}b_{2}\delta_{m}^{2}+b_{s}s^{2}+\mathcal{O}(3)$
(2.4)

where the traceless tidal tensor is

$s_{ij}:=\left(\Delta^{-1}\nabla_{i}\nabla_{j}-\frac{1}{3}\delta_{ij}\right)\delta_{m}~{},\qquad\qquad
s^{2}:=s_{ij}s_{ij}~{}.$ (2.5)

Additional contributions are added (or prohibited) at each order [21, 22] to
yield contributions to $\delta_{g}$ that is compatible with the symmetry in
the system, or the Galilean symmetry in Newtonian dynamics. This is the
essence of the effective field theory (EFT) approach (see, e.g., [26, 27,
28]), and for instance a gravitational potential contribution $\phi$ in eq.
(2.4) is not allowed as a constant shift in $\phi$ everywhere cannot be
measured according to the equivalence principle.

As discussed in the introduction, this standard modeling of galaxy bias works
well on small scales, but it is incomplete in general relativity. First, a
spatial gauge transformation would yield

$\tilde{\delta}_{g}(x)=b_{1}\delta_{m}+\frac{1}{2}b_{2}\delta_{m}^{2}+b_{s}s^{2}-b_{1}{\cal
L}\cdot\nabla\delta_{m}+\mathcal{O}(3)~{},$ (2.6)

with an extra contribution arising, from the transformation of the matter
fluctuation in eq. (2.3). The spatial transformation ${\cal L}^{i}$ is fully
arbitrary, allowed by diffeomorphism symmetry. Hence this term cannot be
absorbed into the existing bias parameters, signaling the breakdown of the EFT
modeling of galaxy bias in Newtonian dynamics from the perspective of
relativistic dynamics. However, mind that there exists a case, in which this
is possible, as the Newtonian limit is part of general relativity. For a
specific transformation from the Lagrangian frame to the Eulerian frame in the
Newtonian dynamics, the spatial transformation ${\cal L}^{i}$ corresponds to
the displacement field $\Psi^{i}$, i.e., ${\cal
L}^{i}\equiv\Psi^{i}=-\Delta^{-1}\nabla^{i}\delta_{m}$, which can be further
decomposed in terms of $\delta_{m}^{2}$, $s^{2}$, and the second-order matter
density $\delta_{m}^{(2)}$, such that the extra term can be absorbed into the
existing bias parameters. In general, however, we emphasize again that the
spatial transformation is fully arbitrary, not limited to this specific case.

Furthermore, the recent analysis of the relativistic effects reveals [80, 33,
81, 54, 55, 82] that there exist numerous relativistic effects in galaxy
clustering such as the gravitational potential $\phi$. Such contributions are
indeed allowed, given the diffeomorphism symmetry in general relativity. In
fact, compatibility with the equivalence principle, or the independence of a
constant shift in $\phi$ is enforced in galaxy clustering [55, 83, 84, 85,
86], not by prohibiting such contributions, but by allowing other relativistic
contributions of $\phi$ that cancel the effect of a constant shift in $\phi$.
It is evident that an EFT approach to modeling galaxy bias in relativistic
dynamics is more involved than the standard EFT approaches in Newtonian
dynamics. This incompatibility in the standard model for galaxy bias comes
from the fact that the diffeomorphism symmetry in general relativity is not
respected in the standard model for galaxy bias, or there is no covariant
description of galaxy bias.

We stress again that Newtonian dynamics and the standard model for galaxy
bias, though incompatible with general relativity, work well in most cases of
our cosmological interest. As discussed, however, the relativistic effects
become important, when the signals of our cosmological interest are also small
relativistic effects such as the primordial non-Gaussianity. The deviation
from a perfect Gaussianity in the primordial fluctuations contains crucial
clues about the early Universe. The primordial non-Gaussianity is often
parametrized [87, 88] as

${\cal R}(x)={\cal R}_{g}(x)+\frac{3}{5}f_{\rm NL}{\cal R}_{g}^{2}(x)~{},$
(2.7)

a local-type deviation from a Gaussian curvature fluctuation ${\cal R}_{g}$ in
the comoving gauge, and the prediction for the standard single-field
inflationary model is of the order $f_{\rm NL}\sim(n_{s}-1)$ [72], while other
non-standard inflationary models often predict larger non-Gaussianities
$f_{\rm NL}\geq 1$ [89, 90, 74, 91]. Constraining the level of primordial non-
Gaussianity in our Universe is one of the key targets in large-scale surveys
such as the Vera C. Rubin Observatory [38], the Euclid mission [39], and the
Nancy Grace Roman Space Telescope [40].

In the presence of a local-type primordial non-Gaussianity, it was shown [35,
36] that the local rms fluctuations are modulated, affecting the probability
of the matter density fluctuation above the threshold $\delta_{c}\simeq
1.686$, and hence the galaxy number density

$\delta_{g}=b~{}\delta_{m}+2f_{\rm
NL}(b-1)\delta_{c}\phi_{p}+\mathcal{O}(2)~{},$ (2.8)

receives additional contribution from the primordial gravitational potential
$\phi_{p}:=(3/5){\cal R}$. Compared to the first term $b~{}\delta_{m}$, this
extra contribution makes the bias parameter scale dependent with respect to
the matter density, as $\phi_{p}\propto\delta_{m}/k^{2}$, boosting the galaxy
power spectrum on large scales, where the bias parameter is believed to be a
scale-independent constant. Measurements of galaxy clustering in the upcoming
surveys would then provide a clear path to best constraining the primordial
non-Gaussianity, given that more large-scale modes are available for galaxy
surveys than in CMB observations (see, e.g., [38, 39, 40, 34]).

However, at first glance, there exist a few unsatisfactory features: the
galaxy number density in eq. (2.8) responds to a shift in the uniform
gravitational potential, and the second-order correction in eq. (2.7) appears
in eq. (2.8) as the first-order correction. A more rigorous derivation for
galaxy clustering was made [37] in Newtonian dynamics by using the spherical
collapse model in the presence of primordial non-Gaussianity [92], and the
two-point correlation function $\xi_{g}$ for galaxy,

$\xi_{g}(x_{1},x_{2})={\nu^{2}\over\sigma^{2}_{R}}~{}\xi_{m}(x_{1},x_{2})+{\nu^{3}\over\sigma^{3}_{R}}~{}\xi_{m}(x_{1},x_{1},x_{2})~{},$
(2.9)

receives an additional contribution from the three-point matter correlation
function $\xi_{m}$ in the squeezed triangular configuration, where
$\sigma_{R}$ is the rms matter fluctuation smoothed by a scale $R$,
$\nu:=\delta_{c}/\sigma_{R}$ is the scaled threshold for the spherical
collapse model, and the bias parameter is then $b=1+\nu/\sigma_{R}$ in this
model [2, 3, 4, 5].

Indeed, the non-Gaussian correction in eq. (2.9) is the next order in
perturbations, compared to the leading two-point correlation function, as
indicated in eq. (2.7) with a non-vanishing $f_{\rm NL}$. While things are in
order in Newtonian dynamics, this implies that the ambiguities in spatial
gauge due to the transformation in eq. (2.3) affect the correction from the
three-point correlation function in eq. (2.9) in relativistic dynamics. Such
relativistic effects may affect the measurements of $f_{\rm NL}$, at the level
of order unity [48, 45, 46, 47, 49]. Note that the prediction for $f_{\rm NL}$
in the standard inflationary model is negligibly small [72], such that $\Delta
f_{\rm NL}\sim$few from the relativistic effects in general relativity might
be misinterpreted as the signal from inflationary models beyond the standard
model. Given the enormous potential to probe the early Universe, the
theoretical description in eqs. (2.8) or (2.9) needs to be further improved in
a relativistic framework for an accurate prediction for $\Delta f_{\rm NL}$
from the relativistic effects.

###  2.3 Spatial gauge choices for galaxy bias in literature

The linear bias relation in eq. (1.1) is expected to be valid with the matter
density fluctuation beyond the linear order in perturbations, and we need a
choice of spatial gauge in addition to a choice of temporal gauge, if we want
to impose the linear bias relation in a specific gauge. Given that a temporal
gauge choice for a proper-time hypersurface is physically justified and also
universally adopted in literature, we discuss the past choices for spatial
gauge and aruge that there is no good physical argument for a specific choice
of spatial gauge so far in literature. The diffeomorphism symmetry in general
relativity puts all choices on a equivalent footing.

The temporal comoving gauge ($v\equiv 0$) with spatial C-gauge was preferred
in [57] (or gauge-I therein), since the equation of motion for a pressureless
medium in this gauge is identical to one in Newtonian dynamics. However, this
correspondence is valid only up to the second order in perturbations [93, 94,
95], or not valid beyond the second order in perturbations. Furthermore, the
solution in this gauge contains not only the Newtonian contributions, but also
the relativistic contributions [45, 46, 47, 96, 97], which weakens the
arguments for this gauge choice. Furthermore, following the argument in [98],
this gauge choice was advocated in [57], as it preserves the “mass
conservation,” which in fact meant that the leading terms in the spatial
expansion of $\delta_{m}$ add up to zero. The “momentum conservation” was used
in [98] to refer to the vanishing sum of the next-leading terms in the spatial
expansion. Unfortunately, this argument for the spatial C-gauge used in [57]
is just a statement about the coordinate properties, rather than a physical
argument for supporting the gauge choice. Moreover, with the relativistic
effects properly considered in this gauge [45, 46, 47, 96, 97], the condition
for the vanishing sum in the leading terms is in fact not satisfied, i.e., the
“mass conservation” argument used in [57] is invalid. Furthermore, the
synchronous gauge ($\alpha\equiv 0$ with spatial B-gauge or gauge-II in [57])
was shown to correspond to one in the Lagrangian frame (see also [64, 65]),
such that this coordinate is also related by a gauge transformation to the
temporal comoving gauge with spatial C-gauge. Similar to the Newtonian
dynamics, these two gauge choices reflect two different views of the same
system with no preference to one or the other.

Another approach in literature that gained popularity is the conformal Fermi
coordinate (CFC; [99, 100, 101, 102, 103, 67]), in which the local metric
along a given geodesic parametrized by a proper time $\tau$ is
$g_{\mu\nu}=a^{2}(\tau)\eta_{\mu\nu}$ and the corrections to the metric in the
neighborhood appear at least quadratically in distances. It provides a useful
framework to interpret Newtonian calculations in general relativity, as it
closely matches Newtonian coordinates in the neighborhood. However, CFC is not
Newtonian either, simply because the metric around the world line deviates
from the Minkowski metric $\eta_{\mu\nu}$ (or Euclidean even). Numerical
simulations can be used to study galaxy bias with the relativistic effects
accounted for [60, 61, 104, 105, 106]. In particular, the light-cone
observations were made [106, 107] in the simulations, and the galaxy bias
parameters obtained in the light-cone observations were compared to the bias
parameters on a hypersurface in the simulations. While these approaches
provide a powerful tool to study the relativistic effects and the light-cone
observables, the focus of this work is galaxy bias in a relativistic
framework, or the relativistic relation between the galaxy and the matter
density distributions, one step before we account for the light propagation.
In both cases, the spatial C-gauge was chosen in the numerical simulations.

In short, Newtonian correspondence in relativistic dynamics is a good
consistency check, but it is not a sufficient condition for validating a model
in general relativity. In the end, general relativity is not Newtonian. Figure
1 highlights this point that most of the matter fluctuations in various gauges
correspond to the Newtonian matter fluctuation on small scales, while they are
different on large scales. We have found no good physical justification for
any specific choice of spatial gauge for galaxy bias. However, as we show
below, a covariant description of galaxy bias resolves the issue of choosing a
gauge condition in a way that any gauge choice is on a equal footing,
consistent with the diffeomorphism symmetry.

##  3 Gauge-invariance of cosmological observables

###  3.1 Cosmological observables

As opposed to the coordinate dependence discussed so far, the observers
utilize a unique coordinate system (up to a trivial rotation) to map the
cosmological sources in the rest frame, i.e., observed angular positions
${\bf\hat{n}}$ and observed redshifts $z$. These numbers and the observer
coordinates are independent of coordinates we use in our theoretical
description of the Universe, i.e., FRW coordinates. Indeed, the mapping of the
observers is not a coordinate, as it does not cover the whole spacetime
manifold, nor is it invertible (see [52] for further discussion). In terms of
galaxy clustering, the observers simply count the number of the observed
galaxies in a small volume defined by a small redshift bin and a small solid
angle in observations, and this number yields the observed galaxy number
density as a function of the observed position $(z,n^{i})$ [80]. Therefore, it
is evident that our theoretical description for the observed galaxy number
density fluctuation should be a scalar invariant under diffeomorphism

$\delta^{\rm obs}_{g}(x^{\mu}_{s})=\tilde{\delta}^{\rm
obs}_{g}(\tilde{x}^{\mu}_{s})~{},$ (3.1)

or gauge-invariant when expressed in terms of the observed position [71] to
properly compare to observations. This is indeed the case [108, 109, 110, 111,
112, 113, 114, 71, 52, 115] as the source galaxy position in a given
coordinate is

$x^{\mu}_{s}=:\bar{x}^{\mu}_{z}+\Delta x^{\mu}~{},$ (3.2)

where the observed position of the source specified by $(z,n^{i})$ in a
coordinate is expressed as $\bar{x}^{\mu}_{z}$ and its deviation from the real
position $x^{\mu}_{s}$ in a coordinate is captured by $\Delta x^{\mu}$. For
instance, the expression up to the second order in perturbations is

$\delta_{g}^{\rm obs}(x_{s}^{\mu})=\delta^{\rm
obs}_{g}(\bar{x}_{z}^{\mu})+\Delta x^{\mu}\partial_{\mu}\delta_{g}^{\rm
obs}+\mathcal{O}(3)~{}.$ (3.3)

With the observed position $\bar{x}^{\mu}_{z}$ fixed as a reference in all
coordinates, the deviation $\Delta x^{\mu}$ gauge transforms, according to the
coordinate transformation in eq. (2.1), and this transformation cancels any
temporal and spatial transformation of $\delta_{g}^{\rm obs}(x^{\mu}_{s})$ in
eq. (2.2) [or Eqs. (1.3) and (2.3)], rendering the observed galaxy number
density fluctuation fully gauge invariant at the observed position.

The theoretical expression for the observed galaxy number density fluctuation
$\delta_{g}^{\rm obs}$ is derived by solving the geodesic equation and
relating the observed position $\bar{x}^{\mu}_{z}$ to the source position
$x^{\mu}$ in a given coordinate, which uniquely determines the volume effect
[116] such as the gravitational lensing [117] and the redshift-space
distortion [118]. However, being a key part of the source effect, the galaxy
bias part in the theoretical description remains ambiguous beyond the linear
order in perturbations due to the lack of covariant galaxy bias model.

###  3.2 A covariant toy model for galaxy bias

It is in fact simple to devise a covariant model for galaxy bias, while the
difficult part lies in building a realistic model. For instance, the simplest
one is to model the galaxy distribution as the matter distribution or a
function of the matter distribution:

$n_{g}(x^{\mu})\equiv{\cal N}\rho_{m}^{b}(x^{\mu})~{},$ (3.4)

where we assumed a power-law function for simplicity and $\cal N$ is a
dimensionful constant. This simple model is covariant, and the bias relation
in eq. (1.1) is valid in all coordinates, with $b$ in the power-law as the
linear bias parameter:

$n_{g}(x^{\mu})={\cal
N}\bar{\rho}_{m}^{b}(t)\left[1+b~{}\delta_{m}(x)+\cdots\right]~{}.$ (3.5)

The evolution of the galaxy number density in this model is, however, related
to the bias parameter $b$, according to eq. (1.4). Hence, the model is not
observationally viable, which illustrates an extra complexity in developing a
covariant model for galaxy bias — one has to account for the time evolution of
the galaxy number density in building a model for galaxy bias. Note that the
time evolution of galaxy bias can be accounted for in the standard approaches
(e.g., [79]), but this aspect is not a necessary condition in the standard
model.

To implement the physical arguments for the proper-time hypersurface into the
model and to match the observational constraint on the evolution of the galaxy
number density, we develop another simple toy model by considering an extra
function $E(t_{p})$ of proper-time $t_{p}$ of the galaxy along the world line:

$n_{g}(x^{\mu})\equiv{\cal N}\rho_{m}^{b}(x^{\mu})E(t_{p})=:{\cal
N}\rho_{m}^{b}(x^{\mu})(1+z_{p})^{b_{t}}~{},$ (3.6)

where we again assumed a simple power-law function for $E(t_{p})$ and
expressed it in terms of the redshift parameter $z_{p}$ in the proper-time
hypersurface, and $\cal N$ is another dimensionful constant of the model. The
functional dependence of the model on $t_{p}$ prefers the proper-time
hypersurface, and the evolution of the galaxy number density is

${d\ln\bar{n}_{g}\over d\ln(1+z)}=3b+b_{t}~{},$ (3.7)

instead of $3b$, which allows freedom to match the observed evolution. The
galaxy number density can be expanded in the proper-time hypersurface as

$n_{g}(x^{\mu})={\cal
N}\bar{\rho}_{m}^{b}(t_{p})(1+z_{p})^{b_{t}}\left[1+b~{}\delta_{p}(x)+\cdots\right]~{},$
(3.8)

and a choice of spatial gauge for galaxy bias is not needed, as it is
covariant in this model, where $\delta_{p}$ is the matter density fluctuation
in the proper-time hypersurface. To illustrate the difference, the question in
the standard approaches is to use the matter density fluctuation with which
gauge conditions:

$\delta_{g}(x)=b~{}\delta_{m}^{\rm
I}(x)~{}?\qquad\qquad\delta_{g}(x)=b~{}\delta_{m}^{\rm
II}(x)~{}?\qquad\qquad\delta_{g}(x)=\cdots?$ (3.9)

where the super-scripts I and II indicate the matter density fluctuation
computed in various choices of gauge conditions, either temporal or spatial.
Note that the matter density fluctuation $\delta_{m}(x)$ is uniquely defined
in each gauge and they are all different. According to the covariant toy
model, the galaxy bias relation in contrast becomes

$\delta_{g}^{\rm I}(x^{\rm I})=b~{}\delta_{p}^{\rm I}(x^{\rm
I})~{},\qquad\qquad\delta_{g}^{\rm II}(x^{\rm II})=b~{}\delta_{p}^{\rm
II}(x^{\rm II}),~{}\qquad\qquad\delta_{g}^{\rm I}(x^{\rm I})=\delta_{g}^{\rm
II}(x^{\rm II})=\cdots~{},$ (3.10)

such that the same bias factor is multiplied by the matter density fluctuation
in each gauge choice and the galaxy fluctuation at a given physical point (not
a coordinate position, $x^{\rm I}\neq x^{\rm II}$) is identical in any gauge
choice. Therefore, any spatial gauge choice is allowed with the same bias
parameter $b$ and no physical justification for a specific choice of spatial
gauge is needed, resolving the issues associated with a spatial gauge choice.

Consequently, when the galaxy two-point correlation
$\langle\delta_{g}\delta_{g}\rangle$ (or for that matter the $N$-point
correlation as well) is computed, the outcome will differ, depending on which
gauge choice is adopted, as shown in Figure 1. However, as discussed in
Section 3.1 or shown in eq. (3.1), the outcome will be identical and
independent of gauge choice, when expressed in terms of the observed position
for a covariant model for galaxy bias.

While these toy models serve the purpose of demonstrating the need to develop
a covariant model, they are not realistic models for galaxy bias. A realistic
covariant model for galaxy bias can be developed by generalizing, for example,
the peak model [4, 11, 24]. As discussed, a temporal gauge choice to the
proper-time hypersurface can be enforced not by choosing a coordinate by hand,
but by adopting a four velocity $u^{\mu}$ of the galaxy sample as the time
direction in the model. The matter overdensity in the proper-time hypersurface
should then satisfy the condition for a collapse $\delta_{m}\geq\delta_{c}$,
and the other condition that the peak position should have a vanishing spatial
derivative and a negative second spatial derivative can be accommodated in a
covariant manner by defining the spatial projection tensor
$\mathcal{H}_{\mu\nu}:=g_{\mu\nu}+u_{\mu}u_{\nu}$ and by using the projected
covariant derivative $\nabla_{\mu}^{(3)}:=\mathcal{H}^{\nu}_{\mu}\nabla_{\nu}$
as a spatial derivative in the proper-time hypersurface, where $g_{\mu\nu}$ is
the metric tensor. In this way it would be rather straightforward to build a
covariant peak model, but of course quantifying the impact of the covariant
peak model compared to the standard model requires a further investigation,
beyond the scope of the current manuscript.

##  4 Conclusions

The standard model for galaxy bias (and many other elements in cosmology as
well) is built on a perturbative framework in Newtonian cosmology, and hence
it lacks a covariant expression for galaxy bias. We have demonstrated that the
reverse-engineered approaches to modeling galaxy bias in literature do not
respect the diffeomorphism symmetry in general relativity.

To mitigate this issue, we have considered a simple, but covariant toy model
in eq. (3.6) that respects the diffeomorphism symmetry in general relativity
and hence resolves all the ambiguities in the standard model for galaxy bias.
One obvious deficiency in the toy model is that galaxies can form even in the
background universe — It is the matter fluctuations in the local neighborhood,
not the absolute value of the matter density that matters for galaxy
formation; in the early Universe, galaxies or even dark matter halos cannot
form, despite the enormous matter density everywhere. Though this
consideration can be put into the model by adjusting the function $E(t_{p})$
in eq. (3.6), a more physical model for galaxy bias remains to be built. In
particular, building a realistic model for galaxy bias, especially in the
presence of the primordial non-Gaussianity, deserves a further in-depth
investigation, and a possibility is to put the peak model for galaxy bias in a
relativistic framework.

Though strictly true in principle, the flaws in the standard approaches pose
no significant issues in practice, as the applications of galaxy clustering in
most cases are well described by the Newtonian dynamics. This work explores
the subtle issues in a regime, where the relativistic effects become
important. Given the level of the measurement uncertainties $\Delta f_{\rm
NL}\approx$ few in measuring the primordial non-Gaussianity in the upcoming
large-scale surveys (e.g., [38, 39, 40, 41, 42, 43, 44]), a covariant
relativistic description of galaxy bias will play a crucial role in extracting
such small relativistic effects that contain yet important cosmological
information [53, 50, 51].

## Acknowledgments

We acknowledge useful discussions with Ermis Mitsou and Julian Adamek. This
work is supported by the Swiss National Science Foundation and a Consolidator
Grant of the European Research Council.

## References

  * Press and Schechter [1974] W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).
  * Kaiser [1984] N. Kaiser, Astrophys. J. Lett. 284, L9 (1984).
  * Politzer and Wise [1984] H. D. Politzer and M. B. Wise, Astrophys. J. Lett. 285, L1 (1984).
  * Bardeen et al. [1986] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * Bond et al. [1991] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Astrophys. J. 379, 440 (1991).
  * Desjacques et al. [2018] V. Desjacques, D. Jeong, and F. Schmidt, Phys. Rep. 733, 1 (2018), 1611.09787.
  * Eisenstein et al. [2005] D. J. Eisenstein, I. Zehavi, D. W. Hogg, R. Scoccimarro, et al., Astrophys. J. 633, 560 (2005), arXiv:astro-ph/0501171.
  * Tegmark et al. [2006] M. Tegmark et al., Phys. Rev. D 74, 123507 (2006), arXiv:astro-ph/0608632.
  * Anderson et al. [2012] L. Anderson, E. Aubourg, S. Bailey, D. Bizyaev, M. Blanton, A. S. Bolton, J. Brinkmann, J. R. Brownstein, et al., Mon. Not. R. Astron. Soc. 427, 3435 (2012), 1203.6594.
  * Aubourg et al. [2015] É. Aubourg, S. Bailey, J. E. Bautista, F. Beutler, et al., Phys. Rev. D 92, 123516 (2015), 1411.1074.
  * Bond and Myers [1996] J. R. Bond and S. T. Myers, Astrophys. J. Suppl. Ser. 103, 1 (1996).
  * Mo and White [1996] H. J. Mo and S. D. M. White, Mon. Not. R. Astron. Soc. 282, 347 (1996), arXiv:9512127.
  * Matsubara [1999] T. Matsubara, Astrophys. J. 525, 543 (1999), arXiv:9906029.
  * Catelan et al. [2000] P. Catelan, C. Porciani, and M. Kamionkowski, Mon. Not. R. Astron. Soc. 318, L39 (2000), astro-ph/0005544.
  * Ma and Fry [2000] C.-P. Ma and J. N. Fry, Astrophys. J. 543, 503 (2000), astro-ph/0003343.
  * Seljak [2000] U. Seljak, Mon. Not. R. Astron. Soc. 318, 203 (2000), astro-ph/0001493.
  * Sheth et al. [2001] R. K. Sheth, H. J. Mo, and G. Tormen, Mon. Not. R. Astron. Soc. 323, 1 (2001), astro-ph/9907024.
  * Scoccimarro et al. [2001] R. Scoccimarro, R. K. Sheth, L. Hui, and B. Jain, Astrophys. J. 546, 20 (2001), astro-ph/0006319.
  * Berlind and Weinberg [2002] A. A. Berlind and D. H. Weinberg, Astrophys. J. 575, 587 (2002), astro-ph/0109001.
  * Cooray and Sheth [2002] A. Cooray and R. Sheth, Phys. Rep. 372, 1 (2002), astro-ph/0206508.
  * McDonald [2006] P. McDonald, Phys. Rev. D 74, 103512 (2006), arXiv:astro-ph/0609413.
  * McDonald and Roy [2009] P. McDonald and A. Roy, J. Cosmol. Astropart. Phys. 8, 20 (2009), 0902.0991.
  * Baldauf et al. [2012] T. Baldauf, U. Seljak, V. Desjacques, and P. McDonald, Phys. Rev. D 86, 083540 (2012), 1201.4827.
  * Desjacques [2013] V. Desjacques, Phys. Rev. D 87, 043505 (2013), 1211.4128.
  * Schmidt et al. [2013] F. Schmidt, D. Jeong, and V. Desjacques, Phys. Rev. D 88, 023515 (2013), 1212.0868.
  * Senatore [2015] L. Senatore, J. Cosmol. Astropart. Phys. 2015, 007 (2015), 1406.7843.
  * Mirbabayi et al. [2015] M. Mirbabayi, F. Schmidt, and M. Zaldarriaga, J. Cosmol. Astropart. Phys. 2015, 030 (2015), 1412.5169.
  * Angulo et al. [2015] R. Angulo, M. Fasiello, L. Senatore, and Z. Vlah, J. Cosmol. Astropart. Phys. 2015, 029 (2015), 1503.08826.
  * Ivanov et al. [2020] M. M. Ivanov, M. Simonović, and M. Zaldarriaga, J. Cosmol. Astropart. Phys. 2020, 042 (2020), 1909.05277.
  * d’Amico et al. [2020] G. d’Amico, J. Gleyzes, N. Kokron, K. Markovic, L. Senatore, P. Zhang, F. Beutler, and H. Gil-Marín, J. Cosmol. Astropart. Phys. 2020, 005 (2020), 1909.05271.
  * Eisenstein et al. [2001] D. J. Eisenstein, J. Annis, J. E. Gunn, A. S. Szalay, A. J. Connolly, R. C. Nichol, N. A. Bahcall, M. Bernardi, S. Burles, F. J. Castander, et al., Astron. J. 122, 2267 (2001), arXiv:astro-ph/0108153.
  * White et al. [2011] M. White, M. Blanton, A. Bolton, D. Schlegel, J. Tinker, A. Berlind, L. da Costa, E. Kazin, Y.-T. Lin, M. Maia, et al., Astrophys. J. 728, 126 (2011), 1010.4915.
  * Yoo [2010] J. Yoo, Phys. Rev. D 82, 083508 (2010), arXiv:1009.3021.
  * Yoo et al. [2012] J. Yoo, N. Hamaus, U. Seljak, and M. Zaldarriaga, Phys. Rev. D 86, 063514 (2012), 1206.5809.
  * Dalal et al. [2008] N. Dalal, O. Doré, D. Huterer, and A. Shirokov, Phys. Rev. D 77, 123514 (2008), 0710.4560.
  * Slosar et al. [2008] A. Slosar, C. Hirata, U. Seljak, S. Ho, and N. Padmanabhan, J. Cosmol. Astropart. Phys. 8, 31 (2008), 0805.3580.
  * Matarrese and Verde [2008] S. Matarrese and L. Verde, Astrophys. J. Lett. 677, L77 (2008), 0801.4826.
  * Stubbs et al. [2004] C. W. Stubbs, D. Sweeney, J. A. Tyson, and LSST Collaboration, in _American Astronomical Society Meeting Abstracts_ (2004), vol. 36 of _Bulletin of the American Astronomical Society_ , p. 108.02.
  * Laureijs et al. [2011] R. Laureijs, J. Amiaux, S. Arduini, J. . Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, et al., ArXiv e-prints (2011), 1110.3193.
  * Green et al. [2012] J. Green, P. Schechter, C. Baltay, R. Bean, D. Bennett, R. Brown, C. Conselice, M. Donahue, et al. (2012), 1208.4012.
  * Cosmic Visions 21 cm Collaboration et al. [2018] Cosmic Visions 21 cm Collaboration, R. Ansari, E. J. Arena, K. Bandura, P. Bull, E. Castorina, T.-C. Chang, Chen, et al., arXiv e-prints arXiv:1810.09572 (2018), 1810.09572.
  * Mueller et al. [2022] E.-M. Mueller, M. Rezaie, W. J. Percival, A. J. Ross, et al., Mon. Not. R. Astron. Soc. 514, 3396 (2022).
  * Cabass et al. [2022] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Phys. Rev. D 106, 043506 (2022), 2204.01781.
  * D’Amico et al. [2022] G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, arXiv e-prints arXiv:2201.11518 (2022), 2201.11518.
  * Bruni et al. [2014a] M. Bruni, J. C. Hidalgo, N. Meures, and D. Wands, Astrophys. J. 785, 2 (2014a), 1307.1478.
  * Bruni et al. [2014b] M. Bruni, J. C. Hidalgo, and D. Wands, Astrophys. J. Lett. 794, L11 (2014b), 1405.7006.
  * Villa et al. [2014] E. Villa, L. Verde, and S. Matarrese, Classical and Quantum Gravity 31, 234005 (2014), 1409.4738.
  * Camera et al. [2015] S. Camera, R. Maartens, and M. G. Santos, Mon. Not. R. Astron. Soc. 451, L80 (2015), 1412.4781.
  * Matarrese et al. [2021] S. Matarrese, L. Pilo, and R. Rollo, J. Cosmol. Astropart. Phys. 2021, 062 (2021), 2007.08877.
  * Castorina and Di Dio [2022] E. Castorina and E. Di Dio, J. Cosmol. Astropart. Phys. 2022, 061 (2022), 2106.08857.
  * Foglieni et al. [2023] M. Foglieni, M. Pantiri, E. Di Dio, and E. Castorina, arXiv e-prints arXiv:2303.03142 (2023), 2303.03142.
  * Mitsou and Yoo [2020] E. Mitsou and J. Yoo, Springer Briefs in Physics arXiv:1908.10757 (2020), 1908.10757.
  * Yoo et al. [2022] J. Yoo, N. Grimm, and E. Mitsou, J. Cosmol. Astropart. Phys. 2022, 050 (2022), 2204.03002.
  * Challinor and Lewis [2011] A. Challinor and A. Lewis, Phys. Rev. D 84, 043516 (2011), arXiv:1105.5292.
  * Jeong et al. [2012] D. Jeong, F. Schmidt, and C. M. Hirata, Phys. Rev. D 85, 023504 (2012), arXiv:1107.5427.
  * Jeong and Schmidt [2014] D. Jeong and F. Schmidt, Phys. Rev. D 89, 043519 (2014), 1305.1299.
  * Yoo [2014a] J. Yoo, Phys. Rev. D 90, 123507 (2014a), 1408.5137.
  * Yoo et al. [2018] J. Yoo, N. Grimm, E. Mitsou, A. Amara, and A. Refregier, J. Cosmol. Astropart. Phys. 4, 029 (2018), 1802.03403.
  * Planck Collaboration et al. [2020] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, et al., Astron. Astrophys. 641, A6 (2020), 1807.06209.
  * Fidler et al. [2015] C. Fidler, C. Rampf, T. Tram, R. Crittenden, K. Koyama, and D. Wands, Phys. Rev. D 92, 123517 (2015), 1505.04756.
  * Fidler et al. [2016] C. Fidler, T. Tram, C. Rampf, R. Crittenden, K. Koyama, and D. Wands, Journal of Cosmology and Astro-Particle Physics 2016, 031 (2016), 1606.05588.
  * Fidler et al. [2017] C. Fidler, T. Tram, C. Rampf, R. Crittenden, K. Koyama, and D. Wands, Journal of Cosmology and Astro-Particle Physics 2017, 022 (2017), 1708.07769.
  * Adamek and Fidler [2019] J. Adamek and C. Fidler, arXiv e-prints p. arXiv:1905.11721 (2019), 1905.11721.
  * Rampf [2014] C. Rampf, Phys. Rev. D 89, 063509 (2014), 1307.1725.
  * Rampf and Wiegand [2014] C. Rampf and A. Wiegand, Phys. Rev. D 90, 123503 (2014), 1409.2688.
  * Koyama et al. [2018] K. Koyama, O. Umeh, R. Maartens, and D. Bertacca, J. Cosmol. Astropart. Phys. 2018, 050 (2018), 1805.09189.
  * Umeh et al. [2019] O. Umeh, K. Koyama, R. Maartens, F. Schmidt, and C. Clarkson, J. Cosmol. Astropart. Phys. 2019, 020 (2019), 1901.07460.
  * Umeh and Koyama [2019] O. Umeh and K. Koyama, J. Cosmol. Astropart. Phys. 2019, 048 (2019), 1907.08094.
  * Stewart and Walker [1974] J. M. Stewart and M. Walker, Proceedings of the Royal Society of London Series A 341, 49 (1974).
  * Bruni et al. [1997] M. Bruni, S. Matarrese, S. Mollerach, and S. Sonego, Classical and Quantum Gravity 14, 2585 (1997), gr-qc/9609040.
  * Yoo and Durrer [2017] J. Yoo and R. Durrer, J. Cosmol. Astropart. Phys. 9, 016 (2017), 1705.05839.
  * Maldacena [2003] J. Maldacena, J. High Energy Phys. 5, 13 (2003), astro-ph/0210603.
  * Creminelli and Zaldarriaga [2004] P. Creminelli and M. Zaldarriaga, J. Cosmol. Astropart. Phys. 10, 6 (2004), arXiv:astro-ph/0407059.
  * Chen et al. [2007] X. Chen, M.-x. Huang, S. Kachru, and G. Shiu, J. Cosmol. Astropart. Phys. 2007, 002 (2007), hep-th/0605045.
  * Senatore and Zaldarriaga [2012] L. Senatore and M. Zaldarriaga, J. Cosmol. Astropart. Phys. 2012, 001 (2012), 1203.6884.
  * Mitsou and Yoo [2022] E. Mitsou and J. Yoo, Physics Letters B 828, 137018 (2022), 2109.13154.
  * Szalay [1988] A. S. Szalay, Astrophys. J. 333, 21 (1988).
  * Fry and Gaztanaga [1993] J. N. Fry and E. Gaztanaga, Astrophys. J. 413, 447 (1993), astro-ph/9302009.
  * Fry [1996] J. N. Fry, Astrophys. J. Lett. 461, L65+ (1996).
  * Yoo et al. [2009] J. Yoo, A. L. Fitzpatrick, and M. Zaldarriaga, Phys. Rev. D 80, 083514 (2009), arXiv:0907.0707.
  * Bonvin and Durrer [2011] C. Bonvin and R. Durrer, Phys. Rev. D 84, 063505 (2011), arXiv:1105.5280.
  * Yoo [2014b] J. Yoo, Class. Quant. Grav. 31, 234001 (2014b), arXiv:1409.3223.
  * Biern and Yoo [2017] S. G. Biern and J. Yoo, J. Cosmol. Astropart. Phys. 026 (2017), 1704.07380.
  * Scaccabarozzi et al. [2018] F. Scaccabarozzi, J. Yoo, and S. G. Biern, J. Cosmol. Astropart. Phys. 10, 024 (2018), 1807.09796.
  * Grimm et al. [2020] N. Grimm, F. Scaccabarozzi, J. Yoo, S. G. Biern, and J.-O. Gong, J. Cosmol. Astropart. Phys. 2020, 064 (2020), 2005.06484.
  * Baumgartner and Yoo [2021] S. Baumgartner and J. Yoo, Phys. Rev. D 103, 063516 (2021), 2012.03968.
  * Komatsu and Spergel [2001] E. Komatsu and D. N. Spergel, Phys. Rev. D 63, 063002 (2001), astro-ph/0005036.
  * Bartolo et al. [2004] N. Bartolo, E. Komatsu, S. Matarrese, and A. Riotto, Phys. Rep. 402, 103 (2004), arXiv:0406398.
  * Seery and Lidsey [2005] D. Seery and J. E. Lidsey, J. Cosmol. Astropart. Phys. 6, 003 (2005), astro-ph/0503692.
  * Byrnes et al. [2006] C. T. Byrnes, M. Sasaki, and D. Wands, Phys. Rev. D 74, 123519 (2006), arXiv:astro-ph/0611075.
  * De Felice and Tsujikawa [2011] A. De Felice and S. Tsujikawa, J. Cosmol. Astropart. Phys. 2011, 029 (2011), 1103.1172.
  * Matarrese et al. [1986] S. Matarrese, F. Lucchin, and S. A. Bonometto, Astrophys. J. Lett. 310, L21 (1986).
  * Hwang and Noh [2005a] J.-C. Hwang and H. Noh, Phys. Rev. D 72, 044011 (2005a), arXiv/0412128.
  * Hwang and Noh [2005b] J.-C. Hwang and H. Noh, Phys. Rev. D 72, 044012 (2005b), gr-qc/0412129.
  * Hwang and Noh [2007] J.-c. Hwang and H. Noh, J. Cosmol. Astropart. Phys. 12, 003 (2007), 0704.2086.
  * Bertacca et al. [2015] D. Bertacca, N. Bartolo, M. Bruni, K. Koyama, R. Maartens, S. Matarrese, M. Sasaki, and D. Wands, Classical and Quantum Gravity 32, 175019 (2015), 1501.03163.
  * Yoo and Gong [2016] J. Yoo and J.-O. Gong, J. Cosmol. Astropart. Phys. 7, 017 (2016), 1602.06300.
  * Peebles [1980] P. J. E. Peebles, _The large-scale structure of the universe_ (Princeton University Press, Princeton, 1980).
  * Pajer et al. [2013] E. Pajer, F. Schmidt, and M. Zaldarriaga, Phys. Rev. D 88, 083502 (2013), 1305.0824.
  * Dai et al. [2015a] L. Dai, E. Pajer, and F. Schmidt, ArXiv e-prints (2015a), 1504.00351.
  * Dai et al. [2015b] L. Dai, E. Pajer, and F. Schmidt, J. Cosmol. Astropart. Phys. 2015, 043 (2015b), 1502.02011.
  * Cabass et al. [2017] G. Cabass, E. Pajer, and F. Schmidt, J. Cosmol. Astropart. Phys. 2017, 003 (2017), 1612.00033.
  * Ip and Schmidt [2017] H. Y. Ip and F. Schmidt, J. Cosmol. Astropart. Phys. 02, 025 (2017), 1610.01059.
  * Adamek et al. [2016a] J. Adamek, D. Daverio, R. Durrer, and M. Kunz, Nature Physics 12, 346 (2016a), 1509.01699.
  * Adamek et al. [2016b] J. Adamek, D. Daverio, R. Durrer, and M. Kunz, J. Cosmol. Astropart. Phys. 2016, 053 (2016b), 1604.06065.
  * Borzyszkowski et al. [2017] M. Borzyszkowski, D. Bertacca, and C. Porciani, Mon. Not. R. Astron. Soc. 471, 3899 (2017), 1703.03407.
  * Lepori et al. [2023] F. Lepori, S. Schulz, J. Adamek, and R. Durrer, J. Cosmol. Astropart. Phys. 2023, 036 (2023), 2209.10533.
  * Yoo and Zaldarriaga [2014] J. Yoo and M. Zaldarriaga, Phys. Rev. D 90, 023513 (2014), 1406.4140.
  * Di Dio et al. [2014] E. Di Dio, R. Durrer, G. Marozzi, and F. Montanari, J. Cosmol. Astropart. Phys. 12, 017 (2014), 1407.0376.
  * Bertacca et al. [2014a] D. Bertacca, R. Maartens, and C. Clarkson, J. Cosmol. Astropart. Phys. 9, 037 (2014a), 1405.4403.
  * Bertacca et al. [2014b] D. Bertacca, R. Maartens, and C. Clarkson, J. Cosmol. Astropart. Phys. 11, 013 (2014b), 1406.0319.
  * Di Dio et al. [2016] E. Di Dio, R. Durrer, G. Marozzi, and F. Montanari, J. Cosmol. Astropart. Phys. 2016, 016 (2016), 1510.04202.
  * Umeh et al. [2017] O. Umeh, S. Jolicoeur, R. Maartens, and C. Clarkson, J. Cosmol. Astropart. Phys. 2017, 034 (2017), 1610.03351.
  * Jolicoeur et al. [2017] S. Jolicoeur, O. Umeh, R. Maartens, and C. Clarkson, J. Cosmol. Astropart. Phys. 2017, 040 (2017), 1703.09630.
  * Magi and Yoo [2022] M. Magi and J. Yoo, J. Cosmol. Astropart. Phys. 2022, 071 (2022), 2204.01751.
  * Yoo [2009] J. Yoo, Phys. Rev. D 79, 023517 (2009), arXiv:0808.3138.
  * Kaiser [1992] N. Kaiser, Astrophys. J. 388, 272 (1992).
  * Kaiser [1987] N. Kaiser, Mon. Not. R. Astron. Soc. 227, 1 (1987).