# Mol-PECO: a deep learning model to predict human olfactory perception from
molecular structures

Author One Author Two
Affiliation  Mengji Zhang
Schoool of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai,
200030, China
Institute of Industrial Science, The University of Tokyo, Meguro-ku, Tokyo,
153-8505, Japan
<EMAIL_ADDRESS>
&Yusuke Hiki
Department of Biosciences and Informatics
Keio University
Kouhoku-ku, Yokohama, 223-8522, Japan
<EMAIL_ADDRESS>
&Akira Funahashi
Department of Biosciences and Informatics
Keio University
Kouhoku-ku, Yokohama, 223-8522, Japan
<EMAIL_ADDRESS>
&Tetsuya J. Kobayashi
Institute of Industrial Science
The University of Tokyo
Meguro-ku, Tokyo, 153-8505, Japan
<EMAIL_ADDRESS>
This work is done as a visiting student in Tokyo University.Correspondence to
Tetsuya J. Kobayashi or Mengji Zhang.

###### Abstract

While visual and auditory information conveyed by wavelength of light and
frequency of sound have been decoded, predicting olfactory information encoded
by the combination of odorants remains challenging due to the unknown and
potentially discontinuous perceptual space of smells and odorants. Herein, we
develop a deep learning model called Mol-PECO (Molecular Representation by
Positional Encoding of Coulomb Matrix) to predict olfactory perception from
molecular structures. Mol-PECO updates the learned atom embedding by
directional graph convolutional networks (GCN), which model the Laplacian
eigenfunctions as positional encoding, and Coulomb matrix, which encodes
atomic coordinates and charges. With a comprehensive dataset of $8,503$
molecules, Mol-PECO directly achieves an area-under-the-receiver-operating-
characteristic (AUROC) of $0.813$ in $118$ odor descriptors, superior to the
machine learning of molecular fingerprints (AUROC of $0.761$) and GCN of
adjacency matrix (AUROC of $0.678$). The learned embeddings by Mol-PECO also
capture a meaningful odor space with global clustering of descriptors and
local retrieval of similar odorants. Our work may promote the understanding
and decoding of the olfactory sense and mechanisms.

##  1 Introduction

Olfaction is one of the essential senses, where the sense of smell is
triggered by the binding of odorant molecules to olfactory receptors and is
shaped by the subsequent neural processing of the received information in the
brain(Sobel et al., 1998; Lapid et al., 2011). Unlike vision and hearing,
however, the prediction of olfactory perception for odorant molecules remains
challenging. On the one hand, some molecules with different functional groups
share identical smells(Sell, 2006). On the other hand, other molecules with
similar structures can produce totally different perceptions(Boesveldt et al.,
2010). The structures and perceptions of the odorant molecules are nonlinearly
and discontinuously related(Keller et al., 2017). Unvieling the quantitative
structure-odor relationship (QSOR) is indispensable for understanding the
coding principle of olfactory information(Keller et al., 2017) and also
essential for predicting and designing smells and flavors for various
applications such as food technologies(Polster and Schieberle, 2017).

Machine learning (ML) is a promising approach to untangle such a complicated
relationship. However, the prediction of olfactory perception from molecular
structures is strongly dependent on molecular representation(Pattanaik and
Coley, 2020). Molecular fingerprints, which encode chemical substructures into
fixed-length vectors, are the major and classical molecular representation,
yet demonstrate limited performances in QSOR due to the inefficient feature
extraction by hand-crafted rules(Rogers and Hahn, 2010; Moriwaki et al.,
2018). To learn a good representation of molecules from data, graph
convolutional networks (GCNs) have been widely applied in molecular
modeling(Kipf and Welling, 2016), e.g., quantum chemistry(Hofstetter et al.,
2022; Yang et al., 2019; Feinberg et al., 2018), biophysics(Wang et al., 2019;
Withnall et al., 2020), and biological side effects(Yang et al., 2019;
Withnall et al., 2020). The conventional GCN models each molecule by the
adjacency matrix, which encodes the chemical bonds as a graph, and performs
the information aggregation among the neighbors prescribed by the adjacency
matrix(Yang et al., 2019; Sanchez-Lengeling et al., 2019). While GCN has been
reported to outperform conventional ML with molecular fingerprints in standard
tasks(Sanchez-Lengeling et al., 2019), it still has technical drawbacks, which
potentially limit the applicability of GCN and the adjacency matrix ((Kreuzer
et al., 2021; Oono and Suzuki, 2019; Topping et al., 2021)) to learning QSOR.
First, the adjacency matrix cannot encode the atomic and global 3D information
of molecules(Pattanaik and Coley, 2020), even though such atomic and 3D
information is the major determinant of binding affinities between odorant
molecules and olfactory receptors(Floriano et al., 2000). Second, graphs do
not have a canonical coordinate representation, which contrasts with sequences
and images for which one-, two-, and three-dimensional lattice coordinates are
canonical. As a result, the GCN employs permutation invariant operations,
e.g., message passing or neighboring aggregation, which then limit its
expressive power to discriminate molecules with different structures((Kreuzer
et al., 2021)) and also induces oversmoothing(Oono and Suzuki, 2019) and
oversquashing(Topping et al., 2021). All of these factors may hamper the
network to efficiently learn the QSOR of odorant molecules with a variety of
structures.

In this work, we develop a deep learning model (Mol-PECO) for QSOR, which aims
at the multi-label classification of olfactory perception from molecular
structure (Figure 1). To address the problems in conventional ML and GCN
approaches, Mol-PECO combines the Coulomb matrix (CM) and Spectral Attention
Network (SAN). CM is a simple global representation of a molecule by Coulombic
forces between atoms in the molecule calculated with the nuclear charges and
the corresponding 3D coordinates(Rupp et al., 2012). CM therefore could encode
more detailed structural information than the adjacency matrix which only
represents the chemical bonds between atoms. SAN is an architecture of graph
attention network (GAT)(Kreuzer et al., 2021), which uses the full Laplacian
spectrum of a molecular graph for a learned positional encoding (LPE).
Eigenfunctions of the graph Laplacian hierarchically describe the global and
local structures in a graph and thereby provides a way to canonically
characterize graphs and to define positional information of nodes
(atoms)(Dwivedi and Bresson, 2020). SAN can also be categorized as an attempt
to extend the expressive power in graph-based architectures by endowing a
canonical coordinate or positional information with graphs by using their
spectral information(Beaini et al., 2021; Kreuzer et al., 2021; Dwivedi and
Bresson, 2020). CM can be naturally combined with SAT by regarding CM as a
weighted adjacency matrix.

Based on the learned representation, Mol-PECO directly predicts $118$ odor
descriptors of perception for each odorant molecule. Mol-PECO achieves area
under the receiver operating characteristic curve (AUROC) of $0.813$ and area
under the precision-recall curve (AUPRC) of $0.181$, whereas the conventional
MLs of molecular fingerprints fail to balance AUROC and AUPRC; the ML method
(cfps-KNN) with the highest AUROC of $0.761$ has low AUPRC of $0.057$ and one
(mordreds-RF) with the highest AUPRC of $0.144$ shows low AUROC of $0.723$.
Thus, Mol-PECO may boost the prediction of QSOR for applications and also
contribute to the understanding of the principle underlying olfactory
information processing.

Figure 1: Overview of Mol-PECO. (a) Typical molecular structures and the
corresponding odor descriptors are shown as examples. (b) The main workflow of
modeling quantitative structure-odor relationship (QSOR). (c) The detailed
model architecture of Mol-PECO and its three features: 1) skip-connection in
graph neural networks to alleviate over-smoothing, 2) fully-connected
molecular representation by Coulomb matrix to suppress over-squashing, and 3)
positional encoding by Laplacian eigenfunctions.

##  2 Dataset:A comprehensive human olfactory perception dataset

In this work, we use the data in which each molecular structure is paired with
multiple odor descriptors (Figure 1 a). The dataset in our work is compiled
from ten expert-labeled sources: Arctander’s dataset ($n=3,102$, (Arctander et
al., 1960)), AromaDb ($n=1,194$, (Kumar et al., 2018)), FlavorDb ($n=525$,
(Garg et al., 2018)), FlavorNet ($n=718$, (Acree, 2004)), Goodscents
($n=6,158$, (Flavor, )), Fragrance Ingredient Glossary ($n=1,135$,
(Association et al., 2003)), Leffingwell’s dataset ($n=3,523$, (Leffingwell,
2001)), Sharma’s dataset ($n=4,006$, (Sharma et al., 2021)), OlfactionBase
($n=5,105$, (Sharma et al., 2022)), and Sigma’s Fragrance and Flavor Catalog
($n=871$, (Corporation, 2011)). These datasets are retrieved from the archive
of https://github.com/pyrfume/pyrfume-data. The data cleaning procedure
includes 1) merging the overlapped molecules, 2) filtering the conflict
descriptors, and 3) filtering the rare descriptors assigned to less than $30$
molecules. After data cleaning, we obtain a comprehensive dataset of $8,503$
molecules and $118$ odor descriptors.

This comprehensive human olfactory perception dataset is multi-labeled, with
every molecule labeled with one or several odor descriptors. For the number of
molecules associated with each odor descriptor, the distribution is
imbalanced: each of $112$ odor descriptors possesses $\leq 800$ molecules
whereas the other $6$ descriptors possess $>800$ molecules (Figure 2 a). For
the number of descriptors associated with each molecule, the distribution is
also skewed, with $8,054$ molecules possessing $\leq 5$ odor descriptors and
$449$ molecules possessing $>5$ odor descriptors (Figure 2 b). For co-
occurrence, descriptors of ’fruity’, ’green’, ’sweet’, ’floral’, and ’woody’
co-occur with almost all the descriptors, while odorless molecules co-occur
with no other molecules (Figure 2 c). The data split is built by second-order
iterative stratification (Szymański and Kajdanowicz, 2017), which aims at
splitting multi-label dataset and preserves the label ratios in each split
with an iterative sampling design. The whole dataset is splitted into train /
validation / test datasets of $6,802$ / $864$ / $837$ pairs, respectively.

Figure 2: The comprehensive human olfactory dataset. (a) Distribution of
molecules across odor descriptors. (b) Distribution of descriptors across
molecules. (c) Co-occurrence matrix of 118 odor descriptors. The heatmap is
demonstrated with logarithm transformation, and the descriptors are ordered
alphabetically.

##  3 Results

In this section, we first introduce Coulomb-GCN, which updates GCN by
replacing the adjacency matrix with CM. After verifying the effectiveness of
Coulomb-GCN, we have Mol-PECO by further replacing the random embedding of
atoms in Coulomb-GCN with positional encoding, in which the spectral
information of the CM is employed to have a structure-aware embedding.

###  3.1 Fully-connected graph by Coulomb matrix is superior to sparse graph
by adjacency matrix

We calculate CM, which models atomic energies with the internuclear Coulomb
repulsion operator (Rupp et al., 2012; Schrier, 2020), and use it as our
molecular representation. In CM, the diagonal entries refer to a polynomial
fit of atomic energies and off-diagonal entries represent the Coulomb
repulsion force between atomic nuclei. Although the adjacency matrix has been
widely used in molecular modeling(Mahmood et al., 2021; Wang et al., 2022), CM
as an emerging molecular representation can have at least two advantages: 1)
CM handles the over-squashing plight by allowing direct paths between distant
nodes in the fully-connected graph representation (Figure 3a); 2) distance by
Frobenius norm between CM and adjacency matrix is $5$–$10$ times smaller than
that between random initialized matrix and adjacency matrix, indicating that
CM is fully-connected while preserving a similarity to adjacency matrix
(Figure 3b).

Figure 3: The motivation and workflow of modeling the Coulomb matrix. (a) An
example of Coulomb matrix as a fully-connected graph for Propargyl alcohol, a
clear colorless liquid with a geranium-like odor. (b) Similarity between
Coulomb matrix and adjacency matrix, indicated by the distances. The distance
is calculated with the Frobenius norm. (c) The workflow of modeling Coulomb
matrix by graph neural networks.

We build a nonlinear map (named Coulomb-GCN) between molecular structures and
human olfactory perception (Figure 3c) by replacing the adjacency matrix in
message passing of GCN with CM. Specifically, starting from random atom
embedding, the learned atom embedding is obtained by message passing on fully-
connected molecular graph with neighbor weights specified by the entries of
CM. The molecular embedding is extracted by sum pooling and fed to a multi-
label classification module to predict $118$ odor descriptors. Considering the
gap between maximal and minimal entries in CM, normalization of entries may
affect the performance. We test Minmax and Frobenius normalizations in a
matrix-wise manner.

We evaluate and compare the prediction accuracy of GCN with adjacency matrix
and those of Coulomb-GCN with the different normalizations of CM (Table 1).
Compared with GCN with adjacency matrix (AUROC of $0.678$), gains in AUROC are
observed in Coulomb-GCN with Frobenius normalization (AUROC of $0.759$) and
minmax normalization (AUROC of $0.713$). Coulomb-GCN with Frobenius
normalization also achieves higher performances in five out of six evaluation
metrics (Table 1): AUROC (improved from $0.678$ to $0.759$), AUPRC (improved
from $0.111$ to $0.143$), specificity (improved from $0.625$ to $0.744$),
precision (improved from $0.079$ to $0.089$), and accuracy (improved from
$0.726$ to $0.780$).

