# Transverse-electric surface plasmon polaritons in periodically modulated
graphene

Zeeshan Ahmad  Sang Soon Oh  Egor A. Muljarov<EMAIL_ADDRESS>School of Physics and Astronomy, Cardiff University, Cardiff CF24 3AA, United
Kingdom

###### Abstract

Transverse-electric (TE) surface plasmon polaritons are unique eigenmodes of a
homogeneous graphene layer that are tunable with the chemical potential and
temperature. However, as their dispersion curve spectrally lies just below the
light line, they cannot be resonantly excited by an externally incident wave.
Here, we propose a way of exciting the TE modes and tuning their peaks in the
transmission by introducing a one-dimensional graphene grating. Using the
scattering-matrix formalism, we show that periodic modulations of graphene
make the transmission more pronounced, potentially allowing for experimental
observation of the TE modes. Furthermore, we propose the use of turbostratic
graphene to enhance the role of the surface plasmon polaritons in optical
spectra.

††preprint: APS/123-QED

##  I Introduction

Surface plasmon polaritons (SPPs) in graphene have been an interesting area of
theoretical and experimental research especially with the possibility of
supporting SPPs with transverse-electric (TE) polarization in a graphene layer
[1]. More recently, TE SPPs were shown to exist in a homogeneous layer of
graphene for an extended frequency range at non-zero temperature [2] using a
complex frequency approach that models an open system with temporal decay.
Although the TE SPP frequency dispersion is very close to the light line due
to the small magnitude of the conductivity of graphene, proportional to the
fine-structure constant [1], TE SPPs cannot be resonantly excited by an
externally incident light because their dispersion curve lies below the light
line.

Graphene is known for its tunability of optical conductivity by application of
suitable gate voltage that induces sufficiently low but finite chemical
potential [3]. This is because the electronic transitions occur near the
K-point [4] where the electronic dispersion is linear and the density of
states vanishes. Devices like optical modulators [5] and polarizers [6], as
well as absorption enhancement devices [7, 8] benefit from this tunability,
which together with the existence of TE SPPs in graphene, provides exciting
prospects for plasmonic applications [9]. In addition, structures with
periodic graphene open the possibilities of generating topological plasmonic
states [10, 11, 12, 13] when a magnetic field is applied. Periodic plasmonic
structures of graphene [14, 15, 16, 17], and even multilayer stacks of
periodic graphene strips [18, 19, 20, 21, 22] have been already studied. The
effect of stacking graphene-dielectric layers on the properties of transverse-
magnetic (TM) SPPs has also been analyzed for up to ten graphene layers [23].
Although periodic graphene grating structure has been studied for TM SPP modes
[14, 16], it has not been proposed for excitation and measurement of TE SPPs
in graphene, to the best of our knowledge.

Various finite-element and finite-difference time-domain methods have been
employed in many studies of plasmonic structures [11, 24, 25, 26, 27, 28, 29,
30] to investigate SPPs. These fully numerical approaches are usually suitable
for treating electromagnetic systems of relatively small sizes, but suffer
from an increase in computation time for larger samples, since they require
discretization of Maxwell’s equations in real space over entire samples. For
larger samples and extended periodic systems, expanding fields into Fourier
harmonics presents a more suitable approach [31]. One can approximate
structural details by expanding the solutions to Maxwell’s equations in the
Fourier space [32], organizing the obtained solutions as coefficients of Bragg
diffraction orders, which is a more appropriate approach for an infinite
system with translational symmetries presented in this work. Although studies
exist that consider periodic modulation of graphene at the atomic scale [33],
these are suitable for periods of the order of nanometers. For a modulation
period of microns order, which is sufficiently large compared to the
electronic wavelengths, we may approximate the optical conductivity of
graphene as isotropic [34] and consider the periodic modulation on the order
of wavelengths of electromagnetic waves.

In this paper, we theoretically show that the TE SPP modes can be excited by
incident light with an angle close to normal incidence with the help of one-
dimensional periodic modulation in graphene. We demonstrate the excitation of
the TE SPP modes as pronounced dips in the zeroth Bragg order transmission
spectra. We also show that the in-plane wave number of the transmission dips
are tunable with the grating period. To enhance the transmission feature of
the graphene grating, we propose the use of multilayer turbostratic graphene
strips [35, 36].

To obtain the transmission spectra showing features due to the TE SPP modes,
we employ the scattering-matrix formalism [37, 38] that is commonly used for
systems periodically modulated in space. The fact that the width of a single-
or multiple-layer graphene is much smaller than the SPP wavelength allows us
to obtain an explicit analytical expression for the scattering matrix. The SPP
modes, which are the eigen solutions of Maxwell’s equations satisfying
outgoing wave boundary conditions, correspond to the poles of the scattering
matrix for the grating system [38, 39] and manifest themselves as peaks and/or
dips in transmission spectra. We have developed an analytical approximation
for the scattering matrix, eigenmodes, and transmission, providing a proof
that the features observed in optical spectra of periodically modulated
graphene are in fact a manifestation of the TE SPP modes, which are unique to
this kind of material.

Figure 1: Schematic of a graphene grating in a perspective view (a) and in the
$xz$-plane side view (b). An infinitesimally thin graphene is periodically
modulated along the $x$ axis at $z=0$ but homogeneous along the $y$ axis.
Arrows labelled $A_{n},B_{n},C_{n},D_{n}$ schematically represent incoming or
outgoing coefficients of diffracted plane waves with an angle of incidence not
necessarily normal.

##  II Theoretical model

In this section, we first derive, in the spirit of the scattering-matrix
formalism [38], the transmission coefficients for a system containing a thin
layer of a periodically modulated dispersive material surrounded by vacuum.
Generalizations of the presented method to two-dimensional (2D) periodic
modulation of a thin material layer and to systems with a substrate and/or
superstrate are straightforward. The developed general formalism is then
applied to a one-dimensional periodic array of single- and multiple-layer
graphene strips, by using the available analytic expression for the
conductivity of graphene. We also present here a useful approximation of the
transmission at normal incidence, valid for the frequencies near the crossing
point of the folded SPP mode dispersion, which allows us to analyse the
optical response of the periodically modulated system in terms of its SPP
modes.

In the following, we use the model of the permittivity which is described in
the entire space as

$\displaystyle\varepsilon(x,z;\omega)=1+\chi(\omega)\Lambda(x)\delta(z)\,,$
(1)

where $\chi(\omega)$ is the two-dimensional (2D) susceptibility of the
infinitesimal layer, $\Lambda(x)=\Lambda(x+d)$ is a periodic function with
period $d$ describing the spatial modulation of the layer, and $\delta(z)$ is
the Dirac delta function approximating the system as an infinitely thin layer.
An illustration of the system, studied in this work, which is a periodic array
of graphene strips, as well as the geometry of the electromagnetic field in
the TE polarization are shown in Fig. 1(a). The use of the $\delta$ function
is justified by the condition that the graphene thickness is much smaller than
the wavelength of the electromagnetic field.

###  II.1 Scattering matrix, transmission, and the eigenmodes of an inifnitely
thin periodically modulated layer

The scattering matrix, transmission and reflection coefficients, as well as
the secular equation determining the dispersion of SPP modes are obtained in
this section by solving Maxwell’s equations,

$\displaystyle\nabla\times\mathbf{E}$ $\displaystyle=i\omega\mathbf{H}\,,$ (2)
$\displaystyle\nabla\times\mathbf{H}$
$\displaystyle=-i\omega\varepsilon\mathbf{E}\,,$ (3)

in which we assume non-magnetic materials and use the units in which the speed
of light $c=1$. We also assume a harmonic time dependence of the
electromagnetic field in the form of $e^{-i\omega t}$ and use the frequency-
and position-dependent permittivity $\varepsilon(x,z;\omega)$ introduced in
Eq. (1). For zero component of the light wave number in the $y$-direction,
solutions of Maxwell’s equations (2) and (3) split into TE and TM
polarizations. We focus below on the TE polarization. The case of TM
polarization for the same system is addressed in Appendix A.3.

For TE-polarized waves propagating in the $x$-direction [Fig. 1(a)] with the
wave number $q$, the magnetic and electric fields in Cartesian coordinates
take the form

$\displaystyle\mathbf{H}=\begin{pmatrix}H_{x}\\\ 0\\\
H_{z}\end{pmatrix}\,,\quad\mathbf{E}=\begin{pmatrix}0\\\ E_{y}\\\
0\end{pmatrix}\,,$ (4)

respectively. The periodic grating profile $\Lambda(x)$ can be written as a
Fourier series,

$\displaystyle\Lambda(x)=\sum_{n}V_{n}e^{ig_{n}x}\,,$ (5)

where

$g_{n}=\frac{2\pi n}{d}$ (6)

are the reciprocal lattice vectors and $V_{n}$ are Fourier coefficients of the
periodic function $\Lambda(x)$. We can take the general electric-field
solution to Eqs. (2) and (3) in the basis of the same Fourier harmonics,

$\displaystyle
E_{y}(x,z)=\sum_{n}e^{i(q+g_{n})x}\times\begin{cases}C_{n}e^{ik_{n}z}+D_{n}e^{-ik_{n}z}&z>0\,,\\\
A_{n}e^{ik_{n}z}+B_{n}e^{-ik_{n}z}&z<0\,,\\\ \end{cases}$ (7)

where $A_{n}$, $B_{n}$, $C_{n}$, and $D_{n}$ are the field coefficients of the
incoming and outgoing waves of the Bragg order $n$, as depicted in Fig. 1(b).
The normal component of the light wave number for each Bragg order is given by

$k_{n}=\sqrt{\omega^{2}-(q+g_{n})^{2}}\,,$ (8)

with the square root in Eq. (8) chosen in such a way that either ${\rm
Re}\,k_{n}>0$ or ${\rm Im}\,k_{n}>0$, for any real frequency $\omega$ and real
wave number $q$, so the term $C_{n}e^{ik_{n}z}$ in Eq. (7) describes,
respectively, either a purely propagating in the $z$-direction or a purely
evanescent wave. For complex frequencies, the choice of the square root is
determined by a proper analytic continuation of Eq. (8) from the real
$\omega$-axis into the complex frequency plane.

Integrating Eqs. (2) and (3) across the infinitesimal grating layer results in
an infinite set of simultaneous equations for the amplitudes $A_{n}$, $B_{n}$,
$C_{n}$, and $D_{n}$:

$\displaystyle k_{n}(C_{n}\\!-\\!D_{n}\\!-\\!A_{n}\\!+\\!B_{n})$
$\displaystyle=-\omega^{2}\chi(\omega)\sum_{m}V_{n-m}(A_{m}\\!+\\!B_{m})\,,$
$\displaystyle C_{n}\\!+\\!D_{n}$ $\displaystyle=A_{n}\\!+\\!B_{n}\,,$ (9)

see Appendix A for a derivation.

Equations (9) can be written in matrix form as

$\displaystyle\begin{pmatrix}\mathbf{C}\\\
\mathbf{B}\end{pmatrix}=\mathbf{\mathbb{S}}\begin{pmatrix}\mathbf{D}\\\
\mathbf{A}\end{pmatrix}\,,$ (10)

where $\mathbf{\mathbb{S}}$ is the scattering matrix, and $\mathbf{A}$,
$\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ are vectors containing the field
coefficients $A_{n}$, $B_{n}$, $C_{n}$, and $D_{n}$, respectively. The
scattering matrix can be expressed as

$\displaystyle\mathbf{\mathbb{S}}=\begin{bmatrix}-\mathbf{W}^{-1}\mathbf{U}&\mathbf{W}^{-1}\\\
\mathbf{W}^{-1}&-\mathbf{W}^{-1}\mathbf{U}\end{bmatrix}\,,$ (11)

where $\mathbf{W}^{-1}$ is the matrix inverse of
$\mathbf{W}=\mathbf{I}+\mathbf{U}$, $\mathbf{I}$ is the identity matrix, and
$\mathbf{U}$ is the matrix describing the coupling between diffraction orders,
with the matrix elements given by

$\displaystyle U_{nm}=\frac{\omega^{2}\chi(\omega)}{2ik_{n}}V_{n-m}\,,$ (12)

see Appendix A for details.

For each Bragg channel $n$, the transmission coefficient $T_{n}$ is given by
the ratio of the energy of the diffracted wave of order $n$ to the energy of
an incoming wave (corresponding to $n=0$), which can be both evaluated from
the normal component of the Poynting vector [38], which yields

$\displaystyle T_{n}=\left|\frac{k_{n}}{k_{0}}\,C^{2}_{n}\right|\,.$ (13)

Here, it is assumed without loss of generality that

$A_{n}=\delta_{n0}\,,\quad D_{n}=0\,,$ (14)

where $\delta_{nm}$ is the Kronecker symbol, and the expansion coefficients
$C_{n}$ in Eq. (13) are found by solving the linear algebraic equations (10)
[or Eq. (9)] with the boundary conditions Eq. (14). Note that the above
definition of the transmission coefficients, Eq. (13), is valid for so-called
“open” Bragg channels only, for which

$\displaystyle\omega^{2}>(q+g_{n})^{2}\,,$ (15)

