# Evidence of Gas Phase Nucleation of Nano Diamond in Microwave Plasma
Assisted Chemical Vapor Deposition

Tanvi Nikhar<EMAIL_ADDRESS>Department of Electrical and Computer
Engineering, Michigan State University, 428 S. Shaw Ln., East Lansing, MI
48824, USA  Sergey V. Baryshev<EMAIL_ADDRESS>Department of Electrical and
Computer Engineering, Michigan State University, 428 S. Shaw Ln., East
Lansing, MI 48824, USA  Department of Chemical Engineering and Materials
Science, Michigan State University, 428 S. Shaw Ln., East Lansing, MI 48824,
USA

###### Abstract

The mechanism of ballas like nano crystalline diamond formation (NCD) still
remains elusive, and this work attempts to analyze its formation in the
framework of activation energy ($E_{\text{a}}$) of NCD films grown from
H${}_{\text{2}}$/CH${}_{\text{4}}$ plasma in a 2.45 GHz chemical vapor
deposition system. The $E_{\text{a}}$ is calculated using the Arrhenius
equation corresponding to the thickness growth rate while using substrate
temperature ($\sim 1000-1300$ K) in all the calculations. While the calculated
values match with the $E_{\text{a}}$ for nano diamond formation throughout the
literature, these values of $\sim$10 kcal/mol are lower compared to $\sim$15 –
25 kcal/mol for standard single crystal diamond (SCD) formation, concluding
thus far, that the energetics and processes involved are different. In this
work, To further investigate this, the substrate preparation and sample
collection method are modified while keeping the growth parameters constant.
Unseeded Si substrates are physically separated from the plasma discharge by a
molybdenum plate with a pinhole drilled in it. Small quantities of a sample
substance are collected on the substrates. The sample is characterized by
electron microscopy and Raman spectroscopy confirming it to be nano diamond,
thus, suggesting that nano diamond self nucleates in the plasma and flows to
the substrate which acts as a mere collection plate. It is hypothesized then,
if nano diamond nucleates in gas phase, gas temperature has to be used in the
Arrhenius analysis. The $E_{\text{a}}$ values for all the nano diamond films
are re-calculated using the simulated gas temperature ($\sim 1500-2000$ K)
obtained from a simple H${}_{\text{2}}$/CH${}_{\text{4}}$ plasma model, giving
values within the range characteristic to SCD formation. A unified growth
mechanism for NCD and SCD is proposed concluding that the limiting reactions
for NCD and SCD formation are the same.

##  I Introduction

The growth of diamond through chemical vapor deposition (CVD) is reported
extensively on an existing diamond structure, such as on a diamond substrate
or on a diamond seed (achieved by various seeding processes performed on a
non-diamond substrate), and a number of studies have been done hypothesizing
the mechanism of adding C atoms from source gas to the base diamond structure
that results in growth Derjaguin and Fedoseev (1975); Harris and Weiner
(1990); Harris (1990); Frenklach (1991); Butler  _et al._ (1993); Wolden  _et
al._ (1994). The precipitation of solid diamond from gaseous carbonaceous
species without having an already defined diamond structure to grow on has
also been reported Spitsyn  _et al._ (1981); Maeda  _et al._ (1992); Pehrsson
_et al._ (1992); Kobayashi  _et al._ (1993); Chai  _et al._ (1994) but
analysis of the mechanism behind it is far less analyzed. To initiate growth
without an existing diamond structure to replicate on, the reactive species
from gas, at some point, need to form the first diamond structure – a nucleus,
on which the replication (epitaxy) process may follow to form a diamond
crystal.

Without an existing diamond seed, this nucleation process is mainly understood
to be originating at favorable nucleation sites on the surface of a substrate,
such as on defects introduced by scratching the substrate with abrasives
(including non-diamond materials that do not leave a diamond seed behind),
carbide forming substrates, etc. Liu and Dandy (1996a) For mirror-polished
surfaces, a commonly used technique is bias-enhanced nucleation, where charged
ions impinge upon the surface implanting C atoms in the upper layers of the
substrate and these act as active sites for nucleation Yugo  _et al._ (1991);
Garcia Poza  _et al._ (2000); Lifshitz  _et al._ (2002). Hence, diamond
nucleation is widely understood to be a surface phenomenon.

In low-pressure CVD systems though, gas phase nucleation of diamond has been
rarely discussed, but there are some reports, however scarce, discussing its
possibility based on experimental observations. Mitura $et$ $al.$ Mitura
(1987) observed agglomerates of diamond particles separate from the rest of
the cohesive film deposited in an RF methane plasma and concluded that these
particles must have formed in the gas phase, separate from the film formed on
the substrate surface. Howard $et$ $al.$ Howard  _et al._ (1990) reported the
synthesis of diamond powder in the gas phase by microwave-assisted combustion
of acetylene in oxygen at low pressures. Frenklach $et$ $al.$ Frenklach  _et
al._ (1989) were successfully able to obtain diamond downstream in a low-
pressure CVD reactor with mixtures of dichloromethane and oxygen hence
definitely confirming that nucleation does not necessarily require a surface
and is happening in the gas phase, although tedious post-treatment of the
collected sample was required to reveal the presence of diamond grains. They
were also able to achieve induced nucleation of diamond by adding diborane to
a mixture of acetylene, hydrogen and argon.Frenklach  _et al._ (1991)

Various theoretical studies as reviewed by Liu $et$ $al.$ Liu and Dandy
(1996b), have reported on the plausibility of gas phase nucleation of diamond.
However, the high nucleation density necessary for the formation of a
continuous film can only be attributed to surface nucleations. Hence, the
importance of gas phase nucleation is not well understood and is recognized at
best as an accompanying process to processes taking place on the substrate
surface. Hwang $et$ $al.$ Hwang  _et al._ (1996) presented a theoretical
confirmation of gas phase nucleation based on the study of carbon atom
potential. But there is not enough experimental proof seen in a standard low-
pressure hydrogen methane plasma CVD.

The hypotheses on the diamond formation mechanism via CVD, include a number of
proposed reactions and models of the growth process to investigate the
energies involved at various steps of these reactions.Wang and Frenklach
(1991); Harris  _et al._ (1991); Frenklach (1992); Garrison  _et al._ (1992);
Frenklach (1994); Kang and Musgrave (2000); Cheesman  _et al._ (2008) A common
thread in all these studies is the treatment of the growth process as a
surface phenomenon. Gas kinetics have mostly been studied separately from
surface kinetics because combining the two would immensely increase the
complexity of the system to be modeled.Frenklach and Wang (1991) These
calculations are done based on the assumption that surface kinetics leading to
deposition can be considered independent of gas kinetics, which are then used
to thermodynamically warrant the proposed mechanisms. This leaves out the
possibility of important processes occurring in the gas phase that may end up
altering the overall kinetics of the process. It begs the question of whether
the growth kinetics can be treated as an independent function of surface
phenomena.

The commonly accepted mechanism for the growth of diamond ($sp^{3}$ phase of
carbon) is by CH${}_{\text{3}}$ radical incorporation. The activation energy
calculated through experimentally obtained data for diamond films with varying
grain sizes reveals that for larger grains having faceted morphology as seen
in microcrystalline diamond films (MCD), the activation energy is $\sim$10
kcal/mol. Barbosa  _et al._ (2014) This energy matches up with the energy
required for CH${}_{\text{3}}$ incorporation into the lattice, thus, making it
the rate-limiting step for such type of diamond growth. On the other hand, for
smaller crystallites of diamond having ballaslike morphology as seen in NCD
films, the activation energy was found to be consistently lower, $\sim$6
kcal/mol. Barbosa  _et al._ (2014) Thus, the formation of ballaslike diamond
morphology was hypothesized to be different from the canonical
CH${}_{\text{3}}$ incorporation mechanism. The reason for the lowering of the
activation energy was not clear and was attributed to various possibilities:
competition of diamond grain with non-diamond phase growth (especially
graphitic); a mechanism led by C${}_{\text{2}}$ dimer incorporation into the
lattice which corresponds to $\sim$6 kcal/mol; high renucleation rate that
probably interrupts $sp^{3}$ grain growth by imposing another rate-limiting
step that is also unknown.

What is concluded so far is that faceted and ballaslike morphologies must
follow different reaction mechanisms owing to their different activation
energies of formation. In this work, we investigate further on this difference
to get a better understanding of the mechanism for the formation of ballaslike
morphology. The activation energies of ballaslike NCD films produced in our
previous work Nikhar  _et al._ (2020) are calculated from the Arrhenius
equation used to approximate the relation of growth rate as a function of
temperature. New syntheses is carried out by depositing ballaslike NCD using
two distinct methods; in one of which plasma was decoupled from the substrate.
From this, a hypothesis of gas phase nucleation of diamond from hydrogen
methane plasma in a conventional low-pressure CVD system is formulated. Aided
by a basic computational model describing the
H${}_{\text{2}}$/CH${}_{\text{4}}$ microwave assisted CVD plasma, a correction
in the calculation of activation energy for the grown films is proposed. The
recalculated energy values are then used to discuss inferences about the
growth mechanism involved in ballaslike morphology formation.

##  II Samples

The E${}_{\text{a}}$ calculations in this study are done for the NCD films
deposited on intrinsic Si(100) substrates in a 2.45 GHz microwave plasma
assisted chemical vapor deposition (MPACVD) system.Asmussen and Reinhard
(2002) The films were comprehensively described and characterized in our
previous work. Nikhar  _et al._ (2020) We particularly focus on the films
synthesized in pure H${}_{\text{2}}$/CH${}_{\text{4}}$ plasma denoted by the
SA series corresponding to 0% nitrogen summarized in Table I of Ref. Nikhar
_et al._ , 2020.

All samples were characterized with a 532 nm probing laser on a Horiba Raman
spectrometer. High-resolution scanning electron microscopy (SEM) was performed
using a JEOL JSM 7500F to study the surface morphology of the samples. The Si
substrates are weighed before and after deposition using a Mettler Toledo XS
105 balance with 10${}^{\text{-5}}$ g readability to find the weight of the
deposited material given by $\Delta w$.

##  III Activation Energy Results

###  III.1 Thickness calculation

The Raman spectra of the NCD films presented in our previous workNikhar  _et
al._ (2020) reveal that as the substrate temperature, $T_{s}$ increases, the D
peak shifts from 1333 cm${}^{\text{-1}}$ to 1350 cm${}^{\text{-1}}$ and the G
peak gets narrower and shifts towards 1590 cm${}^{\text{-1}}$. The evident
phase transformation of the material from ultra-nano-crystalline diamond to
nano-crystalline graphite is used here to approximate the thickness of the
deposited material using the following relation:

$t=\frac{\Delta w\cdot 10^{4}}{S_{a}\cdot\rho}~{}\text{in}~{}[\mu\text{m}]$
(1)

where, $\Delta w$ = weight difference before and after deposition (g), $S_{a}$
= surface area of the film (cm${}^{\text{2}}$), and
$\rho$ = $\begin{cases}\rho_{dia},&\text{if the D peak is at 1333
cm${}^{\text{-1}}$}\\\ \rho_{gra},&\text{if the D peak is shifted to 1350
cm${}^{\text{-1}}$}\end{cases}$ where, $\rho_{dia}$ = 3.51 g/cm3 is the
density of pure diamond, and $\rho_{gra}$ = 2.27 g/cm3 is the density of pure
graphite.

The growth of the films is interpreted in terms of thickness $t$. Since all
the films are deposited over a period of 1 hour, the growth rate is given by
$k$ = thickness/time = $t[\mu\text{m}]$/hour.

###  III.2 Activation energy calculation

The activation energy for the growth is calculated using the Arrhenius
equation, given as:

$k=A\cdot e^{-E_{a}/RT}$ (2)

Hence,

$E_{a}=-RT\hskip 5.0ptln(k/A)$ (3)

where, $R=8.31$ J/mol$\cdot$K = $8.31\times 0.000239$ kcal/mol$\cdot$K is the
ideal gas constant, $T$ = temperature (in Kelvin), $k$ = rate constant, and
$A$ = pre-exponential factor calculated from the Arrhenius plot: ln($k$) vs
1/$T$ whose intercept is ln($A$).

