# Modelling magnetism of C at O and B monovacancies in graphene

T. P. Kaloni1, M. Upadhyay Kahaly1, R. Faccio2,3, and U. Schwingenschlögl1,
Corresponding author. Tel.: +966 544700080. E-mail address:
<EMAIL_ADDRESS>(U. Schwingenschlogl) 1PSE Division, KAUST,
Thuwal 23955-6900, Kingdom of Saudi Arabia
2Crystallography, Solid State and Materials Laboratory (Cryssmat-Lab), DETEMA,
Facultad de Química, Universidad de la República, Gral. Flores 2124, P.O. Box
1157, Montevideo, Uruguay
3Centro NanoMat, Polo Tecnológico de Pando, Facultad de Química, Universidad
de la República, Cno. Aparicio Saravia s/n, 91000, Pando, Canelones, Uruguay

###### Abstract

The presence of defects can introduce important changes in the electronic
structure of graphene, leading to phenomena such as C magnetism. In addition,
vacancies are reactive and permit the incorporation of dopants. This paper
discusses the electronic properties of defective graphene for O and B
decoration. Phonon calculations allow us to address directly the stability of
the systems under study. We show that it is possible to obtain magnetic
solutions with and without dangling bonds, demonstrating that C magnetism can
be achieved in the presence of B and O.

1\. Introduction

C nanostructures have attracted the attention of the scientific community
because of both unique fundamental properties and potential in technological
applications. In particular, the high coherence length and large conductivity
castro  of C materials are promising in electronics and spintronics. On the
other hand, also magnetism has been observed in this class of materials
Palacio , generating huge interest of both experiment and theory han ; Setze ;
rode ; lehtinen ; Flipse1 ; Ugeda ; mousumi . The structural, electronic, and
magnetic properties of defective graphene have been studied in detail and it
has been confirmed that vacancies are the source of magnetism new1 ; new2 ;
new3 ; new4 . It has been demonstrated that point defects lead to notable
paramagnetism but no magnetic ordering is achieved down to liquid helium
temperature nair . In addition, the point defects carry magnetic moments
irrespective of the vacancy concentration. In palacio1  the authors report
that the extended $\pi$-band magnetism reduces to zero in the limit of
monovacancies in graphene. For these reasons, it is important to evaluate the
role of saturation and chemical modifications in the neighborhood of mono and
multivacancies in order to understand how a persistent magnetic ordering can
be achieved.

Graphene samples in general are characterized by the presence of $sp^{2}$
hybridized C atoms, extended pores, and various attached functional groups, as
revealed by scanning electron microscopy, atomic force microscopy, and Raman
spectroscopy. It has been concluded that the appearance of magnetism requires
the presence of defects, adatoms, or topological defects khana . Zagaynova and
coworkers have prepared magnetic C by chemical oxidation mombr . B doping in
general leads to different magnetic responses depending on the amount of
doping. This behavior has been partially addressed by Faccio and coworkers
pardo2 , who have demonstrated that magnetism requires the proximity of B and
a monovacancy. However, when the B atom is too close to the vacancy all the
dangling bonds are reconstructed and magnetism is suppressed. The authors of
kaloni  have addressed the situation of multiple O atoms attached to a
monovacancy in graphene, demonstrating that vacancies are magnetic if and only
if they are metallic and non-magnetic if and only if they are semiconducting.
Metallicity and magnetism thus are simultaneously determined by the presence
or absence of dangling C bonds after the oxidation. In our present work, we
give a theoretical study of the effects of vacancies interacting with
different amounts of B and O atoms, focusing on the structural, electronic,
and magnetic properties. An experimental realization of the proposed
structures in possible along the lines of Ref. carbon .

2\. Methods

We employ density functional theory in the generalized gradient approximation
(Perdew-Burke-Ernzerhof scheme) as implemented in the Quantum-ESPRESSO code
paolo . All calculations are performed with a plane wave cutoff energy of 544
eV. A Monkhorst-Pack $8\times 8\times 1$ k-mesh is used to relax the
structures and a $16\times 16\times 1$ k-mesh to calculate the density of
states (DOS) with a high accuracy. We use a $5\times 5$ supercell of pristine
graphene in our calculations cheng1 . This supercell size has shown in various
studies to be sufficient for functionalization of graphene by atoms, small
molecules, amine and nitrobezene groups, and others new12 ; new13 ; new14 ;
new15 ; new16 ; new17 . Our supercell contains 50 C atoms and has a lattice
constant of $a=12.2$ Å with a 20 Å thick vacuum layer on top. The atomic
positions are relaxed upto an energy convergence of $10^{-7}$ eV and a force
convergence of 0.005 eV/Å. For our vibrational study, the phonon frequencies
and eigenvectors at the $\Gamma$-point are calculated with an energy error of
less than $10^{-14}$ eV. Phonon frequencies at the $\Gamma$-point are
determined by density functional perturbation theory for evaluating the
structural stability Mod. . We note that the numerics break the symmetry of
the dynamical matrix and introduce slight errors in the phonon frequencies of
maximal 15 cm-1, which, however, are not critical.

Figure 1:  Crystal structures for 1 B and 2, 3, or 4 O atoms adsorbed at a
monovacancy in graphene. The yellow, red, and gray spheres represent C, O, and
B, respectively. Figure 2:  Crystal structures for 2 O and 2, 3, 4, or 6 B
atoms adsorbed at a monovacancy in graphene. The yellow, red, and gray spheres
represent C, O, and B, respectively.

In a first step, we create a monovacancy in our $5\times 5$ supercell and
relax the system. Afterwards, we add O and B atoms in the vicinity of the
monovacancy for different starting geometries and relax the system again.
Several prototypical monovacancies are prepared in order to consider various
possibilities with different O and B concentrations. We name these
configurations giving first the number of O atoms and second the number of B
atoms. For example, configuration 1/2 has 1 O and 2 B atoms in the vicinity of
the monovacancy. All the configurations under study are shown in Figs. 1 and
2. Note that we have tested the required supercell size by addressing for the
most stable configurations $10\times 10$ supercells and find no significant
modification in the induced buckling of the graphene sheet after O and B
adsorption.

3\. Results and discussion

At a clean monovacancy the under-coordinated C atom is subject to Jahn-Teller
distortion Jahn  and therefore moves out of the graphene plane by 0.13 Å. We
find a total magnetic moment of 1.35 $\mu_{B}$, where 1 $\mu_{B}$ is due to
the localized dangling bond of one C atom and the remaining moment is carried
by extended $\pi$ states kaloni . The vacancy formation energy, $E_{form}$, is
determined by the expression

$E_{form}=E_{vacancy}-\frac{N-1}{N}E_{graphene},$ (1)

where $E_{vacancy}$ and $E_{graphene}$ are the total energies of defective and
pristine graphene, respectively. Moreover, $N$ is the number of atoms in the
pristine graphene. The obtained value of $E_{form}$ is 7.4 eV, which agrees
well with previous reports pandey ; eggie .

We first study a monovacancy decorated by 1 O atom to which 1 B atom is added.
When forcing the O atom to occupy a position within the graphene plane
connected to 3 C atoms the B atom is released and repelled, leading to
configuration 1/1a. In configuration 1/1b the B atom connects to 1 O atom and
2 C atoms, saturating all dangling bonds. Moreover, configuration 1/1c deals
with B adsorption far from the vacancy, for which we obtain a bridge position
over a C$-$C bond with $d_{B-C}=1.87$ Å. The energetics of the doping process
indicate that the chemical reaction between graphene-oxide and B can be
evaluated as

$E_{form}=E_{vacancy+B+O}-\frac{N-N_{C}}{N}E_{graphene}-\frac{N_{O}}{2}E_{O_{2}}-N_{B}E_{B,bulk},$
(2)

where $E_{vacancy+B+O}$ is the total energy after B adsorption on graphene-
oxide. In addition, $N_{C}$, $N_{O}$, and $N_{B}$ are the numbers of missing C
atoms and additional O and B atoms in the decorated system. It turns out that
configuration 1/1b is energetically favorable with $E_{form}=-6.9$ eV, see
Table I, while for configurations 1/1a and 1/1c we obtain in each case a value
of $-1.7$ eV. By the huge energy difference it is clear that only
configuration 1/1b will be realized. Due to the B adsorption and the induced C
and O reconstructions, this configuration becomes a metal, see Fig. 3(a). The
Dirac cone, which is still due to the C $p_{z}$ orbitals with no contribution
of O or B, splits and shifts above $E_{F}$ by about 0.6 eV. Interestingly, a
magnetic moment of 0.2 $\mu_{B}$ is observed although no dangling bonds are
present.

System |  $E_{form}$ (eV) |  Magnetic moment ($\mu_{B}$)
---|---|---
1/1a |  $-$1.7 | 1.0
1/1b |  $-$6.9 | 0.2
1/1c |  $-$1.7 | 0.0
2/1a |  $-$7.0 | 0.0
2/1b |  $-$7.4 | 0.1
2/1c |  $-$5.3 | 1.0
2/1d |  $-$4.8 | 0.0
2/1e |  $-$6.9 | 0.5∗
3/1a |  $-$4.4 | 0.2
3/1b |  $-$3.4 | 1.0
3/1c |  $-$3.5 | 1.0
4/1a | 1.1 | 0.0
4/1b |  $-$11.3 | 0.0
4/1c | 1.6 | 0.0
2/2 |  $-$8.4 | 1.0
2/3 |  $-$12.8 | 1.0
2/4 |  $-$22.9 | 0.0
2/6 |  $-$34.4 | 0.0
Table 1: Formation energy and magnetic moment. ∗The magnetic moment is located
on the released CO molecule. Figure 3:  Electronic band structure and DOS for
the structures (a) 1/1b, (b) 2/1b, (c) 3/1a, and (d) 4/1b.

Turning to an oxidized monovacancy decorated by 2 O atoms, in configuration
2/1a the B atom establishes strong interaction with O and C, which originally
was under-coordinated. We obtain $E_{form}=-7.0$ eV and no spin polarization
due to saturation of all dangling C bonds. In configurations 2/1b and 2/1c the
B atom is adsorbed at the boundary of the vacancy, leading to reconstructions
with one and two C=O bonds, respectively. We obtain for $E_{form}$ values of
$-7.4$ eV and $-5.3$ eV, where configuration 2/1c is less favorable due to the
presence of a dangling C bond which induces a magnetic moment of 1 $\mu_{B}$.
In configuration 2/1d the B atom is adsorbed far away from the monovacancy,
which yields a surprising geometrical reconstruction: the two C=O bonds change
into one ketone (C=O) and one ether group. However, with $E_{form}=-4.8$ eV
this configuration is not favorable. Finally, in configuration 2/1e the B atom
is located in the graphene plane and a CO molecule is released. We obtain
$E_{form}=-6.9$ eV. In Fig. 3(b) we address the electronic band structure of
the energetically favorable configuration 2/1b. We find that one spin majority
band and one spin minority band cross $E_{F}$, both contributed by the C
$p_{z}$ orbitals. The semiconducting state of a monovacancy decorated by 2 O
atoms is transferred into a metallic state by B adsorption.

For B adsorption at a monovacancy decorated by 3 O atoms, configuration 3/1a
is characterized by three ketone groups, establishing a BO3 unit with an
interatomic distance of $d_{B-O}=1.38$ Å. In configuration 3/1b the B atom is
attached to one ketone group and one C terminated ether group, which changes
the original $sp^{2}$ hybridization into a $sp^{3}$ hybridization, with
interatomic distances of $d_{B-O}=1.35$ Å and $d_{B-C}=1.62$ Å. Since in
configuration 3/1c the B atom substitutes a C atom near the vacancy, the
structure does not change significantly. Configuration 3/1a is clearly
favorable with $E_{form}=-4.4$ eV, as compared to values of $-3.4$ eV and
$-3.5$ eV for configurations 3/1b and 3/1c. For configuration 3/1a we find a
magnetic moment of 0.6 $\mu_{B}$. In addition, Fig. 3(c) shows that this
configuration 3/1a is metallic with two distinct C $p_{z}$ derived bands
crossing $E_{F}$. The observed localized states can be attributed to the
under-coordination of the B atom, which clearly distinguishes this case from
configurations 3/1b and 3/1c.