Table 1: Prediction performances of Coulomb matrix with Minmax and Frobenius
normalizations and adjacency matrix in GCN. The highest score of each metric
is shown in bold.

Representation | Normalization1 | AUROC | AUPRC | Precision | Recall | Specificity | Accuracy
---|---|---|---|---|---|---|---
Adjacency matrix | - | 0.678 | 0.111 | 0.079 | 0.827 | 0.625 | 0.726
Coulomb matrix | Minmax | 0.713 | 0.138 | 0.082 | 0.811 | 0.687 | 0.749
Coulomb matrix | Frobenius | 0.759 | 0.143 | 0.089 | 0.816 | 0.744 | 0.780

  * 1

Normalization refers to the normalization methods for Coulomb matrix.

###  3.2 Directional graph modeling by Laplacian eigenfunctions improves
prediction accuracy

The graph Laplacian and its spectral information enable us to characterize the
global and substructures of graphs(Mohar et al., 1991; Chung, 1997; Spielman,
2012). Specifically, the graph Laplacian is defined as $L=D-A$, where $D$ and
$A$ refer to the degree and adjacency matrices. $L$ is positive semi-definite
with one trivial and the other nontrivial eigenvalues. In this work, the
Laplacian defined by CM acts as an extension of the normal Laplacian
($L=D-W$), where $W$ refers to the weighted matrix (CM in this work) and
possesses the same properties as the graph Laplacian (e.g., symmetric and
positive semi-definite). In particular, the eigenvectors of $L$ provide an
optimal solution to the Laplacian quadratic form
($f^{T}Lf=1/2\sum_{i,j}X(i,j){(f_{i}-f_{j})}^{2}$)(Mohar et al., 1991;
Spielman, 2012; Chung, 1997), encoding the geometric information of graphs.

Figure 4: The motivation and architecture of Mol-PECO. (a) Structural
information carried by the Laplace spectrum of the Coulomb matrix. Low-
frequency eigenvectors, calculated with graph Laplacian, as the input matrix
for positional encoding and $3$ examples, including cyclic and acyclic
molecules, of eigenvalue $\lambda_{i}$ and eigenvector $\phi_{i}$ for
molecular graphs ($i\in\\{1,2,3\\}$). Color indicates the value of each
component (node) of the eigenvectors. (b) The architecture of Mol-PECO. Mol-
PECO learns the positional encoding (LPE) with Transformer of graph Laplacian,
and updates the atom embedding with GCN of Coulomb matrix and LPE.
Specifically, GCN is implemented with skip-connection to relieve over-
smoothing. Coulomb matrix, the fully connected graph representation,
suppresses over-squashing with direct connections between nodes. With the
updated atom embedding, Mol-PECO extracts the molecular embedding with sum
pooling, and predicts 118 odor descriptors with neural networks of molecular
embedding. In this work, $p$ and $d$ is set as $20$ and $32$, respectively.

Given these properties of the Laplacian graph, we use the Laplacian
eigenfunctions of the CM to encode the positional information of molecular
graphs. Typical results of a cyclic odorant (5-pentyloxolan-2-one, flowing
from left to right in $\lambda_{1}$) and acyclic odorants (hexyl
3-methylbutanoate, flowing from left to right in $\lambda_{1}$, and heptyl
pentanoate, flowing from right to left in $\lambda_{1}$) demonstrate the
information carried by low frequency eigenfunctions (Figure 4 a). Combining it
with the Coulomb-GCN, we construct the deep learning framework, named Mol-PECO
(Figure 4b), with the fully-connected molecular representation by CM and the
positional encoding by Laplacian. We choose LPE by Transformer(Kreuzer et al.,
2021) to build the atom embedding. Specifically, LPE concatenates the $p$
lowest eigenvalues and the corresponding eigenvectors as the input matrix
$\Lambda\in R^{p\times 2}$, and learns the encoding with Transformer for every
atom(Kreuzer et al., 2021). We obtain AUROC of 0.796 and AUPRC of 0.153 with
LPE of raw CM. We further perform the experiments for LPE of asymmetric
normalized CM and obtain additional gain of performances by 0.017 and 0.028
for AUROC and AUPRC, respectively.

We compare Mol-PECO with the baseline models (Table 2): the conventional GCN
of graph representations, including the adjacency matrix and the CM, and the
classifiers of fingerprint representations, including Mordreds features
(mordreds)(Moriwaki et al., 2018), bit-based fingerprints (bfps)(Rogers and
Hahn, 2010), and count-based fingerprints (cfps)(Rogers and Hahn, 2010).
Conventional classifiers include $k$-Nearest Neighbor (KNN), random forest
(RF), and gradient boosting (GB). In the fingerprint methods, we first handle
the problem of imbalanced label distribution with Synthetic Minority Over-
sampling Technique (SMOTE)(Chawla et al., 2002), and then perform the
classification. Mol-PECO outperforms the baselines in three out of six
evaluation metrics (Table 2), with AUROC improved from $0.761$ (cfps-KNN) to
$0.813$, AUPRC improved from $0.144$ (mordreds-RF) to $0.181$, and accuracy
improved from $0.780$ (Coulomb-GCN) to $0.808$. Notably, Mol-PECO can balance
AUROC ($0.813$) and AUPRC ($0.181$) whereas the ML method (cfps-KNN) with the
highest AUROC of $0.761$ has a low AUPRC of $0.057$ and one (cfps-RF) with the
highest AUPRC of $0.144$ shows low AUROC of $0.723$. Thus, Mol-PECO boosts the
predictability of QSOR.

Table 2: Performance comparison of Mol-PECO with baseline methods, GCN, and
Coulomb-GCN. The highest scores are shown in bold. The runner-ups are shown
with underlines.

Baseline1 | AUROC2 | AUPRC2 | Precision2 | Recall2 | Specificity2 | Accuracy2
---|---|---|---|---|---|---
cfps-KNN | 0.761 | 0.057 | 0.065 | 0.760 | 0.764 | 0.762
bfps-KNN | 0.758 | 0.055 | 0.064 | 0.748 | 0.769 | 0.759
mordreds-KNN | 0.729 | 0.052 | 0.062 | 0.676 | 0.783 | 0.730
mordreds-RF | 0.723 | 0.144 | 0.241 | 0.483 | 0.964 | 0.723
cfps-RF | 0.689 | 0.137 | 0.258 | 0.418 | 0.961 | 0.690
bfps-RF | 0.671 | 0.119 | 0.227 | 0.381 | 0.962 | 0.672
mordreds-GB | 0.725 | 0.126 | 0.220 | 0.499 | 0.951 | 0.725
cfps-GB | 0.701 | 0.120 | 0.210 | 0.453 | 0.950 | 0.702
bfps-GB | 0.687 | 0.111 | 0.193 | 0.428 | 0.948 | 0.688
adjacency-GCN | 0.678 | 0.111 | 0.079 | 0.827 | 0.625 | 0.726
Coulomb-GCN | 0.759 | 0.143 | 0.089 | 0.816 | 0.744 | 0.780
Mol-PECO-sym | 0.796 | 0.153 | 0.088 | 0.817 | 0.787 | 0.802
Mol-PECO-asym | 0.813 | 0.181 | 0.104 | 0.819 | 0.797 | 0.808

  * 1

Baseline includes conventional classifiers of fingerprint representations and
graph convolutional networks (GCN) of molecular graphs. Fingerprint
representations include count-based fingerprints (cfps), bit-based
fingerprints (bfps), and Mordreds features (mordreds). The conventional
classifiers include $k$-Nearest Neighbor (KNN), random forest (RF), and
gradient boosting (GB). Molecular graph representations include adjacency
matrix (adjacency-GCN) and Coulomb matrix (Coulomb-GCN). Mol-PECO-sym refers
to the performances with LPE of raw Coulomb matrix. Mol-PECO-asym refers to
the performances with LPE of asymmetrically normalized Coulomb matrix.

  * 2

The evaluation metrics are calculated in testing set.

###  3.3 Learned odor space by Mol-PECO

To investigate the learned structure of multiple odors in relation with
descriptors, we perform dimensionality reduction over the output of Mol-PECO’s
penultimate layer to build the latent odor space and evaluate it at global and
local scales. At the global scale, we inspects how appropriately the clusters
of odors in this latent space represent the information of descriptors, while,
at the local scale, we examines whether individual molecules can possess a set
of odor descriptors similar to those of nearby molecules or not.

Figure 5: The odor space built from Mol-PECO and its global and local
properties. Global view of the learned odor space with dimensionality
reduction by t-SNE is shown in (a), (b), and (c) where each red or gray dot
represents one odor molecule with and without the corresponding descriptor,
respectively. The value of AUROC in each panel is the AUROC of the
corresponding descriptor. (a) The descriptors with most molecules (‘odorless’,
‘fruity’, ‘green’, ‘woody’, ‘sweet’, and ‘floral’), (b) those for the
sulfurous compounds (‘alliaceous’, ‘garlic’, and ‘onion’), and (c) those with
similar semantic meaning (‘oily’ and ‘fatty’). (d) Local view of learned odor
space investigated with nearest neighbor retrieval of reference odorless and
odorant molecules. AUROC in (a), (b), and (c) refers to the AUROC of single
odor descriptor. AUROC in (d) refers to the unweighted AUROC of 118 odor
descriptors.

For global structure, we analyze the distribution of the high-frequency
descriptors, the structure-correlated descriptors, and the synonymy
descriptors. For the high-frequency descriptors, we evaluate ‘odorless’,
‘fruity’, ‘green’, ‘woody’, ‘sweet’, and ‘floral’, which are the top six most
assigned descriptors. The odorless molecules show a high AUROC value (AUROC =
0.94) and cluster far away from other molecules (Figure 5a, top left panel),
verifying Mol-PECO’s ability to distinguish odorant and odorless molecules.
Descriptor ‘woody’ localizes on the clusters at the middle of the odor space,
whereas descriptors of ‘fruity’, ‘floral’, and ‘green’ are distributed across
multiple clusters over the learned odor space (Figure 5a), indicating that
’woody’ is a more different characteristic than ‘fruity’, ‘floral’, and
‘green’. As shown in Figure 2c, these three descriptors tend to co-occur. In
the odor space, we can find clusters for different pairs of them, meaning that
Mol-PECO can disentangle the differences between ‘fruity’-‘floral’,
‘floral’-‘green’, and ‘fruity’-‘green’. The descriptor ‘sweet’ achieves a low
AUROC score (AUROC = $0.67$) and the associated molecules spread over the
space without a specific pattern, presumably reflecting its polysemous nature
and strong association with taste.

For structure-correlated descriptors , molecules associated with ‘alliaceous’,
‘garlic’, or ‘onion’ appear in the cluster of ‘sulfurous’ compounds (Figure
5b) and show a very high score (AUROC = $0.97-0.98$), in accordance with
previous studies about sulfur compounds’ olfactory descriptions(Block, 1992;
Miękus et al., 2020). For synonymy descriptors, molecules of ‘oily’ and
‘fatty’ are clustered into neighbors, validating their similar semantics
(Figure 5c).

For local structure, we investigate one odorless molecule (triphenylphosphane)
and one odorant molecule (1-(1-sulfanylpropylsulfanyl)propane-1-thiol with
descriptors of ‘alliaceous’, ‘fruity’, ‘garlic’, ‘green’, ‘onion’, and
‘sulfurous’) as the examples (Figure 5d). We compare the top-5 nearest
molecules searched by Mol-PECO’s embedding and bfps. The top-5 molecules of
Mol-PECO and bfps are calculated with cosine similarity and Tanimoto
similarity, respectively. For the odorless molecule, Mol-PECO retrieves all
neighbors with the odorless descriptor, but bfps retrieves no molecule.
Notably, the molecules fetched by Mol-PECO possess different substructures
compared with the reference (e.g., all with C-Cl bond and top-2 / 3 / 5 with
carbonyl functional group). For the odorant molecule, Mol-PECO retrieves all
neighbors with shared descriptor, and bfps retrieves four. Moreover, all the
molecules retrieved by bfps are open-chain structured, the same with the
reference. In contrast, Mol-PECO retrieves quite different structured
molecules, four of which are cyclic molecules. Both examples would indicate
Mol-PECO’s promising potential in decoding molecules with different structures
but identical smells.

##  4 Discussion

In this work, we develop Mol-PECO for predicting human olfactory perception
from molecular structures. We handle this QSOR problem by improving the graph-
neural-network-based approach from two aspects: the molecular representation
and the graph modeling method. For molecular representation, we use CM, which
is fully-connected and contains 3D conformer information (3D coordinates and
charges of atoms), instead of adjacency matrix. For graph modeling, we use the
positional encoding from the eigenfunctions of Laplacian to make up for GCN’s
deficiency in directional modeling. Mol-PECO improves the olfactory perception
prediction in two stages: starting from adjacency-GCN (AUROC of $0.678$) to
Coulomb-GCN (AUROC of $0.759$) by introduction of CM, following by from
Coulomb-GCN (AUROC of $0.759$) to Mol-PECO with eigenfunctions of the
symmetric (AUROC of $0.796$) and asymmetric Laplacian (AUROC of $0.813$).

