# Exploring chemical compound space with a graph-based recommender system

Elton Ogoshi<EMAIL_ADDRESS>Henrique Ferreira  Center for Natural
and Human Sciences, Federal University of ABC, Santo André, SP, Brazil  João
N. B. Rodrigues  Center for Natural and Human Sciences, Federal University of
ABC, Santo André, SP, Brazil  Gustavo M. Dalpian<EMAIL_ADDRESS>Center for
Natural and Human Sciences, Federal University of ABC, Santo André, SP, Brazil

###### Abstract

With the availability of extensive databases of inorganic materials, data-
driven approaches leveraging machine learning have gained prominence in
materials science research. In this study, we propose an innovative adaptation
of data-driven concepts to the mapping and exploration of chemical compound
space. Recommender systems, widely utilized for suggesting items to users,
employ techniques such as collaborative filtering, which rely on bipartite
graphs composed of users, items, and their interactions. Building upon the
Open Quantum Materials Database (OQMD), we constructed a bipartite graph where
elements from the periodic table and sites within crystal structures are
treated as separate entities. The relationships between them, defined by the
presence of ions at specific sites and weighted according to the thermodynamic
stability of the respective compounds, allowed us to generate an embedding
space that contains vector representations for each ion and each site. Through
the correlation of ion-site occupancy with their respective distances within
the embedding space, we explored new ion-site occupancies, facilitating the
discovery of novel stable compounds. Moreover, the graph’s embedding space
enabled a comprehensive examination of chemical similarities among elements,
and a detailed analysis of local geometries of sites. To demonstrate the
effectiveness and robustness of our method, we conducted a historical
evaluation using different versions of the OQMD and recommended new compounds
with Kagome lattices, showcasing the applicability of our approach to
practical materials design.

###### keywords:

new materials discovery, recommender system

Paulista Faculty of Informatics and Administration, São Paulo, SP, Brazil
Physics Departament, Federal University of Pernambuco, Recife, PE, Brazil
Institute of Physics, University of São Paulo, São Paulo, SP, Brazil

##  1 Introduction

The set of all admissible compounds results from all possible combinations of
atoms, compositions1 and crystal structures2 (ACS3). The combinatorially large
nature of this space makes it extremely difficult to thoroughly explore it. A
naive approach to expand the quantity of available known materials would be a
trial-and-error methodology across all possibilities, but this is is rendered
impractical due to resource limitations, a constraint that holds true even for
in silico exploration. Despite the significant recent advancements in high-
throughput infrastructure and computational power, the most effective strategy
for researchers remains the identification of patterns within known compounds
to effectively restrict this expansive search space.

In the past, there were some attempts to formalize such patterns in the form
of empirical geometrical rules, such as the Goldschmidt tolerance factor 4 for
perovskites, and the Hume-Rothery rules for the formation of substitutional
solid solutions in metallic systems 5. However, when these empirical rules
were formulated, the available dataset of materials was considerably limited
compared to the vast repositories we have today. The advent of comprehensive
materials databases, such as the Open Quantum Materials Database (OQMD) 6,
Aflowlib 7, Materials Project 8, and NOMAD9 has significantly expanded the
compounds’ knowledge base. This broadens the prospects for both finding new
empirical rules and improving the existing ones. For instance, a larger
dataset of compounds allowed Bartel et al 10 to propose a new and more
accurate tolerance factor for perovskites, while Pei et al 11 found new
descriptors that predicted the stability of high-entropy alloys beyond the
Hume-Rothery rules. Yet, these approaches are typically limited to a specific
subset of compounds constrained by a crystal structure.

In recent years, inspired by data-driven methodologies, researchers have
proposed innovative strategies to narrow the search space across a broader set
of compounds. These strategies can be broadly categorized into iterative and
non-iterative methods. Iterative methods employ a cyclical process of
exploration and exploitation, with examples including active learning 12,
genetic algorithms 13, and iterative ion substitution 14. Non-iterative
methods, on the other hand, build models from the entirety of the available
data. Notable such examples include the tensor-based recommender system
proposed by Seko et al 15, 16, and state-of-art graph-based machine learning
models that predict the formation energy of unexplored compounds, thereby
helping to restrict the search space 17, 18, 19.

Another context where one has to cope with an enormous informational space is
the realm of digital life. There users are routinely presented with a vast
array of options, be it media to consume or products to purchase. They have to
parse through all that information and make decisions that are most aligned
with their interests. Some of the most popular information filtering
strategies are known as recommender systems 20 (RS). These attempt to simplify
the decision-making process by first trying to identify items that are
perceived to most closely align with the user’s interests and then suggesting
these to the user.

Collaborative filtering, a method used by RSs 21, makes predictions about the
interests of a user by collecting preferences from many users. This process
can be visualized and modeled using a bipartite graph, where one set of nodes
represents the users, and the other set represents the items. The task of the
RS then becomes a link prediction problem: based on the existing structure of
the graph and the preferences of similar users, which new user-item links are
most likely to form?

In this study, inspired by information filtering ideas, we propose a new way
of looking at chemical compound space and apply recommender system (RS)
principles to the search for new materials. We construct a bipartite graph to
depict the relations between periodic table elements and local crystal
environment geometries (ignorant of chemistry). By then generating an
embedding space from this graph, we are able to explore the affinities between
ions and local environments, and eventually spot new stable materials

In the remainer of this text we will use the expression anonymous motif sites
(AM sites) to refer to the local crystal environment geometries mentioned in
the previous paragraph. We discuss AM sites in detail in the Methods section.
Within our model, a ion-AMsite link means that that ion occupies that AM site
in one or more materials of the set of compounds used to construct the graph.
The weight assigned to each link corresponds to the respective compound’s
thermodynamic stability, as derived from Density Functional Theory (DFT). By
embedding this representation into a vector space, we were able to discern
chemical and structural patterns and develop a RS that recommends new ion-site
occupations. These recommendations may result in the formation of stable
compounds within a desired prototype crystal structure, as we demonstrate in
section Results and Discussion for the case of CsV3Sb522, 23.

Our methodology notably differs from previous machine learning and data-driven
methods that have generated element embeddings from the periodic table24, 25,
26, 27. Unlike those approaches, our embedding for ions is constructed
exclusively from the data of ion-site occupancies in known compounds, without
considering any compositional descriptor or chemical features from the
elements themselves. This unique approach emphasizes the incorporation of site
information within the embedding and, as a consequence, we cannot only delve
into ion-ion relationships within the embedding, but also exploit and analyse
the ion-site relationships.

##  2 Methods

###  2.1 The Link Prediction Task in Bipartite Graphs

Link prediction can be formally described as follows: consider a graph
$G=(V,E^{+})$, where $V$ is the set of nodes and $E^{+}$ the set of edges. In
the context of a recommender system, $G$ is a bipartite graph with $V=U\cup P$
and $U\cap P=\emptyset$. Here, $U$ represents the set of users and $P$ the set
of products, for example. Each edge connects a vertex in $U$ to a vertex in
$P$, with edges $(u,p)\in E^{+}$ representing user-product interactions.

The set of all possible but non-existing edges is denoted by $E^{-}$, such
that ${E^{+}\cap E^{-}=\emptyset}$. The link prediction task can then be
formulated as a binary classification problem, where all allowed edges or
pairs of nodes $(u,p)\in U\times P$ are associated with a label
$y_{up}\in\\{0,1\\}$, indicating the absence or presence of an edge between
$u$ and $p$, respectively.

The objective of the prediction task is to learn a function $f:U\times
P\rightarrow[0,1]$ that approximates the probability $P(y_{up}=1\mid u,p)$ of
an edge existing between nodes $u$ and $p$.

$P(y_{up}=1|u,p)\approx f(u,p)$ (1)

A common approach for link prediction employs node embeddings, that are vector
representations of nodes that capture their local graph structural
information. Let $\phi:V\rightarrow\mathbb{R}^{d}$ be a function mapping each
node $v\in V$ to a $d$-dimensional vector $\phi(v)$. This function outputs
vector representations for both $U$ and $P$ ($\boldsymbol{v}_{u}$ and
$\boldsymbol{v}_{p}$, respectively). The similarity between the embeddings of
two distinct nodes can be computed using a similarity function
$s:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, such as the dot
product or the cosine similarity. This can then be used as a proxy for the
likelihood of a link existing between them, or used as an input to a second
trained model $m$, such as a decision tree or a logistic regression,
responsible for defining the similarity threshold that delineates the existing
edges from the non-existing edges. The sequence of transformations from the
node representations to the final link prediction function $f$ is represented
in Figure 1.

Figure 1: Visualization of the link prediction task for recommender systems
($RS$), encapsulated by function $f$. The goal is to suggest new candidate
connections between a pair of nodes, represented in deep red and deep yellow.
The function $f$ is composed of two components. Firstly, the nodes are
processed by an embedding function $\phi$ which converts each node into a
vector representation. Then, these vectors are used as input for a similarity
function $s$. The output from the similarity function is then processed by a
model $m$, which finally produces the recommendation. Through this process,
the recommender system is capable of suggesting new links between nodes.

To train a link prediction model $f$, a training dataset
$\mathcal{D}={\\{(u_{i},p_{i},y_{u_{i}p_{i}})\\}}_{i=1}^{N}$ of $N$ instances
is constructed from positive examples $(u,p,y_{up}=1)$ drawn from the true
edges $E^{+}$ and negative examples $(u,v,y_{up}=0)$ drawn from the non-
existing edges $E^{-}$. The model parameters are optimized by minimizing a
global loss function $L$ that measures the discrepancy between the predicted
probabilities $f(u,p)$ and the true labels $y_{up}$. The edges from the set of
positive examples are removed from $G$ to avoid data leakage during the
optimization and training of $m$. The dataset $\mathcal{D}$ is partitioned
into training, validation, and test sets, ensuring a balanced distribution of
positive and negative samples. The global loss function can be written as

$L=\sum_{(u,p)\in\mathcal{D}}\mathcal{L}(f(u,p),y_{up})$ (2)

where $f(u,p)$ is the predicted probability of edge existence, $y_{up}$ is the
true label, and $\mathcal{L}$ is the loss function defined by the model $m$.

The learned function $f$ is evaluated on the test set to estimate its
generalization ability to unseen data. Ultimately, $f$ is used to estimate the
probability of all edges $(u,p)\in E^{-}$ for new recommendations.

###  2.2 Building the Recommender System

####  2.2.1 The Data:

We utilized two versions of the OQMD database: version 1.4 from 2020 and
version 1.5 from 2021, denoted as $\text{OQMD}^{\text{past}}$ and
$\text{OQMD}^{\text{current}}$, respectively. Each version was subjected to
the following filters:

  1. 1.