In the case of 4 adsorbed O atoms, three configurations are found to be
stable. Configuration 4/1a has two ketone groups and one peroxide
($-$O$-$O$-$) group and configuration 4/1b one ketone, one peroxide, and one
ether group. The origin of the stability is the fact that the BO3 group itself
is very stable new6 . In the former case the B atom is adsorbed far away from
to the monovacancy, while in the latter it bonds with the 2 O atoms to form a
BO3 group. Configuration 4/1c is similar to configuration 4/1a but with B
substituting C at the boundary of the vacancy. Regarding the energetics,
configuration 4/1b is strongly favorable with $E_{form}=-11.3$ eV, as compared
to values of $1.1$ eV and $1.6$ eV, since a peroxide group is present. The
electronic structure shown in Fig. 3(d) indicates that this configuration is
semiconducting with a band gap of 0.5 eV at the K point, and not magnetic.

We next evaluate the effect of the B concentration on the electronic structure
of oxidized graphene by substituting 2 to 6 B atoms on C sites for a fixed O
concentration (two adsorbed O atoms), see Fig. 2. In configuration 2/2 two B
atoms replace two C atoms at the boundary of the vacancy. The system develops
a magnetic moment of 1 $\mu_{B}$ due to the presence of a dangling C bond. In
configuration 2/3 the substitution of an additional B atom does not change
this situation. The fourth B atom in configuration 2/4 establishes a linear
B$-$C=O bond, which is similar to the bond evolving in configuration 2/6.
$E_{form}$ varies strongly between $-8.4$ eV and $-34.4$ eV due to the
formation of B$-$B dimers. Our results show that there is a strong tendency to
B adsorption at oxidized monovacancies.

We note the appearance of fractional magnetic moments in several cases: 0.1
$\mu_{B}$ in configuration 2/1b and 0.2 $\mu_{B}$ in configurations 1/1b and
3/1a. For the remaining systems the moment is either 0 or 1 $\mu_{B}$, see
Table I, and therefore results from localized dangling bonds. In all systems
with fractional moments there are no contributions of the adsorped O and/or B
atoms to the magnetization. Spin polarization is carried only by C atoms,
where the atoms forming the boundary of the defect contribute less than those
further away. While clean monovacancies in graphene do not give rise to
extended $\pi$-band magnetism, the situation changes under O and B decoration.

Figure 4:  Phonon frequency at the $\Gamma$-point for (a) pristine graphene,
the structures (b) 1/1b, (c) 2/1b, (d) 3/1a, (e) 4/1b, (f) 2/2, (g) 2/4, (h)
2/6, and (i) 2 O atoms adsorbed at a monovacancy.

The calculation of $\Gamma$-point phonons allows us to address the structural
stability after O and B doping. Phonon densities of states are shown as
histograms for the most stable structures (1/1b, 2/1b, 3/1a, 4/1b, 2/2, 2/4,
2/6) in Fig. 4 together with results for pristine graphene and for 2 O atoms
adsorbed at a monovacancy. Note that the reported frequencies refers to the
$\Gamma$-point of the Brillouine zone of our $5\times 5$ supercell and due to
backfolding therefore also include frequencies of other points of the standard
graphene Brillouin zone. This phenomenon is well-known from carbon nanotubes
new7 ; new8 ; new9 ; new10 ; new11 . We find all phonon frequencies to be
positive and therefore conclude that all configurations addressed in Fig. 4
are stable. Pristine graphene shows a two peak phonon spectrum, where the
gross shape is largely maintained under O and B adsorption. This finding is
consistent with recent experimental results of Raman spectroscopy of B doped
graphene new5 . However, there are characteristic differences evident in the
high frequency range beyond 1500 cm-1. Except for configuration 1/1b, the
modes in this energy range are suppressed and the high frequency peak shifts
to the left. This fact cannot be a consequence of O adsorption, since an
oxidized monovacancy without adsorbed B atoms rather comes along with an
enhancement of the high energy modes, see Fig. 4(i). The softening of the G
modes appears to be related to out-of-plane shifts of O atoms and substantial
local distortions in the graphene plane induced by B adsorption.

4\. Conclusion

We have performed first principles calculations to study the structural,
electronic, and magnetic properties of O and B decorated graphene. We have
identified the energetically favorable configurations for a variety of O and B
concentrations and have demonstrated that there exist magnetic solutions with
and without dangling bonds. Since vacancies in graphene are reactive and
permit the incorporation of dopants, our calculations demonstrate that B
doping of oxidized vacancies is a successful approach to induce extended
$\pi$-band magnetism. By controlling the O and B concentrations, it is even
possible to tune the magnetic state. A study of the $\Gamma$-point phonons has
been performed to understand the structural stability of the decorated
monovacancies.

Acknowledgments

We thank KAUST research computing for providing the computational resources
used for this investigation. M. Upadhyay Kahaly thanks SABIC for financial
support. R. Faccio thanks the PEDECIBA, CSIC, and Agencia Nacional de
Investigación e Innovación (ANII) Uruguayan organizations for financial
support.

References

## References

  * (1) Castro Neto AH, Guinea F, Peres NMR, Novoselov KS, Geim AK. The electronic properties of graphene. Rev. Mod. Phys. 2009; 81: 109-162.
  * (2) Makarova T, Palacio F. C-based magnetism: An overview of metal free C-based compounds and materials. Amsterdam: Elsevier. 2005.
  * (3) Esquinazi P, Setzer A, Höhne R, Semmelhack C, Kopelevich Y, Spemann D, et al. Ferromagnetism in oriented graphite samples. Phys. Rev. B 2002; 66: 024429.
  * (4) Esquinazi P, Spemann D, Höhne R, Setzer A, Han KH, T. Butz. Induced magnetic ordering by proton irradiation in graphite. Phys. Rev. Lett. 2003; 91: 227201.
  * (5) Rode AV, Gamaly EG, Christy AG, Gerald JGF, Hyde ST, Elliman RG, et al. Unconventional magnetism in all-carbon nanofoam. Phys. Rev. B 2004; 70: 054407.
  * (6) Lehtinen PO, Foster AS, Ma Y, Krasheninnikov AV, Nieminen RM. Irradiation-induced magnetism in graphite: A density functional study. Phys. Rev. Lett. 2004; 93: 187202.
  * (7) ${\rm\check{C}}$rvenka JS, Katsnelson MI, Flipse CFJ. Room-temperature ferromagnetism in graphite driven by two-dimensional networks of point defects. Nat. Phys. 2009; 5: 840-844.
  * (8) Ugeda MM, Brihuega I, Guinea F, Gómez-Rodríguez JM. Missing atom as a source of carbon magnetism. Phys. Rev. Lett. 2010; 104: 096804.
  * (9) Kaloni TP, Upadhyay Kahaly M, Schwingenschlögl U. Induced magnetism in transition metal intercalated graphitic systems. J. Mater. Chem. 2011; 21: 18681-18685.
  * (10) Singh R, Kroll P. Magnetism in graphene due to single-atom defects: Dependence on the concentration and packing geometry of defects. J. Phys.: Condens. Matter 2009; 21: 196002.
  * (11) Nanda BKR, Sherafati M, Popovic Z, Satpathy S. Electronic structure of the substitutional vacancy in graphene: Density-functional and Green’s function studies. New J. Phys. 2012; 14: 083004.
  * (12) Yang X, Xia H, Qin X, Li W, Dai Y, Liu X et al. Correlation between the vacancy defects and ferromagnetism in graphite. Carbon 2009; 47: 1399.
  * (13) Lim DH, Negreira AS, Wilcox J. DFT Studies on the interaction of defective graphene-supported Fe and Al nanoparticles. J. Phys. Chem. C 2011; 115: 8961.
  * (14) Nair RR, Sepioni M, Tsai IL, Lehtinen O, Keinonen J, Krasheninnikov AV, et al. Spin-half paramagnetism in graphene induced by point defects. Nat. Phys. 2012; 8: 199-202.
  * (15) Palacios JJ, Ynduráin F. Critical analysis of vacancy-induced magnetism in monolayer and bilayer graphene. Phys. Rev. B 2012; 85: 245443.
  * (16) Pardo H, Khane ND, Faccio R, Araújo-Moreira FM, Fernández-Werner L, Makarova T, et al. Raman characterization of bulk ferromagnetic nanostructured graphite. Physica B: Cond. Matt. 2012; 407; 3206-3209.
  * (17) Zagaynova VS, Makarova TL, Bergstrom L, Vasiliev P, Araújo-Moreira FM, Pardo H, et al. Boron doped graphene oxide. Submitted (2012).
  * (18) Faccio R, Fernández-Werner L, Pardo H, Goyenola C, Ventura ON, Mombrú AW. Electronic and structural distortions in graphene induced by carbon vacancies and boron doping. J. Phys. Chem. C 2010; 114: 18961-18971.
  * (19) Kaloni TP, Cheng YC, Faccio R, Schwingenschlögl U. Oxidation of monovacancies in graphene by oxygen molecules. J. Mater. Chem. 2011; 21: 18284-18288.
  * (20) Pardo H, Faccio R, Araujo-Moreira FM, Lima de OF, Mombru WA. Synthesis and characterization of stable room temperature bulk ferromagnetic graphite. Carbon 2006; 44: 565.
  * (21) Giannozzi P, Baroni S, Bonini N, Calandra M, Car R, Cavazzoni C, et al. Quantum ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009; 21: 395502.
  * (22) Cheng YC, Kaloni TP, Zhu ZY, Schwingenschlögl U. Oxidation of graphene in ozone under ultraviolet light. Appl. Phys. Lett. 2012; 101: 073110.
  * (23) Saha KS, Chandrakanth RC, Krishnamurthy HR, Waghmare UV. Mechanisms of molecular doping of graphene: A first-principles study. Phys. Rev. B 2009; 80: 155414.
  * (24) Johari P, Shenoy VB. Modulating optical properties of graphene oxide: Role of prominent functional groups. ACS Nano 2011; 5: 7640.
  * (25) Dai J, Yuan J, Giannozzi P. Gas adsorption on graphene doped with B, N, Al, and S: A theoretical study. Appl. Phys. Lett. 2009; 95: 232105.
  * (26) Yan J-A, XianL, Chou MY. Structural and electronic properties of oxidized graphene. Phys. Rev. Lett. 2009; 103: 086802.
  * (27) Wu M, Cao C, Jiang JZ. Electronic structure of substitutionally Mn-doped graphene. New J. Phys. 2010; 12: 063020.
  * (28) Wu B-R, Yang C-K. Electronic structures of graphane with vacancies and graphene adsorbed with fluorine atoms. AIP Advances 2012; 2: 012173.
  * (29) Baroni S, de Gironcoli S, Dal Corso A. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 2001; 73: 515-562.
  * (30) Jahn A, Teller E. Stability of polyatomic molecules in degenerate electronic states. I. Orbital degeneracy. Proc. R. Soc. Lond. A 1937; 61: 220-235.
  * (31) Banhart F, Kotakoski J, Krasheninnikov AV. Structural defects in graphene. ACS Nano 2011; 5: 26-41.
  * (32) El-Barbary AA, Telling RH, Ewels CP, Heggie MI, Briddon PR. Structure and energetics of the vacancy in graphite. Phys. Rev. B 2003; 68: 144107.
  * (33) Sheng ZH, Gao HJ, Bao WJ, Wang FB, Xia XH. Synthesis of boron doped graphene for reduction reaction in fuel cells. J. Mater. Chem. 2012; 22: 390.
  * (34) Marconcini P and Macucci M. A novel choice of the graphene unit vectors, useful in zone-folding computations. Carbon 2007; 45: 1018.
  * (35) Kahaly MU, Waghmare UV. Vibrational properties of single-wall carbon nanotubes: A first-principles study. J. Nanosci. Nanotechnol. 2007; 7: 1787.
  * (36) Maultzsch J, Reich S, Thomsen C, Dobardzic E, Milosevic I, Damnjanovic M. Phonon dispersion of carbon nanotubes. Solid State Commun. 2002; 121: 471.
  * (37) Saito R, Takeya T, Kimura T, Dresselhaus G, Dresselhaus MS. Raman intensity of single-wall carbon nanotubes. Phys. Rev. B 1998; 57: 4145.
  * (38) Doorn SK, Zheng L, O’Connell MJ, Zhu Y, Huang S, Liu J. Raman spectroscopy and imaging of ultralong carbon nanotubes. Phys. Chem. 2005; 109: 3751.
  * (39) Kim YA, Fujisawa K, Muramatsu H, Hayashi T, Endo M, Fujimori T et al. Raman spectroscopy of boron-doped single-layer graphene. ACS Nano 2012; 6: 6293.


