# Inverse Molecular Design with
Multi-Conditional Diffusion Guidance

Gang Liu, Jiaxin Xu, Tengfei Luo, Meng Jiang
University of Notre Dame
{gliu7, jxu24, tluo<EMAIL_ADDRESS>

###### Abstract

Inverse molecular design with diffusion models holds great potential for
advancements in material and drug discovery. Despite success in unconditional
molecule generation, integrating multiple properties such as synthetic score
and gas permeability as condition constraints into diffusion models remains
unexplored. We introduce multi-conditional diffusion guidance. The proposed
Transformer-based denoising model has a condition encoder that learns the
representations of numerical and categorical conditions. The denoising model,
consisting of a structure encoder-decoder, is trained for denoising under the
representation of conditions. The diffusion process becomes graph-dependent to
accurately estimate graph-related noise in molecules, unlike the previous
models that focus solely on the marginal distributions of atoms or bonds. We
extensively validate our model for multi-conditional polymer and small
molecule generation. Results demonstrate our superiority across metrics from
distribution learning to condition control for molecular properties. An
inverse polymer design task for gas separation with feedback from domain
experts further demonstrates its practical utility.

##  1 Introduction

Diffusion models for molecular graphs are essential for inverse design of
materials and drugs by generating molecules and polymers (macro-molecules)
(thornton2012polymer; wu2018moleculenet), because the models can be
effectively trained to predict discrete graph structures and atom/bond types
in denoising processes (vignac2022digress). Practical inverse designs consider
multiple factors such as molecular synthetic score and various properties
(gebauer2022inverse), known as the task of multi-conditional graph generation.

Existing work converted multiple conditions into a single one and solved the
task as single-condition generation (bilodeau2022generative;
lee2023exploring). However, multi-property relations may not be properly or
explicitly defined (bilodeau2022generative). First, the properties have
diverse scales and units. For example, the synthetic complexity ranges from 1
to 5 (coley2018scscore), while the gas permeability varies widely, exceeding
10,000 in Barrier units (barnett2020designing). This gap makes it hard for
models to balance the conditions. Second, multi-conditions consist of a mix of
categorical and numerical properties. The common practice of addition
(xie2021mars) or multiplication (lee2023exploring) is inadequate for
combination.

Figure 1(a) empirically illustrates the challenges in multi-conditional
generation, i.e., discovering molecules meeting multiple properties. We used a
test set of 100 data points with three properties: synthesizability (Synth.)
(ertl2009estimation), O2 and N2 permeability (O2Perm and N2Perm)
(barnett2020designing). A single-conditional diffusion model generated up to
30 graphs for each condition, resulting in a total of 90 graphs for three
conditions. The 30 graphs in each set were ranked by an Oracle for the polymer
property (see LABEL:add:subsec:motivation-setup). We defined $K$ ($1\leq K\leq
30$) as the highest rank at which a polymer was found in all three sets and
met the multi-property constraints. If not, we set $K$ as 30. Figure 1(a)
shows the frequency distribution of $K$ on the 100 test cases. The median $K$
was 30, indicating that the multiple properties were not met on over half of
the test polymers despite generating a large number of graphs.

(a) Existing work’s limitation: A median rank of 30 showed that on fewer than
half test polymers, the sets of generated graphs from different single
conditions intersected, indicating a failure to generate polymers meeting
multiple properties.

(b) Proposed work: Our idea, multi-conditional guidance for diffusion models,
successfully generated polymers that satisfied multi-property constraints. It
achieved a higher rank than 30 in any set of the single-conditional generated
graphs.

Figure 1: Multi-conditional diffusion guidance in (b) generates polymers of
higher property accuracy than existing work in (a). Explanations are in
Section 1 and details are in LABEL:add:subsec:motivation-setup.

Our idea for multi-conditional inverse design is to project multiple
constraints into representation by _learning_ , thereby guiding the process of
molecule generation in diffusion models. It is distinct from using multi-
conditions in a manual combination or a separate approach. We propose Multi-
Conditional Diffusion (MCD) that has a novel denoising model for predictor-
free diffusion guidance. The denoising model consists of a condition encoder
to convert multiple conditions into a representation. It then uses the
adaptive layer normalization (huang2017arbitrary) to substitute the molecule
representation statistics (e.g., mean and variance) with those of the
condition representation in the Transformer-based structure encoder. Finally,
a multi-layer perceptron (MLP) based structure decoder predicts the denoised
state’s atom and bond probabilities for one-step denoising. To accurately
estimate discrete noise on molecules, we propose a graph-dependent noise model
that aligns with the molecular graph-dependent denoising model. It establishes
a noise transition matrix for four transitions: atom-to-atom, atom-to-bond,
bond-to-atom, and bond-to-bond. Therefore, the noise is applied to the atom-
bond dependent graph, unlike many existing works that apply it separately to
atoms or bonds (jo2022score; vignac2022digress).

Results in Figure 1(b) show that the polymers generated by MCD closely align
with multi-property constraints. For each test case, we have one graph
generated from MCD conditional on three properties. The Oracle determines the
rank of this graph among 30 single-conditionally generated graphs for each
condition. We find the median ranks are 4, 9, and 11, for Synth., O2 Perm, and
N2 Perm, respectively, all much higher than 30. Note that the ranked set of 30
graphs was very competitive because the model was trained on the specific
condition dedicatedly.

In experiments, we evaluate model performance on one polymer and three small
molecule datasets. The polymer dataset includes four numerical conditions for
multi-conditional evaluation. Our model has the lowest average mean absolute
error (MAE), significantly reducing the error by 17.86% compared to the best
baseline. It also excels in small molecule tasks, achieving over 0.9 accuracy
on task-related categorical conditions, notably surpassing the baseline
accuracy of less than 0.6. We also examine the model’s utility in inverse
polymer designs for O2/N2 gas separation, with domain expert feedback
highlighting our model’s practical utility in multi-conditional molecular
design.

##  2 Problem Definition

###  2.1 Multi-Conditional Inverse Molecular Design

A molecular graph $G=(V,E)$ consists of a set of nodes (atoms) $V$ and edges
(bonds) $E$. We follow (vignac2022digress) and define “non-bond” as a type of
edge. There are $N$ atoms and each atom has a one-hot encoding, denoting the
atom type. We represent it as $\mathbf{X}_{V}\in\mathbb{R}^{N\times F_{V}}$,
where $F_{V}$ is the total number of atom types. Similarly, the bond features
are a tensor $\mathbf{X}_{E}\in\mathbb{R}^{N\times N\times F_{E}}$,
representing both the graph structure and $F_{E}$ bond types.

Let $\mathcal{C}=\\{c_{1},c_{2},\dots,c_{M}\\}$ be a set of $M$ numerical and
categorical conditions. The task is: $q(G\mid c_{1},c_{2},\dots,c_{M})\propto
q(G)q(c_{1},c_{2},\dots,c_{M}\mid G)$, where $q$ represents observed
probability. We use a model parameterized by $\theta$ for multi-conditional
molecular generation $p_{\theta}(G\mid\mathcal{C})$. The evaluation involves
both distribution learning $q(G)$ (polykovskiy2020molecular) and condition
control $q(c_{1},c_{2},\dots,c_{M}\mid G)$. We follow previous work in
assuming that there exist different oracle functions $\mathcal{O}$ that can
independently evaluate each conditioned property (gao2022sample):
$q(c_{1},c_{2},\dots,c_{M}\mid G)=\prod_{i=1}^{M}\mathcal{O}_{i}(c_{i}\mid
G)$. Note that the oracles are not used in the training of $p_{\theta}$.

###  2.2 Diffusion Model on Graph Data

Diffusion models consist of forward and reverse diffusion processes
(ho2020denoising). We refer to the forward diffusion process as the diffusion
process following (ho2020denoising). The diffusion process $q(G^{1:T}\mid
G^{0})=\prod_{t=1}^{T}q(G^{t}\mid G^{t-1})$ corrupts molecular graph data
($G^{0}=G$) into noisy states $G^{t}$. As timesteps $T\rightarrow\infty$,
$q(G^{T})$ converges to a stationary distribution $\pi(G)$. The reverse Markov
process $p_{\theta}(G^{0:T})=q(G^{T})\prod_{t=1}^{T}p_{\theta}(G^{t-1}\mid
G^{t})$ is parameterized by neural networks to gradually denoise the latent
states toward the desired data distribution.

#### Diffusion Process

One may perturb $G^{t}$ in a discrete state-space to capture the structural
properties of molecules (vignac2022digress). Two transition matrices
$\mathbf{Q}_{V}\in\mathbb{R}^{F_{V}\times F_{V}}$ and
$\mathbf{Q}_{E}\in\mathbb{R}^{F_{E}\times F_{E}}$ are defined for nodes
$\mathbf{X}_{V}$ and edges $\mathbf{X}_{E}$, respectively (vignac2022digress).
Then, each step $q(G^{t}\mid G^{t-1},G^{0})=q(G^{t}\mid G^{t-1})$ in the
diffusion process is sampled as follows.

$\left\\{\begin{aligned}
q(\mathbf{X}_{V}^{t}\mid\mathbf{X}_{V}^{t-1})&=\operatorname{Cat}\left(\mathbf{X}_{V}^{t};\mathbf{p}=\mathbf{X}^{t-1}_{V}\mathbf{Q}_{V}^{t}\right),\\\
q(\mathbf{X}_{E}^{t}\mid\mathbf{X}_{E}^{t-1})&=\operatorname{Cat}\left(\mathbf{X}_{E}^{t};\mathbf{p}=\mathbf{X}^{t-1}_{E}\mathbf{Q}_{E}^{t}\right),\end{aligned}\right.$
(1)

where $\operatorname{Cat}(\mathbf{X};\mathbf{p})$ denotes sampling from a
categorical distribution with probability $\mathbf{p}$. We remove the
subscript (V/E) when the description applies to both nodes and edges. It is
assumed that the noise $\mathbf{Q}^{i}$ ($i\leq t$) is independently applied
to $\mathbf{X}$ in each step $i$, allowing us to rewrite
$q(\mathbf{X}^{t}\mid\mathbf{X}^{t-1})$ as the probability of the initial
state
$q(\mathbf{X}^{t}\mid\mathbf{X}^{0})=\operatorname{Cat}\left(\mathbf{X}^{t};\mathbf{p}=\mathbf{X}^{0}\bar{\mathbf{Q}}^{t}\right)$,
where $\bar{\mathbf{Q}}^{t}=\prod_{i\leq t}\mathbf{Q}^{i}$.

#### Noise Scheduling

Transition matrices $\mathbf{Q}_{V}$ and $\mathbf{Q}_{E}$ control the noise
applied to atom features and bond features, respectively. vignac2022digress
defined
$\pi(G)=(\mathbf{m}_{X}\in\mathbb{R}^{F_{V}},\mathbf{m}_{E}\in\mathbb{R}^{F_{E}})$
as the marginal distributions of atom types and bond types. The transition
matrix at timestep $t$ is
$\mathbf{Q}^{t}=\alpha^{t}\mathbf{I}+(1-\alpha^{t})\mathbf{1}\mathbf{m}^{\prime}$
for atoms or bonds, where $\mathbf{m}^{\prime}$ denotes the transposed row
vector. Therefore, we have
$\bar{\mathbf{Q}}^{t}=\bar{\alpha}^{t}\mathbf{I}+(1-\bar{\alpha}^{t})\mathbf{1}\mathbf{m}^{\prime}$,
where $\bar{\alpha}^{t}=\prod_{\tau=1}^{t}\alpha^{\tau}$. The cosine schedule
(nichol2021improved) is often chosen for
$\bar{\alpha}^{t}=\cos(0.5\pi(t/T+s)/(1+s))^{2}$.

#### Reverse Process

With the initial sampling $G^{T}\sim\pi(G)$, the reverse process generates
$G^{0}$ iteratively in reversed steps $t=T,T-1,\dots,0$. We use a neural
network to predict the probability $p_{\theta}(\tilde{G}^{0}\mid G^{t})$ as
the product over nodes and edges (austin2021structured; vignac2022digress):

$p_{\theta}(\tilde{G}^{0}\mid G^{t})=\prod_{v\in V}p_{\theta}(v^{t-1}\mid
G^{t})\prod_{e\in E}p_{\theta}(e^{t-1}\mid G^{t})$ (2)

$p_{\theta}(\tilde{G}^{0}\mid G^{t})$ could be combined with $q(G^{t-1}\mid
G^{t},G^{0})$ to estimate the reverse distribution on the graph
$p_{\theta}(G^{t-1}\mid G^{t})$. For example, $p_{\theta}(v^{t-1}\mid G^{t})$
is marginalized over predictions of node types
$\tilde{v}\in\tilde{\mathbf{x}}_{v}$, which applies similarly to edges:

$p_{\theta}(v^{t-1}\mid
G^{t})=\sum_{\tilde{v}\in\tilde{\mathbf{x}}_{v}}q(v^{t-1}\mid\tilde{v},G^{t})p_{\theta}(\tilde{v}\mid
G^{t}).$ (3)

The neural network could be trained to minimize the negative log-likelihood
(vignac2022digress).

$L=\mathbb{E}_{q(G^{0})}\mathbb{E}_{q(G^{t}|G^{0})}\left[-\mathbb{E}_{\mathbf{x}\in
G^{0}}\log p_{\theta}\left(\mathbf{x}\mid G^{t}\right)\right]$ (4)

where $\mathbf{x}\in G^{0}$ denotes the node or edge features. Typically, the
reverse process in diffusion models does not consider molecular properties as
conditions. While there have been efforts to introduce property-related
guidance using additional predictors, the more promising approach of
predictor-free guidance (ho2022classifier), particularly in multi-conditional
generation, remains underexplored.

##  3 Multi-Conditional Diffusion Model

We present the denoising model with multiple predictor-free guidance in Figure
2. A condition encoder is proposed to learn the representation of $M$
conditions. The statistics of this representation like mean and variance are
used to replace the ones from the molecule representation (huang2017arbitrary)
(see Section 3.2). Besides, we introduce a new noise model in the diffusion
process to better fit graph-structured molecules.(see Section 3.1).

Figure 2: Denoising framework with multi-conditional guidance: the model
consists of the condition encoder, structure encoder, and structure decoder.
Details are in Section 3.2.

###  3.1 Graph-dependent Noise Model

In the noise model, the transition probability of a node or an edge should
rely on the joint distribution of nodes and edges in the prior state. However,
as an example shown in Eq. 1, current diffusion models (jo2022score;
vignac2022digress; lee2023exploring) treat node and edge state transitions as
independent, misaligning with the denoising process in Eq. 3. This difference
between the sampling distributions of noise in the diffusion and reverse
processes introduces unnecessary challenges to multi-conditional molecular
generations.

To address it, we use a single matrix $\mathbf{X}_{G}\in\mathbb{R}^{N\times
F_{G}}$ to represent graph features for $G$, where $F_{G}=F_{V}+N\cdot F_{E}$.
This matrix is created by concatenating the node feature matrix
$\mathbf{X}_{V}$ and the flattened edge connection matrix from the last two
dimensions of $\mathbf{X}_{E}$. Each row in $\mathbf{X}_{G}$ contains features
for both nodes and edges, representing all connections and non-connections to
other nodes. We also design a transition matrix $\mathbf{Q}_{G}$ considering
the joint distribution of nodes and edges.
$\mathbf{Q}_{G}\in\mathbb{R}^{F_{G}\times F_{G}}$ is constructed from four
matrices $\mathbf{Q}_{V},\mathbf{Q}_{EV}\in\mathbb{R}^{F_{E}\times
F_{V}},\mathbf{Q}_{E},\mathbf{Q}_{VE}\in\mathbb{R}^{F_{V}\times F_{E}}$,
denoting the transition probability ($\text{``dependent old
state''}\rightarrow\text{``target new state''})$
$\text{node}\rightarrow\text{node};\text{edge}\rightarrow\text{node};\text{edge}\rightarrow\text{edge};\text{node}\rightarrow\text{edge}$,
respectively.

$\mathbf{Q}_{G}=\begin{bmatrix}\mathbf{Q}_{V}&\mathbf{1}_{N}^{\prime}\otimes\mathbf{Q}_{VE}\\\
\mathbf{1}_{N}\otimes\mathbf{Q}_{EV}&\mathbf{1}_{N\times
N}\otimes\mathbf{Q}_{E}\end{bmatrix},$ (5)

where $\otimes$ denotes the Kronecker product, $\mathbf{1}_{N}$,
$\mathbf{1}_{N}^{\prime}$, and $\mathbf{1}_{N\times N}$ represent the column
vector, row vector, and matrix with all 1 elements, respectively. According to
Eq. 5, the first $F_{V}$ columns in $\mathbf{Q}_{G}$ determine the node
feature transitions based on both node features (first $F_{V}$ rows) and edge
features (remaining $N\cdot F_{E}$ rows). Conversely, the remaining $N\cdot
F_{E}$ columns determine the edge feature transitions, depending on the graph
features. Applying it to the graph feature $\mathbf{X}_{G}$, we have a new
diffusion noise model.

$q(\mathbf{X}_{G}^{t}\mid\mathbf{X}_{G}^{t-1})=\operatorname{Cat}\left(\mathbf{X}_{G}^{t};\mathbf{p}=\mathbf{X}^{t-1}_{G}\mathbf{Q}_{G}^{t}\right).$
(6)

#### Choice of $\mathbf{Q}_{VE}$ and $\mathbf{Q}_{EV}$

Similar to the definitions of $\mathbf{m}_{V}$ and $\mathbf{m}_{E}$
(vignac2022digress), we leverage the prior knowledge within the training data
for the formulation of task-specific matrices, $\mathbf{Q}_{EV}$ and
$\mathbf{Q}_{VE}$. Hence, we calculate co-occurrence frequencies of atom and
bond types in the training molecular graphs to derive the marginal atom-bond
co-occurrence probability distribution. For each bond type, each row in
$\mathbf{m}_{EV}$ represents the probability of co-occurring atom types.
$\mathbf{m}_{VE}$ is the transpose of $\mathbf{m}_{EV}$ and has a similar
meaning. Subsequently, we define
$\mathbf{Q}_{EV}=\bar{\alpha}^{t}\mathbf{I}+(1-\bar{\alpha}^{t})\mathbf{1}\mathbf{m}_{EV}^{\prime}$
and
$\mathbf{Q}_{VE}=\bar{\alpha}^{t}\mathbf{I}+(1-\bar{\alpha}^{t})\mathbf{1}\mathbf{m}_{VE}^{\prime}$.

###  3.2 Denoising Generation with Multiple Conditions

We introduce a denoising model to generate molecules under multi-conditions
$\mathcal{C}=\\{c_{1},c_{2},\dots,c_{M}\\}$ without predictors.

#### Predictor-free Guidance

The predictor-free reverse process $\hat{p}_{\theta}(G^{t-1}\mid
G^{t},\mathcal{C})$ aims to generate molecules with a high probability
$q(\mathcal{C}\mid{G}^{0})$. This could be achieved by a linear combination of
the log probability for unconditional and conditional denoising
(ho2022classifier):

$\begin{split}&\hat{p}_{\theta}(G^{t-1}\mid G^{t},\mathcal{C})=\log
p_{\theta}(G^{t-1}\mid G^{t})\\\ &+s\left(\log p_{\theta}(G^{t-1}\mid
G^{t},\mathcal{C})-\log p_{\theta}(G^{t-1}\mid G^{t})\right),\end{split}$ (7)

where $s$ denotes the scale of conditional guidance. Unlike classifier-free
guidance (ho2022classifier), which typically predicts noise, we directly
estimate $p_{\theta}(\tilde{G}^{0}\mid G^{t},\mathcal{C})$. We utilize a
single Graph Transformer model, $f_{\theta}(G^{t},\mathcal{C})$, for both
$p_{\theta}(\tilde{G}^{0}\mid G^{t})$ and $p_{\theta}(\tilde{G}^{0}\mid
G^{t},\mathcal{C})$. Here, $f_{\theta}(G^{t},\mathcal{C}=\emptyset)$ computes
the unconditional probability, substituting original embeddings for a
condition with a learnable dropping embedding. To learn this dropping
embedding, we randomly replace elements in $\mathcal{C}$ with a special
dropping value during training. The Transformer architecture for
$f_{\theta}(G^{t}=\mathbf{X}_{G}^{t},\mathcal{C})$ comprises three components:
the condition encoder, the structure encoder, and the structure decoder. An
overview of the architecture is presented in Figure 3.

Figure 3: The architecture for condition encoder, structure encoder, and
structure decoder.

#### Condition Encoder

We treat the timestep $t$ as a special condition and follow (nichol2021glide)
to obtain a $D$-dimensional representation $\mathbf{t}$ with sinusoidal
encoding. For other task property-related numerical or categorical condition
$c_{i}\in\mathcal{C}$, we apply distinct encoding operations to get their
$D$-dimensional representation. For a categorical condition, we use the one-
hot encoding with a $\operatorname{Linear}$ layer. For a numerical variable,
we introduce a clustering encoding method. This involves defining learnable
centroids, assigning $c_{i}$ to clusters, and transforming the soft assignment
vector of condition values into the representation. It could be implemented
using two $\operatorname{Linear}$ layers and a $\operatorname{Softmax}$ layer
in the middle as:
$\operatorname{Linear}\left(\operatorname{Softmax}\left(\operatorname{Linear}(c_{i})\right)\right)$.
Finally, we could obtain the representation of the condition as
$\mathbf{c}=\sum_{i=1}^{M}\operatorname{encode}(c_{i})$, where
$\operatorname{encode}$ is the specific encoding method based on the condition
type. For numerical conditions, we evaluate our proposed clustering-based
approach against alternatives like direct or interval-based encodings
(liu2023semi). As noted in LABEL:subsec:ablation, the clustering encoding
outperforms the other methods.

#### Structure Encoder

The graph feature is first encoded into the latent space as
$\mathbf{H}=\operatorname{Linear}(\mathbf{X}_{G}^{t})$, where
$\mathbf{H}\in\mathbb{R}^{N\times D}$. We then adapt the standard Transformer
encoder (vaswani2017attention) to design our structure encoder. To learn the
latent representation under conditions, we use an adaptive layer normalization
(denoted by $\operatorname{AdaLN}$) (huang2017arbitrary; peebles2023scalable)
to replace the standard one:
$\mathbf{H}=\operatorname{AdaLN}(\mathbf{H},\mathbf{c}$). For each row
$\mathbf{h}$ in $\mathbf{H}$:

$\operatorname{AdaLN}\left(\mathbf{h},\mathbf{c}\right)=\gamma_{\theta}(\mathbf{c})\odot\frac{\mathbf{h}-\mu\left(\mathbf{h}\right)}{\sigma\left(\mathbf{h}\right)}+\beta_{\theta}(\mathbf{c}),$
(8)

where $\mu(\cdot)$ and $\sigma(\cdot)$ are mean and variance values. $\odot$
indicates element-wise product. $\gamma_{\theta}(\cdot)$ and
$\beta_{\theta}(\cdot)$ are neural network modules in $f_{\theta}(\cdot)$,
each of which consists of two linear layers with $\operatorname{SiLU}$
activation (elfwing2018sigmoid) in the middle. We have a gated variant
$\operatorname{AdaLN}_{gate}$ for residuals:

$\operatorname{AdaLN}_{gate}\left(\mathbf{h},\mathbf{c}\right)=\alpha_{\theta}(\mathbf{c})\odot\operatorname{AdaLN}\left(\mathbf{h},\mathbf{c}\right)$
(9)

We apply the zero initialization for the first layer of
$\gamma_{\theta}(\cdot),\beta_{\theta}(\cdot)$, and $\alpha_{\theta}(\cdot)$
(peebles2023scalable). There are other options to learn the structure
representation from the condition (peebles2023scalable): $\operatorname{In-
Context}$ conditioning adds condition representation to the structure
representation at the beginning of the structure encoder, and
$\operatorname{Cross-Attention}$ calculates cross-attention between the
condition and structure representation. We observe in LABEL:subsec:ablation
that $\operatorname{AdaLN}$ performs best among them.

#### Structure Decoder

The latent representation $\mathbf{H}$ from the final structure encoder layer
is used to predict node probabilities $\tilde{\mathbf{X}}_{V}^{0}$ and edge
probabilities $\tilde{\mathbf{X}}_{E}^{0}$ at $t=0$.

$\tilde{\mathbf{X}}_{G}^{0}=\operatorname{AdaLN}(\operatorname{MLP}(\mathbf{H}),\mathbf{c}).$
(10)

We split the output $\mathbf{X}_{G}$ into atom and bond features
$\tilde{\mathbf{X}}_{V}^{0},\tilde{\mathbf{X}}_{E}^{0}$. The first $F_{V}$
dimensions of $\tilde{\mathbf{X}}_{G}^{0}$ represent node type probabilities,
and the remaining $N\cdot F_{E}$ dimensions cover probabilities for $N$ edge
types associated with the node, as detailed in Section 3.1.

#### Generation to Molecule Conversion

A common way of converting generated graphs to molecules selects only the
largest connected component (vignac2022digress), denoted as MCD-LC in our
model. For MCD, we connect all components by randomly selecting atoms. It
minimally alters the generated structure to more accurately reflect model
performance than MCD-LC.

##  4 Experiment

RQ1: We validate the generative power of MCD compared to baselines from
molecular optimization and diffusion models in LABEL:subsec:generation-main.
RQ2: We study an inverse polymer design for gas separation in
LABEL:subsec:study. RQ3: We conduct further analysis to examine MCD in
LABEL:subsec:ablation.

###  4.1 Experimental Setup

We use datasets with over ten types of atoms and up to fifty nodes in a
molecular graph. We include both numerical and categorical properties for
drugs and materials, offering a benchmark for evaluation across diverse
chemical spaces. Model performance is validated across up to nine metrics,
including distribution coverage, diversity, and condition control capacity for
various properties.

#### Datasets and Input Conditions

We have one polymer dataset (thornton2012polymer) for materials, featuring
three numerical gas permeability conditions: O2Perm, CO2Perm, and N2Perm. For
drug design, we create three class-balanced datasets from MoleculeNet
(wu2018moleculenet): HIV, BBBP, and BACE, each with a categorical property
related to HIV virus replication inhibition, blood-brain barrier permeability,
or human $\beta$-secretase 1 inhibition, respectively. We have two more
numerical conditions for synthesizability from synthetic accessibility (SAS)
and complexity scores (SCS) (ertl2009estimation; coley2018scscore).

Table 1: Generation of 10K Polymers: Results on a numerical synthesizability
score (Synth.) and three numerical properties (gas permeability for O2, N2,
CO2). MAE is calculated between input conditions and generated properties.
Best results are  highlighted.


# Equatorial Lensing in the Balasin-Grumiller Galaxy Model

Marco Galoppo1,∗, Sergio Luigi Cacciatori2,3,†, Vittorio Gorini2,‡ and
Mariarosa Mazza4,+

1School of Physical & Chemical Sciences, University of Canterbury,
Private Bag 4800, Christchurch 8041, New Zealand
2Department of Science and High Technology, Università degli Studi
dell’Insubria,
Via Valleggio 11, 22100, Como, Italy
3INFN, sezione di Milano, Via Celoria 16, 20133, Milano, Italy
4Department of Mathematics, Università degli Studi di Roma Tor Vergata,
Via della Ricerca Scientifica 1, 00133, Rome, Italy

<EMAIL_ADDRESS><EMAIL_ADDRESS>
‡vittorio.gorini<EMAIL_ADDRESS>

###### Abstract

The Balasin-Grumiller model has been the first model employed as an attempt
towards providing a fully general relativistic description of the dynamics of
a disc galaxy. In this paper, we compute the equatorial gravitational lensing
observables of the model. Indeed, our purpose is to investigate the role that
gravitational lensing plays as an observable in distinguishing between the
state-of-the-art galaxy models and the fully general relativistic ones, with
the latter stressing the role of frame-dragging and hence conceivably pointing
to a possible re-weighting of the dark matter content of disc galaxies. We
obtain for the Balasin-Grumiller model the exact formula for the bending angle
of light and we provide a corresponding estimate for the time delay between
images in the equatorial plane. For a reasonable choice for the values of the
parameters of the solution (bulge and scale radiuses, and average rotational
star speeds), the values that we obtain for the bending angle are in agreement
with those observed for typical disc galaxies. On the other hand, the
calculated time delay, which is directly tied to the frame-dragging generated
by the angular momentum of the galaxy, turns out to be some orders of
magnitude larger than the ones measured for the class of galaxies that the
Balasin-Grumiller model would claim to describe. We believe this abnormal
discrepancy to be due to the very nature of the Balasin-Grumiller model.
Namely, it being rigidly rotating, hence providing an unphysical amount of
frame-dragging. Therefore, we conclude that, in spite of its simplicity and
its unquestionable didactical value, the Balasin-Grumiller model is far too
crude to provide an instrument for a reliable general relativistic description
of a disc galaxy and that further work in the fully general relativistic
modelling of galaxies is required to reach a satisfactory stage.

Keywords— Gravitational Lensing, Galaxy Models, General Relativity.

##  1 Introduction

The description of galactic dynamics has been highly contentious in recent
years. Many approaches have been developed to cure the conflict arising with
the observed flat rotation curves when implementing an ordinary Newtonian
description of it [1, 2, 3, 4, 5, 6]. The most common approaches are MOND [7,
8, 9, 10, 11, 12, 13, 14, 15], the MOG theories [16], in which new invariants
are added to the Einstein-Hilbert action, and the Dark Matter (DM) hypothesis
[17]. This hypothesis, namely the assumption of the existence of a non-
baryonic component of mass dominating the matter density in the universe, is
one of the foundations of the widely accepted $\Lambda$CDM cosmological model
(e.g. [18]) and, at the same time, one of the greatest mysteries in the
universe. Indeed, besides justifying the rotation curves of disc galaxies and
the velocity distribution of galaxies in galaxy clusters(GCs) [1, 2, 3, 4, 5,
6, 19], the DM hypothesis has been incredibly successful in interpreting
several other astrophysical and cosmological observables such as the
thermodynamical properties of X-ray emitting gas in GCs, the gravitational
lensing produced by their mass distributions, the features of the two Bullet
Clusters, and the growth of cosmic structures from inhomogeneities in the
matter density content of the early universe [20, 21, 22, 23, 24, 25, 26].
However, despite such remarkable successes, the results of the many
enterprises aimed at the direct detection of DM particles are, to date,
inconclusive [27, 28, 29, 30, 31, 32, 33]. In addition, over the years,
several observations of the universe have, to a larger or lesser extent,
challenged the validity of the $\Lambda$CDM model [34, 35, 36, 37]. In
particular, concerning galaxies, some recent independent observations in
different domains, among which the unexpected, apparent formation of massive
galaxies at very early times spotted by JWST, appear to be in tension with the
standard DM paradigm111Nonetheless, it is fair to say that, concerning the
mentioned ancient galaxies observed by JWST, several explanations have been
put forward, both theoretical as well as based on observations, which
indicates that $\Lambda$CDM might be resilient after all [38, 39, 40, 41][42,
43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59]. Given the
above difficulties confronting the DM hypothesis as well as those affecting
MOND and MOG theories [60], a new approach has recently gained traction. This
approach, which so far has been advocated for the study of the features of
galaxies alone, suggests the full employment of General Relativity (GR),
namely of Einstein equations, to describe galactic dynamics [61, 62, 63, 64,
65, 66, 67, 68, 69]

The core idea behind the proposed approach of fully general relativistic
modelling of galaxies is that the highly nonlinear dynamics of spacetime in GR
carries energy with it, which, in turn, gravitates, thus helping to bind the
galaxies together. Usually, full GR is not considered a viable solution to
galactic dynamics since the velocities of stars and gas in a galaxy are much
smaller than the speed of light, and the gravitational field is considered
”weak” far from the central region. Hence, linearized GR, in particular the
Newtonian limit, is implemented almost universally in the study of galactic
dynamics. Moreover, a recent remarkable paper by Ciotti [70] proves that even
a perturbative implementation of a gravitomagnetic limit would not be a
feasible solution. Thus, if one believes, as we are inclined to do, that the
solution to the issue is to be found in a better implementation of GR, this
solution should lie in a non-perturbative approach to the problem. Indeed,
though in the presence of low velocities and weak gravitational fields, the
linearized approximation of GR is certainly valid everywhere locally, it may
not to be valid anymore globally in spatially extended rotating systems, such
as galaxies. To wit, we expect the angular momentum of a concrete galaxy
featuring a few hundred billion stars in a disc with a radius of tens of kpc
to be so high as to exert such a considerable frame-dragging as to give rise
in the metric to off-diagonal elements of the same order of magnitude of the
diagonal ones.

In this spirit, after the pioneering work carried out in [61, 62], the first
viable full GR model for a galaxy was described by Balasin and Grumiller in
[64]. The Balasin-Grumiller model (BG) eliminated the unphysical behaviour of
the previous solutions and was claimed by the authors to reduce the need for
DM by around 30%. Moreover, in a recent study, Crosta et al [65] have shown
that, after combining the BG model for a portion of the Milky Way with better
modellings of the matter distribution inside the galaxy, one is able to
correctly fit the data from the Gaia DR2 catalogue as well as other Newtonian
models, but with a much smaller set of free parameters. In particular, the
authors state that the rotation curves measured by Gaia agree with the flat BG
prevision. In [70], Ciotti shows that if the Newtonian calculation is
rigorously performed for a thin disc of matter, it too predicts a flat
rotation curve for up to three scale galaxy radii. Hence, Ciotti claims
Crosta’s argument to be at best fallacious, since the Gaia data implemented in
Crosta’s study alone cannot tell which model is the best. Nevertheless,
Newtonian models require a much larger number of free parameters in order to
correctly fit the Gaia data. Furthermore, in a recent communication [66] it
has been anticipated that, based on the recent Gaia DR3 data release, Crosta’s
results seem to be confirmed well beyond the validity of the Newtonian
approximation. Nonetheless, in view of the limitations affecting the BG model,
in [67], the class of full GR models to which BG belongs, namely the most
general class of stationary, axisymmetric dust solutions of Einstein equations
(see [71, 72]), has been further studied as a much more flexible setting than
BG for the behaviour of disc galaxies. The merit of BG lies in being the first
non-perturbative GR model applied to the description of disc galaxies.
Furthermore, it has been the only model of the class attempted to be tested so
far against observational data. As such, it is a good starting point for a
study of gravitational lensing in the class of full GR galaxy models. However,
in spite of its undoubtful pedagogical value and its simplicity, we are fully
aware that BG is a far too crude model to be sufficiently viable for the
description of the dynamics of a disc galaxy, particularly since it implies
unphysical rigid rotation.

Gravitational lensing by galaxies has been employed to study the proposed DM
haloes around galaxies for many years now [73, 74, 75, 76, 77, 78, 79].
Furthermore, microlensing in GCs, and in galaxies, has been utilized to
investigate the distribution of DM in these systems on different scales [80,
81, 82, 83]. Overall, gravitational lensing has been acknowledged as a
powerful tool in the investigation of the DM mystery. Hence, it only sounds
reasonable to compare observed lensing parameters (e.g. deflection angles and
time delays) in disc galaxies to values of the same quantities computed from
full GR galaxy models as tests of the extent of the validity of the model
themselves.

In this paper, we study the equatorial lensing in BG. We fully characterize
the deflection angle and the time delay between the two images and obtain
quantitative results to compare with the observational data. In section 2, we
introduce the mathematical tools necessary to calculate the magnitude of the
bending angle, especially focusing on the construction of the lensing
geometry. In section 3, we define the general class of galaxy models, the
class of observers through which we read their physics and we specialize to
BG. In section 4, we obtain the general formula for the equatorial deflection
angle in BG and compare our predictions for a Milky Way-like galaxy with the
observational data. In section 5, we study the time delay between light rays
following retrograde and prograde motion on the equatorial plane and obtain
quantitative results for a Milky Way-like galaxy. Section 6 is dedicated to a
brief overview of the results and the discussion of future perspectives. In
Appendix A, we discuss the existence of planar motion on the equatorial plane
for photons. Appendix B gives a full account of the numerical methods applied
to obtain the results in this paper.

##  2 Mathematical Background

In this section, we report the mathematical material related to the
calculation of the bending angle of light, which deviates from the standard
techniques. The reader who is not strictly interested in the technical details
can skip this section. Consider a generic spacetime metric

$ds^{2}=g_{00}dt^{2}+2g_{0i}dtdx^{i}+g_{ij}dx^{i}dx^{j}.$ (1)

From this, we construct the generalized optical metric[84] defined by

$\gamma_{ij}\coloneqq-\frac{g_{ij}}{g_{00}}+\frac{g_{0i}g_{0j}}{g_{00}^{2}}.$
(2)

Moreover, we define the shift vectors

$\beta_{i}\coloneqq-\frac{g_{0i}}{g_{00}}.$ (3)

Then, if we consider the motion of a photon

$ds^{2}=0,$ (4)

we can combine (2),(3) and (4) to obtain

$dt=\pm\sqrt{\gamma_{ij}dx^{i}dx^{j}}+\beta_{i}dx^{i},$ (5)

where the sign depends on the specific orientation of the photon. In the case
of stationary, axisymmetric metrics (e.g. the Kerr metric), the $\pm$ sign
indicates prograde and retrograde motions respectively. We also consider

$d\ell=\sqrt{\gamma_{ij}dx^{i}dx^{j}},$ (6)

the arc length in the three-geometry defined by (2). Thus, if the metric under
examination is stationary, from (5) and (6) we can define

$t\coloneqq\int
dt=\int\left(\sqrt{\gamma_{ij}\dot{x}^{i}\dot{x}^{j}}+\beta_{i}\dot{x}^{i}\right)d\ell,$
(7)

where the dot indicates differentiation with respect to $\ell$. We can now use
the generalized Fermat principle [85, 86] to obtain the light’s paths
equations starting from (7). We get

$\displaystyle\ddot{x}^{k}$
$\displaystyle=-{}^{(3)}\Gamma^{k}_{mn}\dot{x}^{m}\dot{x}^{n}+\gamma^{ki}\left(\beta_{m;i}-\beta_{i;m}\right)\dot{x}^{m}$
(8)

where ${}^{(3)}\Gamma^{k}_{mn}$ are the Christoffel symbols of the three-
geometry associated to $\gamma_{ij}$, and $\gamma^{ij}$ are the components of
the inverse three-spatial metric.

Usually, at this point, the standard approach would entail some sort of
integration of the path equations. However, in the specific problem we are
interested in, there is a more suitable approach. This method was first
proposed in [87] and later on generalized in [88, 89, 84]. The strategy rests
on the use of the Gauss-Bonnet Theorem (GBT) which can be summarized through
the formula

$\sum_{i=1}^{n}\int_{C_{i}}k_{g}(\ell)d\ell+\int\int_{R}Kd\sigma+\sum_{j=1}^{p}\theta_{j}=2\pi\chi(R),$
(9)

where $R$ identifies a portion of a two-dimensional surface bounded by the
union of the curves $C_{i}$, $k_{g}$ represents the geodesic curvature along
the boundary of $R$, $K$ is the Gaussian curvature of $R$, $\theta_{j}$ are
the external angles at the intersection points of the boundary curves and
$\chi(R)$ is the Euler characteristic of $R$. In Fig. 1 we report an example
of the geometry considered in the theorem.

Figure 1: Schematic figure for the Gauss-Bonnet Theorem.

To use the GBT in the setting of planar gravitational lensing we start by
defining the angle $\alpha_{\Psi}$ [87] as

$\alpha_{\Psi}\coloneqq\Psi_{O}+\Psi_{L}-\Psi_{S},$ (10)

where $\Psi_{O}$, $\Psi_{L}$ and $\Psi_{S}$ are defined as described in the
lensing geometry reported in Fig. 2, in which $S$ and $O$ indicate the
observer and the source, considered point-like and taken in a near flat region
of spacetime, $L$ indicates the centre of the lens and the curved line joining
$S$ and $O$ is the path of the light ray.

Figure 2: Gravitational lensing geometry: ${}^{S}\nabla^{O}_{L}$.

We apply the GBT to the embedded triangle of Fig. 2 and obtain

$\alpha_{\Psi}=\int_{\partial^{S}\nabla^{O}_{L}}k_{g}(\ell)d\ell+\int\int_{{}^{S}\nabla^{O}_{L}}Kd\sigma.$
(11)

(11) gives us a coordinate-independent way to calculate $\alpha_{\Psi}$.
However, from (10), it is clear that in the case of lenses with strong
gravitational fields, the angle $\Psi_{L}$, and therefore $\alpha_{\Psi}$,
might be ill-defined. Hence, we introduce the angle $\alpha$ [88] defined
through

$\alpha\coloneqq\Psi_{O}+\phi_{OS}-\Psi_{S},$ (12)

where $\phi_{OS}=\phi_{O}-\phi_{S}$ is the coordinate angular separation
between $O$ and $S$. The angle defined in (12) is what we identify as the
bending angle of light since it has been proved to take on the same value as
the canonically defined bending angle in cases in which the latter is known
[87, 88, 89]. We want to move from (12) to an equation for $\alpha$ similar to
(11) for $\alpha_{\Psi}$. To do so, we start by considering the embedded
triangle shown in Fig. 3.

Figure 3: Embedded triangle: ${}^{\infty}\nabla_{L}^{\infty}$. The arc is
taken at spatial infinity.

By means of the GBT applied to the geometry of
${}^{\infty}\nabla_{L}^{\infty}$ we express $\Psi_{L}$ in terms of geometrical
quantities as

$\Psi_{L}=\int_{\partial^{\infty}\nabla_{L}^{\infty}}k_{g}(\ell)d\ell+\int\int_{{}^{\infty}\nabla_{L}^{\infty}}Kd\sigma.$
(13)

Now, since we are considering an asymptotically flat spacetime, and the arc is
at flat infinity, we have that $k_{g}(\ell)\longrightarrow 1/r$ and
$d\ell\longrightarrow rd\phi$. Therefore

$\int_{C_{\infty}}k_{g}(\ell)d\ell=\phi_{OS},$ (14)

and

$\Psi_{L}=\phi_{OS}+\int_{L\left(\uparrow
O\right)}^{\infty}k_{g}(\ell)d\ell+\int_{\infty\left(\downarrow
S\right)}^{L}k_{g}(\ell)d\ell+\int\int_{{}^{\infty}\nabla_{L}^{\infty}}Kd\sigma,$
(15)

where $\left(\uparrow O\right)$ and $\left(\downarrow S\right)$ are used to
underline the path of integration. (15) allows us to express $\phi_{OS}$ in
terms of $\Psi_{L}$ and thus $\alpha$ through $\alpha_{\Psi}$ as

$\alpha=\alpha_{\Psi}-\int_{L\left(\uparrow
O\right)}^{\infty}k_{g}(\ell)d\ell-\int_{\infty\left(\downarrow
S\right)}^{L}k_{g}(\ell)d\ell-\int\int_{{}^{\infty}\nabla_{L}^{\infty}}Kd\sigma.$
(16)

We substitute (11) in (16) and impose the radial lines from the lens to the
source and the observer to be radial geodesics (see Appendix A for the proof
of the validity of this assumption in the case of BG) to obtain the core
formula

$\alpha=\int_{S}^{O}k_{g}(\ell)d\ell-\int\int_{\prescript{\infty}{O}{\square}^{\infty}_{S}}Kd\sigma,$
(17)

where the geometry of $\prescript{\infty}{O}{\square}^{\infty}_{S}$ is showed
in Fig. 4. Finally, we want to highlight how the chosen strategy for the
bending angle calculation allows us to consider an extended lens and does not
force the thin lens approximation.

Figure 4: $\prescript{\infty}{O}{\square}^{\infty}_{S}$ geometry.

##  3 Galaxy Models

BG belongs to the class of galaxy models which attempts at describing disc
galaxies through a stationary, axisymmetric metric expressed in standard
cylindrical coordinates

$\displaystyle ds^{2}=$ $\displaystyle
g_{tt}(r,z)dt^{2}+2g_{t\phi}(r,z)dtd\phi$
$\displaystyle+g_{\phi\phi}(r,z)d\phi^{2}+e^{\mu(r,z)}\left(dr^{2}+dz^{2}\right),$
(18)

where we use the convention c = 1. The metric in (18) is coupled to a dust
energy-momentum tensor of the form

$T_{\mu\nu}=\rho(r,z)u_{\mu}u_{\nu}.$ (19)

The coupling has been worked out in [72, 71]. Thus, we have

$u^{\mu}\partial_{\mu}=\sqrt{-H}\left(\partial_{t}+\Omega\partial_{\phi}\right),$
(20) $\displaystyle
g_{tt}=\frac{\left(H-\eta\Omega\right)^{2}-r^{2}\Omega^{2}}{H},$ (21)
$\displaystyle g_{t\phi}=\frac{r^{2}-\eta^{2}}{H}\ \Omega+\eta,$ (22)
$\displaystyle g_{\phi\phi}=\frac{\eta^{2}-r^{2}}{H},$ (23)
$\displaystyle\mu_{,r}=\frac{1}{2r}\left[g_{tt,r}g_{\phi\phi,r}-g_{tt,z}g_{\phi\phi,z}-\left(g_{t\phi,r}\right)^{2}+\left(g_{t\phi,z}\right)^{2}\right],$
(24)
$\displaystyle\mu_{,z}=\frac{1}{2r}\left[g_{tt,z}g_{\phi\phi,r}-g_{tt,r}g_{\phi\phi,z}-2g_{t\phi,z}g_{t\phi,r}\right],$
(25)

and

$8\pi
G\rho=\frac{\left(\eta_{,r}^{2}+\eta_{,z}^{2}\right)\left[\eta^{2}r^{-2}\left(2-\ell\eta\right)^{2}-r^{2}\ell^{2}\right]}{4\eta^{2}e^{\mu}},$
(26)

where $\eta$ is a function of $r$ and $z$, $H$ is an arbitrary negative
function of $\eta$, $\ell=H^{\prime}/H$ is the logarithmic derivative of $H$
and $\Omega$ is defined as

$\Omega\coloneqq\frac{1}{2}\int\frac{H^{\prime}}{\eta}d\eta.$ (27)

The parameter $\Omega(r,z)$ describes the angular velocity of the dust
referred to the coordinates in use

$\Omega=\frac{d\phi}{dt}$ (28)

The function $\eta(r,z)$ can be implicitly retrieved through

$\mathcal{F}=2\eta+r^{2}\int\ell(\eta)\frac{d\eta}{\eta}-\int\ell(\eta)\eta
d\eta,$ (29)

where, as a consequence of Einstein’s equations, $\mathcal{F}$ satisfies the
harmonic equation

$\mathcal{F}_{,rr}-\frac{1}{r}\mathcal{F}_{,r}+\mathcal{F}_{,zz}=0.$ (30)

From here on out, we refer to this class of galaxy models as the $(\eta,H)$
class.

###  3.1 The ZAMO observers

The physical meaning behind the field $\eta(r,z)$ manifests itself once we
choose the appropriate class of observers for the galaxy. In the case of
stationary, axisymmetric metrics, there exists a natural class of observers in
terms of which to read the physics of the system: the Zero Angular Momentum
Observers (ZAMO) [90, 91]. These are defined by the tetrad

$\displaystyle\bm{e^{0}}=\frac{r}{\sqrt{g_{\phi\phi}}}dt,$ (31)
$\displaystyle\bm{e^{1}}=e^{\mu/2}dr,$ (32)
$\displaystyle\bm{e^{2}}=e^{\mu/2}dz,$ (33)
$\displaystyle\bm{e^{3}}=\sqrt{g_{\phi\phi}}\left(d\phi-\chi dt\right).$ (34)

We define the velocity of the dust in the galaxy as measured by the reference
frame formed by the ZAMO, $v(r,z)$, through

$-e^{0}_{\mu}u^{\mu}\eqqcolon\frac{1}{\sqrt{1-v^{2}}}$ (35)

where $u^{\mu}$ is the four-velocity of the dust. On the other hand, we also
have

$-e^{0}_{\mu}u^{\mu}=\frac{\sqrt{-H}r}{\sqrt{g_{\phi\phi}}}=\frac{1}{\sqrt{1-\left(\eta/r\right)^{2}}}.$
(36)

Thus, the choice of ZAMO as observers and the use of (35) coupled to (36)
allows us to identify $\eta(r,z)$ as

$\eta(r,z)=rv(r,z).$ (37)

Therefore, the field $\eta(r,z)$ is henceforth identified as the product of
the velocity field of the dust, measured in the reference frame built by ZAMO,
times the radial coordinate. Now, using (29) and (37), it is convenient to
rewrite the constraint (30) directly as a functional equation for the velocity
field $v(r,z)$ (see Appendix C of [67]) which uniquely determines $v(r,z)$
once $v(r,0)$ is arbitrarily assigned.

###  3.2 Balasin-Grumiller Galaxy

To get BG we must make the mutually inclusive choices

$\displaystyle H(\eta)=-1,$ (38) $\displaystyle\Omega(r,z)=0,$ (39)

which imply that BG undergoes rigid rotation, a clear drawback of the model.
Indeed, to introduce some degree of differential rotation, we would have to
specify, through an explicit choice, a non-constant $H(\eta)$. From
(21),(22),(23),(38) and (39), we have

$\displaystyle g_{tt}=-1,$ (40) $\displaystyle g_{t\phi}=\eta,$ (41)
$\displaystyle g_{\phi\phi}=r^{2}-\eta^{2},$ (42)

Moreover, we also have to impose reflection symmetry with respect to the
equatorial plane

$\partial_{z}g_{\mu\nu\big{|}z=0}=0.$ (43)

Finally, through an appropriate choice of boundary conditions, the metric has
been fully characterized in [64] by obtaining an analytic expression for
$\eta(r,z)$

$\displaystyle\eta(r,z)=$
$\displaystyle\frac{V_{0}}{2}\sum_{\pm}\left(\sqrt{(z\pm
r_{0})^{2}+r^{2}}-\sqrt{(z\pm R)^{2}+r^{2}}\right)$
$\displaystyle+V_{0}(R-r_{0}),$ (44)

where $V_{0}$ is the asymptotic velocity measured by ZAMO, $R$ is the galaxy
radius and $r_{0}$ is the bulge radius. The other relevant quantities of the
system are obtained through the equations previously discussed for the entire
class of models.

BG achieves a good reconstruction of rotation curves of disc galaxies and, at
the same time, a reduction in DM of approximately 30% [64]. We also highlight
that a recent work [67] provides insight into the possibility that non-rigidly
rotating models could do much better in the latter regard. However, BG is the
simplest model of the $(\eta,H)$ class and the only one which has been solved
entirely. Hence, it is the one we have chosen to investigate the effect of
using a full-GR model for gravitational lensing.

##  4  Equatorial Deflection Angle in the BG Lens

We want to use (17) to calculate the bending angle of light for equatorial
lensing in BG. Thus, we have to express the integrands in terms of
$\gamma_{ij}$ and $\beta_{i}$. We start by calculating the generalized optical
metric and the shift vector for BG. From (2), (3), (18), (40), (41) and (42)
we get

$\gamma_{ij}=\begin{bmatrix}e^{\mu}&0&0\\\ 0&e^{\mu}&0\\\
0&0&r^{2}\end{bmatrix},$ (45) $\displaystyle\beta_{r}=\beta_{z}=0,$ (46)
$\displaystyle\beta_{\phi}=\eta(r,z)=rv(r,z).$ (47)

We can now go through the calculations for the geodesic curvature and the
Gaussian curvature. We start with the latter. We have [92, 93]

$K=\frac{{}^{(3)}R_{r\phi r\phi}}{\gamma},$ (48)

where ${}^{(3)}R_{ijkl}$ are the components of the Riemann tensor associated
with the spatial three-geometry identified by the generalized optical metric
and $\gamma$ is the determinant of $\gamma_{ij}$. As showed in [88], (48) can
be written as

$K=\frac{1}{\sqrt{\gamma}}\left[\partial_{\phi}\left(\frac{\sqrt{\gamma}}{\gamma_{rr}}~{}^{(3)}\Gamma^{\phi}_{rr}\right)-\partial_{r}\left(\frac{\sqrt{\gamma}}{\gamma_{rr}}~{}^{(3)}\Gamma^{\phi}_{r\phi}\right)\right],$
(49)

where ${}^{(3)}\Gamma^{k}_{ij}$ are the Christoffel symbols related to
$\gamma_{ij}$. (49) for BG reduces to

$K=-\frac{1}{\sqrt{\gamma}}\partial_{r}\left(\sqrt{\gamma_{zz}\gamma_{\phi\phi}}~{}^{(3)}\Gamma^{\phi}_{r\phi}\right).$
(50)

We also have

${}^{(3)}\Gamma^{\phi}_{r\phi}=\frac{1}{2\gamma_{\phi\phi}}\partial_{r}\gamma_{\phi\phi},$
(51)

so that

$\displaystyle
K=-\frac{1}{\sqrt{\gamma}}\left(\partial_{r}\sqrt{\gamma_{zz}}\partial_{r}\sqrt{\gamma_{\phi\phi}}+\sqrt{\gamma_{zz}}\partial^{2}_{r}\sqrt{\gamma_{\phi\phi}}\right).$
(52)

Finally, from (45) and (52) we have

$K(r,z)=-\frac{\mu_{,r}(r,z)}{r}.$ (53)

The geodesic curvature is the acceleration of the curve projected along the
vector product between the unit tangent vector to the curve and the normal to
the plane of motion [92, 93]. Hence, we can write

$k_{g}=\epsilon_{ijk}a^{i}e^{j}N^{k},$ (54)

where $a^{i}$ is the acceleration, $N^{k}$ is the normal vector to the
equatorial plane, $e^{j}$ is the unit vector tangent to the light ray’s path
and we define $\epsilon_{ijk}=\sqrt{\gamma}\varepsilon_{ijk}$ where
$\varepsilon_{ijk}$ is the Levi-Civita symbol. From (8), we have

$\displaystyle a^{i}$
$\displaystyle=\gamma^{ia}\left(\beta_{a,b}-\beta_{b,a}\right)e^{b}$
$\displaystyle=\gamma^{ia}\epsilon_{cab}\epsilon^{clm}\beta_{l;m}e^{b}.$ (55)

We substitute (55) in (54) to get

$\displaystyle k_{g}$
$\displaystyle=\epsilon_{ijk}\gamma^{ia}\epsilon_{cab}\epsilon^{clm}\beta_{l;m}e^{b}e^{j}N^{k}$
$\displaystyle=\epsilon^{lmc}\beta_{l;m}e^{b}e_{j}N_{k}\left(\delta^{j}_{b}\delta^{k}_{c}-\delta^{j}_{c}\delta^{k}_{b}\right)$
$\displaystyle=\epsilon^{lmk}\beta_{l;m}N_{k}.$ (56)

$N_{k}$ is then easily found for the equatorial plane to be

$N_{k}=\sqrt{\gamma_{zz}}\delta^{z}_{k}.$ (57)

So that, from (56) and (57) we get

$\displaystyle k_{g}$
$\displaystyle=\sqrt{\gamma_{zz}}\epsilon^{lmz}\beta_{l;m}$
$\displaystyle=\frac{\sqrt{\gamma_{zz}}}{\sqrt{\gamma}}\left(\beta_{\phi;r}-\beta_{r;\phi}\right)$
$\displaystyle=\frac{1}{\sqrt{\gamma_{\phi\phi}\gamma_{rr}}}\beta_{\phi,r}.$
(58)

(58) reduces for BG to

$k_{g}=\frac{e^{\mu/2}}{r}\eta_{,r}.$ (59)

(59) and (53) are the integrands in (17) for the BG. Nonetheless, we still
have to define the differentials and the boundaries of integration. To begin
with, the area differential in the double integral is just

$d\sigma=\sqrt{\gamma}drd\phi=re^{\mu}drd\phi,$ (60)

where $\gamma$ is evaluated on the equatorial plane. On the other hand, the
path differential is given by

$\displaystyle d\ell$ $\displaystyle=\sqrt{\gamma_{ij}dx^{i}dx^{j}}$
$\displaystyle=d\phi\sqrt{\gamma_{rr}r^{\prime}(\phi)^{2}+\gamma_{\phi\phi}}$
$\displaystyle=d\phi\sqrt{e^{\mu}r^{\prime}(\phi)^{2}+r(\phi)^{2}},$ (61)

where we have chosen the plus sign in the second line, $r(\phi)$ indicates the
trajectory of the light ray and $r^{\prime}(\phi)=dr/d\phi$. Finally, from the
geometry in Fig. 4 we can deduce the boundary of integration and, for BG, cast
(17) in the form

$\displaystyle\alpha=$
$\displaystyle\int_{\phi_{S}}^{\phi_{O}}\frac{e^{\mu/2}}{r}\eta_{,r}\sqrt{e^{\mu}r^{\prime
2}+r^{2}}d\phi+\int_{\phi_{S}}^{\phi_{O}}d\phi\int_{r(\phi)}^{\infty}\mu_{,r}e^{\mu}dr$
(62)

where $\phi_{S}$ and $\phi_{O}$ are the angular coordinates of the source and
the observer and in the first integral we must read $r$ as $r(\phi)$. We
highlight that (62) has the upside to have built-in finite-distance
corrections due to its dependence on $\phi_{S}$ and $\phi_{O}$. From the four-
dimensional equation of motion for a photon in a stationary, axisymmetric
spacetime we obtain the equation for the light ray trajectory on the
equatorial plane [84]

$\left(\frac{dr}{d\phi}\right)^{2}=\frac{-g_{tt}g_{\phi\phi}+g_{t\phi}^{2}}{g_{rr}}\cdot\frac{g_{\phi\phi}+2g_{t\phi}b+g_{tt}b^{2}}{\left[g_{tt}b+g_{t\phi}\right]^{2}},$
(63)

which, for BG, becomes

$\left(\frac{dr}{d\phi}\right)^{2}=\frac{r^{2}e^{-\mu}}{\left(b-\eta\right)^{2}}\left[r^{2}-\left(b-\eta\right)^{2}\right],$
(64)

where $b$ is the impact parameter of the light ray, defined as

$b=\frac{L}{E}=\frac{g_{t\phi}+g_{\phi\phi}\left(\frac{d\phi}{dt}\right)}{-g_{tt}-g_{t\phi}\left(\frac{d\phi}{dt}\right)}.$
(65)

In (65), $L$ and $E$ are, respectively, the constants of motion corresponding
to axial symmetry and stationarity, and are consequently interpreted as the
angular momentum and the energy of the photon as measured by a Minkowskian
observer at infinity. (64) is not analytically solvable. However, we underline
how the parameter $\bar{V}_{0}=V_{0}/c$ in BG can play the role of an
expansion parameter for a perturbative approach to the problem. Thus, we
choose to carry out the calculation for $\alpha$ in (62) to the first order in
$\bar{V}_{0}$. From (44), (24), (25), and (53) we conclude that

$K\propto\bar{V}_{0}^{2},$ (66)

and thus that the second integral in the r.h.s. of (62) is negligible to the
first order in the perturbative expansion. (59) and (44) show that the
geodesic curvature is of the first order in $\bar{V}_{0}$. Hence, $\alpha$ is,
at the leading order, of the first order in $\bar{V}_{0}$. To progress in the
calculation we must therefore solve for $r(\phi)$ in (64) to the zeroth order
in $\bar{V}_{0}$. This is equivalent to considering the unperturbed light ray
path for $r(\phi)$. We thus have

$r(\phi)=\frac{b}{\cos(\phi)}.$ (67)

We can then substitute (67) and (44) in (62) and neglect the remaining terms
above the leading order to get

$\displaystyle\alpha$
$\displaystyle=\bar{V}_{0}\int_{\phi_{S}}^{\phi_{O}}\frac{b}{\cos^{2}(\phi)}\left(\frac{1}{\sqrt{r_{0}^{2}+b^{2}/\cos^{2}(\phi)}}-\frac{1}{\sqrt{R^{2}+b^{2}/\cos^{2}(\phi)}}\right)d\phi$
$\displaystyle=\bar{V}_{0}\int_{\phi_{S}}^{\phi_{O}}\frac{\left(B^{2}-A^{2}\right)|\cos(\phi)|}{\sqrt{\left(1+A^{2}\cos^{2}(\phi)\right)\left(1+B^{2}\cos^{2}(\phi)\right)}}\left[\sqrt{1+A^{2}\cos^{2}(\phi)}+\sqrt{1+B^{2}\cos^{2}(\phi)}\right]^{-1}d\phi,$
(68)

where $A=r_{0}/b$ and $B=R/b$. (68) allows us to calculate, to the leading
order in $\bar{V}_{0}$, the deflection angle of light in BG once we have
specified the parameters of the model and the position of source and observer.
In particular, a good approximation to (68) can be obtained if we are
interested in impact parameters for which the condition $B^{2}\ll 1$ is
fulfilled. In this case, we have

$\displaystyle\alpha$
$\displaystyle=\bar{V}_{0}\int_{\phi_{S}}^{\phi_{O}}\frac{1}{|\cos(\phi)|}\left(\frac{1}{\sqrt{1+A^{2}\cos^{2}(\phi)}}-\frac{1}{\sqrt{1+B^{2}\cos^{2}(\phi)}}\right)d\phi\approx\frac{1}{2}\bar{V}_{0}\left(B^{2}-A^{2}\right)\int_{\phi_{S}}^{\phi_{O}}|\cos(\phi)|d\phi.$
(69)

In Fig. 5, we show $\alpha$ as a function of the impact parameter for
scattering outside the galaxy. The plot shows both the exact and approximated
values for $\alpha$ (see Appendix B for the details of the numerical
calculation). The values of the parameters of BG have been taken as the ones
obtained through the likelihood fit to the Gaia DR2 catalogue data in [65].
Moreover, since in the case of the typical cosmological distance between
observer, lens and source for galaxy lensing, the finite distance corrections
have been found to be of the order of $10^{-8}$ arcseconds, we considered
source and observer to be at infinite distance.

Figure 5: Exact and approximated $\alpha$ values at leading order for the
Milky Way as a lens in BG. $R=47.87$ kpc, $r_{0}=0.39$ kpc and ${V}_{0}=263.1$
km/s. Source and observer are taken at infinite distances.

Fig. 5 shows that the model produces deflection angles which match the
observed ones [94, 95, 96, 97]. The standard models used for gravitational
lensing for galaxies predict values of the order of some arcseconds for impact
parameters up to a few hundred kpc and which then rapidly fall down well below
the arcsec before reaching the Mpc milestone [98, 99]. From Fig. 5, we see
that the same behaviour is obtained for the BG model. Indeed, these same
values are those actually observed. Therefore, we conclude that the BG model
predictions match the current observations as to what concerns the deflection
angles. To a certain degree, this result could have been predicted. Indeed,
even though BG reduces the need for DM, the reduction is not so significant as
to greatly move the estimated galaxy mass from the one predicted with current
models [100, 101, 102, 103]. Therefore, as the masses predicted by current DM
galaxy models and BG are in the same ballpark, so should the deflection angles
be.

##  5 Time Delay for the BG LENS

The time delay between the light ray undergoing prograde and, respectively,
retrograde motion is another observable, besides the deflection angle, which
can shed light on the different behaviour of the effects of gravitational
lensing among the standard models, BG and observations. In particular, if the
time delay is calculated for photons which move along two geodesics
characterised by the same impact parameter (in modulus), this observable
proves an optimal probe to gauge the impact of the dragging generated by the
galaxy on the physics of the system. Indeed, the time delay, in this case, can
be attributed almost exclusively to the dragging, as the paths of the light
rays are bound to be approximately symmetrical with respect to the geometry of
the scattering. Thus, we are interested in getting an estimate of

$\Delta t(|b|)=\int_{\gamma_{+|b|}}dt-\int_{\gamma_{-|b|}}dt,$ (70)

where $b$ is the impact parameter, and the integrals are carried over along
$\gamma_{+|b|}$ and $\gamma_{-|b|}$, respectively, the path of the co-rotating
photon and the path of the counter-rotating one. There is a caveat in using
(70). Indeed, with the current coordinates implemented, we get

$\lim_{r\rightarrow\infty}g_{t,\phi}(r,0)=V_{0}(R-r_{0})\equiv\eta_{\infty},$
(71)

which means that the clocks of the observers at infinity are desynchronised
[104]. Therefore, the delay measured by (70) includes the desynchronisation
carried out by the coordinate observers at infinity. To avoid this, we need to
resynchronize the clocks at infinity through the coordinate change

$t\longrightarrow t+\eta_{\infty}\phi.$ (72)

(72) is a well-defined coordinate change insofar that $\phi$ is well-defined
on a 2$\pi$ range and, if so, it is equivalent to the shift
$\eta\longrightarrow\eta-\eta_{\infty}:=\tilde{\eta}$ for the BG metric. Thus,
the observers who use the new coordinates have synchronised clocks at
infinities and are reliable for studying the time delay. At this point, it is
crucial to underline the importance, in GR, to identify correctly the
reference frame one is adopting to compare the theoretical results with those
of the observations. A discussion on possible reference frames for
observations is expounded in [105]. They do not coincide with the observers we
are considering here, and we will not delve into this topic in the present
paper, leaving such an important task for future work including more realistic
models for galaxies. To study (70), we start by writing the equations of
motion, relative to an affine parameter, on the equatorial plane in the 4D BG
metric in which we have substituted $\tilde{\eta}$ in place of $\eta$

$\displaystyle\ddot{t}=\tilde{\eta}_{,r}\dot{r}\dot{\phi}-\frac{\tilde{\eta}\tilde{\eta}_{,r}}{r^{2}}\dot{r}\left(\dot{t}-\tilde{\eta}\dot{\phi}\right)-2\frac{\tilde{\eta}}{r}\dot{r}\dot{\phi},$
(73)
$\displaystyle\ddot{\phi}=-\frac{\tilde{\eta}_{,r}}{r^{2}}\dot{r}\left(\dot{t}-\tilde{\eta}\dot{\phi}\right)-2\frac{\dot{r}}{r}\dot{\phi},$
(74)
$\displaystyle\ddot{r}=-\frac{1}{2}\mu_{,r}\dot{r}^{2}+e^{-\mu}\left[\tilde{\eta}_{,r}\left(\dot{t}-\tilde{\eta}\dot{\phi}\right)\dot{\phi}+r\dot{\phi}^{2}\right].$
(75)

(73),(74), and (75) can be cast in a more suitable form, explicating the path
dependence on the impact parameter $b$, by integrating the first two equations
of motion by means of the constant of motions $E$ and $L$ and through the
defining relation $ds^{2}=0$ along the photon’s path. Thus, we obtain the
system of equations

$\displaystyle\dot{t}=e^{\mu/2}\frac{r|\dot{r}|}{\sqrt{r^{2}-\left(b-\tilde{\eta}\right)^{2}}}\left[1-\frac{\tilde{\eta}}{r^{2}}\left(b-\tilde{\eta}\right)\right],$
(76)
$\displaystyle\dot{\phi}=e^{\mu/2}\frac{r|\dot{r}|}{\sqrt{r^{2}-\left(b-\tilde{\eta}\right)^{2}}}\left(\frac{b-\tilde{\eta}}{r^{2}}\right),$
(77)
$\displaystyle\ddot{r}=\dot{r}^{2}\left[\frac{b-\tilde{\eta}}{r^{2}-\left(b-\tilde{\eta}\right)^{2}}\left(\tilde{\eta}_{,r}+\frac{b}{r}-\frac{\tilde{\eta}}{r}\right)-\frac{1}{2}\mu_{,r}\right].$
(78)

Now, (78) can be exactly integrated. Indeed, we notice that

$\displaystyle\ddot{r}$
$\displaystyle=\dot{r}^{2}\left[\frac{b-\tilde{\eta}}{r^{2}-\left(b-\tilde{\eta}\right)^{2}}\left(\tilde{\eta}_{,r}+\frac{b}{r}-\frac{\tilde{\eta}}{r}\right)-\frac{1}{2}\mu_{,r}\right]$
$\displaystyle=\dot{r}^{2}\left[\frac{1}{2}\log\left(1-\frac{\left(b-\tilde{\eta}\right)^{2}}{r^{2}}\right)_{,r}-\frac{1}{2}\mu_{,r}\right],$
(79)

which implies

$\dot{r}=Ce^{-\mu/2}\sqrt{1-\frac{\left(b-\tilde{\eta}\right)^{2}}{r^{2}}},$
(80)

where $C$ is a constant of integration related to the choice of the affine
parameter.222Notice that the choice of the affine parameter is related to the
normalisation of the photon’s energy. Indeed, $C=\pm 1$ corresponds to a
normalisation of $E$ to $E=1$ as could be directly inferred from the relation
$p_{\mu}p^{\mu}=0$ valid for photons. In particular, its value can always be
taken as $\pm 1$ (the sign defines whether the light ray is moving towards or
away from the lens). Thus, we can write

$\dot{r}=\pm
e^{-\mu/2}\sqrt{1-\frac{\left(b-\tilde{\eta}\right)^{2}}{r^{2}}},$ (81)

By substituting (81) into (76) and (77) we finally obtain

$\displaystyle\dot{t}=1+\frac{\tilde{\eta}}{r^{2}}\left(b-\tilde{\eta}\right),$
(82) $\displaystyle\dot{\phi}=\frac{b-\tilde{\eta}}{r^{2}}.$ (83)

(81), (82) and (83) constitute the first-order system of differential
equations which fully characterises photons’ motion on the equatorial plane
for the BG metric. Using (82) we can write (70) as

$\displaystyle\Delta t(|b|)=$
$\displaystyle\int_{\gamma_{+|b|}}\dot{t}d\lambda-\int_{\gamma_{-|b|}}\dot{t}d\lambda=\int_{\gamma_{+|b|}}\frac{\dot{t}}{\dot{r}}dr-\int_{\gamma_{-|b|}}\frac{\dot{t}}{\dot{r}}dr$
$\displaystyle=$
$\displaystyle\int^{r_{+}}_{D_{SL}}-\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}\right)/r)^{2}}}e^{\frac{\mu}{2}}dr+\int_{r_{+}}^{D_{SO}}\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}\right)/r)^{2}}}e^{\frac{\mu}{2}}dr-\\{r_{+}\longleftrightarrow
r_{-};+|b|\longleftrightarrow-|b|\\}$ $\displaystyle=$ $\displaystyle
2\int_{r_{+}}^{\text{min}(D_{SL},D_{LO})}\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}\right)/r)^{2}}}e^{\frac{\mu}{2}}dr+\int^{\text{max}(D_{SL},D_{LO})}_{\text{min}(D_{SL},D_{LO})}\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}\right)/r)^{2}}}e^{\frac{\mu}{2}}dr$
$\displaystyle-\\{r_{+}\longleftrightarrow
r_{-};+|b|\longleftrightarrow-|b|\\}$ $\displaystyle=$
$\displaystyle\bar{\Delta
t}(\pm|b|,r_{\pm},\min\\{D_{SL},D_{LO}\\})+\hat{\Delta
t}(\pm|b|,D_{SL},D_{LO}),$ (84)