equivalent to $k_{n}$ being real, see Eq. (8). For “closed” channels, the
inequality Eq. (15) is not fulfilled, and the values of $T_{n}$ correspond to
the near-field coefficients for these channels.

All the electromagnetic modes of the system, including the SPP modes and their
dispersion relations, i.e. the dependence of the mode frequency on the wave
number, can be found from Eq. (10) by applying the outgoing boundary
conditions, $\mathbf{A}=\mathbf{D}=0$, for all diffraction orders. This is
equivalent to solving a secular equation

$\displaystyle\text{Det[}\mathbf{\mathbb{S}}^{-1}(\omega;q)\text{]}=0\,,$ (16)

or finding the poles of the scattering matrix $\mathbf{\mathbb{S}}(\omega;q)$
in the complex frequency plane for any real wave number $q$.

Figure 2: Electronic dispersion of graphene near the K-point with non-zero
chemical potential $\mu$ and zero temperature $T$.

###  II.2 Graphene conductivity and normalized frequencies and wave numbers

While the general method introduced in Sec. II.1 is suitable for treating any
thin periodically modulated layer, in what follows, we consider a graphene
grating structure depicted in Fig. 1. In this case, the susceptibility is
given by

$\displaystyle\chi(\omega)=\frac{2iN\sigma(\omega)}{\omega},$ (17)

where $\sigma(\omega)$ is the 2D conductivity of a single homogeneous graphene
layer and $N$ is the number of parallel graphene layers stuck together in such
a way that, on the one hand, the total thickness of the system is much smaller
than the SPP wavelength, but on the other hand, all the layers are
electronically decoupled and periodically modulated with the modulation
function $\Lambda(x)$, the same for all layers. For the graphene strips shown
in Fig. 1, we take this periodic function in the following form

$\displaystyle\Lambda(x)=\sum_{n}\Theta\left(\frac{b}{2}-|x-nd|\right)$ (18)

where $\Theta(x)$ is the Heaviside step function, $d$ is the period, and $b$
is the width of the graphene strips.

In the long-wavelength limit, the single-layer frequency-dependent optical
conductivity of graphene $\sigma(\omega)$ is presented in Refs. [1, 40] and is
shown in Fig. 3(a) for zero temperature. In terms of the notations used in
Ref. [2], it is given by

$\displaystyle\sigma(\omega)=2\pi[\sigma_{\text{intra}}(\omega)+\sigma_{\text{inter}}(\omega)]\,,$
(19)

and full expressions for $\sigma_{\text{intra}}$ and $\sigma_{\text{inter}}$,
the intraband and interband parts of the graphene conductivity, are derived in
Ref. [2]. While the interband part, given by the Drude model, describes the
physical mechanism of conductivity typical for normal metals, the interband
part can be understood from the electronic band structure of graphene for a
finite chemical potential $\mu$, as shown in Fig. 2. The interband part
$\sigma_{\text{inter}}$ comes from optical transitions with energy greater
than $2\mu$, contributing to a dip and a step-like feature at
$\text{Re}\,\Omega=2$ which can be seen in Fig. 3(a), respectively, in the
imaginary and real parts of the total conductivity Eq. (19).

Here and below we use for convenience dimensionless frequencies and wave
numbers, all normalized to the chemical potential $\mu$ as

$\Omega=\frac{\hbar\omega}{\mu}\,,\quad Q=\frac{\hbar cq}{\mu}\,,\quad
G_{n}=\frac{\hbar cg_{n}}{\mu}\,,\quad K_{n}=\frac{\hbar ck_{n}}{\mu}\,.$ (20)

Equation (8) then modifies to $K_{n}=\sqrt{\Omega^{2}-(Q+G_{n})^{2}}$, and
similar changes are made in all other equations. The periodicity of the
structure generates the first Brillouin zone with the range
$Q\in[-G_{1}/2,G_{1}/2]$, as depicted in Fig. 3(b). The conductivity
$\sigma(\omega)$ is then called below $\sigma(\Omega)$, which is a function of
the normalized frequency $\Omega$, but it also depends on a normalized inverse
temperature $\mu\beta=\mu/k_{B}T$, where $k_{B}$ is the Boltzmann constant.

Let us note finally that multilayer graphene structures generally have
interlayer coupling which alters the electronic band structure near the
K-point and thus qualitatively changes the optical conductivity spectrum.
However, it is known that turbostratically stacked graphene layers have
negligible interlayer coupling [41], simply factorizing the graphene
conductivity [36]. Such graphene stacks have already been synthesized [42,
35]. We propose to use turbostratic multilayer graphene strips [43] to enhance
features in transmission due to the TE SPP mode, including the increase of its
linewidth and the frequency gap in SPP mode dispersion. In this work, we
obtain transmission spectra and approximate frequencies of the TE modes for
both single-layer ($N=1$) and turbostratic ($N=10$) graphene grating by using
the suitable number of graphene layers $N$ in Eq. (17).

Figure 3: (a) Conductivity of graphene at zero temperature $T=0$, sharing
frequency axis with dispersion in panel (b). The frequency of interband pole
in conductivity is affected by $\mu$ but not $d$. (b) Real part of the
frequency of the TE mode dispersion of homogeneous graphene structure folded
into the first Brillouin zone, shown for $G_{1}=2$, $2.3$ and $2.6$. The size
of the first Brillouin zone, and consequently the frequency of mode crossing
at $Q=0$, is controlled by $\mu d$. The TE mode frequency is close to the
light line due to small value of the fine-structure constant. (c) The period
of grating $d$ required for normalized frequency $\Omega=2$, various chemical
potentials $\mu$ and corresponding frequencies $\omega$ (legend) as a function
of desired incidence angle $\theta=\sin^{-1}(cq/\omega)$ in degrees, using Eq.
(21).

###  II.3 Controlling mode frequency and wave number

The conductivity spectrum of graphene has a key frequency
$\text{Re}\,\Omega=2$ where the interband dip occurs at low temperatures, see
Fig. 3(a). Since the scattering matrix depends on frequency
$\omega=\mu\Omega/\hbar$ and the SPP mode coupling due to the periodic
modulation of the system depends on the grating period $d$, one can control
this coupling and the SPP mode properties by varying the chemical potential
$\mu$ and the period $d$. In fact, by tuning the grating period $d$, the edge
of the first Brillouin zone $G_{1}/2$ where the bands are folded also changes,
leading to variations of the frequency of the TE SPP folded mode crossing at
$Q=0$, as it is clear from Fig. 3(b). Within and above the micrometer range,
the grating period has a negligible effect on the conductivity spectrum, since
the photon wave numbers in this case are much smaller than the electronic wave
numbers [2] — this regime of a negligible spatial dispersion is called above
the long-wavelength limit of the conductivity. One can therefore independently
choose a particular part of the conductivity spectrum shown in Fig. 3(a) and
the desired normalized coupled mode frequency, where the two TE SPP dispersion
lines (which are close to the light lines) cross. These are, for example,
$\Omega=1.667$ and $\Omega=2$ in Fig. 3(a) and $\Omega=2$, $2.3$, and $2.6$ in
Fig. 3(b). For a better description of this tuning process, we define the
normalized interband detuning $\Delta=G_{1}-2$, between the crossing in
dispersion at $\Omega=G_{1}=2\pi\hbar c/(\mu d)$ and the interband dip in the
imaginary part of conductivity at $\Omega=2$ as illustrated in Fig. 3(a). As
an example, for a chemical potential of $\mu$ = 0.2 eV, the zero detuning
$\Delta=0$ would require the graphene stripe width to be in the order of
micrometers. Since the dispersion of the TE mode is close to the light line
[1, 2], we may approximate the TE mode dispersion as the folded light line and
estimate the grating period required to tune the TE mode to a desired
frequency $\Omega$, for any in-plane wave number $Q$. In Fig. 3(c) we show the
dependence of the grating period $d$ required for the TE SPP mode of the
$n=-1$ diffraction order to reach the frequency $\Omega=2$. For general
$\Omega$ and $Q$ this condition is given by

$\displaystyle d=\frac{\hbar c}{\mu}\,\frac{2\pi}{\Omega+Q}\,.$ (21)

Note that this dependence is shown in Fig. 3(c) in terms of different values
of the chemical potential $\mu$ in eV, unnormalized frequency $\omega$ in Hz,
and the angle of incidence $\theta$ measured from the normal to the graphene
plane and defined as $\sin\theta=Q/\Omega$. Clearly $|\theta|\leqslant
30^{\circ}$ in this case, since $\Omega=G_{1}$ and $|Q|\leqslant G_{1}/2$
within the first Brillouin zone.

###  II.4 Role of the TE SPP modes in the transmission: Analytic approximation
at normal incidence

To reveal the effects of periodic modulation on TE SPP modes and the role of
these modes in the transmission, we focus in this section on the normal
incidence and develop approximate analytical expressions for the transmission
around the frequency Re $\Omega=G_{1}$. To do so, we truncate $\mathbf{W}$ in
the scattering matrix [Eq. (11)] to the lowest three diffraction orders, $n=0$
and $n=\pm 1$:

$\displaystyle\mathbf{W}=\begin{pmatrix}\displaystyle
1+\frac{\Omega\tilde{V}_{0}}{K_{-1}}&\displaystyle\frac{\Omega\tilde{V}_{-1}}{K_{-1}}&\displaystyle\frac{\Omega\tilde{V}_{-2}}{K_{-1}}\\\
\displaystyle\frac{\Omega\tilde{V}_{1}}{K_{0}}&\displaystyle
1+\frac{\Omega\tilde{V}_{0}}{K_{0}}&\displaystyle\frac{\Omega\tilde{V}_{-1}}{K_{0}}\\\
\displaystyle\frac{\Omega\tilde{V}_{2}}{K_{1}}&\displaystyle\frac{\Omega\tilde{V}_{1}}{K_{1}}&\displaystyle
1+\frac{\Omega\tilde{V}_{0}}{K_{1}}\end{pmatrix},$ (22)

where $\tilde{V}_{n}=N\sigma(\Omega)V_{n}$ and
$K_{n}\equiv\sqrt{\Omega^{2}-(Q+G_{n})^{2}}$.

The transmission coefficients are then given by the middle column entries of
$\mathbf{W}^{-1}$,

$\displaystyle T_{\pm 1}=\frac{F_{\pm 1}}{\left|\mathbf{W}\right|}\,,\quad
T_{0}=\frac{F_{0}}{\left|\mathbf{W}\right|}\,,$ (23)

where

$\displaystyle F_{\pm 1}$
$\displaystyle=\frac{\Omega^{2}}{K_{1}K_{-1}}\tilde{V}_{\mp 1}\tilde{V}_{\pm
2}-\left(1+\frac{\Omega}{K_{\mp 1}}\tilde{V}_{0}\right)\frac{\Omega}{K_{\pm
1}}\tilde{V}_{\pm 1}\,,$ (24) $\displaystyle F_{0}$
$\displaystyle=\left(1+\frac{\Omega}{K_{-1}}\tilde{V}_{0}\right)\left(1+\frac{\Omega}{K_{1}}\tilde{V}_{0}\right)-\frac{\Omega^{2}}{K_{1}K_{-1}}\tilde{V}_{2}\tilde{V}_{-2}\,,$
(25)

and the determinant of Eq. (22) is given by

$\left|\mathbf{W}\right|=\left(1+\frac{\Omega}{K_{0}}\tilde{V}_{0}\right)F_{0}+\frac{\Omega}{K_{0}}\tilde{V}_{1}F_{-1}+\frac{\Omega}{K_{0}}\tilde{V}_{-1}F_{1}\,.$
(26)

Note that the modulated TE SPP mode frequencies correspond to the zeros of
this determinant.

For $Q=0$, corresponding to normal incidence, the approximate transmission
coefficients obtained in Eq. (23), can be further simplified by using the fact
that $K_{0}=\Omega$ and $K_{\pm 1}=\sqrt{\Omega^{2}-G_{1}^{2}}$, that results
in

$\displaystyle F_{\pm 1}$ $\displaystyle=$
$\displaystyle-\frac{\Omega\tilde{V}_{1}}{K_{1}}\,\mathcal{D}_{-}(\Omega)\,,$
(27) $\displaystyle F_{0}$ $\displaystyle=$
$\displaystyle\mathcal{D}_{-}(\Omega)\,{\mathcal{D}_{+}}(\Omega)\,,$ (28)
$\displaystyle\left|\mathbf{W}\right|$ $\displaystyle=$
$\displaystyle\mathcal{D}_{-}(\Omega)\,\mathcal{B}(\Omega)\,,$ (29)

where

$\mathcal{D}_{\pm}(\Omega)=1+\frac{\Omega}{K_{1}}\left(\tilde{V}_{0}\pm\tilde{V}_{2}\right)\,,$
(30)

and

$\mathcal{B}(\Omega)=1+\tilde{V}_{0}+\frac{\Omega}{K_{1}}\left[\tilde{V}_{0}+\tilde{V}_{2}+\tilde{V}_{0}\left(\tilde{V}_{0}+\tilde{V}_{2}\right)-2\tilde{V}_{1}^{2}\right]\,.$
(31)