The term $k$ in the expression of $E_{a}$ is the rate constant of the reaction
taking place at temperature $T$. Here, this temperature corresponds to the
reactions that lead to the deposition of the film. A common assumption
throughout literature so far has been that since the deposition of the
nanodiamond particles happens on the substrate surface which eventually builds
up to form the NCD film, the reactions involved in the formation of these
particles must happen on the surface itself. Thus, in all the studies, the
temperature of the substrate measured experimentally is used to calculate the
activation energy for the formation of the films.

For our experiments, the temperature of the substrate $T_{s}$ was measured
using an infrared pyrometer, the values of which are reported as deposition
temperature in Table I in Ref. Nikhar  _et al._ , 2020. The activation
energies calculated using Eq. 3 with $T=T_{s}$ are summarized in Fig. 1 above
the label $T_{s}$ at the 1,000 K mark. Activation energies for the formation
of pure graphiteFedoseev  _et al._ (1979) and SCDKondoh  _et al._ (1993); Kang
and Musgrave (2000) are also plotted for comparison. This result agrees
extremely well with the results of the previous meta-studyBarbosa  _et al._
(2014) that summarized and analyzed data on nanodiamond obtained by different
deposition techniques. Since $E_{a}$ for nanodiamond corresponding to $T_{s}$
has consistently been found lesser than that for SCD, the mechanisms for the
formation of NCD and SCD have been assumed to be different. As summarized by
Barbosa $et$ $al.$, Barbosa  _et al._ (2014) the rate-limiting reaction for
SCD formation is understood to be CH${}_{\text{3}}^{\bullet}$ incorporation
while the limiting reaction for NCD is not known yet.

Figure 1: Activation energies of NCD films calculated using substrate
temperature ($T_{s}$) and gas temperature ($T_{g}$) as compared with the
activation energy of formation of SCD and pure graphite.

The formation of nanodiamond films is attributed to a high renucleation rate.
Based on their observation of an orange glow on the edges of the plasma
discharge, May $et$ $al.$May  _et al._ (2010) suggested that the high
renucleation rate could be due to particulates present in the outer edge of
the plasma. These particulates are speculated to act as nucleation or defect
sites on the growth surface promoting film growth and a high renucleation
rate. A similar orange glow was also observed in the discharges produced in
our previous work on nanodiamond synthesis.Nikhar  _et al._ (2020)

NCD synthesis on unseeded Si substrate was attempted to find whether self-
seeding (from the afore-mentioned particulats) and further coalescence and
continuous film growth was attainable. The results as presented in Ref. Nikhar
_et al._ , 2020 showed that the first hour of deposition resulted in the
formation of colonies mostly separate from each other. After the second hour,
some of the colonies were found to merge together towards the formation of a
continuous film. SEM images and Raman spectra of these colonies confirmed that
the samples obtained on seeded and unseeded substrates were identical. The
synthesis runs were repeated on seeded and unseeded Mo substrates, each
producing identical samples to the ones on Si substrates.

These results provided evidence of self-seeding, independent of the type of
substrate or the substrate pre-treatment process. These results corroborated
the previous hypotheses of diamond gas phase nucleation, where the seeding
diamond particulates are formed in the plasma itself independent of the
substrate. If true, then the temperature term in Eqs. 2 and 3 should be taken
as the gas temperature $T_{g}$ and not the substrate temperature $T_{s}$.

##  IV Experiment

To confirm this hypothesis, in the present work, we further modified our
original setup (Fig. 2 (left)) to physically separate the plasma from the
substrate by placing a Mo disc on top of the substrate removing any direct
contact with the plasma as shown in Fig. 2 (right). An unseeded Si substrate
of $\sim$1 cm $\times$ 1 cm$~{}\times$ 0.3 mm was placed in a pocket of the Mo
substrate holder on top of which a Mo disc of diameter 7 cm and thickness 3.31
mm with a pinhole was placed with the intention of allowing any particles
forming in the plasma to pass through the pinhole onto the substrate while
keeping it physically separate from the plasma.

The diameter of the disc is chosen wide enough so that the entire plasma is
contained on top of it and does not flow off its edge, thus keeping the shape
of the plasma unperturbed. The diameter of the pinhole is chosen in accordance
with the debye length of the plasma near the surface of the disc, ensuring
that the plasma do not reach down to the substrate, hence maintaining physical
separation.

Figure 2: Schematics showing the two setup arrangements used for the
deposition of NCD samples: (left) Original setup: the samples are collected on
a substrate directly exposed to the plasma; (right) Altered setup: the samples
are collected through a pinhole in a Mo stub used to physically separate the
plasma from the substrate.

Figure 3: SEM images of sample collected (a), (b) in the original setup as
compared to samples collected in the altered setup through the pinhole (c),
(d).

Figure 4: Raman spectra comparing samples collected by the original and
altered methods, both showing D and G bands typical to nanodiamond.

Samples are collected with both setups using process conditions: a hydrogen-
rich H2/CH4 mixture was used with the CH4 flow rate was kept at 10 sccm (5% by
volume), while the flow rate of H2 was maintained at 200 sccm.

With the original setup leading to in-plasma collection, using seeded Si as
the substrate, a continuous film of thickness $\sim$500 nm was formed over a
deposition duration of 1 hour. The SEM images and Raman spectrum shown in Fig.
4(a),(b) and 4 respectively, present characteristics of the film typical for
nanodiamond. With the new setup using unseeded Si as the substrate, the out-
of-plasma sample collection through the pinhole was only observed after a
deposition duration of 4 hours, exclusively over the substrate area exposed to
the pinhole opening. A very low yield of the sample is obtained having poor
adherence to the substrate and can be easily removed by a cotton swab. (The
films obtained through conventional in-plasma deposition adhere strongly to
the substrate.) The SEM analysis of the sample collected through the pinhole
revealed the presence of individual colonies scattered only under the area of
the pinhole while the rest of the substrate remains clear of any sample. The
diameter of each colony ranges from 118 to 262 nm with the morphology as shown
in Fig. 4(c),(d) identical to the nanodiamond film shown in Fig. 4(a),(b). The
Raman spectrum shown in Fig. 4 in red is also similar to that of the original
setup shown in blue.

The obtained results even further support the gas phase nucleation of
nanodiamond formation. It indeed calls for replacing $T_{s}$ with $T_{g}$ when
calculating the activation energy in Eq. 3. To obtain $T_{g}$ parameter, a
model of the reactor was built to estimate the basic plasma parameters.

##  V Numerical model of the plasma

####  V.0.1 Model formulation

Numerical modeling of microwave CH${}_{\text{4}}$/H${}_{\text{2}}$ plasma was
carried out in COMSOL with consideration of 2D computational domain with an
axisymmetric problem statement. The corresponding symmetry axis crosses the
center of the quartz dome and the substrate holder at their centers.

CH${}_{\text{4}}$/H${}_{\text{2}}$ plasma chemistry modeling is based on the
solution of species and momentum conservation equations with drift-diffusion
approximation for electrons and positive ions. The electron conservation
equation is written as:

$\frac{\partial}{\partial t}(n_{e})+\nabla\bm{\cdot}\vec{\Gamma_{e}}=R_{e}$
(4)

where $n_{e}$ is the electron density ($m^{-3}$). $\vec{\Gamma_{e}}$ is the
electron flux vector derived from the momentum conservation equation and is
given by:

$\vec{\Gamma_{e}}=-n_{e}\mu_{e}\vec{E}-D_{e}\nabla n_{e}$ (5)

where $\mu_{e}$ and $D_{e}$ are the electron mobility ($m^{2}/Vs$) and
diffusion ($m^{2}/s$) constants respectively. $\vec{E}$ is the electric field.
Suppose that there are $M$ reactions that contribute to the growth/decay of
electron density, then the electron source term $R_{e}$ is given by:

$R_{e}=\sum_{j=1}^{M}x_{j}k_{j}N_{n}n_{e}$ (6)

where $x_{j}$ is the mole fraction of the target species for reaction $j$,
$k_{j}$ is the rate constant ($m^{3}/s$) for reaction $j$, and $N_{n}$ is the
total neutral density ($m^{-3}$).

$\rho_{i}\frac{\partial w_{i}}{\partial
t}=\nabla\bm{\cdot}\vec{\Gamma}_{i}+R_{i}$ (7)

The electron energy balance equation solved to calculate electron temperature
is given by:

$\frac{\partial}{\partial
t}(n_{\varepsilon})+\nabla\bm{\cdot}\vec{\Gamma}_{\varepsilon}+\vec{E}\bm{\cdot}\vec{\Gamma}_{e}=S_{en}-(\vec{u_{e}}\bm{\cdot}\nabla)n_{e}+Q/q$
(8)

where $n_{\varepsilon}$ is the energy density ($V/m^{3}$).
$\vec{\Gamma}_{\varepsilon}$ is the electron energy flux given by:

$\vec{\Gamma}_{\varepsilon}=-\mu_{\varepsilon}n_{\varepsilon}\vec{E}-D_{\varepsilon}\nabla
n_{\varepsilon}$ (9)

where $\mu_{\varepsilon}$ is the electron energy mobility ($m^{2}/Vs$),
$D_{\varepsilon}$ is the electron energy diffusivity ($m^{2}/s$). Suppose that
there are $P$ inelastic electron-neutral collisions, then the energy loss/gain
due to inelastic collions $S_{e}n$ ($V/m^{3}s$) is can be given by:

$S_{e}=\sum_{j=1}^{P}x_{j}k_{j}N_{n}n_{e}\Delta\varepsilon_{j}$ (10)

where $\Delta\varepsilon_{j}$ is the energy loss from reaction j ($V$).
$\vec{u_{e}}$ is the elecron drift velocity vector, and $Q$ is the heat source
($W/m^{3}$).

Ionization, electronic and vibrational state excitation of H${}_{\text{2}}$,
and dissociation cross-sections of H${}_{\text{2}}$ and CH${}_{\text{4}}$
Mankelevich  _et al._ (2008a); Hassouni  _et al._ (2010); Yoon  _et al._
(2008) were specified to get an accurate calculation of electron temperature.
Since CH${}_{\text{4}}$ presence in the mixture is 5%, it was assumed that
electron kinetics and microwave plasma properties are dominated by electron
collisions with H and H${}_{\text{2}}$.

Volume Reaction | Description
---|---
H${}_{\text{2}}$ \+ e → H${}_{\text{2}}$ \+ e | Elastic collisions
H + e → H + e | Elastic collisions
H${}_{\text{2}}$ \+ e → H${}_{\text{2}}$(V) + e | 1-3 vibrational levels excitation
H${}_{\text{2}}$ \+ e → H${}_{\text{2}}$ \+ e | Electronic levels excitation
H${}_{\text{2}}$ \+ e → H${}_{\text{2}}^{\text{+}}$ \+ 2e | Direct ionization
H${}_{\text{2}}$ \+ e → 2H + e | Dissociation
H + e → H${}^{\text{+}}$ \+ 2e | Direct ionization
H${}_{\text{2}}^{\text{+}}$ \+ H${}_{\text{2}}$ → H${}_{\text{3}}^{\text{+}}$ \+ H | H${}_{\text{3}}^{\text{+}}$ ions formation
H${}_{\text{3}}^{\text{+}}$ \+ H → H${}_{\text{2}}^{\text{+}}$ \+ H${}_{\text{2}}$ |
H${}_{\text{3}}^{\text{+}}$ \+ e → 3H |
H${}_{\text{3}}^{\text{+}}$ \+ e → H${}_{\text{2}}$ \+ H |
H${}_{\text{2}}^{\text{+}}$ \+ e → 2H |
2H${}_{\text{2}}$ → 2H+H${}_{\text{2}}$ |
H${}_{\text{2}}$+H → 3H |
2H+H${}_{\text{2}}$ → 2H${}_{\text{2}}$ |
3H → H${}_{\text{2}}$ \+ H |
CH${}_{\text{4}}$ +H → CH${}_{\text{3}}$ \+ H${}_{\text{2}}$ |
CH${}_{\text{3}}$ \+ H${}_{\text{2}}$ → CH${}_{\text{4}}$ \+ H |
CH${}_{\text{3}}$ \+ H → CH${}_{\text{2}}$ \+ H${}_{\text{2}}$ |
CH${}_{\text{2}}$ \+ H${}_{\text{2}}$ → CH${}_{\text{3}}$ \+ H |
CH${}_{\text{2}}$ \+ H → CH + H${}_{\text{2}}$ |
CH+H${}_{\text{2}}$ → CH${}_{\text{2}}$ \+ H |
CH + H → C + H${}_{\text{2}}$ |
C + H${}_{\text{2}}$ → CH${}_{\text{2}}$ \+ H |
CH${}_{\text{3}}$ \+ C → C${}_{\text{2}}$H${}_{\text{2}}$ \+ H |
C${}_{\text{2}}$H${}_{\text{2}}$ \+ H → CH${}_{\text{3}}$ \+ H${}_{\text{2}}$ |
Table 1: Summary of the gas phase reactions.