where $D_{SL}$ is the source-lens distance, $D_{LO}$ is the lens-observer
distance, $r_{\pm}$ are the closest distance to the centre of the lens for a
given choice of $|b|$, respectively, for co-rotating and counter-rotating
photons and $\\{r_{+}\longleftrightarrow r_{-};+|b|\longleftrightarrow-|b|\\}$
indicates the exchanges of $r_{+}$ with $r_{-}$ and of $+|b|$ with $-|b|$ in
the integrals. Moreover, the first $\Delta t$ in (84), $\bar{\Delta t}$, is
defined as the difference between the integral whose lower extreme is $r_{+}$
and the one whose lower extreme is $r_{-}$, whereas the second $\Delta t$,
$\hat{\Delta t}$, is the difference between the remaining two integrals, the
first with $+|b|$, and the second with $-|b|$. (84) has the advantage of
writing $\Delta t(|b|)$ exclusively as the difference of integrals over the
radial coordinate of purely radial terms whose functional form is known.
Therefore, if the extremes of integration are known, the integral can be
computed numerically. To begin with, given $|b|$, $r_{\pm}$ can be obtained by
solving for $r$ the algebraic equations

$r^{2}-(\pm|b|-\tilde{\eta}(r))^{2}=0.$ (85)

(85) can be solved numerically by using the fzero built-in Matlab function
specialized for non-linear root-finding and by taking the positive root close
to $+|b|$. On the other hand, the values of $D_{SL}$ and $D_{LO}$ have to be
picked based on those obtained through observations of galaxy lenses. On
average, it is found that $D_{SL}<D_{LO}$ and, most importantly, that both are
of the order of some Gpc [95]. Therefore, we can infer that $\Delta t$ is of
the order of the first $\Delta t$ in (84), $\bar{\Delta t}$. Indeed, since
$r\geq 1$ Gpc when we calculate the integrals that define $\hat{\Delta t}$, we
are carrying out the integration over an almost flat region of the BG
spacetime. We then expect the contribution to the overall time difference of
the second $\Delta t$ in (84), $\hat{\Delta t}$, to be negligible with respect
to that of yhe first $\Delta t$, $\bar{\Delta t}$. Thus, we infer that

