# Metastable Kitaev Magnets

Faranak Bahrami  Department of Physics, Boston College, Chestnut Hill, MA
02467, USA  Mykola Abramchuk  Department of Physics, Boston College, Chestnut
Hill, MA 02467, USA  Oleg I. Lebedev  Laboratoire CRISMAT, ENSICAEN-CNRS
UMR6508, 14050 Caen, France  Fazel Tafti  Department of Physics, Boston
College, Chestnut Hill, MA 02467, USA<EMAIL_ADDRESS>
###### Abstract

Nearly two decades ago, Alexei Kitaev proposed a model for spin-$1/2$
particles with bond-directional interactions on a two-dimensional honeycomb
lattice which had the potential to host a quantum spin-liquid ground state.
This work initiated numerous investigations to design and synthesize materials
that would physically realize the Kitaev Hamiltonian. The first-generation of
such materials, such as Na2IrO3, $\alpha$-Li2IrO3, and $\alpha$-RuCl3,
revealed the presence of non-Kitaev interactions such as the Heisenberg and
off-diagonal exchange. Both physical pressure and chemical doping were used to
tune the relative strength of the Kitaev and competing interactions; however,
little progress was made towards achieving a purely Kitaev system. Here, we
review the recent breakthrough in modifying Kitaev magnets via topochemical
methods that has led to the second-generation of Kitaev materials. We show how
structural modifications due to the topotactic exchange reactions can alter
the magnetic interactions in favor of a quantum spin-liquid phase.

##  1 Introduction

Recently, the 4d/5d honeycomb layered materials have been vigorously studied
due to their potential in realizing a quantum spin-liquid (QSL) ground state
1, 2, 3, 4, 5, 6, 7, 8. First introduced by Alexei Kitaev in 2006, the Kitaev
model is an exactly solvable theoretical model with bond-dependent Ising
interactions among spin-$1/2$ degrees of freedom on a two-dimensional (2D)
honeycomb lattice, which is described by the Kitaev Hamiltonian:
$\mathcal{H}=-\sum K_{\gamma}\mathbf{S_{i}}^{\gamma}\mathbf{S_{j}}^{\gamma}$
9. The ground state of this system is magnetically frustrated and is predicted
to be a QSL 9. The applications of a Kitaev QSL in quantum information and the
possibility of realizing Majorana fermions have inspired numerous
investigations on quasi-2D honeycomb materials 3, 10, 11, 1, 12.

The first-generation of these compounds, namely Na2IrO3, $\alpha$-Li2IrO3,
Li2RhO3, and $\alpha$-RuCl3, were synthesized using conventional solid state
methods at high temperatures ($T>700$ ∘C). In these materials, heavy
transition metal ions (Ru3+, Rh4+, and Ir4+) are octahedrally coordinated with
oxygen or chlorine atoms (Fig. 1a), and the edge-sharing octahedra create
honeycomb layers (Fig. 1b). The combination of octahedral crystal electric
field (CEF) and strong spin-orbit coupling (SOC) splits the five-fold
degenerate $d$-levels and leaves one electron in the isospin-$1/2$
($J_{\text{eff}}$=$1/2$) state necessary for the Kitaev model (Fig. 1c) 1, 4,
7, 13, 9.

Figure 1:  (a) The bond angle ($\phi$) between edge-shared octahedral units
plays a significant role in tuning the magnetic interactions. (b) Edge-sharing
octahedral units create a honeycomb structure in Kitaev magnets such as
$\alpha$-Li2IrO3 and Na2IrO3. (c) Interplay between CEF and SOC creates the
isospin-1/2 state in the Kitaev magnets.

Finding new Kitaev magnets, beyond the first-generation compounds, has become
a frontier challenge in solid state chemistry. Prior attempts to replace Na
with K in Na2IrO3 or replacing Cl with Br in $\alpha$-RuCl3 have led to other
stable phases with different structures instead of the honeycomb lattice 14,
15. The amount of physical pressure required to substantially tune the
interactions is too high 16 and chemical doping leads to a change of spin
state 17. Therefore, recent success in synthesizing a second-generation of
Kitaev magnets where magnetic interactions can be tuned by topochemical
methods has revitalized the field. In this review, we will first explain the
different types of exchange reactions (partial and complete), then discuss the
interplay between topochemical reactions and magnetism, and finally present
heat capacity and magnetization data to compare the properties of the first
and second-generation Kiteav magnets.

##  2 Topotactic Exchange Reactions

Figure 2:  Synthesis of the second-generation Kitaev magnets from the first-
generation materials through (a) partial and (b) complete exchange reactions.
Both generations have honeycomb layers. The topochemical change of inter-layer
coordination from octahedral to linear modifies the intra-layer Ir-O-Ir bond
angles due to the change of oxygen positions.

The second-generation Kitaev magnets are metastable compounds, i.e. they have
a higher enthalpy of formation and a lower decomposition threshold compared to
stable counterparts 18. Thus, it is impossible to synthesize them with
conventional solid state methods at high temperatures. Instead, they are
stabilized through topochemical reactions from the first-generation compounds
under mild conditions. As shown schematically in Figs. 2 the global symmetries
of the unit cell (space group and honeycomb structure) do not change during a
topochemical reaction. However, the local parameters such as bond lengths and
bond angles are modified efficiently.

Topotactic exchange reactions can be either partial (Fig. 2a) or complete
(Fig. 2b). The most general formulation of a partial exchange reaction is

$\mathrm{2A_{2}MO_{3}+3BX}~{}\rightarrow~{}\mathrm{B_{3}AM_{2}O_{6}+3AX}$ (1)

where the inter-layer A-atoms (typically Li or Na) in a stable honeycomb
structure A2MO3 are exchanged with the B-atoms (typically Cu, Ag, and H) from
a halide, nitrate, or sulfate compound BX. For example, Fig. 2a corresponds to
A=Li, B=Ag, M=Ir, and X=NO3 for the synthesis of Ag3LiIr2O6 from
$\alpha$-Li2IrO3. Replacing the inter-layer Li atoms by H, Cu, or Ag, in
$\alpha$-Li2IrO3 has recently produced H3LiIr2O6, Cu3LiIr2O6, and Ag3LiIr2O6,
respectively 19, 20, 21.

In a complete topotactic exchange reaction, all A-atoms within and between the
layers are replaced by the B-atoms.

$\mathrm{A_{2}MO_{3}+2BX}~{}\rightarrow~{}\mathrm{B_{2}MO_{3}+2AX}$ (2)

For example, Fig. 2b corresponds to A=Na, B=Cu, M=Ir, and X=Cl for the
synthesis of Cu2IrO3 from Na2IrO3. A complete exchange reaction is much less
likely to happen and so far Cu2IrO3 is the only known system in this category
22. It is noteworthy that the copper atoms in Cu2IrO3 are not entirely in a
Cu+ state. Both x-ray absorption and electron energy loss spectroscopy (XAS
and EELS) confirmed a mixed valence of Cu+/Cu2+ with a 1/1 ratio within the
honeycomb layers 23. A mixed valence of copper, induces a mixed valence of
iridium (Ir+3/Ir+4) and leads to magnetic disorder and spin-glass behavior 23,
24.

##  3 Synthesis Details

Figure 3:  (a) After each heat cycle, the powder x-ray pattern of
$\alpha$-Li2IrO3 shows more pronounced peaks, especially between 20 and 30
degrees where the honeycomb Bragg peaks appear. The number of times each
sample has been reheated is shown on the right above its respective pattern.
(b) The x-ray patterns of two second-generation Kitaev systems, H3LiIr2O6
(green) and Ag3LiIr2O6 (gray data, reproduced from 25). The inset shows the
asymmetric broadening of the honeycomb Bragg peaks in Ag3LiIr2O6 due to
stacking faults. In H3LiIr2O6, the honeycomb peaks are hardly discernible due
to high structural disorder.

The first-generation Kitaev magnets are prepared via conventional solid state
reaction at high temperatures (T$\geq 700$ K) either in air or under the flow
of oxygen/argon gas 13, 7, 26. To improve the sample quality and remove
stacking faults, it is necessary to perform successive stages of grinding and
heating. For example, the x-ray patterns in Fig. 3a show that the quality of
$\alpha$-Li2IrO3 samples improve by repeating the heat cycles. Specifically,
the superstructure peaks between 20 and 30 degrees (inset of Fig. 3a) that
represent the honeycomb ordering become more pronounced in each iteration.
Typically, improving the quality of the first-generation compound will improve
the quality of the second-generation material after the exchange reaction 25.

The topotactic cation exchange reaction must be conducted at low temperatures
($T\leq 400$ K) 27, 22, 25, since higher temperatures will decompose the
metastable product. The second-generation Kitaev magnets are prepared by
modifying the inter-layer atoms and the associated chemical bonds, and
therefore they have more stacking faults than their parent compounds 25, 28.
This can be seen in the inset of Fig. 3b that shows an asymmetric broadening
of the honeycomb Bragg peaks in Ag3LiIr2O6. Unlike the solid state reactions,
topotactic exchange cannot be repeated to improve the sample quality. Thus,
getting rid of the stacking faults in these materials remains an open
challenge.

Details of the synthesis procedures for Cu2IrO3 and Ag3LiIr2O6 have previously
been published by Abramchuk and Bahrami _et al._ 22, 25, 19. Here, we present
more details about the synthesis of H3LiIr2O6 based on the earlier work of
Bette _et al._ 29. Polycrystalline samples of H3LiIr2O6 are synthesized using
a modified version of Eq. 1.

$\mathrm{4Li_{2}IrO_{3}+3H_{2}SO_{4}}~{}\rightarrow~{}\mathrm{2H_{3}LiIr_{2}O_{6}+3Li_{2}SO_{4}}$
(3)

After synthesizing a high-quality of $\alpha$-Li2IrO3 (Fig. 3a), approximately
300 mg of the material was added to a 10 ml Teflon-lined steel autoclave
filled with H2SO4 acid (1 M solution) and heated to 120 ∘C for several days.
After completing the reaction, the product was washed with water and the
quality was verified using x-ray diffraction (Fig. 3b).

##  4 Stacking Faults

Figure 4:  HAADF-TEM images from (a) $\alpha$-Li2IrO3 and (b) Ag3LiIr2O6. The
images show an abundance of stacking faults in Ag3LiIr2O6 unlike
$\alpha$-Li2IrO3, due to the weaker inter-layer bonding in the former. The
electron diffraction patterns are presented as insets and reveal fewer
streaking in $\alpha$-Li2IrO3 due to fewer stacking faults compared to
Ag3LiIr2O6.

A comparison between the insets of Figs. 3a and 3b suggests fewer stacking
faults in $\alpha$-Li2IrO3 (sharp and well-separated Bragg peaks from the
honeycomb layers) and considerable stacking faults in Ag3LiIr2O6 (broadened
peaks). The asymmetric broadening of honeycomb peaks is known as the Warren
line shape, which is a signature of stacking disorder 30. The higher amount of
stacking faults in the second-generation Kitaev magnets is due to the inter-
layer chemistry. As seen in Fig. 2, each inter-layer Li atom in
$\alpha$-Li2IrO3 is octahedrally coordinated with three oxygen atoms from the
top and three from the bottom honeycomb layers. In contrast, each Ag atom in
Ag3LiIr2O6 is connected to only one O from the top and one from the bottom
layer in a dumbbell (linear) coordination. The weak dumbbell bonds are
responsible for the larger inter-layer separation in Ag3LiIr2O6 and more
stacking faults compared to $\alpha$-Li2IrO3 31.

Direct lattice imaging with transmission electron microscope (TEM) is a
powerful tool to study the stacking faults. Figures 4a and 4b (reproduced from
Ref. 25) are high angle annular dark field TEM (HAADF-TEM) images of
$\alpha$-Li2IrO3 and Ag3LiIr2O6 samples, respectively. Whereas the stacking
sequence in $\alpha$-Li2IrO3 can be flawless for up to 50 unit cells,
Ag3LiIr2O6 shows a maximum of 5 unit cells stacked without faults (in the form
of twisting between the layers). In H3LiIr2O6, the small size of H atoms and
their high mobility make the chemical bonds even weaker than in Ag3LiIr2O6. As
such, H3LiIr2O6 has the highest degree of stacking faults among the second-
generation Kitaev magnets 29, 27, 28. This is why the honeycomb peaks of
H3LiIr2O6 are not resolved by x-rays (Fig. 3b).

