# Uncoupling System and Environment Simulation Cells for Fast-Scaling Modeling
of Complex Continuum Embeddings

G. Medrano, E. Bainglass, O. Andreussi

###### Abstract

Continuum solvation models are becoming increasingly relevant in condensed
matter simulations, allowing to characterize materials interfaces in the
presence of wet electrified environments at a reduced computational cost with
respect to all atomistic simulations. However, some challenges with the
implementation of these models in plane-wave simulation packages still
persists, especially when the goal is to simulate complex and heterogeneous
environments. Among these challenges is the computational cost associated with
large heterogeneous environments, which in plane-wave simulations has a direct
effect on the basis-set size and, as a result, on the cost of the electronic
structure calculation. Moreover, the use of periodic simulation cells are not
well-suited for modeling systems embedded in semi-infinite media, which is
often the case in continuum solvation models. To address these challenges, we
present the implementation of a double-cell formalism, in which the simulation
cell used for the continuum environment is uncoupled from the one used for the
electronic-structure simulation of the quantum-mechanical system. This allows
for a larger simulation cell to be used for the environment, without
significantly increasing computational time. In this work, we show how the
double-cell formalism can be used as an effective PBC correction scheme for
non-periodic and partially periodic systems. The accuracy of the double-cell
formalism is tested using representative examples with different
dimensionalities, both in vacuum and in a continuum dielectric environment.
Fast convergence and good speedups are observed for all the simulation setups,
provided the quantum-mechanical simulation cell is chosen to completely fit
the electronic density of the system.

##  1 Introduction

Materials simulations have seen an incredible growth in recent years, thanks
to the increasing power of computer hardware and simulation software, together
with the development of new computational infrastructures to handle high-
throughput simulations [1, 2, 3, 4, 5]. Unbiased first-principles simulations
based on density functional theory (DFT) represent the work-horse of most
computational approaches to materials design [6], with many research groups
pushing to improve the accuracy of computational predictions and the extend
the scope of first-principles simulations to more challenging systems and
properties.

Traditional DFT simulations of materials focus on bulk periodic crystalline
structures. For this reason, most condensed-matter and materials simulation
packages rely on basis sets composed of periodic functions (e.g. plane-waves,
PW) which uniformly occupy the simulation cell. This approach provides a
systematically improvable basis set, whose accuracy can be tuned with a single
parameter, usually associated with the kinetic energy of the corresponding
plane-wave. PWs also provide easy access to the solution of the core
electrostatic equations involved in first-principles simulations, with the
Poisson equation in vacuum being readily solved in reciprocal space. Fast
Fourier transforms (FFTs)[7] have proven instrumental to the wide-spread
success and fast scaling of PW-based DFT packages, by providing real-to-
reciprocal space mapping, and vice-versa, at a cost proportional to
$N_{pw}\log(N_{pw})$, with $N_{pw}$ being the number of PWs in the simulation.
However, what is advantageous for bulk crystalline materials can be a source
of artifacts when non-periodic or partially periodic systems need to be
characterized. Traditionally, materials interfaces are simulated using the
slab approach, in which a few-atoms thick slab of the material is described in
a 3-dimensional periodically repeated simulation cell. By increasing the
thickness of the slab and by increasing the cell size in the direction
orthogonal to the material interface, the periodic replicas of the slab are
effectively uncoupled from each other and the results are considered to be
representative of the surface of a semi-infinite bulk material. A similar
strategy can be devised for systems that are only periodic in one dimension,
such as nano-wires or nano-tubes, or that are non-periodic. While adjusting
the cell size can provide an easy solution to periodic boundary conditions
(PBC) artifacts, the number of PWs grows linearly with the cell volume,
leading to a polynomial increase in the cost of the corresponding DFT
calculation. For this reason, multiple alternative strategies to correct for
PBC artifacts have been proposed in the past years, including real-space and
reciprocal-space approaches. Most of these strategies are focused on
simulations of two-dimensional and zero-dimensional (isolated) systems, while
one-dimensional correction schemes are usually less wide-spread.

Unfortunately, as the typical scaling of DFT simulations depend varies as
$N^{3}$, where $N$ is the number of electrons in the system, complex and large
systems with more than a few hundreds of electrons are still out of reach for
systematic DFT studies. To overcome this limitation, a variety of hierarchical
approaches [8, 9] and/or divide and conquer[10, 11, 12] strategies have been
developed in the literature. Without sacrificing the important atomistic
characteristics of the relevant part of the process under investigation (the
system), these methods take advantage of different, possibly simplified,
models to handle the more macroscopic or complex part of the process (the
environment), thus allowing a significant reduction of the computational cost.

In particular, implicit solvation models based on continuum embedding media
have proven to provide a reasonable qualitative accuracy, while significantly
reducing the need for statistical sampling of disordered configurations [13,
14, 15, 16]. Indeed, removing the atomistic details of solvent molecules not
only reduces the number of electrons involved in the DFT simulation, but it
allows to avoid the use of molecular dynamics (MD) to sample the
configurations of the molecules of the liquid. Continuum methods in
computational chemistry usually involve more or less empirical definitions of
the different contributions to solvation free energy, which are then
parameterized on experimental databases generated from solubility data.
Methods such as the Polarizable Continuum Model (PCM)[17] of Tomasi and co-
workers have shown a good success in predicting ground state and response
properties of molecules in solution. More recently, similar approaches have
been translated into condensed-matter simulation packages in order to study
wet, possibly electrified, interfaces of materials[8]. When compared to the
more widespread continuum solvation models in the computational chemistry
community, the models implemented in periodic PW-based simulation packages
rely on a smoothly-varying boundary between the quantum-mechanical system and
the continuum environment, instead of a sharp two-dimensional interface.
Definition of the interface based off the electronic density of the solute[18]
or its atomic positions[19] have been developed and show similar accuracy,
when properly parameterized. While computationally more demanding, the use of
smooth interfaces allows continuum models in condensed matter to more
seamlessly introduce non-local corrections [20, 21] to avoid some of the
typical artifacts of these solvation methods.

Beyond the use of continuum solvation to model bulk neutral solutions,
extensions to handle liquid crystals and diluted ionic solutions were
implemented exploiting the integral equation formalism (IEF) of PCM[22]. Still
following the same strategy, PCM was extended to the study of a single two-
dimensional interface between different dielectric media (liquid-air, liquid-
liquid)[23]. Corni et al. extended the use of the polarizable dielectric model
to handle the electrostatic effects of metal surfaces and nanoparticles on
nearby molecular dyes, thus allowing a more quantum-mechanical
characterization of surface-enhanced spectroscopies[24, 25, 26].

Similar to the PCM literature, more complex and heterogeneous environments can
be introduced in condensed-matter simulations. In particular, a significant
focus has been devoted to study electrolyte solutions[23, 27, 28, 29, 30, 31],
so as to unlock the use of these models for simulations of electrified
interfaces. Going beyond neutral solutions and electrolyte distributions,
Campbell et al. [32, 33] exploited a continuum approach to model the charge
reorganization in a semiconductor substrate.

As the definition of the continuum environment in smooth-interface solvation
models follows a function defined everywhere in the simulation cell, it is
relatively straightforward to introduce more complex and heterogeneous
environments, e.g. by introducing multiple continuum media, each with its own
boundaries and physical properties (e.g. dielectric constant or surface
tension). These media can be used to model substrate effects on the properties
of overlaying materials or molecular compounds. Following this idea, Bononi
and collaborators[34, 35] exploited the possibility to introduce flat two-
dimensional dielectric regions in the simulation cell in order to characterize
in a computationally inexpensive way the effect of an ice substrate on the
absorption properties and photo-degradation of molecular dyes.

Despite the flexibility of continuum embedding approaches in condensed matter
environments, PW-based simulations can pose significant challenges to these
class of multiscale methods. The use of periodic simulation cells is
intrinsically unfit to model systems that are embedded in a semi-infinite
medium, i.e. most solvated systems would typically require non-periodic
simulation cells in the directions where the continuum solvent resides.
Moreover, the fact that the simulation cell size is directly related to the
basis-set size $N_{pw}$ makes simulating large heterogeneous environments
computationally more challenging. Extensions of PBC correction schemes have
been developed to account for continuum dielectric media in simulations of
isolated (0D) and slab (2D) systems[36]. In both continuum electrolyte and
semiconductor models, long-range charge reorganizations are handled
implicitly, by relying on analytical solutions to the electrostatic problems
and/or appropriate boundary conditions. However, simulating the effects of
nano-sized environments, such as plasmonic nanoparticles, micelles, liquid
nanodroplets, would require to use simulations cells, and thus basis-sets,
that are beyond reach of standard computational resources.

The straightforward solution to the aforementioned problem is to uncouple the
simulation cell exploited for the continuum environment from the one used for
the DFT simulation of the quantum-mechanical system. The use of a large
simulation cell for the environment would not cause significant limitations to
the simulation time: the time-consuming part of the continuum environment
calculation is related to computing the electrostatic potentials, by solving
modified and generalized forms of the Poisson or Poisson-Boltzmann equations.
Most of the algorithms implemented rely on iterative strategies in which each
iteration involve the solution of a simple Poisson equation, which can be
performed via FFTs in a time almost linear in the cell size. By keeping the
simualtion cell of the DFT calculation small, the largest overhead in
computational time is avoided. While ideal in theory, in practice this double-
cell formalism requires to fully uncouple the two simulation methods and to
carefully design the mapping between the two simulation cells.