11institutetext: Department of Computer Science, UNC Chapel Hill
11email<EMAIL_ADDRESS>22institutetext: Eshelman School of Pharmacy, UNC
Chapel Hill
22email<EMAIL_ADDRESS>33institutetext: Department of Pharmacology, UNC Chapel Hill
33email<EMAIL_ADDRESS>44institutetext: Joint Department of Biomedical
Engineering at UNC Chapel Hill and NCSU
44email<EMAIL_ADDRESS>
$\dagger$ Authors contributed equally
$\star$ Corresponding authors

# Utilizing Low-Dimensional Molecular Embeddings for Rapid Chemical Similarity
Search

Kathryn E. Kirchoff 11†† James Wellnitz 22†† Joshua E. Hochuli 22†† Travis
Maxfield 22†† Konstantin I. Popov 22 Shawn Gomez 3344⋆⋆ Alexander Tropsha 22⋆⋆

###### Abstract

Nearest neighbor-based similarity searching is a common task in chemistry,
with notable use cases in drug discovery. Yet, some of the most commonly used
approaches for this task still leverage a brute-force approach. In practice
this can be computationally costly and overly time-consuming, due in part to
the sheer size of modern chemical databases. Previous computational
advancements for this task have generally relied on improvements to hardware
or dataset-specific tricks that lack generalizability. Approaches that
leverage lower-complexity searching algorithms remain relatively
underexplored. However, many of these algorithms are approximate solutions
and/or struggle with typical high-dimensional chemical embeddings. Here we
evaluate whether a combination of low-dimensional chemical embeddings and a
k-d tree data structure can achieve fast nearest neighbor queries while
maintaining performance on standard chemical similarity search benchmarks. We
examine different dimensionality reductions of standard chemical embeddings as
well as a learned, structurally-aware embedding—SmallSA—for this task. With
this framework, searches on over one billion chemicals execute in less than a
second on a single CPU core, five orders of magnitude faster than the brute-
force approach. We also demonstrate that SmallSA achieves competitive
performance on chemical similarity benchmarks.

###### Keywords:

Cheminformatics Virtual screening Drug discovery

##  1 Introduction

Searching a database for the nearest neighbors to a query is a task that spans
many fields, from computer graphics [45, 38] to medical diagnostics [18] and
chemistry [32]. In chemistry, the identification and retrieval of similar
compounds plays a vital role in various domains, including drug discovery.
Nearest neighbor search algorithms have been widely employed for this purpose,
allowing researchers to explore large chemical databases to help discover
potential drug candidates. When applied to this task, nearest neighbor
searching, or chemical similarity searching, is often called virtual
screening, where databases are searched for chemicals similar to ones with
known desirable properties, assuming that similar chemicals share similar
properties [13]. This approach is widely adopted and has been used
successfully to discover potent drug candidates [16]. Currently, most chemical
similarity searching methods, like the popular Arthor [2], use a simple,
brute-force algorithm to compare all queries to all database chemicals.

However, the sizes of searchable chemical databases have recently undergone
dramatic expansion. For example, the latest release of the ZINC database [42]
and Enamine’s catalog [17] now number close to 40 billion chemicals each.
Brute-force approaches to nearest neighbor searching struggle to scale to
databases of this size. Often, hundreds of queries need to be searched against
the database, resulting in trillions of calculations. As a result, there is a
growing need for faster search algorithms that can handle the ever-increasing
data volume.

Figure 1: Overview of the similarity search framework. k-d trees are combined
with with low-dimensional chemical embeddings to produce a partitioned
chemical space, which can be quickly queried for nearest neighbors.

Indeed, there exist algorithms that effect faster runtimes for brute-force
searching, one being the k-dimensional tree (k-d tree), which achieves
logarithmic searching complexity [10]. To date, k-d trees have been used to
achieve algorithmic speed-ups in numerous applications. However, this data
structure has yet to be effectively used for the cheminformatics tasks of
large-scale chemical similarity search and virtual screening. This lack of use
could stem from two issues that arise when utilizing k-d trees for these
purposes. First, chemical databases contain tens of billions of entries, and a
k-d tree built on such a database of such size cannot exist purely in memory
on most systems, thus requiring implementations designed for low-memory usage.
Second, k-d trees are only effective at searching low-dimensional data and
provide near-zero benefit as the data dimensionality exceeds about 20 [48], a
problem since nearly all common chemical representations exceed 20 dimensions.

In fact, typical representations of chemicals have dimensionalities in the
thousands [52], which are nowhere near compatible with the k-d tree approach.
There is precedent for reducing high-dimensional objects to lower dimensions
while preserving relative distance [20], leveraging techniques such as random
projection [5] or Principal Components Analysis (PCA) [21]; however, to our
knowledge dimensionality reduction techniques have not yet been applied in
conjunction with k-d trees to allow for faster chemical similarity searching.
Specific to chemistry, chemical representation learning has emerged as a
powerful technique to organize chemicals in a meaningful way at lower
dimensions [51]. However, previous implementations of these techniques still
do not reduce dimensionality low enough, nor are these approaches always
trained to organize the embedding space in a meaningful way for the specific
task of chemical similarity search.

In this work, we propose an effective framework for chemical similarity
searching: combining a meaningful low-dimensional chemical embedding with a
custom k-d tree implementation designed to index billion-sized chemical data
sets with low memory constraints. Further, we show that utilizing
representation learning to generate Small Structurally-Aware embeddings
(SmallSA) provides better performance on similarity searching benchmarks
compared to dimensional reduction methods applied to existing high-dimensional
embeddings.

To our knowledge, this is the first publicly disseminated method that
leverages a k-d tree with low-dimensional chemical embeddings for the task of
rapid chemical similarity searching with applications to drug discovery.

##  2 Background

###  2.1 Nearest Neighbor Searching Algorithms

Nearest neighbor search has been the subject of many algorithmic improvements
aimed at lowering complexity. For our purposes, we will distinguish two main
classes of such improvements: 1) reducing searchable space and 2)
dimensionality reduction techniques.

####  2.1.1 Reducing Searchable Space

##### Indexing Methods.

Given a search database, indexing methods, like the BallTree [33] and
k-dimensional tree (k-d tree) [10], partition it in such a way that individual
queries need not be compared to every element in the database. This
effectively reduces the search space, leading to at best $O(q\log{}n)$
complexity, where $q$ is the number of queries and $n$ is the size of the
space. However, this best-case speed-up only occurs for sufficiently small
data dimension $d$, [43] collapsing to the brute-force complexity of $O(qn)$
beyond about $d=20$ [48].

##### Approximate Nearest Neighbor Search.

Indexing methods are exact, but if absolute correctness is not required there
are also fast approximate nearest neighbor methods [46]. These methods often
involve the use of hashing, as in the case of locality sensitive hashing (LSH)
forests [7] and spectral hashing [50], or the use of clustering [14] to group
together similar portions of the database. Like the indexing methods described
above, these approximate search methods can be described as reducing the
effective search database size, $n$, resulting in sublinear search times.
However, the accuracy of these methods is dependent on the quality of the
grouping, whether hashing clustering or otherwise [46]. Though fast algorithms
exist, effective clustering methods can also have difficulty scaling to
massive data sets [29], bringing into question their utility on billion-sized
chemical datasets.

####  2.1.2 Dimensionality Reduction.

Dimensionality reduction methods, like PCA [21], t-SNE [44] or UMAP [30], can
also be applied to increase search speed by reducing the number of dimensions,
$d$. Accounting for dimensionality, the complexity of brute-force search can
be described as $O(qdn)$, and thus a smaller $d$ can increase efficiency.
Further, if $d$ is reduced to a small enough dimensionality, it can be applied
in tandem with indexing methods, yielding sublinear complexity [6, 23, 25].
Like the approximate methods described above, accuracy is tied to the ability
of the dimensionality reduction to preserve the original spatial relationships
between data points, and not all dimensionality reduction methods can scale to
billion-sized data sets.

###  2.2 Chemical Similarity

In cheminformatics, “chemical similarity searching” amounts to searching a
database for the nearest $n$ neighbors to a given query molecule. “Similarity”
in chemical space is not a well-defined concept, and there is no agreed-upon
definition as to what constitutes a pair of similar chemicals [28]. For
example, a chemist can claim two chemicals are similar due to similar
structure, while another can claim the pair is dissimilar due to different
properties. With this in mind, and to establish a baseline for comparison, we
choose a broadly applicable definition of chemical similarity: two chemicals
are similar if their structures are. In this case, graph edit distance (GED)
can be used as the ground truth metric for chemical similarity.

The graph edit distance (GED) [37] between two graphs is the smallest total
number of edits to one graph that are needed to make it isomorphic to the
second graph. GED is often approximated, as calculating the actual value is
NP-hard [53]. Since a chemical structure can be represented as a graph of
atoms (nodes) and bonds (edges), the distance between any two chemicals can be
measured with approximate GED [11].

###  2.3 Virtual Screening

Chemical similarity searching is often used in early stages of drug discovery,
in a process called “virtual screening” [26]. This practice is based on the
idea that chemicals with similar structures likely share similar biological
activity profiles, a concept known as structure-activity relationship (SAR)
[13]. Given a single chemical with a desirable biological activity, similarity
search can be used to find similar chemicals that are expected to share that
activity profile. Generally, virtual screening is carried out with hundreds of
query molecules against a database of billions. These large chemical databases
are projected to get larger, increasing into the trillions [27].

###  2.4 Related Works

There are some existing chemical similarity search methods designed to
function on large chemical databases, notably Arthor [2, 39] and SmallWorld
[4], both developed by NextMove Software, and SpaceLight [8]. Arthor uses a
brute-force search, but is accelerated through a highly optimized, domain-
specific implementation [40] that, while fast, is still linear in complexity.
Further, it utilizes a hashing technique that is vulnerable to collisions
which can hinder the accuracy. SmallWorld’s approach is based on GEDs between
chemicals, generating a graph of graphs representing the database, where
chemical graphs are connected by an edge if they are $1$ GED apart; given a
graph in the database, the exact GED to nearby graphs can then be quickly
calculated. However, this method requires expensive, high-end hardware to run,
and can only be used in cases where the search query is in the database since
it relies on an extremely expensive pre-calculation of the database.
Spacelight uses an approximate clustering approach that relies on the search
database being composed of a smaller number of chemical building blocks
combined via known chemical reactions, allowing large parts of the database
being searched to be skipped. Not all chemical databases satisfy this
constraint and, even among those that do, the required decomposition is often
proprietary, limiting the utility. Furthermore, all three methods described
here—SmallWorld, Arthor, and Spacelight—are commercial and not open source.