(a) CH${}_{\text{3}}$ density distribution

(b) Gas temperature distribution

(c) $T_{g}$ $\oplus$ CH${}_{\text{3}}^{\bullet}$ distribution

Figure 5: Contour plots representing the 2D spatial distribution of various
plasma parameters across the reactor.

A simplified model of neutral species interaction was developed to describe
the spatial evolution of Cx=1,2Hy=0-3 radicals within the chamber volume and
their fluxes to the surfaces.Dandy and Coltrin (1995); Gorbachev and Leshchev
(2022) The whole set of gas phase reactions considered is presented in Table
1. Surface recombination of charged species was set to 1, and all carbon
containing radicals are lost on the quartz dome walls with unity probability.

Electric field distribution was calculated based on the Helmholtz equation
given as:

$\nabla\times(\nabla\times\vec{E})-k_{0}^{2}\left(\varepsilon_{r}-\frac{j\bm{\cdot}\sigma_{n}}{\omega\bm{\cdot}\varepsilon_{0}}\right)\vec{E}=0,$
(11)

where $k_{0}=\omega\sqrt{\varepsilon_{0}\mu_{0}}$ is the wave number in space.
$\omega$ is the frequency, $\varepsilon_{0}$ is the absolute permittivity and
$\mu_{0}$ is the absolute permeability of free space, $\varepsilon_{r}$ is the
relative permittivity of the medium. $\sigma_{n}$ is the conductivity of the
medium that gives rise to ohmic losses defined by the term
$(j\bm{\cdot}\sigma_{n})/(\omega\bm{\cdot}\varepsilon_{0}$). Electric
conductivity and dielectric permittivity provide coupling between plasma and
Helmholtz equation through the plasma frequency and electron transport
frequency.

A heat transfer equation was used to couple the process of gas heating and
temperature dependent chemical reaction with electromagnetic wave energy
absorption in plasma. The heat balance equation can be written as:

$\rho_{n}C_{p}(\vec{u}\bm{\cdot}\nabla T)+\nabla\bm{\cdot}\vec{q}=Q_{mw}$ (12)

where $\rho_{n}$ is the fluid density $(kg/m^{-3})$, $C_{p}$ is the fluid heat
capacity at a constant pressure $(J/kgK)$, $\vec{u}$ is the velocity vector
responsible for convection in fluid $(m/s)$. In a continuous medium,
conductive heat flux $\vec{q}$ is proportional to the temperature gradient,
hence $\vec{q}=-k\nabla T$, where $k$ is the thermal conductivity $(W/mK)$.
$Q_{mw}$ is the heat source from microwave energy.

Finally, gas inlet and outlet were considered in the model to take into
account the process of CH${}_{\text{4}}$ and H${}_{\text{2}}$ neutrals renewal
due to gas flow. Navie-Stocks equation given below was solved to obtain the
spatial evolution of gas neutrals within the chamber volume.

$\rho_{n}(\vec{u}\bm{\cdot}\nabla)\vec{u}=\nabla\bm{\cdot}(-pI+\tau)+\vec{F}$
(13)

where $p$ is the pressure $(Pa)$, $I$ is the identity tensor, $\tau$ is the
viscous stress tensor $(Pa)$, and $\vec{F}$ is the volume force vector
$(N/m^{3})$.

####  V.0.2 $T_{g}$ and methyl radical distributions

Fig. 5 summarizes the computational results of the plasma model for the input
reactor parameters detailed in Section IV: (a) and (b) are contour plots
representing the 2D spatial distribution of the methyl radical
CH${}_{\text{3}}^{\bullet}$ and gas temperature $T_{g}$, respectively. Noting
that the regions of maximum CH${}_{\text{3}}^{\bullet}$ concentration and
maximum $T_{g}$ do not overlap, Fig. 5(c) is the convolution of distributions
(a) and (b): $T_{g}$ $\oplus$ CH${}_{\text{3}}^{\bullet}$. This is done
because, a priori, the location of the reaction zone (region of the plasma
where diamond nucleation and subsequent grain formation occurs) is not known,
estimates of the reaction zone are made based on:
($i$) the region of maximum CH${}_{\text{3}}^{\bullet}$ concentration from
Fig. 5(a);
($ii$) the region of maximum gas temperature from Fig. 5(b);
($iii$) the region where the concentration of CH${}_{\text{3}}^{\bullet}$ and
gas temperature are both optimal favoring the desired reactions: we estimate
this region as the maximum of the convolution of CH${}_{\text{3}}^{\bullet}$
concentration and $T_{g}$ from Fig. 5(c).

The respective gas temperatures obtained: $T_{g}^{CH_{3}^{\bullet}}$,
$T_{g}^{T}$ and $T_{g}^{\oplus}$ corresponding to the maxima of
CH${}_{\text{3}}^{\bullet}$, $T_{g}$ and $T_{g}$ $\oplus$
CH${}_{\text{3}}^{\bullet}$, respectively, are then used to re-calculate the
activation energy from Eq. 3. The new values are plotted together with the
activation energies calculated using $T_{s}$ in Fig. 1. Each of the four
E${}_{\text{a}}$ values, (obtained using $T_{s}$ and three estimates for
$T_{g}$) is presented here with an error bar that originates from the linear
fitting of the Arrhenius plot ln$(k)$ vs $1/T$ to obtain the term ln$(A)$ in
Eq. 3.

##  VI Discussion

####  VI.0.1 Comparison with SCD formation

$E_{a}$ was reported to range from 13.6 to 22.5 kcal/mol for the formation of
SCDKondoh  _et al._ (1993); Kang and Musgrave (2000) and 0.98 kcal/mol for
graphite.Fedoseev  _et al._ (1979) The $E_{a}$ values for nanodiamond films
calculated using $T_{s}$ range from 7.87 to 11.29 kcal/mol, as shown in Fig.
1. This range matches very well with the values 6–10 kcal/mol reported in the
literature for nanodiamond formation,Barbosa  _et al._ (2014) where $E_{a}$
was always calculated using the substrate temperature. The revisited $E_{a}$
calculation gives a higher range of values: 13.5 $\pm$ 1.5 kcal/mol for
$T_{g}^{CH_{3}^{\bullet}}$, 17.5 $\pm$ 2.5 kcal/mol for $T_{g}^{T}$ and 15
$\pm$ 2 kcal/mol for $T_{g}^{\oplus}$. These values are well within the
characteristic range of SCD growth kinetics.

SCD ($sp^{\text{3}}$ carbon) formation is linked with the presence of
CH${}_{\text{3}}^{\bullet}$ in the plasma coupled with a high atomic H content
to keep CH${}_{\text{3}}^{\bullet}$ saturated with $\sigma$-bonded H
terminations important for subsequent pure sp${}^{\text{3}}$ phase formation
and growth. Butler $et$ $al.$ Mankelevich  _et al._ (2008b) present a
numerical simulation of the plasma used for SCD formation. Because of the high
$T_{g}$ in the plasma core ($\sim$2900 K) achieved by high operating pressures
of 150 Torr, molecular species are broken down to elementary parts and
CH${}_{\text{3}}^{\bullet}$ is distributed as a ring towards the outer cooler
edge of the plasma. The optimal combination is when
CH${}_{\text{3}}^{\bullet}$ concentration is on the order of
10${}^{\text{13}}$ cm${}^{\text{-3}}$ at $T_{g}\sim 1200-1500$ K. This
reaction zone is exactly where the SCD substrate is placed consuming the
incoming flux of CH${}_{\text{3}}^{\bullet}$. Following the hydrogen
abstraction reactions, C atoms are added to the SCD substrate homoepitaxially,
replicating the existing diamond structure.

When compared to the described optical SCD growth conditions, our numerical
simulation results illustrate that under 40 Torr operating condition the core
of the plasma is much cooler with $T_{g}\sim 2000$ K. As in the SCD plasma
case, the region for maximal CH${}_{\text{3}}^{\bullet}$ concentration does
not coincide with the maximal $T_{g}$ region. In fact, the maximal
CH${}_{\text{3}}^{\bullet}$ concentration is found to be at $T_{g}\sim 1400$
K, i.e. lying in the same temperature range as for the SCD plasma case.
Because of the lower core temperature, the high CH${}_{\text{3}}^{\bullet}$
concentration on the order of $~{}10^{14}$ cm-3 is comparatively more evenly
distributed throughout the plasma and the maxima of this distribution is not
just concentrated on the edge. This makes sp${}^{\text{3}}$ formation possible
throughout the plasma giving rise to nucleation of diamond in the gas phase.

Because of the low pressure though, the nucleated grain does not reach the
critical size after which it can keep growing as a stable crystal. Thus, the
process of renucleation begins giving rise to small crystal sized
polycrystalline diamond. If now pressure is increased without changing any
other parameters, the power density changes, the environment is no more in the
metastable regime of diamond growth and the samples thus obtained are
graphitized soot.Nikhar  _et al._ (2020) Instead, if the pressure is increased
and the power is reduced to balance the power density, it creates a suitable
plasma that can sustain the growth of a diamond grain as it grows larger and
hence the renucleation step is eliminated here.

####  VI.0.2 NCD film on SCD substrate

The reasoning of renucleation due to low pressure was also supported by
additional experimental observation. Synthesis parameters for NCD were applied
to deposit a sample on an optical grade SCD substrate. Even though the C atoms
from the plasma have a stable structure to replicate on and grow as a single
crystal, the resulting final product was a highly renucleated NCD film. The
resulting NCD on SCD film is shown in Fig. 6. The SEM micrographs illustrate
ballas-like structure with grain sizes comparable to the ones obtained on Si.
This result is further confirmed by Raman spectroscopy before (SCD substrate)
and after (NCD film) deposition. When the same reactor when operated at higher
pressures (120 to 400 Torrs) and maintained at a constant substrate
temperature of $\sim$1200 K (achieved by decreasing power and active substrate
cooling), where the SCD substrate is directly subject to the ring shaped flux
of high CH${}_{\text{3}}^{\bullet}$ concentration, resulted in high quality
single crystal homoepitaxy.Muehle  _et al._ (2017) These results highlight the
crucial importance of plasma kinetics on the final diamond product (NCD vs.
SCD).

Figure 6: (top) SEM images, and (bottom) Raman spectrum of NCD film deposited
on SCD substrate via in-plasma collection.

####  VI.0.3 Proposed Mechanism

##### Nanodiamond nucleation:

In the gas phase, the particles precipitate from the supersaturated vapor into
diamond crystals. In a H${}_{\text{2}}$/CH${}_{\text{4}}$ plasma, carbon can
condense into both $sp^{\text{2}}$ and $sp^{\text{3}}$ phases from the gas
phase. The preferential formation of diamond over graphite is usually
explained by etching of $sp^{\text{2}}$ phase by atomic H but Hwang $et$ $al.$
Hwang  _et al._ (1996) stated that H etching is not thermodynamically favored
and instead proposed the capillary effect as the major effect behind
$sp^{\text{3}}$ nucleation. The capillarity is expressed by Laplace-Young
equation as

$P_{2}-P_{1}=2\sigma/r,$ (14)

where $P_{2}$ and $P_{1}$ are the pressures inside and outside the nucleus,
respectively, of radius $r$ and surface energy $\sigma$. Given the mass
density of diamond is larger than that of graphite, for the same number of C
atoms precipitating into $sp^{\text{2}}$ and $sp^{\text{3}}$ phase forming a
nucleus, the radius of this nucleus corresponding to $sp^{\text{3}}$ phase
would be smaller than for $sp^{\text{2}}$ phase. Hence, for a given global gas
pressure in the plasma of a low-pressure CVD reactor, $P_{1}$, the local
pressure inside the nucleus, $P_{2}$ of $sp^{\text{3}}$ phase would be higher
than that of $sp^{\text{2}}$ phase making it comparable to the pressure for
making diamond more stable than graphite (local high pressure high
temperature, HPHT condition).

