11institutetext: Osaka University, Japan
22institutetext: University of Alberta, Canada
33institutetext: Nara Institute of Science and Technology
33email<EMAIL_ADDRESS>
# Automated Cancer Subtyping via Vector Quantization Mutual Information
Maximization

Zheng Chen indicates joint first authors.11 Lingwei Zhu††footnotemark:  22
Ziwei Yang 33 Takashi Matsubara 11

###### Abstract

Cancer subtyping is crucial for understanding the nature of tumors and
providing suitable therapy. However, existing labelling methods are medically
controversial, and have driven the process of subtyping away from teaching
signals. Moreover, cancer genetic expression profiles are high-dimensional,
scarce, and have complicated dependence, thereby posing a serious challenge to
existing subtyping models for outputting sensible clustering. In this study,
we propose a novel clustering method for exploiting genetic expression
profiles and distinguishing subtypes in an unsupervised manner. The proposed
method adaptively learns categorical correspondence from latent
representations of expression profiles to the subtypes output by the model. By
maximizing the problem-agnostic mutual information between input expression
profiles and output subtypes, our method can automatically decide a suitable
number of subtypes. Through experiments, we demonstrate that our proposed
method can refine existing controversial labels, and, by further medical
analysis, this refinement is proven to have a high correlation with cancer
survival rates.

###### Keywords:

Cancer Subtypes Information Maximization Clustering

##  1 Introduction

Cancer is by far one of the deadliest epidemiological diseases known to
humans: consider the breast cancer which is the most prevalent (incidence
47.8% worldwide) and the most well-studied cancer in the world [32], the
5-year mortality rate can still reach 13.6% [1]. Its heterogeneity is
considered as the crux of limiting the efficacy of targeted therapies and
compromising treatment outcomes since some tumors that differ radically at the
molecular level might exhibit highly resemblant morphological appearance [22].
Increasing evidence from modern transcriptomic studies has supported the
assumption that each specific cancer is composed of multiple categories (known
as cancer subtypes) [33, 4]. Reliably identifying cancer subtypes can
significantly facilitate the prognosis and personalized treatment [21].
However, currently there is a fierce debate in the cancer community: given
transcriptomic data of one cancer, authoritative resources put that there
might be different number of subtypes from distinct viewpoints, that is, the
fiducial definition of the subtypes is constantly undergoing calibration [12],
suggesting for the majority of cancers the ground-truth labeling remains
partially unavailable and awaits better definition.

In the data science community, the lack of ground truth for the cancer data
can be addressed as a clustering problem [11], in which the clusters give a
hint on the underlying subtypes. Such clustering methods rely crucially on the
quality of the data and suitable representations. Modern subtyping methods
typically leverage molecular transcriptomic expression profiles (expression
profiles in short) which consist of genetic and microRNA (miRNA) expressions
that characterize the cancer properties [26, 21]. However, several dilemmas
exist in the way of fully exploiting the power of expression profiles:

  * •

_High-dimensionality_ : the expression profiles are typically of $>60,000$
dimensions; even after typical preprocessing the dimension can still be
$>10,000$.

  * •

_Scarcity_ : cancer data are scarce and costly. Even for the most well-studied
breast cancer, the largest public available dataset consists of expression
profiles from around only 1500 subjects [30];

  * •

_Dependence_ : expression profiles have complicated dependence: a specific
expression might be under joint control of several genes, and sometimes such
the joint regulation can be circular, forming the well-known gene regulation
network [10].

To extract information from the inherently high-dimensional expression
profiles for tractable grouping [9], traditional methods preprocess the data
via variants of principal components analysis (PCA) or least absolute
shrinkage and selection operator (LASSO) [3] for reducing the dimensionality
of the data. However, expression profiles with such complicated dependence
have already been shown to not perform well with PCA and LASSO [14], since
many seemingly less salient features can play an important role in the gene
regulation network. Motivated by the resurgence of deep learning techniques,
recently the community has seen promising applications leveraging deep
autoencoders (AEs) or variational AEs (VAEs) for compressing the data into a
lower-dimensional latent space that models the underlying genetic regulation
[33]. However, VAEs with powerful autoregressive decoders often ignore the
latent spaces [25, 8], which runs the risk of overfitting [28]. Furthermore,
the latent representation is assumed to be continuous variables (usually
Gaussian) [18, 31], which is at odds with the inherently categorical cancer
subtypes [5]. As a result, those subtyping models might have poor performance
as well as generalization ability.

Aside from feature extraction, another issue concerns the grouping process
itself. Given extracted features from the expression profiles, the above-
mentioned methods usually apply similarity-based clustering algorithms such as
K-means for subsequent grouping. However, such methods require strong
assumptions on the data and are sensitive to representations [27]: one will
have to define a similarity metric for the data (often Euclidean) and find
appropriate transformations (such as logarithm transform) as informative
features. Unsuitable choices of the metric and transformation can greatly
degrade the model performance. Recently, mutual information has been gaining
huge popularity in deep representation learning as a replacement for
similarity metrics [13, 6]: it is the unique measure of relatedness between a
pair of variables invariant to invertible transformations of the data, hence
one does not need to find a _right_ representation [20]. Better yet, if two
genes share more than one bit of information, then the underlying mechanism
must be more subtle than just _on_ and _off_. Such subtlety and more general
dependence can be captured by the mutual information [27].

In this paper, we propose a novel, generally applicable clustering method that
is capable of fully exploiting the expression profiles and outputting sensible
cancer subtyping solutions. Besides tackling the above-mentioned problems in a
unified and consistent manner, the proposed method has an intriguing property
of automatically adjusting the number of groups thanks to its special
architecture, which stands as a sheer contrast to prior methods that
predetermine the number of groups by domain knowledge. Before introducing the
proposed architecture in Section 3, we summarize our contributions as follows:

  * •

(Algorithmic) Inspired by recent work, we propose a novel clustering method
vector quantization regularized information maximization (VQ-RIM) for cancer
subtyping. VQ-RIM maximizes mutual information in the categorical VQ-VAE
model, which results in a combination of VAE reconstruction loss and mutual
information loss. (Section 3.2)

  * •

(Effective) We compare the clustering results of VQ-RIM against existing
ground truth labels (together with controversial labels from the entirety of
labels) on different cancer datasets and find that VQ-RIM concords well with
the ground truth, which verifies the correctness of VQ-RIM. (Section 4.1 and
4.2.2)

  * •

(Medical) Extensive experiments on distinct cancers verify that VQ-RIM
produces subtyping that consistently outperform the controversial labels in
terms of enlarged separation of between-group life expectancies. The clearer
separation suggests VQ-RIM is capable of better capturing the underlying
characteristics of subtypes than controversial labels. We believe such results
are far-reaching in providing new insights into the unsettled debate on cancer
subtyping.(Section 4.2.2)

##  2 Related Work

Feature Extraction for Subtyping.  Building a model suitable for cancer
subtyping is non-trivial as a result of the cancer data scarcity. High
dimensionality and data scarcity pose a great challenge to automated models
for generating reliable clustering results [31]. Conventionally, the problem
is tackled by leveraging classic dimension reduction methods such as PCA [3].
However, since the progress of cancers is regulated by massive genes in a
complicated manner (which themselves are under the control of miRNAs), brute-
force dimension reduction might run the risk of removing informative features
[15]. On the other hand, recently popular AE-based models [33, 21], especially
VAEs, construct the feature space by reconstructing the input through a multi-
dimensional Gaussian posterior distribution in the latent space [31]. The
latent posterior learns to model the underlying causalities, which in the
cancer subtyping context corresponds to modeling the relationship among
expression profiles such as regulation or co-expression [33]. Unfortunately,
recent investigation has revealed that VAEs with powerful autoregressive
decoders easily ignore the latent space. As a result, the posterior could be
either too simple to capture the causalities; or too complicated so the
posterior distribution becomes brittle and at the risk of posterior collapse
[2, 25]. Moreover, the Gaussian posterior is at odds with the inherently
categorical cancer subtypes [5].

In this paper, we propose to leverage the categorical VQ-VAE to address the
aforementioned issues: (i) VQ-VAE does not train its decoder, preventing the
model from ignoring its latent feature space resulting from an over-powerful
decoder; (ii) VQ-VAE learns categorical correspondence between input
expression profiles, latent representations, and output subtypes, which
theoretically suggests better capability of learning more useful features.
(iii) the categorical latent allows the proposed model to automatically set a
suitable number of groups by plugging in mutual information maximization
classifier, which is not available for the VAEs.
Information Maximization for Subtyping.  Cancer subtyping is risk-sensitive
since misspecification might incur an unsuitable treatment modality. It is
hence desired that the clustering should be _as certain as possible for
individual prediction, while keeping subtypes as separated as possible_ [7,
11]. Further, to allow for subsequent analysis and further investigation of
medical experts, it is desired that the method should output probabilistic
prediction for each subject. In short, we might summarize the requirements for
the subtyping decision boundaries as (i) should not be overly complicated;
(ii) should not be located at where subjects are densely populated; (iii)
should output probabilistic predictions. These requirements can be formalized
via the information-theoretic objective as maximizing the mutual information
between the input expression profiles and the output subtypes [19, 29]. Such
objective is problem-agnostic, transformation-invariant, and unique for
measuring the relationship between pairs of variables. Superior performance
over knowledge-based heuristics has been shown by exploiting such an objective
[27].

##  3 Method

###  3.1 Problem Setting

Let $\mathcal{X}$ be a dataset
$\mathcal{X}=\\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N}\\}$, where
$\boldsymbol{x}_{i}\in\mathbb{R}^{d},1\leq i\leq N$ are $d$-dimensional
vectors consisting of cancer expression profiles. For a given
$\boldsymbol{x}$, our goal lies in determining a suitable cancer subtype
$y\in\\{1,2,\dots,K\\}$ given $\boldsymbol{x}$, where $K$ is not fixed
beforehand and needs to be automatically determined. Numeric values such as
$y=1,\dots,K$ do not bear any medical interpretation on their own and simply
represent distinct representations due to the underlying data. It is worth
noting while a label set $\mathcal{Y}$ is available, it comprises a small
subset of ground-truth labels $\mathcal{Y}_{gt}:=\\{y_{gt}\\}$ that have been
medically validated and a larger portion of _controversial_ labels
$\mathcal{Y}_{c}:=\\{y_{c}\\}$, with
$\mathcal{Y}_{gt}=\mathcal{Y}\textbackslash\mathcal{Y}_{c}$. Our approach is
to compare the clustering result $y$ of the proposed method against ground
truth labels $y_{gt}$ to see if they agree well, as a first step of
validation. We then compare $y$ against controversial labels $y_{c}$ and
conduct extensive experiments to verify that the proposed method achieves
improvement upon the subtyping given by $y_{c}$. Our goal is to unsupervisedly
learn a discriminative classifier $D$ which outputs conditional probability
$P(y|\boldsymbol{x},D)$. Naturally it is expected that
$\sum_{k=1}^{K}P(y=k|\boldsymbol{x},D)=1$ and we would like $D$ to be
probabilistic so the uncertainty associated with assigning data items can be
quantitized. Following [28], we assume the marginal class distribution
$P(y|D)$ is close to the prior $P(y)$ for all $k$. However, unlike prior work
[19, 28] we do not assume the amount of examples per class in $\mathcal{X}$ is
uniformly distributed due to the imbalance of subtypes in the data.

###  3.2 Proposed Model

####  3.2.1 Information Maximization.

Given expression profiles of subject $\boldsymbol{x}$, the discriminator
outputs a $K$-dimensional probability logit vector
$D(\boldsymbol{x})\in\mathbb{R}^{K}$. The probability of $\boldsymbol{x}$
belonging to any of the $K$ subtypes is given by the softmax parametrization:

$\displaystyle
P(y=k|\boldsymbol{x},D)=\frac{e^{D_{k}(\boldsymbol{x})}}{\sum_{k=1}^{K}e^{D_{k}(\boldsymbol{x})}},$

where $D_{k}(\boldsymbol{x})$ denotes the $k$-th entry of the vector
$D(\boldsymbol{x})$. Let us drop the dependence on $D$ for uncluttered
notation. It is naturally desired that the individual prediction be as certain
as possible, while the distance between the predicted subtypes as large as
possible. This consideration can be effectively reflected by the mutual
information between the input expression profiles and the output prediction
label. Essentially, the mutual information can be decomposed into the
following two terms:

$\displaystyle\begin{split}\hat{I}(\boldsymbol{x},y)\\!:=\\!\underbrace{-\sum_{k=1}^{K}\\!P(y=k)\log\\!P(y=k)}_{\widehat{\mathcal{H}}\left(P(y)\right)}+\,\alpha\underbrace{\sum_{i=1}^{N}\frac{1}{N}\sum_{k=1}^{K}P(y=k|\boldsymbol{x}_{i})\log
P(y=k|\boldsymbol{x}_{i})}_{-\widehat{\mathcal{H}}\left(P(y|\mathcal{X})\right)}.\end{split}$
(1)

which are the marginal entropy of labels
$\widehat{\mathcal{H}}\left(P(y=k)\right)$ and the conditional entropy
$\widehat{\mathcal{H}}\left(P(y|\mathcal{X})\right)$ approximated by $N$ Monte
Carlo samples $\boldsymbol{x}_{i},i\in\\{1,\dots,N\\}$. $\alpha$ is an
adjustable parameter for weighting the contribution, setting $\alpha=1$
recovers the standard mutual information formulation [19]. This formulation
constitutes the _regularized information maximization (RIM)_ part of the
proposed method. The regularization effect can be seen from the following:

  * $\boldsymbol{\cdot}$

Conditional entropy $\widehat{\mathcal{H}}\left(P(y|\mathcal{X})\right)$
encourages confident prediction by minimizing uncertainty. It effectively
captures the modeling principles that decision boundaries should not be
located at dense population of data [11].

  * $\boldsymbol{\cdot}$

Marginal entropy $\widehat{\mathcal{H}}\left(P(y)\right)$ aims to separate the
subtypes as far as possible. Intuitively, it attempts to keep the subtypes
_uniform_. Maximizing only
$\widehat{\mathcal{H}}\left(P(y|\mathcal{X})\right)$ tends to produce
degenerate solutions by removing subtypes [6, 19], hence
$\widehat{\mathcal{H}}\left(P(y)\right)$ serves as an effective regularization
for ensuring nontrivial solutions.

####  3.2.2 Categorical Latents Generative Feature Extraction.

Recent studies have revealed that performing RIM alone is often insufficient
for obtaining stable and sensible clustering solutions [6, 20, 28]:
Discriminative methods are prone to overfitting spurious correlations in the
data, e.g., some entry $A$ in the expression profiles might appear to have
direct control over certain other entries $B$. The model might naïvely
conclude that the appearance of $B$ shows positive evidence of $A$. However,
such relationship is in general not true due to existence of complicated
biological functional passways: Such pathways have complex (sometimes
circular) dependence between $A$ and $B$ [24]. Since discriminative methods
model $P(y|\boldsymbol{x})$ but not the data generation mechanism
$P(\boldsymbol{x})$ (and the joint distribution $P(\boldsymbol{x},y)$) [11],
such dependence between genes and miRNAs might not be effectively captured by
solely exploiting the discriminator, especially given the issues of data
scarcity and high dimensionality.

A generative model that explicitly captures the characteristics in
$P(\boldsymbol{x})$ is often introduced as a rescue for leveraging RIM-based
methods [13, 23, 28]. Such methods highlight the use of VAEs for modeling the
latent feature spaces underlying input $\mathcal{X}$: given input
$\boldsymbol{x}$, VAEs attempt to compress it to a lower-dimensional latent
$\boldsymbol{z}$, and reconstruct $\tilde{\boldsymbol{x}}$ from
$\boldsymbol{z}$. Recently there has been active research on leveraging VAEs
for performing cancer subtyping [31, 33]. However, existing literature
leverage continuous latents (often Gaussian) for tractability, which is at
odds with the inherently categorical cancer subtypes. Furthermore, VAEs often
ignore the latents which implies the extracted feature space is essential
dismissed and again runs the risk of overfitting [2].

We exploit the recent vector quantization variational auto-encoder (VQ-VAE)
[25] as the generative part of the proposed architecture. The categorical
latents of VQ-VAE are not only suitable for modeling inherently categorical
cancer subtypes, but also avoids the above-mentioned latent ignoring problem
[18]. In VQ-VAE, the latent embedding space is defined as
$\\{e_{i}\\}\in\mathbb{R}^{M\times l}$, where $M$ denotes the number of
embedding vectors and hence a $M$-way categorical distribution. $l<d$ is the
dimension of each latent embedding vector $e_{i},i\in\\{1,\dots,M\\}$. VQ-VAE
maps input $\boldsymbol{x}$ to a latent variable $\boldsymbol{z}$ via its
encoder $z_{e}(\boldsymbol{x})$ by performing a nearest neighbor search among
the embedding vectors $e_{i}$, and output a reconstructed vector
$\tilde{\boldsymbol{x}}$ via its decoder $\boldsymbol{z}_{q}$. VQ-VAE outputs
a deterministic posterior distribution $q$ such that

$\displaystyle q(\boldsymbol{z}=k|\boldsymbol{x})=\begin{cases}1,&\text{ if
}k=\operatorname*{arg\,min}_{j}\left|\left|z_{e}(\boldsymbol{x})-e_{j}\right|\right|_{2}^{2}\\\
0,&\text{ otherwise }\end{cases}$ (2)

The decoder does not possess gradient and is trained by copying the gradients
from the encoder. The final output of the decoder is the log-posterior
probability $\log P(\boldsymbol{x}|\boldsymbol{z}_{q})$ which is part of the
reconstruction loss.

####  3.2.3 Architecture and Optimization.

Figure 1: Overview of the proposed system. $D$ denotes the discriminator, $G$
denotes the generator.

We propose a novel model for clustering expression profiles as shown in Figure
1. The model consists of a discriminator denoted as $D$ that maximizes the
mutual information and a generator $G$ that aims to reconstruct the input via
modeling a categorical underlying latent feature space spanned by
$\\{e_{i}\\}$. $D$ and $G$ are deeply coupled via the latent embeddings
$\boldsymbol{z}$, which is made possible through the fact the decoder of VQ-
VAE does not possess gradients and hence the embedding space can be controlled
by only the encoder and the discriminator. In prior work, the generator is
often architecturally independent from the discriminator and is only weakly
related through loss functions [13, 20, 28]. Intuitively, one can consider the
proposed model attempts to simultaneously minimize reconstruction loss as well
as maximize the mutual information:

$\displaystyle\mathcal{L}:=\underbrace{\widehat{\mathcal{H}}\left(P(y)\right)-{\widehat{\mathcal{H}}\left(P(y|\boldsymbol{z})\right)}-R(\lambda)}_{\mathcal{L}_{D}}+\underbrace{\log
P(\boldsymbol{x}|\boldsymbol{z}_{q})+\left|\left|\text{sg}[\boldsymbol{z}_{e}]-e\right|\right|_{2}+\left|\left|\boldsymbol{z}_{e}-\text{sg}[e]\right|\right|_{2}}_{\mathcal{L}_{G}}$
(3)

where $\mathcal{L}_{D},\mathcal{L}_{G}$ denote the discriminator loss and the
generator loss, respectively. $R(\lambda)$ is a possible regularizer that
controls the weight growth, e.g.
$R(\lambda):=\frac{\lambda}{2}||w^{T}w||^{2}_{2}$, where $w$ denotes the
weight parameters of the model. $\text{sg}[\cdot]$ denotes the stop gradient
operator.

####  3.2.4 Automatically Setting Number of Subtypes.

The proposed model can automatically determine suitable number of subtypes by
exploiting hidden information contained in the expression profiles which is
not available to conventional methods such as K-means relying on prior
knowledge. The automatic subtyping is made possible via the deeply coupled
latents and the discriminator: the multi-layer perceptron in the discriminator
outputs the logarithm of posterior distribution $\log
q(\boldsymbol{z}|\boldsymbol{x})$. However, by definition of Eq. (2) the
posterior is deterministic, which suggests $\log
q(\boldsymbol{z}|\boldsymbol{x})$ must either be $0$ or tend to $-\infty$. The
subsequent softmax layer hence outputs:

$\displaystyle
P(y=k|\boldsymbol{z})=\begin{cases}\frac{q(\boldsymbol{z}=k|\boldsymbol{x})}{\sum_{k=1}^{K}q(\boldsymbol{z}=k|\boldsymbol{x})},&\text{
if
}k=\operatorname*{arg\,min}_{j}\left|\left|z_{e}(\boldsymbol{x})-e_{j}\right|\right|_{2}^{2}\\\
0,&\text{ otherwise }\end{cases}$ (4)

We can set $K$ to a sufficient large integer $\tilde{K}$ initially that covers
the maximum possible number of subtypes. Since the nearest neighbor lookup of
VQ-VAE typically only updates a small number of embeddings $e_{j}$, by Eq. (4)
we see for any unused $e_{i},i\neq j$ the clustering probability is zero,
which suggests the number of subtypes $K$ will finally narrow down to a much
smaller number $K\ll\tilde{K}$.

##  4 Experiments

The expression profile data used in this study were collected from the world’s
largest cancer gene information database Genomic Data Commons (GDC) portal.
All of the used expression data were generated from cancer samples prior to
treatment.

We utilized the expression profiles of three representative types of cancer
for experiments:

  * •

Breast invasive carcinoma (BRCA): BRCA is the most prevalent cancer in the
world. Its expression profiles were collected from the Illumina Hi-Seq
platform and the Illumina GA platform.

  * •

Brain lower grade glioma (LGG): the expression profiles were collected from
the Illumina Hi-Seq platform.

  * •

Glioblastoma multiforme (GBM): the expression profiles were collected from the
Agilent array platform. Results on this dataset are deferred to the appendix.

These datasets consist of continuous-valued expression profiles (feature
length: 11327) of 639, 417 and 452 subjects, respectively. Additional
experimental results and hyperparameters can be seen in Appendix Section A
available at https://arxiv.org/abs/2206.10801.

The experimental section is organized as follows: we first compare the
clustering results with the ground truth labels $\mathcal{Y}_{gt}$ in Section
4.1 to validate the proposed method. We show in Section 4.2 that VQ-RIM
consistenly re-assigns subjects to different subtypes and produces one more
potential subtype with enlarged separation in between-group life expectancies,
which in turn suggests VQ-RIM is capable of better capturing the underlying
characteristics of subtypes. Extensive ablation studies on both the
categorical generator (VQ-VAE) and the information maximizing discriminator
(RIM) are performed to validate the proposed architecture in Section 4.3. We
believe the VQ-RIM subtyping result is far-reaching and can provide important
new insights to the unsettled debate on cancer subtyping.

###  4.1 Ground Truth Comparison

Figure 2: Comparison between $\mathcal{Y}_{gt}$ and the VQ-RIM label $y$ on
the Basal-like subtype of BRCA.

For validating the correctness of VQ-RIM, we show an example in Figure 2,
i.e., the Basal-like cancer subtype of BRCA that has been well-studied and
extensively validated by human experts and can be confidently subtyped, which
can be exploited as the ground-truth labels $\mathcal{Y}_{gt}$.

However, other subtypes lack such well-verified labels and are regarded as the
controversial labels $\mathcal{Y}_{c}$. The left subfigure of Figure 2 shows
the two principal axes of Basal-like expression profiles after PCA. The blue
triangles in the right subfigure indicates the difference between
$\mathcal{Y}_{gt}$ and the VQ-RIM result. It can be seen that VQ-RIM agrees
well with the ground truth.

###  4.2 Controversial Label Comparison

####  4.2.1 Subtype Comparison.

We compare existing controversial labels $\mathcal{Y}_{c}$ with the clustering
results of VQ-RIM in Figure 3. VQ-RIM output sensible decision boundaries that
separated the data well and consistently produced one more subtype than
$\mathcal{Y}_{c}$. As confirmed in Section 4.1, the Basal-like subtype
concorded well with the VQ-RIM Cluster A. On the other hand, other subtypes
exhibited significant differences: controversial labels seem to compactly fit
into a fan-like shape in the two-dimensional visualization. This is owing to
the human experts’ heuristics in subtyping: intuitively, the similarity of
tumors in the clinical variables such as morphological appearance often
renders them being classified into an identical subtype. However, cancer
subtypes are the result of complicated causes on the molecular level. Two main
observations can be made from the BRCA VQ-RIM label: (1) Luminal A was divided
into three distinct clusters C,D,E. Cluster E now occupies the left and right
wings of the fan which are separated by Cluster B and C; (2) A new subtype
Cluster F emerged from Luminal B, which was indistinguishable from Cluster E
if naïvely viewed from the visualization. This counter-intuitive clustering
result confirmed the complexity of cancer subtypes in expression profiles
seldom admits simple representations as was done in the controversial labels.
A similar conclusion holds as well for other datasets such as LGG: IDH mut-
codel was divided into two distinct subtypes (Cluster A, B), among which the
new subtype Cluster A found by VQ-RIM occupied the right wing of IDH mut-
codel. In later subsections, the one more cluster and re-assignment of VQ-RIM
are justified by analyzing the subtype population and from a medical point of
view. Due to page limit, we provide analysis focusing on BRCA only.