##  5 Tuning magnetic interactions with topochemical methods

Figure 5:  Exchange paths for (a) $K$, (b) $J$, and (c) $\Gamma$ terms in Eq.
4. The $d$ and $p$ orbitals are painted in blue and red, respectively. The
numbers show the hopping sequence in the perturbation.

As seen in Fig. 2, the monoclinic unit cell and the honeycomb ordering in the
2D layers remain unchanged before and after exchange reactions. However, the
change of inter-layer coordination from octahedral to dumbbell modifies the
M-O-M bond angles within the honeycomb layers (Fig. 1a and 2). Superexchange
magnetic interactions are sensitive to a change of bond angles and thus
topochemical reactions can be used to tune the magnetic interactions. There
are at least three terms in the magnetic Hamiltonian of the Kitaev materials.

$\mathcal{H}=\sum\limits_{\langle
i,j\rangle\in{\alpha\beta(\gamma)}}\left[-K_{\gamma}S_{i}^{\gamma}S_{j}^{\gamma}+J\textbf{S}_{i}\cdot\textbf{S}_{j}+\Gamma\left(S_{i}^{\alpha}S_{j}^{\beta}+S_{i}^{\beta}S_{j}^{\alpha}\right)\right]$
(4)

The Kitaev term ($K$) favors QSL, the Heisenberg term ($J$) favors AFM
ordering, and the off-diagonal exchange term ($\Gamma$) controls details of
the ordered structure. All three term can be modified via topochemical
reactions as follows.

Figure 5 shows the individual exchange paths for each term in Eq. 4. The
Kitaev term is an indirect exchange interaction with hopping matrix elements
$t_{dpd}$ between the $d_{xz}$, $p_{z}$, and $d_{yz}$ orbitals (Fig. 5a) 32,
33. In addition to the indirect exchange ($K$), Fig. 5b shows a direct
exchange path for the Heisenberg interaction ($J$) with hopping matrix element
$t_{dd}$ between $d_{xy}$ orbitals, leading to $J\sim t_{dd}^{2}/U$ in Eq. 4
34. Finally, a combination of direct and indirect paths in Fig. 5c leads to
the symmetric off-diagonal exchange, $\Gamma\sim t_{dpd}t_{dd}J_{H}/U^{2}$,
where $J_{H}$ is the Hund’s coupling between the $e_{\mathrm{g}}$ and
$t_{\mathrm{2g}}$ orbitals 35, 36. The hopping matrix elements ($t_{dd}$ and
$t_{dpd}$) are tuned by the M-O-M bond angle and the M-M distance which can be
tuned by the exchange reactions. For example, (i) the change of oxygen
positions within the honeycomb layers due to the change of inter-layer
coordinations in Fig. 2 modifies the M-O-M bond angle ($\phi$ in Fig. 1a and
5a) and thereby tunes $t_{dpd}$; (ii) according to theoreical calculations 1,
the Heisenberg interaction is canceled between the opposite paths if the bond
angle $\phi$ is close to $90^{\circ}$ (Fig. 1a and 5a); (iii) the
hybridization between the Ag $d$-orbitals between the layers and O
$p$-orbitals within the layers tunes the ratio of $t_{dpd}/t_{dd}$.

##  6 Magnetic characterization of metastable Kitaev materials

Figure 6:  (a) Heat capacity ($C/T$) plotted as a function of temperature
below 30 K for the first-generation Kitaev magnet $\alpha$-Li2IrO3 and its
second-generation derivatives Ag3LiIr2O6 and H3LiIr2O6. The data for
$\alpha$-Li2IrO3 and Ag3LiIr2O6 are reproduced from Refs. 2, 25. (b) A similar
comparison is made between Na2IrO3 (first-generation) and Cu2IrO3 (second-
generation). The data are reproduced from Ref. 22.

To demonstrate the effect of topochemical modifications on the magnetic
interactions (Eq. 4 and Fig. 5), we compare the heat capacity and magnetic
susceptibility of the first and second-generation Kitaev magnets. The peak in
the heat capacity of $\alpha$-Li2IrO3 in Fig. 6a confirms long-range magnetic
ordering at $T_{N}=15$ K. The order has been characterized as incommensurate
spiral by recent neutron scattering and muon spin relaxation ($\mu$SR)
experiments 5, 8. As seen in Fig. 6a, this peak is shifted to lower
temperatures in Ag3LiIr2O6 and seemingly disappeared in H3LiIr2O6. The
suppression of $T_{N}$ in second-generation compounds Ag3LiIr2O6 and H3LiIr2O6
is a positive sign of approaching the quantum spin-liquid (QSL) phase, where
long-range order is replaced by long-range quantum entanglement 3, 10. A
recent $\mu$SR experiment 25 has shown a similar incommensurate spiral order
in Ag3LiIr2O6; however, the long-range order develops at 8 K in Ag3LiIr2O6,
well below $T_{N}=15$ K in $\alpha$-Li2IrO3. Thus, the topochemical
modification of bond angles seem to strengthen $K$ and weaken $J$ in Eq. 4. A
recent nuclear magnetic resonance (NMR) experiment has shown absence of long-
range order in H3LiIr2O6, which is another promising result toward the
discovery of a QSL phase 27.

A similar trend is observed in Fig. 6b for the first-generation material
Na2IrO3 that shows a peak at $T_{N}=15$ K and its second-generation
counterpart Cu2IrO3 that does not show a peak but seems to have a broad
anomaly below 5 K. Neutron scattering experiments have confirmed a zigzag
antiferromagnetic (AFM) order in Na2IrO3 37. Recent $\mu$SR and NMR
experiments have revealed a coexistence of static and dynamic magnetism below
5 K in Cu2IrO3 but without a long-range order, suggesting proximity to the QSL
phase 23, 38.

Figure 7:  (a) Magnetic susceptibility ($\chi$) plotted as a function of
temperature below 30 K for the first-generation Kitaev magnet $\alpha$-Li2IrO3
and its second-generation derivatives Ag3LiIr2O6 and H3LiIr2O6. The data for
$\alpha$-Li2IrO3 and Ag3LiIr2O6 are reproduced from Refs. 19, 25. (b) A
similar comparison is made between Na2IrO3 (first-generation) and Cu2IrO3
(second-generation). The data for Na2IrO3 and Cu2IrO3 are reproduced from
Refs. 2, 22.

The suppression of magnetic ordering due to topochemcial changes in metastable
Kitaev magnets is also observed in the magnetic susceptibility data. Figure 7a
shows the magnetic susceptibility of $\alpha$-Li2IrO3 (black curve) with a
clear anomaly at $T_{N}=15$ K indicating the incommensurate spiral
antiferromagnetic (AFM) order. The green curve representing Ag3LiIr2O6 shows
two downturns at $T_{F}=14$ K and $T_{N}=8$ K, corresponding to the onsets of
short-range and long-range magnetic orders, respectively 25. The orange curve
representing H3LiIr2O6 does not show any evidence of magnetic ordering. Figure
7b shows a similar trend, where the first-generation material Na2IrO3 orders
at $T_{N}=15$ K and the second-generation material Cu2IrO3 shows a small peak
at 2 K, evidence of short-range spin freezing instead of long-range order.

##  7 Challenges and Opportunities

The above results are exciting; however, they need to be interpreted with
caution. Topotactic exchange reactions increase disorder that has adverse
effects on magnetism. A recent TEM study has shown that the silver atoms in
Ag3LiIr2O6 can enter the honeycomb layers and form small inclusions (up to 50
atoms) that disrupt the magnetic ordering 25. Such a structural disorder can
spuriously hide the long-range order and be misinterpreted as evidence of a
QSL phase. As noted earlier, H3LiIr2O6 is even more disordered compared to
Ag3LiIr2O6 due to the high mobility of the H atoms, which causes bond
randomness and site vacancies within the honeycomb layers 29. Recent
theoretical works show that the absence of magnetic ordering in H3LiIr2O6 may
be due to bond randomness and a large amount of vacancies 39, 40. Thus, the
most important challenge in this field is to optimize the synthesis conditions
for a minimum amount of disorder and to find methods of annealing away the
stacking faults and vacancies.

Metastable Kitaev magnets have opened a new window of opportunity to realizing
the quantum spin-liquid ground state. The Majorana excitations of such
materials will form the building blocks of a solid state quantum computer 41.
Braiding algorithms and logical gates have been theoretically developed for
such computers 42. It remains an open challenge for the solid state chemistry
community to synthesize the appropriate materials for such models. Another
intriguing opportunity is to find unconventional superconductivity in the
Kitaev magnets 43, an exciting theoretical prediction that awaits experimental
discovery.

##  8 ACKNOWLEDGMENTS

This material is based upon work supported by the Air Force Office of
Scientific Research under award number FA2386-21-1-4059.