Here, without loss of generality, we have made use of the fact that the
Fourier coefficients $V_{n}$ of the grating profile function $\Lambda(x)$ are
real, if $\Lambda(x)$ is an even function of $x$; therefore
$\tilde{V}_{n}=\tilde{V}_{-n}$. Thus, for $Q=0$, the transmission coefficients
in Eq. (23) reduce to

$T_{0}(\Omega)=\frac{{\mathcal{D}_{+}}(\Omega)}{\mathcal{B}(\Omega)}\,,\quad
T_{\pm
1}(\Omega)=-\frac{\Omega}{K_{1}}\frac{\tilde{V}_{1}}{\mathcal{B}(\Omega)}\,.$
(32)

As it is clear from Fig. 3(b), the two folded dispersion lines originating
from the unperturbed TE SPP of the homogeneous layer, corresponding to the
$n=1$ and $n=-1$ Bragg channels and crossing at $\Omega\approx G_{1}$, produce
in a periodically modulated system two perturbed TE SPP modes at each $Q$. The
frequencies of these modes can be found from the zeros of the determinant Eq.
(26), which factorizes at $Q=0$, according to Eq. (29). One mode, having
frequency $\Omega_{d}$, is found from the condition
$\mathcal{D}_{-}(\Omega_{d})=0$, which yields

$\frac{\Omega_{d}}{G_{1}}=\frac{1}{\sqrt{\displaystyle
1-\left(\tilde{V}_{0}-\tilde{V}_{2}\right)^{2}}}\,,$ (33)

and is further referred to as dark mode, as it does not manifest itself in the
normal-incidence transmission. In fact, since $F_{\pm 1}$ and $F_{0}$ are
proportional to $\mathcal{D}_{-}(\Omega)$ there is an exact cancellation of
this factor in $T_{\pm 1}$ and $T_{0}$. The other mode, having frequency
$\Omega_{b}$, which we further call bright mode, is obtained from equating to
zero the other factor in $|\mathbf{W}|$ given by Eq. (29),
$\mathcal{B}(\Omega_{b})=0$, which yields

$\frac{\Omega_{b}}{G_{1}}=\frac{1+\tilde{V}_{0}}{\sqrt{\displaystyle\left(1+\tilde{V}_{0}\right)^{2}-\left[(1+\tilde{V}_{0})(\tilde{V}_{0}+\tilde{V}_{2})-\tilde{V}_{1}^{2}\right]^{2}}}\,,$
(34)

This mode does contribute to the transmission, which in turn can be
approximated for the frequencies $\Omega$ close to $\Omega_{b}$ as

$T_{0}(\Omega)\approx\frac{\mathcal{D}_{-}(\Omega_{b})+(\Omega-\Omega_{b})\mathcal{D}_{-}^{\prime}(\Omega_{b})}{(\Omega-\Omega_{b})\mathcal{B}^{\prime}(\Omega_{b})}\,,$
(35)

and a similar approximate expression can be obtained also for $T_{\pm
1}(\Omega)$.

Figure 4: Transmission spectra ($1-T_{0}$) for zeroth diffraction order
(close-to-normal incidence) when a TE polarised light with parallel wave
number $Q$, is incident on the grating with $b/d=0.3$, $\Delta=0$, for (a)
single-layer graphene with
$\sigma/2\pi=\sigma_{\text{intra}}+\sigma_{\text{inter}}$, (b) single-layer
graphene without interband conductivity $\sigma/2\pi=\sigma_{\text{intra}}$,
and (c) turbostratically stacked graphene with
$\sigma/2\pi=N(\sigma_{\text{intra}}+\sigma_{\text{inter}})$, for $N=10$. The
dashed red line indicates $T_{0}=0$. Spectra for $Q\neq 0$ are offset by
powers of $\times 10$ in increasing order of $Q$. Vertical dash-dotted lines
are at TE SPP resonance frequencies for $Q=0$.

Note that Eqs. (33) and (34) are not explicit expressions for, respectively,
$\Omega_{d}$ and $\Omega_{b}$, but should instead be solved self-consistently,
as the matrix elements $\tilde{V}_{n}$, even though they are small, contain
dependence on frequency, which is present in the conductivity
$\sigma(\Omega)$. However, one can use $\Omega=G_{1}$ as a good starting value
of the argument of $\tilde{V}_{n}(\Omega)$ for a quickly converging self-
consistent iterative solution of Eqs. (33) and (34), which is implemented in
this work and illustrated below in Sec. III.

Note also that in the absence of the periodic modulation,
$\tilde{V}_{n}=\delta_{n,0}N\sigma(\Omega)$, and both modes have exactly the
same frequency, given by $\Omega_{b,d}={G_{1}}/{\sqrt{1-\tilde{V}_{0}^{2}}}$,
which coincides with the homogeneous grahene dispersion [2] folded into the
center of the first Brillouin zone.

##  III Results

We now present transmission spectra of periodically modulated graphene at
$\mu\beta=10^{5}$, with $30\%$ filling factor ($b/d=0.3$), for normal and
close-to-normal incidence of a TE-polarized incoming plane wave. Here, the
grating period $d$ is chosen in such a way that $\Delta=0$ (equivalent to
$G_{1}=2$). For this grating period, the crossing point of the folded light
dispersion at $Q=0$ coincides with the interband dip in the imaginary part of
the conductivity, see Figs. 3(a) and (b). It is therefore expected that the
coupling of the TE SPP modes of a homogeneous graphene, which is caused by its
periodic modulation, is maximized, and thus features in transmission which are
due to the TE SPP modes are enhanced.

Fig. 4 shows the transmission spectra in the zeroth diffraction order,
$1-T_{0}$, for single-layer graphene grating in panel (a), single-layer
graphene grating without interband conductivity $\sigma_{\text{inter}}$ in
panel (b), and turbostratic ($N=10$) graphene grating in panel (c), calculated
with scattering matrix truncated to the lowest $21$ diffraction orders
($-10\leqslant n\leqslant 10$). In all three cases, the first ($n=1$) and
minus first ($n=-1$) diffraction-order light lines are evident as dips in
$1-T_{0}$ that move apart from each other almost linearly with $Q$, as $Q$
increases starting from zero. When the contribution of the interband
conductivity $\sigma_{\text{inter}}$ is removed from the total conductivity
$\sigma$, these dips are clearly seen in Fig. 4(b) in addition to a step-like
feature also moving linearly with $Q$ along the minus-first diffraction-order
light line. This step-like feature, also present in some form in the graphene
transmission [Figs. 4(a) and (c)], is caused by opening the minus-first
diffraction channel with increasing frequency in accordance with Eq. (15).
Such features are common to the optical spectra of periodic systems [44, 45].
In fact, as frequency increases to the point that Eq. (15) is satisfied for
$n=-1$, the transmission of the corresponding ($n=-1$) diffraction order
suddenly increases from zero to a finite value [see Fig. 9(a) in Appendix B]
that results in a corresponding abrupt reduction in $T_{0}$, which is observed
in Fig. 4(b).

The most important features in the spectra of periodically modulated graphene
[Figs. 4(a) and (c)] are also highlighted by their comparison with the normal
metal spectra [Fig. 4(b)]. These are peaks in $1-T_{0}$ which appear also
close to the folded light lines and are the manifestation of the TE SPP modes
in periodically modulated graphene. At the frequencies of these TE SPP peaks,
the energy of the incoming TE beam is diverted back into a reflection channel,
as demonstrated in Appendix B, showing in the same and in a much wider
frequency range the transmission in zeroth and other ($n=\pm 1$) diffraction
orders, as well as the absorption, see Figs. 9 and 12 for $\Delta=0$ (results
for non-zero detuning $\Delta$ are also presented in Appendix B).

In the selected frequency range close to $\Omega=G_{1}$, one would expect to
observe in Figs. 4(a) and (c) two TE SPP modes, separated by a frequency gap,
in accordance with Fig. 3(b). This is a general property of planar photonic-
crystal structures which is usually caused by guided-mode folding and their
hybridization due to the coupling between diffraction orders [38, 46, 47].
However, for $Q=0$, we see in Fig. 4(a) only one peak at around
$\Omega-2\approx-1.76\times 10^{-4}$. The second peak is missing for $Q=0$ and
appears only at nonzero values of $Q$. This is, however, in a full agreement
with our analysis of the transmission, presented in Sec. II.4 in terms of the
bright and dark TE SPP modes. The real part of their frequencies, $\Omega_{b}$
and $\Omega_{d}$, are also shown in Fig. 4(a) by vertical dash-dotted lines,
demonstrating a good agreement with the peak positions. The linewidth of the
dark mode is increasing with $Q$, and both modes are deviating from their
positions at $Q=0$. We evaluate the frequency gap between the two TE SPP modes
to be approximately $1.47\times 10^{-4}$. Note that a similar behaviour has
been observed in transmission spectra of a distributed feedback microcavity
structure [46, 47] in which bright and dark modes originate from the
interaction of the folded guided-mode dispersion lines.

A further analysis of the transmission spectrum at $Q=0$ is presented in Fig.
5 which shows $1-T_{0}$ calculated using the approximation Eq. (35) (red line)
and the full scattering matrix truncated up to $21$ diffraction orders (blue
line), demonstrating a good agreement between the two near
$\Omega\approx\Omega_{b}$ (vertical red dashed line). Away from $\Omega_{b}$,
the approximate transmission Eq. (35), based on the first-order Taylor
expansion, naturally deviates from the exact result. This analysis provides a
proof that the observed peaks in the transmission are in fact a manifestation
of the TE SPP modes predicted in a homogeneous graphene layer [1, 2]. However,
in homogeneous systems, these modes spectrally lie below the light line and
thus cannot be excited resonantly. By periodically modulating the homogeneous
graphene layer, the TE SPP modes can be seen in optical spectra.

Figure 5: Approximation of transmission (red solid line) using Eq. (35)
showing consistency with the TE peak for $Q=0$ and $\Delta=0$. Red and black
dashed lines indicate the real part of the frequencies of the bright and dark
TE SPP modes, obtained using Eq. (34) and Eq. (33), respectively. The
frequency gap between the dashed lines is approximately $1.47\times 10^{-4}$.
Figure 6: Transmission spectra ($1-T_{0}$) for zeroth diffraction order when a
TE polarised light is normally incident ($Q=0$) on grating with $b/d=0.3$, for
$\Delta=0$, shown from zero temperature (red) to finite temperature (green),
for (a) single-layer graphene, and (b) turbostratically stacked graphene. Red
and green dashed dotted lines indicate $T_{0}=0$ for zero temperature and
$T=0.3K$ (single) or $T=3K$ (turbostratic), respectively.

Using turbostratic graphene grating as shown in Fig. 4(c), TE mode peaks
coming from interband conductivity can be further pronounced. Similar to the
single-layer case, the main peak reaches its maximum value of 1. However, the
width of the peaks due to TE SPPs and their positions in frequency are both
affected by the use of turbostratic graphene. Note that the TE mode peak at
$Q=0$ has been shifted further away from the light line, about a factor of
$35$ more compared to its previous position for the single-layer case in Fig.
4(a). The peak is also wider in frequency making TE SPP mode easier to detect.
The peak frequency does not change much for $Q$ close to zero, but eventually
starts to shift to lower frequencies as $Q$ is increased (this can be seen in
Fig. 12(c) of in Appendix B). The dark TE SPP mode peak only appears for
$Q\neq 0$ as seen in Fig. 4(c), similar to the single-layer case. However, it
demonstrates a fine structure, consisting of multiple peaks, which are due to
a stronger coupling between the TE SPP modes and higher diffraction orders
more strongly contributing to the scattering matrix.

The effect of change in temperature on the transmission spectrum is presented
in Fig. 6 for $Q=0$. Temperature appears in $\mathbf{\mathbb{S}}$ through the
conductivity of graphene, where the increase in temperature has the smearing
effect on the interband dip in the imaginary part of graphene conductivity
[2]. Figure 6 shows that the TE mode peak is reaching its maximum height of
$1$ at zero temperature in both the single-layer and turbostratic case. For
the single-layer case, increasing temperature smears the peak out, with the
peak maximum reduced. For the turbostratic case, the peak is still reaching
the maximum of $1$ at $3K$ with minimal smearing. This insensitivity of
transmission to temperature suggests further in favour of turbostratic over
single-layer graphene for applications using TE SPPs that aim to be robust
against temperature change.

##  IV Conclusions

We have shown that the complex-frequency transverse-electric surface plasmon
polariton (TE SPP) modes in a graphene grating can be excited by external
electro-magnetic waves with close-to-normal incidence. Specifically, a
significant reduction in the zeroth diffraction-order transmission is observed
near the interband transition frequency of graphene. This reduction can be
explained as a reflection of the energy of the incoming beam due to the
resonant coupling to the TE SPP modes. In addition, we have shown that the
frequency and the in-plane wave number of TE SPP modes seen in transmission
can be tuned by controlling the graphene chemical potential and grating
period. Specifically, we have demonstrated that for zero detuning, that is the
grating period is chosen so that the TE SPP mode energy matches the interband
threshold of the graphene conductivity, the features in the optical spectra
due to the TE SPPs are significantly enhanced. Furthermore, by using
turbostratic graphene layers with negligible interlayer interactions, we can
enhance these features further, making them experimentally more accessible
than in the single-layer graphene.