Among open source tools for chemical similarity searching, there are few
options. ChemMine [12] is a web tool that implements several similarity search
approaches. The most popular, FPSim2 [15], uses the same Baldi bounds [40]
approach as Arthor to partition the database and accelerate brute-force
search. This partitioning scheme can only generate a maximum number of
partitions equal to the dimensionality of the embedding. Additionally the size
of the partitions are unequal, with the largest partitions also likely to be
those that must be searched for an average query. As a results, this approach
is both limited to binary embeddings and suffers from longer run times on
large chemical databases. There are also open source libraries for fast
approximate neighbor search of high dimensional data, like FAISS (Facebook AI
Similarity Search) [19] or Annoy (Approximate Nearest Neighbors Oh Yeah) [1].
However, these implementations struggle on billion-sized databases and require
significant modification to accommodate current large chemical databases on
minimal hardware. These libraries often rely on some form of dimensionality
reduction paired with a partitioning data structure, making their approach
similar to the approach implemented in this work.

##  3 Methods

###  3.1 Overview of Approach

To efficiently search a large database, we encode chemicals into a low-
dimensional embeddings which we then use to construct a k-d tree. Database
embeddings are calculated once at tree construction time, while query
embeddings are calculated on the fly. Query embeddings are then used to search
the k-d tree to find the $k$ nearest neighbors. This process is summarized in
Figure 1.

###  3.2 SmallSA: Learned Embedding

To generate the learned SmallSA embeddings, we opt to utilize the SALSA
framework presented in [24], which is trained to map chemical compounds to an
embedding space that respects GED. We modify the SALSA architecture to our
specific objective, training a model that compresses the chemical space to
dimensionalities below $20$, specifically to $16$ and $8$ dimensions. As a
further deviation from the original SALSA, we train our model on a diverse
sample of $1.4$ million chemical compounds from the Enamine chemical catalog
[17]. Further details on our implementation of the SALSA approach can be found
in appendix 0.A.

###  3.3 Other Embeddings

There are several commonly utilized approaches for embedding chemicals. The
most common embeddings are Extended-Connectivity Fingerprints (ECFPs) and
Extended-Connectivity Count Vectors (ECFCs) [36]. Both embeddings are
topological fingerprints indicating presence of molecular substructures;
however, ECFPs are _binary_ indications of substructures, while ECFCs indicate
_counts_ of substructures. For both fingerprinting methods, we opted for a
radius of two and a resulting 256d vector. Herein, we refer to these versions
as “ECFP-256” and “ECFC-256”. Furthermore, to calculate embedding distances,
we used the Jaccard index (i.e., Tanimoto index in cheminformatics) [41] for
ECFP-256 and Euclidean distance for ECFC-256.

Additionally, we obtained low-dimensional versions of either fingerprint
through two dimensionality reduction methods: (1) sparse random projections
(SparseRP) and (2) PCA. We computed PCA at both eight and 16 dimensions
(“PCA-8” and “PCA-16”), and computed SparseRP at only 16 dimensions
(“SparseRP-16”). Other common dimensionality reduction methods, including
t-Distributed Stochastic Neighbor Embedding (t-SNE) [44] and Uniform Manifold
Approximation and Projection (UMAP) [30], were initially considered but left
out due to incompatibility with large datasets.

###  3.4 k-d Tree Implementation

As available open-source k-d tree implementations are not designed to handle
the scale of large chemical databases, we implement a custom k-d tree based on
a proposal from [9] designed to handle large datasets. Our resulting k-d tree
is able to organize terabytes of chemical data, but requires only tens of
gigabytes of memory to search. Implementation details and source code are
provided at the GitHub link at the end of the document.

###  3.5 Large Benchmarking Dataset

To benchmark the performance of our proposed framework for chemical similarity
searching on billion-sized chemical databases. A database of $1.3$ billion
chemicals randomly sampled from the $33.5$ billion Enamine REAL Space library
[17] was prepared. Additionally, we sampled $100$ random compounds from the
REAL space to serve as the query set. There is no overlap between the training
set used for SmallSA, mentioned previously, and the database or query set.

##  4 Experiments and Analysis

To assess the effectiveness of our proposed framework for chemical similarity
search—utilizing k-d trees with low-dimensional embeddings—we compare
performance of low (SmallSA-8, SmallSA-16, PCA-8, PCA-16, SparseRP-16) and
high (ECFP-256, ECFC-256) dimensional embeddings on three tasks: 1) GED of
retrieved queries to their retrieved hits, 2) performance on the RDKit Virtual
Screening benchmark, and 3) computational speed (and hit quality) with respect
to computing resources. Results of these experiments are summarized in Table
1.

###  4.1 Query–Hit Graph Edit Distance

Figure 2: The average approximate graph edit distance (GED) between query
molecule and nearest neighbors, per method, shown as a function of the number
of neighbors considered. Lower distances are better. Lines are smoothed using
a running average approach for simplicity of analysis.

To assess the quality of retrieved compounds, we evaluate the GEDs of the top
$k$ retrieved compounds ($k$ nearest neighbors) per query, positing that high
quality compounds will be those of low GED from the query. For each method, we
query the same set of $100$ randomly sampled compounds against the 1.3 billion
dataset, and for each query, determine its $100$ nearest neighbors. GED was
approximated (due to expensive exact computation) using a bipartite graph
matching approach [34] with the cost for all edge and node edits set to one.
We plot the resulting GEDs, per k neighbors, for each method in Figure 2.

Results indicate that, for relatively low $k$, six out of seven embedding
techniques are comparably successful at retrieving high quality neighbors. The
exception is SparseRP-16, which generally performed worse, obtaining neighbors
of higher GED compared to all other methods. Overall, the structurally aware
learned embedding, in this case SmallSA-16, performed the best out of all
embeddings. Further, SmallSA-8 performed better than PCA-8, and comparably to
the best high-dimensional method, ECFC-256.

###  4.2 RDKit Virtual Screening Benchmark

To assess performance at the virtual screening task, we implement the RDKit
Virtual Screening benchmark proposed by Riniker and Landrum [35]. Briefly,
this platform works by determining the performance of each similarity
searching method in "simulated" virtual screening against $69$ proteins
(targets) that are targeted by chemicals. Each simulated screening starts with
a set of chemicals with known activity against the target ("actives") and a
set of chemicals with no activity against the target ("decoys"). A random
subset of actives and decoys are used as the database to search against, while
the remainder of actives are used as the query.

Each compound in the database is ranked according to its closest distance to
any of the actives, and this rank is used as a simple nearest neighbor
classification model. The area under the receiver operating curve (AUROC) for
this model can be calculated using the known activity class of the compounds
in the search database. This is related to how similarity search is used in
practice for drug discovery with the key difference being that the true
activity class is not known a priori, as it is in this benchmark. We plot the
resulting AUROCs, per database, for each method in Figure 3.

The results show that SmallSA-based methods (SmallSA-8 and SmallSA-16)
significantly outperform both the high-dimensional and low-dimensional
embeddings. Considered alongside the GED results, this indicates that there
exist low-dimensional chemical embeddings capable of performing at least
equivalently on the virtual screening task compared to the more traditional
high-dimensional embeddings, like ECFC-256. It should be noted that these
results are not meant to imply that SmallSA-based methods are state-of-the-art
among learned representations; indeed, demonstrating this would require a much
broader comparison with other learned representations. Instead, these results
are meant to simply demonstrate that low-dimensional chemical embeddings are
capable of performance on par with some commonplace higher-dimensional
embeddings, implying that similarity search with these small embeddings can be
useful for downstream drug discovery tasks.

Figure 3: AUROC achieved by each embedding on the RDKit virtual screening
benchmark of 69 query targets, grouped by target database. Each database is
indicated on the x-axis. Note that out of the 69 targets, most targets (50)
belong to the ChEMBL database.

###  4.3 Computational Speed

A major goal of our proposed framework is rapid determination of nearest
neighbors utilizing minimal computing resources. To compare speed performance
to other methods, we collected timings for retrieval of the top $100$
neighbors for $100$ queries on the $1.3$ billion database. Queries using low
dimensional embeddings were executed using both the k-d tree framework and
brute force approach, but only brute force for high-dimensional embeddings.

Nearest neighbor queries were executed on a workstation with three consumer-
grade solid states drives in a RAID 5 array and one core of an AMD
Threadripper PRO 3955WX. Due to high computational cost, brute force for high-
dimensional embeddings were carried out on the UNC Longleaf computing cluster.
Timings are reported for single-core performance in Table 1.

Results indicate that combination of low-dimensional embeddings with k-d trees
achieves substantially reduced compute time—speed-ups between $10^{3}$
(16-dimensional) to $10^{5}$ (8-dimensional) fold—compared to the combination
of high dimensional embeddings with a brute-force approach. Practically, this
results in reduced query times, from several hours down to under a minute with
no change in hardware. Further, upon visual inspection, we find that this
substantial speedup does not coincide with a reduction in perceived hit
quality (see example in Figure 4).

We are unable to compare these timings and benchmarks rigorously to the
commercial releases of Arthor, SmallWorld and SpaceLight due to lack of
access; however all methods report run-times ranging between several seconds
to several minutes, depending on the query [8, 3].

Table 1: GED, virtual screening, and timing results. Best is bold, second best is underlined. Runtime is proportional to the dimensionality, not method of dimensionality reduction, therefore methods with the same dimension have identical query times. |  | Query accuracy | Query time (s)
---|---|---|---
Method |  D |  VS AUC |  GED |  Brute-force |  k-d tree
ECFC-256 |  256 |  0.494 |  16.69 |  $24,900\pm 2,200$ |  N/A
ECFP-256 |  |  0.566 |  18.30 |  $13,600\pm{}1,000$
SparseRP-16 |  16 |  0.501 |  15.06 |  $13,000\pm{}70$ |  $37.5\pm{}26$
PCA-16 |  |  0.524 |  16.97
SmallSA-16 |  |  0.692 |  15.74
PCA-8 |  8 |  0.517 |  18.78 |  11,400 ± 98 |  0.606 ± 0.35
SmallSA-8 |  |  0.662 |  17.30*
Figure 4: Example of a query molecule and the hits obtained by select high-
performing embeddings.

##  5 Discussion

In this work, we showcase the feasibility of achieving rapid speed ups in
similarity search, through a novel application combining low-dimensional
embeddings and an efficient, low-memory k-d tree implementation. Further, we
show that utilizing a learned, structurally-aware embedding approach (SmallSA)
effects a low-dimensional embedding that excels at virtual screening tasks,
demonstrating marked improvements over both traditionally high-dimensional
fingerprints as well as their dimension-reduced counterparts.

Alongside our work, others have explored low-dimensional [22, 31] and
structure-aware [24, 47] chemical representations for statistical modeling and
faster, brute-force similarity searching. There is also extensive literature
on methods for sublinear nearest neighbor searching. Ours is the first work,
to our knowledge, to show that combining low-dimensional embeddings and k-d
trees can achieve sublinear chemical similarity searches on billion-sized
datasets while still maintaining sensible results. We do not contend that the
combination specifically of SmallSA and k-d trees is the most optimal setup,
as we did not investigate other learned, low-dimensional embeddings or other
partition-tree based similarity searching algorithms. Instead we point out
that leveraging a meaningful, chemically aware, low-dimensional embedding with
a searching algorithm that excels with low-dimensional data can provide major
speed ups in exact chemical similarity searching.

###  5.1 Requirement for Exactness