Figure 3: PCA visualization of the first two principal axes for BRCA and LGG.
Figure 4:  (a) t-SNE visualization of the sample distribution on BRCA. (b)
t-SNE of the samples with controversial labels. (c) label flows from the
controversial labels (left) to VQ-RIM 5 subtypes (mid) and 6 subtypes (right).

####  4.2.2 Label Flows.

The controversial labels might run the risk of _over-simplifying assignment_
which refers to that in the regions overlapped with several distinct subtypes,
controversial labels put all subjects into one of them without further
identifying their sources. Such assignment can be illustrated by Figure 4.
Here, Figure (4a) plots the sample distribution with darker colors indicating
denser population of samples. It is visible that the samples can be assigned
to five clusters. However, by injecting subtyping label information it is
clear from Figure (4b) that in the lower left corner there existed strong
overlaps of three different subtypes. Controversial labels assigned them to a
single subtype Luminal A. VQ-RIM, on the other hand, was capable of separating
those three subtypes. This separation can be seen from Figure (4c) which
compares the two labeling when setting the number of VQ-RIM subtypes to 5 in
accordance with controversial labels, or to 6 by setting $K$ to a sufficiently
large value and automatically determines the suitable number of subtypes. In
either case, VQ-RIM consistently separated the Luminal A into three distinct
subtypes: (B,C,E) in 5 subtypes case and (C*, D*, E*) in the 6 subtypes case.
In the next subsection, we verify the effectiveness of such finer-grained
subtyping by performing survival analysis.

####  4.2.3 Medical evaluation.

Figure 5: Kaplan-Meier survival analysis within each identified subtype group
(right) compared with original subtyping system (left) as a baseline. The line
in different colors represent patients from different subtypes. P-value was
calculated by Kaplan–Meier analysis with the log-rank test.

To demonstrate the clinical relevance of the identified subtypes, we perform
subtype-specific survival analysis by the Kaplan-Meier (KM) estimate. KM
estimate is one of the most frequently used statistical methods for survival
analysis, which we use to complementarily validate the VQ-RIM labels from a
clinical point of view [16]. KM compares survival probabilities in a given
length of time between different sample groups. The KM estimator is given by:
$\widehat{S}=\prod_{i:t{{}_{i}}<t}\frac{n_{i}-d_{i}}{n_{i}}$, where $n_{i}$ is
the number of samples under observation at time $i$ and $d_{i}$ is the number
of individuals dying at time $i$. The survival analysis graph is plotted
between estimated survival probabilities (on Y-axis) and the time passed after
samples entry into the study (on X-axis), where the survival curve is drawn as
a step function and falls only when a subject dies.

We can compare curves for different subtypes by examining gaps between the
curves in horizontal or vertical direction. A vertical gap means that at a
specific time point, samples belonging to one subtype had a greater fraction
of surviving, while a horizontal one means that it takes longer for these
samples to experience a certain fraction of deaths. The survival curves can
also be compared statistically by testing the null hypothesis, i.e. there is
no difference regarding survival situation among different groups, which can
be tested by classical methods like the log-rank test and the Cox proportional
hazard test.

Figure 5 shows the KM survival analysis graph for BRCA samples, based on the
PAM50 subtyping system and VQ-RIM subtypes. Compared with the PAM50, the
survival curves of VQ-RIM subtypes are more significantly separated. Log-rank
test also shows that there is significant difference in between-group survival
with a smaller $p$-value of 0.0023 compared against the PAM50
($\mathcal{Y}_{c}$). Smaller $p$-values indicate better subtyping results. We
indicate the subtype-specific median survival time with dashed lines. It is
visible that VQ-RIM performed better in identifying subtypes with large median
survival time differences.

###  4.3 Ablation Study

Figure 6: PCA visualization of the first two principal axes for BRCA using
different clustering methods. The number of cluster number is determined by
$\mathcal{Y}_{c}$ (PAM50). Figure 7: The number of clustering is determined by
combining the results from the elbow method, not used for RIM.

In this section, we conduct comprehensive ablation experiments to further
strengthen the effectiveness of VQ-RIM. Specifically, we validate the VQ part
and RIM part respectively by comparing VQ-RIM against the following
combinations:

  * •

we replace the VQ part with AE and VAEs with continuous latent which have been
exploited for subtyping in [31, 33]. The expression profiles are compressed
into continuous latent feature spaces for subsequent RIM clustering.

  * •

we replace the RIM part with existing classic clustering algorithms such as
K-Means, spectral clustering, and Gaussian mixture models [3]. Categorical
latent variables from VQ-VAE are fed into them for subtyping.

Results of the ablation studies can be seen by inspecting Figure 6 row-wise
and column-wise, respectively. All methods applied $\mathcal{Y}_{c}$ (PAM50)
for labeling. By inspecting the results row-wise, a first observation is that
for all clustering methods AE and VAE tended to disperse data points. On the
other hand, it can be seen from VQ-VAE clustering that subjects from distinct
subtypes compactly located in lower-dimensional spaces. Column-wise inspection
indicates that compared to other clustering methods, RIM tended to more
cohesively aggregate the in-group points: by contrast, it is visible from the
VAE row that only Basal-like subjects were cohesively grouped. Among all
compared methods, VQ-RIM stood out as the subjects were located in lower
dimensional spaces from which clear decision boundaries emerged. The clear
separation of VQ-RIM can be attributed to the underlying nearest neighbor
search: such search essentially performed a preliminary grouping on the data,
which greatly facilitated later separation. This observation is consistent
with the recent finding of pseudo-labeling that explicitly performs K-means in
creating preliminary clusters [20]. Besides the aforementioned qualitative
evaluation, we also quantitatively measure in Table 1 the scores of all the
clustering results in Figure 6 by using the three well-accepted metrics:
Normalized Mutual Information (NMI), Sihouette coefficient scores (Silhouette)
and $p$-value of survival analysis [17].

However, in Section 4.2, the labeling $\mathcal{Y}_{c}$ might not be the best
medically even if the clustering result accords well with human intuition. In
Figure 7 we focus on VQ-VAE, and set the number of clusters for RIM to a
sufficiently large value and let RIM automatically determines a suitable
number of subtypes. We term this strategy automatic VQ-VAE in the last row of
Table 1.

For other clustering algorithms, the number of clusters is determined from the
Silhouette coefficient scores and the elbow method [17]. It is visible that
clustering algorithms other than RIM tended to reduce the number of subtypes
for higher scores. By contrast, VQ-RIM produced one more subtype. This
automatic VQ-RIM clustering was superior from a medical perspective since it
achieved greatest subtyping result as demonstrated by the smallest $p$-value
of 0.0023. Furthermore, algorithmically it is better than plain VQ-RIM as it
achieved the highest scores of 0.63 and 0.54 of NMI and Silhouette among all
ablation choices.

Table 1:  Metrics used for measuring the clustering results of Figure 6 and 7.
The top three rows show the number of clustering determined by
$\mathcal{Y}_{c}$ (PAM50), while the last row shows the number of clustering
is automatically determined.

Generator |  Discriminator |  NMI $\uparrow$ |  Silhouette $\uparrow$ |  $p$-value $\downarrow$
---|---|---|---|---
AE + |  K-Means |  0.34 |  0.13 |  0.0861
Spectral clustering |  0.01 |  0.01 |  0.1523
Gaussian mixtures |  0.34 |  0.11 |  0.0734
RIM |  0.31 |  0.13 |  0.0834
VAE + |  K-Mmeans |  0.29 |  0.17 |  0.0812
Spectral clustering |  0.06 |  0.06 |  0.1382
Gaussian mixtures |  0.28 |  0.17 |  0.0899
RIM |  0.33 |  0.22 |  0.0732
VQ-VAE + |  K-Means |  0.33 |  0.29 |  0.0154
Spectral clustering |  0.05 |  0.04 |  0.0194
Gaussian mixtures |  0.44 |  0.24 |  0.0166
RIM |  0.55 |  0.29 |  0.0042
VQ-VAE + (Automatic) |  K-Means |  0.42 |  0.29 |  0.0145
Spectral clustering |  0.06 |  0.05 |  0.0188
Gaussian mixtures |  0.51 |  0.32 |  0.0132
RIM |  0.63 |  0.54 |  0.0023

##  5 Discussion and Conclusion

In this paper we were concerned with the cancer subtyping problem that aimed
to disclose the difference between subtypes within a specific cancer. Existing
literature is having an unsettled debate over the subtyping problem, with
various definition and suitable number of subtypes put forward from distinct
viewpoints. Aiming to aid the medical experts by providing dependable
reference for subtyping, we took a data-scientific standpoint and exploited
genetic expression profiles of cancers without using the controversial labels
current literature has imposed. Such genetic expression profiles featured
scarcity, high dimensionality, and complicated dependence which posed
challenges for both physicians and data scientists. To tackle these problems,
we leveraged information-theoretic principles as well as recent categorical
latent generative modeling techniques that helped in minimizing clustering
confusion and maximizing interpretability. The resultant novel model: Vector
Quantization Regularized Information Maximization (VQ-RIM) can better reveal
the intrinsic difference between cancer genetic expression profiles and based
on which automatically decide a suitable number of subtypes. The experiment on
ground-truth BRCA cancer verified the correctness of VQ-RIM, while more
extensive experiments on multiple authoritative datasets consisting of various
cancers showed the difference between VQ-RIM results and the controversial
labels. By comprehensive analysis from both data scientific and medical views,
we illustrated that the different subtyping result yielded by VQ-RIM
consistently outperformed existing ones in terms of survival analysis, and
contributed important new insights into the unsettled debate.

The future work consists of two interesting directions: (1) to further
validate the effectiveness of VQ-RIM, comprehensive experiments on all
available cancer datasets and comparison with their existing labeling might be
necessary. (2) the VQ-RIM architecture might not only work well with cancer
data but also be generally applicable on radically different data such as
images, voices that inherently exploit discrete nature of the data.

##  6 Acknowledgement

This work was supported by JST Mirai Program (JPMJMI20B8) and JST PRESTO
(JPMJPR21C7), Japan.