$\Delta t\approx\bar{\Delta t}.$ (86)

This is, in fact, confirmed by numerical calculations333See Appendix B for a
complete explanation of the method implemented. of the two $\Delta t$’s,
$\hat{\Delta t}$ and $\bar{\Delta t}$, which show that the ratio of the former
to the latter is of the order of $10^{-7}/10^{-8}$ for reasonable values of
$|b|$ and Milky Way-like galaxies. In particular, the two contributions to
$\Delta t$ have been estimated for $|b|$ values going from 50 to 1000 kpc
(step of 25 kpc) and the same galaxy parameters used for the calculation of
the deflection angle. $D_{SL}$ has been kept fixed at 1.5 Gpc, and $D_{SO}$
has been considered as 1.5, 5 and 10 Gpc. Over the whole physically-informed
parameter space explored, we have found the BG time delays to be orders of
magnitudes larger than the observed ones. Truly, the typical delay time for
galaxy lensing is of the order of months [106, 107, 108], in flat contrast
with the obtained values for BG of several decades. The full behaviour of
$\bar{\Delta t}$ as a function of $|b|$ is reported in Fig. 6.

Figure 6: $\bar{\Delta t}$ as a function of $|b|$. $R=47.87$ kpc, $r_{0}=0.39$
kpc and ${V}_{0}=263.1$ km/s and $D_{SL}=1.5$ Gpc. The negative values
obtained indicate, as expected, that the light ray undergoing retrograde
motion takes more time than the one moving with prograde motion.