While designed to account for large heterogenous environments, the double-cell
formalism can also be used for non-periodic and partially periodic systems as
an effective PBC correction scheme. By increasing the length of the lattice
vectors going along the non-periodic directions in the environment cell,
electrostatics can be fully converged without affecting the basis-set size and
cost of the quantum-mechanical calculation. This general strategy works
equally well for isolated (0D), slab (2D), but also one-dimensional (1D)
systems. In this work, the details of the double-cell methodology are
presented and the approach is tested as a PBC correction scheme for some
representative examples with different dimensionalities, both in vacuum and in
a continuum dielectric environment. Along with this new methodology and in
order to test the accuracy and the speed-up of the double-cell formalism for
PBC corrections, the auxiliary-function countercharge corrections library
(LIBAFCC) of Dabo et al.[37] was coupled to the continuum embedding software,
giving access to an alternative exact correction scheme for all types of
partially periodic systems.

The paper is organized as follows: in Section 2 we review the main details of
the implementation of the double cell formalism, as well as the coupling of
the auxiliary-function countercharge corrections library (LIBAFCC)[37]; in
Section 3 the comptuational details of the DFT calculations as well as how the
double cell calculations are reported. Eventually, in Section 4 we report a
detailed benchmark of the accuracy and performance of the double cell approach
for modelling partially periodic systems, including 0D, 1D, and 2D systems, in
vacuum and in a uniform dielectric medium.

##  2 Methods

###  2.1 Double-Cell Formalism

Figure 1: A schematic view of the double cell formalism with an isolated water
molecule. The green cell is the system cell where all the PWs are introduced
while the blue cell is the environment cell that are used for PBC correction.

The design of a double-cell algorithm is strongly dependent on the numerical
details that are associated with the calling program (e.g. the PW DFT
simulation package), the continuum embedding library, and with their coupling.
In particular, the key aspects that affect the implementation are:

  1. 1.

which scalar fields need to be passed from one cell to the other

  2. 2.

how the scalar fields are stored and parallelized in the two cells

  3. 3.

the way the gridpoints of the two cells are aligned

In the following we will refer to the specific implementation involving
Quantum Espresso as the DFT driver and Environ as the continuum embedding
library. In order to generate the continuum boundary in Environ, information
on the electronic density of the system and/or its atomic positions are
required. The same information is also needed for the calculation of the
electrostatic potential, i.e. the position and valence of the ions as well as
the full electronic density of the system are required inputs for an Environ
calculation. In alternative polarizable dielectric models, such as PCM, the
electrostatic potential of the system in vacuum is passed to the continuum
embedding module, to be used to solve for the polarization charge on the
environment. However, passing the full information of the charge of the solute
allows to have a consistent definition of the electrostatic potential, e.g.
when smeared ionic charges are used to account for core-electrons. In the case
of the double-cell formalism, passing the full charge density of the solute
also allows to seamlessly adapt the calculation to a different simulation
cell. In the presented implementation, the key scalar quantity that needs to
be passed from the system cell into the environment cell is thus the
electronic density of the system.

The results of the continuum embedding calculations are usually reported in
terms of corrections to the corresponding quantities computed in vacuum and
periodic boundary conditions by the host DFT program. In particular,
corrections to the total energy, interatomic forces, and Kohn-Sham potential
are the key outputs of Environ calculations. Of these quantities, only the
latter needs to be mapped from the environment cell into the system cell.

Scalar fields in a 3D simulation cell are discretized in PW simulation
packages in terms of their values on a 3D structured grid, whose spacing is
inversely related to the PW cutoff specified in the simulation input. However,
the values on the 3D grid are stored by Quantum-Espresso as one-dimensional
arrays, by scanning grid-points sequentially along the three cell axis. As
part of the hierarchy of parallelization schemes adopted by Quantum-Espresso,
parallelization on the real-space simulation cell is implemented by dividing
it into slices along the third axis or into sticks along the second and third
axes, depending on the total number of available processors. This allows to
perform operations that only depend on local values of the scalar field in a
fully parallelized way, while still requiring gathering and scattering of the
whole simulation data for real to reciprocal-space transformations (FFTs and
inverse FFTs). Owing to the origin of Environ as an internal plugin of Quantum
Espresso, the same strategies and the same numerical libraries are exploited
by Environ for its internal calculations.

In designing the coupling between system and environment simulation cell, it
is crucial to realize that the quantum-mechanical system lives in a fully
periodic simulation cell: in the DFT calculation the system energy and
properties are not affected by arbitrary translations of the system’s degrees
of freedom. In particular, it is possible for a system to have some of its
atoms or part of its electronic density on opposite sides of the simulation
cell, as the minimum image convention would still give rise to a well-defined
fully connected system. However, when expanding the system cell into the
environment cell, care must be taken to ensure that a fully connected system
is passed to continuum embedding module and that the environment is properly
defined with respect to a meaningful center of the system. Thus, for the
mapping of scalar fields between the two simulation cells, a two step process
was implemented as visualized in Figure 1: first a centering of the system and
second an expansion of the environment cell along the designated directions
(x,y,z).

The first step is accomplished by computing the center of mass of the system
within its original simulation cell, using the minimum image convention to
account for potential splits in neighbouring cells. We then introduce a new
cell with the same cell parameters and the same grid spacing as the system
cell, but with an origin shifted from the Cartesian origin so as to place the
system’s center of mass near its central gridpoint. A one-to-one mapping of
gridpoints, i.e. a change in the integer index of the stored one-dimensional
arrays containing the scalar fields, can be computed on the fly without loss
of information and without affecting the real-space parallelization scheme.

The second step of the double-cell algorithm involves the expansion the system
cell in the different directions in order to generate a larger environment
cell. While arbitrary expansions could in principle be considered, they would
usually result in a environment grid spacing and gridpoint positions that
would not match the initial ones on the system cell. This numerical mismatch
could be easily addressed by interpolation algorithms, which would add an
additional layer in the mapping between the two cells. However, in order to
keep the process as simple and robust as possible, the current implementation
of the double-cell formalism restricts the cell expansion to integer multiples
of the initial cell. An expansion vector
$\mathbf{n}=\left(n_{1},n_{2},n3\right)$ is defined in the input of the
environment module, with $n_{i}$ representing the number of fictitious
replicas added on both sides of the re-centered system cell along its $i$-th
axis. The combination of these two steps creates the desired environment cell,
with lattice vectors $2\mathbf{n}+1$ times as large as the system cell, with
the center of mass of the system located near its central gridpoint.

###  2.2 Auxiliary-Function Correction

When simulating materials using DFT, electron interactions can be hindered by
computational singularities, becoming apparent when looking at partially
periodic systems, arising from the reciprocal-space summation of interaction
contributions. Most simulation software set these divergences to zero,
introducing negligible errors that approach zero as the size of the
computational supercell increases. Real-space electrostatic corrections have
been introduced that can remove these singularities although they can become
costly due to these corrections requiring different interaction potentials per
self-consistent iteration.

Because of this additional computational cost, reciprocal-space countercharge
corrections have been studied in the form of auxiliary-function techniques
which are the basis of the LIBAFCC[37]. These techniques have an added
flexibility with the choice of auxiliary functions which all have different
definitions as well as capabilities in removing PBC errors. The LIBAFCC
library implements exact point-charge reciprocal-space auxiliary functions for
all types of partially periodic systems with the 0D correction giving the same
formulas as the Martyna-Tuckerman[38] (MT) correction scheme.

Following the similarity of the AFCC with the MT approach, we coupled the
library of Li and Dabo to a development version of Environ. The point-charge
correction as computed by the library is mapped in the real-space grid,
centered on the origin. Fast Fourier transforms allow to convert the
correction to reciprocal space and to store it, so that it can be used to
complement the standard vacuum FFT Poisson solver at no extra cost.

##  3 Computational Details

All calculations were performed using Quantum Espresso (QE) v7.1[39, 40]
compiled with the Environ library for continuum embedding effects [18, 39,
41]. While the double-cell implementation has been officially released in
Environ 3.0, for the results reported ion the following we used a local
development version of the library that includes the coupling with the AFCC
library for periodic boundary corrections.

All reported benchmarks, in vacuum and in dielectric environments, involve a
single self-consistent field (SCF) optimization of the electronic density,
starting from a random initial guess. In an effort to keep benchmarks as
consistent as possible, we did not include geometry optimization calculations
in our results. However, convergence of the inter-atomic forces with the newly
implemented algorithms was thoroughly tested, showing an accuracy consistent
with fully periodic simulations and with other PBC correction schemes.

Figure 2: Graphic representation of the two sets of calculations adopted to
investigate convergence of the developed formalism. The double-cell
calculations involve a small system cell (green box) embedded in an expanded
environment cell (blue box). Reference calculations adopt the same expanded
cell for both the system (DFT) and the environment (continuum electrostatic)
calculations.