Based on this explanation, it can be inferred that as the diamond nucleates
and grows to a certain size, the local pressure being inversely proportional
to the radius is no longer high enough to stabilize this larger crystal. The
crystal then stops growing and a new nucleus starts to form, hence the high
renucleation rate in our low-pressure CVD environment. The gas phase growth of
the diamond nanoparticle and the subsequent renucleation depends on the gas
pressure. Lower the gas pressure, smaller the nanodiamond particles and higher
the renucleation rate. For our conditions, the $r\sim 10$ nm. As the operating
pressure drops, nanodiamond stops nucleating and growing due to the decreasing
gas temperature as $T_{g}$ is a strong function of the operating pressure.

##### Nanodiamond particle and film growth:

After subsequent renucleations, the condensed particles can no longer remain
suspended in the plasma and drop down onto the substrate. This is illustrated
by the schematic shown in Fig. 7. In a recent work by Xiong $et$ $al.$,Xiong
_et al._ (2022) it was noted that for nanometer-sized particles, gravity does
not play a significant role. There is trapping in the plasma occurring out of
the balance between the electrostatic force and opposing forces such as the
drag of the flowing neutral gas, the ion drag, and the thermophoretic force.
It was shown that sub-10 nm Si particles formed in a radiofrequency reactor
due to nuclei being temporarily confined in a trap formed by electrostatic
potential. Those particles only escaped after they grew to a size at which the
increasing drag force was imparted by the flowing gas carrying the particles
out. Because, generally, nanoparticles are expected to charge in plasmas a
similar mechanism could be at play in our system. This would explain particles
grow as they renucleate while being suspended in the gas phase until drag
forces can take over to expel the cluster out of the plasma.

Figure 7: Schematic of the proposed nucleation and growth of NCD film.

This is further hypothesized that the main expelling force is of
thermophoretic nature. There were two important references to note: Tsugawa
$et$ $al.$Tsugawa  _et al._ (2010) and Tallaire $et$ $al.$Tallaire  _et al._
(2019) In both cases, very similar sub-100 Torr plasma conditions were
explored. In Tsugawa $et$ $al.$,Tsugawa  _et al._ (2010) nanodiamond film was
obtained at the substrate temperature less than $100\degree$C and in Tallaire
$et$ $al.$,Tallaire  _et al._ (2019) loose nanodiamond film was obtained at
the substrate temperature less than $500\degree$C. In either case, nanodiamond
was drawn onto the substrate when the substrate was cooled significantly
compared to the usual substrate growth temperatures. Low $T_{s}$ would
mitigate the substrate surface kinetics, indirectly highlighting the only
other synthesis route through gas phase kinetics. In the present study, the
substrate could not be cooled to $<500\degree$C and gas and substrate kinetics
were instead divorced by physically separating the plasma and the collecting
substrate through an orifice in the molybdenum plate (see Fig. 2).

Under conventional conditions where NCD or UNCD films, strongly adhering to
the substrate, are grown on a seeded substrate placed into cooler plasma zone
at temperatures 800-900 $\degree$C, both gas and surface kinetics play equally
important role. Particles formed in the gas phase drop on the substrate,
thereby enabling fast renucleation, while their regrowth on the substrate
simultaneously takes place, thereby enhancing overall film growth.

##  VII Conclusion

Barbosa $et$ $al.$, Barbosa  _et al._ (2014) reviewed a consolidation of
reports from which, based on the activation energy analysis, it was concluded
that physico-chemical formation of SCD or MCD (faceted) is different from the
formation of NCD and UNCD (ballaslike). Our work presented new findings that
point out that the formation of any kind of diamond (single crystal, faceted,
ballaslike) is essentially the same and is governed by
CH${}_{\text{3}}^{\bullet}$ reaction as a limiting step. The main difference
is that the formation of single crystal and faceted diamond is dominated by
the substrate surface kinetics, while ballaslike nanodiamond formation takes
place in the plasma , i.e. the nucleation, growth and renucleation is
dominated by the gas phase kinetics. High renucleation happens because the
conditions are not conducive to the formation of a single growing crystal.
When the activation energy of SCD and MCD, calculated using the substrate
temperature, was compared with that of NCD, calculated using the gas
temperature, full agreement between the energies was found. Hence, high
renucleation rate for NCD is not the rate-limiting step.

The comparison of crystal size, morphology, and activation energy illustrated
that the samples produced in the present study are directly comparable with
the ballaslike nanodiamond materials widely reported in the literature. This
warrants the extension of any results and hypotheses from the current work to
nanodiamond materials reported elsewhere. This includes the call for
reconsidering the importance of the gas phase NCD nucleation.

The exact reaction/formation belt in the plasma where nanodiamond particles
(re)nucleate and grow remains unclear. From calculating E${}_{\text{a}}$ using
three gas temperatures estimates, the ideal conditions for nanodiamond
formation and growth are where the zones of CH${}_{\text{3}}^{\bullet}\sim
10^{13}-10^{14}$ cm-3 and $T_{g}\sim 1200-1500$ K overlap. Finding exact
interplays between various reactive species, gas temperature, nucleation and
formation of sp2 and sp3 phases in the gas phase are the subjects of further
studies.

## Acknowledgement

We would like to thank Dr. Mark J. Kushner (University of Michigan) and Dr.
Rebecca Anthony (Michigan State University) for providing valuable insights
into the subject. We are grateful to the U.S. Department of Energy, Office of
Science, Office of Fusion Energy Sciences, Award No. DE-SC0023211, and
National Science Foundation, Division of Chemical, Bioengineering,
Environmental and Transport Systems, Award No. 2333452 for the financial
support.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Derjaguin and Fedoseev (1975) B. V. Derjaguin and D. B. Fedoseev, Scientific American 233, 102 (1975).
  * Harris and Weiner (1990) S. J. Harris and A. M. Weiner, Journal of Applied Physics 67, 6520 (1990).
  * Harris (1990) S. J. Harris, Applied Physics Letters 56, 2298 (1990).
  * Frenklach (1991) M. Frenklach, in  _Proceedings of the Second International Symposium on Diamond Materials_ (Electrochemical Society, 1991) pp. 142–153.
  * Butler  _et al._ (1993) J. E. Butler, R. L. Woodin, L. M. Brown, P. Fallon, A. H. Lettington, and J. W. Steeds, Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences 342, 209 (1993).
  * Wolden  _et al._ (1994) C. Wolden, K. Gleason, and J. Howard, Combustion and Flame 96, 75 (1994).
  * Spitsyn  _et al._ (1981) B. Spitsyn, L. Bouilov, and B. Derjaguin, Journal of Crystal Growth 52, 219 (1981).
  * Maeda  _et al._ (1992) H. Maeda, S. Masuda, K. Kusakabe, and S. Morooka, Journal of Crystal Growth 121, 507 (1992).
  * Pehrsson  _et al._ (1992) P. E. Pehrsson, J. Glesener, and A. Morrish, Thin Solid Films 212, 81 (1992), special Issue: Diamond Films and Related Materials.
  * Kobayashi  _et al._ (1993) K. Kobayashi, T. Nakano, N. Mutsukura, and Y. Machi, Vacuum 44, 1 (1993).
  * Chai  _et al._ (1994) W. Chai, Y. Gu, M. Li, Z. Mai, Q. Li, L. Yuan, and S. Pang, Journal of Crystal Growth 135, 639 (1994).
  * Liu and Dandy (1996a) H. Liu and D. S. Dandy, in  _Diamond Chemical Vapor Deposition: Nucleation and Early Growth Stages_ (Noyes Publications, New Jersey, 1996) Chap. 6, pp. 92–130.
  * Yugo  _et al._ (1991) S. Yugo, T. Kanai, T. Kimura, and T. Muto, Applied Physics Letters 58, 1036 (1991).
  * Garcia Poza  _et al._ (2000) M. Garcia Poza, I. Jiménez, O. Sánchez, C. Gómez-Aleixandre, and L. Vázquez, Phys. Rev. B 61 (2000).
  * Lifshitz  _et al._ (2002) Y. Lifshitz, T. Köhler, T. Frauenheim, I. Guzmann, A. Hoffman, R. Q. Zhang, X. T. Zhou, and S. T. Lee, Science 297, 1531 (2002).
  * Mitura (1987) S. Mitura, Journal of Crystal Growth 80, 417 (1987).
  * Howard  _et al._ (1990) W. Howard, D. Huang, J. Yuan, M. Frenklach, K. E. Spear, R. Koba, and A. W. Phelps, Journal of Applied Physics 68, 1247 (1990).
  * Frenklach  _et al._ (1989) M. Frenklach, R. Kematick, D. Huang, W. Howard, K. E. Spear, A. W. Phelps, and R. Koba, Journal of Applied Physics 66, 395 (1989).
  * Frenklach  _et al._ (1991) M. Frenklach, W. Howard, D. Huang, J. Yuan, K. E. Spear, and R. Koba, Applied Physics Letters 59, 546 (1991).
  * Liu and Dandy (1996b) H. Liu and D. S. Dandy, in  _Diamond Chemical Vapor Deposition: Nucleation and Early Growth Stages_ (Noyes Publications, New Jersey, 1996) Chap. 4, pp. 46–78.
  * Hwang  _et al._ (1996) N. M. Hwang, J. H. Hahn, and D. Y. Yoon, Journal of Crystal Growth 160, 87 (1996).
  * Wang and Frenklach (1991) H. Wang and M. Frenklach, Journal of Applied Physics 70, 7132 (1991).
  * Harris  _et al._ (1991) S. J. Harris, A. M. Weiner, and T. A. Perry, Journal of Applied Physics 70, 1385 (1991).
  * Frenklach (1992) M. Frenklach, The Journal of Chemical Physics 97, 5794 (1992).
  * Garrison  _et al._ (1992) B. J. Garrison, E. J. Dawnkaski, D. Srivastava, and D. W. Brenner, Science 255, 835 (1992).
  * Frenklach (1994) M. Frenklach, MRS Proceedings 339, 255 (1994).
  * Kang and Musgrave (2000) J. K. Kang and C. B. Musgrave, The Journal of Chemical Physics 113, 7582 (2000).
  * Cheesman  _et al._ (2008) A. Cheesman, J. N. Harvey, and M. N. R. Ashfold, The Journal of Physical Chemistry A 112, 11436 (2008).
  * Frenklach and Wang (1991) M. Frenklach and H. Wang, Phys. Rev. B 43, 1520 (1991).
  * Barbosa  _et al._ (2014) D. C. Barbosa, M. R. Baldan, V. J. Trava-Airoldi, and E. J. Corat, Journal of Vacuum Science & Technology B 32, 031808 (2014).
  * Nikhar  _et al._ (2020) T. Nikhar, R. Rechenberg, M. F. Becker, and S. V. Baryshev, Journal of Applied Physics 128, 235305 (2020).
  * Asmussen and Reinhard (2002) J. Asmussen and D. K. Reinhard,  _Diamond Films Handbook_ (New York: Marcel Dekker, 2002) pp. 211–302.
  * Fedoseev  _et al._ (1979) D. Fedoseev, S. Vnukov, and B. Derjaguin, Carbon 17, 453 (1979).
  * Kondoh  _et al._ (1993) E. Kondoh, T. Ohta, T. Mitomo, and K. Ohtsuka, Journal of Applied Physics 73, 3041 (1993).
  * May  _et al._ (2010) P. W. May, J. N. Harvey, N. L. Allan, J. C. Richley, and Y. A. Mankelevich, Journal of Applied Physics 108, 114909 (2010).
  * Mankelevich  _et al._ (2008a) Y. A. Mankelevich, M. N. R. Ashfold, and J. Ma, Journal of Applied Physics 104, 113304 (2008a).
  * Hassouni  _et al._ (2010) K. Hassouni, F. Silva, and A. Gicquel, Journal of Physics D: Applied Physics 43, 153001 (2010).
  * Yoon  _et al._ (2008) J.-S. Yoon, M.-Y. Song, J.-M. Han, S. H. Hwang, W.-S. Chang, B. Lee, and Y. Itikawa, Journal of Physical and Chemical Reference Data 37, 913 (2008).
  * Dandy and Coltrin (1995) D. S. Dandy and M. E. Coltrin, Journal of Materials Research 10, 1993–2010 (1995).
  * Gorbachev and Leshchev (2022) Y. E. Gorbachev and D. V. Leshchev, Fullerenes, Nanotubes and Carbon Nanostructures 30, 76 (2022).
  * Mankelevich  _et al._ (2008b) Y. A. Mankelevich, M. N. R. Ashfold, and J. Ma, Journal of Applied Physics 104, 113304 (2008b).
  * Muehle  _et al._ (2017) M. Muehle, J. Asmussen, M. Becker, and T. Schuelke, Diamond and Related Materials 79, 150 (2017).
  * Xiong  _et al._ (2022) Z. Xiong, S. Lanham, E. Husmann, G. Nelson, M. A. Eslamisaray, J. Polito, Y. Liu, J. Goree, E. Thimsen, M. J. Kushner, and U. R. Kortshagen, Journal of Physics D: Applied Physics 55, 235202 (2022).
  * Tsugawa  _et al._ (2010) K. Tsugawa, M. Ishihara, J. Kim, Y. Koga, and M. Hasegawa, Phys. Rev. B 82, 125460 (2010).
  * Tallaire  _et al._ (2019) A. Tallaire, O. Brinza, M. De Feudis, A. Ferrier, N. Touati, L. Binet, L. Nicolas, T. Delord, G. Hetet, T. Herzig, _et al._ , ACS Applied Nano Materials 2, 5952 (2019).