Such a stark discrepancy between the observed values and those predicted by BG
for an observable which we expect to depend strongly on the frame-dragging, as
is the time delay between prograde and retrograde motion, can be understood
when we consider that BG rotates rigidly. Indeed, this feature is not only
patently unphysical, but it also generates an unreasonably large dragging, the
latter being evidently responsible for such huge and absurd time delays. That
the dragging in BG is so unreasonably strong is also independently highlighted
in [68], where it is shown in particular that in BG the effect of dragging
compensates entirely the special relativistic Doppler shift of photons emitted
by rotating matter. These effects of dragging in BG seem to us an indication
that it might indeed be the dragging due to the angular momentum and not so
much conventional DM responsible for, perhaps, a considerable part of the DM
effect.

##  6 Conclusions and perspectives

When, in 1919, Arthur Eddington and Frank W. Dyson measured for the first time
the gravitational deflection of starlight grazing the Sun [109], the results
left the world breathless: Albert Einstein had been right once more. Indeed,
not only were the paths of light rays affected by gravity, but the angle of
deflection was precisely what Einstein’s theory of General Relativity
predicted [110]. The gravitational lensing effect by our sun was one of the
first experimental verifications of GR. It was crucial in convincing the
physicist’s community of the validity of the theory [111]. Gravitational
lensing was used back in 1919 to decide between the two different theories of
gravity, two different models of the universe.