The 0D and 2D simulations are based off the Environ’s examples distributed
with the library. In particular, the isolated system considered in this study
is an acetamine cation, while a thin Pt (111) slab with an adsorbed carbon
monoxide molecule was used for testing 2D systems. A BN nanoribbon was chosen
as benchmark for the 1D systems. The selected systems are known to display
more pronounced PBC artifacts. However, it is important to stress that the
reported analysis is meant to focus on the consistency of the results with
alternative PBC correction schemes and on the overall computational
performance. For this reason, the physical approximations and the simulation
details that are not relevant for the convergence of the electrostatic energy
have been selected according to the most wide-spread choices in the
literature. All the calculations were performed using the Perdue-Burke-
Ernzerhof generalized gradient approximation (GGA) density functional[42]. For
the 0D and 2D systems, the standard ultra-soft psuedopotentials distributed on
the QE website were adopted. For the 1D system, the pseudopotentials from the
SSSP efficiency library were selected[43, 44]. The Simulations were performed
by only sampling the gamma point sampling of the Brillouin zone.

A dielectric embedding environment was modelled using the self-consistent
interface function of Andreussi et al[18]. A homogeneous static dielectric
permittivity of 100 was selected for the bulk of the environment. The
preconditioned conjugate gradient (PCG) approach of Fisicaro et al [45] was
used for the convergence of the generalized Poisson equation, with a
convergence threshold on the computed electrostatic potential of $5.0e-13$ Ry
[18, 39].

For the acetamine cation (0D system) the wave function cutoff and density
cutoff were chosen to be 30 Ry and 300 Ry, respectively, with a estimated
convergence threshold on the total energy of $5.0e-6$ Ry. A simple cubic
simulation cell was adopted, in order to allow the use of a parabolic
correction to PBCs [46, 47]. For this isolated system, reference simulations
were also performed using the Martyna-Tuckerman reciprocal space correction
[38], both in vacuum and in the dielectric medium.

For the 2D system, an orthorhombic cell was used with the cell parameter in
the z-direction being varied to evaluate the interaction with periodic
replicas. The wave function cutoff and density cutoff were set to 35 Ry and
300 Ry, respectively, with an estimated convergence threshold on the total
energy of $1.0e-6$ Ry. Marzari-Vanderbilt[48] smearing of the band occupations
was adopted with a spread value of 0.03 Ry. A parabolic correction scheme, as
implemented in the Environ library, was adopted for simulations in vacuum and
in the dielectric medium.

Finally, for the 1D system, an ideal boron-nitride (BN) 2D surface, as
obtained from the aflowlib[4] repository, was used to create a 1D ribbon. The
wave function and density cutoffs were set to 60 Ry and 480 Ry, respectively,
with an estimated convergence threshold on the total energy of $1.0e-7$ Ry. A
Gaussian smearing with a spread value of 0.02 Ry was adopted, together with a
local Tomas-Fermi mixing and a mixing parameter of 0.7. Although the parabolic
correction scheme implemented in Environ is only compatible with 0D and 2D
systems, a 2D correction with a large simulation cell in the direction along
the plane of the ribbon was also exploited as reference for the 1D system.

In order to demonstrate the accuracy and performance of the double-cell
formalism, two DFT calculations were done for each unit cell size, as
schematized in Figure 2. The first calculation is performed using double-cell
approach, i.e. having the environment cell defined as an integer expansion of
the system cell in at least one direction. A reference calculation is then
performed by imposing the system and environment cells to be of the same size
and to match the environment cell of the first calculation. As an example, for
a simulation of an isolated (0D) system with a cubic cell of size 10 Bohr, a
double-cell calculation with one replica in every direction,
$\mathbf{n}=\left(1,1,1\right)$, would involve a cubic environment cell with a
30 Bohr side length. The reference calculation would be performed with
identical system and environment cells, both set to a size of 30 Bohr. In the
following section, the energy difference between the double-cell calculation
and the reference are reported together with the speed up in calculation time.

##  4 Results

###  4.1 Isolated Systems (0D)

Figure 3: The acetamine cation used for testing the double-cell formalism in
isolated systems.

As discussed above, one of the uses of the double-cell formalism is to reduce
PBC artifacts in isolated systems, thus overcoming one of the main sources of
artifacts of plane-wave based codes. The system selected for testing the
accuracy and performance of the double-cell strategy is an acetamine cation,
visualized in Figure 3. To simplify the handling of PBC artifacts and to allow
the use of the parabolic point-counter-charge correction, a cubic cell is used
to study this system.

In the top panels of Figure 4 the changes in total energy with respect to the
reference are reported, for both simulations in vacuum and in a dielectric
environment. The considered molecule has overall size of about 8.0 Bohr.
However, from the convergence of the double-cell calculations with respect to
the reference it is clear that system cells smaller than 15 Bohr present
artifacts with respect to the same calculation in the reference cell. As the
electrostatic calculations in the double-cell and the reference are the same,
these small difference must be due to the fact that the double-cell mapping is
affecting the DFT convergence. In fact, inspection of the density of the
molecule reveals that the difference is due to the non-negligible spilling of
the electronic density of the molecule beyond the cell boundaries. While in a
periodic cell this still corresponds to a connected smooth density, the
mapping step of the double-cell formalism breaks the spilling density apart.
This result suggests a general rule for the identification of the minimal
system cell in a double-cell calculation, which needs to be about 10.0 Bohr
larger than the system size.

When looking at the behaviour of the double-cell energy as a function of
system size, it is clear that the double-cell results are fully converged for
cells larger than the identified minimum size. This is consistent with all the
other PBC correction schemes investigated, which also show converged
electrostatic for cell sizes larger than 15 Bohr in vacuum, and 20 Bohr in a
dielectric medium. In the bottom panels of Figure 4 we report the speedup per
SCF cycle of the double-cell formalism when compared with the reference
calculations. This speedup accounts for the relative weight of the FFT-based
electrostatic calculation and the cost of inverting the DFT Hamiltonian in a
single SCF step. The reported results show an average five-fold decrease in
computational time when the system cell is allowed to be smaller than the
environment cell. For the case of dielectric embedding, iterating over the
solution of the Generalized Poisson problem in the dielectric environment
increases the relative costs of the electrostatic calculation with respect to
the SCF step, thus producing a significantly smaller speedup.

While the results in this section aim to provide a comparison between
alternative PBC schemes, it is important to stress that the double-cell
formalism can also be used in combination with any other PBC corrections
scheme reported above: the electrostatic calculation in the expanded
environment cell can be corrected with a real-space term (parabolic
correction) or with any reciprocal-space scheme (Martyna-Tuckerman or AFC)
already implemented in Environ. Those results are not reported in this section
because the simulations are already fully converged without the need for
additional corrections.

Figure 4: The top row plots show the energy difference, in Ry, between an SCF
calculation using the double cell formalism and the corresponding reference
calculation (as visualized in Figure 2). In the middle row, the convergence of
different PBC schemes as a function of DFT system cell size is reported. The
bottom row plots show the time speedup between these simulations. All reported
tests were performed for systems in vacuum (left column) and in a continuum
dielectric medium (right column).

###  4.2 Slab (2D) System

To study 2D systems, systems that should be periodic in two dimensions, we
used a bi-layer of platinum atoms in the (111) orientation and in the presence
a CO molecule in the atop position, as visualized in Figure 5. For this
system, the slab is included in an orthorhombic cell, oriented in the x-y
plane and with the non-periodic axis along third direction. For the double-
cell calculations, the system cell is expanded only in the z direction, while
the system and environment cells are the same in the in-plane directions.
Convergence tests are reported in Figure 6 as a function of the size of the
system cell along the z direction.

The considered system has an extension along the vertical axis of about 11.0
Bohr. When looking at the difference in energy between the double-cell and the
reference calculations, it appears that minimal artifacts are present at the
smallest cell sizes (22-23 Bohr). This suggests that no significant spills in
the electronic density are present for these DFT cells, consistent with our
previous observations for the isolated system. As observed in the 0D case,
double-cell energies are fully converged with respect to the DFT system cell
already for the smallest cells considered. The results are in line with the
convergence

Figure 5: Visualization of the system used for testing the double-cell
algorithm in two dimensions: a Pt (111) surface with two layers of Pt atoms in
the presence of a CO molecule adsorbed in the atop position.

observed for a parabolic (also known as point-counter-charge or dipole
correction) scheme, while AFC90 shows a slower convergence, with exact results
for cells larger than 35 Bohr.

Eventually, in the bottom panels of Figure 6 we report the computational
speedup of using a double cell approach, compared to using the same expanded
cell for both electrostatic and DFT simulations. The relative weight of
electrostatic and inverting the DFT Hamiltonian for this small system provide
a speedup of about 2.5-2.8 in vacuum, that drops to about 1.5 in a dielectric
medium. While the scaling of the FFT-based electrostatic is expected to be
better than the DFT one as a function of the cell size, speedups appear to be
independent of cell size for the considered system. This is probably due to
the relatively small number of electrons in the system and to the other
computational overheads associated with larger FFT grids.

Figure 6: The top row plots show the energy difference, in Ry, between an SCF
calculation using the double cell formalism and the corresponding reference
calculation (as visualized in Figure 2). In the middle row, the convergence of
different PBC schemes as a function of DFT system cell size is reported. The
bottom row plots show the time speedup between these simulations. All reported
tests were performed for systems in vacuum (left column) and in a continuum
dielectric medium (right column).