These results indicate that improvement of the expressive power of GCN can
greatly contribute to learning of non-trivial relationships between molecular
structures and olfactory descriptors (labels) in the QSOR problem. Although
extensions of NN architectures are typically accompanied by an increased cost
of learning, the limited size of odorant molecules up to $400$ daltons in
molecular weight enables the practical application of such extended
architectures. Thus, we believe that the QSOR problem provides a good real-
world task to test and demonstrate the effectiveness of advanced
architectures. In particular, pre-training of graph representation using
unlabeled data may lead to further improvements in accuracy of QSOR
prediction(Li et al., 2020; Hu et al., 2019). In addition, as we found by
experiments, employment of asymmetric Laplacians might also contribute to
further technical advancements.

However, the QSOR problem suffers from ambiguity in the labels (descriptors)
assigned to each molecule and also from low objectivity and consistency of the
labeling conducted by few specialists(Trimmer et al., 2019; Keller et al.,
2007). In particular, individuals can have a quite different sense to highly
ambiguous descriptors such as ‘sweet’(Keller et al., 2007). Thus, consistent
labeling by individuals requires a certain amount of training, resulting in
the difficulty to increase the amount of data. Moreover, for the prediction of
mixed odors, which is important for applications, data acquisition is
prohibitive due to the huge number of possible combinations(Meister, 2015). In
the future, it will be important to apply the approach developed in this paper
to the prediction of more objective and systematically measurable outputs,
such as the chemical properties of odorants(Pannunzi and Nowotny, 2019), the
response of olfactory receptors(Bhandawat et al., 2005; Mainland et al.,
2015), and the response of neural activity(Haddad et al., 2010; Lapid et al.,
2011). Such extensions will lead to more comprehensive and data-oriented
understanding of chemical information coding in the olfactory system.

##  5 Materials and methods

###  5.1 Coulomb matrix and its normalization

As a molecular representation, Coulomb matrix is calculated mainly based on
Coulomb replusion force as follows:

$C_{ij}=\left\\{\begin{array}[]{ll}0.5Z_{i}^{2.4}&\text{ for }i=j\\\
\frac{Z_{i}Z_{j}}{\left|\mathbf{R}_{i}-\mathbf{R}_{j}\right|}&\text{ for
}i\neq j\end{array}\right.,$ (1)

where $C_{ij}$ refers to the entry in $i^{th}$ row and $j^{th}$ column,
$Z_{i}$ refers to the atomic charge, and $R_{i}$ refers to the relative
coordinates.

Considering the gap of minimal and maximal entries in Coulomb matrix, we
perform 2 preprocessing methods to handle it, including matrix-wised Frobenius
normalization and minmax normalization as follows:

$\displaystyle||C||_{F}$
$\displaystyle=\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{n}C_{ij}^{2}},$ (2)
$\displaystyle C_{max}$ $\displaystyle=\max_{i,j}{C_{ij}},$ $\displaystyle
C_{min}$ $\displaystyle=\min_{i,j}{C_{ij}},$

where $||C||_{F}$ refers to the normalization term of Frobenius normalization,
$C_{max}$ and $C_{min}$ refer to the max and min term in minmax normalization.
With the normalization term, we calculated the preprocessed Coulomb matrix as
follows:

$\displaystyle C^{F}$ $\displaystyle=C/(||C||_{F}+\epsilon),$ (3)
$\displaystyle C^{M}$ $\displaystyle=(C-C_{min})/(C_{max}-C_{min}+\epsilon),$

where $C^{F}$ and $C^{M}$ refer to the Frobenuis- and minmax- normalized
matrix, and $\epsilon$ equals to ${10}^{-9}$.

###  5.2 GCN of Coulomb matrix for multi-label classification

We build Coulomb-GCN based on Coulomb matrix with three modules: 1) atom
embedding updating, 2) molecular embedding extraction, and 3) multi-label
classification. In atom embedding updating, we use GCN with the residual
mechanism to learn the molecular embedding as follows:

$H^{(l)}=\sigma(XH^{(l-1)}W_{graph}^{(l-1)})+H^{(l-1)}W_{linear}^{(l-1)},$ (4)

where $H^{(l)}\in R^{n\times d}$ refers to the updated atom embedding in
$l^{th}$ layer, $X\in R^{n\times n}$ refers to Coulomb matrix, $H^{(l-1)}\in
R^{n\times h}$ refers to the updated atom embedding in ${l-1}^{th}$ layer,
$W_{graph}^{(l-1)}\in R^{h\times d}$ refers to the parameters in GCN,
$W_{linear}^{(l-1)}\in R^{h\times d}$ refers to the parameters in the linear
transformation of $H^{(l-1)}$ for residual mechanism, $H^{(0)}$ refers to the
random initialized atom embedding, $n$ refers to the number of atom in
molecule, $d$ refers to the dimension of embedding in $l^{th}$ layer, $h$
refer to the dimension of embedding in ${l-1}^{th}$ layer, and $\sigma$ refers
to SELU activation function(Klambauer et al., 2017). In molecular embedding
extraction, we use the sum-pooling function as follows:

$m_{i}=\sum_{j\in[n]}H_{ji}^{(l)},$ (5)

where $m_{i}$ refers to the $i^{th}$ entry in the molecular embedding of
molecule $m$. In multi-label classification, we use fully-connected layers as
follows:

$y=\sigma(mW_{clf}),$ (6)

where $m\in R^{1\times d}$ refers to the learned molecular embedding,
$W_{clf}\in R^{d\times o}$ refers to the parameters in the fully-connected
layer, and $o$ refers to the number of odor descriptors. Coulomb matrix and
Coulomb-GCN have been implemented by Python (version 3.7.4) with deepchem
(version 2.6.1) and pytorch (version 1.12.1).

###  5.3 Directional graph modeling by Laplacian for multi-label
classification

Mol-PECO differs from Coulomb-GCN only in the atom embedding updating module,
where we build the directional graph modeling with learned positional encoding
(LPE).

The Laplacian matrix ($L^{1}$) of Coulomb matrix ($X$) is calculated as
follows:

$L^{1}=D-X,$ (7)

where $D$ refers to the degree matrix. In LPE with eigen-decomposition of
Laplacian, symmetrical Laplacian matrix is calculated for further spectral
decomposition as follows:

$\displaystyle L^{2}$ $\displaystyle=I-D^{-1/2}XD^{-1/2}$ (8)
$\displaystyle=D^{-1/2}DD^{-1/2}-D^{-1/2}XD^{-1/2}$
$\displaystyle=D^{-1/2}(D-X)D^{-1/2}$ $\displaystyle=D^{-1/2}L^{1}D^{-1/2}.$

where $I$ refers to the identity matrix. With the equations, $L_{2}$ is
calculated with $L_{1}$ divided by $\sqrt{d_{ii}d_{jj}}$. Let
$d_{ii}=\sum_{j=0}^{i-1}{a_{ij}}+\sum_{j=i+1}^{n}{a_{ij}}+a_{ii}$, then the
Laplacian is calculated as follows:

$\displaystyle L^{1}_{ij}$
$\displaystyle=\sum_{j=0}^{i-1}{a_{ij}}+\sum_{j=i+1}^{n}{a_{ij}}+a_{ii}-a_{ii},i=j,$
(9) $\displaystyle L^{1}_{ij}$ $\displaystyle=-a_{ij},i\neq j,$

which indicates setting the diagnoal entries zeros or not has no influence on
$L^{1}$ and $L^{2}$ as $L^{2}=D^{-1/2}L^{1}D^{-1/2}$.

Minimization of quadratic form on graphs acts as the cost function of
link/edge prediction and captures the global structural information.
Naturally, the eigen-decompostion of Laplacian provides the solutions as
follows:

$f^{*}=\min\limits_{f}fL^{2}f^{T}=1/2\sum_{i,j}X(i,j){(f_{i}-f_{j})}^{2},$
(10)

where $f^{*}$ refers to the eigenvectors of $L$. With Laplacian, Mol-PECO
updates the positional encoding in an atom-wise manner proposed in Spectral
Attention Network(Kreuzer et al., 2021). Mol-PECO learns LPE in an atom-by-
atom manner. Specifically, Mol-PECO first performs the linear transformation
of eigenvalues and eigenvectors in a single atom, learns the positional
information by Transformer, and then extracts the atom positional encoding by
sum-pooling as follows:

$\displaystyle W$ $\displaystyle=\Lambda W_{0},$ (11) $\displaystyle K$
$\displaystyle=Transformer(W),$ $\displaystyle e_{i}$
$\displaystyle=\sum_{j\in n}K_{ji},$

where $\Lambda\in R^{p\times 2}$ refers to the $p$ eigenvalues and
eigenvectors for a single atom, $W_{0}\in R^{2\times d}$ refers to the
parameters of the linear transformation in LPE, $K\in R^{p\times d}$ refers to
the updated embedding with Transformer of $W\in R^{p\times d}$, and $e_{i}$
refers to the $i^{th}$ entry of single atom embedding (vector-shaped
representation). Graph Laplacian and Mol-PECO have been implemented by Python
(version 3.7.4) with deepchem (version 2.6.1) and PyTorch (version 1.12.1).

###  5.4 Loss functions

For both Coulomb-GCN and Mol-PECO, we build the classification loss with the
binary cross-entropy loss function and a logarithm regularization term as
follows:

$\displaystyle l^{i}(y_{true}^{i},y_{pred}^{i})$
$\displaystyle=BCE(y_{true}^{i},y_{pred}^{i})+\lvert
log(y_{pred}^{i}+\epsilon)-log(y_{true}^{i}+\epsilon)\rvert,$ (12)
$\displaystyle l(y_{true},y_{pred})$
$\displaystyle=\frac{1}{o}\sum_{i\in[o]}w_{i}l^{i}(y_{true}^{i},y_{pred}^{i}),$

where $y_{true}^{i}$ and $y_{pred}^{i}$ refer to the ground truths and
predictions of $i^{th}$ odor descriptor, $l^{i}$ refers to the loss function
of $i^{th}$ odor descriptor, BCE refers to the binary cross-entropy function,
$w_{i}$ refers to $1-n_{pos}^{i}/n_{tot}$, $n_{pos}^{i}$ refers to the number
of positive samples in $i^{th}$ descriptor, $n_{tot}$ refers to the number of
total samples. The training process has been implemented by Python (version
3.7.4) with PyTorch (version 1.12.1).

###  5.5 Molecular descriptors

We include three classical molecular descriptors in this work as the baseline
molecular representations, including Mordred features (mordred)(Moriwaki et
al., 2018), bit-based Morgan fingerprints (bfps)(Rogers and Hahn, 2010), and
count-based Morgan fingerprints (cfps)(Rogers and Hahn, 2010). Mordred
calculates about 1825 features, including 214 2D and 1611 3D features. Both
bfps and cfps encode the molecule’s topological environments (molecular
fragments), which indicate the presence of atoms and functional groups, into a
vector. Specifically, bfps encodes presence or absence of the molecular
fragments as a binary information, while cfps encodes the number of
atom/functional-group in the topological environment. The calculation of
molecular descriptors has been implemented by Python (version 3.7.4) with
mordred (version 1.2.0) and rdkit-pypi (version 2022.3.4).

###  5.6 Machine learning of molecular descriptors

Machine learning of molecular descriptors used in this work as baselines
includes: 1) Synthetic Minority Over-sampling Technique (SMOTE)(Chawla et al.,
2002) for handling the imbalanced label distribution, and 2) multiple binary
classifiers for predicting multiple odor descriptors. For SMOTE, it performs
minority sampling by generating new minority instances to expand the number of
the minority class. For binary classifiers, we use K- Nearest Neighbor
classifier (KNN), random forest classifier (RF), and gradient boosting
classifier (GB). The KNN acts as a non-parametric classifier and predicts the
label by voting from neighbors. RF acts as the ensemble learning of decision
trees with sample bagging to decrease the variance of model and feature
bagging to decrease the correlation among decision trees. GB builds multiple
weak learners to minimize the differences between the true label and the
predicted value by performing gradient decent. The procedures have been
implemented by Python (version 3.7.4) with imblearn (version 0.9.0) and
scikit-learn (version 1.0.2).

###  5.7 Acknowledgments

We thank the suggestions and discussion with all colleagues in laboratory for
quantitative biology, IIS, Tokyo University and Takahiro G. Yamada, Yusuke
Shibuya, and Mamoru Tomiyama in Funahashi-lab, Keio University. We also thank
Yusuke Ihara and Chiori Ijichi for helpful comments. The work is supported by
China Scholarship Council (CSC) Grant 202106230237 to M. Z., JSPS KAKENHI
Grant Numbers 19H05799 to T.J.K., and JST CREST Grant Number JPMJCR2011 to
A.F. and T.J.K.