# Connections between the Open-boundary Spectrum and Generalized Brillouin
Zone in Non-Hermitian Systems

Deguang Wu  National Laboratory of Solid State Microstructures, Department of
Physics, Nanjing University, Nanjing 210093, China  Jiao Xie  National
Laboratory of Solid State Microstructures, Department of Physics, Nanjing
University, Nanjing 210093, China  Yao Zhou  National Laboratory of Solid
State Microstructures, Department of Physics, Nanjing University, Nanjing
210093, China  Jin An<EMAIL_ADDRESS>National Laboratory of Solid State
Microstructures, Department of Physics, Nanjing University, Nanjing 210093,
China  Collaborative Innovation Center of Advanced Microstructures, Nanjing
University, Nanjing 210093, China

###### Abstract

Periodic-boundary spectrum, open-boundary spectrum, as well as the generalized
Brillouin zone(GBZ) are three essential properties of a one-dimensional non-
Hermitian system. In this paper we illustrate that the deep connections
between them can be revealed by a series of special similar transformations.
This viewpoint closely connects the topological geometry of the open-boundary
spectrum with the GBZ and provides a new efficient numerical method of
calculating them accurately. We further extend these connections to non-
Hermitian systems in the symplectic symmetry class. We show that if just the
open-boundary features of a non-Hermitian system such as the spectrum and the
GBZ, are concerned, the relevant symmetry we should consider is not that of
the original system itself, but that of one which has higher symmetry and is
related to the original system by a similarity transformation.

Introduction. —Recent studies on non-Hermitian systems have revealed dramatic
differences from their Hermitian counterpartshatano1997vortex ; hatano1998non
; bender1998real ; bender2002complex ; rudner2009topological ; esaki2011edge ;
hu2011absence ; lee2016anomalous ; leykam2017edge ; shen2018topological ;
yao2018non ; ezawa2019non ; okuma2019topological ; li2021non ; okugawa2021non
. Non-Hermiticity can arise from optical systems with gain and
lossmakris2008beam ; klaiman2008visualization ; guo2009observation ;
longhi2009bloch ; chong2011p ; regensburger2012parity ; bittner2012p ;
liertzer2012pump ; feng2014single ; hodaei2014parity , open systems with
dissipationcarmichael1993quantum ; rotter2009non ; choi2010quasieigenstate ;
lin2011unidirectional ; lee2014heralded ; malzard2015topologically ;
yamamoto2019theory ; li2019observation , and noninteracting or interacting
electron systems with finite-lifetime quasi-particleszyuzin2018flat ;
shen2018quantum ; yoshida2018non ; zhou2018observation ; papaj2019nodal ;
kimura2019chiral . The central feature of non-Hermitian systems is the
sensitivity of the bulk states to boundary conditions, i.e., the non-Hermitian
skin effectyao2018edge ; lee2019anatomy ; lee2019hybrid ; longhi2019probing ;
okuma2020topological ; yi2020non ; kawabata2020higher ; liu2020helical ;
fu2021non ; okuma2021non ; sun2021geometric ; claes2021skin ; guo2021exact ,
which results in the breakdown of the conventional bulk-boundary
correspondencekunst2018biorthogonal ; herviou2019defining ;
imura2019generalized ; lee2019topological ; wang2019non ; jin2019bulk ;
yang2020non ; zhang2020correspondence ; zhang2020bulk ; wang2020defective ;
xue2021simple ; cao2021non . The efforts of resolving this problem then leads
to the establishment of the concept of the GBZyao2018edge  and the non-Bloch
band theoryyao2018edge ; yao2018non ; yokomizo2019non ; kawabata2020non ;
yokomizo2021non . The topological invariants can be redefined on the GBZ
instead of BZ to give an explanation of the anomalous bulk-boundary
correspondenceyao2018edge ; song2019non ; deng2019non . Non-Hermitian
topological systems can also exhibit unique features without Hermitian
analogsgong2018topological ; yin2018geometrical ; jiang2018topological ;
kawabata2019symmetry ; kawabata2019classification ; zhou2019periodic ;
yokomizo2020topological , which can be attributed to the complex-valued nature
of their spectra. The spectrum of a non-Hermitian system under open-boundary
conditions(OBCs) differs greatly from that under periodic-boundary
conditions(PBCs)okuma2020topological ; zhang2020correspondence . The latter
can also define a set of winding numbers to characterize the spectrum of the
system under the half-infinite boundary conditionsokuma2020topological . The
two spectra, together with the winding numbers, are closely related, and are
playing important roles in a non-Hermitian system.

In this paper, we reveal that there are intrinsic connections between the OBC
spectrum, the PBC spectrum and the GBZ. We find the bridge which connects them
are a series of special similarity transformations. We can thus give an
interpretation of the close relation between the geometry of the OBC spectrum
such as 3-bifurcation or 4-bifurcation states with the structure of the GBZ.
Based on our conclusions, we can also provide a new efficient numerical
approach of calculating accurately the OBC spectrum and the GBZ of non-
Hermitian systems. Next, we generalize the conclusions on the connections from
non-Hermitian systems without symmetry to those with anomalous time-reversal
symmetry. We also find that if just the OBC features of a non-Hermitian system
such as the OBC spectrum and GBZ, are concerned, the original system together
with a series of ones which are related to it by the special similarity
transformations, are actually forming an equivalent class which should be
treated as a whole, even though these systems may have different symmetries.
The open-boundary feature is determined by the symmetry of the system which
owns the highest symmetry in the class.

Correspondence between winding number and the OBC spectrum. —We start with a
generic quasi-1D tight-binding model whose Hamiltonian reads

$\displaystyle\mathcal{H}=\sum_{j}\sum_{m=-M_{2}}^{M_{1}}\sum_{\mu,\nu=1}^{N}T_{m}^{\mu\nu}c_{j+m,\mu}^{\dagger}c_{j,\nu},$
(1)

where $\mu$, $\nu$ represent the degrees of freedom including sublattice,
spin, orbitals, etc. $T_{m}^{\mu\nu}$ are the hopping amplitudes, and
$M_{1}$($M_{2}$) is the range of the hopping to the right(left). Thus, in
terms of complex-valued $\beta$, we have the following non-Hermitian Bloch
Hamiltonian,

$\displaystyle H(\beta)=\sum_{m=-M_{2}}^{M_{1}}{T_{m}{\beta}^{-m}},$ (2)

where $T_{m}$ is the $N\times N$ hopping matrix with
$(T_{m})_{\mu\nu}=T_{m}^{\mu\nu}$.

The energy bands $E_{\alpha}(\beta)(\alpha=1,2,\ldots,N)$ of $H(\beta)$ are
determined by the characteristic equation det$(E-H(\beta))=0$. This is an
algebraic equation for $\beta$ of $(p+q)$th order, with $p=N{M}_{1}$ and
$q=N{M}_{2}$, as the determinant is a polynomial of $E$ and $\beta$, which can
be expressed by
det$(E-H(\beta))=\prod^{N}_{\alpha=1}(E-{E}_{\alpha}(\beta))=a_{-p}\frac{1}{{\beta}^{p}}+\ldots+a_{q}{\beta}^{q}$.
For a fixed eigenenergy $E$, the $p+q$ solutions ${\beta}_{i}(E)$ are assumed
to be numbered as
$|{\beta}_{1}(E)|\leq|{\beta}_{2}(E)|\leq\ldots\leq|{\beta}_{p+q}(E)|$. For a
non-Hermitian system without any symmetry, a state belonging to the OBC
continuum bands must have to obey the condition
$|{\beta}_{p}(E)|=|{\beta}_{p+1}(E)|$yao2018non ; yokomizo2019non , and the
trajectories of ${\beta}_{p}$ and ${\beta}_{p+1}$ give the GBZ.

To illustrate the connections between the OBC spectrum and the GBZ in non-
Hermitian systems, we study a series of systems described by Hamiltonians
$\mathcal{H}_{\rho}$, which are related to the matrix Hamiltonian of
$\mathcal{H}$ by a special similarity transformation,

$\displaystyle\mathcal{H}_{\rho}\equiv S^{-1}\mathcal{H}S,$ (3)

where
$S={\rm{diag}}\left\\{1,\rho,\rho^{2},\ldots,\rho^{L}\right\\}\bigotimes\mathbb{I}$
is a block-diagonal matrix, with $\mathbb{I}$ the $N\times N$ identity matrix
and $L$ the system size.

According to Eq.(2), it is clear that ${H}_{\rho}(\beta)=H(\rho\beta)$, as the
hopping matrix $T_{m}$ becomes $T_{m}\rho^{-m}$ under the similarity
transformation. For open systems, since
${\rm{det}}(S^{-1}\mathcal{H}S-E)={\rm{det}}(\mathcal{H}-E)=0$,
$\mathcal{H}_{\rho}$ and $\mathcal{H}$ share identical OBC spectrum. The unit
circle $|\beta|=1$ on the complex $\beta$-plane is the BZ of $H(\beta)$. The
map of this circle onto the complex $E$-plane by $E_{\alpha}(\beta)$ are $N$
curves $E_{\alpha}(|\beta|=1)(\alpha=1,2,\ldots,N)$, which form the PBC
spectrum of $H(\beta)$ and is denoted for simplicity as $E(|\beta|=1)$.
Because the BZ of the system described by ${H}_{\rho}(\beta)$ exactly
corresponds to $|\beta|=\rho$ of the original system described by $H(\beta)$,
the PBC spectrum of ${H}_{\rho}(\beta)$ can be faithfully given by
$E(|\beta|=\rho)$. For a definite $E$, if we denote the solutions of
det$(H(\beta)-E)=0$ as ${\beta}_{i}(E)$, then the solutions of
det$(H(\rho\beta)-E)=0$ must be $\beta_{i}(E)/\rho$, where $i=$1,
2,$\cdot\cdot\cdot,p+q$. Thus for any positive-valued $\rho$,
${H}_{\rho}(\beta)$ shares with $H(\beta)$ the same condition
$|{\beta}_{p}(E)|=|{\beta}_{p+1}(E)|$ for the OBC continuum bands, resulting
in that ${H}_{\rho}(\beta)$ has exactly the same shape of GBZ to that of
${H}(\beta)$, but with a scale factor $\rho^{-1}$.

Due to the multi-valued nature of $E_{\alpha}(\beta)$, each curve of
$E(|\beta|=\rho)$ is not necessarily closed but all the curves must always
form some closed loops with their number being $\mathcal{N}_{\rho}\leq N$. The
$\mathcal{N}_{\rho}$ closed loops divide the whole $E$-plane into several
areas. Each area can be characterized by a winding number which is the sum of
the phase winding of each loop around any point within this area. For any such
point $E$, this winding number can be formally defined by

$\displaystyle w_{\rho}(E)=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}{\oint}_{|\beta|=1}d\beta\frac{d}{d\beta}{\rm{det}}(H_{\rho}(\beta)-E)$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi\mathrm{i}}{\oint}_{|\beta|=\rho}d\beta\frac{d}{d\beta}{\rm{ln}}({\rm{det}}(H(\beta)-E))$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{p+q}\frac{1}{2\pi\mathrm{i}}\oint_{|\beta|=\rho}d\beta{\rm{ln}}(\beta-\beta_{i}(E))-p$
$\displaystyle=$ $\displaystyle N_{\rho}(E)-p,$ (4)