###  4.3 Wire (1D) System

Figure 7: Visualization of the BN nanoribbon (boron atoms in pink, nitrogen
atoms in blue) investigated to test the double-cell formalism for one-
dimensional systems.

Figure 8: The top row plots show the energy difference, in Ry, between an SCF
calculation using the double cell formalism and the corresponding reference
calculation (as visualized in Figure 2). In the middle row, the convergence of
different PBC schemes as a function of DFT system cell size is reported. The
bottom row plots show the time speedup between these simulations. All reported
tests were performed for systems in vacuum (left column) and in a continuum
dielectric medium (right column).

For the 1D systems, we studied the boron nitride (BN) nanoribbon, as shown in
Figure 7. An orthorhombic cell is adopted, with the ribbon oriented along the
x-axis, making the y and z axes the two non-periodic directions. Also in this
case, only these two directions will be considered in the expansion of the
system cell into the environment cell, while the system and environment cells
have the same size along the x direction. For a sake of simplifying the
discussion, the sizes of the cell axes in the non-periodic directions are kept
identical, allowing to focus our analysis on a single parameter. Convergence
tests are reported in Figure 8 as a function of the such a parameter.

The BN nanoribbon has a width along the y direction of 13 Bohr, with a one-
atom thickness of less than 1 Bohr in the z direction. Convergence of the
double-cell calculations with respect to the reference setup show minor
artifact at the smallest cells considered (10 Bohr). This is consistent with
our understanding of these small-cell artifacts reported for the 0D case and
with the size of the studied system. While a small drift in convergence is
observed at the larger cell sizes, the associated errors are within the
convergence threshold of the performed SCF calculations.

Similarly to the results reported for 2D systems, double-cell energies appear
to be well converged at all the considered system-cell sizes. AFC90 results
show a slightly slower convergence, with some fluctuations for cell sizes
smaller than 20 Bohr. While the parabolic (point-counter-charge) correction is
not implemented for 1D systems, in Figure 8 we report results of using a 2D
scheme applied considering the y axis as the only non-periodic direction of
the system. As the system is very homogeneous along the z direction, we expect
this scheme to capture most of the interaction energy between periodic
replicas of the nanoribbon. Indeed, our simulations show that this parabolic
correction is close to the double-cell results, although some small deviations
are still present for cell sizes as large as 20 Bohr.

For this application, more substantial computational speedups are observed in
using the double-cell formalism. Moreover, a linear increase in the speedup is
observed as a function of the cell size, consistent with the expected
difference in computational scaling between the electrostatic and the DFT
components of the calculations. As for the 0D and 2D cases, the speedup is
significantly reduced for the simulations in a dielectric medium, as the
iterations involved in the solution of the more complex Generalized Poisson
problem increase the weight of the electrostatic calculation.

##  5 Conclusion

We report the development and implementation of a double-cell numerical
algorithm to decouple electrostatic calculations and environment effects from
the underlying DFT simulation. The proposed formalism was tested as a tool to
remove PBC artifacts for partially periodic and non-periodic systems. By
comparing double-cell calculations with a benchmark reference we identified a
source of potential artifacts for very small system cell sizes. Indeed, if the
DFT cell is so small that the electronic density of the system spills across
the cell boundaries, the mapping of the DFT cell into an expanded environment
cell introduces artificial cuts in the electronic density, thus forbidding the
use of the double cell for minimal DFT cell sizes. However, the presented
results should full convergence with cell size for most simulation setups and
system dimensionalities. This fast convergence can be further improved by the
additional use of PBC corrections on the environment cell. While a general
speedup in simulation is observed by allowing the DFT cell to be smaller than
the environment cell, the most impressive result were observed the 1D system,
with speedups of up to 30 times with respect to corresponding single-cell
simulations. The effectiveness and accuracy of the presented approach give us
confidence in the future use of this formalism for simulations with
heterogeneous and nano-scaled continuum environments.