## References

  * Jackeli and Khaliullin 2009 Jackeli, G.; Khaliullin, G. Mott Insulators in the Strong Spin-Orbit Coupling Limit: From Heisenberg to a Quantum Compass and Kitaev Models. _Phys. Rev. Lett._ 2009, _102_ , 017205
  * Singh et al. 2012 Singh, Y.; Manni, S.; Reuther, J.; Berlijn, T.; Thomale, R.; Ku, W.; Trebst, S.; Gegenwart, P. Relevance of the Heisenberg-Kitaev Model for the Honeycomb Lattice Iridates ${A}_{2}{\mathrm{IrO}}_{3}$. _Phys. Rev. Lett._ 2012, _108_ , 127203
  * Takagi et al. 2019 Takagi, H.; Takayama, T.; Jackeli, G.; Khaliullin, G.; Nagler, S. E. Concept and realization of Kitaev quantum spin liquids. _Nature Reviews Physics_ 2019, _1_ , 264–280
  * Plumb et al. 2014 Plumb, K. W.; Clancy, J. P.; Sandilands, L. J.; Shankar, V. V.; Hu, Y. F.; Burch, K. S.; Kee, H.-Y.; Kim, Y.-J. $\alpha-{\mathrm{RuCl}}_{3}$: A spin-orbit assisted Mott insulator on a honeycomb lattice. _Phys. Rev. B_ 2014, _90_ , 041112
  * Williams et al. 2016 Williams, S. C.; Johnson, R. D.; Freund, F.; Choi, S.; Jesche, A.; Kimchi, I.; Manni, S.; Bombardi, A.; Manuel, P.; Gegenwart, P.; Coldea, R. Incommensurate counterrotating magnetic order stabilized by Kitaev interactions in the layered honeycomb $\alpha\text{$-$}{\mathrm{Li}}_{2}{\mathrm{IrO}}_{3}$. _Phys. Rev. B_ 2016, _93_ , 195158
  * Wang et al. 2020 Wang, Y.; Osterhoudt, G. B.; Tian, Y.; Lampen-Kelley, P.; Banerjee, A.; Goldstein, T.; Yan, J.; Knolle, J.; Ji, H.; Cava, R. J.; Nasu, J.; Motome, Y.; Nagler, S. E.; Mandrus, D.; Burch, K. S. The range of non-Kitaev terms and fractional particles in $\alpha$-RuCl3. _npj Quantum Materials_ 2020, _5_ , 1–8
  * Todorova and Jansen 2011 Todorova, V.; Jansen, M. Synthesis, Structural Characterization and Physical Properties of a New Member of Ternary Lithium Layered Compounds – Li2RhO3. _Zeitschrift für anorganische und allgemeine Chemie_ 2011, _637_ , 37–40
  * Choi et al. 2019 Choi, S.; Manni, S.; Singleton, J.; Topping, C. V.; Lancaster, T.; Blundell, S. J.; Adroja, D. T.; Zapf, V.; Gegenwart, P.; Coldea, R. Spin dynamics and field-induced magnetic phase transition in the honeycomb Kitaev magnet $\alpha\text{$-$}{\mathrm{Li}}_{2}{\mathrm{IrO}}_{3}$. _Phys. Rev. B_ 2019, _99_ , 054426
  * Kitaev 2006 Kitaev, A. Anyons in an exactly solved model and beyond. _Annals of Physics_ 2006, _321_ , 2–111
  * Knolle and Moessner 2019 Knolle, J.; Moessner, R. A Field Guide to Spin Liquids. _Annual Review of Condensed Matter Physics_ 2019, _10_ , 451–472
  * Chaloupka et al. 2010 Chaloupka, J.; Jackeli, G.; Khaliullin, G. Kitaev-Heisenberg Model on a Honeycomb Lattice: Possible Exotic Phases in Iridium Oxides A2IrO3. _Physical Review Letters_ 2010, _105_ , 027204
  * Kim et al. 2008 Kim, B. J.; Jin, H.; Moon, S. J.; Kim, J.-Y.; Park, B.-G.; Leem, C. S.; Yu, J.; Noh, T. W.; Kim, C.; Oh, S.-J.; Park, J.-H.; Durairaj, V.; Cao, G.; Rotenberg, E. Novel ${J}_{\mathrm{eff}}=1/2$ Mott State Induced by Relativistic Spin-Orbit Coupling in Sr2IrO4. _Physical Review Letters_ 2008, _101_ , 076402
  * Mehlawat et al. 2017 Mehlawat, K.; Thamizhavel, A.; Singh, Y. Heat capacity evidence for proximity to the Kitaev quantum spin liquid in A2IrO3 (A = Na, Li). _Physical Review B_ 2017, _95_ , 144406
  * Weber et al. 2017 Weber, D.; Schoop, L. M.; Wurmbrand, D.; Nuss, J.; Seibel, E. M.; Tafti, F. F.; Ji, H.; Cava, R. J.; Dinnebier, R. E.; Lotsch, B. V. Trivalent Iridium Oxides: Layered Triangular Lattice Iridate K0.75Na0.25IrO2 and Oxyhydroxide IrOOH. _Chemistry of Materials_ 2017, _29_ , 8338–8345
  * Merlino et al. 2004 Merlino, S.; Labella, L.; Marchetti, F.; Toscani, S. Order-disorder transformation in RuBr3 and MoBr3: A two-dimensinal Ising model. _Chemistry of Materials_ 2004, _16_ , 3895–3903
  * Bhattacharjee et al. 2012 Bhattacharjee, S.; Lee, S.-S.; Kim, Y. B. Spin–orbital locking, emergent pseudo-spin and magnetic order in honeycomb lattice iridates. _New Journal of Physics_ 2012, _14_ , 073015
  * Cao and Schlottmann 2018 Cao, G.; Schlottmann, P. The challenge of spin–orbit-tuned ground states in iridates: a key issues review. _Reports on Progress in Physics_ 2018, _81_ , 042502
  * Aykol et al. 2018 Aykol, M.; Dwaraknath, S. S.; Sun, W.; Persson, K. A. Thermodynamic limit for synthesis of metastable inorganic materials. _Science Advances_ 2018, _4_ , eaaq0148
  * Bahrami et al. 2019 Bahrami, F.; Lafargue-Dit-Hauret, W.; Lebedev, O. I.; Movshovich, R.; Yang, H.-Y.; Broido, D.; Rocquefelte, X.; Tafti, F. Thermodynamic Evidence of Proximity to a Kitaev Spin Liquid in ${\mathrm{Ag}}_{3}{\mathrm{LiIr}}_{2}{\mathrm{O}}_{6}$. _Phys. Rev. Lett._ 2019, _123_ , 237203
  * Geirhos et al. 2020 Geirhos, K.; Lunkenheimer, P.; Blankenhorn, M.; Claus, R.; Matsumoto, Y.; Kitagawa, K.; Takayama, T.; Takagi, H.; Kézsmárki, I.; Loidl, A. Quantum paraelectricity in the Kitaev quantum spin liquid candidates ${\mathrm{H}}_{3}\mathrm{LiI}{\mathrm{r}}_{2}{\mathrm{O}}_{6}$ and ${\mathrm{D}}_{3}\mathrm{LiI}{\mathrm{r}}_{2}{\mathrm{O}}_{6}$. _Phys. Rev. B_ 2020, _101_ , 184410
  * Roudebush et al. 2013 Roudebush, J. H.; Andersen, N. H.; Ramlau, R.; Garlea, V. O.; Toft-Petersen, R.; Norby, P.; Schneider, R.; Hay, J. N.; Cava, R. J. Structure and Magnetic Properties of Cu3Ni2SbO6 and Cu3Co2SbO6 Delafossites with Honeycomb Lattices. _Inorganic Chemistry_ 2013, _52_ , 6083–6095
  * Abramchuk et al. 2017 Abramchuk, M.; Ozsoy-Keskinbora, C.; Krizan, J. W.; Metz, K. R.; Bell, D. C.; Tafti, F. Cu2IrO3: A New Magnetically Frustrated Honeycomb Iridate. _Journal of the American Chemical Society_ 2017, _139_ , 15371–15376
  * Kenney et al. 2019 Kenney, E. M. et al. Coexistence of static and dynamic magnetism in the Kitaev spin liquid material Cu2IrO3. _Physical Review B_ 2019, _100_ , 094418
  * Choi et al. 2019 Choi, Y.; Lee, C.; Lee, S.; Yoon, S.; Lee, W.-J.; Park, J.; Ali, A.; Singh, Y.; Orain, J.-C.; Kim, G.; Rhyee, J.-S.; Chen, W.-T.; Chou, F.; Choi, K.-Y. Exotic Low-Energy Excitations Emergent in the Random Kitaev Magnet Cu2IrO3. _Physical Review Letters_ 2019, _122_ , 167202
  * Bahrami et al. 2021 Bahrami, F.; Kenney, E. M.; Wang, C.; Berlie, A.; Lebedev, O. I.; Graf, M. J.; Tafti, F. Effect of structural disorder on the Kitaev magnet ${\mathrm{Ag}}_{3}{\mathrm{LiIr}}_{2}{\mathrm{O}}_{6}$. _Phys. Rev. B_ 2021, _103_ , 094427
  * O’Malley et al. 2008 O’Malley, M. J.; Verweij, H.; Woodward, P. M. Structure and properties of ordered Li2IrO3 and Li2PtO3. _Journal of Solid State Chemistry_ 2008, _181_ , 1803–1809
  * Kitagawa et al. 2018 Kitagawa, K.; Takayama, T.; Matsumoto, Y.; Kato, A.; Takano, R.; Kishimoto, Y.; Bette, S.; Dinnebier, R.; Jackeli, G.; Takagi, H. A spin–orbital-entangled quantum liquid on a honeycomb lattice. _Nature_ 2018, _554_ , 341–345
  * 28 Tsirlin, A. A.; Gegenwart, P. Kitaev Magnetism through the Prism of Lithium Iridate. _physica status solidi (b)_ _n/a_ , 2100146
  * Bette et al. 2017 Bette, S.; Takayama, T.; Kitagawa, K.; Takano, R.; Takagi, H.; Dinnebier, R. E. Solution of the heavily stacking faulted crystal structure of the honeycomb iridate H3LiIr2O6. _Dalton Transactions_ 2017, _46_ , 15216–15227
  * Balzar 1993 Balzar, D. X-Ray Diffraction Line Broadening: Modeling and Applications to High-Tc Superconductors. _Journal of Research of the National Institute of Standards and Technology_ 1993, _98_ , 321–353
  * Abramchuk et al. 2018 Abramchuk, M.; Lebedev, O. I.; Hellman, O.; Bahrami, F.; Mordvinova, N. E.; Krizan, J. W.; Metz, K. R.; Broido, D.; Tafti, F. Crystal Chemistry and Phonon Heat Capacity in Quaternary Honeycomb Delafossites: Cu[Li1/3Sn2/3]O2 and Cu[Na1/3Sn2/3]O2. _Inorganic Chemistry_ 2018, _57_ , 12709–12717
  * Rau et al. 2016 Rau, J. G.; Lee, E. K.-H.; Kee, H.-Y. Spin-Orbit Physics Giving Rise to Novel Phases in Correlated Systems: Iridates and Related Materials. _Annual Review of Condensed Matter Physics_ 2016, _7_ , 195–221
  * Kim and Kee 2016 Kim, H.-S.; Kee, H.-Y. Crystal structure and magnetism in $\alpha$-RuCl3: An ab initio study. _Physical Review B_ 2016, _93_ , 155143
  * Winter et al. 2016 Winter, S. M.; Li, Y.; Jeschke, H. O.; Valentí, R. Challenges in design of Kitaev materials: Magnetic interactions from competing energy scales. _Physical Review B_ 2016, _93_ , 214431
  * Rau et al. 2014 Rau, J. G.; Lee, E. K.-H.; Kee, H.-Y. Generic Spin Model for the Honeycomb Iridates beyond the Kitaev Limit. _Phys. Rev. Lett._ 2014, _112_ , 077204
  * Rusnacko et al. 2019 Rusnacko, J.; Gotfryd, D.; Chaloupka, J. Kitaev-like honeycomb magnets: Global phase behavior and emergent effective models. _Physical Review B_ 2019, _99_ , 064425
  * Choi et al. 2012 Choi, S. K.; Coldea, R.; Kolmogorov, A. N.; Lancaster, T.; Mazin, I. I.; Blundell, S. J.; Radaelli, P. G.; Singh, Y.; Gegenwart, P.; Choi, K. R.; Cheong, S.-W.; Baker, P. J.; Stock, C.; Taylor, J. Spin Waves and Revised Crystal Structure of Honeycomb Iridate Na2IrO3. _Physical Review Letters_ 2012, _108_ , 127204
  * Takahashi et al. 2019 Takahashi, S. K.; Wang, J.; Arsenault, A.; Imai, T.; Abramchuk, M.; Tafti, F.; Singer, P. M. Spin Excitations of a Proximate Kitaev Quantum Spin Liquid Realized in Cu2IrO3. _Physical Review X_ 2019, _9_ , 031047
  * Knolle et al. 2019 Knolle, J.; Moessner, R.; Perkins, N. B. Bond-Disordered Spin Liquid and the Honeycomb Iridate H3LiIr2O6: Abundant Low-Energy Density of States from Random Majorana Hopping. _Physical Review Letters_ 2019, _122_ , 047202
  * Kao et al. 2021 Kao, W.-H.; Knolle, J.; Halász, G. B.; Moessner, R.; Perkins, N. B. Vacancy-Induced Low-Energy Density of States in the Kitaev Spin Liquid. _Physical Review X_ 2021, _11_ , 011034
  * Frolov 2021 Frolov, S. Quantum computing’s reproducibility crisis: Majorana fermions. _Nature_ 2021, _592_ , 350–352
  * Field and Simula 2018 Field, B.; Simula, T. Introduction to topological quantum computation with non-Abelian anyons. _Quantum Science and Technology_ 2018, _3_ , 045004
  * You et al. 2012 You, Y.-Z.; Kimchi, I.; Vishwanath, A. Doping a spin-orbit Mott insulator: Topological superconductivity from the Kitaev-Heisenberg model and possible application to Na2/Li2IrO3. _Physical Review B_ 2012, _86_ , 085145, Publisher: American Physical Society


# Superconductivity in silicon: a review

Monika Moun  Goutam Sheet<EMAIL_ADDRESS>Department of Physical
Sciences, Indian Institute of Science Education and Research (IISER) Mohali,
Sector 81, S. A. S. Nagar, Manauli, PO: 140306, India.

###### Abstract

Silicon, one of the most abundant elements found on Earth, has been an
excellent choice of the semiconductor industry for ages. Despite it’s
remarkable applications in modern semiconductor-based electronic devices, the
potential of cubic silicon in superconducting electronics remained a challenge
because even heavily doped silicon crystals do not superconduct under normal
conditions. It is apparent that if superconductivity can be introduced in
cubic silicon, that will bring a breakthrough in low-dissipation electronic
circuitry. Motivated by this, attempts have been made by several research
groups to induce superconductivity in silicon through a number of different
routes. Some of the other structural phases of silicon like $\beta$-Sn and
simple hexagonal are, however, known to display superconductivity. In the
present review article, various theoretical and experimental aspects of
superconductivity in silicon are discussed. Superconductivity in different
phases and different structural forms of silicon are also reviewed. We also
highlight the potential of superconducting phases of silicon for technological
applications in super-conducting nano-electronics.

###### pacs:

73.23.Ad, 73.63.Rt, 74.10.+v, 74.45.+c, 74.78.Na

## I. Introduction