## Acknowledgments

Z.A acknowledgments financial support of EPSRC under the DTP scheme. S.S.O.
acknowledges the support of Sêr Cymru II Rising Star Fellowship [80762-CU 145
(East)], which is partly funded by the European Regional Development Fund
(ERDF) via the Welsh Government.

##  Appendix A Scattering-matrix approach

###  A.1 Solving Maxwell’s equations for an infinitely thin periodic layer

We start from Maxwell’s equations having the form

$\displaystyle\nabla\times\mathbf{E}$
$\displaystyle=-\partial_{t}\mathbf{H}\,,$ (36)
$\displaystyle\nabla\times\mathbf{H}$
$\displaystyle=\partial_{t}\mathbf{D}\,,$ (37)

where we assume non-magnetic materials, i.e. the permeability $\mu=1$
everywhere in space, and use the units in which the speed of light $c=1$.
Assuming a harmonic dependence on time of the electro-magnetic field in the
form of $e^{-i\omega t}$, we obtain

$\displaystyle\nabla\times\mathbf{E}$ $\displaystyle=i\omega\mathbf{H}\,,$
(38) $\displaystyle\nabla\times\mathbf{H}$
$\displaystyle=-i\omega\varepsilon(\omega)\mathbf{E}\,.$ (39)

For the TE polarization, the magnetic and electric fields can be expressed in
terms of their Cartesian components as:

$\displaystyle\mathbf{H}=\begin{pmatrix}H_{x}\\\ 0\\\
H_{z}\end{pmatrix}\,,\quad\mathbf{E}=\begin{pmatrix}0\\\ E_{y}\\\
0\end{pmatrix}\,.$ (40)

Maxwell’s equations then take the form

$\displaystyle\begin{pmatrix}-\partial_{z}E_{y}\\\ 0\\\
\partial_{x}E_{y}\end{pmatrix}$ $\displaystyle=+i\omega\begin{pmatrix}H_{x}\\\
0\\\ H_{z}\end{pmatrix}\,,$ (41)
$\displaystyle\begin{pmatrix}\partial_{y}H_{z}\\\
\partial_{z}H_{x}-\partial_{x}H_{z}\\\ -\partial_{y}H_{x}\end{pmatrix}$
$\displaystyle=-i\omega\varepsilon\begin{pmatrix}0\\\ E_{y}\\\
0\end{pmatrix}\,.$ (42)

These can be rewritten as

$\displaystyle-\partial_{z}E_{y}$ $\displaystyle=i\omega H_{x}\,,$ (43)
$\displaystyle\partial_{x}E_{y}$ $\displaystyle=i\omega H_{z}\,,$ (44)
$\displaystyle\partial_{z}H_{x}-\partial_{x}H_{z}$
$\displaystyle=-i\omega\varepsilon E_{y}.$ (45)

As defined in Eq. (1), the permittivity in the whole space is given by
$\varepsilon(x,z;\omega)=1+\chi(\omega)\Lambda(x)\delta(z)$, where
$\chi(\omega)$ is the 2D susceptibility of a thin (graphene) layer and
$\Lambda(x)$ is a periodic function describing its spatial modulation, which
can be represented by its Fourier series as

$\displaystyle\Lambda(x)=\sum_{n}V_{n}e^{ig_{n}x}\,.$ (46)

Integrating Eqs. (43) and (45) over $z$ across the point $z=0$, we obtain

$\displaystyle-\left.E_{y}\right|_{z=0^{+}}+\left.E_{y}\right|_{z=0^{-}}$
$\displaystyle=0\,,$ (47)
$\displaystyle\left.H_{x}\right|_{z=0^{+}}-\left.H_{x}\right|_{z=0^{-}}$
$\displaystyle=-i\omega\chi(\omega)\Lambda(x)\left.E_{y}\right|_{z=0}\,,$ (48)

where $0^{+(-)}$ is a positive (negative) infinitesimal. For the TE
polarization, it is easier to eliminate $H_{x}$ and write equations in terms
of the $E_{y}$ component having the following general form

$\displaystyle
E_{y}(x,z)=\sum_{n}e^{i(q+g_{n})x}\times\begin{cases}C_{n}e^{ik_{n}z}+D_{n}e^{-ik_{n}z}&z>0\,,\\\
A_{n}e^{ik_{n}z}+B_{n}e^{-ik_{n}z}&z<0\,.\\\ \end{cases}$ (49)

in which $g_{n}$ and $k_{n}$ are defined, respectively, by Eqs. (6) and (8),
and $q$ is the wave number along $x$.

Substituting Eq. (43) into Eq. (48), we obtain

$\displaystyle\left.\partial_{z}E_{y}\right|_{z=0^{+}}$
$\displaystyle-\left.\partial_{z}E_{y}\right|_{z=0^{-}}=-\omega^{2}\chi(\omega)\Lambda(x)\left.E_{y}\right|_{z=0},$
(50)

which after using the series Eqs. (46) and (49) becomes

$\displaystyle\sum_{n}$ $\displaystyle
ik_{n}(C_{n}-D_{n}-A_{n}+B_{n})e^{i(q+g_{n})x}$
$\displaystyle=-\omega^{2}\chi(\omega)\sum_{mn}e^{iqx}e^{ig_{m+n}x}V_{m}(A_{n}+B_{n}).$
(51)

By shifting the indices of summation $n\xrightarrow{}n-m$ in the right-hand
side of Eq. (51), we find

$\displaystyle\sum_{n}$ $\displaystyle
ik_{n}(C_{n}-D_{n}-A_{n}+B_{n})e^{i(q+g_{n})x}$
$\displaystyle=-\omega^{2}\chi(\omega)\sum_{m,n}e^{iqx}e^{ig_{n}x}V_{m}(A_{n-m}+B_{n-m})\,,$
(52)

and then equating the coefficients at $e^{i(q+g_{n})x}$ in Eq. (52), we obtain

$\displaystyle ik_{n}$ $\displaystyle(C_{n}-D_{n}-A_{n}+B_{n})$
$\displaystyle=-\omega^{2}\chi(\omega)\sum_{m}V_{m}(A_{n-m}+B_{n-m})\,.$ (53)

Finally, shifting the indices again in the summation in Eq. (53), by
substituting $m\xrightarrow{}-m+n$, we arrive at

$ik_{n}(C_{n}\\!-\\!D_{n}\\!-\\!A_{n}\\!+\\!B_{n})\\!=\\!-\omega^{2}\chi(\omega)\sum_{m}V_{n-m}(A_{m}\\!+\\!B_{m})\,.$
(54)

In addition to this, we find from the continuity of the electric field, Eq.
(47), that

$\displaystyle C_{n}+D_{n}=A_{n}+B_{n}\,.$ (55)

An infinite set of simultaneous equations given by Eqs. (54) and (55) presents
a general solution of Maxwell’s equation in TE polarization for the infinitely
thin periodically modulated layer.

Figure 7: Transmission spectra for single layer graphene ($N=1$), for the
minus-first diffraction order (a),(b), zeroth order (c),(d), and first order
(e),(f), for the period $d$ adjusted so that $\Delta=-1$ ($G_{1}=1$). Panels
(b) and (d) zoom in the spectral features near $\Omega=G_{1}$ in (a) and (c),
respectively, for $Q=0$. Panel (g) is absorption of open channels calculated
using Eq. (64), with waterfall increment of $+0.02$. $A(\Omega<2)=0$ for zero
temperature. The gray solid lines correspond to the forbidden frequency range
where Eq. (15) is not satisfied, and the curve has the meaning of the near-
field amplitude. Figure 8: As Fig. 7 but for $\Delta=-2+1.667$
($G_{1}=1.667$), illustrating suppression of near-field and transmission
amplitudes. Figure 9: As Fig. 7 but for $\Delta=0$ ($G_{1}=2$), illustrating
the amplification of the near-field and transmission amplitudes. Figure 10:
Transmission spectra for turbostratic ($N=10$) graphene, for the minus first
diffraction order (a),(b), zeroth order (c),(d), and first order (e),(f) for
detuning $\Delta=-1$ ($G_{1}=1$). Panels (b) and (d) zoom in the spectral
features near $\Omega=G_{1}$ in (a) and (c), respectively, for $Q=0$.
Similarly, panel (f) is zoom in of panel (e) at $Q$ close to but not zero.
Gray regions and gold curves indicate forbidden frequency regions. Panel (g)
is absorption spectrum obtained using Eq. (64), with waterfall increment of
$+0.2$. Figure 11: As Fig. 10 but for $\Delta=-2+1.667$ ($G_{1}=1.667$).
Figure 12: As Fig. 10 but for $\Delta=0$ ($G_{1}=2$). Figure 13: As Fig. 10
but for $\Delta=0.5$ ($G_{1}=2.5$).

###  A.2 Deriving the scattering matrix

To build a scattering matrix, let us rearrange Eqs. (54) and (55) in the
following way

$\displaystyle C_{n}+B_{n}$
$\displaystyle+\sum_{m}\frac{\omega^{2}\chi(\omega)}{ik_{n}}V_{n-m}B_{m}$
$\displaystyle=A_{n}+D_{n}-\sum_{m}\frac{\omega^{2}\chi(\omega)}{ik_{n}}V_{n-m}A_{m}$
(56) $\displaystyle C_{n}-B_{n}$ $\displaystyle=A_{n}-D_{n}$ (57)

Thus, in terms of $2\times 2$ matrices of diffraction-order blocks, the above
equations can be written as

$\displaystyle\begin{bmatrix}\mathbf{I}&\mathbf{I}+2\mathbf{U}\\\
\mathbf{I}&-\mathbf{I}\end{bmatrix}\begin{bmatrix}\mathbf{C}\\\
\mathbf{B}\end{bmatrix}=\begin{bmatrix}\mathbf{I}&\mathbf{I}-2\mathbf{U}\\\
-\mathbf{I}&\mathbf{I}\end{bmatrix}\begin{bmatrix}\mathbf{D}\\\
\mathbf{A}\end{bmatrix}\,,$ (58)

where vectors $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ have
components $A_{n}$, $B_{n}$, $C_{n}$, and $D_{n}$, respectively, matrix
$\mathbf{U}$ has the matrix elements

$\displaystyle U_{nm}=\frac{\omega^{2}\chi(\omega)}{2ik_{n}}V_{n-m}\,,$ (59)

and $\mathbf{I}$ is the identity matrix.

The scattering matrix, relating the amplitudes of incoming and outgoing waves,
is then given by

$\displaystyle\mathbf{\mathbb{S}}=\begin{bmatrix}\mathbf{I}&\mathbf{I}+2\mathbf{U}\\\
\mathbf{I}&-\mathbf{I}\end{bmatrix}^{-1}\begin{bmatrix}\mathbf{I}&\mathbf{I}-2\mathbf{U}\\\
-\mathbf{I}&\mathbf{I}\end{bmatrix}.$ (60)

To further simplify it, we can write the inverse of the square matrix in Eq.
(60) in terms of its blocks:

$\displaystyle\mathbf{\mathbb{S}}=\frac{1}{2}\begin{bmatrix}\mathbf{W}^{-1}&\mathbf{W}^{-1}(\mathbf{I}+2\mathbf{U})\\\
\mathbf{W}^{-1}&-\mathbf{W}^{-1}\end{bmatrix}\begin{bmatrix}\mathbf{I}&\mathbf{I}-2\mathbf{U}\\\
-\mathbf{I}&\mathbf{I}\end{bmatrix},$ (61)

where $\mathbf{W}=\mathbf{I}+\mathbf{U}$. Finally, performing the matrix
multiplication, we obtain

$\displaystyle\mathbf{\mathbb{S}}=\begin{bmatrix}-\mathbf{W}^{-1}\mathbf{U}&\mathbf{W}^{-1}\\\
\mathbf{W}^{-1}&-\mathbf{W}^{-1}\mathbf{U}\end{bmatrix}.$ (62)

###  A.3 TM polarization

We note that TM polarization can be treated in exactly the same manner as TE
polarization solved above, leading to the same result for the scattering
matrix, provided that matrix $\mathbf{U}$ is replaced with
$\mathbf{U}^{\text{TM}}$ having the following matrix elements:

$\displaystyle U_{nm}^{\text{TM}}=-\frac{1}{2}ik_{n}\chi(\omega)V_{n-m}\,.$
(63)

Also, in TM polarization, the vectors $\mathbf{A}$, $\mathbf{B}$,
$\mathbf{C}$, and $\mathbf{D}$ have the same meaning of the expansion
coefficients in Eq. (49), but for $H_{y}(x,z)$ instead, which in this case is
the only non-zero component of the magnetic field.

##  Appendix B Transmission and absorption spectra

Here we present the transmission spectra, calculated using Eq. (13), for
grating with the filling factor $b/d=0.3$ for the first Brillouin zone
$Q\in[0,G_{1}/2]$ and for different grating periods such that $\Delta=-1$,
$-2+1.667$, and 0. We also show the absorption defined as

$A(\Omega)=1-\sum_{n\in\\{\Omega>|Q+G_{n}|\\}}\left[T_{n}(\Omega)+R_{n}(\Omega)\right]\,,$
(64)