## References

  * Sobel et al. (1998) Noam Sobel, Vivek Prabhakaran, John E Desmond, Gary H Glover, RL Goode, Edith V Sullivan, and John DE Gabrieli.  Sniffing and smelling: separate subsystems in the human olfactory cortex.  _Nature_ , 392(6673):282–286, 1998.
  * Lapid et al. (2011) Hadas Lapid, Sagit Shushan, Anton Plotkin, Hillary Voet, Yehudah Roth, Thomas Hummel, Elad Schneidman, and Noam Sobel.  Neural activity at the human olfactory epithelium reflects olfactory perception.  _Nature neuroscience_ , 14(11):1455–1461, 2011\.
  * Sell (2006) Charles S Sell.  On the unpredictability of odor.  _Angewandte Chemie International Edition_ , 45(38):6254–6261, 2006.
  * Boesveldt et al. (2010) Sanne Boesveldt, Mats J Olsson, and Johan N Lundström.  Carbon chain length and the stimulus problem in olfaction.  _Behavioural brain research_ , 215(1):110–113, 2010.
  * Keller et al. (2017) Andreas Keller, Richard C Gerkin, Yuanfang Guan, Amit Dhurandhar, Gabor Turu, Bence Szalai, Joel D Mainland, Yusuke Ihara, Chung Wen Yu, Russ Wolfinger, et al.  Predicting human olfactory perception from chemical features of odor molecules.  _Science_ , 355(6327):820–826, 2017.
  * Polster and Schieberle (2017) Johannes Polster and Peter Schieberle.  Structure–odor correlations in homologous series of mercaptoalkanols.  _Journal of Agricultural and Food Chemistry_ , 65(21):4329–4340, 2017.
  * Pattanaik and Coley (2020) Lagnajit Pattanaik and Connor W Coley.  Molecular representation: going long on fingerprints.  _Chem_ , 6(6):1204–1207, 2020.
  * Rogers and Hahn (2010) David Rogers and Mathew Hahn.  Extended-connectivity fingerprints.  _Journal of chemical information and modeling_ , 50(5):742–754, 2010.
  * Moriwaki et al. (2018) Hirotomo Moriwaki, Yu-Shi Tian, Norihito Kawashita, and Tatsuya Takagi.  Mordred: a molecular descriptor calculator.  _Journal of cheminformatics_ , 10(1):1–14, 2018\.
  * Kipf and Welling (2016) Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Hofstetter et al. (2022) Albert Hofstetter, Lennard Böselt, and Sereina Riniker.  Graph-convolutional neural networks for (qm) ml/mm molecular dynamics simulations.  _Physical Chemistry Chemical Physics_ , 24(37):22497–22512, 2022.
  * Yang et al. (2019) Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, et al.  Analyzing learned molecular representations for property prediction.  _Journal of chemical information and modeling_ , 59(8):3370–3388, 2019.
  * Feinberg et al. (2018) Evan N Feinberg, Debnil Sur, Zhenqin Wu, Brooke E Husic, Huanghao Mai, Yang Li, Saisai Sun, Jianyi Yang, Bharath Ramsundar, and Vijay S Pande.  Potentialnet for molecular property prediction.  _ACS central science_ , 4(11):1520–1530, 2018\.
  * Wang et al. (2019) Xiaofeng Wang, Zhen Li, Mingjian Jiang, Shuang Wang, Shugang Zhang, and Zhiqiang Wei.  Molecule property prediction based on spatial graph embedding.  _Journal of chemical information and modeling_ , 59(9):3817–3828, 2019.
  * Withnall et al. (2020) Michael Withnall, Edvard Lindelöf, Ola Engkvist, and Hongming Chen.  Building attention and edge message passing neural networks for bioactivity and physical–chemical property prediction.  _Journal of cheminformatics_ , 12(1):1–18, 2020\.
  * Sanchez-Lengeling et al. (2019) Benjamin Sanchez-Lengeling, Jennifer N Wei, Brian K Lee, Richard C Gerkin, Alán Aspuru-Guzik, and Alexander B Wiltschko.  Machine learning for scent: Learning generalizable perceptual representations of small molecules.  _arXiv preprint arXiv:1910.10685_ , 2019.
  * Kreuzer et al. (2021) Devin Kreuzer, Dominique Beaini, Will Hamilton, Vincent Létourneau, and Prudencio Tossou.  Rethinking graph transformers with spectral attention.  _Advances in Neural Information Processing Systems_ , 34:21618–21629, 2021.
  * Oono and Suzuki (2019) Kenta Oono and Taiji Suzuki.  Graph neural networks exponentially lose expressive power for node classification.  _arXiv preprint arXiv:1905.10947_ , 2019.
  * Topping et al. (2021) Jake Topping, Francesco Di Giovanni, Benjamin Paul Chamberlain, Xiaowen Dong, and Michael M Bronstein.  Understanding over-squashing and bottlenecks on graphs via curvature.  _arXiv preprint arXiv:2111.14522_ , 2021.
  * Floriano et al. (2000) Wely B Floriano, Nagarajan Vaidehi, William A Goddard III, Michael S Singer, and Gordon M Shepherd.  Molecular mechanisms underlying differential odor responses of a mouse olfactory receptor.  _Proceedings of the National Academy of Sciences_ , 97(20):10712–10716, 2000.
  * Rupp et al. (2012) Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O Anatole Von Lilienfeld.  Fast and accurate modeling of molecular atomization energies with machine learning.  _Physical review letters_ , 108(5):058301, 2012\.
  * Dwivedi and Bresson (2020) Vijay Prakash Dwivedi and Xavier Bresson.  A generalization of transformer networks to graphs.  _arXiv preprint arXiv:2012.09699_ , 2020.
  * Beaini et al. (2021) Dominique Beaini, Saro Passaro, Vincent Létourneau, Will Hamilton, Gabriele Corso, and Pietro Liò.  Directional graph networks.  In _International Conference on Machine Learning_ , pages 748–758. PMLR, 2021.
  * Arctander et al. (1960) Steffen Arctander et al.  Perfume and flavor materials of natural origin.  _Perfume and Flavor Materials of Natural Origin._ , 1960.
  * Kumar et al. (2018) Yogesh Kumar, Om Prakash, Himanshu Tripathi, Sudeep Tandon, Madan M Gupta, Laiq-Ur Rahman, Raj K Lal, Manoj Semwal, Mahendra Pandurang Darokar, and Feroz Khan.  Aromadb: A database of medicinal and aromatic plant’s aroma molecules with phytochemistry and therapeutic potentials.  _Frontiers in plant science_ , 9:1081, 2018.
  * Garg et al. (2018) Neelansh Garg, Apuroop Sethupathy, Rudraksh Tuwani, Rakhi Nk, Shubham Dokania, Arvind Iyer, Ayushi Gupta, Shubhra Agrawal, Navjot Singh, Shubham Shukla, et al.  Flavordb: a database of flavor molecules.  _Nucleic acids research_ , 46(D1):D1210–D1216, 2018.
  * Acree (2004) T Acree.  Flavornet and human odor space.  _http://www. flavornet. org/flavornet. html_ , 2004.
  * (28) Fragrance Flavor.  Food, and cosmetics ingredients information. the good scents company.
  * Association et al. (2003) International Fragrance Association et al.  Gc/ms quantitation of potential fragrance allergens in fragrance compounds, 2003.
  * Leffingwell (2001) John C Leffingwell.  Olfaction ii.  _Leffingwell reports_ , 1(4):1–26, 2001.
  * Sharma et al. (2021) Anju Sharma, Rajnish Kumar, Shabnam Ranjta, and Pritish Kumar Varadwaj.  Smiles to smell: decoding the structure–odor relationship of chemical compounds using the deep neural network approach.  _Journal of Chemical Information and Modeling_ , 61(2):676–688, 2021.
  * Sharma et al. (2022) Anju Sharma, Bishal Kumar Saha, Rajnish Kumar, and Pritish Kumar Varadwaj.  Olfactionbase: a repository to explore odors, odorants, olfactory receptors and odorant–receptor interactions.  _Nucleic Acids Research_ , 50(D1):D678–D686, 2022\.
  * Corporation (2011) Sigma-Aldrich Corporation.  _Aldrich Chemistry 2012-2014: Handbook of Fine Chemicals_.  Sigma-Aldrich, 2011.
  * Szymański and Kajdanowicz (2017) Piotr Szymański and Tomasz Kajdanowicz.  A network perspective on stratification of multi-label data.  In _First International Workshop on Learning with Imbalanced Domains: Theory and Applications_ , pages 22–35. PMLR, 2017.
  * Schrier (2020) Joshua Schrier.  Can one hear the shape of a molecule (from its coulomb matrix eigenvalues)?  _Journal of Chemical Information and Modeling_ , 60(8):3804–3811, 2020.
  * Mahmood et al. (2021) Omar Mahmood, Elman Mansimov, Richard Bonneau, and Kyunghyun Cho.  Masked graph modeling for molecule generation.  _Nature communications_ , 12(1):1–12, 2021.
  * Wang et al. (2022) Yuyang Wang, Jianren Wang, Zhonglin Cao, and Amir Barati Farimani.  Molecular contrastive learning of representations via graph neural networks.  _Nature Machine Intelligence_ , 4(3):279–287, 2022.
  * Mohar et al. (1991) Bojan Mohar, Y Alavi, G Chartrand, and OR Oellermann.  The laplacian spectrum of graphs.  _Graph theory, combinatorics, and applications_ , 2(871-898):12, 1991.
  * Chung (1997) Fan RK Chung.  _Spectral graph theory_ , volume 92.  American Mathematical Soc., 1997.
  * Spielman (2012) Daniel Spielman.  Spectral graph theory.  _Combinatorial scientific computing_ , 18, 2012.
  * Chawla et al. (2002) Nitesh V Chawla, Kevin W Bowyer, Lawrence O Hall, and W Philip Kegelmeyer.  Smote: synthetic minority over-sampling technique.  _Journal of artificial intelligence research_ , 16:321–357, 2002.
  * Block (1992) Eric Block.  The organosulfur chemistry of the genus allium–implications for the organic chemistry of sulfur.  _Angewandte Chemie International Edition in English_ , 31(9):1135–1178, 1992.
  * Miękus et al. (2020) Natalia Miękus, Krystian Marszałek, Magdalena Podlacha, Aamir Iqbal, Czesław Puchalski, and Artur H Świergiel.  Health benefits of plant-derived sulfur compounds, glucosinolates, and organosulfur compounds.  _Molecules_ , 25(17):3804, 2020.
  * Li et al. (2020) Pengyong Li, Jun Wang, Yixuan Qiao, Hao Chen, Yihuan Yu, Xiaojun Yao, Peng Gao, Guotong Xie, and Sen Song.  Learn molecular representations from large-scale unlabeled molecules for drug discovery.  _CoRR_ , abs/2012.11175, 2020.  URL https://arxiv.org/abs/2012.11175.
  * Hu et al. (2019) Weihua Hu, Bowen Liu, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay S. Pande, and Jure Leskovec.  Pre-training graph neural networks.  _CoRR_ , abs/1905.12265, 2019.  URL http://arxiv.org/abs/1905.12265.
  * Trimmer et al. (2019) Casey Trimmer, Andreas Keller, Nicolle R Murphy, Lindsey L Snyder, Jason R Willer, Maira H Nagai, Nicholas Katsanis, Leslie B Vosshall, Hiroaki Matsunami, and Joel D Mainland.  Genetic variation across the human olfactory receptor repertoire alters odor perception.  _Proceedings of the National Academy of Sciences_ , 116(19):9475–9480, 2019.
  * Keller et al. (2007) Andreas Keller, Hanyi Zhuang, Qiuyi Chi, Leslie B Vosshall, and Hiroaki Matsunami.  Genetic variation in a human odorant receptor alters odour perception.  _Nature_ , 449(7161):468–472, 2007.
  * Meister (2015) Markus Meister.  On the dimensionality of odor space.  _Elife_ , 4:e07865, 2015.
  * Pannunzi and Nowotny (2019) Mario Pannunzi and Thomas Nowotny.  Odor stimuli: not just chemical identity.  _Frontiers in physiology_ , 10:1428, 2019.
  * Bhandawat et al. (2005) Vikas Bhandawat, Johannes Reisert, and King-Wai Yau.  Elementary response of olfactory receptor neurons to odorants.  _Science_ , 308(5730):1931–1934, 2005.
  * Mainland et al. (2015) Joel D Mainland, Yun R Li, Ting Zhou, Wen Ling L Liu, and Hiroaki Matsunami.  Human olfactory receptor responses to odorants.  _Scientific data_ , 2(1):1–9, 2015.
  * Haddad et al. (2010) Rafi Haddad, Tali Weiss, Rehan Khan, Boaz Nadler, Nathalie Mandairon, Moustafa Bensafi, Elad Schneidman, and Noam Sobel.  Global features of neural activity in the olfactory system form a parallel code that predicts olfactory behavior and perception.  _Journal of Neuroscience_ , 30(27):9017–9026, 2010.
  * Klambauer et al. (2017) Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter.  Self-normalizing neural networks.  _CoRR_ , abs/1706.02515, 2017.  URL http://arxiv.org/abs/1706.02515.

##  Appendix A Appendix

###  A.1 Parameter tuning of Mol-PECO

We optimize Mol-PECO with different number of Transformer layers in learned
positional encoding and choose $4$ as the optimized parameter (Table 3).

Table 3: Scores of Mol-PECO with different number of Transformer layers and
the layer number optimization. The highest scores are shown in bold.