Silicon has been a material of choice for semiconductor technology for the
ultimate control that has been achieved on electrical properties.1; 2; 3; 4
Furthermore, silicon displays outstanding thermal stability and is abundantly
found on Earth. Functionality of silicon can be further enhanced by
incorporating other novel physical properties such as magnetism and
superconductivity.5; 6; 7; 8 The electrical conductivity of crystalline
silicon could be greatly altered in a control way by doping with both
electrons and holes, a dissipationless superconducting state has been a far
reaching property at reasonable temperature10; 9 Theoretical investigation
suggested that superconductivity can be obtained in doped elemental
semiconductors like silicon and germanium under certain conditions.11; 12; 13
Following such studies, experimentalists all over the world have been trying
to investigate a superconducting phase in silicon. It was seen that
superconductivity could possibly appear in silicon under very high pressure.14
The cubic diamond phase of silicon has also manifested superconductivity at
extreme conditions i.e. low temperatures and heavy doping beyond the
solubility limit.15; 16 Recently, a superconducting phase in moderately doped
crystalline silicon was studied under non-superconducting metallic point
contacts depicting a high critical temperature of 11 K.17 Besides the bulk
form of silicon, it’s different structural forms including silicene and nano-
structured silicon have shown potential for superconducting applications with
relatively higher transition temperatures. This review article provides an
overview of theoretical and experimental investigations in the emerging field
of superconductivity in different phases of silicon. The review article is
structured as follows: various phases of silicon including bulk
superconducting phases, superconductivity in doped silicon and it’s different
structural forms will be discussed. Superconducting phase under point-contacts
in doped silicon will also be reviewed followed by the author’s perspective
and potential applications of superconductivity in silicon.

## II. Various phases of silicon

The pioneering work of Minomura and Drickamer discussed various pressure
induced phase transitions in different semiconductors including silicon, GaAs,
GaSb, InP, InAs, AlSb and germanium.18 For all the semiconductors, a phase
transition to a metallic state was observed at high applied pressures of the
order of a few GPa. In silicon, the metallic phase was obtained at applied
pressure of nearly 19-20 GPa. While the cubic diamond phase of silicon has
been the most commonly investigated of all structures, it’s high-pressure
phases have also been examined both experimentally and theoretically.19; 20;
21 The phase diagram of silicon manifests several metastable phases also under
the ambient conditions (Figure 1). Under applied pressure of upto $\sim$ 11
GPa, the cubic diamond (Si-I) structure of bulk silicon is preserved.22 It
undergoes phase change from cubic diamond phase (Si-I) to body centered
tetragonal ($\beta$-Sn, Si-II) at applied pressure of slightly above 11 GPa.
The $\beta$-Sn(Si-II) phase changes to an orthorhombic Imma phase and
transforms to a simple hexagonal (sh, Si-V) at a pressure of about 15 GPa.
Structural phase change from simple hexagonal phase (sh) to a hexagonal closed
packed (hcp) occurs at a pressure around 40 GPa with an intermediate
orthorhombic Cmca phase (Si-VI) as revealed from x-ray diffraction studies and
ab-initio calculations.23 A face centered cubic (fcc) phase of silicon is
obtained at a pressure of 79 GPa, which remains stable upto the highest
investigated pressure of 248 GPa.24

Figure 1: Schematic representation of various low and high pressure structural
phases of silicon23, reproduced with permission, 2016, American Institute of
Physics

Interestingly, all the phases of silicon above pressure of nearly 12 GPa are
metallic in nature. This semiconductor-metal transition between cubic diamond
phase to high pressure metallic phase is irreversible. Apart from the stable
high-pressure silicon phases, a series of metastable phases has also been
produced as a result of releasing pressure instead of recovering the cubic
diamond phase under different conditions. When pressure is slowly released
from $\beta$-Sn(Si-II) phase, a metastable body centered cubic phase bc8 (Si-
III) has been observed with an intermediate phase called rhombohedral (r8)
phase.25 Two other metastable phases of silicon Si-VIII and Si-IX consisting
of large tetragonal cells have been obtained under rapid decompression. Si-IX
phase was examined to be the mixture of two structures having space groups
P42/m and P-4.26 Two more silicon phases were theoretically predicted, ST12
phase with distorted tetragonal cells27 and a body centered tetragonal (bct)
phase28. Figure 1 depicts various low and high pressure structural phases of
silicon.

## III. Bulk superconducting phases of silicon

In this section, we discuss the high pressure structural phases of silicon in
which superconductivity was observed. The theoretical findings of structural,
vibrational and electronic properties of silicon at high-pressure performed by
Chang and Cohen indicated that a pressure-sensitive soft-transverse-acoustic
phonon mode could result in the successful phase-transition from
$\beta$-Sn(Si-II) to simple hexagonal (sh) to hexagonal closed packed (hcp)
structure.29 The high density of states at the Fermi level of simple hexagonal
(sh) phase of pressurized silicon, and strong covalent bonding in the
transverse direction might lead to local field effects which could improve the
electron-electron interaction in sh silicon. This may facilitate the
appearance of a superconducting phase in the sh phase of silicon. The
observation of superconductivity in sh phase of Si was also suggested by Needs
and Martin from self-consistent density functional theory.30 Experimentally,
the resistance vs temperature ($R-T$) measurements down to liquid helium
temperatures were done by Wittig on $\beta$-Sn(Si-II) and a discontinuity in
the resistance-temperature curve was reported at 6.7 K.31 This was attributed
to a superconducting transition in the $\beta$-Sn(Si-II) phase of bulk
silicon. Il’ina and Istkevich observed a similar discontinuity at 7 K for
$\beta$-Sn(Si-II).32 Following these studies, superconductivity was also
theoretically verified and experimentally examined for $\beta$-Sn(Si-II) and
sh metallic phases of silicon at low temperatures and high pressures by Chang
and Cohen.33 For the $\beta$-Sn(Si-II) phase, the measured critical
temperature was 6.3 K at a pressure of 12 GPa which was lower than the
previous reports. Critical temperature was observed to go to maximum of 8.2 K
at a pressure of 15.2 GPa for the sh phase of silicon.

## IV. Superconductivity in boron doped silicon

In 2006, the discovery of superconductivity in doped cubic phase of silicon by
Bustarret et al. opened a new era of research on such superconducting
transition in doped silicon based systems.16 Until this discovery,
superconductivity had only been obtained in high-pressure structural phases of
silicon. Bustarret et al. showed that substitutional doping of boron achieved
using gas immersion laser doping technique beyond it’s equilibrium solubility
in silicon could induce superconductivity. Critical temperature and critical
field for the boron-doped (B concentration: 6 x 1020) silicon was found to be
0.35 K and 0.4 T, respectively as measured by electrical resistivity and
magnetic susceptibility measurements. Ab initio calculations were performed to
calculate electron-phonon coupling, superconducting transition temperature and
the microscopic pairing mechanism in B-doped silicon. The electron-phonon
coupling constant was found to be 0.28 which well matched with the previous
studies on pressurised silicon. The study also suggested that the observed
transition in silicon could be accounted by standard BCS mechanism and even
higher critical temperature could be possible via high boron doping in
silicon.

Following this experimental observation, Bourgeois et al. investigated the
vibrational and electronic modes along with electron- phonon coupling in cubic
phase of silicon with varying doping elements using a first-principles
approach.34 Their investigation provides evidences that the as-discovered
superconductiivity in highly doped silicon could be explicated on the basis of
a standard phonon-mediated BCS mechanism. In addition, authors suggested a
possible way to increase the critical temperature by one order of magnitude
via incorporating aluminium doping in silicon rather than boron.

## V. Superconductivity in silicene

The low-dimensional systems for applications in nanoscale superconducting
devices have come up recently as a possible alternative with the possibility
of higher $T_{c}$.35 In fact, 2D materials can be ideal candidates for such
application as it has been predicted theoretically that in highly doped 2D
materials, critical temperature could go higher than the boiling point of
liquid nitrogen.36 Silicene is the monolayer counterpart of silicon having
honeycomb lattice structure like graphene and has proved it’s potential to be
compatible with current nano-technology. Several semiconductors from group IV
have been found to be superconductors upon heavy boron doping. In these strong
covalent materials such as silicon and diamond, charge doping occurs into
$sp^{3}$-bonding bands near the Fermi level having large electron-phonon
coupling.37 Considering the case of 2D materials, for instance, pristine
graphene has $sp^{2}$ bonding but hole-doped graphene favors $sp^{3}$ bonding
and thus, enhanced electron-phonon coupling.36 Silicene can exhibit even
stronger electron-phonon coupling owing to low-buckled honeycomb structure
which may lead to partial $sp^{3}$ hybridization.38 In that case, silicene may
also superconduct.

Figure 2: (a) Conductance spectra acquired on silicene and silver substrate
(Ag (111)), respectively. (b) Conductance spectra curves taken along silicene
film at various positions, exhibiting the uniformity of the superconducting
gap. (c) Temperature dependent conductance spectra on silicene film.40,
reproduced with permission, 2013, American Institute of Physics

Experiments show that the charge carriers in silicene on Ag (111) behave as
massless Dirac Fermions.39 Following this study, a superconducting gap was
observed experimentally using scanning tunneling spectroscopy in silicene on
Ag(111) surface with the relatively higher $T_{c}$ values ranging from 35 – 40
K.40 The superconductor-metal transition was investigated by temperature-
dependent transport and tunneling experiments shown in figure 2. The gap
observed in spectrum appeared precisely at the Fermi energy strongly hinting
to the possibility of the gap being a superconducting gap. The amplitude of
the gap was measured to be 35 meV. It was proposed that the charge transfer
between Ag and silicene along with strong electron-phonon coupling led to the
high $T_{C}$ superconductivity in silicene/Ag. However, the knowledge of the
superconducting pairing mechanism is still lacking, hence the mechanism behind
superconductivity in silicene remains an open problem. Following the
experimental observations, Wan et al. theoretically studied the
superconductivity in silicene upon application of strain using first
principles calculations.38 The group suggested that electron-doped silicene
could serve as a good 2D superconductor with increased density of states under
biaxial tensile strain and the critical temperature could go higher above 10 K
by 5 % applied strain. The enhancement in the transition temperature was
attributed to the strong interaction between out of plane phonon modes
transverse to the silicene plane and additional strain induced electronic
states near the Fermi level.

## VI. Superconductivity in nanostructured silicon

$p$-type silicon quantum wells on $n$-type silicon have also shown potential
towards enhanced superconducting properties in silicon. Bagraev et al.
fabricated heavily boron doped $p$-type silicon quantum wells on the $n$-type
Si(100) with high mobility.41; 42 The group proposed that the formation of
oxide overlayers on silicon surfaces could promote the generation of the self-
interstitial sites and vacancies known as “microdefects”. Longitudinal and
lateral quantum wells between microdefect layers were formed depending on the
oxide layer thickness. Although both longitudinal and lateral silicon quantum
wells could be used as electrical and optical active microcavities for
nanoelectronic and optoelectronic applications, the dangling bonds present at
the interfaces could hinder the practical applications. Therefore, short-time
boron diffusion has been suggested for the passivation of silicon surface as
well as the created defects during oxidation process of the silicon wafers.
This aids the transformation of the microdefect arrays into $\delta$-barriers
as observed from the infrared Fourier spectroscopy, cyclotron resonance and
scanning tunneling microscopy (STM) techniques. Possible superconducting gap
of 0.022 eV and high superconducting transition temperature of 150 K was
observed.43

Figure 3: (a) Proposed mechanism for the boron dipole centers (b) Dependence
of temperature and magnetic field on static magnetic susceptibility in
$\delta$ barriers confined silicon quantum wells.43, reproduced with
permission, 2006, Elsevier (c-d) Multiple Andreev reflection (MAR) signatures
with the (c) forward bias and (d) reverse bias applied to the silicon quantum
wells47, reproduced with permission, 2008, Elsevier

The same group investigated the detailed hole transport in the silicon
sandwich nanostructures using tunneling spectroscopy.44 They showed that the
multiple Andreev reflection features of holes in the silicon quantum wells
could provide the mechanism behind the proximity effect phenomenon in silicon
sandwich nanostructures. Moreover, the oscillations of current-voltage
characteristics at temperatures above and below $T_{c}$ were identical
reflecting the presence of interaction between the single hole tunneling and
Cooper pairs. Their observation could lead to controlled fabrication of nano-
dimensional Josephson junctions and hybrid systems suggesting new
possibilities in the field of superconducting nanoelectronics.45; 46

## VII. Superconductivity under point contacts on heavily doped silicon

Recently, evidence of superconductivity was found in moderately doped silicon
crystals under point contacts with non-superconducting metallic tips with
reasonably high critical temperature($\sim 11K$).17