Non-unary materials with negative formation energy ($E_{f}<0.0\ \text{eV}$)
were selected. For unary materials, no restrictions were applied as the energy
above hull $E_{\text{hull}}$ is equivalent to $E_{f}$.

  2. 2.

Only unary, binary, ternary, and quaternary structures with 24 or fewer sites
in the unit cell were considered. This step was applied to try to eliminate
structures with a high number of low-symmetry non-equivalent occupied sites,
such as in amorphous phases.

  3. 3.

Duplicate entries were removed, i.e., data entries reported by OQMD as
repeated for the same compound (indicated by the column duplicate_of_id).

After applying these filters, $\text{OQMD}^{\text{past}}$ contained 307,912
entries and $\text{OQMD}^{\text{current}}$ contained 451,522 entries. The
energy above hull $E_{\text{hull}}$ and the crystal structure were extracted
for each compound in both databases.

Our first objective is to build a recommender system using
$\text{OQMD}^{\text{past}}$ and then validate it by verifying two different
things: (i) how much do the recommendations for ion-site occupations restrict
the search space; (ii) how many of the thermodynamically stable coumpounds
from $(\text{OQMD}^{\text{current}}-\text{OQMD}^{\text{past}})$ can be found
in this restricted search space. For simplicity, we carried out the latter
validation experiment on a subset of cubic halide perovskites, rather than
using the entire set of compounds $\mathcal{C}$. This process, often referred
to as ”historical validation”, also helped identify the optimal set of
hyperparameters for the model and the graph $G$ construction.

Following this validation step, we built a recommender system using
$\text{OQMD}^{\text{current}}$ data and generated recommendations for ion-site
occupations. These recommendations resulted in the proposal of new compounds
in a crystal structure prototype. Subsequently, we conducted Density
Functional Theory (DFT) calculations, using the same parameters of OQMD, and
utilized OQMD’s convex hull constructions to calculate the recommended
compounds’ $E_{\text{hull}}$. This allowed us to infer about the thermodynamic
stability of the new compounds and from it validate the RS recommendations.

####  2.2.2 Anonymous Motifs (AMs)

Here we present the definition of an Anonymous Motif (AM), since to the best
of our knowledge, a formal mathematical definition for AMs does not exist in
the literature of crystallography.

AM is defined solely by the geometric arrangement and connectivity of atomic
positions in space, disregarding the chemical identity of the atoms and,
consequently, the compound’s stoichiometry. Given a compound’s crystal
structure, we first disregard its composition and reduce it, when possible, to
a primitive cell $p$. Then, we use spglib28 to analyze the symmetry of the
resulting structure and of its sites. With the symmetry analysis data, we then
define AM as the following tuple:

$AM=(sg,n_{s},Wy,l)$ (3)

where:

  * •

$sg$ is the space group number (integer) of $p$

  * •

$n_{s}$ is the number of sites in $p$ (integer)

  * •

$Wy$ is a list of non-equivalent sites in $p$, given by their respective
Wyckoff positions, represented as tuples of the form $(n_{i},sym_{i})$, where
$n_{i}$ is the occupancy count of site $i$ (integer) and $sym_{i}$ is the
Hermann-Mauguin symbol for the point group symmetry of site $i$ (string) in
$p$.

  * •

$l$ is a subgroup number (integer), assigned after using pymatgen’s
StructureMatcher29 to further differentiate the $p$’s with the same
$(sg,N,Wy)$ but with different coordination environments. This is a ”dummy
index”, as it is used purely to differentiate otherwise identical AM labels
that represent different coordination environments.

As an illustrative example, we can consider the representation for cubic
perovskites: ${(221,5,[(1,\text{m-3m}),(1,\text{m-3m}),(3,4/\text{mmm})],0)}$.
The three labels in $Wy$ represent the cuboctahedral A-site, the octahedral
B-site, and the linear X-site, respectively (as shown in Figure 2. As the
chemical identities are disregarded, both simple cubic perovskites and double
perovskites would be grouped together into the same AM. By using AMs, one can
reduce the complexity of the materials space by identifying common
coordination patterns and geometrical motifs shared among different compounds.
They provide a standardized way of categorizing materials based on their
crystal structures and occupied sites.

Figure 2: Illustrative example of the Anonymous Motif label of cubic
perovskites. The elements of the $Wy$ list represents the A, B, and X sites.

####  2.2.3 AM sites

We define $\mathcal{C}$ as the set of all unique compounds $c$ from the OQMD
data. Let’s denote the set of all AMs found for all compounds in $\mathcal{C}$
as $AM^{\mathcal{C}}$. Now, we can define the set $S^{\mathcal{C}}$ of all
sites of all AMs in $AM^{\mathcal{C}}$.

$S^{\mathcal{C}}=\\{(sg,n_{s},n_{i},sym_{i},l)\mid\forall AM\in
AM^{\mathcal{C}},(n_{i},sym_{i})\in Wy\\}$ (4)

As an example, the AM site $(227,2,2,-43m,0)$ groups together all sites of
silicon, zincblende, chalcopyrite, stannite and kesterite structures, as shown
in detail in Figure 3. For statistical reasons, we only considered AMs with
more than 5 compounds.

Figure 3: AM site grouping example. As only the lattice geometry is taken into
account by ignoring compositions and stoichiometries, the occupied sites of
all the presented compounds get mapped to the same (227, 2, 2, -43m, 0) AM
site.

####  2.2.4 Defining the Bipartite Graph $G^{T}$

We denote the set of all ions in the chemical compositions of all compounds in
$\mathcal{C}$ as $I$. The bipartite graph is constructed with vertices divided
into two disjoint sets, $I$ and $S^{\mathcal{C}}$, where each edge connects a
vertex in $I$ to one in $S^{\mathcal{C}}$. The graph is represented as
$G^{T}=(I,S^{\mathcal{C}},E^{+},W^{T})$, where $E^{+}$ is the set of existing
edges, and $W^{T}:E^{+}\rightarrow\mathbb{R}$ is a Boltzmann-like weight
function defining the edge weights.

The link $(i,s)$ is added to $E^{+}$ if the occupation of an AM site $s$ by an
ion $i$ occurs in a compound $c$. The set of all links is thus defined as:

$E^{+}=\\{(i,s)\mid i\text{ occupies }s\text{ in compound }c,\forall
c\in\mathcal{C}\\}$ (5)

A weight $w^{c,T}(i,s)$ is defined as being proportional to the energy above
the hull of the compound $c$, $E_{\text{hull}}^{c}$, using a Boltzmann factor
and introducing a parametric temperature $T$:

$w^{c,T}(i,s)=\exp(-E_{\text{hull}}^{c}/kT)$ (6)

We created three graphs $G^{T}$ with three parameteric temperatures: $0\
\text{K}$, $100\ \text{K}$, and $300\ \text{K}$. The objective was to define a
threshold for compounds $c\in\mathcal{C}$: a 0K parametric temperature
includes only compounds in the convex hull (here defined as the compounds with
$E_{\text{hull}}\leq 30\text{ meV/atom}$), while T $>$ 0 K parametric
temperatures attribute a small weight even for compounds with large
$E_{\text{hull}}$. We wanted to investigate if adding these compounds into the
model would improve the recommendations. Larger temperatures penalizes less
the weights of the graph links associated with compounds that are
thermodynamically unstable. Figure 4 illustrates this concept on a
hypothetical binary convex hull.

Figure 4: Illustration of the effect of the temperature $T$ on the selection
of the compounds $\mathcal{C}$ used for building $G^{T}$, and on the weights
$W^{T}(i,s)$ between the ion $i$ and the site $s$.

Each edge’s total weight $W^{T}(i,s)$ is computed iteratively over all
compounds $c^{\prime}$ of the set of compounds $\mathcal{C}^{\prime}$, where
$\mathcal{C}^{\prime}$ is the subset of all compounds that are formed by the
occupation of $s$ by $i$:

$W^{T}(i,s)=\sum_{c^{\prime}\in\mathcal{C}^{\prime}}w^{c^{\prime},T}(i,s)$ (7)

We truncated the values of $w^{c^{\prime},T}$: if
$w^{c^{\prime},T}(i,s)<0.001$, then we considered that $(i,s)$ is non-
existent.

This workflow assigns larger weights $W^{T}$ to common occupations, such as
the X-site in cubic perovskites ${(\text{O},(221,5,3,(3,4/\text{mmm}),0))}$
that occur in numerous oxide cubic perovskites, per comparison with the less
common occupations like ${(\text{Ca},(221,5,3,(3,4/\text{mmm}),0))}$, found in
the anti-perovskite $\text{Ca}_{3}\text{SnO}$. The weight $W^{T}(i,s)$
quantifies the frequency of the $(i,s)$ occupation in $\mathcal{C}$, weighted
by the thermodynamic stability of the compounds featuring $(i,s)$.

Figure 5 pictorially illustrates the structure of the bipartite graph $G^{T}$
built from the set of compounds $\mathcal{C}$.

Figure 5: Illustration of the bipartite graph
$G^{T}=(I,S^{\mathcal{C}},E^{+},W^{T}$)

.

####  2.2.5 The Embedding for Ions and AM sites

In order to generate embeddings that are encoded into vector representations
for each type of node (ions and AM sites) we conduct random walks across the
graph. DeepWalk 30, a widely-used algorithm, generates these random walks by
starting at a randomly selected node and choosing the next node uniformly at
random from its neighbors. This process is repeated for a fixed number of
steps to generate a single random walk, and is performed multiple times to
generate a set of random walks for the entire graph. The premise of DeepWalk
is that nodes frequently co-occurring on these walks likely have similar roles
or positions in the graph, hence the learned embeddings will capture these
similarities. Once the random walks are generated, they are input into the
Word2Vec 31 algorithm to learn the embeddings. These embeddings, which capture
the local and global structure of the graph, can then be used for downstream
tasks such as link prediction or node classification.

This methodology was adapted for a weighted graph. A random walk $rw$ in the
weighted bipartite graph $G$ can be defined as a sequence of alternating
vertices from the sets $I$ and $S^{\mathcal{C}}$, where the probability of
transitioning from one vertex to another depends on the normalized weights of
the edges $W^{T}(i,s)$ linking to the current graph node.

First, we calculate the transition probabilities for each ion and AM site at
parametric temperature $T$. For each ion $i\in I$ and each AM site $s\in
S^{\mathcal{C}}$, the sum of weights $W^{T}(i,s)$ for all edges connected to
it is computed. Then, the transition probability $P(i,s)$ is given as follows:

$P(i,s)=\frac{W^{T}(i,s)}{\sum_{s\in S^{\mathcal{C}}}W^{T}(i,s)}$ (8)

A single $rw$’s in $G$ can be expressed as a sequence of alternating vertices
from $I$ and $S$, $rw=(i_{1},s_{1},i_{2},s_{2},\dots,i_{n},s_{n})$, where
$i_{k}\in I$, $s_{k}\in S^{\mathcal{C}}$, and the choice of $s_{k}$ given
$i_{k}$ is determined by the transition probability $P(i_{k},s_{k})$, and the
choice of $i_{k+1}$ given $s_{k}$ is determined by the transition probability
$P(s_{k},i_{k+1})$. A total of 150 random walks per node, with a walk length
of 50 nodes was used.

The set of all $rw$ in the graph is defined as $RW$. The Word2Vec function
takes $RW$ as input and a set of hyperparameters $\theta$ and outputs an
embedding function $\phi_{\theta}=\text{Word2Vec}(RW,\theta)$.

The Figure 6a illustrates the embedding creation process.

Figure 6: a) Graphical representation of the embedding creation workflow:
given a graph $G^{T}$, the set of weighted random walks is made and then used
as input for Word2Vec. The result is an embedding function $\phi_{t}heta$. b)
Results obtained from the embedding’s vector representations of ions and AM
sites.