## References

  * [1] Cancer today - iarc. https://gco.iarc.fr/today/home, accessed: 2021
  * [2] Alemi, et al.: Fixing a broken ELBO. In: ICML-18. pp. 159–168 (2018)
  * [3] Alexe, G., et al.: Analysis of breast cancer progression using principal component analysis and clustering. Journal of biosciences pp. 1027–1039 (2007)
  * [4] Bair, E., Tibshirani, R.: Machine learning methods applied to dna microarray data can improve the diagnosis of cancer. SIGKDD Explor. Newsl. p. 48–55 (2003)
  * [5] Berger, et al.: A comprehensive pan-cancer molecular study of gynecologic and breast cancers. Cancer cell pp. 690–705 (2018)
  * [6] Boudiaf, M., Ziko, I., Rony, J., Dolz, J., Piantanida, P., Ben Ayed, I.: Information maximization for few-shot learning. In: NIPS-20. pp. 2445–2457 (2020)
  * [7] Bridle, J., Heading, A., MacKay, D.: Unsupervised classifiers, mutual information and 'phantom targets. In: NIPS. pp. 1096–1101 (1991)
  * [8] Chen, X., Kingma, D.P., Salimans, T., Duan, Y., Dhariwal, P., Schulman, J., Sutskever, I., Abbeel, P.: Variational lossy autoencoder. In: ICLR-17 (2017)
  * [9] Dagdia, et al.: Rough set theory as a data mining technique: A case study in epidemiology and cancer incidence prediction. In: ECML-PKDD. pp. 440–455 (2018)
  * [10] Ghasedi Dizaji, K., Wang, X., Huang, H.: Semi-supervised generative adversarial network for gene expression inference. In: KDD ’18. p. 1435–1444 (2018)
  * [11] Grandvalet, Y., Bengio, Y.: Semi-supervised learning by entropy minimization. In: NIPS-04. p. 529–536 (2004)
  * [12] Heiser, L.M., Sadanandam, A., et al.: Subtype and pathway specific responses to anticancer compounds in breast cancer. PNAS (8), 2724–2729 (2012)
  * [13] Hu, W., et al.: Learning discrete representations via information maximizing self-augmented training. In: ICLR-17. p. 1558–1567 (2017)
  * [14] Jiang, B., Ding, C., Bin, L.: Covariate-correlated lasso for feature selection. In: ECML-PKDD ’14. pp. 595–606 (2014)
  * [15] Jolliffe, I.T., Cadima, J.: Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences (2065), 20150202 (2016)
  * [16] Kaplan, E.L., Meier, P.: Nonparametric estimation from incomplete observations. Journal of the American statistical association (282), 457–481 (1958)
  * [17] Kaufman, L., Rousseeuw, P.J.: Finding groups in data: an introduction to cluster analysis. John Wiley & Sons (2009)
  * [18] Kingma, D.P., Welling, M.: Auto-Encoding Variational Bayes. In: ICLR-14. pp. 1–9 (2014)
  * [19] Krause, A., Perona, P., Gomes, R.: Discriminative clustering by regularized information maximization. In: NIPS-10. p. 775–783 (2010)
  * [20] Liang, et al.: Do we really need to access the source data? Source hypothesis transfer for unsupervised domain adaptation. In: ICML-20. pp. 6028–6039 (2020)
  * [21] Liang, C., Shang, M., Luo, J.: Cancer subtype identification by consensus guided graph autoencoders. Bioinformatics (24), 4779–4786 (2021)
  * [22] Liu, C.C., et al.: LRP6 overexpression defines a class of breast cancer subtype and is a target for therapy. PNAS (11), 5136–41 (2010)
  * [23] Löwe, S., O' Connor, P., Veeling, B.: Putting an end to end-to-end: Gradient-isolated learning of representations. In: NIPS-19. pp. 1–13 (2019)
  * [24] Mangan, S., Alon, U.: Structure and function of the feed-forward loop network motif. PNAS (21), 11980–11985 (2003)
  * [25] van den Oord, A., Vinyals, O., Kavukcuoglu, K.: Neural discrete representation learning. In: NIPS-17. p. 6309–6318 (2017)
  * [26] de Ronde, J., Wessels, L., Wesseling, J.: Molecular subtyping of breast cancer: ready to use? The Lancet. Oncology (4), 306–307 (2010)
  * [27] Slonim, N., Atwal, G.S., Tkačik, G., Bialek, W.: Information-based clustering. Proceedings of the National Academy of Sciences, PNAS (51), 18297–18302 (2005)
  * [28] Springenberg, J.T.: Unsupervised and semi-supervised learning with categorical generative adversarial networks. In: ICLR-16. pp. 1–12 (2016)
  * [29] Tschannen, M., Djolonga, J., Rubenstein, P.K., Gelly, S., Lucic, M.: On mutual information maximization for representation learning. In: ICLR-20. pp. 1–12 (2020)
  * [30] Weinstein, et al.: The cancer genome atlas pan-cancer analysis project. Nature genetics (10), 1113–1120 (2013)
  * [31] Withnell, et al.: XOmiVAE: an interpretable deep learning model for cancer classification using high-dimensional omics data. Briefings in Bioinformatics (6) (2021)
  * [32] Xi, J., et al.: Tolerating data missing in breast cancer diagnosis from clinical ultrasound reports via knowledge graph inference. In: KDD ’21. p. 3756–3764 (2021)
  * [33] Yang, B., Xin, T.T., Pang, S.M., Wang, M., Wang, Y.J.: Deep Subspace Mutual Learning for cancer subtypes prediction. Bioinformatics (21), 3715–3722 (2021)

##  Appendix 0.A Parameter Settings and Training

The details of the parameter settings are shown in Table 2. To make the utmost
of the model, a grid search of hyperparameters was implemented in this work to
seek the best combination. Note that the optimal settings (values in Table 2)
were used both for all ablation studies. The experiments were conducted on a
server with an NVIDIA GeForce RTX 3090Ti GPU.

Table 2: Parameter setting of experiments Model Parameter |  Value
---|---
#Embedding |  64
Dimension of encoder |  512
Dimension of embedding |  64
Commitment cost |  1
Dropout rate |  0.5
Training settings |  Value
#Training epoch |  200
Batch size |  32
Optimizer |  AdamW
Learning rate |  $e\times 10^{-5}$
Figure 8: PCA visualization of the first two principal axes for LGG.

This work is conducted a pre-training-to-fine-tuning strategy in the training
process. Specifically, during the pre-training phase the all cancer sample
from 4 cancers was used to initialize the model parameters of VQ-VAE. The
$Adam$ with a biggish learning rate of $10^{-4}$ was utilized for the
optimization. In the fine-tuning phase, we implemented RIM on pre-trained VQ-
VAE to conduct the subtyping tasks for different cancer datasets. Here, the
$AdamW$ optimizer was used to meticulously optimize the loss of subtyping.

##  Appendix 0.B Complementary Results

This section basically shows the results similar to Section 4 of the main
paper. Figure 8 is a reminiscent of Figure 3 in the main paper and shows a
comparison between labeling $\mathcal{Y}_{c}$ and the number of subtypes
determined automatically by VQ-RIM. Figure 9 plots the density population
density and label flows from $\mathcal{Y}_{c}$ and VQ-RIM which has one more
subtype. Figure 10 shows a KM analysis indicating the VQ-RIM result achieved
clearer separation in the life expectancy with one more subtype. The same
conclusion holds for other datasets as well. For simplicity we only show the
clustering result for GBM in Figure 11.

From Figure 8, IDH mut-codel was divided into two distinct subtypes (Cluster
A, B), among which the new subtype Cluster A found by VQ-RIM occupied the
right wing of IDH mut-codel. In later subsections, the one more cluster and
re-assignment of VQ-RIM are justified by analyzing the subtype population and
from a medical point of view.

Figure 9: (a) t-SNE visualization of the sample distribution on LGG. (b) t-SNE
of the samples with controversial labels. (c) label flows from the
controversial IDH cluster labels (left) to VQ-RIM 4 subtypes (right).

Figure (9a) plots the sample distribution with darker colors indicating denser
population of samples. It is visible that the samples can be assigned to three
clusters. From Figure (9b) that we can find in the lower left corner there
existed a overlaps of two different subtypes. Controversial labels assigned
them to a single subtype IDH mut-codel. VQ-RIM was capable of separating those
two subtypes. This separation can be seen from Figure (9c) which compares the
two labeling when setting the number of VQ-RIM subtypes to 4 by setting $K$ to
a sufficiently large value and automatically determines the suitable number of
subtypes. In either case, VQ-RIM consistently separated the IDH mut-codel into
two distinct subtypes: (A* and B*) in 4 subtypes case.

Figure 10: Kaplan-Meier survival analysis within each identified subtype group
(right) compared with original subtyping system (left) as a baseline. The line
in different colors represent patients from different subtypes. P-value was
calculated by Kaplan–Meier analysis with the log-rank test. Figure 11: PCA
visualization of the first two principal axes for GBM. It illuminates the
proposed method has better clustering performance than others.

Figure 10 shows the KM survival analysis graph for LGG samples, based on the
IDH subtyping system and VQ-RIM subtypes. Compared with the IDH, the survival
curves of VQ-RIM subtypes are more significantly separated. Log-rank test also
shows that there is significant difference in between-group survival with a
smaller $p$-value of 0.0124 compared against the IDH with 0.0525.


[3,4]Danfeng Li [1]Alexandre Gloter

1]Laboratoire de Physique des Solides, CNRS, Université Paris-Saclay, Orsay,
91400, France 2]Synchrotron SOLEIL, L’Orme des Merisiers, BP 48 St Aubin, Gif
sur Yvette, 91192, France 3]Department of Physics, City University of Hong
Kong, Kowloon, Hong Kong

4]City University of Hong Kong Shenzhen Research Institute, Shenzhen,
Guangdong, 518057, China

5]CPHT, Ecole Polytechnique, Palaiseau cedex, 91128, France

6]Department of Quantum Matter Physics, University of Geneva, Geneva,
Switzerland

7]Electron Spectrometry and Microscopy Laboratory (LSME), Institute of Physics
(IPHYS), Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne,
Switzerland

8]LCPMR, Sorbonne Université, CNRS, Paris, 75005, France

9]These authors contributed equally to this work

# Emergent electronic landscapes in a novel valence-ordered nickelate with
tri-component nickel coordination