Figure 4: (a) $R-T$ plot depicting superconducting transition for point
contacts in intermediate regime shown in (c). Conductance spectra obtained on
As-doped Si(111) showing (b) ballistic regime (c) intermediate regime (d)
thermal transport regime (e) conductance spectra for thermal regime depicting
multiple critical current peak features. (g) $R-T$ curves at different
magnetic fields (f) $H–T$ phase diagram obtained from (g).17, reprinted with
permission, 2018, American Institute of Physics

The superconducting phase was examined in As-doped Si (111) using point
contact spectroscopy and different transport regimes i.e ballistic,
intermediate and thermal regime were investigated. The comparative estimate of
the contact diameter with the mean free path of electrons decide the transport
regimes.48 Ballistic transport regime is said to occur in the condition when
the contact diameter becomes smaller than electron’s elastic mean free path.
On the other hand, in thermal regime, the contact diameter is larger than the
electron’s inelastic mean free path. For a point-contact between metallic tip
and the sample, according to the Wexler’s formula, total resistance is sum of
both thermal or Maxwell’s resistance $(R_{M}=\rho(T)/2r)$ and ballistic or
Sharvin’s resistance $(R_{S}=\frac{2h/e^{2}}{{(rk_{F})}^{2}})$ where
$\rho(T)$, r, $k_{F}$, h and e represent the resistivity of the sample, radius
of the tip, Fermi radius, Planck’s constant and electronic charge,
respectively. In the intermediate regime, the contact diameter is such that
both ballistic resistance and thermal resistance become comparable. The
conductance spectra reveals different signatures for each kind of transport.
If the conductance spectra ($dI/dV$ vs $V$) shows symmetric two peaks about
V=0, it is the signature of Andreev reflection, known to appear in ballistic
regime of superconducting point contacts.49 Next comes the ”intermediate
regime” where the conductance spectrum displays distinct dips accompanying
suppressed Andreev reflection dominated symmetric double peak feature. In the
thermal regime, strong conductance dips are observed with no signature of
Andreev reflection associated conductance peaks. The spectroscopic data shown
in figure 5 confirms the presence of a superconducting phase in doped-
crystalline silicon under point contacts. The superconducting phase was
obtained with a considerable high $T_{c}$ 11 K and superconducting energy gap
of 2 meV.17 shown in figure 4. Subsequently, the superconducting phase under
point contacts was examined on both $n$-doped (As dopants) and $p$-doped
(B-dopants) facets of silicon crystals using various metallic tips exhibiting
universal behavior of mesoscopic superconductivity in silicon crystals.50
Before this, point-contact spectroscopy was been used to understand
semiconductor-metal junctions by several groups. Tanaka et al. suggested that
superconductivity in silicon could be obtained in rapidly quenched granules of
silicon without applying high pressure.51 Conductance spectra was acquired in
both Au-Si and Si-Si point contact junctions. It was proposed that the
observation of symmetric negative $dI/dV$ spike features at non-zero voltages
was due to the growing quenched metallic phase granules (figure 5). These
granules led to lattice condition similar to the high pressure superconducting
phase of silicon and possibly created a microjunction network in the vicinity
of point contact. This study followed the similar investigations performed in
different materials such as GaAs, PbSe etc. In the model as proposed by Shaw,
the sample fabrication process of point contacts and $p-n$ junctions involves
heat treatment which may create a micro-Josephson junction’s network.52 These
consist of superconducting granules in the boundary of the tunnel junction,
with the condition of presence of superconducting element in the materials
comprising of the junction (for example, Pb in PbSe junction, Ga in GaP).

Figure 5: (a) Conductance spectra of Au-Si point contact junctions (b)
Conductance spectra of silicon-silicon point contact junction.51, reproduced
with permission, 1974, Elsevier

## VIII. Future outlook

Superconducting devices have potential applications in various fields ranging
from sensors and detectors to quantum information technology. However, all
superconducting future advancements will require complex framework and the
adopted technology must be perfectly managed in terms of reliability and
reproducibility. In this quest, superconducting silicon can address the
challenges owing to the established fabrication of high quality silicon based
functional devices. Superconductivity in silicon can provide an avenue towards
designing superconducting nanomechanical resonators and Superconducting
Quantum Interference Device (SQUID) on chip that could work as electro-
magnetic sensors. Duvauchelle et al. made an attempt to form SQUID from a
monolayer superconducting silicon fabricated using the GILD technique.53 The
device consists of two nano-sized bridges and in the fabricated device,
magnetic flux could be modulated at low values of magnetic field and
temperature. More practical devices can be implemented in future such as
Josephson Field Effect Transistors comprising of superconducting drain and
source contacts and current flowing through a channel modulated by an
electrostatic gate. Josephson junction comprising of two superconducting
materials with an insulating layer sandwiched between them, is the basic
building block of most of the superconducting devices. Superconductor-metal-
superconductor junction, an alternative to Josephson junctions, was proposed
theoretically by De Gennes54 and experimentally investigated by Clarke.55 When
a normal metal and superconductor are electrically contacted, Cooper pairs can
penetrate from the superconducting side to the normal metal, due to the
proximity effect, a commonly known phenomenon in superconductivity.56
Researchers have studied a special category of this kind of junctions by
replacing the metal by a degenerately doped semiconductor.57 The advantage of
the superconductor-semiconductor-superconductor junctions lies in the
modulation of the coupling strength between the superconductors by changing
the doping/gate voltage and hence, carrier concentration in the semiconductor,
which results in the variable supercurrent. In the case of semiconductor-
superconductor junction, the Schottky barrier present at the semiconductor-
superconductor interface limits the seamless flow of Cooper pairs from
superconductor side to semiconductor side. In this respect,
silicon(semiconductor)-silicon(superconductor) junctions can provide Schottky
barrier free contacts. This can be made possible using the gas immersion laser
doping (GILD) technique in which highly doped regions can be fabricated with
abrupt interfaces. The ultimate control over supercurrent by varying the
doping level in the semiconductor can lead to applications such as low power,
ultrafast electronic devices. Another important application of
superconductivity in silicon can find in energy saving devices by reducing
energy losses to a considerable extent. For modern day technology, there
exists a huge requirement of various energy saving devices by limiting the
energy losses to an extent with high efficiency. Especially, superconducting
silicon can be a potential material for such devices including superconducting
electric motors and generators, magnetic energy storage systems and so
forth.58

## IX. Conclusions

On account of the experimental and theoretical research efforts, the field of
superconductivity in various structural forms of silicon has progressed
enormously. Although advancement in structural phases prediction in silicon
has identified a number of superconducting phases, an effective route towards
discovering superconductivity in moderately doped silicon is required. On the
other hand, discovery of superconductivity with high critical temperature is
crucial if the promise of silicon as superconductor is to be realized for
industrial purpose. Indeed, different structural forms of silicon have shown
promising superconducting properties, but thorough experimental investigation
is lacking. For example, nano-dimensional silicon quantum wells and Josephson
junctions have presented the possibility of improved superconducting
properties with relatively higher critical temperature. Silicene, the
monolayer conterpart of silicon, can also serve as a promising candidate for
miniaturized future superconductor technology due to strong electron-phonon
coupling as predicted theoretically and verified experimentally with the
possibility to achieve higher critical temperatures. For the fundamental
understanding of superconductivity in silicon at nanoscale, point contact
spectroscopy is an effective tool. Tip-induced superconductivity in silicon
with different kinds of non-superconducting metallic probes can open a way
towards more fundamental understanding of superconductivity in silicon with
higher critical temperature. At last, an attempt is made to provide an
overview of the developments and ongoing research in the field of
superconductivity in silicon. The future research in this field will lead to
advanced understanding of such exotic superconducting properties of silicon,
especially towards high critical temperature to make this material more
feasible for superconducting device applications.

## Acknowledgements

GS acknowledges the support of Swarnajayanti fellowship awarded by the
Department of Science and Technology (DST), Govt. of India (grant number
DST/SJF/PSA01/2015-16). MM acknowledges the support of NPDF-SERB fellowship
(PDF/2020/002122).

## References