Many of the open-source libraries for similarity searching of large databases
leverage approximate searching algorithms. This means they can fail to return
the true nearest neighbor. This is a concern for chemists, as often times the
goal of a chemical similarity search is to find a small number of compounds
most similar to a query for purchase and experimental assessment. Thus,
assuming a good structurally-aligned representation of chemicals, a similarity
searching tool should be designed to prioritize exactness. Leveraging an exact
implementation from available open source methods for high-dimensional
similarity searching often results in either brute force searches (e.g.
’IndexFlatL2’ in FAISS [19]) or dramatically slower run time. Instead, we
chose to implement a k-d tree, which provides exact performance while
providing an efficient sublinear search, assuming dimensionality of embeddings
is below $20$.

###  5.2 Additional Applications

There are additional cheminformatics tasks that can be accelerated using the
framework reported here. The k-d tree can execute range queries to return all
chemicals within some bounds in embedding space. Quantitative structure-
activity relationship (QSAR) models often use tree-based algorithms (e.g.
decision trees) for predicting the biological activity of chemicals. A learned
tree-based model built on the same low-dimensional embeddings used in this
work can then encode the regions of chemical space where desirable chemicals
exist. Thus, range queries on a k-d tree can quickly filter out predicted
active compounds from billion-sized chemical databases according to a learned
model. Fine-tuning of other statistical models, like deep neural networks,
could also make them compatible with this approach if they are built or
constricted to partition a space. Such an approach would allow QSAR models to
be more efficiently utilized for screening of massive chemical databases.

#### Data and Code Availability.

Source code and dataset information can be found on GitHub.

#### Acknowledgments.

We would like to thank the Research Computing group at the University of North
Carolina at Chapel Hill for providing computational resources and support.
This work was supported by the National Institutes of Health (Grant
R01GM140154). JW is supported by the National Institute of General Medical
Sciences of the National Institutes of Health under Award Number T32GM135122.
This work was supported by the following grants to SMG from the National
Institutes of Health - U24DK116204 (NIDDK), U01CA238475 (NCI), R01CA233811
(NCI), U01CA274298 (NCI). T.M. was supported by the National Institute of
General Medical Sciences of the NIH under Award Number T32GM086330. The
funders had no role in study design, data collection and analysis, decision to
publish, or preparation of the manuscript.

## Appendix

##  Appendix 0.A SALSA Embedding Procedure

In learning molecular representations, chemicals are often initially expressed
as alphanumeric sequences, so-called SMILES strings [49] which effectively
describe molecular graphs (that is, atoms connected by bonds). SMILES strings
enable the use of sequence-to-sequence autoencoders for learning molecular
representations. However, autoencoders trained solely on SMILES are
insufficient to learn molecular representations that capture the structural
(graph-to-graph) similarities between molecules, resulting in disorganized
embeddings that potentially hinder downstream applications.

The recently published molecular representation model, SALSA [24], seeks to
remedy this limitation of SMILES-based autoencoders by explicitly enforcing
structural awareness onto learned representations. This is accomplished by
modifying a SMILES-based autoencoder to additionally learn a contrastive
objective of mapping structurally similar molecules to similar
representations. Here, “structurally similar” is defined as any two molecules
having a graph edit distance (GED) of one. This contrastive objective
necessitates a dataset comprised of pairs of 1-GED molecules. For their
implementation, the SALSA authors chose to use a curated chemical dataset from
ChEMBL, and then, from this dataset, generated 1-GED neighbor compounds (or,
mutants) for each ChEMBL anchor compound.

The standard operational definition of GED permits a wide array of
modifications that may drastically change chemical and biological properties
of small molecules. The SALSA methodology operates on a more conservative
notion of graph edits aimed at maintaining relative chemical similarity
between an anchor compound and mutant compounds. To obtain 1-GED mutants,
SALSA defines three node-level (i.e. atom-level) transformations:

  1. 1.

Addition: Append a new node and corresponding edge to an existing node

  2. 2.

Substitution: Change the atom type of an existing node

  3. 3.

Deletion: Remove a singly-attached node and its corresponding edge

Notably, these transformations neither break nor create new cycles in the
molecular graphs, and they do not disconnect previously connected graphs.
There are additional nuances to the definition of graph edits, and we refer
readers to the referenced publication for full details.

We extend the SALSA framework starting with a dataset of approximately 1.3
million chemicals diversity-sampled from the Enamine REAL library of
approximately 40 billion chemicals. Additionally, we modified the
hyperparameters of the model to output a chemical embedding whose
dimensionality would enable the fast indexing search methods described
previously. Specifically, we chose embeddings of dimensions $8$ and $16$.