where

$R_{n}(\Omega)=\left|\frac{K_{n}}{K_{0}}B_{n}^{2}\right|$ (65)

is the reflection coefficient of the $n$-th Bragg channel. Note that the
summation in Eq. (64) is performed over open Bragg channels only, satisfying
the condition Eq. (15), as indicated in the formula.

###  B.1 Single-layer case

Figures 7, 8, and 9 show transmission and absorption spectra for single-layer
($N=1$) graphene grating for $\Delta=-1$, $-2+1.667$, and $0$, respectively,
corresponding to $G_{1}=1$, $1.667$, and $2$.

The SPP modes manifest themselves as peaks in the transmission spectra in
Figs. 7, 8, and 9 in panels (a), (c), and (e), and as dips in panel (g), all
moving almost linearly with $Q$. The absorption $A$ in panel (g) is zero for
$\Omega<2$ and is about $10^{-2}$ otherwise, which correlates with the
interband absorption having a sharp frequency cut-off at $\Omega=2$ for zero
temperature. In the zeroth-order transmission $1-T_{0}$, the TE SPP mode peaks
are lost due to the stronger interband absorption for frequencies above
$\Omega>2$, instead appearing as dips in the absorption spectrum $A$, for such
frequencies. Since the mode features are sharp in frequency, we show zoom in
of the TE mode transmission at $Q=0$ for minus-first [panels (b)] and zeroth
[panels (d)] diffraction orders. First order transmission for small but finite
$Q$ is shown in panels (f).

For $\Delta=-1$ (Fig. 7), the TE mode dip (not peak) for $Q=0$ appears at
frequency very close to $G_{1}$, see panel (b). When $Q$ is slightly increased
[panel (f)], the TE mode dip moves away from $G_{1}$ increasing in frequency
within the open channel region (not gray).

For $\Delta=-2+1.667$ (Fig. 8), the TE mode peak can be also seen moving to
higher frequencies for $Q$ slightly increasing from $Q=0$ [panel (b)] to
$Q\neq 0$ [panel (f)]. The peak width is very small, about $10^{-10}$ in this
case due to the vanishing graphene conductivity near this frequency of
$\Omega\approx 1.667$ [2].

For $\Delta=0$ (Fig. 9), the TE mode peak is clearly enhanced due to the dip
in the imaginary part of the conductivity. The TE SPP peak in the coefficient
of the evanescent near field (gray area) moves to lower frequencies for
increasing $Q$ from panel (b) to (f), in agreement with the observations in
Fig. 4(a).

###  B.2 Turbostratic case

Figures 10, 11, 12, and 13 show transmission and absorption spectra for
turbostratic ($N=10$) graphene grating for $\Delta=-1$, $-2+1.667$, $0$, and
$0.5$ respectively, corresponding to $G_{1}=1$, $1.667$, $2$ and $2.5$.
Features related to the mode peaks and interband peaks are qualitatively the
same here as in the single layer case in Sec. B.1 which can be seen by
comparing panels (a), (c), (e), and (g) of all transmission spectra in this
Appendix.

The factor of $N=10$ multiplying the conductivity increases the absorption by
the same factor to the order of $10^{-1}$. The zoom in panels (b), (d), and
(f) show effects of this factor in amplitude of peaks and the frequency of
peaks and dips moving further away from $G_{1}$ by a factor of $10^{2}$,
compared to the single-layer case. Importantly, the zoom in frequency near
$\Omega=2$ in Figs. 12(b), (d), and (f) shows evidence of enhancement in first
and minus-first diffraction orders and a suppression to zero of the
transmission ($1-T_{0}=1$) in the zeroth diffraction order.

## References

  * Mikhailov and Ziegler [2007] S. A. Mikhailov and K. Ziegler, New electromagnetic mode in graphene, Physical Review Letters 99, 016803 (2007).
  * Ahmad  _et al._ [2021] Z. Ahmad, E. A. Muljarov, and S. S. Oh, Extended frequency range of transverse-electric surface plasmon polaritons in graphene, Physical Review B 104, 085426 (2021).
  * Lee  _et al._ [2012] S. H. Lee, M. Choi, T.-T. Kim, S. Lee, M. Liu, X. Yin, H. K. Choi, S. S. Lee, C.-G. Choi, S.-Y. Choi, X. Zhang, and B. Min, Switching terahertz waves with gate-controlled active graphene metamaterials, Nature Materials 11, 936 (2012).
  * Zhan  _et al._ [2012] D. Zhan, J. Yan, L. Lai, Z. Ni, L. Liu, and Z. Shen, Engineering the electronic structure of graphene, Advanced Materials 24, 4055 (2012).
  * Bahadori-Haghighi  _et al._ [2020] S. Bahadori-Haghighi, R. Ghayour, and A. Zarifkar, Tunable graphene–dielectric metasurfaces for terahertz all-optical modulation, Journal of Applied Physics 128, 044506 (2020).
  * Bao  _et al._ [2011] Q. Bao, H. Zhang, B. Wang, Z. Ni, C. H. Y. X. Lim, Y. Wang, D. Y. Tang, and K. P. Loh, Broadband graphene polarizer, Nature Photonics 5, 411 (2011).
  * Maji and Das [2018] P. S. Maji and R. Das, Absorption enhancement in monolayer graphene using tamm plasmon polaritons, OSA Continuum 1, 392 (2018).
  * Guo  _et al._ [2019] Y. Guo, S. Wang, Y. Zhou, C. Chen, J. Zhu, R. Wang, and Y. Cai, Broadband absorption enhancement of graphene in the ultraviolet range based on metal-dielectric-metal configuration, Journal of Applied Physics 126, 213103 (2019).
  * He  _et al._ [2013] X. Y. He, J. Tao, and B. Meng, Analysis of graphene te surface plasmons in the terahertz regime, Nanotechnology 24, 345203 (2013).
  * Kuzmin  _et al._ [2018] D. A. Kuzmin, I. V. Bychkov, V. G. Shavrov, and V. V. Temnov, Plasmonics of magnetic and topological graphene-based nanostructures, Nanophotonics 7, 597 (2018).
  * Jin  _et al._ [2017] D. Jin, T. Christensen, M. Soljačić, N. X. Fang, L. Lu, and X. Zhang, Infrared topological plasmons in graphene, Physical Review Letters 118, 245301 (2017).
  * Slipchenko  _et al._ [2021] T. M. Slipchenko, J.-M. Poumirol, A. B. Kuzmenko, A. Y. Nikitin, and L. Martín-Moreno, Interband plasmon polaritons in magnetized charge-neutral graphene, Communications Physics 4, 110 (2021).
  * Taillefumier  _et al._ [2011] M. Taillefumier, V. K. Dugaev, B. Canals, C. Lacroix, and P. Bruno, Graphene in a periodically alternating magnetic field: An unusual quantization of the anomalous hall effect, Physical Review B 84, 085427 (2011).
  * Gonçalves  _et al._ [2016] P. A. D. Gonçalves, E. J. C. Dias, Y. V. Bludov, and N. M. R. Peres, Modeling the excitation of graphene plasmons in periodic grids of graphene ribbons: An analytical approach, Phys. Rev. B 94, 195421 (2016).
  * Cui  _et al._ [2021] L. Cui, J. Wang, and M. Sun, Graphene plasmon for optoelectronics, Reviews in Physics 6, 100054 (2021).
  * Rhouma  _et al._ [2022] M. B. Rhouma, B. Guizal, P. Bonnet, F. Paladian, and K. Edee, Semi-analytical model for the analysis of a magnetically biased 1d subwavelength graphene-strip-grating, Opt. Continuum 1, 1144 (2022).
  * Xiong  _et al._ [2019] L. Xiong, C. Forsythe, M. Jung, A. S. McLeod, S. S. Sunku, Y. M. Shao, G. X. Ni, A. J. Sternbach, S. Liu, J. H. Edgar, E. J. Mele, M. M. Fogler, G. Shvets, C. R. Dean, and D. N. Basov, Photonic crystal for graphene plasmons, Nature Communications 10, 4780 (2019).
  * Emani  _et al._ [2015] N. K. Emani, D. Wang, T.-F. Chung, L. J. Prokopeva, A. V. Kildishev, V. M. Shalaev, Y. P. Chen, and A. Boltasseva, Plasmon resonance in multilayer graphene nanoribbons, Laser & Photonics Reviews 9, 650 (2015).
  * Li  _et al._ [2019] G. Li, V. Semenenko, V. Perebeinos, and P. Q. Liu, Multilayer graphene terahertz plasmonic structures for enhanced frequency tuning range, ACS Photonics 6, 3180 (2019).
  * Rodrigo  _et al._ [2017] D. Rodrigo, A. Tittl, O. Limaj, F. J. G. de Abajo, V. Pruneri, and H. Altug, Double-layer graphene for enhanced tunable infrared plasmonics, Light: Science & Applications 6, e16277 (2017).
  * Yan  _et al._ [2012] H. Yan, X. Li, B. Chandra, G. Tulevski, Y. Wu, M. Freitag, W. Zhu, P. Avouris, and F. Xia, Tunable infrared plasmonic devices using graphene/insulator stacks, Nature Nanotechnology 7, 330 (2012).
  * He  _et al._ [2014] M.-D. He, G. Zhang, J.-Q. Liu, J.-B. Li, X.-J. Wang, Z.-R. Huang, L. Wang, and X. Chen, Plasmon resonances in a stacked pair of graphene ribbon arrays with a lateral displacement, Optics Express 22, 6680 (2014).
  * Gong  _et al._ [2020] S. Gong, L. Wang, Y. Zhang, Z. Yang, X. Li, Q. Wen, Z. He, S. Liang, L. Yuan, C. Yu, Z. Feng, Z. Yang, and X. Zhang, Ultra-extraordinary optical transmission induced by cascade coupling of surface plasmon polaritons in composite graphene–dielectric stack, Optics Express 28, 30502 (2020).
  * Haddouche and Cherbi [2017] I. Haddouche and L. Cherbi, Comparison of finite element and transfer matrix methods for numerical investigation of surface plasmon waveguides, Optics Communications 382, 132 (2017).
  * Maier  _et al._ [2017] M. Maier, D. Margetis, and M. Luskin, Dipole excitation of surface plasmon on a conducting sheet: Finite element approximation and validation, Journal of Computational Physics 339, 126 (2017).
  * Jung  _et al._ [2007] J. Jung, T. Søndergaard, and S. I. Bozhevolnyi, Theoretical analysis of square surface plasmon-polariton waveguides for long-range polarization-independent waveguiding, Phys. Rev. B 76, 035434 (2007).
  * Shibayama  _et al._ [2021] J. Shibayama, J. Yamauchi, and H. Nakano, Frequency-dependent FDTD analyses of terahertz plasmonic devices, in  _2020 International Symposium on Antennas and Propagation (ISAP)_ (2021) pp. 457–458.
  * Zang  _et al._ [2020] T. Zang, H. Zang, Z. Xi, J. Du, H. Wang, Y. Lu, and P. Wang, Asymmetric excitation of surface plasmon polaritons via paired slot antennas for angstrom displacement sensing, Phys. Rev. Lett. 124, 243901 (2020).
  * Rana  _et al._ [2015] M. Rana, B. Hossain, R. Islam, and Y. G. Guo, Surface plasmon polariton propagation modeling for graphene parallel pair sheets using FDTD, in  _2015 IEEE International Conference on Applied Superconductivity and Electromagnetic Devices (ASEMD)_ (2015) pp. 179–180.
  * Olkkonen [2010] J. Olkkonen, FDTD scattered field formulation for scatterers in stratified dispersive media, Opt. Express 18, 4380 (2010).
  * Weiss  _et al._ [2011] T. Weiss, N. A. Gippius, S. G. Tikhodeev, G. Granet, and H. Giessen, Derivation of plasmonic resonances in the fourier modal method with adaptive spatial resolution and matched coordinates, J. Opt. Soc. Am. A 28, 238 (2011).
  * Anttu and Xu [2011] N. Anttu and H. Q. Xu, Scattering matrix method for optical excitation of surface plasmons in metal films with periodic arrays of subwavelength holes, Phys. Rev. B 83, 165431 (2011).
  * Maksimova  _et al._ [2012] G. M. Maksimova, E. S. Azarova, A. V. Telezhnikov, and V. A. Burdov, Graphene superlattice with periodically modulated dirac gap, Physical Review B 86, 205422 (2012).
  * Nikitin  _et al._ [2012] A. Y. Nikitin, F. Guinea, and L. Martin-Moreno, Resonant plasmonic effects in periodic graphene antidot arrays, Applied Physics Letters 101, 151119 (2012).
  * Wei  _et al._ [2019] C. Wei, R. Negishi, Y. Ogawa, M. Akabori, Y. Taniyasu, and Y. Kobayashi, Turbostratic multilayer graphene synthesis on cvd graphene template toward improving electrical performance, Japanese Journal of Applied Physics 58, SIIB04 (2019).
  * Baek  _et al._ [2017] I. H. Baek, J. M. Hamm, K. J. Ahn, B. J. Kang, S. S. Oh, S. Bae, S. Y. Choi, B. H. Hong, D.-I. Yeom, B. Min, O. Hess, Y. U. Jeong, and F. Rotermund, Boosting the terahertz nonlinearity of graphene by orientation disorder, 2D Materials 4, 025035 (2017).
  * Whittaker and Culshaw [1999] D. M. Whittaker and I. S. Culshaw, Scattering-matrix treatment of patterned multilayer photonic structures, Phys. Rev. B 60, 2610 (1999).
  * Tikhodeev  _et al._ [2002] S. G. Tikhodeev, A. L. Yablonskii, E. A. Muljarov, N. A. Gippius, and T. Ishihara, Quasiguided modes and optical properties of photonic crystal slabs, Phys. Rev. B 66, 045102 (2002).
  * Bykov and Doskolovich [2013] D. A. Bykov and L. L. Doskolovich, Numerical methods for calculating poles of the scattering matrix with applications in grating theory, Journal of Lightwave Technology 31, 793 (2013).
  * Falkovsky [2008] L. A. Falkovsky, Optical properties of graphene, Journal of Physics: Conference Series 129, 012004 (2008).
  * Min and MacDonald [2008] H. Min and A. H. MacDonald, Electronic structure of multilayer graphene, Progress of Theoretical Physics Supplement 176, 227 (2008).
  * Garlow  _et al._ [2016] J. A. Garlow, L. K. Barrett, L. Wu, K. Kisslinger, Y. Zhu, and J. F. Pulecio, Large-area growth of turbostratic graphene on Ni(111) via physical vapor deposition, Scientific Reports 6, 19804 (2016).
  * Negishi  _et al._ [2017] R. Negishi, K. Yamamoto, H. Kitakawa, M. Fukumori, H. Tanaka, T. Ogawa, and Y. Kobayashi, Synthesis of very narrow multilayer graphene nanoribbon with turbostratic stacking, Applied Physics Letters 110, 201901 (2017).
  * Bolotovskii and Lebedev [1968] B. Bolotovskii and A. Lebedev, On threshold phenomenain classical electrodynamics, Sov. Phys. JETP 26, 784 (1968).
  * Wojcik  _et al._ [2021] C. C. Wojcik, H. Wang, M. Orenstein, and S. Fan, Universal behavior of the scattering matrix near thresholds in photonics, Phys. Rev. Lett. 127, 277401 (2021).
  * Fujita  _et al._ [1998] T. Fujita, Y. Sato, T. Kuitani, and T. Ishihara, Tunable polariton absorption of distributed feedback microcavities at room temperature, Phys. Rev. B 57, 12428 (1998).
  * Yablonskii  _et al._ [2001] A. Yablonskii, E. Muljarov, N. Gippius, S. Tikhodeev, T. Fujita, and T. Ishihara, Polariton effect in distributed feedback microcavities, Journal of the Physical Society of Japan 70, 1137 (2001).