## References

  * (1) W. N. Ye, and Y. Xiong, Review of silicon photonics: history and recent advances. Journal of Modern Optics 60, 1299-1320 (2013).
  * (2) D. H. Kim, J-H Ahn, W. M. Choi, H-S Kim, T-H Kim, J. Song, Y. Y. Huang, Z. Liu, C. Lu, and J. A. Rogers, Stretchable and foldable silicon integrated circuits.Science 320, 507 (2008).
  * (3) M. Ieong, B. Doris, J. Kedzierski, K. Rim, and M. Yang, Silicon device scaling to the sub-10-nm regime. Science 306, 2057 (2004).
  * (4) J. E. Duvauchelle, A. Francheteau, C. Marcenat, F. Chiodi, D. Debarre, K. Hasselbach, J. R. Kirtley, and F. Lefloch, Silicon superconducting quantum interference device. Appl. Phys. Lett. 107, 072601 (2015).
  * (5) B. -C. Min, K. Motohashi, C. Lodder, and R. Jansen, Tunable spin-tunnel contacts to silicon using low-work-function ferromagnets Nature Materials 5, 817-822 (2006).
  * (6) M. P. Delmo, S. Yamamoto, S. Kasai, T. Ono, and K. Kobayashi, Large positive magnetoresistive effect in silicon induced by the space-charge effectNature 457, 1112 (2009).
  * (7) D. Erskine, P. Y. Yu, K. J. Chang, and M. L. Cohen, Superconductivity and Phase Transitions in Compressed Si to 45 GPaPhys. Rev. Lett. 57, 2741-2744 (1986).
  * (8) R. J. Cava, Super siliconNature 444, 427 (2006).
  * (9) P. Dai, Y. Zhang, and M. P. Sarachik, Critical conductivity exponent for Si:B Phys. Rev. Lett. 66, 1914-1917 (1991).
  * (10) C. Kittel, Introduction to Solid State Physics Wiley, (2004).
  * (11) M. L. Cohen, Superconductivity in Many-Valley Semiconductors and in Semimetals* Phys. Rev. 134, A511-A522 (1964).
  * (12) D. Connétable, V. Timoshevskii, B. Masenelli, J. Beille, J. Marcus, B. Barbara, A. M. Saitta, G.-M. Rignanese, P. Mélinon, S. Yamanaka, and X. Blase, Superconductivity in Doped sp3 Semiconductors: The Case of the Clathrates Phys. Rev. Lett. 91, 247001 (2003).
  * (13) L. Boeri, J. Kortus, and O. K. Andersen, Three-Dimensional MgB2-Type Superconductivity in Hole-Doped Diamond. Phys. Rev. Lett. 93, 237002 (2004).
  * (14) T. H. Lin, W. Y. Dong, K. J. Dunn, and C. N. J. Wagner,Pressure-induced superconductivity in high-pressure phases of Si Phys. Rrv. B 33, 7820 (1986).
  * (15) E. Antoncik, The influence of the solubility limit on diffusion of phosphorus and arsenic into silicon Appl. Phys. A 58, 117-123 (1994).
  * (16) E. Bustarret, C. Marcenat, P. Achatz, J. Kačmarčik, F. Lévy, A. Huxley, L. Ortéga, E. Bourgeois, X. Blase, D. Débarre, and J. Boulmer Superconductivity in doped cubic silicon Nature 444, 465-468 (2006).
  * (17) A. Sirohi, S. Gayen, M. Aslam, and G. Sheet, Mesoscopic superconductivity above 10 K in silicon point contacts Appl. Phys. Lett. abc, 12345 (2018).
  * (18) S. Minomura, and H. G. Drickamer, Pressure Induced Phase Transitions in silicon, Germanium and some III-V compounds J. Phys. Chem. Solids Pergamon Press 23, 451-456 (1962).
  * (19) H. Olijnyk, S. K. Sikka, and W. B. Holzapfel, Structural phase transitions in Si and Ge under pressures up to 50 GPa Phys. Lett. 103, 137-140 (1984).
  * (20) J. Z. Hu, and I. L. Spain, Phases of silicon at high pressure Solid State Commun. 51, 263 (1984).
  * (21) S. J. Duclos, Y. K. Vohra, and A. L. Ruoff, hcp to fcc transition in silicon at 78 GPa and studies to 100 GPa Phys. Rev. Lett. 58, 775 (1987).
  * (22) A. Mujica, A. Rubio, A. Munõz, and R. J. Needs, High-pressure phases of group-IV, III–V, and II–VI compounds Rev. Mod. Phys. 75, 863 (2003).
  * (23) S. Wippermann, Y. He, M. Vorovos, and G. Galli, Novel silicon phases and nanostructures for solar energy conversion Appl. Phys. Rev. 3, 040807 (2016).
  * (24) S. J. Duclos, Y. K. Vohra, and A. L. Ruoff, Experimental study of the crystal stability and equation of state of Si to 248 Gpa Phys. Rev. B 41, 17 (1990).
  * (25) R. H. Wentorf, and J. S. Kasper, Two new forms of silicon Science 139, 338 (1963).
  * (26) M. C. Nguyen, X. Zhao, Y. Wanga, C.-Z. Wang, and K.-M. Ho, Genetic algorithm prediction of crystal structure of metastable Si-IX phase Solid state communications 182, 14-16 (2014).
  * (27) S. J. Clark, G. J. Ackland, and J. Crain, Theoretical study of high-density phases of covalent semiconductors. II. Empirical treatment Phys. Rev. B 49, 5341 (1994).
  * (28) Y. Fujimoto, T. Koretsune, S. Saito, T. Miyake, and A. Oshiyama, A new crystalline phase of four-fold coordinated silicon and germanium New J. Phys.  10, 083001 (2008).
  * (29) K. J. Chang and M. L. Cohen, Structural and electronic properties of the high-pressure hexagonal phases of Si Phys. Rev. B 30, 9 (1984).
  * (30) R. J. Needs, R. M. Martin, Transition from $\beta-tin$ to simple hexagonal silicon under pressure Phys. Rev. B 30, 9 (1984).
  * (31) J. Wittig, For the superconductivity of germanium and silicon under high pressure magazine for physics 195, 215-227 (1966).
  * (32) G. Martinez, J. M. Mignot, G. Chouteau, K. J. Chang, M. M. Dacorogna, and M. L. Cohen, Superconductivity of SiliconPhysica Scripta. T13, 226-229 (1986).
  * (33) K. J. Chang, M. M. Dacorogna, M. L. Cohen, J. M. Mignot, G. Chouteau and G. Martinez, Superconductivity in High-Pressure Metallic Phases of Si Phys. Rev. Lett. 54, 21 (1985).
  * (34) E. Bourgeois and X. Blasea, Superconductivity in doped cubic silicon: An ab initio study Appl. Phys. Lett. 90, 142511 (2007).
  * (35) X. Blase, E. Bustarret, C. Chapelier, T. Klein, and C. Marcenat, Superconducting group-IV semiconductors Nature Mater. 8, 375 (2009).
  * (36) G. Savini, A. C. Ferrari, and F. Giustino, First-principles prediction of doped graphane as a high-temperature electron-phonon superconductor Phys. Rev. Lett. 105, 037002 (2010).
  * (37) K. W. Lee and W. E. Pickett, Superconductivity in boron-doped diamond Phys. Rev. Lett. 93, 237003 (2004).
  * (38) W. Wan, Y. Ge, F. Yang and Y. Yao, Phonon-mediated superconductivity in silicene predicted by first-principles density functional calculations EPL 104, 36001 (2013).
  * (39) B. Feng, Z. J. Ding, S. Meng, Y. G. Yao, X. Y. He, P. Cheng, L. Chen, and K. H. Wu, Evidence of silicene in honeycomb structures of silicon on Ag (111) Nano Lett. 12, 3507 (2012).
  * (40) L. Chen, B. Feng, and K. Wua, Observation of a possible superconducting gap in silicene on Ag(111) surface Appl. Phys. Lett. 102, 081602 (2013).
  * (41) N. T. Bagraev, A. Shelykh, V. K. Ivanov and L. E. Klyachkin, Spin depolarization in quantum wires polarized spontaneously in zero magnetic field Phys. Rev. B 70, 155315 (2004).
  * (42) N. T. Bagraev, L. E. Klyachkin, A. A. Koudryavtsev, A. M. Malyarenko, and V. V. Romanov, Superconducting Properties of Silicon Nanostructures Semiconductors 43, 1441-1454 (2009).
  * (43) N. T. Bagraev, W. Gehlhoff, L. E. Klyachkin, A. M. Malyarenko, V. V. Romanov, and S. A. Rykov, Superconductivity in silicon nanostructures Physica C 437, 21-24 (2006).
  * (44) N. T. Bagraev, L. E. Klyachkin, A. A. Koudryavtsev, A. M. Malyarenko, G. A. Oganesyan, and D. S. Poloskin, Quantum Supercurrent and Andreev Reflection in Silicon Nanostructures Semiconductors 43, 1455–1465 (2009).
  * (45) C. W. J. Beenakker and H. van Houten, Josephson current through a superconducting quantum point contact shorter than the coherence length Phys. Rev. Lett. 66, 3056 (1991).
  * (46) P. Jarillo-Herrero, J. van Dam, and L. Kouwenhoven, Quantum supercurrent transistors in carbon nanotubes Nature 439, 953-956 (2006).
  * (47) N. T. Bagraev, W. Gehlhoff, L. E. Klyachkin, A. A. Kudryavtsev, A. M. Malyarenko, G. A. Oganesyan, D. S. Poloskin, and V. V. Romanov, Spin-dependent transport of holes in silicon quantum wells confined by superconductor barriers Physica C 468, 840–843 (2008).
  * (48) G. E. Blonder, M. Tinkham, and T. M. Klapwijk, Transition from metallic to tunneling regimes in superconducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion Phys. Rev. B 25, 4515-4532 (1982).
  * (49) G. Sheet, S. Mukhopadhyay, and P. Raychaudhuri, Role of critical current on the point-contact Andreev refection spectra between a normal metal and a superconductor Phys. Rev. B 69, 134507 (2004).
  * (50) M. Moun, A. Sirohi, and G. Sheet, Universality of Interfacial Superconductivity in Heavily Doped Silicon ACS Appl. Electron. Mater. 3, 1594-1600 (2021).
  * (51) K. Tanaka, Possibility of superconductivity in rapidly quenched silicon Physics Letters 47A, 4 (1974).
  * (52) R. W. Shaw, Superconducting Origin of Spikes in the Differential Conductance of p‐n and Point Contact Tunnel Junctions Journal of Applied Physics 40, 4464 (1969).
  * (53) J. E. Duvauchelle, A. Francheteau, C. Marcenat, F. Chiodi, D. Débarre, K. Hasselbach, J. R. Kirtley, and F. Lefloch, Silicon superconducting quantum interference device Applied Phys. Lett. 107, 072601 (2015).
  * (54) P. G. de Gennes, Boundary Effects in Superconductors Rev. Mod. Phys. 36, 225 (1964).
  * (55) J. Clarke, Supercurrents in lead-copper-lead sandwiches Proc. R. Soc. Lond. Ser. A: Math. Phys. Sci. 308, 447 (1969).
  * (56) T. M. Klapwijk, Proximity Effect From an Andreev Perspective Journal of Superconductivity: Incorporating Novel Magnetism 17, 5 (2004).
  * (57) T. Schäpers, Superconductor/Semiconductor Junctions Springer (2001).
  * (58) T. R. Schneider, Energy applications of superconductivity Annu. Rev. energy environ 16, 533-555 (1991).


# Sum rules for energy deposition eigenchannels in scattering systems

Alexey Yamilov  Physics Department, Missouri University of Science &
Technology, Rolla, Missouri, USA  Corresponding author<EMAIL_ADDRESS>Nicholas Bender  School of Applied and Engineering Physics, Cornell
University, Ithaca, New York, USA  Hui Cao  Department of Applied Physics,
Yale University, New Haven, Connecticut, USA

###### Abstract

In a random-scattering system, the deposition matrix maps the incident
wavefront to the internal field distribution across a target volume. The
corresponding eigenchannels have been used to enhance the wave energy
delivered to the target. Here we find the sum rules for the eigenvalues and
eigenchannels of the deposition matrix in any system geometry: including two
and three-dimensional scattering systems, as well as narrow waveguides and
wide slabs. We derive a number of constraints on the eigenchannel intensity
distributions inside the system as well as the corresponding eigenvalues. Our
results are general and applicable to random systems of arbitrary scattering
strength as well as different types of waves including electromagnetic waves,
acoustic waves, and matter waves.

††journal: ol

## Introduction

Wavefront shaping opened a new frontier for coherent manipulation of wave
propagation in complex media [1, 2, 3, 4]. The core idea is rooted in the
determinism of coherent wave propagation in static, linear scattering media
[5]. Although the incident wavefront can be optimized via an iterative
optimization procedure [6], solving the eigenvalue problem of a linear
operator (matrix) is a more predictive approach [7]. Recently, a number of
matrices have been introduced in order to manipulate quantities such as
transmittance, reflectance, dwell-time, spatial distribution etc, see Refs.
([8, 5, 9, 4]) for review. Control and optimization of wave energy inside a
scattering medium requires an access to the internal field distribution [10].
To find the ultimate limit of energy delivery to a target buried deep inside a
diffusive medium, we recently introduced a deposition matrix that relates the
incident wavefront to the internal field distribution across the target [11].
The maximal eigenvalue of the matrix gives the largest possible energy
enhancement, while the corresponding eigenvector gives the optimal incident
wavefront. However, the properties of deposition eigenvalues and eigenchannels
remain essentially unknown.

In this work, we theoretically obtain a series of sum rules satisfied by the
deposition eigenchannels and their eigenvalues. In Ref. [11] a deposition
eigenchannel is decomposed by the transmission eigenchannels to illustrate the
incoherent (i.e. intensity summation) and coherent (interference)
contributions. Here, we demonstrate that these contributions also obey
rigorous sum rules. These relationships not only provide physical insights
into the deposition eigenchannels, but also play an important role in
utilizing such channels for the most efficient energy delivery in experiments.

Figure 1: Schematic depiction of the slab (a,b) and the waveguide scattering-
medium geometries. The deposition matrix relates the input to the internal
fields inside a scattering medium, panels (a,c). The transmission matrix
$\mathcal{T}_{mn}$ relates the input and output degrees of freedom, panels
(b,d).

## Deposition matrix

Fig. 1 illustrates two geometries we consider for energy delivery into a
linear scattering system: (a,b) wide slab with open (leaky) boundary, (c,d)
narrow waveguide with closed (reflecting) sidewall. The deposition matrix (DM)
$\mathcal{Z}$ is introduced for a target region of arbitrary size, shape, and
depth [11]. It relates an orthonormal set of input waves to the corresponding
spatial field distributions within the target region. As shown in Fig. 1(a) of
a wide slab, the orthogonal input modes have distinct wavevectors. While in
Fig. 1(c) the input waves are a complete set of waveguide modes. The total
number of input modes is $N$. With a unit flux of incident light in the $n$-th
mode, the complex field distribution throughout the scattering system is
$E^{(0)}_{n}(\vec{r})$. We sample the fields uniformly across a target region
of volume $\mathcal{V}$ centered at $\vec{r}_{D}$, see Figs. 1(a). The field
at the $m$-th sampling point $\vec{r}_{m}$ is
$E^{(0)}_{n}(\vec{r}_{m};\vec{r}_{D})$, where $m=\\{1,...,M\\}$. The volume
$\mathcal{V}/{M}$ covered by each sampling point is much smaller than
$\lambda^{3}$, where $\lambda$ is the wavelength. The deposition matrix of
dimension $M\times N$ is defined as

$\mathcal{Z}_{mn}(\vec{r}_{D})\equiv\left[\epsilon(\vec{r}_{m})\,\frac{\mathcal{V}}{M}\right]^{1/2}\,E^{(0)}_{n}(\vec{r}_{m};\vec{r}_{D})\,,$
(1)

where $\epsilon(\vec{r})$ is the spatially-varying dielectric constant.

The deposition eigenchannels are obtained from the singular value
decomposition (SVD),
$\mathcal{Z}_{mn}(\vec{r}_{D})=\sum_{\alpha=1}^{N}U_{m\alpha}^{(D)}(\vec{r}_{D})\,\zeta_{\alpha}^{1/2}(\vec{r}_{D})\,\left[V_{\alpha
n}^{(D)}(\vec{r}_{D})\right]^{*}$. The incident wavefront of the $\alpha$-th
eigenchannel is
$E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})=\sum_{n=1}^{N}E^{(0)}_{n}(\vec{r})\,V_{n\alpha}^{(D)}(\vec{r}_{D})$.