####  2.2.6 Defining the function $f$

Our RS is encoded in the function $f(i,s)$. The embedding function
$\phi_{\theta}$ maps a node (either ion $i\in I$ or a AM site $s\in S$) to a
$d$-dimensional vector representation. These representations are then
processed by a similarity function, with the output fed into a binary
classification model. We employ cosine similarity as the similarity metric and
a decision tree (DT) as the model.

For training, we adapted the dataset $\mathcal{D}$ creation methodology for
weighted graphs. Edges $E^{+}$ were sampled by their weight, with higher
weight edges having a greater selection chance. Once selected, the edge was
removed from the original $G^{T}$ graph to prevent data leakage during
embedding. The edge’s weight was then used as a multiplier for instance
repetition in the dataset. As an example, a positive edge with a weight of 10
would be repeated 10 times. Negative examples were uniformly sampled from
$E^{-}$ until an approximate 50% ratio of positive and negative examples was
achieved.

The dataset was divided into 80% for training, 10% for validation, and 10% for
testing. Optimal hyperparameters for each temperature were determined using
the Optuna python library 32. We optimized the hyperparameters $\theta$ of the
embedding function $\phi$, as these have been shown to significantly impact
recommender systems 33. The accuracy results for the optimal hyperparameters
are presented in the Supplementary Information.

After training, the learned function $f$ was used to predict the likelihood of
edges $(i,s)\in E^{-}$, between node pairs in the graph, and then use these
new reccomended occupations for building new compounds.

##  3 Results and Discussion

Establishing an embedding that encapsulates the structure of bipartite graphs
$G^{T}$ facilitated two distinct analyses concerning the relations among the
graph entities. The first one pertains to the relationships between ions,
while the second one pertains to the relationship between the ions and the AM
sites.

The first was done by correlating the ion-ion embedding cosine similarities to
ionic substitution in crystal structures. This approach enabled the validation
of chemical trends within the continuous embedding space created solely with
ion-site occupation data.

The second analysis involved the study of proximate AM sites surrounding the
ions within the embedding space. This approach produced two distinct types of
outcomes. Firstly, it enabled the verification of the overarching local
geometries of the AM sites in association with the corresponding ion.
Secondly, it allowed the exploration of novel ion-site occupancies informed by
the trained RS function $f$.

Figure 6b illustrates the results obtained from the embedding and further
discussed below.

###  3.1 Embedding Chemical Trends

After optimizing the hyperparameters of Word2Vec, an embedding was obtained
from $G^{0}$, exclusively considering compounds on the hull and built with
$\text{OQMD}^{\text{current}}$. We applied the uniform manifold approximation
and projection (UMAP) technique to reduce the dimensionality of the embedding
34 for a two-dimensional visualization (Figure 7).

Figure 7: Embedding of $G^{0}$ built with $\text{OQMD}^{\text{current}}$. In
light gray are the AM sites. The colored points represent the set of 89 ions,
each color representing its corresponding family in the periodic table. UMAP
parameters: metric: cosine, n_neighbours: 250, min_dist: 0.5, densmap: True,
random_state: 0.

The embedding shows that, with some few distinct exceptions, the ions from the
same family in the periodic table usually occupy the same group of sites, as
they are close together in the embedding space projection of Figure 7. Within
the same cluster family, there is also a larger proximity between ions with
similar atomic radii. For instance, within the halogen family cluster,
Chlorine (Cl) is more closely related to Bromine (Br) than to Iodine (I).

A similar pattern is also observed within transition metal clusters. Ions from
the 5d period display a greater likeness to 4d period ions than those from the
3d period. This can be exemplified by 5d transition metals (Hf, Ta, W, Re, Os,
Ir, Pt), which showcase similar atomic radii to 4d metals within the same
group (Zr, Nb, Mo, Tc, Ru, Rh, Pd). This phenomena can be explained by the
lanthanide contraction effect35. As one moves through the lanthanide series
(f-block elements) on the periodic table, the atomic radii tend to decrease
due to the poor shielding of the increasing nuclear charge by the 4f
electrons, resulting in a higher effective nuclear charge and a smaller atomic
radius. This effect continues into the subsequent d-block elements, causing
the 5d elements to have atomic radii similar to their 4d counterparts, despite
being in a new period.

The embedding projection also illuminates trends related to the ionic
character of ions. Anions, such as oxygen, nitrogen, carbon, hydrogen,
halogens, and chalcogens, largely congregate on the left side, while cations
typically occupy the right side. This pattern supports the notion that cations
and anions tend to inhabit distinct sites in crystal structures.

The visualization also underscores a sizable and dense cluster surrounding
oxygen (O). This clustering likely stems from the original ICSD 36 database
that serves as the foundation for OQMD. ICSD primarily focuses on published
crystallographic research, which heavily features oxides due to their
prevalence and to their wide-ranging applications.

A unique trend is apparent with ions in the second period, or the first row of
p-block elements (B, C, N, O, F), as these elements behave differently from
their same-family heavier counterparts. This divergence arises from their
increased electronegativity and smaller atomic radii, affecting bond polarity
and lengths.

Certain ion-ion similarities spanning across different periodic families can
be observed, often due to shared oxidation states. For instance, Ytterbium
(Yb) and Europium (Eu) from the lanthanide series show a striking resemblance
to alkaline earth metals, likely attributable to their shared +2 oxidation
state. Thallium (Tl), despite belonging to a different group, demonstrates
notable proximity to Sodium (Na) in the embedding, an occurrence that can be
explained by Thallium’s ability to adopt a +1 oxidation state. In addition,
Scandium (Sc) and Yttrium (Y), which typically present a +3 oxidation state,
are found to be closely aligned with the lanthanides, further emphasizing the
influence of overlapping oxidation states on ion-ion relationships.

Supplementary to the UMAP visualization, a dendrogram depicting the
hierarchical clustering of the 89 ions within the embedding space is provided
in the Supplementary Information.

###  3.2 Embedding Local Geometry Trends

In further analysis, we generated heatmaps to provide insights into the local
geometric configurations of the sites. Utilizing the CrystalNN method37, 38,
we computed scores representing the likelihood of each type of local geometry
for every site. Subsequently, we divided the 2D UMAP space and computed the
average score for each segment. To smooth out the data distribution, we
applied a Gaussian fitting across the 2D space, focusing on the dimension
representing the geometry of interest.

The results for two distinct types of local geometries—linear and square
coplanar—are presented in Figure 8. The color variations surrounding the ions
in these maps correspond to the degree to which ions occupy stable compound’s
sites with the given geometry.

Figure 8: Heatmaps indicating local geometry coordination for linear and
square coplanar configurations. The color scale corresponds to the density of
sites within the given area of the UMAP 2D-transformed embedding space that
conform to the respective geometry.

Something important to note is that, as the color encodes a site average over
the area, a light color does not necessarily imply that the ion is not capable
of occupying a site with the given local geometry. It implies instead that,
statistically in the data set, that geometry is less common than other
geometries.

The heatmap for linear local geometry shows that sites with this geometry tend
to concentrate around anions such as halogens, H and O. Cu and Ag also present
a certain concentration around them, which could be attributed to their
capacity of having a less common oxidation state of 1+. An example of such
occupation is found in the delaffosite class of compounds 39.

The square planar coordination is not as common as other geometries but is
typically found with certain transition metals that present a 2+ oxidation
state. For instance, in certain crystal structures, Pd, Pt, and Cu can adopt
square planar geometry. The heatmap confirms this and suggests that Ir and Zn,
which can also present 2+ oxidation state, can also adopt this geometry. One
prominent such family is that of the copper oxide high-temperature
superconductors 40, 41, where the square planar geometry and the Cu 2+
oxidation state are believed to be two crucial ingredients for the distinctive
properties of these materials

The heatmaps for other local geometries—cuboctahedral, octahedral,
tetrahedral, trigonal planar, and single bound—are presented in the
Supplementary Information.

###  3.3 Ion-site Recommendations

####  3.3.1 Historical Validation on Cubic Halide Perovskites

Utilizing the methodology outlined in the previous section, we trained a RS
using $\mathcal{C}=\text{OQMD}^{\text{past}}$. We then validated the
methodology using the set of new cubic halide perovskites in
$(\text{OQMD}^{\text{current}}-\text{OQMD}^{\text{past}})$. Our objective was
to assess the impact of the parametric temperature $T$ on two key aspects of
RSs: the degree to which it restricts the search space; and the number of
thermodynamically stable compounds from
$(\text{OQMD}^{\text{current}}-\text{OQMD}^{\text{past}})$ that can be
identified within this restricted search space by our framework.

Halide perovskites are important since they are a promising class of materials
for optoelectronics and photovoltaics.42 They possess unique crystal
structures and exceptional optoelectronic properties, including high
absorption coefficients and long charge carrier diffusion lengths. Halide
perovskite solar cells have achieved impressive power conversion efficiencies,
surpassing traditional thin-film technologies.43 These materials can be
fabricated using low-cost techniques and integrated into various devices such
as LEDs and photodetectors. However, challenges remain in their stability, and
ongoing research aims to improve their performance and understand their
fundamental properties.

A comparison to a brute-force approach needs to be done when evaluating a RS.
A cubic halide perovskites brute-force search has a space size of 11,200
possibilities, provided we consider four halogens {F, Cl, Br, I}, only non-
metalic ions, and ions from the 1st and 2nd row for the A and B sites, while
excluding noble gases, radioactive ions Additionally, for the A site, we
eliminated the option of p-orbital ions.