where $N_{\rho}(E)$ is the number of $\beta_{i}(E)$ inside the circle
$|\beta|=\rho$.

We consider the crossing points of these closed loops, which can be self-
intersection ones between the same loop or intersection ones between different
loops. One of our central results is a theorem: For a non-Hermitian system
without any symmetry, for a definite $\rho$, if there exists a crossing point
$E_{c}$ which is on the boundary between two areas with winding number
$w_{\rho}=+1$ and $w_{\rho}=-1$ respectively, then $E_{c}$ belongs to the OBC
spectrum. Conversely, for any “normal” point $E_{s}$ on the OBC spectrum,
there must be one and only one $\rho$ so that the PBC spectrum
$E(|\beta|=\rho)$ has a crossing point exactly at $E_{s}$, with their two
adjacent areas possessing winding numbers $\pm 1$.

Figure 1: Connections between the PBC spectrum and OBC spectrum via winding
numbers. The solid(dotted) lines denote the spectra under PBCs(OBCs). (a)-(b)
are for a single-band model with different $\rho$, where the nonzero
parameters are $T_{-2}=-1$, $T_{-1}=\mathrm{i}$, $T_{1}=-0.5\mathrm{i}$,
$T_{2}=0.5$, while (c)-(d) are for a two-band model, where the nonzero
parameters are $T_{-1}=\mathbb{I}+2\sigma_{3}$,
$T_{0}=-\mathbb{I}+\mathrm{i}\sigma_{2}$, $T_{1}=2\mathbb{I}+\sigma_{3}$, with
$\sigma_{i}(i=1,2,3)$ Pauli matrices. The colored areas are possessing winding
number $w_{\rho}(E)=\pm 1$.

In Fig.1 we illustrate the essential points of the above result by using both
a single-band and two-band model as examples. In the single-band model, even
one closed loop give multiple self-crossing points which are located on the
OBC spectrum, while in the two-band model, for a definite $\rho$, the PBC
spectrum $E(|\beta|=\rho)$ give simultaneously the self-crossing points and
intersection points between the two loops, which also belong to the
corresponding OBC spectrum.

Now we give the proof. Let $E_{c}$ be a crossing point of a PBC spectrum
$E(|\beta|=\rho)$. This crossing point means that there must be two different
$\beta_{c}$ and $\beta^{{}^{\prime}}_{c}$ on the circle $|\beta|=\rho$( see
Fig.2(a)), satisfying
$E_{\alpha}(\beta_{c})=E_{\gamma}(\beta^{{}^{\prime}}_{c})=E_{c}$ for two
different bands $\alpha$ and $\gamma$, or
$E_{\alpha}(\beta_{c})=E_{\alpha}(\beta^{{}^{\prime}}_{c})=E_{c}$ for the same
band $\alpha$. Adjacent to $E_{c}$, there are four areas with winding numbers
$(w+1)/w/(w-1)/w$ respectively, as shown in Fig 2.(b). We have to prove that
if $w=0$, $E_{c}$ is on the OBC spectrum. Consider a small curve starting from
point $A$ to point B, which is passing through $E_{c}$. The area where $A/B$
sits possesses winding number $(w+1)/(w-1)$. Due to Eq.( 4), the number of
$\beta_{i}(E_{A/B})$ inside the circle $|\beta|=\rho$ should be
$(p+w+1)/(p+w-1)$. If $A$ and $B$ are sufficiently close to $E_{c}$, among the
$p+q$ image curves of the map $\beta_{i}(E)$ of the trajectory
$E_{A}\rightarrow E_{c}\rightarrow E_{B}$ onto the $\beta$-plane, only two
small curves would intersect the circle $|\beta|=\rho$ at $\beta=\beta_{c}$
and $\beta=\beta^{{}^{\prime}}_{c}$. Moreover, these two curves must be
pointing outwards, as shown in Fig 2.(a), in order to reduce the winding
number by two. Therefore, the number of $\beta_{i}(E_{c})$ inside the circle
should be $p+w-1$. If $w=0$, then $\beta_{c}$ and $\beta^{{}^{\prime}}_{c}$
must be the $p$th and $(p+1)$th $\beta$ solutions of $E_{c}$ and then their
equal absolute value means that $E_{c}$ is on the OBC spectrum.

Conversely, consider a point $E_{s}$ on the OBC spectrum. This $E_{s}$
generally corresponds to two different $\beta_{p}(E_{s})$,
$\beta_{p+1}(E_{s})$ on the GBZ. Since
$|\beta_{p}(E_{s})|=|\beta_{p+1}(E_{s})|=\rho_{s}$, the PBC spectrum
$E(|\beta|=\rho_{s})$ would inevitable form a crossing point at $E_{s}$.
Obviously, the number of $\beta_{i}(E_{s})$ inside the circle
$|\beta|=\rho_{s}$ is $p-1$, and for the four areas adjacent to $E_{s}$,
similar to the situation as shown in Fig.2(b), their winding numbers must be
$(+1)/0/(-1)/0$. Thus we complete the proof.

Figure 2: Illustration of the theorem: (a) $|\beta|=\rho$ and (b) a crossing
point $E_{c}$ of the corresponding PBC spectrum $E(|\beta|=\rho)$, where the
four areas adjacent to $E_{c}$ have been marked by their winding numbers. The
spectrum on the $E$-plane is the image of the circle on the $\beta$-plane,
under the mapping $E(\beta)$. The two solid dots $\beta_{c}$,
$\beta^{{}^{\prime}}_{c}$ on the circle in (a) are the preimages of $E_{c}$ in
(b), and accordingly, the two small arcs passing through the two points
corresponds to the trajectory of from A to B in (b).

Therefore, any crossing point of a PBC spectrum of $H_{\rho}(\beta)$ with $\pm
1$ winding numbers in its two adjacent areas must be on the OBC spectrum. This
theorem thus provides an accurate and effective new method to numerically
obtain the OBC spectrum as well as the GBZ for 1D non-Hermitian systemsSupp1 .
Compared with the previous worksyokomizo2019non ; yang2020non , this numerical
method is much easier to handle.

Obviously, not every crossing point is on the OBC spectrum, because some
crossing points are not adjacent to the areas with winding numbers $\pm 1$.
These energy states $E$ generally require $|\beta_{i}(E)|=|\beta_{i+1}(E)|$
with $i\neq p$, and thus the two $\beta$ values should be on the auxiliary GBZ
introduced in Ref.yang2020non  in obtaining the GBZ.

An exception to the theorem is when the GBZ is a perfect circle $|\beta|=r$.
This can happen even in multi-band models. In this particular situation, the
PBC spectrum $E(|\beta|=r)$ covers the whole OBC spectrum, and the other PBC
spectra $E(|\beta|=\rho)$ with $\rho\neq r$ give no contributions to the OBC
spectrumSupp2 .

N-bifurcation states, topological geometry of the OBC spectrum and their
connection with the GBZ. —According to the theorem, nearly any state on the
OBC spectrum can be a crossing point of a PBC spectrum for a certain $\rho$.
These are the “normal” points mentioned in the theorem, which nearly
constitute the whole OBC spectrum and are characterized by their unique
feature that each has two and only two extending directions along the
spectrum. However, there always exist several other energy states on the OBC
spectrum, which are the 3-bifurcation or 4-bifurcation statesSupp3 (see the
solid dots in Fig.3(a) and (c)), and end-point states. These particular states
can also be the crossing points of some PBC spectra, but with different
winding number distributions in their adjacent areas. Take a 3-bifurcation
state $E$ as an example. We find it obeys,

$\displaystyle|\beta_{p-1}(E)|=|\beta_{p}(E)|=|\beta_{p+1}(E)|,$
$\displaystyle{\rm{or}}\qquad|\beta_{p}(E)|=|\beta_{p+1}(E)|=|\beta_{p+2}(E)|,$
(5)

i.e., each 3-bifurcation state has not two, but three $\beta$ solutions on the
GBZ. More generally, a $n$-bifurcation state $E$ must have exactly $n$ $\beta$
solutions on the GBZ. In this opinion, all the energy states can be viewed in
an universal way. The end points on the OBC spectrum can be viewed as
“1-bifurcation” ones and the “normal” points as the 2-bifurcation ones. The
single $\beta$ on the GBZ corresponding to an end-point state $E$ is actually
a doubly degenerate $\beta$ satisfying $\beta_{p}=\beta_{p+1}$, as discussed
in Ref.yang2020non . These behavior are exhibited in Fig.3. In the one-band
model with $p=q=4$, for the 3-bifurcation states, we have
$|\beta_{3}|=|\beta_{4}|=|\beta_{5}|$, and for the 4-bifurcation state, we
have $|\beta_{4}|=|\beta_{5}|=|\beta_{6}|=|\beta_{7}|$. In the two-band model
with $p=q=2$, for the four 3-bifurcation states, we have
$|\beta_{1}|=|\beta_{2}|=|\beta_{3}|$ for $A_{1}$ and $A_{2}$, while
$|\beta_{2}|=|\beta_{3}|=|\beta_{4}|$ for $B_{1}$ and $B_{2}$.

To prove Eq.(5), for a 3-bifurcation state $E$, assume its counterpart on the
GBZ is a set of $n$ $\beta_{i}(E)$ points, where $i=j+1,j+2,...j+n$, obeying
$|\beta_{j+1}(E)|=|\beta_{j+2}(E)|=...=|\beta_{j+n}(E)|$, with $j<p$ and
$j+n>p$. Consider three small arcs along the three branches connecting to $E$
(see Fig.3(c)). Because any other state on the arcs is a normal one which
corresponds to two different $\beta$s on the GBZ, the three small arcs must be
mapped onto the GBZ to the six small arcs departing from the $n$
$\beta_{i}(E)$. But departing from any $\beta_{i}(E)$, there exist only two
extending directions to accommodate two small arcs. Then $n$ must be three and
we complete the proof.

Figure 3: N-bifurcation states on the OBC spectrum and the GBZ. (a) and (c):
3-bifurcation and 4-bifurcation states, on the OBC spectrum, together with the
PBC spectra passing through them. (b) and (d): GBZ, where the solid circles
are the corresponding three(four) $\beta$s for the 3(4)-bifurcation states.
Here (a)-(b) are for a single-band model with the nonzero parameters
$T_{-4}=0.8$, $T_{\pm{3}}=\mp{1}$, $T_{\pm{2}}=\pm{1}$, $T_{-1}=2$, $T_{1}=3$,
$T_{4}=0.3$, while (c)-(d) are for a two-band model with the nonzero
parameters are $T_{-1}=1.5\mathbb{I}+\sigma_{1}+0.5\sigma_{3}$,
$T_{0}=\mathrm{i}\sigma_{2}$,
$T_{1}=0.5\mathbb{I}+\sigma_{1}-2.5\sigma_{3}$,where the PBC spectrum passing
through the 3-bifurcation state $A_{1}(B_{1})$ is also passing through
$A_{2}(B_{2})$. The dotted loops in (b) and (d) are the circles $|\beta|=\rho$
corresponding to the PBC spectra in (a) and (c). (e) All possible winding
number distributions adjacent to a 3-bifurcation or 4-bifurcation state.