#Transformer layers | AUROC1 | AUPRC1 | Precision1 | Recall1 | Specificity1 | Accuracy1
---|---|---|---|---|---|---
1 | 0.809 | 0.170 | 0.103 | 0.826 | 0.789 | 0.807
2 | 0.798 | 0.147 | 0.100 | 0.817 | 0.779 | 0.798
3 | 0.807 | 0.163 | 0.090 | 0.831 | 0.779 | 0.805
4 | 0.813 | 0.181 | 0.104 | 0.819 | 0.797 | 0.808
5 | 0.807 | 0.181 | 0.089 | 0.819 | 0.780 | 0.800
6 | 0.802 | 0.162 | 0.093 | 0.814 | 0.790 | 0.802

  * 1

The evaluation metrics are calculated with the validation set.

The training process of the optimized parameter (the number of Transformer
layer = $4$) is monitored by the loss curves (Figure 6a). The decreased trends
of the loss curves are consistent, indicating that Mol-PECO has found the
right bias/variance tradeoff in the training and validation sets. Mol-PECO
chooses the checkpoint of $554$ epochs with minimal loss of $0.140$ as the
final model for testing (Figure 6b).

Figure 6: The training process with the optimized parameter. (a) The
unweighted loss during training. (b) The unweighted AUROC during training.

After training and validation, we obtain the detailed performances of $118$
odor descriptors in $6$ evaluation metrics (Figure 7).

Figure 7: The detailed performances of Mol-PECO in 6 evaluation metrics,
including (a) AUROC, (b) AUPRC, (c) Precision, (d) Recall, (e) Specificity,
and (f) accuracy. The x-axis refers to the ratio of positive samples. Each dot
refers to one odor descriptor.

###  A.2 Supplementary materials for the learned odor space

Despite of the high-frequency descriptors, the structure-correlated
descriptors, and the synonymy descriptors, we also investigate the alcohol-
related descriptors, the fruit-related descriptors, and the synonym
descriptors in Figure 8. For alcohol-related descriptors, ‘ethereal’ and
‘winey’ possess the most molecules, leading to an obvious cluster in nearby
locations (Figure 8a). For the fruit-related descriptors, ‘apple’, ‘banana’,
‘pear’, and ‘pineapple’ are clustered together (Figure 8b). For synonymy
descriptors, partial ‘roasted’ and full ‘cooked’ are clustered into neighbors,
validating their similar semantics (Figure 8c).

Figure 8: Global view of learned odor space with dimensionality reduction by
t-SNE in (a) the alcohol-related descriptors (‘alcoholic’, ‘fermented’,
‘ethereal’, ‘winey’, ‘rum’, and ‘brandy’), (b) the fruit-related descritors
(‘apple’, ‘banana’, ‘pear’, and ‘pineapple’), and (c) the molecules with
similar semantic meaning (‘roasted’ and ‘cooked’).


# Gain assisted controllable fast light generation in cavity magnomechanics

Sanket Das  Department of Physics, Indian Institute of Technology Guwahati,
Guwahati-781039, Assam, India.  Subhadeep Chakraborty  Centre for Quantum
Engineering Research and Education, TCG Centres for Research and Education in
Science and Technology, Sector V, Salt Lake, Kolkata 70091, India  Tarak N.
Dey<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology Guwahati, Guwahati-781039, Assam, India.

###### Abstract

We study the controllable output field generation from a cavity
magnomechanical resonator system that consists of two coupled microwave
resonators. The first cavity interacts with a ferromagnetic yttrium iron
garnet (YIG) sphere providing the magnon-photon coupling. Under passive
cavities configuration, the system displays high absorption, prohibiting
output transmission even though the dispersive response is anamolous. We
replace the second passive cavity with an active one to overcome high
absorption, producing an effective gain in the system. We show that the
deformation of the YIG sphere retains the anomalous dispersion. Further,
tuning the exchange interaction strength between the two resonators leads to
the system’s effective gain and dispersive response. As a result, the
advancement associated with the amplification of the probe pulse can be
controlled in the close vicinity of the magnomechanical resonance.
Furthermore, we find the existence of an upper bound for the intensity
amplification and the advancement of the probe pulse that comes from the
stability condition. These findings may find potential applications for
controlling light propagation in cavity magnomechanics.

##  I Introduction

Cavity magnonics Rameshti  _et al._ (2022); Schlawin  _et al._ (2022), has
become an actively pursued field of research due to its potential application
in quantum information processing Lachance-Quirion  _et al._ (2019); Hisatomi
_et al._ (2016). The key constituent to such systems is a ferrimagnetic
insulator with high spin density and low damping rate. It also supports
quantized magnetization modes, namely, the magnons Zhang  _et al._ (2016);
Tabuchi  _et al._ (2014). With strongly coupled magnon-photon modes, cavity
magnonics is an excellent platform for studying all the strong-coupling cavity
QED effects Zhang  _et al._ (2014). Besides originating from the shape
deformation of the YIG, the magnon can also couple to a vibrational or phonon
mode Zhang  _et al._ (2016). This combined setup of magnon-photon-phonon
modes, namely the cavity magnomechanics, has already demonstrated
magnomechanically induced transparency Zhang  _et al._ (2016), magnon-induced
dynamical backaction Potts  _et al._ (2021), magnon-photon-phonon entanglement
Li  _et al._ (2018); Li and Gröblacher (2021), squeezed state generation Li
_et al._ (2019), magnomechanical storage and retreival of a quantum state
Sarma  _et al._ (2021).

Recently, $\mathcal{PT}$-symmetry drew extensive attention to elucidate the
dynamics of a coupled system characterized by gain and loss Bender and
Boettcher (1998); Bender  _et al._ (2002). Here, $\mathcal{P}$ stands for the
parity operation, that results in an interchange between the two constituent
modes of the system. The time reversal operator $\mathcal{T}$ takes $i$ to
$-i$. $\mathcal{PT}$-symmetry demands the Hamiltonian is commutative with the
joint $\mathcal{PT}$ operators i.e., $[H,PT]=0$. This system possesses a
spectrum of entirely real and imaginary eigenvalues that retain
distinguishable characteristics El-Ganainy  _et al._ (2018). The point
separating these two eigenvalues is the exceptional point (EP) Heiss (2012)
where the two eigenvalues coalesce, and the system degenerates. A natural
testbed for $\mathcal{PT}$-symmetric Hamiltonian is optical as well as quantum
optical systems Makris  _et al._ (2008); Musslimani  _et al._ (2008); Guo  _et
al._ (2009) which already led to the demonstration of some of the exotic
phenomena, like nonreciprocal light propagation Peng  _et al._ (2014),
unidirectional invisibility Longhi (2011); Feng  _et al._ (2013), optical
sensing and light stopping Goldzak  _et al._ (2018). Very recently, a
tremendous effort has been initiated to explore non-Hermitian physics in
magnon assisted hybrid quantum systems. The second-order exceptional point is
detected in a two-mode cavity-magnoic system, where the gain of the cavity
mode is accomplished by using the idea of coherent perfect absorption Zhang
_et al._ (2017). The concept of Anti-$\mathcal{PT}$ symmetry has been realized
experimentally Zhao  _et al._ (2020), where the adiabatic elimination of the
cavity field produces dissipative coupling between two magnon modes. Beyond
the unique spectral responses, these non-Hermitian systems can manipulate the
output microwave field transmission Wang and Gou (2023); Kong  _et al._
(2019). The underlying mechanism behind such an application is magnetically
induced transparencyUllah  _et al._ (2020); Zhang  _et al._ (2016), where the
strong magnon-photon coupling produces a narrow spectral hole inside the probe
absorption spectrum. Further studies in this direction establish the
importance of the weak magnon-phonon coupling to create double transmission
windows separated by an absorption peak. Moreover, manipulating the absorption
spectrum is also possible by varying the amplitude and phase of the applied
magnetic field Li  _et al._ (2020).

It is well established over the past decade that optomechanically induced
transparency (OMIT) Agarwal and Huang (2010); Safavi-Naeini  _et al._ (2011);
Weis  _et al._ (2010) is an essential tool for investigating slow light Hau
_et al._ (1999) and light storage Fiore  _et al._ (2011, 2013) in cavity. In
addition, incorporating $\mathcal{PT}$-symmetry in optomechanical systems,
provides a better controllability of light transmission Jing  _et al._ (2015);
Li  _et al._ (2016) and produces subluminal to superluminal light conversion.
Nonetheless, their proposals may find experimental challenges as the gain of
the auxiliary cavity can lead the whole system to instability Vanderhaegen
_et al._ (2021). An eminent advantage of the magnomechanical system over the
optomechanical system is that it offers strong hybridization between the
magnon-photon mode. The magnomechanical systems offer better tunability as an
external magnetic field can vary the magnon frequency. Exploiting these
advantages, a $\mathcal{PT}$-symmetry-like magnomechanical system can be
constructed by resonantly driving the YIG sphere to an active magnon mode Huai
_et al._ (2019). The controllable sideband generation with tunable group delay
can be feasible by changing the power of the control field.

This paper investigates a controllable advancement and transmission of the
microwave field from a coupled cavity magnomechanical system. Optical coupling
between a passive cavity resonator containing YIG sphere and a gain-assisted
auxiliary cavity can form a coupled cavity resonator. An external drive has
been used to deform the YIG sphere’s shape, resulting in the magnon-phonon
interaction in the passive cavity. We show how the gain of the auxiliary
cavity helps to overcome absorptive behaviour in our hybrid system. As a
result, the output microwave field amplifies at the resonance condition.
Moreover, the weak magnon-phonon interaction exhibits anomalous dispersion
accompanied by a gain spectrum, demonstrating superluminal light. We also
examine how the slope of the dispersion curve can be controlled by tuning the
photon hopping interaction strength between the two cavities.

The paper is organized as follows. In Section II, a theorical model for the
coumpound cavity magnomechanical system with $\mathcal{PT}$-symmetric
resonator is described. The Heisenserg equations of motion to govern the
expectation values of operators of every system are derived in this Section.
In Section III.1, we analyse the stability criteria of the model system and
examine the effect of the auxiliary cavity gain on the absorptive and
dispersive response of the system in Section III.2. Section III.3 discusses
the output probe field transmission. Further, the group velocity of the
optical probe pulse has been studied analytically and verified numerically in
Section III.4. Finally, we draw our conclusions in Section IV.

##  II Theoretical Model

Recently, there has been a growing interest in realizing a gain in different
components of cavity magnonics systems Zhang  _et al._ (2017); Huai  _et al._
(2019). In this work, we investigate the effect of medium gain on the probe
response and its transmission. The system under consideration is a hybrid
cavity magnomechanical system that consists of two coupled microwave cavity
resonators. One of the resonators is passive and contains a YIG sphere inside
it. We refer to this resonator as a cavity magnomechanical (CMM) resonator.
Applying a uniform bias magnetic field to the YIG sphere excites the magnon
mode. The magnon mode, in turn, couples with the cavity field by the magnetic-
dipole interaction.

Figure 1: The schematic diagram of a hybrid cavity magnomechanical system. The
system consists of two coupled microwave cavities. One of them is passive, and
another one is active. The passive cavity contains a ferromagnetic YIG sphere
inside it. The applied bias magnetic field produces the magnetostrictive
interaction between magnon and phonon. The coupling rates between the magnon-
photon and magnon-phonon are $g_{ma}$ and $g_{mb}$, respectively. Strong
control field of frequency $\omega_{l}$ and a weak probe field of frequency
$\omega_{p}$ are applied to the passive cavity.

Nonetheless, the external bias magnetic field results in shape deformation of
the YIG sphere, leading to the magnon-phonon interaction. The second resonator
(degenerate with the first one) is coupled to the first resonator via optical
tunnelling at a rate $J$. Two input fields drive the first resonator. The
amplitude of the control, $\varepsilon_{l}$, and probe fields,
$\varepsilon_{p}$, are given by
$\varepsilon_{i}=\sqrt{P_{i}/\hslash\omega_{i}},(i\in l,p)$ with
$P_{i}~{}\textrm{and}~{}\omega_{i}$ being the power and frequency of the
respective input fields. The Hamiltonian of the combined system can be written
as

$\displaystyle
H=\hslash\omega_{c}a_{1}^{\dagger}a_{1}+\hslash\omega_{c}a_{2}^{\dagger}a_{2}+\hslash\omega_{m}m^{\dagger}m+\hslash\omega_{b}b^{\dagger}b$
$\displaystyle+\hslash J(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1})+\hslash
g_{ma}(a_{1}^{\dagger}m+a_{1}m^{\dagger})$ $\displaystyle+\hslash
g_{mb}m^{\dagger}m(b^{\dagger}+b)+i\hslash\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{l}(a_{1}^{\dagger}e^{-i\omega_{l}t}-a_{1}e^{i\omega_{l}t})$
$\displaystyle+i\hslash\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{p}(a_{1}^{\dagger}e^{-i\omega_{p}t}-a_{1}e^{i\omega_{p}t}),$
(1)