$\text{OQMD}^{\text{past}}$ has a total of 135 cubic halide perovskites, 63 of
which are thermodynamically stable ($E_{\text{hull}}\leq 30\text{meV/atom}$).
$\text{OQMD}^{\text{current}}$ has 432 additional cubic halide perovskites, 50
of which are thermodynamically stable. In the context of a binary
classification problem, these 50 compounds constitute the set of positive
instances, while the remaining 382 instances are defined as the negative class
of new thermodynamically unstable compounds. Since the positive and negative
sets are unbalanced and false positives are undesired in the context of a RS,
the precision was the chosen metric to evaluate this task. Precision is a
measure of relevancy — it quantifies the proportion of recommended items that
are actually relevant to the user. The results for parametric temperatures of
$T=\\{0K,100K,300K\\}$ are presented in Table 1.

Table 1: Comparison of different models with different temperatures to the brute-force approach. The search space is the number of possible combinations given the ion-site suggestions made by the RS. The precision is given within $(\text{OQMD}^{\text{current}}-\text{OQMD}^{\text{past}})$ by the ratio of the thermodynamically stable compounds (true positives) to the total number of compounds suggested by the recommender (true positives + false positives). | $G^{0}$ | $G^{100}$ | $G^{300}$ | Brute-force
---|---|---|---|---
A search space | 9 | 9 | 12 | 50
B search space | 11 | 12 | 15 | 57
ABX3 search space | 392 (3.5%) | 428 (3.8%) | 708 (6.3%) | 11,200 (100%)
In $(\text{OQMD}^{\text{current}}-\text{OQMD}^{\text{past}})$ | 49 | 73 | 109 | 432
$E_{\text{hull}}\leq 30\text{ meV/atom}$ | 13 | 23 | 28 | 50
Precision | 26.5% | 31.5% | 25.7% | 11.6%

The results indicate that $G^{100}$ provides the optimal balance between
restricting the search space and offering the largest proportion of
thermodynamically stable recommendations. These findings suggest that
including all compounds and weighting them by their $E_{\text{hull}}$ is the
most effective strategy. The comparatively poorer performance of $T=300K$
suggests that a higher temperature underpenalizes non-stable compounds.

Another trend that can be attributed to the parametric temperature is on the
ABX3 search space reduction: a higher $T$ resulted in a higher number of
possibilites. Such result can be attributed to the higher connectivity of ions
to AM sites in the resulting $G^{T}$ graphs, resulting from $W^{T}(i,s)$.

It is worth noting that there are compounds that were recommended that are not
present in $(\text{OQMD}^{\text{current}}-\text{OQMD}^{\text{past}})$ and thus
remain to be verified with DFT calculations. The list of the recommendations
for A and B sites are given in the Supplementary Information.

####  3.3.2 Recommendations for Kagome Lattice Prototype Structure

Using the validated methodology, optimal hyperparameters $\theta$, and the
best parametric temperature $T=100K$, we trained a recommender system with
$\text{OQMD}^{\text{current}}$. This was then used to recommend ion-site
occupations for the sites of the class of compounds AB3C5.

Their crystal structure presents a layered structure of Kagome lattices
separated by A-ions in a honeycomb lattice. The Kagome lattice is made of
three B-ions coordinated with one C-ion, sandwiched between two honeycomb
lattices of C. The structure presents a symmetry given by the space group
$P6/mmm$.

Ortiz et al first synthesized this prototype structure in 2019 with the
compositions KV3Sb5 , RbV3Sb5, CsV3Sb523. Since then, many characterizations
of these compounds were made, revealing properties such as chiral charge
ordering 44, large anomalous Hall effect 45, non-trivial topological band
structures and superconducting ground states 46, 22.

Figure 9: Prototype Crystal structure of compound CsV3Sb5, with stoichiometry
AB3C5. The ions A and B are represented in blue and red, respectively. The C
ion occupies two non-equivalent types of sites, represented in green and in
yellow.

The prototype crystal structure maps to the four non-equivalent sites
$(191,9,[(1,6/mmm),(1,6/mmm),(3,mmm),(4,3m)],0)$ AM label, as shown in Figure
9. In this prototype, A occupies the first site in the $Wy$ list $(1,6/mmm)$,
B occupies the third site $(3,mmm)$, and C occupies both the second and the
fourth sites $(1,6/mmm)$ and $(4,3m)$. For simplicity, we truncated the
recommendations to a AB3C5 stoichiometry, using the intersection of
recommendations of C for both sites (represented by yellow and green atoms in
Figure 9).

We input these sites into the RS, which recommended ions to occupy them. Known
compositions for this structure in $\text{OQMD}^{\text{current}}$ were
filtered out, and the remaining recommendations were ranked according to the
sum of ion-site distances in the embedding space. The top 10 recommendations
were then verified using DFT calculations, and the results are presented in
Table 2. All ion-site recommendations are given in the Supplementary
Information.

Table 2: Thermodynamical stability results for the top 10 compounds built with the ion-site suggestions of the RS built with T = 100K and $\text{OQMD}^{\text{current}}$. A | B | C | Composition |  $\Sigma$(ion-site distances) |  $E_{\text{hull}}$ (meV/atom)
---|---|---|---|---|---
Rb | Ti | Sb |  RbTi3Sb5 | 0.125 | 0
K | Ti | Sb |  KTi3Sb5 | 0.128 | 13
Cs | Ti | Sb |  CsTi3Sb5 | 0.134 | 0
Rb | Ti | Bi |  RbTi3Bi5 | 0.173 | 5
K | Ti | Bi |  KTi3Bi5 | 0.176 | 23
Cs | Ti | Bi |  CsTi3Bi5 | 0.182 | 0
Rb | Hf | Sb |  RbHf3Sb5 | 0.219 | 182
K | Hf | Sb |  KHf3Sb5 | 0.223 | 165
Tl | Ti | Sb |  TlTi3Sb5 | 0.225 | 102
Cs | Hf | Sb |  CsHf3Sb5 | 0.228 | 139

The results show that 6 out of the 10 recommendations are indeed stable (with
$E_{\text{hull}}<30\text{ meV/atom}$). These 6 compounds were also reported by
Jiang et al 47 as having a $E_{\text{hull}}<30\text{ meV/atom}$, but at the
cost of a total of 1386 high-throughput DFT calculations performed. Notably,
these materials were entirely new in $\text{OQMD}^{\text{current}}$.

Another significant finding is the novel recommendation of the ion family Ti,
Hf, and Zr for the B site occupancy. This type of occupation was not
previously present in $\text{OQMD}^{\text{current}}$, making it a non-trivial
recommendation that could not have been made simply by visually or intuitively
filling the gaps in the periodic table.

##  4 Conclusion

Leveraging data from OQMD, we successfully constructed a graph-based
recommender system that incorporates information on the thermodynamic
stability of compounds along with symmetry and local geometry of sites. The
developed embedding enabled us to discern chemical and local geometrical
trends through ion-ion similarities, and to extract ion-site recommendations
from the ion-site relationships.

Our experiments with a parametric temperature that weights the
$E_{\text{hull}}$ of compounds led us to conclude that considering not only
the compounds on the hull, but also those above the hull, can be beneficial to
the identification of new compounds.

By employing a recommender system with optimally tuned hyperparameters and
parametric temperature, we managed to propose novel compounds that are not
present in the current version of OQMD. This demonstrates the potential of our
approach to facilitate the discovery of new materials.

##  5 Data Availability

The RS trained with $\text{OQMD}^{\text{current}}$, along with the graphs
$G^{100}$ and $G^{0}$ (as the data used to build them), as well as the ion and
site embeddings, can be downloaded from our GitHub repository at
https://github.com/simcomat/ionic-sub-RS. A detailed explanation of the
process for obtaining recommendations for ion-site occupations and new
compounds is provided in a Jupyter notebook available on the same page.

##  6 Competing Interests

The Authors declare no Competing Financial or Non-Financial Interests.

##  7 Author Contributions

E.O. conceived the work, built the model, and ran the DFT calculations and
wrote the first version of the paper. H.F., J.N.B.R. and G.M.D. helped with
refinements of the model. G.M.D. supervised the work. All authors contributed
with the final version of the paper.