## References

  * [1] Sebastiaan P. Huber, Spyros Zoupanos, Martin Uhrin, Leopold Talirz, Leonid Kahle, Rico Häuselmann, Dominik Gresch, Tiziano Müller, Aliaksandr V. Yakutovich, Casper W. Andersen, and et al.  Aiida 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance.  Scientific Data, 7(1), 2020.
  * [2] Martin Uhrin, Sebastiaan P. Huber, Jusong Yu, Nicola Marzari, and Giovanni Pizzi.  Workflows in aiida: Engineering a high-throughput, event-based engine for robust and modular computational workflows.  Computational Materials Science, 187:110086, 2021.
  * [3] Shyue Ping Ong, William Davidson Richards, Anubhav Jain, Geoffroy Hautier, Michael Kocher, Shreyas Cholia, Dan Gunter, Vincent L. Chevrier, Kristin A. Persson, Gerbrand Ceder, and et al.  Python materials genomics (pymatgen): A robust, open-source python library for materials analysis.  Computational Materials Science, 68:314–319, 2013.
  * [4] Stefano Curtarolo, Wahyu Setyawan, Gus L.W. Hart, Michal Jahnatek, Roman V. Chepulskii, Richard H. Taylor, Shidong Wang, Junkai Xue, Kesong Yang, Ohad Levy, and et al.  Aflow: An automatic framework for high-throughput materials discovery.  Computational Materials Science, 58:218–226, 2012.
  * [5] Claudia Draxl and Matthias Scheffler.  The nomad laboratory: From data sharing to artificial intelligence.  Journal of Physics: Materials, 2(3):036001, 2019.
  * [6] Nicola Marzari, Andrea Ferretti, and Chris Wolverton.  Electronic-structure methods for materials design.  Nature Materials, 20(6):736–749, 2021.
  * [7] David Padua.  Fftw.  Encyclopedia of Parallel Computing, page 671–671, 2011.
  * [8] Oliviero Andreussi and Giuseppe Fisicaro.  Continuum embeddings in condensed-matter simulations.  International Journal of Quantum Chemistry, 119(1), 2018.
  * [9] Hai Lin and Donald G. Truhlar.  Qm/mm: What have we learned, where are we, and where do we go from here?  ChemInform, 38(22), 2007.
  * [10] Wenhui Mi, Xuecheng Shao, Alessandro Genova, Davide Ceresoli, and Michele Pavanello.  Eqe 2.0: Subsystem dft beyond gga functionals.  Computer Physics Communications, 269:108122, 2021.
  * [11] Christoph R. Jacob and Johannes Neugebauer.  Subsystem density-functional theory.  Wiley Interdisciplinary Reviews: Computational Molecular Science, 4(4):325–362, 2014.
  * [12] Xuecheng Shao, Wenhui Mi, and Michele Pavanello.  Gga-level subsystem dft achieves sub-kcal/mol accuracy intermolecular interactions by mimicking nonlocal functionals.  Journal of Chemical Theory and Computation, 17(6):3455–3461, 2021\.
  * [13] Modesto Orozco and F. Javier Luque.  Theoretical methods for the description of the solvent effect in biomolecular systems.  Chemical Reviews, 101(1):203–204, 2001.  PMID: 11749364.
  * [14] Christopher J. Cramer and Donald G. Truhlar.  Implicit solvation models: Equilibria, structure, spectra, and dynamics.  Chemical Reviews, 99(8):009–2665, 1999.
  * [15] Jacopo Tomasi and Maurizio Persico.  Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent.  Chemical Reviews, 94(7):2027–2094, 1994.
  * [16] Jacopo Tomasi, Benedetta Mennucci, and Roberto Cammi.  Quantum mechanical continuum solvation models.  Chemical Reviews, 105(8):2999–3094, 2005.  PMID: 16092826.
  * [17] Jacopo Tomasi and Maurizio Persico.  Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent.  Chemical Reviews, 94(7):2027–2094, 1994.
  * [18] Oliviero Andreussi, Ismaila Dabo, and Nicola Marzari.  Revised self-consistent continuum solvation in electronic-structure calculations.  The Journal of Chemical Physics, 136(6):064102, 2012.
  * [19] Giuseppe Fisicaro, Luigi Genovese, Oliviero Andreussi, Sagarmoy Mandal, Nisanth N. Nair, Nicola Marzari, and Stefan Goedecker.  Soft-sphere continuum solvation in electronic-structure calculations.  Journal of Chemical Theory and Computation, 13(8):3829–3845, 2017\.  PMID: 28628316.
  * [20] Oliviero Andreussi, Nicolas Georg Hörmann, Francesco Nattino, Giuseppe Fisicaro, Stefan Goedecker, and Nicola Marzari.  Solvent-aware interfaces in continuum solvation.  Journal of Chemical Theory and Computation, 15(3):1996–2009, 2019\.
  * [21] Matthew Truscott and Oliviero Andreussi.  Field-aware interfaces in continuum solvation.  The Journal of Physical Chemistry B, 123(16):3513–3524, 2019.
  * [22] E. Cancès, B. Mennucci, and J. Tomasi.  A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics.  The Journal of Chemical Physics, 107(8):3032–3041, 1997.
  * [23] Luca Frediani, Roberto Cammi, Stefano Corni, and Jacopo Tomasi.  A polarizable continuum model for molecules at diffuse interfaces.  The Journal of Chemical Physics, 120(8):3893–3907, 2004.
  * [24] S. Corni and J. Tomasi.  Enhanced response properties of a chromophore physisorbed on a metal particle.  The Journal of Chemical Physics, 114(8):3739–3751, 2001.
  * [25] Oliviero Andreussi, Stefano Corni, Benedetta Mennucci, and Jacopo Tomasi.  Radiative and nonradiative decay rates of a molecule close to a metal particle of complex shape.  The Journal of Chemical Physics, 121(20):10190–10202, 2004.
  * [26] Oliviero Andreussi, Alessandro Biancardi, Stefano Corni, and Benedetta Mennucci.  Plasmon-controlled light-harvesting: Design rules for biohybrid devices via multiscale modeling.  Nano Letters, 13(9):4475–4484, 2013.
  * [27] Ravishankar Sundararaman, William A. Goddard, and Tomas A. Arias.  Grand canonical electronic density-functional theory: Algorithms and applications to electrochemistry.  The Journal of Chemical Physics, 146(11):114104, 2017.
  * [28] Ravishankar Sundararaman and Kathleen Schwarz.  Evaluating continuum solvation models for the electrode-electrolyte interface: Challenges and strategies for improvement.  The Journal of Chemical Physics, 146(8):084111, 2017.
  * [29] Ryosuke Jinnouchi and Alfred B. Anderson.  Electronic structure calculations of liquid-solid interfaces: Combination of density functional theory and modified poisson-boltzmann theory.  Physical Review B, 77(24), 2008.
  * [30] Francesco Nattino, Matthew Truscott, Nicola Marzari, and Oliviero Andreussi.  Continuum models of the electrochemical diffuse layer in electronic-structure calculations.  The Journal of Chemical Physics, 150(4):041722, 2019.
  * [31] Ravishankar Sundararaman, Kendra Letchworth-Weaver, Kathleen A. Schwarz, Deniz Gunceler, Yalcin Ozhabes, and T.A. Arias.  Jdftx: Software for joint density-functional theory.  SoftwareX, 6:278–284, 2017.
  * [32] Quinn Campbell and Ismaila Dabo.  Quantum-continuum calculation of the surface states and electrical response of silicon in solution.  Physical Review B, 95(20), 2017.
  * [33] Quinn Campbell and Ismaila Dabo.  Erratum: Quantum-continuum calculation of the surface states and electrical response of silicon in solution [phys. rev. b 95, 205308 (2017)].  Physical Review B, 96(3), 2017.
  * [34] Ted Hullar, Fernanda C. Bononi, Zekun Chen, Danielle Magadia, Oliver Palmer, Theo Tran, Dario Rocca, Oliviero Andreussi, Davide Donadio, Cort Anastasio, and et al.  Photodecay of guaiacol is faster in ice, and even more rapid on ice, than in aqueous solution.  Environmental Science: Processes &; Impacts, 22(8):1666–1677, 2020.
  * [35] Fernanda C. Bononi, Zekun Chen, Dario Rocca, Oliviero Andreussi, Ted Hullar, Cort Anastasio, and Davide Donadio.  Bathochromic shift in the uv–visible absorption spectra of phenols at ice surfaces: Insights from first-principles calculations.  The Journal of Physical Chemistry A, 124(44):9288–9298, 2020.
  * [36] Oliviero Andreussi and Nicola Marzari.  Electrostatics of solvated systems in periodic boundary conditions.  Physical Review B, 90(24), 2014.
  * [37] Yanli Li and Ismaila Dabo.  Electronic levels and electrical response of periodic molecular structures from plane-wave orbital-dependent calculations.  Phys. Rev. B, 84:155127, Oct 2011.
  * [38] Glenn J. Martyna and Mark E. Tuckerman.  A reciprocal space based method for treating long range interactions in ab initio and force-field-based calculations in clusters.  The Journal of Chemical Physics, 110(6):2810–2821, 1999.
  * [39] P Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, D Ceresoli, M Cococcioni, and et al.  Advanced capabilities for materials modelling with quantum espresso.  Journal of Physics: Condensed Matter, 29(46):465901, 2017.
  * [40] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, and et al.  Quantum espresso: A modular and open-source software project for quantum simulations of materials.  Journal of Physics: Condensed Matter, 21(39):395502, 2009.
  * [41] Edan Bainglass, Medrano Gabriel, Matthew Truscott, Quinn Campbell, Ismaila Dabo, Iurii Timrov, Nicola Marzari, and Oliviero Andreussi.  Environ 3.0: A continuum solvation package [unpublished manuscript].
  * [42] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Physical Review Letters, 77(18):3865–3868, 1996.
  * [43] Gianluca Prandini, Antimo Marrazzo, Ivano E. Castelli, Nicolas Mounet, and Nicola Marzari.  Precision and efficiency in solid-state pseudopotential calculations.  npj Computational Materials, 4(1), 2018.
  * [44] Kurt et al. Lejaeghere.  Reproducibility in density functional theory calculations of solids.  Nature, 351(6280), 2016.
  * [45] G. Fisicaro, L. Genovese, O. Andreussi, N. Marzari, and S. Goedecker.  A generalized poisson and poisson-boltzmann solver for electrostatic environments.  The Journal of Chemical Physics, 144(1):014103, 2016.
  * [46] Oliviero Andreussi and Nicola Marzari.  Electrostatics of solvated systems in periodic boundary conditions.  Physical Review B, 90(24), 2014.
  * [47] Ismaila Dabo, Boris Kozinsky, Nicholas E. Singh-Miller, and Nicola Marzari.  Electrostatics in periodic boundary conditions and real-space corrections.  Physical Review B, 77(11), 2008.
  * [48] Nicola Marzari and David Vanderbilt.  Maximally localized generalized wannier functions for composite energy bands.  Physical Review B, 56(20):12847–12865, 1997.


# A topological photonic ring-resonator for on-chip channel filters

Linpeng Gu1 Qingchen Yuan1 Qiang Zhao2 Yafei Ji1 Ziyu Liu2 Liang Fang1 Xuetao
Gan1∗, Jianlin Zhao1<EMAIL_ADDRESS><EMAIL_ADDRESS>1Key
Laboratory of light-field manipulation and information acquisition, Ministry
of Industry and Information Technology, and Shaanxi Key Laboratory of Optical
Information Technology, School of Science, Northwestern Polytechnical
University, Xi’an 710072, China  2Qian Xuesen laboratory of Space Technology,
China Academy of Space Technology, Beijing 100094, China

A topologically protected ring-resonator formed in valley photonic crystals is
proposed and fabricated on a silicon slab. The unidirectional transmission and
robustness against structure defects of its resonant modes are illustrated.
Coupled with topological waveguides, the topological ring is functioned as
notch and channel-drop filters. The work opens up a new avenue for developing
advanced chip-integrated photonic circuits with attributes of topological
photonics.

##  I Introduction

Recently, topological photonics has attracted remarkable attention due to its
unique properties wang2009observation ; shalaev2019robust ; lu2014topological
; khanikaev2017two ; he2019silicon ; gong2020topological . Photonic analogs of
topological quantum systems can be easily realized via gyro photonic crystal
(PC) with an available promising platform wang2008reflection . A variety of
structures have been designed to emulate the quantum Hall effect, quantum
spin- and valley-Hall effects in electronics by introducing the concept of
pseudo-spin degrees of freedom haldane2008possible ; raghu2008analogs ;
ma2016all ; chen2017valley ; dong2017valley ; wu2015scheme ;
yang2018visualization . Besides, according to the $bulk-boundary$
correspondence, the topologically protected edge waves existing at the
interface between two photonic topological insulators exhibit excellent
robustness, defect immunity, backscatter suppression and unidirectional
transmission properties shalaev2019robust ; he2019silicon ;
khanikaev2013photonic ; chen2014experimental ; ma2015guiding . Based on this,
boundary topological channels have been constructed using two-dimensional
topological photonic crystals (2D-TPC) and potentially employed in a number of
promising scenarios, like high-efficiency non-reciprocal lasing with pure mode
and low threshold bahari2017nonreciprocal ; smirnova2020room ; shao2020high ;
zhang2020low , topological photonic routing he2019silicon ;
gong2020topological ; liu2020z2 ; wang2020valley ; tian2020dispersion ,
topological all-optical logic gates he2019topological , robust delay lines
ma2016all ; hafezi2011robust , and so on. In addition, due to the topological
protection characteristics, light can smoothly transport but sharply turn
along the topological channel. Thus, a compact optical topological cavity with
an arbitrary shape can be designed gong2020topological ; ma2016all ;
bahari2017nonreciprocal , which provides the possibility for constructing a
large-scale robust photonic circuit on a single chip.