Here, we are proposing to use again gravitational lensing in this fashion,
albeit on a smaller philosophical scale, to choose once more between models:
models of galaxies. The observational data available today allows us to put
quite strict bounds on the predictions of new models of galaxies. If a model
produces estimates in direct contradiction with the data, it should be
discarded.

The calculations put forth in the present work show that the BG equatorial
lensing produces values for the deflection angle in agreement with those
measured for typical gravitational lensing by galaxies [94, 95, 96, 97].
However, the model produces a time delay which is orders of magnitudes away
from those currently observed [106, 107, 108]. These results can be
theoretically understood when one considers that, of the two chosen lensing
observables, only the second one is strongly sensitive to the totally
unphysical property of BG: the rigid rotation of the galaxy. Indeed, suppose
we were to slice the galaxy into concentric, small, circular rings. Then, the
condition of rigid rotation could be loosely translated to mimic the presence
of extreme friction between these rings, which forces them to move at the same
rotational velocity. In such a situation, the light ray would be subject to an
enormous drag throughout its path. Therefore, the continuous dragging action
would add up, over the paths of the photons, to huge differences in the travel
time for co-rotating and counter-rotating light rays. On the other hand, the
same dragging comes into play as a second-order effect when we look at light-
bending angles, the latter being instead mostly determined by the mass of the
lens. Thus, since BG predicts masses for Milky Way-like galaxies of the right
order of magnitudes as the actual ones [100, 101, 102, 103], it is not
surprising that the predicted deflection angles do indeed match the observed
ones.