# Higher-order protection of quantum gates: Hamiltonian engineering
coordinated with dynamical decoupling

P. Z. Zhao  Center for Quantum Technologies, National University of Singapore,
Singapore 117543, Singapore  Sirui Liu  Department of Chemistry, National
University of Singapore, Singapore 117543, Singapore  Department of Physics,
National University of Singapore, Singapore 117551, Singapore  Joint School of
National University of Singapore and Tianjin University, International Campus
of Tianjin University, Binhai New City, Fuzhou 350207, China  Jiangbin Gong
<EMAIL_ADDRESS>Center for Quantum Technologies, National University of
Singapore, Singapore 117543, Singapore  Department of Physics, National
University of Singapore, Singapore 117551, Singapore  Joint School of National
University of Singapore and Tianjin University, International Campus of
Tianjin University, Binhai New City, Fuzhou 350207, China

###### Abstract

Dynamical decoupling represents an active approach towards the protection of
quantum memories and quantum gates. Because dynamical decoupling operations
can interfere with system’s own time evolution, the protection of quantum
gates is more challenging than that of quantum states. In this work, we put
forward a simple but general approach towards the realization of higher-order
protection of quantum gates. The central idea of our approach is to engineer
(hence regain the control of) the quantum gate Hamiltonian in coordination
with higher-order dynamical decoupling sequences originally proposed for the
protection of quantum memories. In our computational examples presented for
illustration, the required engineering can be implemented by only quenching
the phase of an external driving field at particular times.

Introduction.— A crucial challenge in implementing quantum computation is to
overcome the decoherence induced by the interaction between a quantum system
and its environment. As an active approach, dynamical decoupling (DD) aims to
average out system-environmental interaction and hence protects either quantum
memories or quantum gates Viola ; Viola1999 . For quantum memories, periodic
DD (PDD) Viola2003  eliminates the effect of system-environment interaction to
the first-order (e.g., in the sense of a Magnus expansion), whereas
concatenated DD (CDD) Khodjasteh  and Uhrig DD (UDD) Uhrig ; Yang ; West2010
can offer higher-order elimination of decoherence Gong ; Gong2010 .
Experimental advances on DD-based quantum memory protection have been highly
successful Biercuk ; Du ; Ryan ; deLange ; Wang ; Alvarez , with recent
benchmarking studies of DD carried out on actual noisy quantum computers
Ezzell . However, most of these experimental progresses on quantum memory
protection have not been translated to the protection of quantum gates,
because DD operations often interfere with quantum gate operations, e.g., DD
pulses tend to freeze system’s own time evolution.

One solution towards the DD protection of quantum gates started with certain
driving Hamiltonian commutable with DD operations Byrd ; Lidar2008 ; West ,
the so-called “non-interference” condition. This stimulating route however
needs extra physical qubit resources to encode logical qubits, the
implementation of which is often based on particular dynamical symmetries or
specific physical interactions to meet the non-interference requirement Xu ;
Zhao ; Wu . The second solution utilizes gate operations inserted in between
DD pulses. This approach typically extends the total gate time and can only
realize first-order protection Suter ; Zhang ; Souza . To further fill in the
gap between quantum memory protection and quantum gate protection by DD, the
dynamically corrected gate (DCG) was proposed as a sophiticated composite
quantum gate constructed from primitive building blocks without encoding
Viola2009 ; ViolaPRA ; Ng . DCG may offer higher-order protection with
concatenate design Viola2010 , but the total number of the required gate
operations will be large and thus experimentally challenging. It remains to
see if it is possible to reduce the overhead in the realization of higher-
order protection of quantum gates.

This work aims to show that higher-order protection of quantum gates can be
directly achieved on general physical qubits, with straightforward driving
Hamiltonians that are piecewise continuous and at no cost of the total gate
time. The central idea of our higher-order approach lies in the engineering of
the time dependence of system’s driving Hamiltonian Viola2003 , in
coordination with higher-order DD sequences originally designed for quantum
memory protection, such as CDD, UDD, and even nested UDD. Because the
realization of DD pulses already requires to design explicit and rapid time
dependence of system’s Hamiltonian, engineering some additional time
dependence of system’s driving Hamiltonian does not present a challenge when
compared with the task of qubit encoding. As seen below, the required
engineering of system’s driving Hamiltonian is to act against the unwanted
impact of DD pulses on the system. This being done, only the system-
environment interaction part is suppressed to higher orders and we recover our
control over system’s driving Hamiltonian. The DD protection of quantum
memories and quantum gates can thus be unified and applied on an equal
footing. Our theoretical proposal is verified by our computational examples
below.

Approach.— Consider a quantum system coupled to its environment with the total
Hamiltonian $H(t)=H_{\rm S}(t)+H_{\rm E}+H_{\rm I}$. Here, $H_{\rm S}(t)$ is
our system’s driving Hamiltonian, $H_{\rm E}$ is the environment Hamiltonian
and $H_{\rm I}$ is the system-environment interaction. The most general
interaction between a qubit and its environment is in the form $H_{\rm
I}=\sigma_{x}\otimes{E_{x}}+\sigma_{y}\otimes{E_{y}}+\sigma_{z}\otimes{E_{z}}$,
where ${E}_{x(y,z)}$ represent the associated environment operators. To
illustrate the underlying principle of our idea, we start with a one-qubit
case under PDD pulses, where

$H_{\rm S}(t)=\Omega_{x}(t)\sigma_{x}+\Omega_{y}(t)\sigma_{y},$ (1)

with $\Omega_{x}(t)$ and $\Omega_{y}(t)$ being time-dependent system
parameters.

Let us now consider a periodic sequence of fast and strong symmetrizing pulses
with the DD operations $\\{\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3}\\}$
($\sigma_{0}\equiv{I}$, $\sigma_{1}\equiv{\sigma_{x}}$,
$\sigma_{2}\equiv{\sigma_{y}}$ and $\sigma_{3}\equiv{\sigma_{z}}$ with $I$
corresponding to no pulse applied), over a duration of $4\tau$, applied to the
native evolution governed by the driving Hamiltonian $H_{\rm S}(t)$. Assuming
that the applied pulses are instantaneous (that is, much shorter than the
total gate time $4\tau$), the time evolution operator is given by

$\displaystyle\mathcal{U}=$
$\displaystyle\prod^{3}_{k=0}\sigma_{k}\mathcal{T}e^{-i\int^{(k+1)\tau}_{k\tau}\mathcal{H}(t)dt}\sigma_{k}$
$\displaystyle=$
$\displaystyle\prod^{3}_{k=0}\mathcal{T}e^{-i\left[\int^{(k+1)\tau}_{k\tau}\sigma_{k}H_{\rm
S}(t)\sigma_{k}dt+\sigma_{k}H_{\rm I}\sigma_{k}\tau+H_{\rm E}\tau\right]}$
$\displaystyle=$
$\displaystyle\prod^{3}_{k=0}\mathcal{T}e^{-i\int^{(k+1)\tau}_{k\tau}\sigma_{k}H_{\rm
S}(t)\sigma_{k}dt}e^{-i\left(\sum^{3}_{l=0}\sigma_{l}H_{\rm
I}\sigma_{l}+4H_{\rm E}\tau\right)}+\mathcal{O}(\tau^{2}),$ (2)

where the term with the smallest $k$ is placed at the right most and
$\mathcal{T}$ denotes the time ordering. Note that
$\sum^{3}_{l=0}\sigma_{l}H_{\rm I}\sigma_{l}=0$, so we have

$\displaystyle\mathcal{U}=$
$\displaystyle\prod^{3}_{k=0}\mathcal{T}e^{-i\int^{(k+1)\tau}_{k\tau}\sigma_{k}H_{\rm
S}(t)\sigma_{k}dt}{\otimes}e^{-i4H_{\rm E}\tau}+\mathcal{O}(\tau^{2}).$ (3)

Clearly then, the first-order effect arising from $H_{\rm I}$ is eliminated.
It is also seen that the time evolution of the quantum system is strongly
influenced by the DD pulses. Indeed, during each time interval of duration
$\tau$, the effective Hamiltonian [see the exponentials in Eq. (3)] is updated
by the DD pulses.

We now unravel a key observation behind our Hamiltonian engineering approach.
Denote the effective Hamiltonians altered by DD pulses associated with the
time interval $[k\tau,(k+1)\tau]$ as $H_{k}(t)\equiv\sigma_{k}H_{\rm
S}(t)\sigma_{k}$. Using Eq. (1), one obtains
$H_{0}(t)=\Omega_{x}(t)\sigma_{x}+\Omega_{y}(t)\sigma_{y}$,
$H_{1}(t)=\Omega_{x}(t)\sigma_{x}-\Omega_{y}(t)\sigma_{y}$,
$H_{2}(t)=-\Omega_{x}(t)\sigma_{x}+\Omega_{y}(t)\sigma_{y}$, and
$H_{3}(t)=-\Omega_{x}(t)\sigma_{x}-\Omega_{y}(t)\sigma_{y}$. To fight against
these changes to the system Hamiltonian imposed by DD pulses, we can quench
the parameters of $H_{\rm S}(t)$ during the respective time intervals.
Specifically for $t\in[0,\tau]$, $\Omega_{x}(t)=\Omega(t)\cos\varphi$ and
$\Omega_{y}(t)=\Omega(t)\sin\varphi$, for $t\in(\tau,2\tau]$,
$\Omega_{x}(t)=\Omega(t)\cos\varphi$ and
$\Omega_{y}(t)=-\Omega(t)\sin\varphi$, for $t\in(2\tau,3\tau]$,
$\Omega_{x}(t)=-\Omega(t)\cos\varphi$ and
$\Omega_{y}(t)=\Omega(t)\sin\varphi$, and finally, for $t\in(3\tau,4\tau]$,
$\Omega_{x}(t)=-\Omega(t)\cos\varphi$ and
$\Omega_{y}(t)=-\Omega(t)\sin\varphi$. This can be experimentally implemented
because $\Omega(t)$ is identified as the Rabi frequency of the qubit
Hamiltonian and $\varphi$ represents a phase of the driving field coupling the
two levels of the qubit. The proposed quenching can be realized by simply
adjusting the phase parameter of the driving field, namely, the phase
parameter should be chosen to be $\varphi$ for $t\in(0,\tau]$, $-\varphi$ for
$t\in(\tau,2\tau]$, $-\varphi+\pi$ for $t\in(2\tau,3\tau]$ and $\varphi+\pi$
for $t\in(3\tau,4\tau]$. With phase parameters quenched in this manner, one
has $H_{k}(t)=\Omega(t)(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})$. That
is, the effective Hamiltonians appearing in the exponentials in Eq. (3) become
the same Hamiltonian $\Omega(t)(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})$.
If we further require $\int^{4\tau}_{0}\Omega(t)dt=\theta/2$, Eq. (3)
indicates