## References

  * [1] ANNOY library. https://github.com/spotify/annoy, accessed: 2017-08-01
  * [2] NextMove software | Arthor, https://www.nextmovesoftware.com/arthor.html
  * [3] NextMove software | Arthor, {https://www.nextmovesoftware.com/talks/Sayle_EvolutionVsRevolution_ICCS_201805.pdf}
  * [4] NextMove software | SmallWorld, https://www.nextmovesoftware.com/smallworld.html
  * [5] Achlioptas, D.: Database-friendly random projections: Johnson-lindenstrauss with binary coins 66(4), 671–687. https://doi.org/10.1016/S0022-0000(03)00025-4
  * [6] Agrawal, R., Faloutsos, C., Swami, A.: Efficient similarity search in sequence databases pp. 69–84. https://doi.org/10.1007/3-540-57301-1_5
  * [7] Bawa, M., Condie, T., Ganesan, P.: LSH forest: self-tuning indexes for similarity search pp. 651–660. https://doi.org/10.1145/1060745.1060840, https://doi.org/10.1145/1060745.1060840
  * [8] Bellmann, L., Penner, P., Rarey, M.: Topological similarity search in large combinatorial fragment spaces 61(1), 238–251. https://doi.org/10.1021/acs.jcim.0c00850, https://doi.org/10.1021/acs.jcim.0c00850, publisher: American Chemical Society
  * [9] Bentley, J.: Multidimensional Binary Search Trees in Database Applications. IEEE Transactions on Software Engineering SE-5(4), 333–340 (Jul 1979). https://doi.org/10.1109/TSE.1979.234200
  * [10] Bentley, J.L.: Multidimensional binary search trees used for associative searching. Communications of the ACM 18(9), 509–517 (Sep 1975). https://doi.org/10.1145/361002.361007
  * [11] Birchall, K., Gillet, V.J., Harper, G., Pickett, S.D.: Training similarity measures for specific activities: Application to reduced graphs 46(2), 577–586. https://doi.org/10.1021/ci050465e, https://doi.org/10.1021/ci050465e, publisher: American Chemical Society
  * [12] Cao, Y., Jiang, T., Girke, T.: Accelerated similarity searching and clustering of large compound sets by geometric embedding and locality sensitive hashing. Bioinformatics 26(7), 953–959 (Feb 2010). https://doi.org/10.1093/bioinformatics/btq067, http://dx.doi.org/10.1093/bioinformatics/btq067
  * [13] Crum-Brown, A., Fraser, T.R.: The connection of chemical constitution and physiological action. Trans R Soc Edinb 25(1968-1969), 257 (1865)
  * [14] Deng, Z., Zhu, X., Cheng, D., Zong, M., Zhang, S.: Efficient kNN classification algorithm for big data 195, 143–148. https://doi.org/10.1016/j.neucom.2015.08.112, https://www.sciencedirect.com/science/article/pii/S0925231216001132
  * [15] Félix, E., Dalke, A., Landrum, G., Bushuiev, R.: chembl/fpsim2: 0.5.0 (2023). https://doi.org/10.5281/ZENODO.10041218, https://zenodo.org/doi/10.5281/zenodo.10041218
  * [16] Gahbauer, S., Correy, G.J., Schuller, M., Ferla, M.P., Doruk, Y.U., Rachman, M., Wu, T., Diolaiti, M., Wang, S., Neitz, R.J., Fearon, D., Radchenko, D.S., Moroz, Y.S., Irwin, J.J., Renslo, A.R., Taylor, J.C., Gestwicki, J.E., von Delft, F., Ashworth, A., Ahel, I., Shoichet, B.K., Fraser, J.S.: Iterative computational design and crystallographic screening identifies potent inhibitors targeting the nsp3 macrodomain of sars-cov-2. Proceedings of the National Academy of Sciences 120(2), e2212931120 (2023). https://doi.org/10.1073/pnas.2212931120, https://www.pnas.org/doi/abs/10.1073/pnas.2212931120
  * [17] Grygorenko, O.O.: Enamine Ltd.: The Science and Business of Organic Chemistry and Beyond. European Journal of Organic Chemistry 2021(47), 6474–6477 (Dec 2021). https://doi.org/10.1002/ejoc.202101210, https://onlinelibrary.wiley.com/doi/10.1002/ejoc.202101210
  * [18] Gupta, D., Loane, R., Gayen, S., Demner-Fushman, D.: Medical image retrieval via nearest neighbor search on pre-trained image features (arXiv:2210.02401). https://doi.org/10.48550/arXiv.2210.02401, http://arxiv.org/abs/2210.02401, version: 1
  * [19] Johnson, J., Douze, M., Jégou, H.: Billion-scale similarity search with GPUs. IEEE Transactions on Big Data 7(3), 535–547 (2019)
  * [20] Johnson, W.B., Lindenstrauss, J., Schechtman, G.: Extensions of lipschitz maps into banach spaces 54(2), 129–138. https://doi.org/10.1007/BF02764938, https://doi.org/10.1007/BF02764938
  * [21] Jolliffe, I.: Principal Component Analysis. In: Balakrishnan, N., Colton, T., Everitt, B., Piegorsch, W., Ruggeri, F., Teugels, J.L. (eds.) Wiley StatsRef: Statistics Reference Online. Wiley, first edn. (Sep 2014). https://doi.org/10.1002/9781118445112.stat06472
  * [22] Karlova, A., Dehaen, W., Svozil, D.: Molecular Fingerprint VAE p. 6
  * [23] Keogh, E., Chakrabarti, K., Pazzani, M., Mehrotra, S.: Dimensionality reduction for fast similarity search in large time series databases 3(3), 263–286. https://doi.org/10.1007/PL00011669, https://doi.org/10.1007/PL00011669
  * [24] Kirchoff, K.E., Maxfield, T., Tropsha, A., Gomez, S.M.: Salsa: Semantically-aware latent space autoencoder (2023)
  * [25] Korn, P., Sidiropoulos, N.D., Faloutsos, C., Siegel, E.L., Protopapas, Z.: Fast and effective similarity search in medical tumor databases using morphology: Multimedia storage and archiving systems 2916, 116–129. https://doi.org/10.1117/12.257282
  * [26] Lavecchia, A., Giovanni, C.D.: Virtual screening strategies in drug discovery: A critical review 20(23), 2839–2860, https://www.eurekaselect.com/article/53238
  * [27] Lyu, J., Irwin, J.J., Shoichet, B.K.: Modeling the expansion of virtual screening libraries 19(6), 712–718. https://doi.org/10.1038/s41589-022-01234-w, https://www.nature.com/articles/s41589-022-01234-w, number: 6 Publisher: Nature Publishing Group
  * [28] Maggiora, G., Vogt, M., Stumpfe, D., Bajorath, J.: Molecular Similarity in Medicinal Chemistry: Miniperspective. Journal of Medicinal Chemistry 57(8), 3186–3204 (Apr 2014). https://doi.org/10.1021/jm401411z
  * [29] Mahdi, M.A., Hosny, K.M., Elhenawy, I.: Scalable clustering algorithms for big data: A review 9, 80015–80027. https://doi.org/10.1109/ACCESS.2021.3084057, conference Name: IEEE Access
  * [30] McInnes, L., Healy, J., Melville, J.: UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction (2018). https://doi.org/10.48550/ARXIV.1802.03426
  * [31] Nasser, M., Salim, N., Saeed, F., Basurra, S., Rabiu, I., Hamza, H., Alsoufi, M.A.: Feature Reduction for Molecular Similarity Searching Based on Autoencoder Deep Learning. Biomolecules 12(4), 508 (Mar 2022). https://doi.org/10.3390/biom12040508
  * [32] Nikolova, N., Jaworska, J.: Approaches to Measure Chemical Similarity - a Review. QSAR & Combinatorial Science 22(9-10), 1006–1026 (Dec 2003). https://doi.org/10.1002/qsar.200330831
  * [33] Omohundro, S.: Five Balltree Construction Algorithms. Tech. Rep. TR-89-063, International Computer Science Institute (1989)
  * [34] Riesen, K., Bunke, H.: Approximate graph edit distance computation by means of bipartite graph matching 27(7), 950–959. https://doi.org/10.1016/j.imavis.2008.04.004, https://www.sciencedirect.com/science/article/pii/S026288560800084X
  * [35] Riniker, S., Landrum, G.A.: Open-source platform to benchmark fingerprints for ligand-based virtual screening. J. Cheminform. 5(1), 26 (May 2013)
  * [36] Rogers, D., Hahn, M.: Extended-Connectivity Fingerprints. Journal of Chemical Information and Modeling 50(5), 742–754 (May 2010). https://doi.org/10.1021/ci100050t
  * [37] Sanfeliu, A., Fu, K.S.: A distance measure between attributed relational graphs for pattern recognition SMC-13(3), 353–362. https://doi.org/10.1109/TSMC.1983.6313167, conference Name: IEEE Transactions on Systems, Man, and Cybernetics
  * [38] Sankaranarayanan, J., Samet, H., Varshney, A.: A fast all nearest neighbor algorithm for applications involving large point-clouds. Computers & Graphics 31(2), 157–174 (Apr 2007). https://doi.org/10.1016/j.cag.2006.11.011
  * [39] Swamidass, S.J., Baldi, P.: Bounds and algorithms for fast exact searches of chemical fingerprints in linear and sublinear time 47(2), 302–317. https://doi.org/10.1021/ci600358f, https://pubs.acs.org/doi/10.1021/ci600358f
  * [40] Swamidass, S.J., Baldi, P.: Bounds and algorithms for fast exact searches of chemical fingerprints in linear and sublinear time. Journal of Chemical Information and Modeling 47(2), 302–317 (2007). https://doi.org/10.1021/ci600358f, https://doi.org/10.1021/ci600358f, pMID: 17326616
  * [41] Tanimoto, T.T.: An elementary mathematical theory of classification and prediction by T.T. Tanimoto. International Business Machines Corporation New York (1958)
  * [42] Tingle, B.I., Tang, K.G., Castanon, M., Gutierrez, J.J., Khurelbaatar, M., Dandarchuluun, C., Moroz, Y.S., Irwin, J.J.: ZINC-22-A Free Multi-Billion-Scale Database of Tangible Compounds for Ligand Discovery. Journal of Chemical Information and Modeling 63(4), 1166–1176 (Feb 2023). https://doi.org/10.1021/acs.jcim.2c01253
  * [43] Toth, C.D., O’Rourke, J., Goodman, J.E.: Handbook of Discrete and Computational Geometry. CRC Press, google-Books-ID: 9mlQDwAAQBAJ
  * [44] van der Maaten, L., Hinton, G.: Visualizing data using t-SNE. Journal of Machine Learning Research 9(86), 2579–2605 (2008)
  * [45] Wald, I., Havran, V.: On building fast kd-Trees for Ray Tracing, and on doing that in O(N log N). In: 2006 IEEE Symposium on Interactive Ray Tracing. pp. 61–69. IEEE, Salt Lake City, UT, USA (Sep 2006). https://doi.org/10.1109/RT.2006.280216
  * [46] Wang, J., Shen, H.T., Song, J., Ji, J.: Hashing for similarity search: A survey (arXiv:1408.2927). https://doi.org/10.48550/arXiv.1408.2927, http://arxiv.org/abs/1408.2927
  * [47] Wang, Y., Wang, J., Cao, Z., Farimani, A.B.: Molclr: Molecular contrastive learning of representations via graph neural networks. CoRR abs/2102.10056 (2021), https://arxiv.org/abs/2102.10056
  * [48] Weber, R., Schek, H.J., Blott, S.: Analysis and Performance Study for Methods in High-Dimensional Spaces p. 12
  * [49] Weininger, D.: SMILES, a chemical language and information system. 1. introduction to methodology and encoding rules 28(1), 31–36. https://doi.org/10.1021/ci00057a005, https://doi.org/10.1021/ci00057a005, publisher: American Chemical Society
  * [50] Weiss, Y., Torralba, A., Fergus, R.: Spectral hashing 21
  * [51] Wigh, D.S., Goodman, J.M., Lapkin, A.A.: A review of molecular representation in the age of machine learning. WIREs Computational Molecular Science 12(5) (Sep 2022). https://doi.org/10.1002/wcms.1603
  * [52] Yang, J., Cai, Y., Zhao, K., Xie, H., Chen, X.: Concepts and applications of chemical fingerprint for hit and lead screening 27(11), 103356. https://doi.org/10.1016/j.drudis.2022.103356
  * [53] Zeng, Z., Tung, A.K.H., Wang, J., Feng, J., Zhou, L.: Comparing stars: on approximating graph edit distance 2(1), 25–36. https://doi.org/10.14778/1687627.1687631, https://doi.org/10.14778/1687627.1687631


# Observation of mode-mixing in the eigenmodes of an optical microcavity:
supplemental document

###### Abstract

This supplement contains four sections: (i) it presents equations for the
roundtrip operator $M$, (ii) it links the evolution of the intracavity field
to the dynamics of the coupled cavity modes, (iii) it presents the analysis of
AFM measurements of the mirror shape and links this to the shape-induced mode
coupling, and (iv) it presents measurements on another cavity to demonstrate
the generic behaviour of mode coupling.

##  1 Roundtrip operator $M$

The evolution of the optical field in a cavity can be described by the
roundtrip operator $M$, which transforms the forward-propagating intra-cavity
field $|\psi^{+}\rangle$ into $M|\psi^{+}\rangle$ after a roundtrip [1]. If we
expand the field $|\psi^{+}\rangle=\sum_{i}c_{i}|\psi_{i}^{+}\rangle$ in a set
of orthogonal basis states $|\psi^{+}_{i}\rangle$, we can write
$|\psi^{+}\rangle$ as a vector and $M$ as a matrix. Finding the resonance
conditions and eigenmodes of the cavity now boils down to finding the
eigenvalues and eigenvectors of the roundtrip matrix.

We consider a plano-concave cavity with a large smooth planar mirror and a
smaller concave mirror. The reflection at the large planar mirror results in
eigenmodes with a flat wavefront at this mirror and forward and backward
propagating fields $|\psi^{+}\rangle$ and $|\psi^{-}\rangle$ that are mirror
imaged copies of each other. This symmetry allows one to write the round-trip
operator as a product $M=PB$, where $P$ describes the roundtrip propagation
and $B$ describes the reflection from the concave mirror.

Kleckner et al. [1] present an equation for the elements $B_{s,t}$ of the
reflection matrix, which couple a forward-propagating field
$|\psi_{t}^{+}\rangle$ into a backward-propagating field
$|\psi_{s}^{-}\rangle$. Unfortunately, their Eq. (5) is the complex conjugate
of the correct expression

$B_{s,t}=\langle\psi^{-}_{s}|B|\psi^{+}_{t}\rangle=\int_{0}^{r_{max}}\int_{0}^{2\pi}rdrd\theta\,\,\psi_{s}^{-*}\psi_{t}^{+}e^{-2ikz_{\rm
mirror}}\,,$ (S1)

where the integral, in polar coordinates $r$ and $\theta$, runs over the
finite extent of the mirror. This equation simply states that the reflected
field is a copy of the input field with a modified phase front and size,
assuming perfect mirror reflection. The function $z_{\rm mirror}(r,\theta)$
describes the shape of the concave mirror, where $z>0$ points towards the
planar mirror such that $z_{\rm mirror}(r,\theta)=r^{2}/(2R)$ for a
paraboloidal concave mirror with radius $R>0$.

It is instructive to discuss some properties of the roundtrip operator
$M=r_{1}r_{2}\exp(i2kL)B$, where $r_{1}$ and $r_{2}$ are the reflection
coefficients of the mirrors. The associated roundtrip matrix is symmetric,
i.e. $M_{ij}=M_{ji}$, even if the roundtrip phase delays of the two modes -
hidden in the phases of $\psi^{\pm}$ \- differs. The roundtrip operator is
unitary and the system is lossless (energy conserving) when mirror B is large
($r_{max}\rightarrow\infty$) and when both mirrors are perfect reflectors
($r_{1}=r_{2}=1$). Scattering is described by the off-diagonal elements of
$M$, which quantify the (energy-conserving) scattering between the different
modes per roundtrip.

The matrix elements of $M$ and $B$ depend on the choice of basis states. For
our geometry, the obvious choice is the set of Laguerre-Gaussian (LG) modes
with waist $w_{0}$ at the plane mirror and wavefront radius $R_{\rm
mode}=L+z_{0}^{2}/L=R$ at the concave mirror (Rayleigh range
$z_{0}=\frac{1}{2}kw_{0}^{2}$). The modes in this set are

$\psi_{p\ell}^{\pm}(r,\theta)=|\psi_{p\ell}(r)|\,\exp[\pm
i(\frac{kr^{2}}{2R}-\chi_{p\ell})+i\ell\theta]\,,$ (S2)

at $z=L$, where $p$ and $\ell$ are the radial and angular quantum number that
together label the modes. The Gouy phase $\chi_{p\ell}=(2p+|\ell|+1)\chi_{0}$,
with $\chi_{0}=\arcsin{\sqrt{L/R}}=\arctan(L/z_{0})$, is the single-pass phase
lag that LG-modes acquire with respect to a plane wave. The amplitude function
$|\psi_{p\ell}(r)|=C_{p\ell}x^{|\ell|/2}L_{p}^{\ell}(x)\exp(-x/2)$, with
$x=2r^{2}/w^{2}$ and normalization constant $C_{p\ell}$, describes the usual
LG-mode profiles. The angular momentum $\ell$ does not change in a coordinate
system with fixed orientation.
Next, we will calculate the roundtrip matrix for a rotational-symmetric
cavity, where $\ell$ is conserved, with $r_{1}=r_{2}\approx 1$. The on-
diagonal elements are

$\displaystyle M_{aa}$ $\displaystyle=$
$\displaystyle\exp(2ikL-2i\chi_{a})\int_{0}^{r_{max}}\pi
dr^{2}\,|\psi_{a}^{+}(r)|^{2}\,e^{-2ik\Delta z(r)}$ (S3)
$\displaystyle\approx$
$\displaystyle\exp(2ikL-2i\chi_{a}-2i\Delta\chi_{a})\times$
$\displaystyle\left(1-\int_{0}^{r_{max}}\pi dr^{2}\,2k^{2}[\Delta
z(r)]^{2}\,|\psi_{a}^{+}(r)|^{2}-\int_{r_{max}}^{\infty}\pi
dr^{2}\,|\psi_{a}^{+}(r)|^{2}\right)\,,$

where $\Delta z=z_{\rm mirror}-z_{\rm mode}$. The first equation combines Eq.
(S1) with $M=\exp(i2kL)B$ and with the Gouy phase lag $\chi_{a}$ of the
paraxial modes. The second equation is based on the Taylor expansion
$\exp(-2ik\Delta z)\approx 1-2ik\Delta z-2k^{2}[\Delta z]^{2}$. The imaginary-
valued term in this expansion results in a shift $\Delta\chi_{a}$ of the
resonance condition. The quadratic term in this expansion quantifies the
scattering loss induced by height variations. The final term in Eq. (S3)
quantifies the clipping loss at a mirror of radius $r_{max}$. Both losses are
amplitude losses. The middle term in Eq. (S3) thus yields the familiar
expression $(4\pi\sigma/\lambda)^{2}$ for the intensity loss, where
$\sigma^{2}=\langle\Delta z^{2}\rangle$ is the mean-square height difference
between mirror and phase front.

The coupling between modes $a$ and $b$ is described by the off-diagonal
element

$\displaystyle M_{ab}=M_{ba}=\int_{0}^{r_{max}}\pi
dr^{2}\,|\psi_{a}^{+}(r)\psi_{b}^{+}(r)|e^{-2ik\Delta z(r)}\approx\hskip
96.73936pt$ $\displaystyle\hskip 28.45274pt-\int_{r_{max}}^{\infty}\pi
dr^{2}\,|\psi_{a}^{+}(r)\psi_{b}^{+}(r)|-2ik\int_{0}^{r_{max}}\pi
dr^{2}\,\Delta z(r)\,|\psi_{a}^{+}(r)\psi_{b}^{+}(r)|\,,$ (S4)

where $\Delta z=z_{\rm mirror}-(z_{\rm a}+z_{\rm b})/2$ and where we again
assumed $k\Delta z\ll 1$. Modes with different angular momentum $\ell$ do not
coupling in a rotational-symmetric cavity.

We distinguish between two forms of coupling: dissipative and conservative.
Dissipate coupling is associated with energy loss and, in our system,
originates from clipping at the mirror. It is described by the first term in
Eq. (1) and produces a real-valued $M_{ab}$. Conservative coupling is energy
conserving and, in our system, originates from deformations of the wavefront
upon reflection. It is described by the second term in Eq. (1) and produces an
imaginary-valued $M_{ab}$.

Equations (S3) and (1) show that the coupling $M_{ab}$ is linear in $k\Delta
z$, whereas the scattering loss in $M_{aa}$ is quadratic in $k\Delta z$. The
relation between these two quantities is as follows. The relative intensity
loss by scattering from mode $a$ to mode $b$ is $|M_{ba}|^{2}$ per roundtrip,
making the scattering loss to all discrete modes equal to
$\sum_{b}|M_{ba}|^{2}$. This sum must be equal to the total scattering loss,
$1-M_{aa}^{2}$, from mode $a$. This relation is a consequence of the
conservative character of the operator $M$ in the large-mirror-limit
($r_{max}\rightarrow\infty$).

For our cavity with large concave mirror, dissipative coupling is expected to
be small and $M_{ab}$ will thus be (mainly) imaginary-valued and described by
the second term in Eq. (1). To support this statement we calculate the
expected clipping losses for three modes at a typical cavity length
$L=L_{3}=5.8\leavevmode\nobreak\ \mu$m, where the $p=0$ and $p=3$ modes are
frequency degenerate. For this calculation we also need the dimensions of the
concave mirror (radius of curvature $R=23.8\leavevmode\nobreak\ \mu$m, radius
$r_{max}=5.0\leavevmode\nobreak\ \mu$m) and the calculated size of the
fundamental mode $w_{z}=1.65\leavevmode\nobreak\ \mu$m, with
$I(r)\propto\exp(-2r^{2}/w_{z}^{2})$, at this mirror. The calculated
(amplitude) clipping loss of the fundamental ($p=\ell=0$) mode is order
$10^{-8}$ and hence negligible. The calculated (amplitude) clipping loss of
the ($p=3,\ell=0$) mode is $6\times 10^{-3}$, which limits its modal finesse
to $F\approx 250$, to be compared with an observed finesse of $F\approx 5000$
for the fundamental mode. But even a loss $\gamma_{b}=6\times 10^{-3}$ is
relatively small as compared to a typical detuning $\varphi_{a}-\varphi_{b}$
and can hence be neglected (see below).

We end this section with a disclaimer. The description presented above is
incomplete because the LG-modes are only the paraxial solutions of the full
propagation. Three non-paraxial corrections should be added to the propagator
$P$ when the opening angles are too large, i.e. when the modes are compact
($w_{0}\approx\lambda$) [2, 3, 4, 5]. These non-paraxial corrections lead to
small modifications of the resonance conditions and to an extra coupling
between the paraxial modes. The two scalar non-paraxial corrections, which are
the only ones we will consider, will modify the modal phase front from $z_{\rm
mode}=r^{2}/(2R)$ to $z_{\rm mode}\approx(1-\delta)r^{2}/(2R)+br^{4}/R^{3}$
with positive $\delta\ll 1$. For typical cavity lengths not too close to
concentric, $b<0$ and the non-paraxial phase front is “flatter" than
paraboloidal. The virtual absence of higher-order terms like $r^{6}$ and
$r^{8}$ (only even orders due to rotation symmetry) limits the coupling range.
More precisely, non-paraxial effects can only couple the fundamental $p=0$
mode to the $p=1$ and $p=2$ higher-order modes. For a paraboloidal reference
mirror, the predicted non-paraxial coupling between the $p=0$ and $p=2$ mode
is

$M_{ba}=M_{ab}=-\frac{3i}{2kR}\left(\frac{R-\frac{4}{3}L}{R-L}\right)\,.$ (S5)

For a cavity with $R=23.8\leavevmode\nobreak\ \mu$m we calculate
$M_{20}=-0.004\,i$ for $L=L_{2}=11.3\leavevmode\nobreak\ \mu$m. For our
mirrors, which are “steeper than paraboloidal", this extra non-paraxial
contribution increases the total coupling between the $\ell=0$ modes, but only
slightly.

##  2 Dynamic operator $K$ and coupled modes

Resonances in optical cavities behave like harmonic oscillators [6, 7]. They
have a Lorentzian spectrum under continuous excitation and the resonant field
decays exponential after transient excitation. These resonances are determined
by the condition that the intra-cavity field $|\psi_{i}\rangle$ reproduces
itself after a roundtrip, such that
$M|\psi_{i}\rangle=\exp[i(\varphi-\varphi_{i})-\gamma_{i}]|\psi_{i}\rangle$
where $\varphi=2kL$ and where $\varphi_{i}$ and $\gamma_{i}$ determine the
resonance condition and damping of the mode.

Below, we will demonstrate the equivalence between cavity resonances and
harmonic oscillators by treating the optical cavity as a multi-transverse-mode
Fabry-Perot interferometer. But we will start from a more general perspective.
The dynamics of any multi-mode system can be described by its dynamic operator
$K$, which is related to our roundtrip operation via the relation
$M\equiv\exp[i(\varphi-K)]$. This rewrite is useful because the relevant
dynamics is determined by $(1-M)\approx i(K-\varphi)$ instead of $M$ (see
below) and because the expansion highlights the relevant, close-to-resonance,
modes. The dynamic matrix can be written as $K={\Omega}-i{\Gamma}$, with
Hermitian matrices ${\Omega}$ and ${\Gamma}$, where ${\Omega}$ describes the
conservative dynamics while ${\Gamma}$ describes the dissipative dynamics [6].
In the absence of losses, the roundtrip operator $M=\exp(-iK)=\exp(-i\Omega)$
is unitary (energy-conserving ) and the dynamic matrix $K$ is Hermitian. Just
like in quantum mechanics, where the evolution operator $U=\exp(-iHt/\hbar)$
is unitary and the Hamiltonian $H$ is Hermitian.

As an aside, we note that this general approach has among others been used to
describe the dynamics of piano strings as a set of coupled harmonic
oscillators [8]. The description of that system also involves concepts like
the dynamics matrix $K$ and conservative versus dissipative coupling, which
Weinreich [8] calls reactive and resistive coupling, in analogy with
electronics.

Next, we treat the optical cavity as a multi-transverse-mode Fabry-Perot
interferometer. In a typical experiment, we illuminate the cavity with an
input field $|\psi_{in}\rangle$ and observe the output field
$|\psi_{out}\rangle$ as a function of cavity length or optical frequency. The
output field can be calculated by repeated operation of the roundtrip operator
and summation over the infinite series of reflections to find

$\left[1-M\right]|\psi_{out}\rangle\approx
i(K-\varphi)|\psi_{out}\rangle=t_{1}t_{2}|\psi_{in}\rangle\,.$ (S6)

For the single-mode case, $M$ and $K$ are scalars and the solution is trivial.
We then write $K=\varphi_{i}-i\gamma_{i}$, to find that the output field
indeed varies as a complex Lorentzian, i.e. as $\propto
1/[\gamma_{i}+i(\varphi_{i}-\varphi)]$, where $\varphi_{i}$ is the roundtrip
phase lag of mode $i$ and $(K-\varphi)$ is taken modulo $2\pi$. The modal loss
$\gamma_{i}=\gamma_{r}+\Delta\gamma_{i}$ combines the reflection loss
$\gamma_{r}\approx 1-r_{1}r_{2}$ with potential extra loss $\Delta\gamma_{i}$,
The resonance will repeat itself when $\varphi=2kL$ changes by $2\pi$. The
ratio between this repetition and the FWHM resonance width $2\gamma_{i}$
defines the modal finesse $F_{i}=\pi/\gamma_{i}$.

Next, we analyze the case for two close-to-resonance modes $a$ and $b$, which
together define the field
$|\psi\rangle=\alpha|\psi_{a}\rangle+\beta|\psi_{b}\rangle$, and write Eq.
(S6) as

$\begin{pmatrix}i(\varphi_{a}-\varphi)+\gamma_{a}&-M_{ab}\\\
-M_{ba}&i(\varphi_{b}-\varphi)+\gamma_{b}\end{pmatrix}\begin{pmatrix}\alpha\\\
\beta\end{pmatrix}_{out}=t_{1}t_{2}\begin{pmatrix}\alpha\\\
\beta\end{pmatrix}_{in}\,.$ (S7)

Equation (S7) describes the physics of two coupled harmonic oscillators.

We solve Eq. (S7) for the case of small mixing $M_{ab}$ and dominant
excitation of the less-lossy mode ($\alpha_{in}\gg\beta_{in}$ with
$\gamma_{a}<\gamma_{b}$). The off-diagonal element will then admix a small
fraction

$\left(\frac{\beta}{\alpha}\right)_{\rm
out}\approx\frac{M_{ba}}{\gamma_{ba}+i\Delta\varphi}\,,$ (S8)

of mode $b$ into mode $a$, where $\Delta\varphi=\varphi_{b}-\varphi_{a}$ is
the detuning between the resonances (= twice the difference of their single-
pass Gouy phases) and $\gamma_{ba}=\gamma_{b}-\gamma_{a}$ is the excess
roundtrip loss of mode $b$. The associated change $\Delta\lambda_{a}$ in the
eigenvalue $i\varphi_{a}+\gamma_{a}+\Delta\lambda_{a}$ of the mixed
(dominantly $a$) mode is

$\Delta\lambda_{a}\approx\frac{M^{2}_{ba}}{\gamma_{ba}+i\Delta\varphi}\,,$
(S9)

as $M_{ab}=M_{ba}$. At large detuning, $|\Delta\varphi|>\gamma_{ba}$, the
change $\Delta\lambda_{a}$ is mainly imaginary-valued. For our case of
conservative coupling, with imaginary-valued $M_{ab}=M_{ba}$, the mixing
fraction $(\beta/\alpha)_{\rm out}$ is mainly real-valued and the imaginary-
valued change in $\Delta\lambda_{a}$ results in a frequency shift of the two
resonance away from each other (avoided crossing). At smaller detuning,
$\Delta\lambda_{a}$ also acquires a real-valued component, which corresponds
to a change in loss.

##  3 AFM data yield mirror profile and coupling parameter

We have measured the profile of our concave mirror with atomic force
microscopy (AFM). Figure 1(a) shows that our mirror has a diameter $d\approx
10\leavevmode\nobreak\ \mu$m and a mirror indentation depth
$D=0.64(3)\leavevmode\nobreak\ \mu$m. It also shows that the central region
$r<2.5\leavevmode\nobreak\ \mu$m is approximately rotationally symmetric and
that deformations mainly occur for $r>2.5\leavevmode\nobreak\ \mu$m, in
particular in the lower region of the figure. These deformations are most
likely due to dust and might not even affect the optical data much.

(a)

(b)

Figure S1: (a) Height profile of our concave mirror as measured with AFM. (b)
Height difference $\Delta z=z-r^{2}/(2R)$ of mirror relative to reference
paraboloid. The red cross indicates the center of this reference.

Figure 1(b) shows the area in the white box of figure 1(a), in an expanded
view and after subtraction of a reference paraboloid $z_{\rm ref}=r^{2}/(2R)$
with $R=23.8\leavevmode\nobreak\ \mu$m from the height profile. This figure
shows that the central region of the mirror sticks out, which makes the
effective radius of curvature at the center larger than the reference
curvature. It also highlights the deformations that we mentioned earlier.

Figure S2: Rotation-averaged mirror height profile. The blue curve show the
AFM measurements. The orange curve shows a reference parabola $z=r^{2}/(2R)$
with $R=23.8\leavevmode\nobreak\ \mu$m. The green curve shows a spherical
reference surface with the same $R$. Figure S3: Height difference $\Delta z$
between the rotational-average AFM data and the ($R=23.8\leavevmode\nobreak\
\mu$m) reference parabola. The red, purple, and green curves show the mode
profile of the fundamental modes (dashed) and the nearly-degenerate higher-
order ($p=4,3,$ and 2) modes around the three degenerate cavity lengths.

Figure S2 show the rotational-average mirror height profile of the AFM data in
blue. It also shows the height profile of a reference parabola (orange) and a
reference circle (green), with the same radius of curvature
$R=23.8\leavevmode\nobreak\ \mu$m as deduced from the average transverse mode
splittings in the optical transmission spectra (Fig. 2a in main text). The AFM
data show that the mirror is much flatter than the reference in the center and
steeper than the reference towards the edges.

Figure S3 compares the rotational-average height profile of the mirror with
the ($R=23.8\leavevmode\nobreak\ \mu$m) reference parabola. It shows (i) the
height difference $\Delta z(r)=z_{\rm mirror}(r)-z_{\rm ref}(r)$ in blue and
(ii) the amplitude profiles $\psi(r)$ of the fundamental mode (dashed) and the
nearly-degenerate $p=4,3$ or 2 modes (solid) in red, purple, and green. Note
how the width $w_{z}$ of the fundamental mode at the concave mirror increases
from $p=4$ to 3 to 2 because this width $w_{z}(L)$ increases monotonically
with cavity length $L$. The mode profiles are calculated at the resonance
lengths $L_{4}=3.3\leavevmode\nobreak\ \mu$m, $L_{3}=5.8\leavevmode\nobreak\
\mu$m, and $L_{4}=11.3\leavevmode\nobreak\ \mu$m.

The data depicted in Fig. S3 allows one to calculate the coupling between the
depicted LG-modes, using the equation

$M_{p0}=(-1)\times\matrixelement{\psi_{p0}}{-2ik\Delta
z}{\psi_{00}}=\int_{0}^{\infty}dx\,\,2ik\Delta z(x)\,L_{p}(x)\exp(-x)\,,$
(S10)

where $x\equiv 2r^{2}/w_{z}^{2}$. The factor $(-1)\times$ accounts for the
different phase delay that the modes experience upon propagation from the flat
to the curved mirror. By defining the height difference $\Delta z=z_{\rm
mirror}-z_{\rm ref}$ with respect to a reference paraboloid, we implicitly
assume that the modal phase fronts have this paraboloidal shape; this
assumption proofs to be pretty good. Inspection of Fig. S3 shows that the
positive/negative mirror deformation $\Delta z$ generally overlaps with the
positive/negative part of the higher-order mode. This suggests a positive
coupling parameter $M_{p0}$ since the fundamental mode is positive over the
full range off the mirror. The calculated coupling parameters are
$M_{40}/i=0.015$, $M_{30}/i=0.049$, and $M_{20}/i=0.042$.

The analysis presented above was performed for LG-modes that are defined by
the condition that their radius of curvature is $R=23.8\leavevmode\nobreak\
\mu$m at the concave mirror. A disadvantage of this “fixed-radius" criterion
is that the curvature of the fundamental mode is typically not perfectly
matched to the curvature of the mirror. The fundamental matrix element
$M_{00}$ is therefore not as large as it can be if that radius of curvature
could be freely probed.

Alternatively, we can choose the radius of curvature of the fundamental mode
based on the criterion that it should provide the best match over the modal
region. This “best-match" criterion translates in a two-parameter fit for the
best paraboloid $z_{\rm ref}(x)=a+bx$ in the new coordinate $x\equiv
2r^{2}/w_{z}^{2}$. This fit minimizes the rms height difference $\sigma$
defined by

$\sigma^{2}=\matrixelement{\psi_{00}}{\Delta
z^{2}}{\psi_{00}}=\int_{0}^{\infty}dx\,\,[\Delta z(x)]^{2}\,\exp(-x)\,.$ (S11)

using linear regression. The fit parameter $b$ yields the effective radius
$R_{\rm eff}$ of the “best-match" fundamental mode. The calculated effective
radii at the three degenerate points are $R_{\rm eff}=27.0\leavevmode\nobreak\
\mu$m for $p=4$, $R_{\rm eff}=24.8\leavevmode\nobreak\ \mu$m for $p=3$ and
$R_{\rm eff}=22.7\leavevmode\nobreak\ \mu$m for $p=2$. This is consistent with
the earlier observation that the mirror height profile is “flatter at the
center and steeper towards the edges". It is also consistent with the value
$R=23.8\leavevmode\nobreak\ \mu$m deduced from fits of the optical
transmission spectra.

The “best-match" criterion yields a new beam size $w_{z}$ and a new set of LG-
modes. For this new set of LG modes, we again use Eq. (S10) to calculate the
coupling to find: $M_{40}/i=0.015$, $M_{30}/i=0.052$ and $M_{20}/i=0.053$.
These values are different from the ones presented earlier, but not
dramatically different. The rms fluctuations in the height profile are
$\sigma\approx 4$ nm, which is significantly more than the value
$\sigma\approx 2$ nm expected from the measured Finesse around 2000. A
possible explanation is that the cavity reshapes the fundamental mode in a
more optimal way then we expect from the current theory.

##  4 Repeated measurements on other microcavities

We have repeated the measurements of mode coupling with other microcavities.
We present the full analysis of one, to show that the behaviour is general.
This mirror is from a different batch than the mirror presented in the main
text.

The radius of curvature, $R$, and penetration depth, $L_{D}$, of the mirror
are again determined from the transverse-mode spectrum. The distance between
the fundamental and higher-order transverse mode is plotted in Fig. S4, which
resembles Fig. 2a in the main text. The fit yields
$R=17.4(1)\leavevmode\nobreak\ \mu$m and $L_{D}=0.04(3)\leavevmode\nobreak\
\mu$m.

The AFM measurement of this mirror (not shown) indicate that its indentation
depth $D\approx 0.9\leavevmode\nobreak\ \mu$m is larger than the value
$D\approx 0.6$ reported for the mirror in the main text. The AFM data also
shows that this second mirror is not as clean and has several defects as high
as $0.5\leavevmode\nobreak\ \mu$m next to the mirror. In practice, this means
that the smallest cavity length we can reach has a longitudinal mode number
$q=5$ instead of $q=3$ in the main text. The measured finesse of this
microcavity is also lower and the mirror coatings has a higher transmission.

Figure S5 shows the measured angular width of the fundamental mode, deduced
from single-Gaussian fits. The green line shows the theoretical behaviour of
the uncoupled fundamental modes; the fitted values are increased by 6%, to
correct an error in the magnification. The data follows the theoretical curve
throughout the full cavity length, but strongly deviates at the three points
of frequency degeneracy. These deviations have the same form as the ones
presented in the main text. This means that the coupling has the same nature,
and that the effect occurs for different micromirrors. This suggests that the
effect is general, rather then specific for each mirror.

Next, we use a two-mode fit for the CCD images to find the mode mixing
$c_{p0}$, shown in figure S6. The complex coupling constant is fitted, which
yield $M_{30}/i=0.021(2)$ and $M_{20}/i=0.028(1)$. These values are of the
same order as those in the main text, where we found $M_{30}/i=0.034(3)$ and
$M_{20}/i=0.029(1)$ for a different cavity.

Finally we have measured the finesse of the fundamental modes. Figure S7 shows
that the finesse decreases gradually with cavity length and shows prominent
features around the degenerate points labeled $N=4$ and $N=3$. The arrow in
Fig. S7 points at a significant drop in finesse around the coupling with the
$N=4$ ($p=2$) transverse modes. The observed asymmetry of this drop is
consistent with the positive sign of coupling deduced from Fig. S6; at
positive detuning, the size of the mode decreases at the flat mirror and
increases at the curved mirror, which results in more clipping losses and a
reduced finesse. The finesses around the $N=6$ and $N=8$ degenerate points
show no clear drop, indicating that mirror boundaries play no role here. The
asymmetry that we observed around the $N=3$ degenerate point is outside the
scope of this paper.

Figure S4: Transverse mode splitting versus mirror position for all transverse
modes we observed. Figure S5: The Gaussian opening angle $\theta_{0}$ of the
fundamental mode obtained from CCD images. The green line shows the theory for
the uncoupled system with a radius of curvature $R=17.4\leavevmode\nobreak\
\mu$m. All data points are shifted up by 6%, to correct for a misalignment
error in the imaging setup.  Figure S6: Mode mixing $c_{p0}=\beta/\alpha$ of
modes $p=4,3$, and $2$ into the fundamental mode. The $p=3$ (green) and $p=2$
(red) mixing is fitted to find their complex coupling parameters $M_{p0}$. The
$p=4$ mixing, around $L=2.4\leavevmode\nobreak\ \mu$m is not prominent enough
to allow such fitting. Figure S7: Finesse as a function of cavity length. The
vertical lines indicate the points of frequency degeneracy, with the even
modes in orange (from left to right N=8, 6, 4) and odd modes in green (N=7, 5,
3).

## References

  * [1] D. Kleckner, W. T. M. Irvine, S. S. R. Oemrawsingh, and D. Bouwmeester, “Diffraction-limited high-finesse optical cavities,” Physical Review A 81, 043814 (2010).
  * [2] A. L. Cullen and P. K. Yu, “Complex source-point theory of the electromagnetic open resonator,” Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 366, 155–171 (1979).
  * [3] C. Erickson, “High-order modes in a spherical Fabry-Perot resonator,” IEEE Trans. Micro. Theory and Techn. 23, 218–223 (1975).
  * [4] P. Yu and K. Luk, “High-order azimuthal modes in the open resonator,” Elec. Lett. 19, 539–541 (1983).
  * [5] K.-M. Luk, “Improvement in the resonant formula of a spherical Fabry–Perot resonator with unequal mirrors,” Journal of the Optical Society of America A 3, 3 (1986).
  * [6] R. Spreeuw, N. VanDruten, M. Beijersbergen, E. Eliel, and J. Woerdman, “Classical realization of a strongly driven two-level system,” Phys. Rev. Lett. 65, 2642–2645 (1990).
  * [7] H. Haus, _Waves and fields in optoelectronics_ (Prentice-Hall, 1984).
  * [8] G. Weinreich, “Coupled piano strings,” J. of the Acoustic Soc. of Am. 62, 1474–1484 (1977).