As such, the conflict between the theoretical prediction for BG and the
measured values for the time delay disproves BG as a viable model for disc
galaxies, even though it seems to reproduce the rotation curves of the Milky
Way correctly [65]. Nonetheless, there are a couple of points which are still
worthy of consideration. The galaxy lensings actually observed are related to
face-on or tilted galaxies [96]. Hence, to be entirely sure of our
conclusions, we should investigate the BG lensing in those conditions, most
likely through a numerical approach. Moreover, the analysis put forth in this
work was carried out by implementing the same coordinates chosen to describe
the galaxy in the original work of Balasin and Grumiller [64]. In particular,
a $2\pi$ periodicity was assumed for the $\phi$ coordinate, which thus implies
a local asymptotic flatness at spatial infinity, but not a global one. As some
possible problems with this choice of coordinates, away from the equatorial
plane, were already pointed out in [64], we believe that further investigation
is required to confute any concern regarding their use. We are currently
investigating this issue.

We conclude by mentioning some future perspectives on the studies of disc
galaxies. As stated, the limiting condition of rigid rotation for the BG
galaxy has been proven to create conflict with the observational data when we
look at gravitational lensing. The $(\eta,H)$ class of models investigated in
[67, 68], which generalises BG by dropping rigidity and allowing for
differential rotation, may move in the right direction concerning the
agreement between theoretical predictions and the results of observations. In
addition, it has also been argued in [67, 68] that the $(\eta,H)$ class of
models may allow for a further reduction of the need for DM in a disc galaxy.
Therefore, as reasons to move on from the rigid rotation approximation begin
to pile-up, the need to explore further the whole $(\eta,H)$ class grows. In
particular, a solution for the $(\eta,H)$ class with differential rotation
might help in understanding the dynamics of the so called ultra-diffuse
galaxies (UDG), namely of fluffy galaxies whose star density is very low and
which are spread over vast distances and which appear to contain very little,
or no DM at all [57, 58, 59]. However, as indicated in [67], something more
than the simple $(\eta,H)$ framework is required to model young disc galaxies,
as the former can only be applied to the description of old dust galaxies.
Hence, in future works, we plan to move beyond the dust approximation. Indeed,
we aim to introduce pressure terms that might help in better modelling for
active, young galaxies and their bulge. Nonetheless, as it stands, we have
shown the potential of the lensing observables for the selection of full GR
models for galaxies, and we plan to make use of them in future works. In
particular, if the current techniques employed in the calculations were to
result too limited for the new scenarios, we aim to modify them appropriately.
Finally, as already mentioned, we plan to numerically investigate non-
equatorial gravitational lensing in BG as well as in the full $(\eta,H)$
model.

##  Appendix A Planar motion and radial geodesic in the equatorial plane for
BG.

To prove the existence of planar motion on the equatorial plane for photons in
the BG model we consider (8) for $\ddot{z}$ and we impose the initial
condition $\dot{z}=0$

$e^{\mu}\ddot{z}=\frac{1}{2}\gamma_{rr,z}\dot{r}^{2}+\beta_{\phi,z}\dot{\phi}=\frac{1}{2}e^{\mu}\mu_{,z}\dot{r}^{2}+\eta_{,z}\dot{\phi}.$
(87)

The condition $\ddot{z}=0$ is then equivalent to

$\frac{1}{2}e^{\mu}\mu_{,z}\dot{r}^{2}+\eta_{,z}\dot{\phi}=0.$ (88)

However, given the condition (43) of the BG model, (88) is automatically
satisfied on the equatorial plane.

As for what concerns the existence of radial geodesics on the equatorial plane
for the three-geometry described by (2), the result is immediately verified if
we consider (8) for $\ddot{z}$ and $\ddot{\phi}$ with the initial conditions
$\dot{z}=0$ and $\dot{\phi}=0$ whilst keeping in mind (43).

##  Appendix B Numerical strategy for the estimation of the time delay in the
equatorial plane for BG.

We intend to gauge numerically the values of $\bar{\Delta
t}(\pm|b|,r_{\pm},D_{SL})$ for different choices of $|b|$ and of $D_{SL}$
(where, w.l.o.g., we have implicitly taken $D_{SL}<D_{LO}$). As by the results
of section 5, $\bar{\Delta t}$ is given by

$\bar{\Delta
t}=2\underbrace{\int_{r_{+}}^{D_{SL}}\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}\right)/r)^{2}}}e^{\frac{\mu}{2}}dr}_{I_{+}}-2\underbrace{\int_{r_{-}}^{D_{SL}}\frac{\left[1-\tilde{\eta}\left(|b|+\tilde{\eta}\right)/r^{2}\right]}{\sqrt{1-(\left(|b|+\tilde{\eta}\right)/r)^{2}}}e^{\frac{\mu}{2}}dr}_{I_{-}}.$
(89)