$\displaystyle\mathcal{U}=e^{-i\theta(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})/2}\otimes
e^{-i4H_{E}\tau}+\mathcal{O}(\tau^{2}).$ (4)

As such, to the first order of $\tau$, system-environment coupling is
eliminated and yet we still implement a gate operation, namely, rotating the
qubit along an axis in the $xy$ plane by angle $\theta$. Remarkably, this is
achieved without requiring that $H_{\rm S}(t)$ commutes with the DD
operations.

We are now ready to engineer the time dependence of the qubit Hamiltonian for
higher-order protection of quantum gates. As an example, we consider the
second-order protection of a quantum gate by CDD with two layers. To that end,
we use the time evolution protected by the DD operations
$\\{\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3}\\}$ as the basic building
block and then nest it with a second layer of DD operations
$\\{\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3}\\}$. In doing so, we examine
how the time dependence of $H_{\rm S}(t)$ should be engineered. The evolution
operator with two layers of DD operations is given by

$\displaystyle\mathcal{U}=$
$\displaystyle\prod^{3}_{m=0}\sigma_{m}\left[\prod^{3}_{k=0}\mathcal{T}e^{-i\int^{(k+1)\tau/4}_{k\tau/4}\sigma_{k}H_{\rm
S}(t)\sigma_{k}dt}\right]\sigma_{m}\otimes{e^{-i4H_{E}\tau}}+\mathcal{O}(\tau^{3})$
$\displaystyle=$
$\displaystyle\prod^{3}_{m=0}\prod^{3}_{k=0}\mathcal{T}e^{-i\int^{(k+1)\tau/4}_{k\tau/4}\sigma_{m}\sigma_{k}H_{\rm
S}(t)\sigma_{k}\sigma_{m}dt}\otimes{e^{-i4H_{E}\tau}}+\mathcal{O}(\tau^{3}),$
(5)

where all the terms in the product should be arranged from right to left with
increasing $k$ or $m$. We can analogously define the function appearing in the
exponential as the effective Hamiltonian
$H_{mk}(t)\equiv\sigma_{m}\sigma_{k}H_{\rm S}(t)\sigma_{k}\sigma_{m}$, with
the explicit form easily found by use of Eq. (1). Just like the previous
example showcasing first-order protection, here we can also quench the phase
parameter $\varphi$ without changing the Rabi frequency parameter $\Omega(t)$
during time intervals as well as subintervals such that
$H_{mk}(t)=\Omega(t)(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})$. With
$\int^{4\tau}_{0}\Omega(t)dt=\theta/2$, we have

$\displaystyle\mathcal{U}=e^{-i\theta(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})/2}\otimes
e^{-i4H_{E}\tau}+\mathcal{O}(\tau^{3}),$ (6)

indicating that our one-qubit gate with second-order protection is achieved.
It is straightforward to extend our approach to cases of higher-order
protection of quantum gates by increasing the concatenation level of CDD.

Let us now turn to more efficient higher-order protection of single-qubit
quantum gates, using the same qubit Hamiltonian $H_{\rm S}(t)$ introduced
above. According to the UDD theory proposed for quantum memory protection, the
pure dephasing term $\sigma_{z}\otimes{E_{z}}$ can be filtered by applying
nonequidistant decoupling operator $\sigma_{x}$. Likewise, the longitudinal
relaxation term $\sigma_{x}\otimes{E_{x}}+\sigma_{y}\otimes{E_{y}}$ can be
suppressed by applying $\sigma_{z}$ operations at special UDD timings Uhrig ;
Yang ; West2010 . Building upon this, in order to average out the general
system-environmental interaction
$H_{I}=\sigma_{x}\otimes{E_{x}}+\sigma_{y}\otimes{E_{y}}+\sigma_{z}\otimes{E_{z}}$,
we apply $n$ decoupling operations $\sigma_{x}$ to the quantum system at times
$t_{j}=4\tau\sin^{2}[j\pi/(2n+2)]$, and during each time interval
$[t_{j-1},t_{j}]$, we further apply $n$ decoupling operations $\sigma_{z}$ at
times $t^{j}_{k}=t_{j}\sin^{2}[k\pi/(2n+2)]$, where $j,\ k=1,...,n$ and $n$ is
taken as an even number to facilitate our engineering approach. Such a nested
UDD protocol yields the following time evolution operator:

$\displaystyle\mathcal{U}=$ $\displaystyle\mathcal{U}({4\tau-
t_{n}})\prod^{n}_{j=1}\sigma_{x}\mathcal{U}(t_{j}-t_{j-1})$ $\displaystyle=$
$\displaystyle\prod^{n+1}_{j=1}\sigma^{j-1}_{x}\mathcal{U}(t_{j}-t_{j-1})\sigma^{j-1}_{x}$
(7)

where $t_{0}\equiv 0$, $t_{n+1}\equiv 4\tau$, all the terms in the product
from right to left should have an increasing $j$, and
$\mathcal{U}(t_{j}-t_{j-1})$, the unitary evolution operator associated with
the interval $[t_{j-1},t_{j}]$, is expressed by the subinterval evolution
operators as

$\displaystyle\mathcal{U}(t_{j}-t_{j-1})=$
$\displaystyle\mathcal{U}({t_{j}-t^{j}_{n}})\prod^{n}_{k=1}\sigma_{z}\mathcal{U}(t^{j}_{k}-t^{j}_{k-1})$
$\displaystyle=$
$\displaystyle\prod^{n+1}_{k=1}\sigma^{k-1}_{z}\mathcal{U}(t^{j}_{k}-t^{j}_{k-1})\sigma^{k-1}_{z}$
(8)

with $t^{j}_{0}\equiv t_{j-1}$ and $t^{j}_{n+1}\equiv t_{j}$ and all the terms
in the product from right to left should have an increasing $k$. The UDD
theory can guarantee that the effect of system-environment interaction is
suppressed up to the $n$th order. That is,

$\displaystyle\mathcal{U}=$
$\displaystyle\prod^{n+1}_{j=1}\prod^{n+1}_{k=1}\mathcal{T}e^{-i\int^{t^{j}_{k}}_{t^{j}_{k-1}}\sigma^{j-1}_{x}\sigma^{k-1}_{z}H_{\rm
S}(t)\sigma^{k-1}_{z}\sigma^{j-1}_{x}dt}$
$\displaystyle\otimes{e}^{-i4H_{E}\tau}+\mathcal{O}(\tau^{n+1}).$ (9)

To find out how to engineer the time dependence of $H_{\rm S}(t)$, we define
effective Hamiltonians
$H^{j}_{k}(t)\equiv\sigma^{j-1}_{x}\sigma^{k-1}_{z}H_{\rm
S}(t)\sigma^{k-1}_{z}\sigma^{j-1}_{x}$, i.e., those terms appearing in the
exponentials in the equation above. Further using Eq. (1), one obtains
$H^{j}_{k}(t)=(-1)^{k-1}\Omega_{x}(t)\sigma_{x}+(-1)^{k+j}\Omega_{y}(t)\sigma_{y}$,
explicitly depicting the impact of our nested UDD operations on the system
Hamiltonian. To combat against these unwanted changes to the system
Hamiltonian, we can just adjust only the phase parameter of the qubit driving
field, such that for $t\in[t^{j}_{k-1},t^{j}_{k}]$,
$\Omega_{x}(t)=(-1)^{k-1}\Omega(t)\cos\varphi$ and
$\Omega_{y}(t)=(-1)^{k+j}\Omega(t)\sin\varphi$. The evident outcome of this
phase quenching approach is that all the effective Hamiltonians will be the
same, namely,
$H^{j}_{k}(t)=\Omega(t)(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})$.
Assuming $\int^{4\tau}_{0}\Omega(t)dt=\theta/2$, then one arrives at

$\displaystyle\mathcal{U}=e^{-i\theta(\cos\varphi\sigma_{x}+\sin\varphi\sigma_{y})/2}\otimes
e^{-i4H_{E}\tau}+\mathcal{O}(\tau^{n+1}),$ (10)

which is a quantum gate protected by nested UDD to the $n$th order. To our
knowledge, such a possibility of higher-order UDD protection of quantum gates
is never shown before. In addition, the quantum gate time is still the same as
in the bare case without protection. To see that more explicitly, let $\Omega$
be the average Rabi frequency during $[0,4\tau]$, namely,
$\Omega=[\int^{4\tau}_{0}\Omega(t)dt]/(4\tau)$. Then for all the cases
illustrated above, the total gate time $4\tau$ remains to be
$\theta/(2\Omega)$.

Figure 1: Quantum gate fidelity as a function of the system-environment
coupling strength $\epsilon$. (a) Fidelity $F$ with the first-order protection
(blue), second-order protection (red), and without protection (black) for
$\epsilon\in[0,0.5\Omega]$. (b) More results of fidelity $F$ under the second-
order protection with CDD for $\epsilon\in[0,2\Omega]$.

Numerical results.— To verify our engineering concept, let us consider a
simple form of system-environment interaction of the following Heisenberg
coupling form

$\displaystyle H_{\rm
I}=\epsilon(\sigma_{x}^{s}\sigma_{x}^{b}+\sigma_{y}^{s}\sigma_{y}^{b}+\sigma_{z}^{s}\sigma_{z}^{b}),$
(11)

with the first spin represents the qubit Hamiltonian and the second spin
represents a quantum bath, and $\epsilon$ being the coupling strength. We also
assume that there is no driving field acting solely on the environment spin.
Admittedly, this $H_{\rm I}$ is extremely simple as compared with a real
quantum bath, but it does capture the essence that the qubit system is subject
to both dephasing and population relaxation and it suffices to illustrate how
our engineering approach leads to higher-order protection. $H_{\rm I}$ may
also model some residue interaction between a qubit under gate operation and
its surrounding qubit.

Table 1: Computational results of $[\epsilon,F_{\mathrm{bare}},F_{\mathrm{PDD}},F_{\mathrm{CDD}}]$, namely, quantum gate fidelity without protection, with first order PDD and second-order CDD vs the system-environment coupling strength $\epsilon$. $\epsilon/\Omega$ |  | $0.05$ |  | $0.1$ |  | $0.2$ |  | $0.3$ |  | $0.4$ |  | $0.5$ |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$F_{\mathrm{bare}}$ |  | $98.75\%$ |  | $95.05\%$ |  | 81.47% |  | $63.57\%$ |  | 47.45% |  | $38.67\%$ |
$F_{\mathrm{PDD}}$ |  | $99.90\%$ |  | $99.61\%$ |  | 98.39% |  | $96.31\%$ |  | 93.37% |  | $89.60\%$ |
$F_{\mathrm{CDD}}$ |  | $99.99\%$ |  | $99.98\%$ |  | $99.92\%$ |  | $99.82\%$ |  | 99.67% |  | $99.49\%$ |

The performance of our higher-order protection is characterized by the
fidelity $F=\langle\phi|\rho|\phi\rangle$. Here, $|\phi\rangle$ is the target
output state and $\rho$ is the real output state that is obtained by first
solving the Liouville equation $\dot{\varrho}(t)=-i[H(t),\varrho(t)]$ and then
performing a partial trace over the environment spin. We take the initial
state as $(|0\rangle+|1\rangle)/\sqrt{2}$, the parameters of the driving
Hamiltonian as $\Omega=10\pi\ \mathrm{MHz}$ and $\varphi=\pi/4$, and the
evolution time as $0.05\mathrm{\mu{s}}$. The ideal quantum gate under these
choices are hence $U=(\sigma_{x}+\sigma_{y})/\sqrt{2}$ and the target output
state is
$|\phi\rangle=[\exp(-i\pi/4)|0\rangle+\exp(i\pi/4)|1\rangle]/\sqrt{2}$. To
illustrate how higher-order protection can be achieved, we tune the system-
environment coupling strength $\epsilon$ over a wide range, namely,
$\epsilon\in[0,0.5\Omega]$. We then plot the fidelity $F$ vs the coupling
strength $\epsilon$ for the quantum gate under PDD and CDD with first-order
and second-order protection, as shown by the blue and red lines in Fig. 1(a).
As a comparison, we also plot the gate fidelity without any DD protection, as
shown by the black line in Fig. 1(a).