On the other hand, silicon photonics has achieved rapid development in recent
years owing to the CMOS-compatibility. Silicon’s high refractive index can
support highly confined optical mode even with a small bend. However, the
backscatter caused by surface or side-wall roughnesses of silicon photonic
devices has always been a trouble for the photonic circuits ji2017ultra ;
li2016backscattering . For example, silicon ring-resonators composed by strip
waveguides are important for wavelength division multiplexing (WDM) based on
their narrowband resonant modes little1997microring ; bogaerts2012silicon .
However, if there is a light backscatter arisen from the surface roughnesses,
the forward incident propagating mode would be coupled by the backward
associated propagating mode, giving rise to the resonance-splitting and the
distorted Lorentzian-shaped spectrum li2016backscattering . As a consequence,
it is desirable to involve the topological robustness into the silicon
photonics to be immune from the structure defect and suppress the light
backscatter.

In this letter, we propose and experimentally demonstrate a topological
photonic ring-resonator fabricated on a silicon chip. It is designed in valley
photonic crystals (VPCs) with a triangular-loop topological edge. The
unidirectional transmission and robustness against the structure defects of
the resonant modes in the ring are illustrated. By coupling the topological
ring with topological waveguides, wavelength-selective notch and channel-drop
filters are illustrated for the potential WDM. This result may extend
applications of topological photonics into silicon photonic integrated
circuits.

##  II Simulation

Figures 1(a) and (b) display schematic diagrams of the proposed topological
ring-resonator designed in 2D-TPCs, which is coupled with a topological
waveguide. The arrows in Fig. 1(a) indicate the direction of the energy flow
of light in the propagating mode. The basic units of the 2D-TPC include two
kinds of honeycomb PC lattices with air-holes in a silicon slab, which are
marked by VPC1 and VPC2 in the inset of Fig. 1(b). Both of the PC lattices
contain three large air-holes and three small air-holes, but they are
carefully arranged with broken inversion symmetry. The corresponding Berry
curvature and topological invariant valley Chern indices $C_{V}$ of the two
VPCs are thus calculated to be opposite, i.e., $C_{V}$<0 for VPC1 and
$C_{V}$>0 for VPC2 he2019silicon . Detailed bulk band for VPC1 and VPC2 has
been calculated and plotted in Fig. 1(c) by using the finite-difference time-
domain (FDTD) method with a lattice constant of $a=410$ nm. A wide TE-like
bandgap appears when two inequivalent air-holes of the unit cell are adjusted
with diameters of $d_{1}=174$ nm and $d_{2}=81$ nm, respectively. According to
the $bulk-boundary$ correspondence, the interface between VPC1 and VPC2 will
exhibit topologically protected robust propagations. In Fig. 1(d), we present
the dispersion relationship of the valley-dependent edge states for the
staggered interface when VPC1 and VPC2 are constructed together. The edge
state is momentum-valley locked with a positive group velocity at K’ valley
and a negative group velocity at K valley. Only one surface mode can propagate
in the forward (K’ valley-locked) or backward (K valley-locked) direction
according to the valley-dependence.

We choose this valley pair of 2D-TPC because it opens a TE-like bandgap at
K(K’) point other than the $\Gamma$ center in the $k$-space. Thus, the
topologically protected edge mode could be well confined in the silicon PC
slab considering the K and K’ valleys are both under the light-cone (the grey
area) as displayed in Fig. 1(c). This property is beneficial for lossless
transmission and for improving the quality of the topological ring-based
devices. Therefore, many efforts have been made to construct the high quality
($Q$) factor ring cavity yamaguchi2019gaas ; ma2019topological ;
noh2020experimental ; zeng2020electrically ; barik2020chiral ;
mehrabad2020chiral  using the VPCs and successfully contribute to the
applications in integrated sources like lasers and emitters yamaguchi2019gaas
; noh2020experimental ; zeng2020electrically ; mehrabad2020chiral . Here, the
topological ring-resonator is formed by folding the topological interface into
a triangular loop, considering the geometric $C_{3v}$ group symmetry. To
facilitate the coupling of the resonant modes in the ring, a topological
waveguide at the interface of VPC1 and VPC2 is designed. Note, because of the
valley-selected unidirectional propagation of the topological edge mode, the
waveguide is connected with the ring directly to satisfy the effective
coupling.

For the light (for example, with a normalized amplitude $E_{0}$) enters from
the input port, as shown in Fig. 1(a), only one topological edge mode
transports in the forward direction, as marked by the yellow arrows. The
backscatter is suppressed owing to the valley ’locked’ feature of the
interface. At the connection point between the topological ring cavity and the
waveguide, the lattice is carefully arranged as shown in Fig. 1(b). The TPC
waveguide is directly connected to the corner of the triangular cavity without
a gap. As a result, the incident waveguide can be topologically connected to
the waveguide of the cavity with effective coupling. Then the incident light
would propagate and split topologically into two routes by following the
valley-dependent rule as presented in Fig. 1(a). One part ($\alpha{E_{0}}$)
passes directly through the connection point. The other part ($\beta{E_{0}}$)
turns over a sharp $120^{\circ}$ bending and propagates forward in a clockwise
(CW) direction in the ring instead of spreading along the $60^{\circ}$
direction. The light is ’locked’ to transport in this specific direction
during its propagation in the triangular loop. After a round, under the same
valley-dependent rule, part of the light in the ring will merge to the through
port and the remaining part continues to participate in a new cycle at the
connection point.

After several cycles, the light in the ring could form stable resonant modes
if its wavelengths meet the resonant condition. The resonant modes are
confined in the ring-resonator, which results in periodic dips in the
transmission spectrum of the output waveguide. The whole process could be
described by adopting the traditional superimposition of a multi-beam
interference. By assuming a lossless balanced, i.e. a 50:50 beam splitting
process at the branch considering the topological property and structural
symmetry, the final transmission of the topological ring can be written as

$T=\frac{0.5+a_{l}^{2}-\sqrt{2}a_{l}{\rm
cos}\delta}{1+0.5a_{l}^{2}-\sqrt{2}a_{l}{\rm cos}\delta}$ (1)

where $a_{l}$ and $\delta=6{\pi}n_{eff}L/\lambda$ are the round trip loss
coefficient and phase shift, respectively. $n_{eff}$ is the effective guiding
index of the propagation mode, $L$ is the side length of the triangular ring
and $\lambda$ is the operation wavelength. When the phase shift $\delta$ is a
multiple of $2\pi$, the light is on resonance and shows a notch in the
transmission spectrum. Note, the resonant mode is topologically protected and
shows an opposite routing direction (CW direction) from the ring-resonator
composed by a strip waveguide, which has the counter-clockwise (CCW) direction
and couples with a bus-waveguide via the evanescent field bogaerts2012silicon
.

Figure 1: (a) & (b) Schematic images of a triangular ring-resonator using the
VPCs on a silicon slab. (c) Bulk photonic bands for VPC1 and VPC2 as they
share the same band structure. (d) Dispersion of the valley-dependent edge
states for the staggered interface between VPC1 and VPC2. (e) Simulated
normalized transmission spectrum of a VPC triangular ring-resonator. Inset
gives the electric field distribution of the ring-resonator at the resonant
wavelength of 1525.2 nm.

Using the FDTD method, we conduct the simulations to analyse the performance
of the topological waveguide-ring structure shown in Fig. 1(a). With the same
lattice parameters in the calculation above, a fine-tuned bandgap is desired
to match with the telecom band, which is compatible with the laser wavelength
in our experiment. In Fig. 1(e), we obtain the simulated normalized
transmission spectrum of the proposed structure. It consists of a series of
transmission notches at the resonant wavelengths of the triangular ring-
resonator. The free spectrum range (FSR) is measured as 10.7 nm in
correspondence to the circumference of the cavity ($3L=105a$). By fitting the
resonant dips with the standard Lorentzian function, the linewidth of the
resonant mode is estimated as 1.26 nm at 1525.2 nm. Therefore, the loaded $Q$
factor is correspondingly calculated as 1,210 with $Q=\lambda/\Delta\lambda$
owing to the large loss induced by the coupling coefficient indicated in Eq.
1. The inset of Fig. 1(c) gives the electric field distribution of the on-
resonance state in the waveguide-ring structure at 1525.2 nm, showing a good
resonant suppression in the waveguide output.

To verify the light propagation behaviour proposed in Fig. 1(a), which is
determined by the topological property of the interface between VPC1 and VPC2,
we further simulate the energy flow of the propagation mode. Figure 2(a) gives
the poynting vectors of the propagating mode in the whole topological
waveguide-ring structure at the resonant wavelength of 1525.2 nm. The incident
light is confined in the interface corresponding to the K’ valley-locked
transmission. A clear way of the energy flow transporting along the CW
direction can be observed throughout the whole ring-resonator. Coherent
interference occurs at the waveguide-ring connection point (region A, shown in
Fig. 2(b)), where no light propagates to the output. The enlarged image in
region B (shown in Fig. 2(c)) gives the vanishing field overlapping between
the two valley states along the VPC channel, which could greatly resist the
defects. In region C, i.e., the corner of the ring triangle, which is zoomed
in Fig. 2(d), the flow shows anti-scatter transmission properties. Light
transports smoothly along a sharp acute angle with little scatterings due to
the topologically protected surroundings of the VPCs he2019silicon .