where the first four terms of the Hamiltonian describe the free energy
associated with each system’s constituents. The constituents of our model are
characterized by their respective resonance frequencies: $\omega_{c}$ for the
cavity mode, $\omega_{m}$ for the magnon mode, $\omega_{b}$ for the phonon
mode. The annihilation operators for the cavity, magnon and phonon modes are
represented by $a_{i}$, ($i=1,2$), $m$ and $b$, respectively. The fifth term
signifies the photon exchange interaction between the two cavities with
strength, $J$. The sixth term of the Hamiltonian corresponds to the
interaction between the magnon and photon modes, characterized by a coupling
rate $g_{ma}$. The interaction between the magnon and phonon modes is
described by the seventh term of the Hamiltonian and the coupling rate between
magnon and phonon mode is $g_{mb}$. Finally, the last two terms arise due to
the interaction between the cavity field and two input fields. The cavity,
magnon and phonon decay rates are characterized by
$\kappa_{1},\kappa_{m}~{}\textrm{and}~{}\kappa_{b}$, respectively. The
coupling between the CMM resonator and the output port is given by
$\eta_{a}=\kappa_{c_{1}}/2\kappa_{1}$, where $\kappa_{c_{1}}$ is the cavity
external decay rate. In particular, we will consider the CMM resonator to be
working in the critical-coupling regime where $\eta_{a}$ is $1/2$. At this
point, it is convenient to move to a frame rotating at $\omega_{l}$. Following
the transformation $H_{rot}=RHR^{\dagger}+i\hslash\left({\partial R}/{\partial
t}\right)R^{\dagger}$ with
$R=e^{i\omega_{l}(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2}+m^{\dagger}m)t}$,
the Hamiltonian in Eq. (II) can be rewritten as

$\displaystyle H_{rot}$
$\displaystyle=\hbar\Delta_{a}(a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2})+\hbar\Delta_{m}m^{\dagger}m+\hbar\omega_{b}b^{\dagger}b$
$\displaystyle+\hbar J(a_{1}^{\dagger}a_{2}+a_{2}^{\dagger}a_{1})+\hbar
g_{ma}(a_{1}^{\dagger}m+a_{1}m^{\dagger})$ $\displaystyle+\hslash
g_{mb}m^{\dagger}m(b^{\dagger}+b)+i\hbar\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{l}(a_{1}^{\dagger}-a_{1})$
$\displaystyle+i\hbar\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{p}(a_{1}^{\dagger}e^{-i\delta
t}-h.c),$ (2)

where $\Delta_{a}=\omega_{c}-\omega_{l}$ ($\Delta_{m}=\omega_{m}-\omega_{l}$)
and $\delta=\omega_{p}-\omega_{l}$ are, respectively, the cavity (magnon) and
probe detuning. The mean response of the system can be obtained by the
Heisenberg - Langevin equation as
$\langle\mathcal{\dot{O}}\rangle=i/\hbar\langle[H_{rot},\mathcal{O}]\rangle+\langle
N\rangle$. Further, we consider the quantum fluctuations $(N)$ as white noise.
Then starting form Eq. II, the equations of motion of the system can be
expressed as

$\displaystyle\langle\dot{a_{1}}\rangle$
$\displaystyle=(-i\Delta_{a}-\kappa_{1})\langle a_{1}\rangle-ig_{ma}\langle
m\rangle-iJ\langle a_{2}\rangle$
$\displaystyle+\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{l}+\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{p}e^{-i\delta
t},$ $\displaystyle\langle\dot{m}\rangle$
$\displaystyle=(-i\Delta_{m}-\kappa_{m})\langle m\rangle-ig_{ma}\langle
a_{1}\rangle$ $\displaystyle-ig_{mb}\langle m\rangle(\langle
b^{\dagger}\rangle+\langle b\rangle),$ $\displaystyle\langle\dot{b}\rangle$
$\displaystyle=(-i\omega_{b}-\kappa_{b})\langle b\rangle-ig_{mb}\langle
m^{\dagger}\rangle\langle m\rangle,$ $\displaystyle\langle\dot{a_{2}}\rangle$
$\displaystyle=(-i\Delta_{a}+\kappa_{2})\langle a_{2}\rangle-iJ\langle
a_{1}\rangle,$ (3)

where $\kappa_{2}$ and $\kappa_{b}$ respectively denote the gain of the second
resonator and phonon damping rates. We note that $\kappa_{2}>0$ corresponds to
a coupled passive-active CMM resonators system and $\kappa_{2}<0$ describes a
passive-passive coupled CMM resonators system. Assuming the control field
amplitude $\varepsilon_{l}$ to be larger than the probe field
$\varepsilon_{p}$, each operator expectation values
$\langle\mathcal{O}(t)\rangle$ can be decomposed into its steady-state values
$\mathcal{O}_{s}$ and a small fluctuating term $\delta\mathcal{O}(t)$. The
steady-state values of each operator are

$\displaystyle a_{1s}$
$\displaystyle=\frac{(-i\Delta_{a}+\kappa_{2})(-ig_{ma}m_{s}+\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{l})}{(i\Delta_{a}+\kappa_{1})(-i\Delta_{a}+\kappa_{2})-J^{2}},$
(4a) $\displaystyle m_{s}$
$\displaystyle=\frac{-ig_{ma}a_{1s}}{i\Delta_{m}^{\prime}+\kappa_{m}},$ (4b)
$\displaystyle b_{s}$
$\displaystyle=\frac{-ig_{mb}|m_{s}|^{2}}{i\omega_{b}+\kappa_{b}},$ (4c)
$\displaystyle a_{2s}$
$\displaystyle=\frac{iJa_{1s}}{(-i\Delta_{a}+\kappa_{2})}.$ (4d)

While the fluctuating parts of Eq. II can be expressed as

$\displaystyle\delta\dot{a}_{1}$
$\displaystyle=-\left(i\Delta_{a}+\kappa_{1}\right)\delta a_{1}-iJ\delta
a_{2}-ig_{ma}\delta m$
$\displaystyle+\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{p}e^{-i\delta t},$
$\displaystyle\delta\dot{m}$
$\displaystyle=-(i\Delta^{\prime}_{m}+\kappa_{m})\delta m-ig_{ma}\delta
a_{1}-iG\delta b-iG\delta b^{\dagger},$ $\displaystyle\delta\dot{b}$
$\displaystyle=-\left(i\omega_{b}+\kappa_{b}\right)\delta b-iG\delta
m^{\dagger}-iG^{*}\delta m,$ $\displaystyle\delta\dot{a}_{2}$
$\displaystyle=-\left(i\Delta_{a}-\kappa_{2}\right)\delta a_{2}-iJ\delta
a_{1},$ (5)

where $\Delta_{m}^{\prime}=\Delta_{m}+g_{mb}(b_{s}+b_{s}^{*})$ is the
effective magnon detuning and $G=g_{mb}m_{s}$ is the enhanced magnon-phonon
coupling strength. For simplicity, we express these fluctuation equations as

$\displaystyle i\frac{d}{dt}|\psi\rangle$
$\displaystyle=H_{eff}|\psi\rangle+F,$ (6)

where the fluctuation vector $|\psi\rangle=(\delta a_{1},\delta
a_{1}^{\dagger},\delta a_{2},\delta a_{2}^{\dagger},\delta b,\delta
b^{\dagger},\delta m,\delta m^{\dagger})^{T}$, input field
$F=(\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{p}e^{-i\delta
t},\sqrt{2\eta_{a}\kappa_{1}}\varepsilon_{p}e^{i\delta t},0,0,0,0,0,0)^{T}$.
Next, we adopt the following ansatz to solve Eq. II:

$\displaystyle\delta a_{1}(t)$ $\displaystyle=A_{1+}e^{-i\delta
t}+A_{1-}e^{i\delta t},$ $\displaystyle\delta m(t)$
$\displaystyle=M_{+}e^{-i\delta t}+M_{-}e^{i\delta t}$ $\displaystyle\delta
b(t)$ $\displaystyle=B_{+}e^{-i\delta t}+B_{-}e^{i\delta t},$
$\displaystyle\delta a_{2}(t)$ $\displaystyle=A_{2+}e^{-i\delta
t}+A_{2-}e^{i\delta t}.$ (7)

Here $A_{i+}~{}\textrm{and}~{}A_{i-}$ correspond to the $i^{th}$ cavity
generated probe field amplitude and the four-wave mixing field amplitude,
respectively. By considering
$h_{1}=-i\Delta_{a}+i\delta-\kappa_{1},~{}h_{2}=-i\Delta_{a}-i\delta-\kappa_{1},~{}h_{3}=-i\Delta_{a}+i\delta+\kappa_{2},~{}h_{4}=-i\Delta_{a}-i\delta+\kappa_{2},~{}h_{5}=-i\omega_{b}+i\delta-\kappa_{b},~{}h_{6}=-i\omega_{b}-i\delta-\kappa_{b},~{}h_{7}=-i\Delta_{m}^{\prime}+i\delta-\kappa_{m},~{}h_{8}=-i\Delta_{m}^{\prime}-i\delta-\kappa_{m}$,
we obtain $A_{1+}$ which corresponds to the output probe field amplitude from
the CMM resonator as

$\displaystyle A_{1+}(\delta)=\frac{C(\delta)}{D(\delta)},$ (8)

where

$\displaystyle
C(\delta)=-\sqrt{2\eta_{a}\kappa_{a}}\varepsilon_{p}h_{3}(h_{5}h_{7}h_{6}^{*}(J^{2}h_{8}^{*}+h_{4}^{*}(g_{ma}^{2}+h_{2}^{*}h_{8}^{*}))$
$\displaystyle+|G|^{2}(h_{5}-h_{6}^{*})(J^{2}(h_{7}-h_{8}^{*})-h_{4}^{*}(g_{ma^{2}}+h_{2}^{*}(h_{8}^{*}-h_{7}^{*})))),$
$\displaystyle
D(\delta)=h_{5}h_{6}^{*}(g_{ma}^{2}h_{3}+h_{7}(h_{1}h_{3}+J^{2}))$
$\displaystyle(J^{2}h_{8}^{*}+h_{4}^{*}(g_{ma}^{2}+h_{2}^{*}h_{8}^{*}))+|G|^{2}(h_{5}-h_{6}^{*})$
$\displaystyle(J^{2}(g_{ma}^{2}h_{3}+(h_{1}h_{3}+J^{2})(h_{7}-h_{8}^{*}))-h_{4}^{*}$
$\displaystyle((h_{1}h_{3}+J^{2})(g_{ma}^{2}-h_{2}^{*}(h_{7}-h_{8}^{*}))-h_{2}^{*}h_{3}g_{ma}^{2})).$
(9)

The output field from the CMM resonator is obtained by the cavity input-output
relation

$\displaystyle\varepsilon_{out}$
$\displaystyle=\sqrt{2\eta_{a}\kappa_{1}}\langle
a_{1}\rangle-\varepsilon_{l}-\varepsilon_{p}e^{-i\delta t}.$ (10)

By substituting Eq. II into Eq. 10, we obtain the normalized output probe
field intensity from the CMM resonator as

$\displaystyle T=|t_{p}|^{2}$
$\displaystyle=\left|\frac{\sqrt{2\eta_{a}\kappa_{1}}A_{1+}}{\varepsilon_{p}}-1\right|^{2}.$
(11)

In order to numerically simulate the transmitted output probe field spectrum,
we use the following experimentally realizable set of parameter values Nair
_et al._ (2022); Zhang  _et al._ (2016). The degenerate microwave cavities of
frequency $\omega_{c}/2\pi=7.86$ GHz. The decay rate of the first cavity is
$\kappa_{1}/\pi=3.35$ MHz. The spin density $\rho=4.22\times 10^{27}$ m-3 and
the diameter of the YIG sphere $D=25$ $\mu$m. It results in $3\times 10^{16}$
number of spins ($N_{m}$) present in the YIG sphere. The phonon mode has
frequency $\omega_{b}/2\pi=11.42$ MHz with decay rate $\kappa_{b}/\pi=300$ Hz,
and the magnon-phonon coupling strength $g_{mb}/2\pi$ is $1$ Hz. The Kittel
mode frequency of the YIG sphere is $\omega_{m}=\gamma_{e}B_{0,i}$, with
gyromagnetic ratio, $\gamma_{e}/2\pi=28$ GHz/T and $B_{0,i}$ is the input bias
magnetic field amplitude. The magnon decay rate is $\kappa_{m}=3.52$ MHz.
Magnon-photon coupling strength $g_{ma}=\gamma_{e}B_{vac}\sqrt{5N_{m}}/2$ can
be controlled by changing the vacuum magnetic field amplitude as
$B_{vac}=\sqrt{2\pi\hbar\omega_{c}/V}$.

##  III Results

###  III.1 Stability Analysis

Initially we consider the two coupled cavities which are operating under a
balanced gain-loss condition. The Hamiltonian describing such coupled
resonator system ($g_{ma}=g_{mb}=0$) can be written as

$\displaystyle H_{cav}$ $\displaystyle=\hbar(\Delta_{a}-i\kappa_{1})\delta
a_{1}^{\dagger}\delta a_{1}+\hbar(\Delta_{a}+i\kappa_{1})\delta
a_{2}^{\dagger}\delta a_{2}$ $\displaystyle+\hbar J(\delta
a_{1}^{\dagger}\delta a_{2}+\delta a_{2}^{\dagger}\delta a_{1}).$ (12)