$\zeta_{\alpha}(\vec{r}_{D})$ and $V_{n\alpha}^{(D)}(\vec{r}_{D})$ are the
eigenvalue and eigenvector of matrix
$\mathcal{Z}^{\dagger}(\vec{r}_{D})\mathcal{Z}(\vec{r}_{D})$, respectively.
Using the unitarity of matrices $U^{(D)}(\vec{r}_{D})$ and
$V^{(D)}(\vec{r}_{D})$ as well as the definition of $\mathcal{Z}(\vec{r}_{D})$
in Eq. (1), we get an explicit relationship for the eigenvalue [12]

$\zeta_{\alpha}(\vec{r}_{D})=\frac{\mathcal{V}}{M}\,\sum_{m=1}^{M}\epsilon(\vec{r}_{m})\left|E_{\alpha}^{(D)}(\vec{r}_{m};\vec{r}_{D})\right|^{2}.$
(2)

This relationship reveals that the deposition eigenvalue is equal to the total
energy inside the target region with a coherent excitation of
$E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})$.

For comparison, Fig. 1(b,d) shows the transmission matrix $\mathcal{T}$ that
maps the incident fields to the transmitted fields. When a monochromatic light
with a unit flux in the $n$-th mode is incident to the scattering system,
$\left|\mathcal{T}_{mn}\right|^{2}$ is the amount of flux carried away by the
$m$-th outgoing mode in transmission. The transmission eigenchannels are
obtained from the SVD
$\mathcal{T}_{mn}=\sum_{\alpha=1}^{N}U_{m\alpha}^{(T)}\cdot\tau_{\alpha}^{1/2}\cdot\left[V_{\alpha
n}^{(T)}\right]^{*}$. The input and output wavefronts of the $\alpha$-th
transmission eigenchannel are given by $V_{\alpha n}^{(T)}$ and
$U_{m\alpha}^{(T)}$, and the transmittance by eigenvalue $\tau_{\alpha}$.

## Sum rules for deposited energy

Singular value decomposition is a linear transformation of the input basis
that yields an orthonormal set defined by the unitary matrix
$V^{(D)}(\vec{r}_{D})$ or $V^{(T)}$ for the deposition and transmission,
respectively. This orthonormal property of the input eigenvectors
$V_{n\alpha}^{(D)}(\vec{r}_{D})$ and $V_{n\alpha}^{(T)}$ has a profound impact
on the distribution of the wave intensity inside the scattering medium.
Summing the intensity of all eigenchannels at any given point $\vec{r}$ gives
[12]

$\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})\right|^{2}=\sum_{n=1}^{N}\left|E_{\alpha}^{(T)}(\vec{r})\right|^{2}=\sum_{n=1}^{N}\left|E^{(0)}_{n}(\vec{r})\right|^{2}.$
(3)

The sum of intensities by subsequently exciting the system with individual
deposition or transmission eigenchannels is equal to that with any orthogonal
set of input modes. The above relationships lead to two remarkable properties.
First, the physical quantity being preserved by SVD transformation is the
intensity, and not the field. Because Eq. (3) holds at every position, and
thus can be multiplied by $\epsilon(\vec{r})$ from both sides, it can be
interpreted as point-wise (i.e. local) conservation of energy. Secondly, the
above relationship is not statistical – it does not involve any statistical
averaging over an ensemble of disorder realizations, instead, it holds for
every realization.

Figure 2: (a) Numerically computed probability density of the deposition (red
circle) and transmission (green cross) eigenvalues for a slice target at the
depth $z_{D}/L=1/2$ in a 2D disordered waveguide [12]. (b) Intensity
participation ratio for deposition and transmission eigenchannels at different
depths, see text.

The sum rule in Eq. (3) can be exploited to estimate the intensity enhancement
when combined with a known probability density function (PDF) of the
corresponding eigenvalues. To illustrate this point, lets consider
transmission eigenchannels of a diffusive system and use a numerical
simulation of 2D disordered waveguide in Fig. 1(d). The PDF of the
transmission eigenvalues in this case is the celebrated bimodal distribution
[13], predicting that only $g\ll N$ eigenchannels have $\tau_{\alpha}\sim 1$
while the rest have $\tau_{\alpha}\sim 0$, c.f. Fig. 2a, where $g$ is the
dimensionless conductance. In contrast, for such system, the PDF of deposition
eigenvalues has been predicted in Ref. [11] to exhibit a long tail toward the
largest value of $\zeta_{max}(\vec{r}_{D})/\langle\zeta(\vec{r}_{D})\rangle$,
see Fig. 2a. Consequently, fewer deposition eigenchannels with very large
$\zeta(\vec{r})$ provide a higher enhancement at the target region, compared
to the high-transmission eigenchannels. This effect can be quantified using
the intensity participation ratio ${\cal
P}^{(D)}(z_{D})=\left(\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(D)}(z_{D};z_{D})\right|^{2}\right)^{2}/\left(\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(D)}(z_{D};z_{D})\right|^{4}\right)$
for deposition eigenchannels, and ${\cal
P}^{(T)}(z)=\left(\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(T)}(z)\right|^{2}\right)^{2}/\left(\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(T)}(z)\right|^{4}\right)$
for transmission eigenchannels [14]. Fig. 3b indeed shows a lower
participation ratio ${\cal P}^{(D)}$ than ${\cal P}^{(T)}$, reflecting a
smaller number of deposition eigenchannels with major contributions. This, in
combination with the sum rule in Eq. (3), leads to a higher cross-section
integrated intensity at the target depth $z_{D}$, see also Fig. 3a.
Furthermore, the extraordinarily large intensities of a few deposition
eigenchannels, in turn, allow for the minimal intensity to be even below that
of the lowest-transmission eigenchannel, c.f. Fig. 3b. All of these
conclusions are indeed supported by Ref. [11], demonstrating that a
combination of the knowledge of the eigenvalue PDF and the constraints imposed
by Eq. (3) has important implications for local energy density enhancement by
the eigenchannels.

Figure 3: Deposition eigenchannels (red) in disordered waveguide geometry [12]
with $\alpha=1,N$ for $z_{D}/L=1/2$ and their incoherent components (blue) are
compared to the corresponding TEs (green). Random input intensity profile
(dashed line) is shown for reference. Around the target region, the
enhancement/suppression of the local intensity of the DEs above/below that of
the maximum/minimum transmission eigenchannels is observed. The incoherent
part, on the other hand, is below/above the maximum/minimum transmission
eigenchannels.

We also find the sum rule for the deposition eigenvalues. Summing over
$\zeta_{\alpha}(\vec{r}_{D})$ and recalling Eq. (3), we obtain

$\displaystyle\sum_{n=1}^{N}\zeta_{\alpha}(\vec{r}_{D})$ $\displaystyle=$
$\displaystyle\frac{\mathcal{V}}{M}\,\sum_{m=1}^{M}\epsilon(\vec{r}_{m})\left[\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(D)}(\vec{r}_{m};\vec{r}_{D})\right|^{2}\right]$
(4) $\displaystyle=$
$\displaystyle\frac{\mathcal{V}}{M}\,\sum_{m=1}^{M}\epsilon(\vec{r}_{m})\left[\sum_{\alpha=1}^{N}\left|E_{n}^{(T)}(\vec{r}_{m})\right|^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{\mathcal{V}}{M}\,\sum_{m=1}^{M}\epsilon(\vec{r}_{m})\left[\sum_{n=1}^{N}\left|E_{n}^{(0)}(\vec{r}_{m})\right|^{2}\right]$

This equation shows that the sum of all deposition eigenvalues is equal to sum
of local energy (sampled over $M$ points $\left\\{\vec{r}_{m}\right\\}$)
excited by all deposition or transmission eigenchannels or any orthogonal set
of input waves.

Therefore, the sum of all deposition eigenvalues gives the total energy within
the target region excited by all input degrees of freedom. Moreover, the
degree of control of the energy delivery via the deposition eigenchannels is
determined by the PDF of deposition eigenvalue $P(\zeta)$. In diffusive
systems, $P(\zeta)$ is amenable to a theoretical description within framework
of filtered random matrix (FRM) theory[15, 11]. The maximal and minimal
energies over the target are set by the range of $P(\zeta(\vec{r}_{D}))$.

## Interference of transmission eigenchannels

At first glance, the intensity sum rule in Eq. (3) might be misconstrued to
mean that the intensity profile of a channel in one basis, e.g. the deposition
eigenchannel, is a linear superposition of the intensity pattern in another
basis, e.g. of the transmission eigenchannel. Such a conclusion is markedly
wrong, as seen by comparing $\left|E^{(0)}_{n}(\vec{r})\right|^{2}$ and
$\left|E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})\right|^{2}$, and recalling that
the former set is composed of the statistically identical intensity profiles.
The fault of the above argument lies in neglecting the interference between
channels inside the scattering system. Indeed, we can express

$E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})=\sum_{\beta=1}^{N}E_{\beta}^{(T)}(\vec{r})d_{\beta\alpha}(\vec{r}_{D}),$
(5)

where the decomposition coefficient
$d_{\beta\alpha}(\vec{r}_{D})=\sum_{n^{\prime}=1}^{N}\left[V_{\beta
n^{\prime}}^{(T)}\right]^{*}V_{n^{\prime}\alpha}^{(D)}(\vec{r}_{D})$
represents a projection of input vector of the $\alpha$-th deposition
eigenchannel onto that of the $\beta$-th transmission eigenchannel. Using
these decomposition coefficients, the intensity pattern of a deposition
eigenchannel can be expressed in terms of two distinct terms:

$\begin{split}\left|E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})\right|^{2}=&\sum_{\beta=1}^{N}\left|E_{\beta}^{(T)}(\vec{r})\right|^{2}|d_{\beta\alpha}|^{2}\\\
+&\sum_{\beta\neq\beta^{\prime}}^{N}d_{\beta\alpha}\,d^{*}_{\alpha\beta^{\prime}}E_{\beta}^{(T)}(\vec{r})\left[E_{\beta^{\prime}}^{(T)}(\vec{r})\right]^{*}.\end{split}$
(6)

The first term is an incoherent sum of TE intensity patterns, whereas the
second term is the result of interference between different transmission
eigenchannels inside the scattering medium. The numerical simulation results
in Fig. 3a,b illustrate that the interference contributions can be positive or
negative in order to enhance or suppress the energy in the target. Their
existence is essential to make the energy higher than the largest-transmission
eigenchannel or lower than the smallest-transmission eigenchannel.

Figure 4: Coherent contribution in Eq. (6) for each deposition eigenchannel
$\alpha$, computed at $z=z_{D}$ in numerical model [12], is shown in the main
plot. A small number of large positive contributions is balanced by a large
number of smaller negative contributions in accordance with Eq. (8). The inset
shows that, when negative, the coherent contribution (red) cannot exceed the
incoherent one (blue) in absolute value, see Eq. (9). In both plots symbols /
lines represent one disorder realization / statistically averaged results
respectively.

The quantity $|d_{\beta\alpha}|^{2}$ represents the incoherent (i.e. a real
positive intensity) contribution due to $\beta$-th transmission eigenchannel
to the $\alpha$-th deposition eigenchannel. We obtain the sum rule [12]

$\sum_{\beta=1}^{N}|d_{\beta\alpha}|^{2}=1=\sum_{\alpha=1}^{N}|d_{\beta\alpha}|^{2}.$
(7)

This is a non-trivial result, since $|d_{\beta\alpha}|^{2}$ cannot be
interpreted as a weight coefficient due to presence of the interference term
in Eq. (6). Furthermore, it leads to an important constraint on this very
term. By summing both sides of that equation and using Eqs. (3,7), we get

$\sum_{\alpha=1}^{N}\left[\sum_{\beta\neq\beta^{\prime}}d_{\beta\alpha}\,d^{*}_{\alpha\beta^{\prime}}E_{\beta}^{(T)}(\vec{r})\left[E_{\beta^{\prime}}^{(T)}(\vec{r})\right]^{*}\right]\equiv
0.$ (8)

As shown in Ref. [11], the interference contribution can be quite large, even
dominant in some cases. The relationship in Eq. (8) states that the sum of
interference contributions to all deposition eigenchannels is, in fact, zero.
To illustrate this point, we plot [12] in Fig. 4 the coherent contributions at
the target depth $z_{D}=L/2$ for all deposition eigenchannels. Each
contribution is normalized by
$(1/N)\sum_{n=1}^{N}\left\langle\left|E^{(0)}_{n}(z_{D})\right|^{2}\right\rangle$,
which represents the unoptimized intensity at the target depth. A small number
of deposition eigenchannels have large positive contributions. In contrast,
the number of small negative contributions is large, to ensure the sum is
equal to 0 in Eq. (8).