To prove the robustness of the topological ring-resonator, an air-hole is
intentionally missed at the side-wall of the ring to imitate certain scatters
or the presence of the defect. Figure 2(e) presents the same area of region B
in case with an air-hole defect. The missing-hole defect has been highlighted
with blue dotted curve. Compared with that shown in Fig. 2(c), the flow goes
through the defect with no apparent scatterings along the VPC channel. In Fig.
2(f), more air-holes are missed while the directionality of the flow is still
maintained well. In addition, in Fig. 2(g), we obtain the simulated resonance
dips of the proposed VPC triangular waveguide-ring structures with and without
air-hole defects. For the case with one air-hole defect in region B, as shown
in Fig. 2(e), there is only a slight shift of the resonance dip without
decrease of $Q$ factor. Similar defects in regions A and C are also
considered, showing an ignored impact on the resonance dip as well. For the
case shown in Fig. 2(f), where three air-holes are missed, the resonance shift
become larger. However, when the air-hole defect with the same size of
$d_{2}=81$ nm is dug out from a conventional silicon ring resonator composed
by a strip waveguide with a width of 500 nm, an obvious shift of the resonance
dip is observed, as shown in Fig. 2(h). Moreover, the backscattering induced
by this air-hole defect is enhanced in the ring, resulting in an unwanted
splitting of the resonance li2016backscattering . Therefore, we can conclude
that the intervalley scattering aroused by the defects is naturally suppressed
and this ring-resonator could be designed to against certain perturbations. It
is beneficial for future applications to construct backscattering-free
photonic circuits.

Figure 2: (a) Poynting vector of the resonant mode propagating in the ring-
resonator, showing a CW rotation direction. (b) Zoomed image of region A in
(a) at the waveguide-ring coupling region. (c) Zoomed image of region B in (a)
along the VPC channel. (d) Zoomed image of region C in (a) at the corner of
the triangular ring-resonator showing the anti-scatter transmission. (e) & (f)
Zoomed image of region B along the VPC channel with one and three air-hole
defects, respectively, showing defect immune in comparison with (c). (g)
Simulated resonance dips in the proposed topological ring-resonator with and
without air-hole defects in different regions shown in (a). (h) Simulated
resonance dips in a conventional ring-resonator (with a radius of 4 $\upmu$m)
composed by a strip waveguide with and without an air-hole defect. The inset
is a schematic of the ring-resonator with an air-hole defect.

##  III Experiment Results

To experimentally verify the proposed topological ring-resonator and
demonstrate its function of wavelength-selective filtering, we fabricate it on
a silicon-on-insulator with a 220 nm top silicon layer and 2 $\upmu$m buried
oxide layer. Electron beam lithography combined with inductively coupled
plasma etching are used to accomplish the fabrication. The lattice constant
and diameters of the air-holes of the VPCs are designed appropriately with the
parameters mentioned above. The side length of the ring is about $L=35a$. To
facilitate the light coupling into and out of the topological waveguide-ring,
strip waveguides with tapered width are designed at the two ends of the
topological waveguide. In addition, the strip waveguides have grating couplers
(GCs) at the ends to vertically couple with optical fibers for laser incidence
and transmitted power measurement. Figure 3(a) displays an optical microscope
image of the fabricated device. To clearly characterize the topological
photonic structure, a zoomed scanning electron microscope (SEM) image is shown
in Fig. 3(b) as well. The false colors are added to highlight the regions of
the inequivalent VPC1 and VPC2.

A narrowband tunable laser is coupled into one of the grating couplers, which
transmits along the tapered waveguide and couple into the topological
waveguide and ring. The transmission power of the whole device is monitored
from the other grating coupler. By tuning the laser wavelength, the
transmission spectrum of the device is obtained, as shown in Fig. 3(c). The
transmittance has been normalized by dividing the maximum value obtained from
a single TPC waveguide with the same parameters. Periodic transmission notches
with the Lorentzian lineshapes are observed. The non-flat transmission
background is attributed to the wavelength-dependent coupling efficiency of
the grating couplers. The measured FSR is about 10.7 nm with an extinction
ratio (ER) close to 20 dB, which is consistent with the simulation result
shown in Fig. 1(c). The $Q$ factors are estimated around 1,050, which are
close to the simulation results as well. The obtained transmission dips with a
narrow bandwidth at the resonant wavelengths indicate their abilities of all-
pass notch filter, i.e., light meeting the resonant wavelengths of the ring-
resonator will be filtered out and could not transport from the ’Through’
port.

In addition, similar to the discussion in Fig. 2, we intentionally miss air-
hole defects around the VPC channels in the fabricated devices, as displayed
in Fig. 3(d). Corresponding transmission spectra are plotted in Fig. 3(c) as
well. Though the missed air-hole defects and the inevitable fabrication
errors, the resonant modes of the ring-resonator are not perturbed, showing
barely unchanged resonant wavelengths and $Q$ factors. It indicates the
topologically protected robust transport and defect immunity of the proposed
topological ring-resonator.

Figure 3: (a) Optical microscope image and (b) zoomed SEM image of the
fabricated topological ring-resonator. (c) Measured transmission spectra of
the fabricated devices with and without the missed air-hole. (d) Zoomed SEM
images around the VPC channels with missed air-holes.

The above experiments indicate the possibility of an all-pass notch filter
using the proposed topological ring-resonator, as well as its good defect
immunity property. Next, we illustrate that this structure could also be
employed to achieve a channel-drop filter. As shown in Fig. 4(a), a second
topological channel waveguide is added at another corner of the triangular
ring. To satisfy the topologically protected guiding mode at the VPC
interfaces and effectively couple with the resonant mode in the ring-
resonator, this channel waveguide has an angle of 60 ∘ with respective to the
two sides of the triangular ring and connects with the corner of the ring-
resonator directly. For the on-resonance light incident from the input
topological channel waveguide, it couples into the ring-resonator and
circulates to constructively build up in energy. The spectrum from the direct
output end of the input waveguide (’Through’ port) will exhibit transmission
dips owing to the destructive interference and function as the notch filter,
as discussed above. The resonant mode circulating in the ring will couple into
the second topological channel waveguide, which is similar to the coupling
process between the ring and the input topological waveguide discussed in Fig.
1. Differently, coherent transmission peaks will come from the output of the
second channel waveguide without the interference with the primary input wave.
Consequently, the output end of the second topological channel waveguide could
be considered as a ’Drop’ port for functioning as a channel-drop filter to
select the optical signal at specific wavelengths, which is similar to the
definition in the ring-resonator composed by strip waveguide
little1997microring .

With the parameters of the VPCs employed in Fig. 3, we fabricate the channel-
drop filter, as shown in Figs. 4(a) and (b). The transmission spectra from the
’Through’ port and ’Drop’ port are both measured and plotted in Fig. 4(c).
Similar as the results shown in Fig. 3(c), the ’Through’ port of the filter
presents periodic transmission dips at the resonant wavelengths. The FSR
remains consistent since the ring length is unchanged. However, the linewidths
of the resonant dips are broadened and the ERs are reduced. This could be
attributed to the addition of the second channel waveguide that reduces the
confinement of the resonant mode in the ring via the extra waveguide-coupling
loss. Correspondingly, the spectrum from the ’Drop’ port shows transmission
peaks at the resonant wavelengths with the ER of more than 10 dB, i.e., light
at specific wavelengths are selected via the topological ring-resonator.

Figure 4: (a) Optical microscope image and (b) zoomed SEM image of a channel-
drop filter. (c) Measured transmission spectra from the ’Through’ port and
’Drop’ port, respectively.

##  IV Conclusion

In conclusion, we have successfully demonstrated a topological ring-resonator
fabricated on a silicon chip. It is constructed by folding a topological edge
between two VPCs into a triangular loop, which supports resonant modes with
unidirectional transmission and robustness against structure defects.
Considering the ring-resonator is one of the most important on-chip photonic
elements, the realization of the topological ring-resonator indicates the
topological photonics could provide a new platform immune from material and
structure defects for developing advanced photonic integrated circuits. As an
example, we present notch and channel-drop filters by coupling the topological
ring-resonator with topological waveguides. Although the realization of nano-
scale air-holes in VPCs is currently a challenge as it requires high-
resolution nano-fabrication technique, the topological protection promises the
proposed filters to be applied in more complex integrated devices with higher
performance such as cascaded high-order filters and low crosstalk WDM devices.

## Fundings

the Key Research and Development Program 2017YFA0303800, in part by the
National Natural Science Foundation 91950119, 11634010, 61775183, 61905196, in
part by the Key Research and Development Program in Shaanxi Province of China
2020JZ-10, in part by the Fundamental Research Funds for the Central
Universities 310201911cx032, 3102019JC008, and in part by the seed Foundation
of Innovation and Creation for Graduate Students in Northwestern Polytechnical
University CX2020208.

## Acknowledgment

The authors would thank the Analytical & Testing Center of NPU for the
assistances of device fabrication.