Aravind Raji  Zhengang Dong  Victor Porée  Alaska Subedi  Xiaoyan Li  Bernat
Mundet  Lucia Varbaro  Claribel Domínguez  Marios Hadjimichael  Bohan Feng
Alessandro Nicolaou  Jean-Pascal Rueff<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [  [  [  [  [

###### Abstract

The metal-hydride-based ‘topochemical reduction’ process has produced novel
thermodynamically unstable phases across various transition metal oxide series
with unusual crystal structures and non-trivial ground states. Here, by such
an oxygen (de-) intercalation method we synthesis a novel samarium nickelate
with ordered nickel valences associated with tri-component coordination
configurations. This structure, with a formula of Sm9Ni9O22 as revealed by
four-dimensional scanning transmission electron microscopy, emerges from the
intricate planes of {303}${}_{\text{pc}}$ ordered apical oxygen vacancies.
X-ray spectroscopy measurements and ab-initio calculations show the
coexistence of square-planar, pyramidal and octahedral Ni sites with mono-,
bi- and tri-valences. It leads to an intense orbital polarization, charge-
ordering, and a ground state with a strong electron localization marked by the
disappearance of ligand-hole configuration at low-temperature. This new
nickelate compound provides another example of previously inaccessible
materials enabled by topotactic transformations and presents a unique platform
where mixed Ni valence can give rise to exotic phenomena.

###### keywords:

Topotactic transitions, Nickelates, 4D-STEM, RIXS, Charge order, Oxygen
vacancies, Carrier localization

## Introduction

The on-demand design of transition-metal oxides (TMO) with emerging properties
is imparted by the multivalent nature of the transition-metal ions and the
accessible complexity of lattice structures [1, 2]. Varying the structure and
constituent elements in TMO allows access to a range of competing ground
states and exotic electronic landscapes, in which the interplay of the
transition-metal orbitals and the connecting oxygen network often plays an
essential role [3]. One phenomenal example is the cuprate family, where
different coordinational Cu-O frameworks lead to distinctive $T_{\textrm{c}}$
($T_{\textrm{c}}$, superconducting transition temperature) values [4]. In
addition, the fast developments in topotactic soft-chemistry approaches have
enabled the synthesis of unconventional TMO with unusual valence states and/or
oxygen connectivity [5, 6, 7]. One of the recent seminal examples is the
synthesis of a family of superconducting low-valence nickelates [8, 9, 10, 11]
by a topotactic reduction process, which displays intriguing analogies and
distinctions to the high-$T_{\textrm{c}}$ cuprates and is the focus of the
current extensive research.

Nickel oxides or nickelates form a large family of structures with diverse
crystalline configurations and electronic properties. Perovskite rare-earth
nickelates RNiO3 have a rich phase diagram generated by changing the rare-
earth cation. They exhibit metal-to-insulator, magnetic and structural
transitions [12, 13, 14], and upon stoichiometry alteration such as oxygen
deintercalation or protonization can become reconfigurable neuromorphic
devices [15]. Upon a strong topotactic oxygen reduction, RNiO2 infinite-layer
nickelates can be obtained, which have recently been discovered as
superconducting after hole-doping [8]. While the origin of most of these
properties is not yet fully understood, the interplay of spin, charge and
orbital degrees of freedom is known to drive these exotic properties, along
with the associated crystalline configurations.

In the case of infinite-layer nickelate superconductors, a structurally-
perfect infinite-layer thin film is essential to achieve superconductivity
[16]. The topotactic reduction method for infinite-layer nickelate synthesis
involves metal hydride, e.g. CaH2, to remove the apical oxygens from a
perovskite nickelate. This often results in a mixture of secondary phases with
partial removal of apical oxygen [17, 18]. This is due to the partial
instability attributed to the infinite-layer structure that does not host Ni2+
(3d8), but a Ni1+ in a 3d9 configuration [19]. Along with superconductivity,
the observed charge ordering (CO) [20, 21, 22, 23] in these systems now
appears to be connected to a structural reordering due to oxygen intercalation
in the apical sites [24, 25]. Nevertheless, the structure of the charge-
ordered systems cannot be directly related to known phases, such as
brownmillerite [17] or octahedra and square planar bearing R3Ni3O7 phase where
Ni2+ is absent [26]. The complex transport properties and notably the unusual
dramatic resistivity increase upon oxygen deintercalation of the perovskite
nickelates RNiO3 is also intriguing [27].

This concomitance indicates the possible existence of an unknown stable
intermediate phase with multivalent, strong electron localisation and charge
ordered ground state between the perovskite ABO3 and the infinite-layer ABO2.

Here, we demonstrate that the topochemical reduction of a perovskite SmNiO3
thin film grown on an orthorhombic NdGaO3 substrate stabilizes to such an
intermediate phase whose stoichiometry is derived as Sm9Ni9O22. This
ubiquitous phase existing between the much pronounced perovskite ABO3 and the
infinite-layer ABO2 nickelates hosts intriguing properties associated with
multi-valent and multi-component coordinated Ni sites. Intriguing properties
of this phase include a strong-orbital polarization, temperature dependent
carrier localization, and charge ordering. The observed structural and
electronic properties from microscopy and spectroscopy are validated by ab-
initio simulations, contributing to a coherent analysis and uncovering some
regions of its possibly rich phase-diagram.

## Results

### A new phase from the {303}${}_{\text{pc}}$ apical oxygen vacancy ordering

We begin with a macroscopic comparison between the parent perovskite SmNiO3
and the reduced Sm9Ni9O22 thin-films, using x-ray diffraction (XRD) and
transport measurements. As shown in Fig. 1a, we find that the reduction of
SmNiO3 results in a scaling down of the out-of-plane (o-o-p) lattice parameter
from 3.78 Å to 3.63 Å, indicative of a structural transition to the reduced
Sm9Ni9O22 phase. The thin-film quality is maintained, as evidenced by the
intense Bragg peaks and finite-size oscillations for both phases.
Additionally, in the transport measurements shown in Fig.1b, SmNiO3 evidences
a semiconducting behavior with an increase in conductivity of 4 orders of
magnitude between 75K and room-temperature, and a metal-to-insulator
transition around 375K. This is in agreement with previous reports on SmNiO3
thin-films [28]. On the other hand, the resistivity of Sm9Ni9O22 was above the
measurement range $>2\times 10^{10}$ $\mu$$\Omega$.cm making it highly
insulating with respect to its parent phase.

The microscopic characterization begins with a structural comparison between
the parent perovskite SmNiO3 and the reduced Sm9Ni9O22 thin-films, using high-
angle annular dark-field (HAADF)-STEM imaging, geometrical phase analysis
(GPA) [29], and Fourier transform (FT) analysis. As shown in Fig.1c & d, the
HAADF-STEM image shows homogeneous thin films for both SmNiO3 and Sm9Ni9O22,
with a coherent growth obtained on the (110)${}_{\text{o}}$ cut NdGaO3
substrate (where subscript ’o’ stands for orthorhombic unit cell in Pbnm space
group with a = 5.433 Å, b = 5.504 Å, c= 7.716 Å ). The low magnification
HAADF-STEM image of both are given in Supplementary Information, Fig. S1 & S2.

Figure 1: Structural and transport comparison of SmNiO3 with Sm9Ni9O22 by XRD,
HAADF-STEM imaging, GPA, and FT. (a) XRD comparison of both. The 001 and 002
peaks for SmNiO3 occur at lower 2$\theta$ values compared to the 001 and 002
peaks of Sm9Ni9O22, indicating an out-of-plane compression of the pseudocubic
(pc) unit cell upon reduction. (b) Transport measurement of parent SmNiO3 in
increasing and decreasing temperatures. The resistivity curve of Sm9Ni9O22 is
not included in the plot since it was significantly higher, and is beyond the
detection range of the equipment ($>2\times 10^{10}$ $\mu\Omega$.cm). (c)
HAADF-STEM image of SmNiO3 thin film. (d) HAADF-STEM image of Sm9Ni9O22 thin
film, showing in good contrast the stripes coinciding with the
(303)${}_{\text{pc}}$ direction. (e & f) Maps of the o-o-p strain by GPA in
(e) SmNiO3 and (f) Sm9Ni9O22. Both samples show uniform o-o-p compression. The
Sm9Ni9O22 exhibits a more negative $\epsilon_{zz}$ compared to SmNiO3,
corresponding to a lower o-o-p parameter. (g & h) The corresponding FT of the
HAADF image of (g) SmNiO3, (h) Sm9Ni9O22. The NGO substrate is indexed with an
orthorhombic unit cell and the super-structure spots associated to the tilt
pattern of the NGO are encircled in red. The super-structure spots associated
to the (303)${}_{\text{pc}}$ stripes in Sm9Ni9O22 are encircled in yellow and
are indexed within a pseudo-cubic unit cell.

Fig.1e& f show the maps of o-o-p strain from GPA analysis. The SmNiO3 films
have a uniform o-o-p reduction of 3% as compared to the substrate,
corresponding to an o-o-p parameter of 3.78 Å, while the Sm9Ni9O22 films have
a uniform o-o-p parameter of 3.63 Å. Such uniformity indicates a high
crystalline quality devoid of any defects, which is also visible at larger
field of view (see Supplementary Information, Fig. S1). Furthermore, the
Sm9Ni9O22 have faint but robust contrasts corresponding to
(303)${}_{\text{pc}}$ (subscript pc refers to the pseudocubic unit cell)
oriented stripes which result in superstructure spots in the Fourier transform
(FT) shown in Fig.1h at Q = ($\pm\frac{1}{3},0,\pm\frac{1}{3}$),
($\pm\frac{2}{3},0,\pm\frac{2}{3}$), ($\pm\frac{1}{3},0,\mp\frac{2}{3}$),
($\pm\frac{2}{3},0,\mp\frac{1}{3}$) r.l.u. Here we define the reciprocal
lattice units (r.l.u.), with in plane components a = b = 3.86 Å, and the o-o-p
lattice constant c = 3.63 Å for the Sm9Ni9O22.

Complementary 4D-STEM measurements were carried out to further structurally
characterize this new phase (Sm9Ni9O22) and determine the local structural
modifications at the atomic level. Fig.2 shows the HAADF and 4D-STEM results
obtained along different crystallographic orientations of the Sm9Ni9O22. In
Fig.2a,e are the Sm9Ni9O22 HAADF images at higher magnification, corresponding
to a projection of the NGO substrate along the [001]${}_{\text{o}}$ zone axis.
This is denoted as the 0∘ projection, or [010]${}_{\text{pc}}$ zone axis, for
the new Sm9Ni9O22 phase. The HAADF image shows the extended and reduced Sm-Sm
distances, and one third of the Sm column exhibit a doubling, with all these
reconstructions running periodically along the {303}${}_{\text{pc}}$ stripes.

Figure 2: 4D-STEM analysis of Sm9Ni9O22 in 0∘, 90∘, and 45∘ orientations that
corresponds to [010], [100], and [110]${}_{\text{pc}}$ viewing axis. (a) HAADF
image of Sm9Ni9O22 viewed along [010], showing doubled rare-earth sites and
elongated - shortened Sm-Sm distances. (b) The COM image of the same region,
obtained by 4D-STEM analysis. (c) The corresponding dCOM image, with a very
good oxygen phase contrast, that shows Vo line in-between the extended Sm-Sm
distance, hence along the stripes. (d) A relaxed structural model that is
identical to the dCOM image in (c). (e) A magnified HAADF image of region
highlighted in (a), indicating the location of each atom, and showing the
extended and reduced Sm-Sm distances. (f) A magnified dCOM image of region
highlighted in (c), indicating the location of each atom, and periodic Vo. (g)
HAADF image of Sm9Ni9O22 viewed along [100], showing features similar as in
the view along [010]. (h) The corresponding dCOM image. (i) HAADF image of
Sm9Ni9O22 viewed along [110]. (j) The corresponding dCOM image. (k) A low
magnification HAADF image of Sm9Ni9O22 viewed along [100], that shows stripes
along (033)${}_{\text{pc}}$ and (0$\bar{3}$3)${}_{\text{pc}}$directions. (l)
The corresponding FT showing spots at Q = (0,
$\pm\frac{1}{3},\pm\frac{1}{3}$), (0, $\pm\frac{2}{3},\pm\frac{2}{3}$), (0,
$\pm\frac{1}{3},\mp\frac{2}{3}$), (0, $\pm\frac{2}{3},\mp\frac{1}{3}$) r.l.u.
and Q = (0, $\pm\frac{1}{3},\pm\frac{2}{3}$), (0,
$\pm\frac{2}{3},\pm\frac{1}{3}$), (0, $\pm\frac{1}{3},\mp\frac{1}{3}$), (0,
$\pm\frac{2}{3},\mp\frac{2}{3}$) r.l.u. (m) The inverse-FT image showing the
two domains obtained from the highlighted spots in (l).

4D-STEM data collected simultaneously with the HAADF images are shown in
Fig.2b,c,f. Figure 2b corresponds to the displacement intensity of the center-
of-mass (COM) of the diffraction pattern measured at each pixel, and donut
shapes reveal the presence of atomic columns at their centers, even for oxygen
atoms [30, 31]. The figures 2c,f are obtained by the divergence of the center-
of mass (dCOM) approximating a projected charge-density image, with extrema
located at the atomic positions. These images reveal unambiguously that
columns of apical oxygen vacancies (Vo) are ordered along the stripes. This is
causing the rare-earth super-structure, since the Sm atoms are moving away
from the Vo sites giving separations between the Sm atoms of $\approx$ 4.20 Å
as shown in Fig.2e (while a reduced Sm-Sm distance of $\approx$ 3.10 Å is
measured when columns of apical oxygen are preserved). Such local variations
in the A-cation distances, arising due to the coexistence of different
coordinations of B sites have been observed in studies of cobaltite[32],
ferrite[33], and manganite[34] systems.

A structural model showing a similar atomic reconstruction is shown in Fig.2d.
This was obtained by an ab-initio crystal relaxation of a cubic SmNiO3, after
introducing a single family of (303)${}_{\text{pc}}$ planes of apical Vo. The
presence of such planes of Vo leads to a mixture of octahedral and pyramidal
Ni sites, and it forms the n = 3 member of a family made of
(n0n)${}_{\text{pc}}$ planes of Vo, with a formula AnBnO3n-1. Interestingly,
infinite-layer nickelates, with the formula RNiO2 form the n = 1 member. The
structural details of such AnBnO3n-1 series as predicted by ab-initio
calculations are in Supplementary Information, Fig.S3 and S4.

The STEM images measured along a projection rotated by 90° compared to Fig.2a
(corresponding to a substrate [$\bar{1}$10]${}_{\text{o}}$ zone axis and a
[100]${}_{\text{pc}}$ Sm9Ni9O22 zone axis ) as shown in Fig.2g and h indicate
that the vacancy ordering is taking place along two perpendicular faces. While
the observation along the 0° projection of the Sm9Ni9O22 was only showing
(303)${}_{\text{pc}}$ planes (e.g. and not a ($\bar{3}$03)${}_{\text{pc}}$),
the Fig.2k-m shows different domains in 90∘ projections, where the stripes
cpincide with the (033)${}_{\text{pc}}$ and the (0$\bar{3}$3)${}_{\text{pc}}$
planes. This indicates that 0∘ and 90∘ projections are not exactly
symmetrical, which is most likely due to the orthorhombicity of the substrate,
that causes an asymmetry between the two in-plane directions.

Nevertheless, locally both orientations reveal very similar cation and oxygen
atomic reconstruction of the Sm9Ni9O22. This essentially identifies that this
phase is obtained by intersecting apical Vo planes of (303)${}_{\text{pc}}$
and (033)${}_{\text{pc}}$ orientations. Such a possible tetragonal symmetry is
also confirmed by the observation along 45∘ projection as shown in Fig.2i and
j. The apical oxygen removal indicatively reduces the o-o-p parameter, as is
the case for infinite-layer nickelates, where the complete removal of apical
oxygen results in an o-o-p parameter of $\approx$ 3.30 Å [8, 35, 36, 37]. The
observed o-o-p parameter of 3.63 Å in Sm9Ni9O22 is thus associated with these
oxygen deficiencies.

### Structural and charge ordering from ab-initio calculations

From these microscopic results, it is evident that this new phase cannot be
simply explained by one family of Vo planes, and an intersecting
(303)${}_{\text{pc}}$ and (033)${}_{\text{pc}}$ planes of apical vacancies are
needed. Corresponding simulations involved the structural relaxations of a
6x6x6 cubic SmNiO3 supercell with (303)${}_{\text{pc}}$ and
(033)${}_{\text{pc}}$ ordered apical Vo. The calculations resulted in a
pseudo-tetragonal structure with an o-o-p parameter of $c=3.63$ Å. This 6x6x6
tetragonal supercell was then symmetrized, giving an 80-atom structure with
the formula Sm18Ni18O44, as shown in Fig.3a, which can be further shortened as
Sm9Ni9O22. More details are given in Methods and Supplementary Information,
Fig S6 and table S1. This structure essentially identifies three types of Ni
sites: NiO6 octahedral (Oh), NiO4 square-planar (Sq) and NiO5 pyramidal (Py)
in the ratio 1:2:6 in the lattice. In this regard, it is very different from
the A3B3O7 [26], which have the Vo in the periodicity $\frac{1}{3}$ from the
parent ABO3 structure and is built with apical Vo lines running only along one
face. It is here also noted that an alternative A3B3O7 structure can also been
obtained by mixing apical Vo (303)${}_{\text{pc}}$ planes and equatorial Vo
(033)${}_{\text{pc}}$ planes, but the resulting structure would not match our
observation (See Supplementary Information, Fig. S5).

Figure 3: Structural and electronic refinement by DFT. (a) Relaxed structural
model of Sm9Ni9O22, with Ni sites in a mixed coordination of pyramidal, square
planar and octahedral. The axis indicated in this figure is the local axis
with respect to the structure symmetry. (b) The polyhedral distribution viewed
along [010]${}_{\text{pc}}$ and a 10∘ rotation along [100], with the area
matching with that of the dCOM image. (c) The polyhedral distribution viewed
along [100]${}_{\text{pc}}$ and a 10∘ rotation along [010]. (d) The polyhedral
distribution viewed along [110]. From viewing along the three projections, the
different polyhedra appear to have a periodic distribution. (e) Calculated 3d
orbital occupancy of different Ni sites for U = 0 to 10 eV by DFT+U. (f)
Calculated net spin polarization of different Ni sites for U = 0 to 10 eV by
DFT+U. (g) A demonstration of the d-orbital occupancy of different Ni sites.
Ni in the pyramidal and octahedral sites are in a 3d8 configuration, while Ni
in the square planar sites are in a 3d9 configuration.

On the grounds of differently coordinated Ni sites in the lattice, a mixed
valence is expected for the Sm9Ni9O22 composition. Ab-initio calculations give
evidence of different Ni electronic configurations for the different site
coordination as indicated by Fig.3e-g. Calculations have been done including a
Hubbard U parameter ranging from 0 to 10 eV as shown in Fig.3e,f. At U = 0 eV,
differences of total charges on the Ni sites are already present, but become
more segregated between the square and the octahedra or the pyramid when U
increases. The Ni in square planar sites are in a 3d9 configuration,
indicating Ni1+ and a spin = 1/2 (typical occupation matrix is given in
Supplementary Information, Table S2). The Ni in the pyramidal and octahedral
sites are closer to a 3d8 configuration and a spin = 1. A Sm9Ni9O22
composition would imply 17 charges to be distributed over the nickel sites and
a possible distribution would then be 2 sq(Ni+, d9), 6 Py(Ni2+, d8) and a
single Oh (Ni3+, 3d8$\underline{L}$ with $\underline{L}$ referring to a hole
in the ligand). The exact nature of the Ni3+ ground state (i.e. the respective
weight of 3$d^{7}$ and 3$d^{8}\underline{L}$ in the ground state) is difficult
to address solely based from the DFT+U level where, in brief, electron
localization and correlation are brutally parameterized. We also assume that a
ligand hole will reside on the octahedral sites by analogy from the electronic
structure of the parent SmNiO3 phase. In the metallic state, it is known to
host octahedral Ni3+ with primary covalent ground state (a.k.a, as negative
charge transfer), id est, primarily composed by a 3d8$\underline{L}$. In the
insulating state, the bond disproportionates with the presence of collapsed
$d^{8}\underline{L}^{2}$ and expanded $d^{8}$ octahedral sites [38]. These
charge distribution and orbital physics are summarized in figure 3g, and they
are further probed by spectroscopic techniques involving x-ray absorption
spectroscopy (XAS), resonant inelastix x-ray scattering (RIXS) and hard x-ray
photoemission spectroscopy (HAXPES).

### Orbital polarization and charge localization by spectroscopic measurements

The schematic of the XAS and RIXS measurements is briefly indicated in Fig.4a.
The XAS Ni-L3 edge of the parent perovskite SmNiO3 thin film measured at 25K
is given in Fig.4b which shows a typical splitting associated to bond
disproportionation and the presence of collapsed $d^{8}\underline{L}^{2}$ and
expanded $d^{8}$ octahedra for the insulating state [38]. A more complete set
of the XAS and RIXS data at the Ni-L edges measured for the parent SmNiO3 thin
film at RT (300K) and 25K can be seen in Supplementary Information, Fig.S7.
They are in agreement with previously reported spectroscopic work on
perovskite nickelate thin films [39], and they show a fully developed
insulating behavior at 25K, and a partial metal-to-insulator transition
occurring at RT (300K) as revealed by a blurring of the XAS Ni-L3 edges
together with a blue shift of the localized dd-orbital excitations in the RIXS
[39].

Figure 4: Electronic structure analysis by XAS and RIXS. (a) The experimental
geometry showing the incident x-ray polarization $\sigma$ and $\pi$, incident
x-ray angle $\theta_{i}$ and the scattering angle 2$\Theta$.(b & c) XAS
partial fluorescence yield (PFY) spectra at 25 K at Ni-L3 edge of (b) SmNiO3
(c) Sm9Ni9O22. (d & e) Ni-L3 edge RIXS map at 25 K of Sm9Ni9O22 in (d)
$\sigma$ \- incident polarization and (e) $\pi$ \- incident polarization of
the photons. (f) Difference of RIXS map between the two polarizations. (g & h)
Calculated Ni-L3 edge RIXS map of Ni2+ pyramidal sites, in (g) $\sigma$
polarization and (h) $\pi$ polarization.(i) Difference of calculated RIXS map
between the two polarizations.

The Ni-L3 edge XAS for the Sm9Ni9O22 sample (Fig.4c) occurs at a lower
excitation energy than that of SmNiO3 indicating a reduced valence toward Ni2+
and Ni1+ as it was inferred from the structural investigation. In complement,
a Ni 2p${}_{\frac{3}{2}}$ and Ni 2p${}_{\frac{1}{2}}$ core level bulk
sensitive HAXPES measurements between Sm9Ni9O22 and SmNiO3 are given in
(Supplementary Information, Fig. S13) that show a similar valence trend. In
the XAS, a strong dichroism is observed at the Ni-L edges of the Sm9Ni9O22
with a noticeable o-o-p higher-intensity at the low-energy side of the Ni-L3
main edge (at 852.5 eV). From atomic multiplet calculations (shown in
Supplementary Information, Fig.S14) [40], this strong dichroism is considered
to be stemming from the Ni2+ in the pyramidal sites, which are the most
populated in the lattice. By projecting the XAS final-state 2p53d9 in terms of
crystal field configuration, the features at 852.5 eV and 853.3 eV are
identified respectively with 3d${}_{z^{2}-r^{2}}$ and 3d${}_{x^{2}-y^{2}}$
symmetry of the pyramidal sites explaining the sign and the hierarchy of the
XLD signal [40, 41]. An in-plane contribution is also expected from the XAS
spectra of the square planar 3d${}_{x^{2}-y^{2}}$ state of the Ni1+ and it is
at the origin of the shoulder observed at 852.3 eV for the $\sigma$
polarization. Nevertheless, the overall shape of the XAS-XLD is dominated by
the Ni2+ in the pyramidal site with only a smaller contribution of the Ni1+
site in the square plane. These multi-valent features of the Sm9Ni9O22
structure are made clearer in the RIXS maps. The data for the $\sigma$ and
$\pi$ polarization, along with their linear dichroism are displayed in Fig.4
d-f. The RIXS data show a similar dichroism as observed in XAS, e.g. with a
larger contribution in the $\pi$ polarization at $\approx$ 852.5 eV. The
calculated RIXS maps obtained for a Ni2+ within a pyramid is in a fair
agreement with the experimental maps (Fig.4g-i). A noticeable difference on
comparing Fig4f and i is the experimental intensity at an energy loss of 1.7
eV obtained for an excitation energy of $\approx$ 852.3 eV, that is not
present in the calculation done for Ni2+. This excitation energy corresponds
to Ni1+ and its higher intensity for a $\sigma$ polarization is consistent
with a square planar symmetry (see calculation in figure S15), in accordance
with 3d${}_{x^{2}-y^{2}}$ $\longrightarrow$ 3dxy, 3dxz/3dyz d-d transitions as
also discussed by Rossi et al. [42] in the case of an IL-nickelate. The RIXS
signatures have also been investigated at RT for some excitation energies. The
dichroic nature of Sm9Ni9O22 is rather similar at room temperature, as shown
in (Supplementary Information, Fig. S7 and S8) but some d-d transitions
experience blue shift of $\approx$ 0.2 eV indicating electronic modification
between low to room temperature.

Figure 5: Temperature dependent spectroscopic studies at the O-K edge by XAS
and RIXS maps (in energy loss) of Sm9Ni9O22. (a-c) Measurements at T = 25 K,
(a) O-K edge PFY XAS spectra in $\sigma$ and $\pi$ incident polarizations, (b
& c) RIXS map at O-K edge for (b) $\sigma$ incident photon polarization and
(c) $\pi$ incident photon polarization. (d-f) Measurements at T = 300 K, (d)
O-K edge PFY XAS spectra in $\sigma$ and $\pi$ incident polarizations, (e & f)
RIXS map at O-K edge for (e) $\sigma$ incident photon polarization and (f)
$\pi$ incident photon polarization.

In complement, similar spectroscopic studies have been carried out at the O-K
edge for the parent and reduced phase, both at low and room temperature. O-K
edge XAS has been previously used to investigate the electronic transition in
a perovskite nickelate [43]. Here in the parent phase, along with the XAS, the
electronic transition is further characterized by the RIXS maps that also
indicates a band gap opening at low temperature (Supplementary Information,
Figure S9 and S10).

The O-K edge XAS and RIXS maps of Sm9Ni9O22 shown in Fig.5a-f and in
Supplementary Information, Fig.S11 & S12 are very different from the parent
SmNiO3. The pre-peak of the O-K is very weak which also confirms a reduced
valency and co-valency as compared to the Ni3+ of SmNiO3. The position of the
pre-peak also confirms a band gap opening of $\approx$ 3 eV for the reduced
Sm9Ni9O22. While being of small intensity, the pre-peak shows a strong
incident photon polarization dependence. At low temperature, this pre-peak
intensity is successively stronger in energy for a $\sigma$, a $\pi$ and then
a $\sigma$ polarization that could be interpreted by a combination of
hybridization of the O-2p orbital with the unoccupied 3d${}_{x^{2}-y^{2}}$ of
the Ni1+ and then with the 3d${}_{z^{2}}$ and 3d${}_{x^{2}-y^{2}}$ of the
Ni2+. This goes hand in hand with the respective energy positions of the
calculated unoccupied density of states of these orbitals as shown in
(Supplementary Information, Fig.S16).

The XAS O-K pre-edge at room temperature seems smeared out, and exhibits a
weaker dichroism. While observing the RIXS maps, a strong fluorescence line
appears at room temperature causing the blurring of the O-K XAS pre-edge
feature. This fluorescence line occurs for an incident photon energy
corresponding to Ni-O hybridization and with the smallest energy loss,
involving energy bands at the upper part of the valence band. A disappearance
with temperature of this fluorescence line reveals an electronic phase
transition and somehow a charge localization at low temperature for the
Sm9Ni9O22. Due to the weak dichroic character of this band, the localisation
might concern the Ni3+ in the octahedra whose electron configuration in the
parent phase is also prone to charge localization [44]. Futhermore, additional
electrons associated to oxygen deficient perovskite nickelates have also been
predicted to strongly localize on octahedral sites, explaining the dramatic
decrease in the electrical conductivity by several orders of magnitude in them
[27].

## Discussion

Here, the topotactic reduction of a parent perovskite SmNiO3 thin-film grown
on NdGaO3 results in a homogeneous thin-film of Sm9Ni9O22 phase, with a
tetragonal symmetry and a reduced o-o-p parameter of 3.63 Å. This structure
emerges from families of {303}${}_{\text{pc}}$ apical oxygen vacancy (Vo)
ordering giving a formula A9B9O22 as revealed by 4D-STEM. A single family
plane will give A3B3O8 phase with only octahedral and pyramidal Ni sites. Here
we show that the A9B9O22 phase emerges from Vo families along two
orientations, giving rise to square-planar, pyramidal and octahedral Ni sites.
Square-planar environment occurs at the intersection of orthogonal families
and harbours a Ni1+ indicating a possible route to IL-structure. It is to note
that such (303)${}_{\text{pc}}$ phases have also been reported as structural
defects, constrained in quasi-2D nano-sheet on the top of IL-NNO2 [24]. In
fact, such (303)${}_{\text{pc}}$ defect might be present as nanodomains,
presumably along a single orientation, or overlap along different orientation
and order at longer range as we showed here, enabling phase engineering in the
ABO2-ABO3 system. In such case, this A9B9O22 phase could be a member of a
potentially new family of nickel oxides, that forms by {n0n}${}_{\text{pc}}$
ordered Vo. Overlapping (n0n)${}_{\text{pc}}$ and (0nn)${}_{\text{pc}}$ apical
oxygen vacancy planes results in a general formula
A${}_{n^{2}}$B${}_{n^{2}}$O${}_{3n^{2}-2n+1}$ that is indeed spanning from the
IL-ABO2 with only square planar to the perovskite ABO3 with only octahedral Ni
sites. In-between, octahedral, square-planar and pyramidal sites are in the
ratio $(n-2)^{2}:2:(4n-6)$. Such {n0n}${}_{\text{pc}}$ planes could be seen as
an analogous case of the Ruddlesden-Popper (RP) faults that is also present
from the individual 2D defects and scale to a fully developed RP phase.

The stability of this family is still to be understood and being an
intermediate phase, it will depend on epitaxial strain and redox condition.
This is in line with our observations indicating that the orthorhombic
substrate plays a role, as strain in thin films is generally known to control
the creation of Vo, even for the stabilization of superconducting IL-nickelate
films [45]. Additionally, the choice of the reducing agent such as CaH2 [8],
alumina [11] and atomic hydrogen [45] can influence the selective
stabilization of each phase. The A9B9O22 phase is also a neighbour of other
members of the oxygen deficient nickelate series, such as the brownmillerite
ABO2.5 [17] and the A3B3O7 [26]. The chemical dynamics behind the selective
stabilization of these phases is still not understood. Possibly, this new
phase A9B9O22 might have been present along the samples that reported the
A3B3O7 [35, 36], or it might have been misinterpreted as that.

Interestingly, ab-initio calculations and spectroscopic studies coherently
indicate that the Sm9Ni9O22 phase consists of tri-component Ni sites, with a
varying valence from Ni1+ to Ni3+. Each component is associated with a unique
symmetry, charge, orbital, and spin distribution, that collectively express a
multivalent character with different degrees of orbital polarisation at each
site, as evidence by our spectroscopic studies. The tri-components are
distributed in the lattice with a periodicity determined by the Vo ordering,
contributing to a commensurate charge order for this system with a periodicity
($\frac{1}{3}$, 0, $\frac{1}{3}$) r.l.u. This periodicity is strikingly
similar with that of the CO reported in IL-nickelate thin films [22, 21, 20,
23]. This possibly indicates the existence of this phase in IL-nickelate thin-
films as a defect, and this is further strengthened from recent studies
reporting apical oxygen ordering in IL-nickelate thin-films that exhibited
similar CO [24, 25]. It questions how intrinsic is the ($\frac{1}{3}$, 0)
r.l.u. CO in IL-NdNiO2.

The parent perovskite SmNiO3 exhibits a metal-to-insulator (MIT) around 375K
as evidenced by our transport measurements, and this is similar to what was
previously reported [28] in these. However the novel system Sm9Ni9O22 has a
complex multivalent electronic landscape, with an intriguing carrier
localization contributing to an enhanced resistivity. This is in contrast with
the case of oxide semiconductors, where Vo cause an enhanced conductivity.
However, this is in line with previous reports on perovskite nickelates, where
the Vo cause localization of electrons at the octahedral Ni sites [27], with
the subsequent filling of the ligand hole resulting in a 3$d^{8}$
configuration. The oxygen deficient nickelate discussed in that report might
correspond to this novel phase. Furthermore, RIXS evidence an additional
charge localization at low temperature, presumably associated to the charge
configuration at octahedral Ni site. It is worth noting that the Ni sites at
the parent phase octahedra have known charge fluctuation 3$d^{8}\underline{L}$
$\longleftrightarrow$ 3$d^{8}$ \+ 3$d^{8}\underline{L}^{2}$ playing a role in
the temperature driven MIT. This might be hindered for an isolated octahedron,
giving rise to different charge fluctuation (e.g. 3$d^{8}\underline{L}$
$\longleftrightarrow$ 3$d^{7}$), making the phase diagram of this system more
complex.

On a broader perspective, this is another example of the idea of selective
oxygen incorporation driving stabilization of phases/interfaces that have
metallic, semiconducting, insulating, and even superconducting (with doping)
properties [8, 46, 47, 48]. The carrier distribution in the lattice, which is
associated to the distribution of oxygen atoms, has an identical effect as
that of hydrogenated perovskite nickelate devices proposed for neuromorphic
computing applications [15]. This was demonstrated for an oxygen deficient
NdNiO3-δ nickelate device, that could potentially have this A9B9O22 phase,
where multiple resistance states were achieved by adjusting the electronic
bias [27]. Such findings uncover the possible multitude of applications of
this phase on par with other strongly correlated electron systems.

## Acknowledgements

A.G., A.N., and J.-P.R. acknowledges financing from EU-H2020 IMPRESS under
grant agreement (n. 101094299). A.R. acknowledges financing from LABEX
NanoSaclay and EU-H2020 ESTEEM-3 under grant agreement (n. 823717) for the
doctoral funding. A.G. and X.L. also acknowledges financing from EU-H2020
ESTEEM-3 under grant agreement (n. 823717). Nion UltraSTEM at LPS Orsay and
the FIB at C2N, University of Paris-Saclay was accessed in the TEMPOS project
framework (ANR 10-EQPX-0050). Z.D., B.F. and D.L. acknowledge the funding
support from the Natural Science Foundation of China (No. 12174325) and Guang
Dong Basic and Applied Basic Research Foundation (2023A1515011352). The work
performed in Hong Kong is supported by the Research Grants Council of Hong
Kong through the Early Career Scheme (CityU 21301221) and the General Research
Fund (CityU 11309622). L.V., C.D. and M.H. acknowledge support by the Swiss
National Science Foundation - division II (grant no. 200020_179155 and
200020_207338), and by the European Research Council under the European Union
Seventh Framework Programme (FP7/2007–2013)/ERC Grant Agreement no. 319286
(Q-MAC). We acknowledge SOLEIL Synchrotron for provision of beamtime under the
proposals 20220630 and 20221574. We thank Daniele Preziosi, Alberto Zobelli,
Zhizhong Li, Nicolas Jaouen and Jean-Marc Triscone for fruitful discussions.

## Methods

### Samples

Before growth, the (110)${}_{\text{o}}$-oriented NdGaO3 single-crystalline
substrates were annealed at 1150 ∘C for 10 hours in flowing oxygen to obtain
sharp step-and-terrace surface morphology. Then the pristine SmNiO3 epitaxial
thin films were grown by off-axis radio-frequency magnetron sputtering, at a
substrate temperature of 460 ∘C and a total pressure of 0.18 millibars with an
oxygen/argon mixture of 1:3 [49]. After growth, the samples were cut into
pieces with the size of 5$\times$2.5 mm2 and loosely wrapped in the aluminum
foil. The wrapped samples were then sealed under vacuum in Pyrex glass tubes,
embedded with CaH2 powder of $\approx$ 0.1 g. The sealed glass tubes were
heated in a tube furnace to 230-240 ∘C for 1h and cooled down to room
temperature at a rate of 10 ∘C/min.

### High-resolution HAADF-STEM Imaging and 4D-STEM

Cross sectional transition electron microscopy (TEM) lamellae were prepared
using a focused ion beam (FIB) technique (at C2N, University of Paris-Saclay,
France). Before FIB lamellae preparation, around 20 nm of amorphous carbon was
deposited on top for protection. The HAADF imaging and 4D-STEM was carried out
in a NION UltraSTEM 200 C3/C5-corrected scanning transmission electron
microscope (STEM). The experiments were done at 200 keV with a probe current
of approximately 10 pA and convergence semi-angles of 30 mrad. A MerlinEM
(Quantum Detectors Ltd) in a 4×1 configuration (1024 × 256) has been installed
on a Gatan ENFINA spectrometer mounted on the microscope [50]. The EELS
spectra are obtained using the full 4 × 1 configuration and the 4D-STEM by
selecting only one of the chips (256 × 256 pixels). For 4D-STEM, the EELS
spectrometer is set into non-energy dispersive trajectories and we have used
6-bit detector mode that gives a diffraction pattern with a good signal to
noise ratio without compromising much on the scanning speed.

### Ab-initio Calculations

The first principles calculation were performed using the density functional
theory [51] as implemented in the Quantum ESPRESSO package [52, 53]. The
exchange-correlation functional was approximated by the generalized gradient
approach [54]. Planewave cutoffs of 40 and 400 Ry were used basis-set and
charge density expansions, respectively. Structural relaxations were done
without considering a Hubbard U parameter, while the electronic structure were
calculated with it. All DFT+U calculations were performed with a J = 1.2 eV,
making Ueffective = U - J. All structural relaxations begin with creating a
supercell from a cubic SmNiO3, with $a=3.80$ Å. Then, the oxygen atoms along
the corresponding family of Vo planes were removed. For the Sm9Ni9O22 phase,
after removing the two orthogonal families of apical Vo, it results in a
960-atom $6\times 6\times 6$ supercell. The resulting pseudo-tetragonal o-o-p
parameter of $c=3.63$ Å was obtained by minimizing the total force in a 6x6x6
supercells (the in-plane pseudo-tetragonal parameter $a=b=3.90$ Å, was
frozen). This tetragonal supercell was then reduced using symmetry to an
monoclinic 80-atom cell. Thus obtained structure has the space group $C2/m$
and lattice parameters of $a=9.02$, $b=16.09$, $c=6.47$ Å and
$\beta=87.57^{\circ}$. The structural parameters with atomic positions are
given in Table S1. The lattice vectors of the tetragonal 960-atom phase
($\vec{a_{T}},\ \vec{b_{T}},\ \vec{c_{T}}$) expressed in the symmetrized
monoclinic parameters ($\vec{a_{m}},\ \vec{b_{m}},\ \vec{c_{m}}$) is:
$\vec{a_{T}}=\vec{a_{m}}+\vec{b_{m}}+2\vec{c_{m}},\
\vec{b_{T}}=\vec{a_{m}}-\vec{b_{m}}+2\vec{c_{m}},\
\vec{c_{T}}=2\vec{a_{m}}-2\vec{c_{m}}$.

### XAS and RIXS Measurements

The x-ray absorption spectroscopy (XAS) and resonant inelastic x-ray
scattering (RIXS) measurements of the Ni-L edges and the O-K edge were
performed at the AERHA instrument of the SEXTANTS beamline at the SOLEIL
synchrotron. A detailed optical scheme of the spectrometer that involves
dispersing the x-ray photons as a function of their energy onto a two-
dimensional detector is given in [55], and the experimental geometry in
Fig.4a. The measurements were performed with the a and c sample axes in the
horizontal scattering plane. Measurements at the Ni-L and O-K edges were taken
with $\theta_{i}$=20∘ (where $\theta_{i}$ is the angle between the incoming
x-ray and the sample surface), and linear dichroism was obtained by switching
the incident x-ray polarization between $\sigma$ and $\pi$. RIXS was measured
at 2$\Theta$=85∘, with an overall full-width at half-maximum (FWHM) energy
resolution of 170$\pm$10 meV for the O-K edge and 300$\pm$10 meV for the Ni-L
edge. The measurements were done at T = 25K and T = 300K. The two-dimensional
RIXS data are binned in the isoenergetic direction to form spectra and the
pixel to energy conversion is performed by measuring the position of the
elastic line of the spectrometer while changing the beamline energy.

### HAXPES Measurements

The measurements were carried out at the GALAXIES beamline at the SOLEIL
synchrotron [56] on the HAXPES endstation [57] using a photon energy of 3000
eV, with an incidence angle of 10∘. The bulk sensitivity is defined from the
SESSA simulations [58], that gives a probing depth of around 10 nm for 10∘
incidence. About 95% of the detected signal will be from the elements within
these estimated probing depth. The synchrotron operated with a ring current of
450 mA, giving an intensity of 3.4 × 1013 photons/s at 3000 eV, which was then
reduced using a built-in filter to 5% of the original intensity. The
photoelectrons were detected using a SCIENTA Omicron EW4000 HAXPES
hemispherical analyzer, and a Shirley background [59] was removed prior to
fitting the core levels spectra.

## References

  *   * [1] Dagotto, E.  Complexity in strongly correlated electronic systems.  _Science_ 309, 257–262 (2005).
  * [2] Ahn, C. _et al._ Designing and controlling the properties of transition metal oxide quantum materials.  _Nature materials_ 20, 1462–1468 (2021).
  * [3] Cheong, S.-W.  The exciting world of orbitals.  _Nature Materials_ 6, 927–928 (2007).
  * [4] Keimer, B., Kivelson, S. A., Norman, M. R., Uchida, S. & Zaanen, J.  From quantum matter to high-temperature superconductivity in copper oxides.  _Nature_ 518, 179–186 (2015).
  * [5] Hayward, M.  Topochemical reactions of layered transition-metal oxides.  _Semiconductor Science and Technology_ 29, 064010 (2014).
  * [6] Sanjaya Ranmohotti, K., Josepha, E., Choi, J., Zhang, J. & Wiley, J. B.  Topochemical manipulation of perovskites: Low-temperature reaction strategies for directing structure and properties.  _Advanced Materials_ 23, 442–460 (2011).
  * [7] Wu, Y., Li, D., Wu, C.-L., Hwang, H. Y. & Cui, Y.  Electrostatic gating and intercalation in 2d materials.  _Nature Reviews Materials_ 8, 41–53 (2023).
  * [8] Li, D. _et al._ Superconductivity in an infinite-layer nickelate.  _Nature_ 572, 624–627 (2019).
  * [9] Zeng, S. _et al._ Phase diagram and superconducting dome of infinite-layer Nd1-xSrxNiO2 thin films.  _Physical Review Letters_ 125, 147003 (2020).
  * [10] Pan, G. A. _et al._ Superconductivity in a quintuple-layer square-planar nickelate.  _Nature materials_ 21, 160–164 (2022).
  * [11] Wei, W., Vu, D., Zhang, Z., Walker, F. J. & Ahn, C. H.  Superconducting Nd1-xEuxNiO2 thin films using in situ synthesis.  _Science Advances_ 9, eadh3327 (2023).
  * [12] Kumah, D. P. _et al._ Tuning the structure of nickelates to achieve two-dimensional electron conduction.  _Adv. Mater_ 26, 1935–1940 (2014).
  * [13] Zhang, Z. _et al._ Perovskite nickelates as electric-field sensors in salt water.  _Nature_ 553, 68–72 (2018).
  * [14] Mercy, A., Bieder, J., Íñiguez, J. & Ghosez, P.  Structurally triggered metal-insulator transition in rare-earth nickelates.  _Nature communications_ 8, 1677 (2017).
  * [15] Zhang, H.-T. _et al._ Reconfigurable perovskite nickelate electronics for artificial intelligence.  _Science_ 375, 533–539 (2022).
  * [16] Lee, K. _et al._ Aspects of the synthesis of thin film superconducting infinite-layer nickelates.  _Apl Materials_ 8, 041107 (2020).
  * [17] Wu, Y.-M. _et al._ Topotactically induced oxygen vacancy order in nickelate single crystals.  _Physical Review Materials_ 7, 053609 (2023).
  * [18] Osada, M., Fujiwara, K., Nojima, T. & Tsukazaki, A.  Improvement of superconducting properties in La1-xSrxNiO2 thin films by tuning topochemical reduction temperature.  _Physical Review Materials_ 7, L051801 (2023).
  * [19] Subedi, A.  Possible structural quantum criticality tuned by rare-earth ion substitution in infinite-layer nickelates.  _Physical Review Materials_ 7, 024801 (2023).
  * [20] Krieger, G. _et al._ Charge and spin order dichotomy in NdNiO2 driven by the capping layer.  _Physical Review Letters_ 129, 027002 (2022).
  * [21] Ren, X. _et al._ Symmetry of charge order in infinite-layer nickelates.  _arXiv preprint arXiv:2303.02865_ (2023).
  * [22] Tam, C. C. _et al._ Charge density waves in infinite-layer NdNiO2 nickelates.  _Nature Materials_ 21, 1116–1120 (2022).
  * [23] Rossi, M. _et al._ A broken translational symmetry state in an infinite-layer nickelate.  _Nature Physics_ 18, 869–873 (2022).
  * [24] Raji, A. _et al._ Charge distribution across capped and uncapped infinite-layer neodymium nickelate thin films.  _(Accepted, In production in Small doi:10.1002/smll.202304872), arXiv preprint arXiv:2306.10507_ (2023).
  * [25] Parzyck, C. _et al._ Absence of $3a_{0}$ charge density wave order in the infinite layer nickelates.  _arXiv preprint arXiv:2307.06486_ (2023).
  * [26] Moriga, T. _et al._ Reduction of the perovskite-type LnNiO3 (Ln=Pr, Nd) to Ln3Ni3O7 with monovalent nickel ions.  _Solid State Ionics_ 74, 211–217 (1994).
  * [27] Kotiuga, M. _et al._ Carrier localization in perovskite nickelates from oxygen vacancies.  _Proceedings of the National Academy of Sciences_ 116, 21992–21997 (2019).
  * [28] Torriss, B., Margot, J. & Chaker, M.  Metal-insulator transition of strained SmNiO3 thin films: Structural, electrical and optical properties.  _Scientific reports_ 7, 40915 (2017).
  * [29] Hÿtch, M., Snoeck, E. & Kilaas, R.  Quantitative measurement of displacement and strain fields from hrem micrographs.  _Ultramicroscopy_ 74, 131–146 (1998).
  * [30] Yang, C., Wang, Y., Sigle, W. & van Aken, P. A.  Determination of grain-boundary structure and electrostatic characteristics in a SrTiO3 bicrystal by four-dimensional electron microscopy.  _Nano Letters_ 21, 9138–9145 (2021).
  * [31] Jia, Q. & Gloter, A.  Real-space observation of potential reconstruction at metallic/insulating oxide interface.  _Advanced Materials Interfaces_ 10, 2202165 (2023).
  * [32] Kim, Y.-M. _et al._ Probing oxygen vacancy concentration and homogeneity in solid-oxide fuel-cell cathode materials on the subunit-cell level.  _Nature Materials_ 11, 888–894 (2012).
  * [33] Hirai, K. _et al._ Anisotropic in-plane lattice strain relaxation in brownmillerite SrFeO2.5 epitaxial thin films.  _Journal of Applied Physics_ 114 (2013).
  * [34] Parsons, T. G., D’Hondt, H., Hadermann, J. & Hayward, M. A.  Synthesis and structural characterization of La1-xAxMnO2.5 (A= Ba, Sr, Ca) phases: mapping the variants of the brownmillerite structure.  _Chemistry of Materials_ 21, 5527–5538 (2009).
  * [35] Wang, B.-X. _et al._ Synthesis and characterization of bulk Nd1-xSrxNiO2 and Nd1-xSrxNiO3.  _Physical Review Materials_ 4, 084409 (2020).
  * [36] Wei, W. _et al._ Solid state reduction of nickelate thin films.  _Physical Review Materials_ 7, 013802 (2023).
  * [37] Krieger, G. _et al._ Synthesis of infinite-layer nickelates and influence of the capping-layer on magnetotransport.  _Journal of Physics D: Applied Physics_ 56, 024003 (2022).
  * [38] Green, R. J., Haverkort, M. W. & Sawatzky, G. A.  Bond disproportionation and dynamical charge fluctuations in the perovskite rare-earth nickelates.  _Physical Review B_ 94, 195127 (2016).
  * [39] Bisogni, V. _et al._ Ground-state oxygen holes and the metal–insulator transition in the negative charge-transfer rare-earth nickelates.  _Nature Communications_ 7, 13017 (2016).
  * [40] Stavitski, E. & De Groot, F. M.  The CTM4XAS program for EELS and XAS spectral shape analysis of transition metal L edges.  _Micron_ 41, 687–694 (2010).
  * [41] Ikeno, H., de Groot, F. M., Stavitski, E. & Tanaka, I.  Multiplet calculations of L2,3 x-ray absorption near-edge structures for 3d transition-metal compounds.  _Journal of Physics: Condensed Matter_ 21, 104208 (2009).
  * [42] Rossi, M. _et al._ Orbital and spin character of doped carriers in infinite-layer nickelates.  _Physical Review B_ 104, L220505 (2021).
  * [43] Meyers, D. _et al._ Pure electronic metal-insulator transition at the interface of complex oxides.  _Scientific Reports_ 6, 1–8 (2016).
  * [44] Mikheev, E. _et al._ Tuning bad metal and non-fermi liquid behavior in a Mott material: Rare-earth nickelate thin films.  _Science Sdvances_ 1, e1500797 (2015).
  * [45] Lee, K. _et al._ Linear-in-temperature resistivity for optimally superconducting (Nd, Sr)NiO2.  _Nature_ 619, 288–292 (2023).
  * [46] Chen, L. _et al._ Roles of point defects in thermal transport in perovskite barium stannate.  _The Journal of Physical Chemistry C_ 122, 11482–11490 (2018).
  * [47] Suntivich, J. _et al._ Design principles for oxygen-reduction activity on perovskite oxide catalysts for fuel cells and metal–air batteries.  _Nature Chemistry_ 3, 546–550 (2011).
  * [48] Liu, Z. _et al._ Dominant role of oxygen vacancies in electrical properties of unannealed LaAlO3/SrTiO3 interfaces.  _Journal of Applied Physics_ 115 (2014).
  * [49] Catalano, S. _et al._ Electronic transitions in strained SmNiO3 thin films.  _APL Materials_ 2 (2014).
  * [50] Tencé, M. _et al._ Electron energy-loss spectroscopy using MerlinEM-Medipix3 detector.  _Microscopy and Microanalysis_ 26, 1940–1942 (2020).
  * [51] Kohn, W. & Sham, L. J.  Self-consistent equations including exchange and correlation effects.  _Physical Review_ 140, A1133 (1965).
  * [52] Giannozzi, P. _et al._ Quantum espresso: a modular and open-source software project for quantum simulations of materials.  _Journal of Physics: Condensed matter_ 21, 395502 (2009).
  * [53] Giannozzi, P. _et al._ Advanced capabilities for materials modelling with Quantum ESPRESSO.  _Journal of Physics: Condensed matter_ 29, 465901 (2017).
  * [54] Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Physical Review Letters_ 77, 3865 (1996).
  * [55] Chiuzbăian, S. G. _et al._ Design and performance of AERHA, a high acceptance high resolution soft x-ray spectrometer.  _Review of Scientific Instruments_ 85, 043108 (2014).
  * [56] Rueff, J.-P. _et al._ The GALAXIES beamline at SOLEIL synchrotron: Inelastic X-ray scattering and photoelectron spectroscopy in the hard X-ray range.  _J. Synchrotron Rad._ 22, 175 (2015).
  * [57] Céolin, D. _et al._ Hard X-ray photoelectron spectroscopy on the GALAXIES beamline at the SOLEIL synchrotron.  _Journal of Electron Spectroscopy and Related Phenomena_ 190, 188–192 (2013).
  * [58] Smekal, W., Werner, W. S. & Powell, C. J.  Simulation of electron spectra for surface analysis (SESSA): a novel software tool for quantitative Auger-electron spectroscopy and X-ray photoelectron spectroscopy.  _Surface and Interface Analysis_ 37, 1059–1067 (2005).
  * [59] Shirley, D. A.  High-resolution x-ray photoemission spectrum of the valence bands of gold.  _Physical Review B_ 5, 4709 (1972).

See pages - of Supplementary.pdf