The eigenvalues of $H_{cav}$ are
$\lambda_{\pm}=\Delta_{a}\pm\sqrt{J^{2}-\kappa_{1}^{2}}$. Note that the above
Hamiltonian remains invariant under the simultaneous parity
$\mathcal{P}:a_{1}\leftrightarrow a_{2}$ and time-reversal operation
$\mathcal{T}:i\rightarrow-i$ operations, and, its eigenvalues are entirely
real and complex for $J>\kappa_{1}$ and $J<\kappa_{2}$. The point
$J=\kappa_{1}$, which marks this transition from $\mathcal{PT}$ symmetric to
the $\mathcal{PT}$ breaking phase, is known as the exceptional point (EP). One
must understand the competitive behaviour between the inter-cavity field
coupling and the loss/gain rates to get insight into this transition. For
$J>\kappa_{1}$ the intracavity field amplitudes can be coherently exchanged
and thus give rise to a coherent oscillation between the field amplitudes.
However, for $J<\kappa_{1}$ the intracavity field can not be transferred to
the other one, resulting in a strong field localization or in other words
exponential growth. A quick look at Eq. 4(a) also suggests such gain-induced
dynamic instability in $a_{1}$ at $J=\kappa_{1}$ for $\Delta_{a}=0$. This
situation becomes more complicated in the presence of magnon-photon coupling.
Now, the combined system ($g_{ma},g_{mb}\neq 0$) ceases to become
$\mathcal{PT}$ symmetric. However, the effect of an additional gain cavity
($\kappa_{2}>0$) can be understood by analyzing the stability diagram of the
whole system. In the following, we derive the stability condition by invoking
the Routh-Hurwitz criterion which requires all the eigenvalues of $H_{eff}$
have negative real parts. The magenta region of Fig. 2 suggests that when
$g_{ma}$ is small the instability threshold remains close to the
$J=\kappa_{1}$ (the conventional EP for a binary $\mathcal{PT}$ symmetric
system). However, with increasing $g_{ma}$ the system reaches instability at a
larger exchange interaction $J$. Such a restriction over the choice of the
photon exchange rate parameter $J$ will be followed throughout this paper.

Figure 2: The stable and unstable regions are determined as a function of
normalized evanescent coupling strength ($J/\kappa_{1}$) and the cavity-magnon
coupling strength ($g_{ma}$) when the loss of the CMM-resonator is perfectly
balanced by the gain of the auxiliary cavity $(\kappa_{1}=\kappa_{2})$. We
consider the control field intensity to be $10$ mW. The other parameters are
$\omega_{c}=2\pi\times 7.86$ MHz, $\omega_{b}=2\pi\times 11.42$ MHz,
$\Delta_{a}=\Delta_{m}^{\prime}=\omega_{b}=2\pi\times 11.42$ MHz,
$\kappa_{1}=\kappa_{2}=\pi\times 3.35$ MHz, $\kappa_{m}=\pi\times 1.12$ MHz,
$\kappa_{b}=\pi\times 300$ Hz, and $g_{mb}=2\pi$ Hz.

###  III.2 Absorption and dispersion spectrum

The magnomechanical system under consideration corresponds to the level
diagram of Fig. 3. Application of a probe field excites the passive cavity
mode and allows the transition between $|1\rangle~{}\textrm{and}~{}|2\rangle$.
The exchange interaction, $J$, couples two degenerate excited states
$|2\rangle$ and $|5\rangle$. The presence of the strong control field
distributes the population between the two states $|2\rangle$ and $|3\rangle$.
The magnon-phonon coupling, $g_{mb}$, couples both the metastable ground
states $|3\rangle$ and $|4\rangle$. Here, we consider both the microwave
cavities to be passive ($\kappa_{2}<0$) under a weak magnon-photon coupling
strength, $g_{ma}$. In this situation, the magnon-photon hybridization is
insignificant. The absorptive and dispersive response can be quantified by the
real and imaginary components of $(t_{p}+1)$ that will be presented as
$\alpha$ and $\beta$, respectively. In Fig. 4(a), we present the absorptive
response of the system as a function of normalized probe detuning. The black-
solid curve depicts a broad absorption spectrum of the probe field when the
exchange interaction is much weaker than the cavity decay rate.

Figure 3: Level diagram of the model system.
$|n_{i}\rangle$,$|m\rangle~{}\textrm{and}~{}|b\rangle$ represents the photon
number state of $i^{th}$ cavity, magnon mode and phonon mode, respectively.
The application of strong control field to the CMM resonator couples
$|n_{1}+1,n_{2},m,b\rangle$ and $|n_{1},n_{2},m+1,b\rangle$, whereas, the
presence of a weak probe field increases the photon number of CMM resonator by
unity. $g_{mb}$ couples
$|n_{1},n_{2},m+1,b\rangle~{}\textrm{and}~{}|n_{1},n_{2},m,b+1\rangle$. The
hopping interaction between the two cavities directly couples
$|n_{1}+1,n_{2},m,b\rangle\leftrightarrow|n_{1},n_{2}+1,m,b\rangle$.

One can explain it by considering the level diagram of Fig. 3, where the
initial population stays in the ground state $|1\rangle$. Applying a probe
field transfers the population from the ground state to the excited state,
$|2\rangle$. In addition, the weak magnon-photon coupling (with respect to
$\kappa_{1}$) restricts a significant transition from $|2\rangle$ to
$|3\rangle$. As a result, it allows the transfer of a fraction of the excited
state’s population by invoking the exchange interaction $J$.

Figure 4: (a) Absorption and (b) dispersion spectrum of the model system. The
slope of the dispersion curve is shown in the inset. Here we consider both the
microwave cavities are passive, with identical decay rates
$(\kappa_{1}=-\kappa_{2})$. The magnon-photon couplong strength, $g_{ma}$ is
taken as $2$ MHz. All the other parameters are mentioned earlier.

The increase in the exchange interaction strength causes a gradual decrease in
$|2\rangle$’s population. It reduces the absorption coefficient around the
resonance condition except for $\delta=\omega_{b}$. This phenomenon is shown
by the red-dashed and the blue dotted-dashed curve of Fig. 4(a). We observe a
narrow absorption peak inside the broad absorption peak for $J=0.4$
$\kappa_{1}$. The sharp absorption peak, exactly at $\delta=\omega_{b}$,
occurs due to the magnomechanical resonance. Further increasing the exchange
interaction virtually cuts off the population distribution from $|2\rangle$ to
$|3\rangle$. As a result, the effect of magnon-phonon resonance also
decreases, and the absorption peak at $\delta=\omega_{b}$ eventually
diminishes. In Fig. 4(b), we present the dispersion spectrum as a function of
normalized detuning $\delta/\omega_{b}$. For the time being, we neglect the
effect of magnomechanical coupling and observe the occurrence of anomalous
dispersion around $\delta=\omega_{b}$ for $J=0.4$ $\kappa_{1}$. Further,
increasing the exchange interaction strength more significant than the cavity
decay rate can alter the dispersive response from anomalous to normal, as
shown by the red-dashed and blue-dot-dashed curves. In the inset of Fig. 4(b),
we plot the slope of the temporal dispersion $d\beta/d\delta$ at the extreme
vicinity of the magnon-phonon resonance condition. The positive values of the
slope of the temporal dispersion signify anomalous dispersion due to the
magnomechanical coupling. However, the steepness of the dispersion curve can
be reduced by increasing the exchange interaction strength, as shown by the
red-dashed and blue-dot-dashed curves. Note that this dispersion curve is
accompanied by absorption. Output transmission of the probe field is
prohibited in the presence of huge absorption. Therefore, reducing absorption
or introducing the gain to the system is mandatory for observing the group
velocity phenomena.

Figure 5: (a) Absorption and (b) dispersion spectrum of the model system. The
slope of the dispersion curve is shown in the inset. Here we consider the
second cavity as a gain cavity, with $\kappa_{2}=\kappa_{1}$. All the other
parameters are the same as before.

To achieve reasonable transmission at the output, we replace the auxiliary
passive cavity with an active one where the second cavity’s gain
($\kappa_{2}>0$) completely balances the first cavity’s loss. In this
scenario, the stability criterion for the hybrid system allows us to consider
the exchange interaction strength $J$ to be greater than $1.053$ $\kappa_{1}$
for $g_{ma}=2$ MHz. We present the absorptive response of the model system in
Fig. 5(a). The black solid curve of Fig. 5(a) illustrates the occurrence of a
double absorption peak spectrally separated by a broad gain regime. The
graphical nature is determined by the roots of $D(\delta)$, which are, in
general, complex. The real parts of the roots determine the spectral peak
position, and the imaginary parts correspond to their widths. To illustrate
this, we consider $J=1.30$ $\kappa_{1}$ with all other parameters remaining
the same as earlier. The real parts of the root of $D(\delta)$ present two
distinct normal mode positions at $\delta/\omega_{b}$ values
$0.88~{}\textrm{and}~{}1.12$. The other two normal modes are spectrally
located at the same position $\delta/\omega_{b}=1$. The interference between
these two normal modes becomes significant while approaching the stability
bound as depicted by the red dashed and blue dot-dashed curve of Fig. 5(a). In
turn, it reduces the overall gain of the composite system. Further, we
investigate the effect of a gain-assisted auxiliary cavity on the medium’s
dispersive response in Fig. 5(b). For $J=1.30$ $\kappa_{1}$, the two
absorption peaks produce two distinct anomalous dispersion regions separated
by a broad normal dispersive window. Weakening the exchange interaction
strength reveals prominent normal dispersion around the resonance condition
except for $\delta=\omega_{b}$, and the window shrinks. In the inset of Fig.
5(b), we present the slope of the dispersive response due to the
magnomechanical resonance. The black solid curve of Fig. 5(b) suggests the
occurrence of anomalous dispersion at the magnon-phonon resonance condition.
Moreover, one can increase the steepness of the dispersion curve by simply
approaching the instability threshold, as delineated by the red-dashed and
blue-dotted-dashed curve of the inset of Fig. 5(b). In the consecutive
section, we will discuss how the change in the dispersion curve can produce
controllable group velocity of the light pulses through the medium and
investigate the role of the exchange interaction.

###  III.3 Output probe transmission

The output probe intensity from the system depends on its absorptive response.
Equation 11 dictate the transmission of the probe field and is presented in
Fig. 6. For Fig. 6(a), we consider both the microwave cavities as passive ones
with identical decay rates, i.e., $\kappa_{1}=-\kappa_{2}$. The black solid
curve shows a broad absorptive response for $J=0.40$ $\kappa_{1}$. Increasing
the exchange interaction strength causes gradual enhancement in the output
probe transmission, as delineated by the red-dashed and blue dot-dashed curve
of Fig. 6(a), and the absorption window splits into two parts. A precise
observation confirms the presence of extremely weak transmission dip exactly
at $\delta=\omega_{b}$ for all the three exchange interaction strengths under
consideration.

Figure 6: Exchange interaction $J$ dependent normalized output probe
transmission is plotted as a function of normalized detuning between the
control and the probe field when (a) both the cavities are passive ones, and
(b) one is active and another one is passive.

In Fig. 6(b), we present the advantage of using a gain-assisted auxiliary
cavity along with a CMM resonator to obtain a controllable amplification of
the output probe field. We begin our discussion considering the photon hopping
interaction, $J=1.30$ $\kappa_{1}$. The black solid curve of Fig. 6(b)
estimates the normalized probe transmission of $6.03$. Here the normalization
is done with respect to the input probe field intensity. By decreasing the
parameter $J$, we approach the unstable region and observe the occurrence of a
double transmission peak separated by a sharp and narrow transmission dip. The
amplitude of the double transmission peak demonstrates the probe pulse
amplification by a factor of $830$, as presented by the blue dotted-dashed
curve. However, an explicit observation suggests the output probe field
amplification by a factor of $67$ at the resonance condition
$\delta=\omega_{b}$. The physics behind the probe field amplification can be
well understood as: Introduction of gain to the second cavity compensates a
portion of losses in the first cavity through $J$. This leads to an enhanced
field amplitude in the first cavity. In the presence of moderate magnon-photon
coupling it increases the effective magnon-photon coupling strength. Hence, we
observe a higher transmission at the two sidebands but also find a large
transmission dip at $\delta=\omega_{b}$.

###  III.4 Group delay

Controllable group delay has gained much attention due to its potential
application in quantum information processing and communication. The
dispersive nature of the medium is the key to controlling the group delay of
the light pulse under the assumption of low absorption or gain. The pulse with
finite width in the time domain is produced by superposing several independent
waves with different frequencies centered around a carrier frequency
$(\omega_{s})$. The difference in time between free space propagation and a
medium propagation for the same length can create a group delay. The
analytical expression for the group delay can be constructed by considering
the envelope of the optical pulse as

$f(t_{0})=\int_{-\infty}^{\infty}\tilde{f}(\omega)e^{-i\omega
t_{0}}d\omega,\\\ $

where $\tilde{f}(\omega)$ corresponds to the envelope function in the
frequency domain. Accordingly, the reflected output probe pulse can be
expressed as