Figure 1 illustrates that while the first-order protection yields gate
fidelities much better than that of the bare gate, the second-order protection
further dramatically improves the performance of the quantum gate. In
particular, for the worse case $\epsilon=2\Omega$ (the system-bath coupling
strength is twice of the Rabi frequency), the gate fidelity with second-order
protection can still reach more than $92\%$, shown in Fig. 1(b). To appreciate
the DD improvement more quantitatively, we present in table 1 our
computational findings based on our simple model. As seen from table 1 , for a
moderate coupling strength $\epsilon=0.05\Omega$, the gate fidelities under
first-order and second-order protection can reach $F_{\mathrm{PDD}}=99.90\%$
and $F_{\mathrm{CDD}}=99.99\%$, respectively, both being much higher than the
fidelity $F_{\mathrm{bare}}=98.75\%$ without DD protection. Notably, for the
case $\epsilon=0.5\Omega$, i.e., the coupling strength is approximately half
of the Rabi frequency of gate Hamiltonian such that the fidelity of the bare
gate is as low as $38.67\%$ (the gate is basically totally destroyed), the
first-order protection can boost the fidelity back to $90\%$ (which may not
suffice for an actual quantum computation application), but the second-order
protection can further improve the fidelity to $99.49\%$. These numerical
findings do not represent at all what one can achieve on an actual quantum
computing platform, but indicating the promise of higher-order protection of
gate operations, now made possible by our engineering approach.

Figure 2: Quantum gate fidelity as a function of the system-environment
coupling strength over a wide range of $\epsilon\in[0,2\Omega]$ under nested
UDD protection. Black, blue and red lines respectively correspond to the
protection under the orders of $n=2$, $4$ and $6$.

Next let us check the performance of our approach when using nested UDD to
protect the quantum gate to the $n$th order according to our protocol outlined
above. Figure 2 presents the corresponding gate fidelity for a wide range of
$\epsilon\in[0,2\Omega]$, with $n=2$, $4$ and $6$. The gate fidelity is seen
to be greatly enhanced when we increase the order of protection. We further
present our numerics in table 2. As seen from table 2, for
$\epsilon\in[0,2\Omega]$, the gate fidelity with $4$th and $6$th protection
always exceed $98\%$. Remarkably, the fidelity with $6$th order protection can
be as high as $99.54\%$ even when the system-environment coupling strength is
twice of the Rabi frequency of the system Hamiltonian.

Table 2: Computational results of $[\epsilon,F_{\rm UDD-2},F_{\rm UDD-4},F_{\rm UDD-6}]$, namely, quantum gate fidelity under second-order, fourth-order, and sixth-order UDD protection vs the system-environment coupling strength $\epsilon$, with a nested UDD scheme. $\epsilon/\Omega$ |  | $0.2$ |  | $0.6$ |  | $1.0$ |  | $1.4$ |  | $2.0$ |
---|---|---|---|---|---|---|---|---|---|---|---
$F_{\rm UDD-2}$ |  | $99.57\%$ |  | $94.96\%$ |  | $84.51\%$ |  | 70.18% |  | $50.28\%$ |
$F_{\rm UDD-4}$ |  | $99.9998\%$ |  | $99.995\%$ |  | $99.93\%$ |  | 99.63% |  | $98.14\%$ |
$F_{\rm UDD-6}$ |  | $99.9999\%$ |  | $99.996\%$ |  | $99.97\%$ |  | 99.88% |  | $99.54\%$ |

Two-qubit gates.— Having computationally verified that the concept underlying
our engineering approach is valid, we finally discuss how to protect a
nontrivial two-qubit gate, which is important for universal quantum
computation under DD protection. As an example, we show below how to engineer
a two-qubit Hamiltonian to realize a nontrivial two-qubit gate under second-
order protection. Let us now start to consider a system Hamiltonian different
from all single-qubit cases discussed above, namely,

$\displaystyle H_{\rm
S}=J_{1}(t)\sigma_{z}^{(1)}+J_{2}(t)\sigma_{z}^{(2)}+J(t)\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)},$
(12)

where $J_{1}$ and $J_{2}$ are the parameters of local operations respectively
acting on the first and second qubits, and $J$ is the coupling parameter
between the two qubits. We assume each qubit is subject to its own dephasing
and population relaxation. To suppress the decoherence to the second order, we
exploit the time evolution protected by the PDD operations
$\\{\sigma_{0}^{(1)}\otimes\sigma_{0}^{(2)},\sigma_{1}^{(1)}\otimes\sigma_{1}^{(2)},\sigma_{2}^{(1)}\otimes\sigma_{2}^{(2)},\sigma_{3}^{(1)}\otimes\sigma_{3}^{(2)}\\}$
as the basic building block and then nest it into a second layer of the same
DD operations. As a result, the time evolution operator over a period of
$4\tau$ is found to be

$\displaystyle\mathcal{U}=$
$\displaystyle\Bigg{[}\prod^{3}_{m=1}\sigma_{m}^{(1)}\otimes\sigma_{m}^{(2)}\prod^{3}_{k=1}\mathcal{T}e^{-i\int^{(k+1)\tau/4}_{k\tau/4}\sigma_{k}^{(1)}\otimes\sigma_{k}^{(2)}H_{\rm
S}(t)\sigma_{k}^{(1)}\otimes\sigma_{k}^{(2)}dt}$
$\displaystyle\times\sigma_{m}^{(1)}\otimes\sigma_{m}^{(2)}\Bigg{]}\otimes{e^{-i4H_{E}\tau}}+\mathcal{O}(\tau^{3})$
$\displaystyle=$
$\displaystyle\prod^{3}_{m=1}\prod^{3}_{k=1}\mathcal{T}e^{-i\int^{(k+1)\tau/4}_{k\tau/4}\sigma_{m}^{(1)}\sigma_{k}^{(1)}\otimes\sigma_{m}^{(2)}\sigma_{k}^{(2)}H_{\rm
S}(t)\sigma_{k}^{(1)}\sigma_{m}^{(1)}\otimes\sigma_{k}^{(2)}\sigma_{m}^{(2)}dt}$
$\displaystyle\otimes{e^{-i4H_{E}\tau}}+\mathcal{O}(\tau^{3}),$ (13)

Again we may define the effective Hamiltonians
$H_{mk}(t)\equiv[\sigma_{m}^{(1)}\sigma_{k}^{(1)}]\otimes[\sigma_{m}^{(2)}\sigma_{k}^{(2)}]H_{\rm
S}(t)[\sigma_{k}^{(1)}\sigma_{m}^{(1)}]\otimes[\sigma_{k}^{(2)}\sigma_{m}^{(2)}]$.
Using the expression of such effective Hamiltonians, it can be observed that
the impact of DD pulses on the system Hamiltonian will not flip the sign of
the qubit-qubit interaction term
$J(t)\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}$. In other words, only the sign
of the single-qubit terms $J_{1}(t)\sigma_{z}^{(1)}$ and
$J_{2}(t)\sigma_{z}^{(2)}$ may be changed. When this change does occur, we can
just reverse the direction of the driving field to cancel the unwanted action
of the CDD pulses. As an example where $|J_{1}(t)|=|J_{2}(t)|=J(t)$, our
quench protocol coordinated with DD pulses will lead to
$H_{mk}=-J(t)[\sigma^{(1)}_{z}+\sigma^{(2)}_{z}-\sigma^{(1)}_{z}\sigma^{(2)}_{z}]$.
For $\int^{4\tau}_{0}J(t)dt=\theta/2$, we obtain

$\displaystyle\mathcal{U}=e^{-i\theta(\sigma_{z}^{(1)}\otimes{I}+{I}\otimes\sigma_{z}^{(2)}-\sigma_{z}^{(1)}\sigma_{z}^{(2)})/2}\otimes
e^{-i4H_{E}\tau}+\mathcal{O}(\tau^{3}).$ (14)

Ignoring an unimportant global phase factor $\exp(-i\theta/2)$, this gate
under second-order protection yields

$\displaystyle U=|00\rangle\langle 00|+|01\rangle\langle 01|+|10\rangle\langle
10|+e^{i\theta}|11\rangle\langle 11|,$ (15)

which is a controlled-phase gate. Combined with protected single-qubit gates
illustrated above, we have thus conceptually shown that in principle it is
possible to execute universal quantum computation under higher-order
protection with Hamiltonian engineering.

Discussions.— We have proposed a simple approach towards the realization of
higher-order protection of quantum gates. The central idea is to quench a gate
Hamiltonian to combat against the influence of DD pulses on the time evolution
of the system. This way, on the one hand we manipulate the quantum states by
DD to average out system-bath interaction to high orders; on the other hand we
regain our control on quantum evolution and hence achieve high-order
protection of quantum gates. This work thus indicates that highly efficient
schemes such as UDD for quantum memory protection can be integrated with
quantum gates with higher-order protection. Though in this work we only
illustrated our idea using general-purpose PDD, CDD and UDD sequences, our
approach can exploit other DD designs as well. In particular, if there is
spectral information about the environment noise, then the time intervals
between DD sequences can be further optimized according to the noise spectrum
Gong11  and as a result, quantum gate protection using our concept can
directly benefit from such optimization that can further drastically improve
the performance of decoherence suppression. Our Hamiltonian engineering
approach in coordination with DD hence makes the measurement of environment
spectrum much more relevant than before to the realization of high-fidelity
quantum gates. Numerics presented in this work are not aimed to simulate a
realistic situation, but to verify that our approach is conceptually correct.
Indeed, as a price to implement our design, we have assumed instantaneous DD
control pulses as compared with the total gate time. It will be of great
interest to examine the actual performance of our simple design on real
physical platforms, including noisy intermediate scale quantum computers.

Acknowledgments.— J.G. is grateful to Lorenza Viola for technical discussions
and for her highly constructive comments on the first version of our
manuscript. This work was supported by the National Research Foundation,
Singapore and A*STAR under its CQT Bridging Grant.

## References

  * (1) L. Viola, E. Knill, and S. Lloyd, Phys. Rev. Lett. 82, 2417 (1999).
  * (2) L. Viola, S. Lloyd, and E. Knill, Phys. Rev. Lett. 83, 4888 (1999).
  * (3) L. Viola and E. Knill, Phys. Rev. Lett. 90, 037901 (2003).
  * (4) K. Khodjasteh and D. A. Lidar, Phys. Rev. Lett. 95, 180501 (2005).
  * (5) G. S. Uhrig, Phys. Rev. Lett. 98, 100504 (2007).
  * (6) W. Yang and R. B. Liu, Phys. Rev. Lett. 101, 180403 (2008).
  * (7) J. R. West, B. H. Fong, and D. A. Lidar, Phys. Rev. Lett. 104, 1130501 (2010).
  * (8) M. Mukhtar, T. B. Saw, W. T. Soh, and J. Gong, Phys. Rev. A 81, 012331 (2010).
  * (9) M. Mukhtar, W. T. Soh, T. B. Saw, and J. Gong, Phys. Rev. A 82, 052338 (2010).
  * (10) M. J. Biercuk, H. Uys, A. P. VanDevender, N. Shiga, W. M. Itano, and J. J. Bollinger, Nature 458, 996 (2009).
  * (11) J. F. Du, X. Rong, N. Zhao, Y. Wang, J. H. Yang, and R. B. Liu, Nature 461, 1265 (2009).
  * (12) C. A. Ryan, J. S. Hodges, and D. G. Cory, Phys. Rev. Lett. 105, 200402 (2010).
  * (13) G. deLange, Z. H. Wang, D. Risté, V. V. Dobrovitski, and R. Hanson, Science 330, 60 (2010).
  * (14) Y. Wang, X. Rong, P. B. Feng, W. J. Xu, B. Chong, J. H. Su, J. Gong, and J. F. Du, Phys. Rev. Lett. 106, 040501 (2011).
  * (15) G. A. Álvarez, A. M. Souza, and D. Suter, Phys. Rev. A 85, 052324 (2012).
  * (16) N. Ezzell, B. Pokharel, L. Tewala, G. Quiroz, and D. A. Lidar, arXiv: 2207.03670.
  * (17) M. S. Byrd and D. A. Lidar, Phys. Rev. Lett. 89, 047901 (2002).
  * (18) D. A. Lidar, Phys. Rev. Lett. 100, 160506 (2008).
  * (19) J. R. West, D. A. Lidar, B. H. Fong, and M. F. Gyure, Phys. Rev. Lett. 105, 230503 (2010).
  * (20) G. F. Xu and G. L. Long, Phys. Rev. A 90, 022323 (2014).
  * (21) P. Z. Zhao , X. Wu, and D. M. Tong, Phys. Rev. A 103, 012205 (2021).
  * (22) C. F. Wu, C. F. Sun, J. L. Chen, and X. X. Yi , Phys. Rev. Appl. 19, 034069 (2023).
  * (23) A. M. Souza, G. A. Álvarez, and D. Suter, Phys. Rev. A 86, 050301(R) (2012).
  * (24) J. F. Zhang, A. M. Souza, F. D. Brandao, and D. Suter, Phys. Rev. Lett. 112, 050502 (2014).
  * (25) D. Suter and G. A. Álvarez, Rev. Mod. Phys. 88, 041001 (2016).
  * (26) K. Khodjasteh and L. Viola, Phys. Rev. Lett. 102, 080501 (2009).
  * (27) K. Khodjasteh and L. Viola, Phys. Rev. A 80, 032314 (2009).
  * (28) H. K. Ng, D. A. Lidar, and J. Preskill , Phys. Rev. A 84, 012305 (2011).
  * (29) K. Khodjasteh, D. A. Lidar, and L. Viola, Phys. Rev. Lett. 104, 090501 (2010).
  * (30) Y. Pan, Z. R. Xi, and J. Gong, J. Phys. B 44, 175501 (2011).