For a generic $n$-bifurcation state $E$, since there are $n$ successive
$\beta_{i}(E)$ solutions with the same absolute value located on the GBZ, the
PBC spectrum $E(|\beta|=\rho_{0})$ is passing through $E$ exactly $n$ times,
where $\rho_{0}=|\beta_{p}(E)|$. The local areas adjacent to $E$ would be
divided into $2n$ parts. A $n$-bifurcation state generally has relatively
definite winding number distributions in these adjacent areas. For a
3-bifurcation or 4-bifurcation energy state, all possible winding number
distributions are shown explicitly in Fig.3(e). The $n$ $\beta_{i}(E)$
solutions on the GBZ can also be viewed in another way: They are the
intersection points between the circle $|\beta|=\rho$ and the GBZ, as seen in
Fig.3(b) and Fig.3(d). So starting with a circle $|\beta|=\rho$, with $\rho$
increasing gradually from $0$ to infinity, one has a set of intersection
points for each $\rho$ with the GBZ. Each set determines one or several energy
states $E$ on the OBC spectrum and the number of the points in the set
indicates the topological geometry of this state: One only means an end-point
state, two means a “normal” state or two end-point states, and three means a
3-bifurcation state, or a “normal” state together with an end-point state, or
three end-point states, etc. Thus the $\rho$ parameter of the PBC spectrum
having crossing points on the open spectrum has actually a finite range:
$|\beta|_{\rm{min}}\leq\rho\leq|\beta|_{\rm{max}}$, where
$|\beta|_{\rm{max}/\rm{min}}$ is the maximum/minimum value of $|\beta|$ on the
GBZ. Now we can also give the non-Hermitian skin effect a new viewpoint. Under
OBCs, any state belonging to the continuum bands generally will be localized
at end. However, if we assign each state $E$ on the OBC spectrum an unique
open system described by $H_{\rho}(\beta)$, with $\rho=|\beta_{p}(E)|$, then
this state with respect to $H_{\rho}(\beta)$ must be extended, since its
$\beta$ solutions $\beta^{{}^{\prime}}_{i}(E)=\beta_{i}(E)/\rho$ obey:
$|\beta^{{}^{\prime}}_{p}(E)|=|\beta^{{}^{\prime}}_{p+1}(E)|=1$. Thus from
this opinion, the whole OBC spectrum of the continuum bands can still be
viewed as being extended without skin effect, not for the original system
$H(\beta)$, but for a series of representatives $H_{\rho}(\beta)$ with
$\rho=|\beta_{p}(E)|$.

Extension to systems with anomalous time-reversal symmetry. —Now we extend our
results to systems with symmetry. We focus on the non-Hermitian systems in the
symplectic symmetry class, which has the anomalous time-reversal
symmetryyi2020non ; kawabata2020non ; kawabata2019symmetry  defined by:
$U_{T}\mathcal{H}^{t}U^{-1}_{T}=\mathcal{H}$, where $U_{T}$ is a unitary
matrix satisfying $U_{T}{U_{T}}^{*}=-1$. This symmetry indicates
$U_{T}H^{t}(\beta^{-1})U^{-1}_{T}=H(\beta)$, so the standard non-Bloch band
theory breaks down and should be modified according to Ref.kawabata2020non .
The characteristic equation det$[H(\beta)-E]=0$ has $2p$ numbered solutions
$\beta_{i}(E)$ in total( $p=q=NM_{1}$), obeying
$\beta_{p+i}(E)=\beta_{p+1-i}^{-1}(E)$ due to Kramers
degeneracykawabata2020non , which results in a new condition for the GBZ:
$|\beta_{p-1}(E)|=|\beta_{p}(E)|$ or $|\beta_{p+1}(E)|=|\beta_{p+2}(E)|$Supp4
.

Analogous discussion leads us to a similar conclusion for a non-Hermitian
system in the symplectic symmetry class: The crossing points of the PBC
spectrum $E(|\beta|=\rho)$ with its two adjacent areas possessing winding
numbers $\pm 2/0$ belong to the OBC spectrum. Conversely, for any “normal”
point $E_{s}$ on the OBC spectrum, these must exist two $\rho_{1}$ and
$\rho_{2}$ with $\rho_{1}\rho_{2}=1$ and $\rho_{1}\leq\rho_{2}$ so that the
PBC spectrum $E(|\beta|=\rho_{1})$ ( $E(|\beta|=\rho_{2})$) has a crossing
point exactly at $E_{s}$, with its two adjacent areas possessing winding
numbers $-2/0$(+2/0).

Figure 4: A non-Hermitian system in the symplectic symmetry class. (a) OBC
spectrum, PBC spectra for a two-band model with the non-zero parameters
$T_{\pm{1}}=2\mathbb{I}\pm(1+1.5\mathrm{i})\sigma_{1}\pm\sigma_{2}\pm
2\mathrm{i}\sigma_{3}$,
$T_{\pm{2}}=3\mathrm{i}\mathbb{I}\pm(-0.6+0.2\mathrm{i})\sigma_{1}\pm\sigma_{2}\pm
0.8\sigma_{3}$. Here the crossing points of the PBC spectrum in (a) have two
adjacent areas with winding number $0$ or $\pm 2$ and so they belong to the
OBC spectrum. (c) A 3-bifurcation state on the OBC spectrum and the PBC
spectra passing through it, and (d) its counterpart on the GBZ, where the
nonzero parameters are $T_{\pm{1}}=2\mathbb{I}\pm 1.5\mathrm{i}\sigma_{1}\pm
1.5\mathrm{i}\sigma_{2}\pm\sigma_{3}$, $T_{\pm{2}}=\mathrm{i}\mathbb{I}\pm
0.7\mathrm{i}\sigma_{1}\pm 0.7\mathrm{i}\sigma_{2}\pm 0.5\sigma_{3}$. In both
models, $U_{T}=\sigma_{2}$.

Now we give our explanation. Due to the symmetry
$U_{T}H^{t}(\beta^{-1})U^{-1}_{T}=H(\beta)$, the $\mathcal{N}_{\rho}$ closed
loops of the PBC spectrum $E(|\beta|=\rho)$ exactly coincide with those of the
PBC spectrum $E(|\beta|=\rho^{-1})$, but the two sets of loops are pointing in
opposite directions. So each state on the OBC spectrum will be experienced
twice by the PBC ones. This also indicates that any area on the $E-$plane has
opposite winding number with respect to $E(|\beta|=\rho)$ and
$E(|\beta|=\rho^{-1})$. Since here a “normal” state on the OBC spectrum
corresponds to four $\beta$s on the GBZ, an end-point state should have two
degenerate points on the GBZ obeying: $\beta_{p-1}=\beta_{p}$ or
$\beta_{p}=\beta_{p+1}$. The counterpart of a 3-bifurcation state should have
two sets of $\beta$ points on the GBZ, with each set being composed of the
three of them, which obey:

$\displaystyle|\beta_{p-2}(E)|=|\beta_{p-1}(E)|=|\beta_{p}(E)|,$
$\displaystyle{\rm{or}}\qquad|\beta_{p+1}(E)|=|\beta_{p+2}(E)|=|\beta_{p+3}(E)|.$
(6)

In this case, there still exist two PBC spectra for $\rho_{1}$ and $\rho_{2}$
with $\rho_{1}\rho_{2}=1$, each of which is passing through the 3-bifurcation
state three times and is overlapping exactly with each other. We illustrate
these results by a two-band model in Fig.4, where $p=4$, and the 3-bifurcation
state obeys: $|\beta_{2}|=|\beta_{3}|=|\beta_{4}|$, or
$|\beta_{5}|=|\beta_{6}|=|\beta_{7}|$.

It should be noted that $H_{\rho}(\beta)$ with $\rho\neq 1$ generally have no
symmetry, but they share identical OBC spectrum with $H(\beta)$ since they are
related to each other by similarity transformations. What’s remarkable is that
$H_{\rho}(\beta)$ share with $H(\beta)$ the same GBZ conditions, even if they
possess different symmetries. This is because, as mentioned before, for the
same state with eigenenergy $E$, the $\beta_{i}(E)$ solutions of
det$(E-H_{\rho}(\beta))=0$ is proportional to those of det$(E-H(\beta))=0$. To
make this point more clear, we consider a class of Hamiltonians related by the
special similarity transformations,
$\mathcal{C}_{H}=\\{H_{\rho}(\beta)=H(\rho\beta)\mid 0<\rho<\infty\\}$. If the
original system we start with is described by $H(\beta)$ which does not have
any symmetry, and simultaneously, one member $H_{\rho}(\beta)$ of the class
has certain symmetry which leads to a modified GBZ conditions, then $H(\beta)$
must have the same GBZ conditions. Therefore, if only the OBC features of a
non-Hermitian system such as the spectrum and the GBZ, are concerned, the
relevant symmetry we should consider is not that of the original system
$H(\beta)$ itself, but the one belonging to $\mathcal{C}_{H}$ which has the
highest symmetry. So to extract the properties of an open non-Hermitian
systems, the class $\mathcal{C}_{H}$ should be treated as a whole. This
viewpoint could be extended straightforwardly to higher dimensions.

This work is supported by NSFC under Grants No.11874202.