Again we stress that the above relationships apply for every disorder
configuration and do not require any statistical averaging. Furthermore,
because the left hand side of Eq. (6) is a positively defined quantity, we
note that

$\sum_{\beta=1}^{N}\left|E_{\beta}^{(T)}(\vec{r})\right|^{2}|d_{\beta\alpha}|^{2}\geq-\sum_{\beta\neq\beta^{\prime}}^{N}d_{\beta\alpha}\,d^{*}_{\alpha\beta^{\prime}}E_{\beta}^{(T)}(\vec{r})\left[E_{\beta^{\prime}}^{(T)}(\vec{r})\right]^{*}.$
(9)

It illustrates that when the interference contribution does become negative,
i.e. for a low-deposition eigenhcannel, it cannot exceed in absolute value the
incoherent contribution, i.e. it cannot become dominant, see the inset in Fig.
4. However, such restriction does not apply for the high-deposition
eigenchannel with positive coherent contribution, which can and, in fact, does
become dominant in a diffusive medium for the deposition depth $z_{D}<L/2$,
c.f. Fig. 4 and Ref. [11]. We note that such larger positive interference
contributions are related to the PDF of the deposition eigenvalues. Since the
PDF $P(\zeta)$ has a long tail at large $\zeta$ in diffusive systems, the
number of positive contributions is small. Consequently, the sum rule in Eq.
(8) dictates that such contributions must be large in order to balance the
numerous negative ones.

## Conclusions

Targeted delivery of electromagnetic energy inside a random-scattering system
has important applications in imaging, optogenetics, photothermal therapy,
etc. The deposition eigenchannels accomplish the goal of delivering maximal or
minimal amount of energy to a target region of arbitrary size, shape and
depth. Little is known about the spatial structure of the deposition
eignchannels, albeit some progress has been made in understanding the spatial
distribution of transmission eigenchannels [14, 16, 17]. On the other hand,
the PDF of deposition eigenvalues has been predicted by the filtered random
matrix theory [11]. The sum rules, presented in this work, establish a
connection between the spatial structure of the deposition eigenchannels and
the eigenvalues. They represent the rigorous constraints for any disorder
realization of 2D and 3D scattering systems in both waveguide and slab
geometry.

## Acknowledgements

## Funding

We acknowledge Arthur Goetschy, Chia-Wei Hsu, Hasan Yılmaz for stimulating
discussions. This work was supported by the National Science Foundation under
Grant Nos. DMR-1905442, DMR-1905465.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

All data underlying the results is presented in this paper.

See Supplement 1 for supporting content.

## References

  * [1] P. A. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Nat. Photon. 6, 283 (2012).
  * [2] H. Yu, J. Park, K. Lee, J. Yoon, K. Kim, S. Lee, and Y. K. Park, Curr. Appl. Phys. 15, 632 (2015).
  * [3] J. Kubby, S.Gigan, and M. Cui, _Wavefront shaping for biomedical imaging_ (Cambridge University Press, 2019).
  * [4] S. Gigan, O. Katz, H. B. de Aguiar, E. R. Andresen, A. Aubry, J. Bertolotti, E. Bossy, D. Bouchet, J. Brake, S. Brasselet, Y. Bromberg, H. Cao, T. Chaigne, Z. Cheng, W. Choi, T. Čižmár, M. Cui, V. R. Curtis, H. Defienne, M. Hofer, R. Horisaki, R. Horstmeyer, N. Ji, A. K. LaViolette, J. Mertz, C. Moser, A. P. Mosk, N. C. Pégard, R. Piestun, S. Popoff, D. B. Phillips, D. Psaltis, B. Rahmani, H. Rigneault, S. Rotter, L. Tian, I. M. Vellekoop, L. Waller, L. Wang, T. Weber, S. Xiao, C. Xu, A. Yamilov, C. Yang, and H. Yılmaz, arXiv 2111.14908 (2021).
  * [5] S. Rotter and S. Gigan, Rev. Mod. Phys. 89, 015005 (2017).
  * [6] I. M. Vellekoop, Opt. Express 23, 12189 (2015).
  * [7] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Phys. Rev. Lett. 104, 100601 (2010).
  * [8] R. Horstmeyer, H. Ruan, and C. Yang, Nat. Photon. 9, 563 (2015).
  * [9] S. Yoon, M. Kim, M. Jang, Y. Choi, W. Choi, S. Kang, and W. Choi, Nature Reviews Physics 2, 141 (2020).
  * [10] X. Cheng and A. Z. Genack, Opt. Lett. 39, 6324 (2014).
  * [11] N. Bender, A. Yamilov, A. Goetschy, H. Yilmaz, C. W. Hsu, and H. Cao, Nat. Phys. 18, 309 (2022).
  * [12] Supplementary Material.
  * [13] O. N. Dorokhov, Solid State Comm. 51, 381 (1984).
  * [14] M. Davy, Z. Shi, J. Park, C. Tian, and A. Z. Genack, Nat. Comm. 6, 6893 (2015).
  * [15] A. Goetschy and A. D. Stone, Phys. Rev. Lett. 111, 063901 (2013).
  * [16] O. S. Ojambati, H. Yilmaz, A. Lagendijk, A. P. Mosk, and W. L. Vos, New J. Phys. 18, 043032 (2016).
  * [17] M. Koirala, R. Sarma, H. Cao, and A. Yamilov, Phys. Rev. B 96, 054209 (2017).
  * [18] R. Sarma, A. Yamilov, S. Petrenko, Y. Bromberg, and H. Cao, Phys. Rev. Lett. 117, 086803 (2016).
  * [19] N. Bender, A. Yamilov, H. Yilmaz, and H. Cao, Phys. Rev. Lett. 125, 165901 (2020).

## Supplement 1: Sum rules for energy deposition eigenchannels in scattering
systems

Alexey Yamilov1, Nicholas Bender2, Hui Cao3
1Physics Department, Missouri University of Science & Technology, Rolla,
Missouri, USA
2School of Applied and Engineering Physics, Cornell University, Ithaca, New
York, USA
3Department of Applied Physics, Yale University, New Haven, Connecticut, USA

### Derivation of Eq. (3)

Taking advantage of the unitarity of $V^{(D)}$ matrices,
$\sum_{n=1}^{N}\left[V_{\beta
n}^{(D)}(\vec{r}_{D})\right]^{*}V_{n\alpha}^{(D)}(\vec{r}_{D})=\delta_{\alpha\beta}$,
where $\delta_{\alpha\beta}$ is the Kronecker delta, we obtain

$\begin{split}&\sum_{\alpha=1}^{N}\left|E_{\alpha}^{(D)}(\vec{r};\vec{r}_{D})\right|^{2}\\\
&=\sum_{\alpha,nn^{\prime}=1}^{N}E^{(0)}_{n}(\vec{r})\,V_{n\alpha}^{(D)}(\vec{r}_{D})\,\left[V_{\alpha
n^{\prime}}^{(D)}(\vec{r}_{D})\right]^{*}\,\left[E^{(0)}_{n^{\prime}}(\vec{r})\right]^{*}\\\
&=\sum_{nn^{\prime}=1}^{N}E^{(0)}_{n}(\vec{r})\,\delta_{nn^{\prime}}\,\left[E^{(0)}_{n^{\prime}}(\vec{r})\right]^{*}=\sum_{n=1}^{N}\left|E^{(0)}_{n}(\vec{r})\right|^{2}.\end{split}$
(S1)

A similar relationship between intensities of TE
$\left|E_{\alpha}^{(T)}(\vec{r})\right|^{2}$ and
$\left|E^{(0)}_{n}(\vec{r})\right|^{2}$ can also be obtained analogously
completing the derivation of Eq. (3).

### Derivation of Eq. (4)

Singular values $\zeta_{\alpha}(\vec{r}_{D})$ can be extracted from
$\mathcal{Z}^{\dagger}(\vec{r}_{D})\mathcal{Z}(\vec{r}_{D})$ matrix. Using the
unitarity of $U^{(D)}(\vec{r}_{D})$, $V^{(D)}(\vec{r}_{D})$ matrices and the
definition of $\mathcal{Z}(\vec{r}_{D})$ in Eq. (1), we can find an explicit
relationship for the eigenvalue

$\begin{split}&\zeta_{\alpha}(\vec{r}_{D})=\sum_{nn^{\prime}=1}^{N}\sum_{m=1}^{M}\left[V_{\alpha
n^{\prime}}^{(D)}(\vec{r}_{D})\right]^{*}\mathcal{Z}_{n^{\prime}m}^{*}(\vec{r}_{D})\mathcal{Z}_{mn}(\vec{r}_{D})V_{n\alpha}^{(D)}(\vec{r}_{D})\\\
&=\frac{\mathcal{V}}{M}\,\sum_{nn^{\prime}=1}^{N}\sum_{m=1}^{M}\epsilon(\vec{r}_{m})\left[V_{\alpha
n^{\prime}}^{(D)}(\vec{r}_{D})\right]^{*}\left|E_{n}(\vec{r}_{m};\vec{r}_{D})\right|^{2}V_{n\alpha}^{(D)}(\vec{r}_{D})\\\
&=\frac{\mathcal{V}}{M}\,\sum_{m=1}^{M}\epsilon(\vec{r}_{m})\left|E_{\alpha}^{(D)}(\vec{r}_{m};\vec{r}_{D})\right|^{2}.\end{split}$
(S2)

We stress that this relationship is exact and holds in every disorder
realization.

### Derivation of Eq. (7)

Recalling the definition
$d_{\beta\alpha}(\vec{r}_{D})=\sum_{n^{\prime}=1}^{N}\left[V_{\beta
n^{\prime}}^{(T)}\right]^{*}V_{n^{\prime}\alpha}^{(D)}(\vec{r}_{D})$, the
quantity $|d_{\beta\alpha}|^{2}$ can be evaluated as

$\begin{split}\sum_{\beta=1}^{N}|d_{\beta\alpha}|^{2}&\equiv\sum_{nn^{\prime}\beta=1}^{N}V_{n\beta}^{(T)}\left[V_{\beta
n^{\prime}}^{(T)}\right]^{*}V_{n^{\prime}\alpha}^{(D)}(\vec{r}_{D})V_{\alpha
n}^{(D)*}(\vec{r}_{D})\\\
&\equiv\sum_{nn^{\prime}=1}^{N}\delta_{nn^{\prime}}V_{n^{\prime}\alpha}^{(D)}(\vec{r}_{D})\left[V_{\alpha
n}^{(D)}(\vec{r}_{D})\right]\\\
&\equiv\sum_{n=1}^{N}V_{n\alpha}^{(D)}(\vec{r}_{D})\left[V_{\alpha
n}^{(D)}(\vec{r}_{D})\right]^{*}\equiv 1,\end{split}$ (S3)

where we used the unitarity of the $V$ matrices. Summation over the other
index $\alpha$ yieslds the same result.

### Numerical simulation in the waveguide geometry

The results presented in Figs. 2,3,4 are obtained in two-dimensional waveguide
geometry, see Refs. [11, 17, 18, 19] for the detailed description of the
numerical model. The parameters of in the simulation are chosen to represent a
diffusive transport with no loss. Specifically, system length is $L/\ell\simeq
15$, where $\ell$ is the transport mean free path. The width of the waveguide
$W/L=0.3$, the number of modes in the waveguide $N=56$. The value of the
dimensionless conductance $g\simeq 5.3$.

Transmission eigenchannels are obtained using singular value decomposition of
the transmission matrix $\mathcal{T}$. The deposition matrix
$\mathcal{Z}_{mn}(\vec{r}_{D})$ is defined using Eq. (1) by sampling fields at
all spatial points $E^{(0)}_{n}(y_{m},z_{D})$ at a specific depth $z_{D}$. $y$
and $z$ are the transverse and the longitudinal coordinates respectively.
Subsequently, the deposition eigenchannels are obtained using singular value
decomposition of the rectangular matrix $\mathcal{Z}_{mn}(z_{D})$. The results
presented in Figs. 3,4 are obtained for the target region at the middle of the
disordered region $z_{D}=L/2$.

All numerical results in Figs. 2-4 are statistically averaged over an ensemble
of $1000$ random disorder realizations.