Both $I_{\pm}$ are improper integrals as their integrands diverge at the lower
ends. As a consequence, estimating numerically $\bar{\Delta t}$ asks for some
extra care. Various approaches can be adopted to this aim; see, e.g., [112].
Above all, one could proceed by operating a change of variable that either
eliminates the singularity of the integrands or introduces a weight function.
The first option allows us to use routine quadrature rules, e.g., trapezoidal
or Simpson’s, while the second one opens to the implementation of specialized
high-precision formulas, e.g., Gauss-Jacobi quadratures (see, e.g., [113] for
more details on quadrature rules). Disposing of these strategies comes at the
price of expressing $r$ in (89) in terms of the newly introduced variable,
which, due to the presence of $\tilde{\eta}$, could result far from trivial.

A technique known as ‘subtracting out the singularity’ is a valuable
alternative to change of variable and it consists in approximating the
singular integrand with a function that shares the same kind of singularity,
being at the same time easier to be integrated than the original one. A
modification of this approach naturally applies to the numerical calculation
of (89) by observing that $\bar{\Delta t}$ is a difference of two integrals
whose integrands share the same kind of singularity at $r_{\pm}$. Indeed, if
we reduce both $I_{\pm}$ to the same integration interval by the following
translation

$r\longrightarrow y=r-r_{\pm},$ (90)

where $r_{+}$ is chosen for the first integral and $r_{-}$ for the second one,
we can cast (89) into the form

$\displaystyle\bar{\Delta t}=$ $\displaystyle
2\int_{0}^{D_{SL}-r_{-}}\left[\frac{\left[1+\tilde{\eta}_{+}\left(|b|-\tilde{\eta}_{+}\right)/(y+r_{+})^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}_{+}\right)/(y+r_{+}))^{2}}}e^{\frac{\mu_{+}}{2}}-\frac{\left[1-\tilde{\eta}_{-}\left(|b|+\tilde{\eta}_{-}\right)/(y+r_{-})^{2}\right]}{\sqrt{1-(\left(|b|+\tilde{\eta}_{-}\right)/(y+r_{-}))^{2}}}e^{\frac{\mu_{-}}{2}}\right]dy$
$\displaystyle+2\int_{D_{SL}-r_{-}}^{D_{SL}-r_{+}}\frac{\left[1+\tilde{\eta}_{+}\left(|b|-\tilde{\eta}_{+}\right)/(y+r_{+})^{2}\right]}{\sqrt{1-(\left(|b|-\tilde{\eta}_{+}\right)/(y+r_{+}))^{2}}}e^{\frac{\mu_{+}}{2}}dy,$
(91)

where $\tilde{\eta}_{\pm}$ and $\mu_{\pm}$ indicate the functions
$\tilde{{\eta}}$ and $\mu$ in which the change of variables of (90) have been
carried out, and for the determination of the integration bounds, we have used
that $r_{+}<r_{-}$. (91) allows for a reliable numerical estimation of
$\bar{\Delta t}$ as in the first integral, the integrand does not diverge
anymore when we value it on the lower bound of integration, and clearly, the
second integral does not pose any challenge to numerical integration.
Therefore, we can reliably apply standard quadrature rules for the numerical
integration of (91) and the estimation of $\bar{\Delta t}$.

A further alternative to both change of variable and subtraction of the
singularity consists in relying on adaptive quadratures. This is done by
recursively splitting the integration interval into pieces and applying a
basic formula to each piece. The interval is split in a way adapted to the
behavior of the integrand (long subintervals where the integrand is easy to
approximate, short ones where it is difficult). Proceeding in this manner, a
prescribed accuracy can be reached while facing automatically the singularity
of the integrand. Based on these considerations and on the fact that Matlab is
equipped with two built-in functions [114], integral (double precision, non
symbolic) and vpaintegral (variable precision, symbolic), that implement
numerical integration by means of global adaptive quadratures, we compute
‘almost exactly’ $\alpha$ in (69) by vpaintegral with a relative tolerance of
$10^{-12}$ and provide satisfactory estimates of $\bar{\Delta t}$ by integral
with default relative tolerance of $10^{-6}$.

For the sake of completeness, we conclude by mentioning that another way of
proceeding is to integrate by parts a recasting of (89) obtained by
multiplying and dividing the integrands of $I_{\pm}$ by the derivative of
$r^{2}-(|b|\pm\tilde{\eta})^{2}$. This brings to writing $I_{\pm}$ as follows

$\displaystyle I_{+}$
$\displaystyle=\int_{r_{+}}^{D_{SL}}\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{(r+\frac{{\rm
d}\tilde{\eta}}{{\rm d}r}(|b|-\tilde{\eta}))}re^{\frac{\mu}{2}}\frac{{\rm
d}}{{\rm d}r}\sqrt{r^{2}-(|b|-\tilde{\eta})^{2}}dr$
$\displaystyle=\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{(r+\frac{{\rm
d}\tilde{\eta}}{{\rm
d}r}(|b|-\tilde{\eta}))}re^{\frac{\mu}{2}}\sqrt{r^{2}-(|b|-\tilde{\eta})^{2}}\bigg{|}_{r=D_{SL}}-\int_{r_{+}}^{D_{SL}}\sqrt{r^{2}-(|b|-\tilde{\eta})^{2}}\frac{{\rm
d}}{{\rm
d}r}\left(\frac{\left[1+\tilde{\eta}\left(|b|-\tilde{\eta}\right)/r^{2}\right]}{(r+\frac{{\rm
d}\tilde{\eta}}{{\rm d}r}(|b|-\tilde{\eta}))}re^{\frac{\mu}{2}}\right)dr$ (92)
$\displaystyle I_{-}$
$\displaystyle=\int_{r_{-}}^{D_{SL}}\frac{\left[1-\tilde{\eta}\left(|b|+\tilde{\eta}\right)/r^{2}\right]}{(r-\frac{{\rm
d}\tilde{\eta}}{{\rm d}r}(|b|+\tilde{\eta}))}re^{\frac{\mu}{2}}\frac{{\rm
d}}{{\rm d}r}\sqrt{r^{2}-(|b|+\tilde{\eta})^{2}}dr$
$\displaystyle=\frac{\left[1-\tilde{\eta}\left(|b|+\tilde{\eta}\right)/r^{2}\right]}{(r-\frac{{\rm
d}\tilde{\eta}}{{\rm
d}r}(|b|+\tilde{\eta}))}re^{\frac{\mu}{2}}\sqrt{r^{2}-(|b|+\tilde{\eta})^{2}}\bigg{|}_{r=D_{SL}}-\int_{r_{-}}^{D_{SL}}\sqrt{r^{2}-(|b|+\tilde{\eta})^{2}}\frac{{\rm
d}}{{\rm
d}r}\left(\frac{\left[1-\tilde{\eta}\left(|b|+\tilde{\eta}\right)/r^{2}\right]}{(r-\frac{{\rm
d}\tilde{\eta}}{{\rm d}r}(|b|+\tilde{\eta}))}re^{\frac{\mu}{2}}\right)dr$ (93)

where the resulting integrals are not affected by any numerical criticality
and can then be treated by basic quadrature rules. Note that, by leveraging on
Matlab symbolic tools, we are not asked to explicitly compute any of the
derivatives appearing in (92)-(93).