The authors acknowledge the National Laboratory for Scientific Computing
(LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer.
The authors also thank the Brazilian funding agencies FAPESP (grants
17/02317-2 and 18/11641-0), CNPq (309384/2020-6) and CAPES for financial
support. This research has been conducted within the scope of the Materials
Informatics INCT (CNPq).

## References

  * Davies et al. 2016 Davies, D. W.; Butler, K. T.; Jackson, A. J.; Morris, A.; Frost, J. M.; Skelton, J. M.; Walsh, A. Computational screening of all stoichiometric inorganic materials. _Chem_ 2016, _1_ , 617–627
  * Oganov and Glass 2006 Oganov, A. R.; Glass, C. W. Crystal structure prediction using ab initio evolutionary techniques: Principles and applications. _The Journal of chemical physics_ 2006, _124_ , 244704
  * Zunger 2018 Zunger, A. Inverse design in search of materials with target functionalities. _Nature Reviews Chemistry_ 2018, _2_
  * Goldschmidt 1926 Goldschmidt, V. M. Die gesetze der krystallochemie. _Naturwissenschaften_ 1926, _14_ , 477–485
  * Callister and Rethwisch 2018 Callister, W. D.; Rethwisch, D. G. _Materials science and engineering: An introduction_ ; Wiley, 2018; p 96
  * Saal et al. 2013 Saal, J. E.; Kirklin, S.; Aykol, M.; Meredig, B.; Wolverton, C. Materials Design and Discovery with High-Throughput Density Functional Theory: The Open Quantum Materials Database (OQMD). _JOM_ 2013, _65_ , 1501–1509
  * Curtarolo et al. 2012 Curtarolo, S.; Setyawan, W.; Wang, S.; Xue, J.; Yang, K.; Taylor, R. H.; Nelson, L. J.; Hart, G. L.; Sanvito, S.; Buongiorno-Nardelli, M.; Mingo, N.; Levy, O. AFLOWLIB.ORG: A distributed materials properties repository from high-throughput ab initio calculations. _Computational Materials Science_ 2012, _58_ , 227–235
  * Jain et al. 2013 Jain, A.; Ong, S. P.; Hautier, G.; Chen, W.; Richards, W. D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G.; Persson, K. A. Commentary: The Materials Project: A materials genome approach to accelerating materials innovation. _APL Materials_ 2013, _1_ , 011002
  * Draxl and Scheffler 2019 Draxl, C.; Scheffler, M. The NOMAD laboratory: from data sharing to artificial intelligence. _Journal of Physics: Materials_ 2019, _2_ , 036001
  * Bartel et al. 2019 Bartel, C. J.; Sutton, C.; Goldsmith, B. R.; Ouyang, R.; Musgrave, C. B.; Ghiringhelli, L. M.; Scheffler, M. New tolerance factor to predict the stability of perovskite oxides and halides. _Science Advances_ 2019, _5_ , 1–10
  * Pei et al. 2020 Pei, Z.; Yin, J.; Hawk, J. A.; Alman, D. E.; Gao, M. C. Machine-learning informed prediction of high-entropy solid solution formation: Beyond the Hume-Rothery rules. _npj Computational Materials_ 2020, _6_ , 50
  * Montoya et al. 2020 Montoya, J. H.; Winther, K. T.; Flores, R. A.; Bligaard, T.; Hummelshøj, J. S.; Aykol, M. Autonomous intelligent agents for accelerated materials discovery. _Chemical Science_ 2020, _11_ , 8517–8532
  * Lyakhov et al. 2013 Lyakhov, A. O.; Oganov, A. R.; Stokes, H. T.; Zhu, Q. New developments in evolutionary structure prediction algorithm USPEX. _Computer Physics Communications_ 2013, _184_ , 1172–1182
  * Wang et al. 2021 Wang, H.-C.; Botti, S.; Marques, M. A. Predicting stable crystalline compounds using chemical similarity. _npj Computational Materials_ 2021, _7_ , 12
  * Seko et al. 2018 Seko, A.; Hayashi, H.; Kashima, H.; Tanaka, I. Matrix-and tensor-based recommender systems for the discovery of currently unknown inorganic compounds. _Physical Review Materials_ 2018, _2_ , 013805
  * Hayashi et al. 2022 Hayashi, H.; Seko, A.; Tanaka, I. Recommender system for discovery of inorganic compounds. _npj Computational Materials_ 2022, _8_ , 217
  * Reiser et al. 2022 Reiser, P.; Neubert, M.; Eberhard, A.; Torresi, L.; Zhou, C.; Shao, C.; Metni, H.; van Hoesel, C.; Schopmans, H.; Sommer, T., et al. Graph neural networks for materials science and chemistry. _Communications Materials_ 2022, _3_ , 93
  * Chen and Ong 2022 Chen, C.; Ong, S. P. A universal graph deep learning interatomic potential for the periodic table. _Nature Computational Science_ 2022, _2_ , 718–728
  * Choudhary et al. 2023 Choudhary, K.; DeCost, B.; Major, L.; Butler, K.; Thiyagalingam, J.; Tavazza, F. Unified graph neural network force-field for the periodic table: solid state applications. _Digital Discovery_ 2023, _2_ , 346–355
  * Ko et al. 2022 Ko, H.; Lee, S.; Park, Y.; Choi, A. A survey of recommendation systems: recommendation models, techniques, and application fields. _Electronics_ 2022, _11_ , 141
  * Su and Khoshgoftaar 2009 Su, X.; Khoshgoftaar, T. M. A survey of collaborative filtering techniques. _Advances in artificial intelligence_ 2009, _2009_
  * Ortiz et al. 2020 Ortiz, B. R.; Teicher, S. M.; Hu, Y.; Zuo, J. L.; Sarte, P. M.; Schueller, E. C.; Abeykoon, A. M.; Krogstad, M. J.; Rosenkranz, S.; Osborn, R., et al. Cs v 3 sb 5: A z 2 topological kagome metal with a superconducting ground state. _Physical Review Letters_ 2020, _125_ , 247002
  * Ortiz et al. 2019 Ortiz, B. R.; Gomes, L. C.; Morey, J. R.; Winiarski, M.; Bordelon, M.; Mangum, J. S.; Oswald, I. W.; Rodriguez-Rivera, J. A.; Neilson, J. R.; Wilson, S. D., et al. New kagome prototype materials: discovery of KV 3 Sb 5, RbV 3 Sb 5, and CsV 3 Sb 5. _Physical Review Materials_ 2019, _3_ , 094407
  * Jha et al. 2018 Jha, D.; Ward, L.; Paul, A.; Liao, W.-k.; Choudhary, A.; Wolverton, C.; Agrawal, A. Elemnet: Deep learning the chemistry of materials from only elemental composition. _Scientific reports_ 2018, _8_ , 1–13
  * Zhou et al. 2018 Zhou, Q.; Tang, P.; Liu, S.; Pan, J.; Yan, Q.; Zhang, S.-C. Learning atoms for materials discovery. _Proceedings of the National Academy of Sciences_ 2018, _115_ , E6411–E6417
  * Schütt et al. 2018 Schütt, K. T.; Sauceda, H. E.; Kindermans, P.-J.; Tkatchenko, A.; Müller, K.-R. Schnet–a deep learning architecture for molecules and materials. _The Journal of Chemical Physics_ 2018, _148_ , 241722
  * Antunes et al. 2022 Antunes, L. M.; Grau-Crespo, R.; Butler, K. T. Distributed representations of atoms and materials for machine learning. _npj Computational Materials_ 2022, _8_ , 44
  * Togo and Tanaka 2018 Togo, A.; Tanaka, I. Spglib: a software library for crystal symmetry search. _arXiv preprint arXiv:1808.01590_ 2018,
  * Ong et al. 2013 Ong, S. P.; Richards, W. D.; Jain, A.; Hautier, G.; Kocher, M.; Cholia, S.; Gunter, D.; Chevrier, V. L.; Persson, K. A.; Ceder, G. Python Materials Genomics (pymatgen): A robust, open-source python library for materials analysis. _Computational Materials Science_ 2013, _68_ , 314–319
  * Perozzi et al. 2014 Perozzi, B.; Al-Rfou, R.; Skiena, S. Deepwalk: Online learning of social representations. Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining. 2014; pp 701–710
  * Mikolov et al. 2013 Mikolov, T.; Chen, K.; Corrado, G.; Dean, J. Efficient estimation of word representations in vector space. _arXiv preprint arXiv:1301.3781_ 2013,
  * Akiba et al. 2019 Akiba, T.; Sano, S.; Yanase, T.; Ohta, T.; Koyama, M. Optuna: A next-generation hyperparameter optimization framework. Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining. 2019; pp 2623–2631
  * Caselles-Dupré et al. 2018 Caselles-Dupré, H.; Lesaint, F.; Royo-Letelier, J. Word2vec applied to recommendation: Hyperparameters matter. Proceedings of the 12th ACM Conference on Recommender Systems. 2018; pp 352–356
  * McInnes et al. 2018 McInnes, L.; Healy, J.; Melville, J. Umap: Uniform manifold approximation and projection for dimension reduction. _arXiv preprint arXiv:1802.03426_ 2018,
  * Housecroft and Sharpe 2018 Housecroft, C. E.; Sharpe, A. G. _Inorganic Chemistry_ ; Pearson, 2018; pp 809, 1036
  * Hellenbrandt 2004 Hellenbrandt, M. The Inorganic Crystal Structure Database (ICSD)—Present and Future. _Crystallography Reviews_ 2004, _10_ , 17–22
  * Zimmermann and Jain 2020 Zimmermann, N. E.; Jain, A. Local structure order parameters and site fingerprints for quantification of coordination environment and crystal structure similarity. _RSC Advances_ 2020, _10_ , 6063–6081
  * Pan et al. 2021 Pan, H.; Ganose, A. M.; Horton, M.; Aykol, M.; Persson, K. A.; Zimmermann, N. E.; Jain, A. Benchmarking coordination number prediction algorithms on inorganic crystal structures. _Inorganic chemistry_ 2021, _60_ , 1590–1603
  * Shannon et al. 1971 Shannon, R. D.; Prewitt, C. T.; Rogers, D. B. Chemistry of noble metal oxides. II. Crystal structures of platinum cobalt dioxide, palladium cobalt dioxide, coppper iron dioxide, and silver iron dioxide. _Inorganic Chemistry_ 1971, _10_ , 719–723
  * Bednorz and Müller 1986 Bednorz, J. G.; Müller, K. A. Possible highTc superconductivity in the Ba-La-Cu-O system. _Zeitschrift für Physik B Condensed Matter_ 1986, _64_ , 189–193
  * Proust and Taillefer 2019 Proust, C.; Taillefer, L. The Remarkable Underlying Ground States of Cuprate Superconductors. _Annual Review of Condensed Matter Physics_ 2019, _10_ , 409–429
  * Walsh 2021 Walsh, A. Atomistic models of metal halide perovskites. _Matter_ 2021, _4_ , 3867–3873
  * Tian et al. 2020 Tian, J.; Xue, Q.; Yao, Q.; Li, N.; Brabec, C. J.; Yip, H.-L. Inorganic Halide Perovskite Solar Cells: Progress and Challenges. _Advanced Energy Materials_ 2020, _10_ , 2000183
  * Jiang et al. 2021 Jiang, Y.-X.; Yin, J.-X.; Denner, M. M.; Shumiya, N.; Ortiz, B. R.; Xu, G.; Guguchia, Z.; He, J.; Hossain, M. S.; Liu, X., et al. Unconventional chiral charge order in kagome superconductor KV3Sb5. _Nature materials_ 2021, _20_ , 1353–1357
  * Yang et al. 2020 Yang, S.-Y.; Wang, Y.; Ortiz, B. R.; Liu, D.; Gayles, J.; Derunova, E.; Gonzalez-Hernandez, R.; Šmejkal, L.; Chen, Y.; Parkin, S. S., et al. Giant, unconventional anomalous Hall effect in the metallic frustrated magnet candidate, KV3Sb5. _Science advances_ 2020, _6_ , eabb6003
  * Ortiz et al. 2021 Ortiz, B. R.; Sarte, P. M.; Kenney, E. M.; Graf, M. J.; Teicher, S. M.; Seshadri, R.; Wilson, S. D. Superconductivity in the Z 2 kagome metal KV 3 Sb 5\. _Physical Review Materials_ 2021, _5_ , 034801
  * Jiang et al. 2022 Jiang, Y.; Yu, Z.; Wang, Y.; Lu, T.; Meng, S.; Jiang, K.; Liu, M. Screening promising CsV3Sb5-like kagome materials from systematic first-principles evaluation. _Chinese Physics Letters_ 2022, _39_ , 047402


# Absence of gravitational polarization mechanism
in the canonical bimetric theory

Luc Blanchet<EMAIL_ADDRESS>$\mathcal{G}\mathbb{R}\varepsilon{\mathbb{C}}\mathcal{O}$, Institut
d’Astrophysique de Paris,
UMR 7095, CNRS, Sorbonne Université,
98bis boulevard Arago, 75014 Paris, France  Lavinia Heisenberg
<EMAIL_ADDRESS>Institut für Theoretische Physik, ETH Zürich,
Wolfgang-Pauli-Str. 27, 8093 Zürich, Switzerland  Institute for Theoretical
Physics, Heidelberg University, Philosophenweg 16, 69120 Heidelberg, Germany
François Larrouturou<EMAIL_ADDRESS>Deutsches Elektronen-
Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

###### Abstract

Motivated by a version of the “Dipolar Dark Matter” model, that aims at a
relativistic completion of the phenomenology of MOND, we investigate the
gravitational polarization mechanism in the canonical bimetric theory with an
effective matter coupling. We explicitly show the fundamental obstacle why
such theories cannot achieve a consistent gravitational polarization, and thus
fail to recover the MONDian phenomenology at low energies.

###### pacs:

95.35.+d, 04.50.Kd

††preprint: DESY-23-010

##  I Introduction

Despite passing all current tests with flying colors (the last in date being
the discovery of gravitational waves [1, 2]), our theory of General Relativity
(GR) suffers from some “darkness”, _i.e._ , yet ambiguous gravitational
phenomena that can be explained either by invoking non-gravitational
ingredients, or by “beyond-GR” completions [3]. The paramount instance is that
stars orbiting in periphery of galaxies do not move accordingly to GR
predictions [4, 5], but instead behave as if the visible mass of the galaxy
was only a small fraction of the total mass. The usual way to address this
“missing mass” problem is to postulate the existence of some unknown non-
relativistic form of matter, insensitive to electromagnetic forces (otherwise
it would have already been seen), coined “cold dark matter” [6]. If this
paradigm is a cornerstone of current physics, there is still no evidence of
this new particle, despite decades of experimental efforts [7, 8, 9].

Another way to explain the missing mass problem could be to consider that in
some regime, the laws of gravity deviate from GR. This path has been initiated
by the 40-years old seminal work of Milgrom [10], setting a very simple
modification of Newton’s law that accounts for the motion of stars in galaxies
and provides a natural explanation for the Tully-Fischer relation.
Nevertheless, the MOND (for MOdified Newtonian Dynamics) approach is yet only
phenomenological, despite great efforts to embed it within a viable
fundamental theory (see _e.g._ [11, 12, 13], and [14] for a comprehensive
review on the MONDian approach). An alternative possibility to build such a
relativistic completion is to exploit the dielectric analogy of MOND [15],
which has led to several “Dipolar Dark Matter” (DDM) models [16, 17, 18, 19,
20, 21]. In the present work, we focus on a subclass of such models [19, 20,
21] that are based on a bimetric extension of GR in which the two metrics are
supposed to host opposite “gravitational charges” (see also [22] for a
different type of MONDian bimetric theory).

Thus the DDM phenomenology relies on a mechanism of gravitational polarization
associated with opposite gravitational charges. It was indeed claimed in [20],
that in the non-relativistic limit, using the canonical bimetric theory, the
matter present in each sector feel forces with opposite signs, leading to the
desired polarization. Nevertheless, as shown in the present work, due to a
missing term in [20], unfortunately such polarisation mechanism cannot happen
in this bimetric framework, which ruins this DDM model. However previously
proposed DDM models [16, 17] in standard GR are still viable.

The plan of this paper is as follows. The DDM model used, which is an
extension of the one studied in [20, 21], is presented in Sec. II. The proof
of the impossibility of a gravitational polarization mechanism in this general
set up is then exposed in Sec. III. Concluding remarks are presented in Sec.
IV.

##  II Model

This work investigates an attempt for a Dipolar Dark Matter (DDM) model,
generalizing the one used in previous works (see _e.g._ [21]). Starting from
the two dynamical metrics $\\{g_{\mu\nu},f_{\mu\nu}\\}$,111Through this work,
we use a mostly positive signature for both metrics, denote space-time indices
with greek letters and, once a proper 3+1 decomposition achieved, purely
spatial indices with latin letters. For simplicity, we set $G=c=\hbar=1$,
unless otherwise specified. our action reads

$\mathcal{S}=\int\\!\\!\mathrm{d}^{4}x\,\mathcal{L}_{\text{grav}}+\int\\!\\!\mathrm{d}^{4}x\,\mathcal{L}_{\text{mat}}+\int\\!\\!\mathrm{d}^{4}x\,\mathcal{L}_{\text{bar}}\,,$
(1)

where $\mathcal{L}_{\text{mat}}$ describes the dark matter, coupled to both
$g_{\mu\nu}$ and $f_{\mu\nu}$, and where $\mathcal{L}_{\text{bar}}$ is the
usual baryonic matter, coupled to the metric $g_{\mu\nu}$. The gravitational
sector is given by

$\mathcal{L}_{\text{grav}}=\frac{M_{g}^{2}}{2}\sqrt{-g}\,R\left[g\right]+\frac{M_{f}^{2}}{2}\sqrt{-f}\,R\left[f\right]+M_{\text{eff}}^{2}m^{2}\,\mathcal{V}\left[g,f\right]\,,$
(2)

where $R$ denotes the usual Ricci scalar, $M_{g}$ and $M_{f}$ are the Planck
scales associated with each metrics and $M_{\text{eff}}^{2}\,m^{2}$ is an
unspecified scale. As for the interaction term $\mathcal{V}\left[g,f\right]$,
we will naturally use the ghost-free potential function of de Rham, Gabadadze
and Tolley [23, 24], parametrized by five arbitrary constants
$\\{\alpha_{n}\\}$ as

$\mathcal{V}\left[g,f\right]=\sqrt{-g}\sum_{n=0}^{4}\alpha_{n}\,e_{n}\left(X\right)\,,$
(3)

where the square-root matrix $X^{\mu}_{\ \nu}$ is defined as $X^{\mu}_{\
\rho}X^{\rho}_{\ \nu}=g^{\mu\rho}f_{\rho\nu}$, together with its inverse
$Y^{\mu}_{\ \nu}$, such that $Y^{\mu}_{\ \rho}Y^{\rho}_{\
\nu}=f^{\mu\rho}g_{\rho\nu}$, and $e_{n}(X)$ represent the elementary
symmetric polynomials

$\displaystyle e_{0}(X)$ $\displaystyle=1\,,$ (4) $\displaystyle e_{1}(X)$
$\displaystyle=\bigl{[}X\bigr{]}\,,$ $\displaystyle e_{2}(X)$
$\displaystyle=\frac{1}{2}\bigl{(}\bigl{[}X\bigr{]}^{2}-\bigl{[}X^{2}\bigr{]}\bigr{)}\,,$
$\displaystyle e_{3}(X)$
$\displaystyle=\frac{1}{6}\bigl{(}\bigl{[}X\bigr{]}^{3}-3\bigl{[}X\bigr{]}\bigl{[}X^{2}\bigr{]}+2\bigl{[}X^{3}\bigr{]}\bigr{)}\,,$
$\displaystyle e_{4}(X)$
$\displaystyle=\frac{1}{24}\bigl{(}\bigl{[}X\bigr{]}^{4}-6\bigl{[}X\bigr{]}^{2}\bigl{[}X^{2}\bigr{]}+3\bigl{[}X^{2}\bigr{]}^{2}+8\bigl{[}X\bigr{]}\bigl{[}X^{3}\bigr{]}-6\bigl{[}X^{4}\bigr{]}\bigr{)}\,,$

where the brackets denote the trace operation. The interaction term (3) is a
generalization of the one taken in previous works, corresponding to the choice
$\alpha_{n}=\alpha^{4-n}\beta^{n}$, where $\alpha$ and $\beta$ enter the
effective metric (6), see _e.g._ Eq. (2.8) of [20]. The Lagrangian (2) with
(3) is nothing but the gravitational sector of the canonical ghost-free
bimetric theory of Hassan and Rosen [25, 26], where $\alpha_{0}$ (resp.
$\alpha_{4}$) plays the role of a bare cosmological constant for the metric
$g_{\mu\nu}$ (resp. $f_{\mu\nu}$). This gravitational sector is symmetric
under the exchange of the metrics $g_{\mu\nu}\leftrightarrow f_{\mu\nu}$
together with the relabeling
$\\{M_{g},\alpha_{n}\\}\leftrightarrow\\{M_{f},\alpha_{4-n}\\}$.

As for the matter sector of our model, in addition to the baryonic matter
$\mathcal{L}_{\text{bar}}$ coupled to $g_{\mu\nu}$ in the usual way, we will
consider the following Lagrangian

$\mathcal{L}_{\text{mat}}=\mathcal{L}^{g}_{\text{mat}}[J_{g}^{\mu};g_{\mu\nu}]+\mathcal{L}^{f}_{\text{mat}}[J_{f}^{\mu};f_{\mu\nu}]+\mathcal{L}^{\text{eff}}_{\text{mat}}[A_{\mu};g^{\text{eff}}_{\mu\nu}]+\mathcal{L}^{\text{int}}_{\text{mat}}[J_{*g}^{\mu},J_{*f}^{\mu},A_{\mu}]\,,$
(5)

that describes three matter fields: (i) particles coupled to $g_{\mu\nu}$ and
described by a conserved222With the conserved current we can describe
particles without interaction or a perfect fluid with zero temperature. Adding
a dependence of the specific entropy we could also describe the case with
finite temperature, which we will not consider in this work. current
$J_{g}^{\mu}$ (such that $\nabla^{g}_{\mu}J_{g}^{\mu}=0$); (ii) particles
coupled to $f_{\mu\nu}$ and described by a conserved current $J_{f}^{\mu}$
(such that $\nabla^{f}_{\mu}J_{f}^{\mu}=0$); (iii) an Abelian $U(1)$ vector
field $A_{\mu}$ coupled to the effective metric [27]

$g^{\text{eff}}_{\mu\nu}=\alpha^{2}g_{\mu\nu}+2\alpha\beta\,g_{\mu\rho}X^{\rho}_{\
\nu}+\beta^{2}f_{\mu\nu}\,,$ (6)

where $\alpha$ and $\beta$ are arbitrary constants. The last term entering the
matter sector (5) is an interaction between the particles and the $U(1)$ field

$\mathcal{L}^{\text{int}}_{\text{mat}}=\bigl{(}J_{*g}^{\mu}-J_{*f}^{\mu}\bigr{)}A_{\mu}\,.$
(7)

where $J_{*g}^{\mu}\equiv\sqrt{-g}\,J_{g}^{\mu}$ and
$J_{*g}^{\mu}\equiv\sqrt{-f}\,J_{f}^{\mu}$ are metric independent matter
currents.

Note that the interaction term (7) induces a supplementary, although indirect,
coupling between the two metrics, and can thus reintroduce the ghost that was
suppressed by the interaction term (3) [27, 28]. Restricting the interaction
term to the case $\alpha_{n}=\alpha^{4-n}\beta^{n}$ as was done previously, it
has been shown that the ghost is not excited for energies beyond the
decoupling limit of the theory [21]. In our more general framework (3), it
could be that the theory is plagued by a ghost at energies well below the
decoupling limit. Nevertheless, as our aim is to prove that those theories are
not well suited for a MONDian phenomenology, and as such a ghost would
definitively not play a role in this low-energy analysis, we will consider the
larger class of theories with (7) and arbitrary coefficients
$\\{\alpha_{n}\\}$.

##  III Proof of the absence of polarization mechanism

Let us now turn to the study of the weak-field (non-relativistic) limit of the
theory (1), investigated along the lines of [20]. In order to reach a
consistent weak-field limit, we assume that the Universe is empty in average,
_i.e._ that the matter only enters at perturbative order in the action.
Moreover, seeking for a viable gravitational polarization mechanism requires
to perturb both metrics around a common background metric, say
$\bar{g}_{\mu\nu}$. Note that, as we will see later on, an empty background is
not necessarily a flat one, as the interaction term will play the role of an
effective cosmological constant. Thus $\bar{g}_{\mu\nu}$ is not the
Minkowskian metric, but rather a de Sitter (or anti-de Sitter) one. We then
perturb each metric as

$\displaystyle g_{\mu\nu}$
$\displaystyle=(\bar{g}_{\mu\nu}+h_{\mu\nu})^{2}\equiv(\bar{g}_{\mu\rho}+h_{\mu\rho})\bar{g}^{\rho\lambda}(\bar{g}_{\lambda\nu}+h_{\lambda\nu})\,,$
(8) $\displaystyle f_{\mu\nu}$
$\displaystyle=(\bar{g}_{\mu\nu}+\ell_{\mu\nu})^{2}\equiv(\bar{g}_{\mu\rho}+\ell_{\mu\rho})\bar{g}^{\rho\lambda}(\bar{g}_{\lambda\nu}+\ell_{\lambda\nu})\,.$

Up to quadratic order in the metric perturbations, the Lagrangian becomes333We
do not write here the zeroth order in perturbation, naturally irrelevant for
our study. From now on, the barred quantities are associated with the
background metric $\bar{g}_{\mu\nu}$, with which the indices are operated.

$\displaystyle\mathcal{L}^{(2)}=\sqrt{-\bar{g}}\,\biggl{\\{}$
$\displaystyle-\Bigl{(}M_{g}^{2}h^{\mu\nu}+M_{f}^{2}\ell^{\mu\nu}\Bigr{)}\overline{G}_{\mu\nu}-\Bigl{(}M_{g}^{2}\Lambda_{g}h^{\mu\nu}+M_{f}^{2}\Lambda_{f}\ell^{\mu\nu}\Bigr{)}\bar{g}_{\mu\nu}$
(9) $\displaystyle-
M_{g}^{2}h^{\mu\nu}\overline{\mathcal{E}}^{\rho\lambda}_{\mu\nu}h_{\rho\lambda}-M_{f}^{2}\ell^{\mu\nu}\overline{\mathcal{E}}^{\rho\lambda}_{\mu\nu}\ell_{\rho\lambda}+h_{\mu\nu}\Bigl{(}T_{\text{bar}}^{\mu\nu}+T_{g}^{\mu\nu}\Bigr{)}+\ell_{\mu\nu}\,T_{f}^{\mu\nu}$
$\displaystyle+\frac{M_{\text{eff}}^{2}\,m^{2}}{2}\Bigl{[}\Gamma_{0}h^{\mu\nu}h^{\rho\lambda}+2\Gamma_{1}h^{\mu\nu}\ell^{\rho\lambda}+\Gamma_{2}\ell^{\mu\nu}\ell^{\rho\lambda}\Bigr{]}\bigl{(}\bar{g}_{\mu\nu}\bar{g}_{\lambda\rho}-\bar{g}_{\mu\rho}\bar{g}_{\lambda\nu}\bigr{)}\biggr{\\}}\,,$

where $\overline{G}_{\mu\nu}$ is the Einstein tensor associated to
$\bar{g}_{\mu\nu}$, and $\overline{\mathcal{E}}^{\rho\sigma}_{\mu\nu}$ is the
Lichnerowicz operator

$\displaystyle-2\overline{\mathcal{E}}^{\rho\lambda}_{\mu\nu}h_{\rho\lambda}$
$\displaystyle=\overline{\Box}\bigl{(}h_{\mu\nu}-\bar{g}_{\mu\nu}h\bigr{)}+\overline{\nabla}_{\mu}\\!\overline{\nabla}_{\nu}h-2\overline{\nabla}_{(\mu}H_{\nu)}+\bar{g}_{\mu\nu}\overline{\nabla}_{\rho}H^{\rho}$
(10)
$\displaystyle-2\overline{C}_{\mu\rho\lambda\nu}h^{\rho\lambda}-\frac{2}{3}\,\overline{R}\,\Bigl{(}h_{\mu\nu}-\frac{1}{4}\bar{g}_{\mu\nu}h\Bigr{)}-\overline{G}_{\rho(\mu}h^{\rho}_{\nu)}\,.$

Here $\overline{C}_{\mu\nu\rho\lambda}$ is the Weyl tensor, we denoted the
trace as $h\equiv h^{\rho\sigma}\bar{g}_{\rho\sigma}$ and posed
$H_{\rho}\equiv\overline{\nabla}^{\sigma}h_{\rho\sigma}$. The stress-energy
tensors are defined as usual by

$T^{\mu\nu}_{\text{bar}}=\frac{2}{\sqrt{-g}}\frac{\delta\mathcal{L}_{\text{bar}}}{\delta
g_{\mu\nu}}\,,\qquad
T^{\mu\nu}_{g}=\frac{2}{\sqrt{-g}}\frac{\delta\mathcal{L}_{\text{mat}}}{\delta
g_{\mu\nu}}\,,\qquad
T^{\mu\nu}_{f}=\frac{2}{\sqrt{-f}}\frac{\delta\mathcal{L}_{\text{mat}}}{\delta
f_{\mu\nu}}\,.$ (11)

They are first-order quantities (as we required an empty Universe) that obey
conservation laws with respect to the background metric at leading order,
namely

$\overline{\nabla}_{\nu}\Bigl{(}T^{\mu\nu}_{\text{bar}}+T^{\mu\nu}_{g}\Bigr{)}=\mathcal{O}(2)\,,\qquad\overline{\nabla}_{\nu}T^{\mu\nu}_{f}=\mathcal{O}(2)\,.$
(12)

Here we neglect the stress-energy tensor associated with the internal vector
field $A_{\mu}$, which is coupled to the effective metric (6) and is second
order $\mathcal{O}(2)$ in perturbation, see (3.5) in [20]. Finally, we have
introduced the effective cosmological constants

$\Lambda_{g}=-\frac{M_{\text{eff}}^{2}\,m^{2}}{M_{g}^{2}}\big{(}\alpha_{0}+3\alpha_{1}+3\alpha_{2}+\alpha_{3}\big{)}\,,\qquad\Lambda_{f}=-\frac{M_{\text{eff}}^{2}\,m^{2}}{M_{f}^{2}}\big{(}\alpha_{1}+3\alpha_{2}+3\alpha_{3}+\alpha_{4}\big{)}\,,$
(13)

as well as the combinations

$\Gamma_{0}\equiv\alpha_{0}+2\alpha_{1}+\alpha_{2}\,,\qquad\Gamma_{1}\equiv\alpha_{1}+2\alpha_{2}+\alpha_{3}\,,\qquad\Gamma_{2}\equiv\alpha_{2}+2\alpha_{3}+\alpha_{4}\,.$
(14)

The Lichnerowicz operator considered in Ref. [20] had a missing term: Eq.
(3.10) there missed the last term in (10), _i.e._
$-\overline{G}_{\rho(\mu}h^{\rho}_{\nu)}$, proportional to the Einstein tensor
of the background. This term plays a major role when establishing the relation
(18) below, and in turn the Eq. (20), instead of the relation (3.13) found in
[20]. Unfortunately, this invalidates the model [20] for the proposed
application to MONDian phenomenology.

At the linear level, the variation of the perturbed Lagrangian (9) gives the
field equations

$\overline{G}_{\mu\nu}=-\Lambda_{g}\,\bar{g}_{\mu\nu}=-\Lambda_{f}\,\bar{g}_{\mu\nu}\,,$
(15)

which can be fulfilled only for a common cosmological constant in the
background, $\Lambda_{g}=\Lambda_{f}$. This equality is nothing but the
condition that the background metric $\bar{g}_{\mu\nu}$ is indeed the same for
both sectors. When translated in terms of the notation (14), this condition
gives a generalisation of Eq.(3.8) of [20]:

$\frac{\Gamma_{0}+\Gamma_{1}}{M_{g}^{2}}=\frac{\Gamma_{1}+\Gamma_{2}}{M_{f}^{2}}\,.$
(16)

Turning now to the quadratic level, the Einstein field equations for the two
metric perturbations $h_{\mu\nu}$ and $\ell_{\mu\nu}$ read

$\displaystyle
2M_{g}^{2}\,\overline{\mathcal{E}}^{\mu\nu}_{\rho\lambda}h^{\rho\lambda}+M_{\text{eff}}^{2}\,m^{2}\biggl{[}\Gamma_{0}\bigl{(}h^{\mu\nu}-h\,\bar{g}^{\mu\nu}\bigr{)}+\Gamma_{1}\bigl{(}\ell^{\mu\nu}-\ell\,\bar{g}^{\mu\nu}\bigr{)}\biggr{]}=T_{\text{bar}}^{\mu\nu}+T_{g}^{\mu\nu}\,,$
(17) $\displaystyle
2M_{f}^{2}\,\overline{\mathcal{E}}^{\mu\nu}_{\rho\lambda}\ell^{\rho\lambda}+M_{\text{eff}}^{2}\,m^{2}\biggl{[}\Gamma_{1}\bigl{(}h^{\mu\nu}-h\,\bar{g}^{\mu\nu}\bigr{)}+\Gamma_{2}\bigl{(}\ell^{\mu\nu}-\ell\,\bar{g}^{\mu\nu}\bigr{)}\biggr{]}=T_{f}^{\mu\nu}\,.$

As usual in bigravity, we observe that the mass eigenstates are different from
the kinetic ones. Using the expression (10) as well as the background
condition (15) which says that $\Lambda\equiv\Lambda_{g}=\Lambda_{f}$, the
divergence of the Lichnerowicz operator reads

$2\overline{\nabla}_{\nu}\left(\overline{\mathcal{E}}^{\mu\nu}_{\rho\lambda}h^{\rho\lambda}\right)=2\overline{G}_{\rho\lambda}\overline{\nabla}^{[\mu}h^{\rho]\lambda}+\overline{\nabla}_{\rho}\Bigl{(}\overline{G}^{\lambda[\mu}h^{\rho]}_{\lambda}\Bigr{)}=\Lambda\,\overline{\nabla}_{\rho}\bigl{[}h^{\rho\mu}-h\,\bar{g}^{\rho\mu}\bigr{]}\,,$
(18)

together with the mirror relation obtained by replacing
$h^{\rho\sigma}\rightarrow\ell^{\rho\sigma}$. Therefore, using the
conservation laws (12) of the stress-energy tensors, we find that the
divergence of both field equations (17) boil down to the same constraint,
independent of the background relation (16),

$\Gamma_{1}\overline{\nabla}_{\rho}\big{[}h^{\rho\mu}-\ell^{\rho\mu}-h\,\bar{g}^{\rho\mu}+\ell\,\bar{g}^{\rho\mu}\big{]}=0\,.$
(19)

Imposing $\Gamma_{1}=0$, the two sectors would decorrelate at first order, as
clear from Eq. (17), and thus no gravitational polarization mechanism could be
possible, by definition. Therefore, we are forced to impose

$\overline{\nabla}_{\rho}\big{[}h^{\rho\mu}-\ell^{\rho\mu}\big{]}=\overline{\nabla}^{\mu}\big{[}h-\ell\big{]}\,.$
(20)

We now investigate the phenomenological consequences of this constraint in the
non-relativistic limit $c\to\infty$. In this regime, we consider Newtonian-
like systems, with sizes significantly smaller that the effects of a
reasonable cosmological constant. We can thus discard the effect of the
effective cosmological constant $\Lambda$, and take a Minkowskian ansatz for
$\bar{g}_{\mu\nu}$, which allows a proper 3+1 decomposition of the components.
On this flat background, we parametrize the metric perturbations by single
Newtonian potentials $U_{g}$ and $U_{f}$ and two “PPN” constants $\gamma_{g}$
and $\gamma_{f}$ as

$\displaystyle
h^{00}=\frac{U_{g}}{c^{2}}+\mathcal{O}\left(\frac{1}{c^{4}}\right)\,,\qquad
h^{0i}=\mathcal{O}\left(\frac{1}{c^{3}}\right)\,,\qquad
h^{ij}=\gamma_{g}\,\frac{U_{g}}{c^{2}}\,\delta_{ij}+\mathcal{O}\left(\frac{1}{c^{4}}\right)\,,$
(21)
$\displaystyle\ell^{00}=\frac{U_{f}}{c^{2}}+\mathcal{O}\left(\frac{1}{c^{4}}\right)\,,\qquad\ell^{0i}=\mathcal{O}\left(\frac{1}{c^{3}}\right)\,,\qquad\ell^{ij}=\gamma_{f}\,\frac{U_{f}}{c^{2}}\,\delta_{ij}+\mathcal{O}\left(\frac{1}{c^{4}}\right)\,.$

The $ij$ components of the field equations (17) set $\gamma_{g}=\gamma_{f}=1$,
and the constraint equation (20) then gives

$\partial_{i}\big{[}U_{g}-U_{f}\big{]}=0\,.$ (22)

However, the gravitational polarization mechanism discussed in [20] required
that the gradients of $U_{g}$ and $U_{f}$ have opposite signs: as the
particles described by $J_{g}^{\mu}$ and $J_{f}^{\mu}$ feel the force
respectively due to the potentials $U_{g}$ and $U_{f}$, they will be displaced
in the opposite direction. The $U(1)$ field $A_{\mu}$ coupling both sectors
would then create an effective polarization of the spacetime medium,
reproducing the main feature of the model [15]. It is thus clear that the
constraint (22) is incompatible with such a mechanism.

##  IV Conclusion

This study has proven that, contrarily to what was claimed in the work [20],
the framework based on canonical bimetric theory is unable to produce a
consistent gravitational polarization mechanism. This shortage unfortunately
spoils its capacity to account for an effective MONDian phenomenology, and
thus, to be a proper relativistic completion of the MOND paradigm. Note that
another DDM model, based on standard GR [16, 17], is still viable in this
respect; this model has recently been simulated numerically on galactic scales
[29]. Nevertheless, there is still hope that a viable mechanism could be
achieved in a ghost-free bigravity framework, as DDM was based on the Hassan-
Rosen model, but other bimetric theories exist, for instance the newly
proposed minimal bigravity [30], which has the advantage of having stable
cosmological solutions as well as escaping the Higuchi ghost present in usual
bigravity scenarii.

###### Acknowledgements.

F.L. would like to thank the Institut d’Astrophysique de Paris for hosting him
during both early and late stages of this project.

## References

  * Abbott  _et al._ [2016] B. Abbott _et al._ (LIGO Scientific Collaboration and VIRGO Collaboration), Observation of gravitational waves from a binary black hole merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2017] B. Abbott _et al._ (LIGO Scientific Collaboration and VIRGO Collaboration), Gw170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Heisenberg [2019] L. Heisenberg, A systematic approach to generalisations of General Relativity and their cosmological implications, Phys. Rept. 796, 1 (2019), arXiv:1807.01725 [gr-qc] .
  * Bosma [1981] A. Bosma, 21-cm line studies of spiral galaxies. ii. the distribution and kinematics of neutral hydrogen in spiral galaxies of various morphological types, Astron. J. 86, 1825 (1981).
  * Rubin  _et al._ [1982] V. Rubin, J. Ford, N. Thonnard, and D. Burstein, Rotational properties of 23 sb galaxies, Astrophys. J. 261, 439 (1982).
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Particle dark matter: Evidence, candidates and constraints, Phys. Rept. 405, 279 (2005), hep-ph/0404175 .
  * Aprile  _et al._ [2020] E. Aprile _et al._ (XENON), Excess electronic recoil events in XENON1T, Phys. Rev. D 102, 072004 (2020), arXiv:2006.09721 [hep-ex] .
  * Albert  _et al._ [2020] A. Albert _et al._ (ANTARES, IceCube), Combined search for neutrinos from dark matter self-annihilation in the Galactic Center with ANTARES and IceCube, Phys. Rev. D 102, 082002 (2020), arXiv:2003.06614 [astro-ph.HE] .
  * Aad  _et al._ [2021] G. Aad _et al._ (ATLAS), Search for new phenomena in events with two opposite-charge leptons, jets and missing transverse momentum in pp collisions at $\sqrt{\mathrm{s}}$ = 13 TeV with the ATLAS detector, JHEP 04, 165, arXiv:2102.01444 [hep-ex] .
  * Milgrom [1983] M. Milgrom, A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis, Astrophys. J. 270, 365 (1983).
  * Bekenstein [2004] J. Bekenstein, Relativistic gravitation theory for the modified Newtonian dynamics paradigm, Phys. Rev. D 70, 083509 (2004), astro-ph/0403694 .
  * Bruneton and Esposito-Farèse [2007] J.-P. Bruneton and G. Esposito-Farèse, Field-theoretical formulations of MOND-like gravity, Phys. Rev. D 76, 124012 (2007), arXiv:0705.4043 [gr-qc] .
  * Deffayet  _et al._ [2014] C. Deffayet, G. Esposito-Farese, and R. P. Woodard, Field equations and cosmology for a class of nonlocal metric models of MOND, Phys. Rev. D 90, 064038 (2014), [Addendum: Phys.Rev.D 90, 089901 (2014)], arXiv:1405.0393 [astro-ph.CO] .
  * Famaey and McGaugh [2012] B. Famaey and S. McGaugh, Modified Newtonian dynamics (MOND): Observational phenomenology and relativistic extensions, Living Rev. Rel. 15, 10 (2012), arXiv:1112.3960 [astro-ph.CO] .
  * Blanchet [2007] L. Blanchet, Gravitational polarization and the phenomenology of MOND, Class. Quant. Grav. 24, 3529 (2007), astro-ph/0605637 .
  * Blanchet and Le Tiec [2008] L. Blanchet and A. Le Tiec, Model of dark matter and dark energy based on gravitational polarization, Phys. Rev. D 78, 024031 (2008), astro-ph/0804.3518 .
  * Blanchet and Le Tiec [2009] L. Blanchet and A. Le Tiec, Dipolar dark matter and dark energy, Phys. Rev. D 80, 023524 (2009), arXiv:0901.3114 [astro-ph] .
  * Bernard and Blanchet [2015] L. Bernard and L. Blanchet, Phenomenology of dark matter via a bimetric extension of general relativity, Phys. Rev. D 91, 103536 (2015), arXiv:1410.7708 [astro-ph] .
  * Blanchet and Heisenberg [2015a] L. Blanchet and L. Heisenberg, Dark Matter via Massive (bi-)Gravity, Phys. Rev. D 91, 103518 (2015a), arXiv:1504.00870 [gr-qc] .
  * Blanchet and Heisenberg [2015b] L. Blanchet and L. Heisenberg, Dipolar Dark Matter with Massive Bigravity, J. Cosmo. Astropart. 12, 026 (2015b), arXiv:1505.05146 [gr-qc] .
  * Blanchet and Heisenberg [2017] L. Blanchet and L. Heisenberg, Dipolar Dark Matter as an Effective Field Theory, Phys. Rev. D 96, 083512 (2017), arXiv:1701.07747 [gr-qc] .
  * Milgrom [2009] M. Milgrom, Bimetric MOND gravity, Phys. Rev. D 80, 123536 (2009), arXiv:0912.0790 [gr-qc] .
  * de Rham  _et al._ [2011] C. de Rham, G. Gabadadze, and A. Tolley, Resummation of Massive Gravity, Phys. Rev. Lett. 106, 231101 (2011), arXiv:1011.1232 [hep-th] .
  * de Rham [2014] C. de Rham, Massive Gravity, Living Rev.Rel. 17, 7 (2014), arXiv:1401.4173 [hep-th] .
  * Hassan and Rosen [2012a] S. Hassan and R. Rosen, Bimetric gravity from ghost-free massive gravity, J. High Energy Phys. 1202, 126, arXiv:1109.3515 [hep-th] .
  * Hassan and Rosen [2012b] S. Hassan and R. Rosen, Confirmation of the secondary constraint and absence of ghost in massive gravity and bimetric gravity, J. High Energy Phys. 1204, 123, arXiv:1111.2070 [hep-th] .
  * de Rham  _et al._ [2015] C. de Rham, L. Heisenberg, and R. Ribeiro, On coupling to matter in massive (bi-)gravity, Class. Quantum Grav. 32, 035022 (2015), arXiv:1408.1678 [hep-th] .
  * de Rham  _et al._ [2014] C. de Rham, L. Heisenberg, and R. Ribeiro, Ghosts & matter couplings in massive (bi-& multi-)gravity, Phys. Rev. D 90, 124042 (2014), arXiv:1409.3834 [hep-th] .
  * Stahl  _et al._ [2022] C. Stahl, B. Famaey, G. Thomas, Y. Dubois, and R. Ibata, Dipolar dark matter simulations on galaxy scales with the RAMSES code, Mon. Not. Roy. Astron. Soc. 517, 498 (2022), arXiv:2209.07831 [astro-ph.CO] .
  * De Felice  _et al._ [2021] A. De Felice, F. Larrouturou, S. Mukohyama, and M. Oliosi, Minimal Theory of Bigravity: construction and cosmology, JCAP 04, 015, arXiv:2012.01073 [gr-qc] .