$\displaystyle f^{R}(t_{0})=$
$\displaystyle\int_{-\infty}^{\infty}t_{p}(\omega)\tilde{f}(\omega)e^{-i\omega
t_{0}}d\omega,$ (13) $\displaystyle=$ $\displaystyle
e^{-i\omega_{s}t_{0}}\int_{-\infty}^{\infty}t_{p}(\omega_{s}+\delta)\tilde{f}(\omega_{s}+\delta)e^{-i\delta
t_{0}}d\delta,$ $\displaystyle=$ $\displaystyle
t_{p}(\omega_{s})e^{-i\omega_{s}\tau_{g}}f(t_{0}-\tau_{g}).$ (14)

This expression can be obtained by expanding $t_{p}(\omega_{s}+\delta)$ in the
vicinity of $\omega_{s}$ by a Taylor series and keeping the terms upto first
order in $\delta$. An expression for time-delay is obtained as Safavi-Naeini
_et al._ (2011); Das and Dey (2022)

$\displaystyle\tau_{g}={\textrm{R}e}\left[\frac{-i}{t_{p}(\omega_{s})}\left(\frac{dt_{p}}{d\omega}\right)\bigg{|}_{\omega_{s}}\right],$
(15)

which can be further simplified as

$\displaystyle\tau_{g}=\frac{{\left(\alpha(\omega_{s})-1\right)}{\frac{d\beta}{d\omega}}\bigg{|}_{\omega_{s}}-{\beta(\omega_{s})}{\frac{d\alpha}{d\omega}}\bigg{|}_{\omega_{s}}}{|t_{p}(\omega_{s})|^{2}}.$
(16)

From Eq. 16, the slope of the medium’s absorption and dispersion curves
determine the probe pulse propagation delay or advancement. However, Fig. 5(b)
suggests that the value of $\beta$ is negligibly small near the
magnomechanical resonance. Hence, the group delay depends on the first term of
the numerator of Eq. 16.

Figure 7: Time delay of the probe pulse for different evanescent coupling
strength $J$ have been plotted against the normalized probe detuning
$\delta/\omega_{b}$, while the control power is $10$ mW. All other parameters
are taken as the same as in Fig. 5.

In Fig. 7, we examine the effect of photon-photon exchange interaction on the
probe pulse propagation delay when both cavities operate under balanced gain-
loss condition. The system produces anomalous dispersion accompanied by a gain
response. The black solid curve of Fig. 7 depicts the probe pulse advancement
of $2.4$ $\mu$s for the photon hopping interaction strength, $J=1.3$
$\kappa_{1}$. Moreover, one can enhance the effective gain and the slope of
the anomalous dispersion curve by approaching the instability threshold. That,
in turn, brings out the super luminosity of the output probe pulse,
characterized by the advancement of $17.9$ $\mu$s as shown by the blue dotted-
dashed curve of Fig. 7.

To verify the above results, we consider a Gaussian probe pulse with a finite
width around the resonance condition, i.e., $\delta=\omega_{b}$, and
numerically integrate it by using Eq. 13. The shape of the input envelope is
considered as,

$\tilde{f}(\omega)=\frac{\varepsilon_{p}}{\sqrt{\pi\Gamma^{2}}}e^{-\frac{(\omega-\omega_{s})^{2}}{\Gamma^{2}}},$
(17)

where $\Gamma$ is the spectral width of the optical pulse. We consider
$\Gamma$ to be $7.17$ kHz, such that the Gaussian envelope is well-contained
inside the gain-window around the resonance condition ($\delta=\omega_{b}$),
as depicted in Fig. 5(a). The dispersive, absorptive as well as gain response
of the system can be demonstrated by examining the effect of the probe
transmission coefficient ($t_{p}$) on the shape of the input envelope. The
gain of the auxiliary cavity manipulates the probe transmission coefficient in
such a way that it amplifies the intensity of the output probe pulse. The
black solid curve depicts the output probe pulse amplification of $6.2$ for
photon-hopping interaction strength, $J=1.30$ $\kappa_{1}$. A decrease in the
$J$ value gradually enhances the effective gain in the system.

Figure 8: The relative intensity of the output probe pulse is plotted against
the normalized time ($\Gamma$t) for different photon-photon exchange
interaction strength when both cavities are operating under balanced gain-loss
condition.

It amplifies the output probe transmission as presented by the red dashed and
blue dotted-dashed curves of Fig. 8. We observe that the output field
amplification can reach to a factor of 65.3 while considering the exchange
interaction strength to be $1.07$ $\kappa_{1}$. Further decreasing the
exchange interaction will lead to dynamical instability in our model system.
Interestingly, the temporal width of the probe pulse is almost unaltered
during the propagation through the magnon-assisted double cavity system. This
numerical result agrees with our analytical results for the output probe
transmission, as shown in Fig. 6(b). Moreover, the importance of the photon-
photon exchange interaction on the probe pulse propagation advancement can be
observed from the inset of Fig. 8. The peak separation between the input pulse
$(t=0)$ and the output pulse for $J=1.30$ $\kappa_{1}$ estimates the probe
pulse advancement of $2.34$ $\mu$s. The red dashed, and blue dashed-dotted
curve of the inset estimates the probe pulse advancement of $8.75$ $\mu s$ and
$13.30$ $\mu s$ for $J=1.10$ $\kappa_{1}~{}\textrm{and}~{}1.07$ $\kappa_{1}$,
respectively.

##  IV Conclusion

In conclusion, we have theoretically investigated the controllable output
field transmission from a critically coupled cavity magnomechanical system. We
drive the first cavity with a YIG sphere inside it, establishing the magnon-
photon coupling. The photon exchange interaction connects the second microwave
cavity with the first. An external magnetic field induces the deformation
effect of the YIG sphere. In this study, the interaction between the magnon
and photon modes lies under the weak coupling regime. The medium becomes
highly absorbent when both cavities are passive, and the output probe
transmission is prohibited. We introduce a gain to the auxiliary cavity to
overcome this situation. It is noteworthy that the instability threshold must
be close to the conventional exceptional point for a binary
$\mathcal{PT}$-symmetric system. At the magnomechanical resonance, the
auxiliary cavity produces an effective gain associated with anomalous
dispersion. Further, decreasing the photon exchange interaction strength
causes gradual enhancement of the effective gain and the steepness of the
dispersion spectrum. As a result, we observe a controllable superluminal
microwave pulse propagation associated with amplification by a factor of 67.
By studying the propagation dynamics of a Gaussian probe pulse of width $7.17$
kHz, we confirm that the numerical study is consistent with the analytical
results. Our study may find potential applications in weak signal sensing and
communication purposes in the newly emerging field of cavity magnomechanics.

## References

  * Rameshti  _et al._ (2022) B. Z. Rameshti, S. V. Kusminskiy, J. A. Haigh, K. Usami, D. Lachance-Quirion, Y. Nakamura, C.-M. Hu, H. X. Tang, G. E. Bauer, and Y. M. Blanter, Physics Reports 979, 1 (2022).
  * Schlawin  _et al._ (2022) F. Schlawin, D. M. Kennes, and M. A. Sentef, Applied Physics Reviews 9, 011312 (2022).
  * Lachance-Quirion  _et al._ (2019) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Applied Physics Express 12, 070101 (2019).
  * Hisatomi  _et al._ (2016) R. Hisatomi, A. Osada, Y. Tabuchi, T. Ishikawa, A. Noguchi, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. B 93, 174427 (2016).
  * Zhang  _et al._ (2016) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Science Advances 2, e1501286 (2016), https://www.science.org/doi/pdf/10.1126/sciadv.1501286 .
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 113, 083603 (2014).
  * Zhang  _et al._ (2014) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * Potts  _et al._ (2021) C. A. Potts, E. Varga, V. A. S. V. Bittencourt, S. V. Kusminskiy, and J. P. Davis, Phys. Rev. X 11, 031053 (2021).
  * Li  _et al._ (2018) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. Lett. 121, 203601 (2018).
  * Li and Gröblacher (2021) J. Li and S. Gröblacher, Quantum Science and Technology 6, 024005 (2021).
  * Li  _et al._ (2019) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. A 99, 021801(R) (2019).
  * Sarma  _et al._ (2021) B. Sarma, T. Busch, and J. Twamley, New Journal of Physics 23, 043041 (2021).
  * Bender and Boettcher (1998) C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998).
  * Bender  _et al._ (2002) C. M. Bender, D. C. Brody, and H. F. Jones, Phys. Rev. Lett. 89, 270401 (2002).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Physics 14, 11 (2018).
  * Heiss (2012) W. D. Heiss, Journal of Physics A: Mathematical and Theoretical 45, 444016 (2012).
  * Makris  _et al._ (2008) K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Phys. Rev. Lett. 100, 103904 (2008).
  * Musslimani  _et al._ (2008) Z. H. Musslimani, K. G. Makris, R. El-Ganainy, and D. N. Christodoulides, Phys. Rev. Lett. 100, 030402 (2008).
  * Guo  _et al._ (2009) A. Guo, G. J. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. A. Siviloglou, and D. N. Christodoulides, Phys. Rev. Lett. 103, 093902 (2009).
  * Peng  _et al._ (2014) B. Peng, Ş. K. Özdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Nature Physics 10, 394 (2014).
  * Longhi (2011) S. Longhi, Journal of Physics A: Mathematical and Theoretical 44, 485302 (2011).
  * Feng  _et al._ (2013) L. Feng, Y.-L. Xu, W. S. Fegadolli, M.-H. Lu, J. E. B. Oliveira, V. R. Almeida, Y.-F. Chen, and A. Scherer, Nature Materials 12, 108 (2013).
  * Goldzak  _et al._ (2018) T. Goldzak, A. A. Mailybaev, and N. Moiseyev, Phys. Rev. Lett. 120, 013901 (2018).
  * Zhang  _et al._ (2017) D. Zhang, X.-Q. Luo, Y.-P. Wang, T.-F. Li, and J. Q. You, Nature Communications 8, 1368 (2017).
  * Zhao  _et al._ (2020) J. Zhao, Y. Liu, L. Wu, C.-K. Duan, Y.-x. Liu, and J. Du, Phys. Rev. Appl. 13, 014053 (2020).
  * Wang and Gou (2023) F. Wang and C. Gou, Opt. Lett. 48, 1164 (2023).
  * Kong  _et al._ (2019) C. Kong, H. Xiong, and Y. Wu, Phys. Rev. Appl. 12, 034001 (2019).
  * Ullah  _et al._ (2020) K. Ullah, M. T. Naseem, and Ö. E. Müstecaplioğlu, Phys. Rev. A 102, 033721 (2020).
  * Li  _et al._ (2020) X. Li, W.-X. Yang, T. Shui, L. Li, X. Wang, and Z. Wu, Journal of Applied Physics 128, 233101 (2020), https://doi.org/10.1063/5.0028395 .
  * Agarwal and Huang (2010) G. S. Agarwal and S. Huang, Phys. Rev. A 81, 041803(R) (2010).
  * Safavi-Naeini  _et al._ (2011) A. H. Safavi-Naeini, T. P. M. Alegre, J. Chan, M. Eichenfield, M. Winger, Q. Lin, J. T. Hill, D. E. Chang, and O. Painter, Nature 472, 69 (2011).
  * Weis  _et al._ (2010) S. Weis, R. Rivière, S. Deléglise, E. Gavartin, O. Arcizet, A. Schliesser, and T. J. Kippenberg, Science 330, 1520 (2010), https://www.science.org/doi/pdf/10.1126/science.1195596 .
  * Hau  _et al._ (1999) L. V. Hau, S. E. Harris, Z. Dutton, and C. H. Behroozi, Nature 397, 594 (1999).
  * Fiore  _et al._ (2011) V. Fiore, Y. Yang, M. C. Kuzyk, R. Barbour, L. Tian, and H. Wang, Phys. Rev. Lett. 107, 133601 (2011).
  * Fiore  _et al._ (2013) V. Fiore, C. Dong, M. C. Kuzyk, and H. Wang, Phys. Rev. A 87, 023812 (2013).
  * Jing  _et al._ (2015) H. Jing, Ş. K. Özdemir, Z. Geng, J. Zhang, X.-Y. Lü, B. Peng, L. Yang, and F. Nori, Scientific Reports 5, 9663 (2015).
  * Li  _et al._ (2016) W. Li, Y. Jiang, C. Li, and H. Song, Scientific Reports 6, 31095 (2016).
  * Vanderhaegen  _et al._ (2021) G. Vanderhaegen, C. Naveau, P. Szriftgiser, A. Kudlinski, M. Conforti, A. Mussot, M. Onorato, S. Trillo, A. Chabchoub, and N. Akhmediev, Proceedings of the National Academy of Sciences 118, e2019348118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2019348118 .
  * Huai  _et al._ (2019) S.-N. Huai, Y.-L. Liu, J. Zhang, L. Yang, and Y.-x. Liu, Phys. Rev. A 99, 043803 (2019).
  * Nair  _et al._ (2022) J. M. P. Nair, D. Mukhopadhyay, and G. S. Agarwal, Phys. Rev. B 105, 214418 (2022).
  * Das and Dey (2022) S. Das and T. N. Dey, J. Opt. Soc. Am. B 39, 859 (2022).