## References

  * [1] Naomichi Hatano and David R Nelson.  Vortex pinning and non-hermitian quantum mechanics.  Physical Review B, 56(14):8651, 1997.
  * [2] Naomichi Hatano and David R Nelson.  Non-hermitian delocalization and eigenfunctions.  Physical Review B, 58(13):8384, 1998.
  * [3] Carl M Bender and Stefan Boettcher.  Real spectra in non-hermitian hamiltonians having p t symmetry.  Physical Review Letters, 80(24):5243, 1998.
  * [4] Carl M Bender, Dorje C Brody, and Hugh F Jones.  Complex extension of quantum mechanics.  Physical Review Letters, 89(27):270401, 2002.
  * [5] Mark S Rudner and LS Levitov.  Topological transition in a non-hermitian quantum walk.  Physical Review Letters, 102(6):065703, 2009.
  * [6] Kenta Esaki, Masatoshi Sato, Kazuki Hasebe, and Mahito Kohmoto.  Edge states and topological phases in non-hermitian systems.  Physical Review B, 84(20):205128, 2011.
  * [7] Yi Chen Hu and Taylor L Hughes.  Absence of topological insulator phases in non-hermitian p t-symmetric hamiltonians.  Physical Review B, 84(15):153101, 2011.
  * [8] Tony E Lee.  Anomalous edge state in a non-hermitian lattice.  Physical Review Letters, 116(13):133903, 2016.
  * [9] Daniel Leykam, Konstantin Y Bliokh, Chunli Huang, Yi Dong Chong, and Franco Nori.  Edge modes, degeneracies, and topological numbers in non-hermitian systems.  Physical Review Letters, 118(4):040401, 2017.
  * [10] Huitao Shen, Bo Zhen, and Liang Fu.  Topological band theory for non-hermitian hamiltonians.  Physical Review Letters, 120(14):146402, 2018.
  * [11] Shunyu Yao, Fei Song, and Zhong Wang.  Non-hermitian chern bands.  Physical Review Letters, 121(13):136802, 2018.
  * [12] Motohiko Ezawa.  Non-hermitian boundary and interface states in nonreciprocal higher-order topological metals and electrical circuits.  Physical Review B, 99(12):121411, 2019.
  * [13] Nobuyuki Okuma and Masatoshi Sato.  Topological phase transition driven by infinitesimal instability: Majorana fermions in non-hermitian spintronics.  Physical Review Letters, 123(9):097701, 2019.
  * [14] Tianyu Li, Jia-Zheng Sun, Yong-Sheng Zhang, and Wei Yi.  Non-bloch quench dynamics.  Physical Review Research, 3(2):023022, 2021.
  * [15] Ryo Okugawa, Ryo Takahashi, and Kazuki Yokomizo.  Non-hermitian band topology with generalized inversion symmetry.  Physical Review B, 103(20):205205, 2021.
  * [16] Konstantinos G Makris, R El-Ganainy, DN Christodoulides, and Ziad H Musslimani.  Beam dynamics in p t symmetric optical lattices.  Physical Review Letters, 100(10):103904, 2008.
  * [17] Shachar Klaiman, Uwe Günther, and Nimrod Moiseyev.  Visualization of branch points in p t-symmetric waveguides.  Physical Review Letters, 101(8):080402, 2008.
  * [18] A Guo, GJ Salamo, D Duchesne, R Morandotti, M Volatier-Ravat, V Aimez, GA Siviloglou, and DN Christodoulides.  Observation of p t-symmetry breaking in complex optical potentials.  Physical Review Letters, 103(9):093902, 2009.
  * [19] Stefano Longhi.  Bloch oscillations in complex crystals with p t symmetry.  Physical Review Letters, 103(12):123601, 2009.
  * [20] YD Chong, Li Ge, and A Douglas Stone.  P t-symmetry breaking and laser-absorber modes in optical scattering systems.  Physical Review Letters, 106(9):093902, 2011.
  * [21] Alois Regensburger, Christoph Bersch, Mohammad-Ali Miri, Georgy Onishchukov, Demetrios N Christodoulides, and Ulf Peschel.  Parity–time synthetic photonic lattices.  Nature, 488(7410):167–171, 2012.
  * [22] S Bittner, B Dietz, Uwe Günther, HL Harney, M Miski-Oglu, A Richter, and F Schäfer.  P t symmetry and spontaneous symmetry breaking in a microwave billiard.  Physical Review Letters, 108(2):024101, 2012.
  * [23] M Liertzer, Li Ge, A Cerjan, AD Stone, Hakan E Türeci, and S Rotter.  Pump-induced exceptional points in lasers.  Physical Review Letters, 108(17):173901, 2012.
  * [24] Liang Feng, Zi Jing Wong, Ren-Min Ma, Yuan Wang, and Xiang Zhang.  Single-mode laser by parity-time symmetry breaking.  Science, 346(6212):972–975, 2014.
  * [25] Hossein Hodaei, Mohammad-Ali Miri, Matthias Heinrich, Demetrios N Christodoulides, and Mercedeh Khajavikhan.  Parity-time–symmetric microring lasers.  Science, 346(6212):975–978, 2014.
  * [26] HJ Carmichael.  Quantum trajectory theory for cascaded open systems.  Physical Review Letters, 70(15):2273, 1993.
  * [27] Ingrid Rotter.  A non-hermitian hamilton operator and the physics of open quantum systems.  Journal of Physics A: Mathematical and Theoretical, 42(15):153001, 2009.
  * [28] Youngwoon Choi, Sungsam Kang, Sooin Lim, Wookrae Kim, Jung-Ryul Kim, Jai-Hyung Lee, and Kyungwon An.  Quasieigenstate coalescence in an atom-cavity quantum composite.  Physical Review Letters, 104(15):153601, 2010.
  * [29] Zin Lin, Hamidreza Ramezani, Toni Eichelkraut, Tsampikos Kottos, Hui Cao, and Demetrios N Christodoulides.  Unidirectional invisibility induced by p t-symmetric periodic structures.  Physical Review Letters, 106(21):213901, 2011.
  * [30] Tony E Lee and Ching-Kit Chan.  Heralded magnetism in non-hermitian atomic systems.  Physical Review X, 4(4):041001, 2014.
  * [31] Simon Malzard, Charles Poli, and Henning Schomerus.  Topologically protected defect states in open photonic systems with non-hermitian charge-conjugation and parity-time symmetry.  Physical Review Letters, 115(20):200402, 2015.
  * [32] Kazuki Yamamoto, Masaya Nakagawa, Kyosuke Adachi, Kazuaki Takasan, Masahito Ueda, and Norio Kawakami.  Theory of non-hermitian fermionic superfluidity with a complex-valued interaction.  Physical Review Letters, 123(12):123601, 2019.
  * [33] Jiaming Li, Andrew K Harter, Ji Liu, Leonardo de Melo, Yogesh N Joglekar, and Le Luo.  Observation of parity-time symmetry breaking transitions in a dissipative floquet system of ultracold atoms.  Nature communications, 10(1):1–7, 2019.
  * [34] Alexander A Zyuzin and A Yu Zyuzin.  Flat band in disorder-driven non-hermitian weyl semimetals.  Physical Review B, 97(4):041203, 2018.
  * [35] Huitao Shen and Liang Fu.  Quantum oscillation from in-gap states and a non-hermitian landau level problem.  Physical Review Letters, 121(2):026403, 2018.
  * [36] Tsuneya Yoshida, Robert Peters, and Norio Kawakami.  Non-hermitian perspective of the band structure in heavy-fermion systems.  Physical Review B, 98(3):035141, 2018.
  * [37] Hengyun Zhou, Chao Peng, Yoseob Yoon, Chia Wei Hsu, Keith A Nelson, Liang Fu, John D Joannopoulos, Marin Soljačić, and Bo Zhen.  Observation of bulk fermi arc and polarization half charge from paired exceptional points.  Science, 359(6379):1009–1012, 2018.
  * [38] Michał Papaj, Hiroki Isobe, and Liang Fu.  Nodal arc of disordered dirac fermions and non-hermitian band theory.  Physical Review B, 99(20):201107, 2019.
  * [39] Kazuhiro Kimura, Tsuneya Yoshida, and Norio Kawakami.  Chiral-symmetry protected exceptional torus in correlated nodal-line semimetals.  Physical Review B, 100(11):115124, 2019.
  * [40] Shunyu Yao and Zhong Wang.  Edge states and topological invariants of non-hermitian systems.  Physical Review Letters, 121(8):086803, 2018.
  * [41] Ching Hua Lee and Ronny Thomale.  Anatomy of skin modes and topology in non-hermitian systems.  Physical Review B, 99(20):201103, 2019.
  * [42] Ching Hua Lee, Linhu Li, and Jiangbin Gong.  Hybrid higher-order skin-topological modes in nonreciprocal systems.  Physical Review Letters, 123(1):016805, 2019.
  * [43] Stefano Longhi.  Probing non-hermitian skin effect and non-bloch phase transitions.  Physical Review Research, 1(2):023013, 2019.
  * [44] Nobuyuki Okuma, Kohei Kawabata, Ken Shiozaki, and Masatoshi Sato.  Topological origin of non-hermitian skin effects.  Physical Review Letters, 124(8):086801, 2020.
  * [45] Yifei Yi and Zhesen Yang.  Non-hermitian skin modes induced by on-site dissipations and chiral tunneling effect.  Physical Review Letters, 125(18):186802, 2020.
  * [46] Kohei Kawabata, Masatoshi Sato, and Ken Shiozaki.  Higher-order non-hermitian skin effect.  Physical Review B, 102(20):205118, 2020.
  * [47] Chun-Hui Liu, Kai Zhang, Zhesen Yang, and Shu Chen.  Helical damping and dynamical critical skin effect in open quantum systems.  Physical Review Research, 2(4):043167, 2020.
  * [48] Yongxu Fu, Jihan Hu, and Shaolong Wan.  Non-hermitian second-order skin and topological modes.  Physical Review B, 103(4):045420, 2021.
  * [49] Nobuyuki Okuma and Masatoshi Sato.  Non-hermitian skin effects in hermitian correlated or disordered systems: Quantities sensitive or insensitive to boundary effects and pseudo-quantum-number.  Physical Review Letters, 126(17):176601, 2021.
  * [50] Xiao-Qi Sun, Penghao Zhu, and Taylor L Hughes.  Geometric response and disclination-induced skin effects in non-hermitian systems.  arXiv preprint arXiv:2102.05667, 2021.
  * [51] Jahan Claes and Taylor L Hughes.  Skin effect and winding number in disordered non-hermitian systems.  Physical Review B, 103(14):L140201, 2021.
  * [52] Cui-Xian Guo, Chun-Hui Liu, Xiao-Ming Zhao, Yanxia Liu, and Shu Chen.  Exact solution of non-hermitian systems with generalized boundary conditions: size-dependent boundary effect and fragility of skin effect.  arXiv preprint arXiv:2102.03781, 2021.
  * [53] Flore K Kunst, Elisabet Edvardsson, Jan Carl Budich, and Emil J Bergholtz.  Biorthogonal bulk-boundary correspondence in non-hermitian systems.  Physical Review Letters, 121(2):026808, 2018.
  * [54] Loic Herviou, Jens H Bardarson, and Nicolas Regnault.  Defining a bulk-edge correspondence for non-hermitian hamiltonians via singular-value decomposition.  Physical Review A, 99(5):052118, 2019.
  * [55] Ken-Ichiro Imura and Yositake Takane.  Generalized bulk-edge correspondence for non-hermitian topological systems.  Physical Review B, 100(16):165430, 2019.
  * [56] Jong Yeon Lee, Junyeong Ahn, Hengyun Zhou, and Ashvin Vishwanath.  Topological correspondence between hermitian and non-hermitian systems: Anomalous dynamics.  Physical Review Letters, 123(20):206404, 2019.
  * [57] Huaiqiang Wang, Jiawei Ruan, and Haijun Zhang.  Non-hermitian nodal-line semimetals with an anomalous bulk-boundary correspondence.  Physical Review B, 99(7):075130, 2019.
  * [58] L Jin and Z Song.  Bulk-boundary correspondence in a non-hermitian system in one dimension with chiral inversion symmetry.  Physical Review B, 99(8):081103, 2019.
  * [59] Zhesen Yang, Kai Zhang, Chen Fang, and Jiangping Hu.  Non-hermitian bulk-boundary correspondence and auxiliary generalized brillouin zone theory.  Physical Review Letters, 125(22):226402, 2020.
  * [60] Kai Zhang, Zhesen Yang, and Chen Fang.  Correspondence between winding numbers and skin modes in non-hermitian systems.  Physical Review Letters, 125(12):126402, 2020.
  * [61] Zhicheng Zhang, Zhesen Yang, and Jiangping Hu.  Bulk-boundary correspondence in non-hermitian hopf-link exceptional line semimetals.  Physical Review B, 102(4):045412, 2020.
  * [62] Xiao-Ran Wang, Cui-Xian Guo, and Su-Peng Kou.  Defective edge states and number-anomalous bulk-boundary correspondence in non-hermitian topological systems.  Physical Review B, 101(12):121116, 2020.
  * [63] Wen-Tan Xue, Ming-Rui Li, Yu-Min Hu, Fei Song, and Zhong Wang.  Simple formulas of directional amplification from non-bloch band theory.  Physical Review B, 103(24):L241408, 2021.
  * [64] Yang Cao, Yang Li, and Xiaosen Yang.  Non-hermitian bulk-boundary correspondence in a periodically driven system.  Physical Review B, 103(7):075126, 2021.
  * [65] Kazuki Yokomizo and Shuichi Murakami.  Non-bloch band theory of non-hermitian systems.  Physical Review Letters, 123(6):066404, 2019.
  * [66] Kohei Kawabata, Nobuyuki Okuma, and Masatoshi Sato.  Non-bloch band theory of non-hermitian hamiltonians in the symplectic class.  Physical Review B, 101(19):195147, 2020.
  * [67] Kazuki Yokomizo and Shuichi Murakami.  Non-bloch band theory in bosonic bogoliubov–de gennes systems.  Physical Review B, 103(16):165123, 2021.
  * [68] Fei Song, Shunyu Yao, and Zhong Wang.  Non-hermitian topological invariants in real space.  Physical Review Letters, 123(24):246801, 2019.
  * [69] Tian-Shu Deng and Wei Yi.  Non-bloch topological invariants in a non-hermitian domain wall system.  Physical Review B, 100(3):035102, 2019.
  * [70] Zongping Gong, Yuto Ashida, Kohei Kawabata, Kazuaki Takasan, Sho Higashikawa, and Masahito Ueda.  Topological phases of non-hermitian systems.  Physical Review X, 8(3):031079, 2018.
  * [71] Chuanhao Yin, Hui Jiang, Linhu Li, Rong Lü, and Shu Chen.  Geometrical meaning of winding number and its characterization of topological phases in one-dimensional chiral non-hermitian systems.  Physical Review A, 97(5):052115, 2018.
  * [72] Hui Jiang, Chao Yang, and Shu Chen.  Topological invariants and phase diagrams for one-dimensional two-band non-hermitian systems without chiral symmetry.  Physical Review A, 98(5):052116, 2018.
  * [73] Kohei Kawabata, Ken Shiozaki, Masahito Ueda, and Masatoshi Sato.  Symmetry and topology in non-hermitian physics.  Physical Review X, 9(4):041015, 2019.
  * [74] Kohei Kawabata, Takumi Bessho, and Masatoshi Sato.  Classification of exceptional points and non-hermitian topological semimetals.  Physical Review Letters, 123(6):066405, 2019.
  * [75] Hengyun Zhou and Jong Yeon Lee.  Periodic table for topological bands with non-hermitian symmetries.  Physical Review B, 99(23):235112, 2019.
  * [76] Kazuki Yokomizo and Shuichi Murakami.  Topological semimetal phase with exceptional points in one-dimensional non-hermitian systems.  Physical Review Research, 2(4):043045, 2020.
  * [77] See Supplemental Material at [URL will be inserted by publisher] for (i) numberical method of calculating the OBC spectrum and GBZ, (ii) Non-Hermitian systems with their GBZ being a perfect circle, (iii) proof of correspondence between winding number and OBC spectrum in anomalous time-reversal symmetry systems, (iv) reality of the 3-bifurcation states.
  * [78] We have provided some models as examples of this exceptional case in Supplemental Material [77].
  * [79] For a simple model, we have a detailed proof of the existence of the 3-bifuraction states on the open-boundary condition spectrum, see Supplemental Material [77].
  * [80] The proof is similar to the previous one in proving the theorem, see Supplemental Material [77] for detail.