## References

  * [1] F. Zwicky “Die Rotverschiebung von extragalaktischen Nebeln” In _Helvetica Physica Acta_ 6, 1933, pp. 110
  * [2] F. Zwicky “On the Masses of Nebulae and of Clusters of Nebulae” In _Astrophysical Journal_ 86, 1937, pp. 217
  * [3] J.H. Oort “Some Problems Concerning the Structure and Dynamics of the Galactic System and the Elliptical Nebulae NGC 3115 and 4494” In _Astrophysical Journal_ 91, 1940, pp. 273
  * [4] Rubin V. C. and W.. Jr “Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions” In _Astrophysical Journal_ 159, 1970, pp. 379
  * [5] Rubin V. C. and W.. Jr “Extended rotation curves of high-luminosity spiral galaxies. IV – Systematic dynamical properties, SA through SC” In _The Astrophysical Journal Letters_ 225, 1978, pp. L107–L111
  * [6] Y. Sofue and V.. Rubin “Rotation curves of spiral galaxies” In _Ann.Rev.Astron.Astrophys._ 39, 2001, pp. 137
  * [7] B. Famaey and S.. McGaugh “Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions” In _Living Reviews in Relativity_ 15, 2012, pp. 10
  * [8] S.. McGaugh “The Third Law of Galactic Rotation” In _Galaxies_ , 2014
  * [9] S. Trippe “The ‘Missing Mass Problem’ in Astronomy and the Need for a Modified Law of Gravity” In _Zeitschrift für Naturforschung A_ , 2014
  * [10] S.. McGaugh “A tale of two paradigms: the mutual incommensurability of $\Lambda$CDM and MOND” In _Canadian Journal of Astrophysics_ , 2015
  * [11] X. Wu and P. Kroupa “Galactic rotation curves, the baryon-to-dark-halo-mass relation and space–time scale invariance” In _Monthly Notices of the Royal Astronomical Society_ 446, 2015
  * [12] P. Kroupa “Galaxies as simple dynamical systems: observational data disfavor dark matter and stochastic star formation” In _Canadian Journal of Physics_ 93, 2015
  * [13] H. Haghi, A.. Bazkiaei, A. Zonoozi and P. Kroupa “Declining rotation curves of galaxies as a test of gravitational theory” In _Monthly Notices of the Royal Astronomical Society_ 458, 2016
  * [14] M. Haslbauer et al. “The High Fraction of Thin Disk Galaxies Continues to Challenge $\Lambda$CDM Cosmology” In _The Astrophysical Journal_ 925, 2022
  * [15] N. Wittenburg et al. “Hydrodynamical structure formation in Milgromian cosmology” In _Canadian Journal of Astrophysics_ 523, 2023
  * [16] S.D. S. and V.K. Oikonomou “Modified gravity theories on a nutshell: Inflation, bounce and late-time evolution” In _Physics Report_ 692, 2017, pp. 1
  * [17] Bing-Lin Young “A survey of dark matter and related topics in cosmology” In _Frontiers of Physics_ 12, 2016
  * [18] Planck Collaboration “Planck 2018 results. I. Overview and the cosmological legacy of Planck” In _A &AP_ 641, 2020
  * [19] Edvige Corbelli and Paolo Salucci “The extended rotation curve and the dark matter halo of M33” In _Monthly Notices of the Royal Astronomical Society_ 311, 2000, pp. 441–447
  * [20] Allen “Cosmological Parameters from Observations of Galaxy Clusters” In _ARAA_ 409, 2011
  * [21] A. Taylor “Gravitational lens magnification and galaxy cluster masses” In _Evolution of Large Scale Structure: From Recombination to Garching_ , 1999, pp. 226
  * [22] Douglas Clowe et al. “A Direct Empirical Proof of the Existence of Dark Matter*” In _The Astrophysical Journal_ 648.2, 2006, pp. L109
  * [23] Maruša Bradač et al. “Revealing the Properties of Dark Matter in the Merging Cluster MACS J0025.4–1222*” In _The Astrophysical Journal_ 687.2, 2008, pp. 959
  * [24] M.J. Jee, A. Mahdavi and H. Hoekstra “A Study of the Dark Core in A520 with Hubble Space Telescope: The Mystery Deepens” In _The Astrophysical Journal_ 747, 2012, pp. 96
  * [25] C. Conroy, Risa H. Wechsler and Andrey V. Kravtsov “Modeling Luminosity-dependent Galaxy Clustering through Cosmic Time” In _The Astrophysical Journal_ 647.1, 2006, pp. 201–214
  * [26] Oliver Piattella “Lecture Notes in Cosmology” Springer International Publisher, 2018
  * [27] Gianfranco Bertone and Dan Hooper “History of dark matter” In _Reviews of Modern Physics_ 90, 2018
  * [28] Jianglai Liu, Xun Chen and Xiangdong Ji “Current status of direct dark matter detection experiments” In _Nature Physics_ 13, 2017
  * [29] Teresa Marrodán Undagoitia and Ludwig Rauch “Dark matter direct-detection experiments” In _Journal of Physics G: Nuclear and Particle Physics_ 43, 2015
  * [30] E. al. Aprile “Excess electronic recoil events in XENON1T” In _Phys. Rev. D_ 102, 2020
  * [31] J. Aalbers “First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment”, e-print arXiv:hep-ex: 2207.03764, 2022
  * [32] Simon Viel “Dark matter search results from DEAP-3600 at SNOLAB” In _PoS_ ICHEP2020, 2021, pp. 655
  * [33] Yue et al. Meng “Dark Matter Search Results from the PandaX-4T Commissioning Run” In _Phys. Rev. Lett._ 127, 2021
  * [34] C Krishnan et al. “Does Hubble tension signal a breakdown in FLRW cosmology?” In _Classical and Quantum Gravity_ 38, 2021
  * [35] Elcio Abdalla “Cosmology intertwined: A review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies” In _Journal of High Energy Astrophysics_ 34, 2022, pp. 49–211
  * [36] Kroupa, P. et al. “Local-Group tests of dark-matter concordance cosmology - Towards a new paradigm for structure formation” In _A &A_ 523, 2010
  * [37] P.. Aluri “Is the Observable Universe Consistent with the Cosmological Principle?”, e-print arXiv:astro-ph: 2207.05765, 2023
  * [38] Avishai Dekel “Efficient formation of massive galaxies at cosmic dawn by feedback-free starbursts” In _Monthly Notices of the Royal Astronomical Society_ 523, 2023
  * [39] B.. Keller “Can Cosmological Simulations Reproduce the Spectroscopically Confirmed Galaxies Seen at z $\leq$ 10?” In _The Astrophysical Journal Letters_ 943, 2023
  * [40] F. Prada “Confirmation of the standard cosmological model from red massive galaxies $\sim$600 Myr after the Big Bang”, e-print arXiv:astro-ph: 2304.11911, 2023
  * [41] J. McCaffrey “No Tension: JWST Galaxies at z$>$10 Consistent with Cosmological Simulations”, e-print arXiv:astro-ph: 2304.13755, 2023
  * [42] David Harvey, F. Courbin, J.. Kneib and Ian G. McCarthy “A detection of wobbling brightest cluster galaxies within massive galaxy clusters” In _Monthly Notices of the Royal Astronomical Society_ 472.2, 2017, pp. 1972–1980
  * [43] Elena Asencio et al. “The distribution and morphologies of Fornax Cluster dwarf galaxies suggest they lack dark matter” In _Monthly Notices of the Royal Astronomical Society_ 515, 2022, pp. 2981–3013
  * [44] Oliver Müller, Marcel S. Pawlowski, Helmut Jerjen and Federico Lelli “A whirling plane of satellite galaxies around Centaurus A challenges cold dark matter cosmology” In _Science_ 359, 2018, pp. 534–537
  * [45] Michal Bílek, Ingo Thies, Pavel Kroupa and Benoit Famaey “Are disks of satellites comprised of tidal dwarf galaxies?” In _Galaxies_ 9, 2021, pp. 100
  * [46] Müller, Oliver “The dwarf galaxy satellite system of Centaurus A” In _A &A_, 2019, pp. A18
  * [47] Sanjaya Paudel, Suk-Jin Yoon and Rory Smith “A Corotating Group of Dwarf Galaxies around NGC 2750 as a Centaurus A Analog” In _The Astrophysical Journal Letters_ 917.2, 2021, pp. L18
  * [48] D. Lynden-Bell “Dwarf Galaxies and Globular Clusters in High Velocity Hydrogen Streams” In _Monthly Notices of the Royal Astronomical Society_ 174, 1976, pp. 695–710
  * [49] Heesters, Nick et al. “Flattened structures of dwarf satellites around massive host galaxies in the MATLAS low-to-moderate density fields” In _A &A_ 654, 2021, pp. A161
  * [50] Marco Castellano “Early Results from GLASS-JWST. III. Galaxy Candidates at z $\sim$9–15$*$” In _The Astrophysical Journal Letters_ 938, 2022, pp. L15
  * [51] Rohan P. Naidu “Two Remarkably Luminous Galaxy Candidates at z $\approx$ 10–12 Revealed by JWST” In _The Astrophysical Journal Letters_ 940, 2022, pp. L14
  * [52] Ivo Labbé “A population of red candidate massive galaxies $\sim$600 Myr after the Big Bang” In _Nature_ , 2023
  * [53] Christopher C Lovell et al. “Extreme value statistics of the halo and stellar mass distributions at high redshift: are JWST results in tension with $\Lambda$CDM?” In _Monthly Notices of the Royal Astronomical Society_ 518.2, 2022, pp. 2511–2520
  * [54] Michael Boylan-Kolchin “Stress testing $\Lambda$CDM with high-redshift galaxy candidates” In _Nature Astronomy_ , 2023
  * [55] Clara Giménez-Arteaga “Spatially Resolved Properties of Galaxies at 5 ¡ z ¡ 9 in the SMACS 0723 JWST ERO Field” In _The Astrophysical Journal_ 948, 2023
  * [56] A. Amruth, T. Broadhurst and J. Lim “Einstein rings modulated by wavelike dark matter from anomalies in gravitationally lensed images” In _Nature Astronomy_ 7, 2023
  * [57] P. Dokkum “A galaxy lacking dark matter” In _Nature_ 555, 2018, pp. 629
  * [58] S. Trujillo-Gomez, J… Kruijssen and M. Reina-Campos “The emergence of dark matter-deficient ultra-diffuse galaxies driven by scatter in the stellar mass–halo mass relation and feedback from globular clusters” In _Mon. Not. R. Astron. Soc._ 510, 2021, pp. 3356
  * [59] Demao Kong et al. “The Odd Dark Matter Halos of Isolated Gas-rich Ultradiffuse Galaxies” In _The Astrophysical Journal_ 936, 2022
  * [60] D.. Rodrigues, V. Marra, A. Popolo and Z. Davari “Absence of a fundamental acceleration scale in galaxies” In _Nature Astronomy_ 2, 2018, pp. 668
  * [61] F.I. Cooperstock and S. Tieu “General relativity resolves galactic rotation without exotic dark matter”, e-print arXiv:astro-ph/0507619, 2005
  * [62] F.I. Cooperstock and S. Tieu “Galactic dynamics via general relativity: a compilation and new developments” In _Int. J. Mod. Phys._ 22, 2007, pp. 2293
  * [63] J. Carrick and Fred I. Cooperstock “General relativistic dynamics applied to the rotation curves of galaxies” In _Astrophys. Space Sci_ 337, 2012, pp. 321
  * [64] H. Balasin and D. Grumiller “Non-Newtonian behavior in weak field general relativity for extended rotating sources” In _Int. J. Mod. Phys. D_ 17, 2008, pp. 475
  * [65] M. Crosta, M. Giammaria, M.G. Lattanzi and E. Poggio “CDM and geometry-driven Milky Way rotation curve models with Gaia DR2” In _Mon. Not. R. Astron. Soc._ 496, 2020, pp. 2107
  * [66] William Beordo et al. “Geometry-driven and dark-matter-sustained Milky Way rotation curves with Gaia DR3” Manuscript submitted for publication, 2023
  * [67] D. Astesiano, S.. Cacciatori, V. Gorini and F. Re “Towards a full general relativistic approach to galaxies” In _Eur. Phys. J. C_ 82, 2022, pp. 554
  * [68] F. Re et al. “Re-weighting dark matter in disc galaxies: a new general relativistic observational test”, e-print arXiv:astro-ph:2204.05143, 2022
  * [69] D. Astesiano and M.. Ruggiero “Galactic dark matter effects from purely geometrical aspects of general relativity” In _Phys. Rev. D_ 106, 2022
  * [70] L. Ciotti “On the Rotation Curve of Disk Galaxies in General Relativity” In _The Astrophysical Journal_ 936, 2022
  * [71] H. Stephani et al. “Exact Solutions of Einstein’s Field Equations” Cambridge University Press, 2003
  * [72] J.. Islam “Rotating Fields in General Relativity” Cambridge University Press, 2009
  * [73] C.. Keeton “Cold Dark Matter and Strong Gravitational Lensing: Concord or Conflict?” In _The Astrophysical Journal_ 561, 2001, pp. 46
  * [74] R.. Metcalf and P. Madau “Compound Gravitational Lensing as a Probe of Dark Matter Substructure within Galaxy Halos” In _The Astrophysical Journal_ 563, 2001, pp. 9
  * [75] L.. Moustakas and R.. Metcalf “Detecting dark matter substructure spectroscopically in strong gravitational lenses” In _Mon. Not. R. Astron. Soc._ 339, 2003, pp. 607
  * [76] H. Hoekstra, H… Yee and M.. Gladders “Properties of Galaxy Dark Matter Halos from Weak Lensing” In _The Astrophysical Journal_ 606, 2004, pp. 67
  * [77] R. Mandelbaum et al. “Ellipticity of dark matter haloes with galaxy–galaxy weak lensing” In _Mon. Not. R. Astron. Soc._ 370, 2006, pp. 1008
  * [78] S.. Suyu “Disentangling baryons and dark matter in the spiral gravitational lens B1933+ 503” In _The Astrophysical Journal_ 750, 2012, pp. 10
  * [79] M. Velander “CFHTLenS: the relation between galaxy dark matter haloes and baryons from weak gravitational lensing” In _Mon. Not. R. Astron. Soc._ 437, 2013, pp. 2111
  * [80] C. Alcock “Experimental Limits on the Dark Matter Halo of the Galaxy from Gravitational Microlensing” In _Phys. Rev. Lett._ 74, 1995, pp. 2867
  * [81] A.. Tuntsov, G.. Lewis, R.. Ibata and J.. Kneib “Detecting compact dark matter in galaxy clusters via gravitational microlensing: A2218 and A370” In _Mon. Not. R. Astron. Soc._ 353, 2004, pp. 853
  * [82] M. Ricotti and A. Gould “A new probe of dark matter and high-energy universe using microlensing” In _The Astrophysical Journal_ 707, 2009, pp. 979
  * [83] F. Iocco, M. Pato, G. Bertone and P. Jetzer “Dark Matter distribution in the Milky Way: microlensing and dynamical constraints” In _J. of Cosm. and Astropar. Phys._ 2011, 2011, pp. 29
  * [84] T. Ono, A. Ishihara and H. Asada “Gravitomagnetic bending angle of light with finite-distance corrections in stationary axisymmetric spacetimes” In _Phys. Rev. D_ 96, 2017
  * [85] F. Giannoni, A. Masiello and P. Piccione “The Fermat principle in general relativity and applications” In _J. Math. Phys._ 43, 2002, pp. 44
  * [86] V.. Frolov “Generalized Fermat’s principle and action for light rays in a curved spacetime” In _Phys. Rev. D_ 88, 2013
  * [87] G.. Gibbons and M.. Werner “Applications of the Gauss-Bonnet theorem to gravitational lensing” In _Class. Quantum Grav._ 25, 2008
  * [88] M.. Wernre “Gravitational lensing in the Kerr-Randers optical geometry” In _Gen. Rel. Grav._ 44, 2012, pp. 3047
  * [89] A. Ishihara et al. “Gravitational bending angle of light for finite distance and the Gauss-Bonnet theorem” In _Phys. Rev. D_ 94, 2016
  * [90] J.. Bardeen “A variational principle for rotating stars in General Relativity” In _The Astrophysical Journal_ 162, 1970, pp. 71
  * [91] J.. Bardeen, W.. Press and S.. Teukolsky “Rotating Black Holes: Locally Nonrotating Frames, Energy Extraction, and Scalar Synchrotron Radiation” In _The Astrophysical Journal_ 178, 1972, pp. 347
  * [92] M.. Carmo “Differential geometry of curves and surfaces” Dover Books on Mathematics, 1976
  * [93] D.d W. Henderson “Differential Geometry: A Geometric Introduction” Prentice Hall, 1998
  * [94] L.C. Parker “The Masses and Shapes of Dark Matter Halos from Galaxy-Galaxy Lensing in the CFHT Legacy Survey” In _The Astrophysical Journal_ 669, 2007, pp. 21
  * [95] T.Treu “The Sloan Lens ACS Survey. IX. Colors, Lensing, and Stellar Masses of Early-Type Galaxies” In _The Astrophysical Journal_ 705, 2012, pp. 1099
  * [96] A. Sonnenfeld “The SL2S galaxy-scale lens sample. III. Lens models, surface photometry, and stellar masses for the final sample” In _The Astrophysical Journal_ 777, 2013, pp. 97
  * [97] T. Verdugo “Characterizing SL2S galaxy groups using the Einstein radius” In _Astronomy and Astrophysics_ 571, 2014, pp. 21
  * [98] S. Mollerach and E. Roule “Gravitational Lensing and Microlensing” World Scientific, 2002
  * [99] F. Garzón “Astrophysical Applications of Gravitational Lensing” Cambridge University Press, 2016
  * [100] J.. Ostriker, P… Peebles and A. Yahil “The size and mass of galaxies, and the mass of the universe” In _Astrophysical Journal_ 193, 1974
  * [101] D. Zaritsky, R. Smith, C. Frenk and S… White “Satellites of spiral galaxies” In _Astrophysical Journal_ 405, 1993
  * [102] D. Zaritsky and S… White “The massive halos of spiral galaxies” In _Astrophysical Journal_ 435, 1994
  * [103] Sofue Yoshiaki “Rotation and mass in the Milky Way and spiral galaxies” In _Astronomical Society of Japan_ 69, 2016
  * [104] L.. Landau and E.. Lifschits “The Classical Theory of Fields” Volume 2, Course of Theoretical Physics Oxford: Pergamon Press, 1975
  * [105] L. O. Costa, José Natário, F. Frutos-Alfaro and Michael Soffel “Reference frames in General Relativity and the galactic rotation curves”, e-print arXiv:gr-qc 2303.17516, 2023
  * [106] A.. Biggs “Time delay for the gravitational lens system B0218+357” In _Mon. Not. R. Astron. Soc._ 304, 1999, pp. 349
  * [107] J.. Cohn, C.. Kochanek, B.. McLeod and C.. Keeton2 “Constraints on Galaxy Density Profiles from Strong Gravitational Lensing: The Case of B1933+503” In _The Astrophysical Journal_ 554, 2001, pp. 1216
  * [108] M. Oguri, A. Taruya, Y. Suto and E.. Turner “Strong Gravitational Lensing Time Delay Statistics and the Density Profile of Dark Halo” In _The Astrophysical Journal_ 568, 2002, pp. 48
  * [109] E. John and G. Clark “Relativity and eclipses: the British eclipse expeditions of 1919 and their predecessors” In _Hist. Stud. Phys. Sci._ 11, 1980
  * [110] Albert Einstein “Collected Papers of Albert Einstein Vol. 3” Princeton University Press, 1994
  * [111] D. Valls-Gabaud “A Conceptual and Cultural History of Gravitational Lensing” In _Cosmology Across Cultures_ 409, 2008
  * [112] Lawrence F. Shampine, Richard C. Allen and Steven Pruess “Fundamentals of Numerical Computing” Wiley, 1997
  * [113] David Ronald Kincaid and Elliott Ward Cheney “Numerical analysis: mathematics of scientific computing” Brooks/Cole Publishing Company, 1991
  * [114] L.F. Shampine “Vectorized adaptive quadrature in MATLAB” In _Journal of Computational and Applied Mathematics_ 211.2, 2008, pp. 131–140