## References

  * (1) Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, “Observation of unidirectional backscattering-immune topological electromagnetic states,” _Nature_ , vol. 461, no. 7265, pp. 772–775, 2009.
  * (2) M. I. Shalaev, W. Walasik, A. Tsukernik, Y. Xu, and N. M. Litchinitser, “Robust topologically protected transport in photonic crystals at telecommunication wavelengths,” _Nature Nanotechnology_ , vol. 14, no. 1, pp. 31–34, 2019.
  * (3) L. Lu, J. D. Joannopoulos, and M. Soljačić, “Topological photonics,” _Nature Photonics_ , vol. 8, no. 11, pp. 821–829, 2014.
  * (4) A. B. Khanikaev and G. Shvets, “Two-dimensional topological photonics,” _Nature Photonics_ , vol. 11, no. 12, pp. 763–773, 2017.
  * (5) X.-T. He, E.-T. Liang, J.-J. Yuan, H.-Y. Qiu, X.-D. Chen, F.-L. Zhao, and J.-W. Dong, “A silicon-on-insulator slab for topological valley transport,” _Nature Communications_ , vol. 10, no. 1, pp. 1–9, 2019.
  * (6) Y. Gong, S. Wong, A. J. Bennett, D. L. Huffaker, and S. S. Oh, “Topological insulator laser using valley-hall photonic crystals,” _ACS Photonics_ , vol. 7, no. 8, pp. 2089–2097, 2020.
  * (7) Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, “Reflection-free one-way edge modes in a gyromagnetic photonic crystal,” _Physical Review Letters_ , vol. 100, no. 1, p. 013905, 2008.
  * (8) F. Haldane and S. Raghu, “Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry,” _Physical Review Letters_ , vol. 100, no. 1, p. 013904, 2008.
  * (9) S. Raghu and F. D. M. Haldane, “Analogs of quantum-hall-effect edge states in photonic crystals,” _Physical Review A_ , vol. 78, no. 3, p. 033834, 2008\.
  * (10) T. Ma and G. Shvets, “All-si valley-hall photonic topological insulator,” _New Journal of Physics_ , vol. 18, no. 2, p. 025012, 2016.
  * (11) X.-D. Chen, F.-L. Zhao, M. Chen, and J.-W. Dong, “Valley-contrasting physics in all-dielectric photonic crystals: Orbital angular momentum and topological propagation,” _Physical Review B_ , vol. 96, no. 2, p. 020202, 2017.
  * (12) J.-W. Dong, X.-D. Chen, H. Zhu, Y. Wang, and X. Zhang, “Valley photonic crystals for control of spin and topology,” _Nature Materials_ , vol. 16, no. 3, pp. 298–302, 2017.
  * (13) L.-H. Wu and X. Hu, “Scheme for achieving a topological photonic crystal by using dielectric material,” _Physical Review Letters_ , vol. 114, no. 22, p. 223901, 2015.
  * (14) Y. Yang, Y. F. Xu, T. Xu, H.-X. Wang, J.-H. Jiang, X. Hu, and Z. H. Hang, “Visualization of a unidirectional electromagnetic waveguide using topological photonic crystals made of dielectric materials,” _Physical Review Letters_ , vol. 120, no. 21, p. 217401, 2018.
  * (15) A. B. Khanikaev, S. H. Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, and G. Shvets, “Photonic topological insulators,” _Nature Materials_ , vol. 12, no. 3, pp. 233–239, 2013.
  * (16) W.-J. Chen, S.-J. Jiang, X.-D. Chen, B. Zhu, L. Zhou, J.-W. Dong, and C. T. Chan, “Experimental realization of photonic topological insulator in a uniaxial metacrystal waveguide,” _Nature Communications_ , vol. 5, no. 1, pp. 1–7, 2014.
  * (17) T. Ma, A. B. Khanikaev, S. H. Mousavi, and G. Shvets, “Guiding electromagnetic waves around sharp corners: topologically protected photonic transport in metawaveguides,” _Physical Review Letters_ , vol. 114, no. 12, p. 127401, 2015.
  * (18) B. Bahari, A. Ndao, F. Vallini, A. El Amili, Y. Fainman, and B. Kanté, “Nonreciprocal lasing in topological cavities of arbitrary geometries,” _Science_ , vol. 358, no. 6363, pp. 636–640, 2017.
  * (19) D. Smirnova, A. Tripathi, S. Kruk, M.-S. Hwang, H.-R. Kim, H.-G. Park, and Y. Kivshar, “Room-temperature lasing from nanophotonic topological cavities,” _Light: Science & Applications_, vol. 9, no. 1, pp. 1–8, 2020\.
  * (20) Z.-K. Shao, H.-Z. Chen, S. Wang, X.-R. Mao, Z.-Q. Yang, S.-L. Wang, X.-X. Wang, X. Hu, and R.-M. Ma, “A high-performance topological bulk laser based on band-inversion-induced reflection,” _Nature Nanotechnology_ , vol. 15, no. 1, pp. 67–72, 2020.
  * (21) W. Zhang, X. Xie, H. Hao, J. Dang, S. Xiao, S. Shi, H. Ni, Z. Niu, C. Wang, K. Jin _et al._ , “Low-threshold topological nanolasers based on the second-order corner state,” _Light: Science & Applications_, vol. 9, no. 1, pp. 1–6, 2020.
  * (22) W. Liu, M. Hwang, Z. Ji, Y. Wang, G. Modi, and R. Agarwal, “Z2 photonic topological insulators in the visible wavelength range for robust nanoscale photonics,” _Nano Letters_ , vol. 20, no. 2, pp. 1329–1335, 2020.
  * (23) M. Wang, W. Zhou, L. Bi, C. Qiu, M. Ke, and Z. Liu, “Valley-locked waveguide transport in acoustic heterostructures,” _Nature Communications_ , vol. 11, no. 1, pp. 1–6, 2020.
  * (24) Z. Tian, C. Shen, J. Li, E. Reit, H. Bachman, J. E. Socolar, S. A. Cummer, and T. J. Huang, “Dispersion tuning and route reconfiguration of acoustic waves in valley topological phononic crystals,” _Nature Communications_ , vol. 11, no. 1, pp. 1–10, 2020.
  * (25) L. He, W. Zhang, and X. Zhang, “Topological all-optical logic gates based on two-dimensional photonic crystals,” _Optics Express_ , vol. 27, no. 18, pp. 25 841–25 859, 2019.
  * (26) M. Hafezi, E. A. Demler, M. D. Lukin, and J. M. Taylor, “Robust optical delay lines with topological protection,” _Nature Physics_ , vol. 7, no. 11, pp. 907–912, 2011.
  * (27) X. Ji, F. A. Barbosa, S. P. Roberts, A. Dutt, J. Cardenas, Y. Okawachi, A. Bryant, A. L. Gaeta, and M. Lipson, “Ultra-low-loss on-chip resonators with sub-milliwatt parametric oscillation threshold,” _Optica_ , vol. 4, no. 6, pp. 619–624, 2017.
  * (28) A. Li, T. Van Vaerenbergh, P. De Heyn, P. Bienstman, and W. Bogaerts, “Backscattering in silicon microring resonators: a quantitative analysis,” _Laser & Photonics Reviews_, vol. 10, no. 3, pp. 420–431, 2016.
  * (29) B. E. Little, S. T. Chu, H. A. Haus, J. Foresi, and J.-P. Laine, “Microring resonator channel dropping filters,” _Journal of Lightwave Technology_ , vol. 15, no. 6, pp. 998–1005, 1997.
  * (30) W. Bogaerts, P. De Heyn, T. Van Vaerenbergh, K. De Vos, S. Kumar Selvaraja, T. Claes, P. Dumon, P. Bienstman, D. Van Thourhout, and R. Baets, “Silicon microring resonators,” _Laser & Photonics Reviews_, vol. 6, no. 1, pp. 47–73, 2012.
  * (31) T. Yamaguchi, Y. Ota, R. Katsumi, K. Watanabe, S. Ishida, A. Osada, Y. Arakawa, and S. Iwamoto, “Gaas valley photonic crystal waveguide with light-emitting inas quantum dots,” _Applied Physics Express_ , vol. 12, no. 6, p. 062005, 2019.
  * (32) J. Ma, X. Xi, and X. Sun, “Topological photonic integrated circuits based on valley kink states,” _Laser & Photonics Reviews_, vol. 13, no. 12, p. 1900087, 2019.
  * (33) W. Noh, H. Nasari, H.-M. Kim, Q. Le-Van, Z. Jia, C.-H. Huang, and B. Kanté, “Experimental demonstration of single-mode topological valley-hall lasing at telecommunication wavelength controlled by the degree of asymmetry,” _Optics Letters_ , vol. 45, no. 15, pp. 4108–4111, 2020.
  * (34) Y. Zeng, U. Chattopadhyay, B. Zhu, B. Qiang, J. Li, Y. Jin, L. Li, A. G. Davies, E. H. Linfield, B. Zhang _et al._ , “Electrically pumped topological laser with valley edge modes,” _Nature_ , vol. 578, no. 7794, pp. 246–250, 2020.
  * (35) S. Barik, A. Karasahin, S. Mittal, E. Waks, and M. Hafezi, “Chiral quantum optics using a topological resonator,” _Physical Review B_ , vol. 101, no. 20, p. 205303, 2020.
  * (36) M. J. Mehrabad, A. P. Foster, R. Dost, E. Clarke, P. K. Patil, A. M. Fox, M. S. Skolnick, and L. R. Wilson, “Chiral topological photonics with an embedded quantum emitter,” _Optica_ , vol. 7, no. 12, pp. 1690–1696, 2020.