# El Niño Modoki thus far can be mostly predicted
more than 10 years ahead of time

X. San Liang  X.S. Liang<EMAIL_ADDRESS>http://www.ncoads.org/
Nanjing Institute of Meteorology, Nanjing 210044, China  Fen Xu  Nanjing
Institute of Meteorology, Nanjing 210044, China  Yineng Rong  Nanjing
Institute of Meteorology, Nanjing 210044, China

(January 4, 2021)

###### Abstract

The 2014-2015 “Monster”/“Super” El Niño failed to be predicted one year
earlier due to the growing importance of a new type of El Niño, El Niño
Modoki, which reportedly has much lower forecast skill with the classical
models. In this study, we show that, so far as of today, this new El Niño
actually can be mostly predicted at a lead time of more than 10 years. This is
achieved through tracing the predictability source with an information flow-
based causality analysis, which is rigorously established from first
principles in the past decade. We show that the information flowing from the
solar activity 45 years ago to the sea surface temperature results in a causal
structure resembling the El Niño Modoki mode. Based on this, a
multidimensional system is constructed out of the sunspot number series with
time delays of 22-50 years. The first 25 principal components are then taken
as the predictors to fulfill the prediction, which through information flow-
based causal deep learning reproduces rather accurately the events 12 years in
advance.

El Niño Modoki; Causal pattern; Information flow; Solar activity; Causal AI

Revised version submitted to Nature’s Scientific Reports

Submission history:
12/16/2019, Nature, MS#2019-12-18780, rejected 12/19/2019
01/17/2020, PNAS, MS#2020-01004, rejected 02/11/2020
02/19/2020, Nature Communications, MS#NCOMMS-20-06318-T, rejected 02/21/2020
03/06/2020, Science Advances, MS#abb6053, rejected 03/07/2020
04/17/2020, Geophys. Res. Lett., MS#2020GL088370, rejected 04/25/2020
08/12/2020, Nature Climate Change, MS#NCLIM-20082068, rejected 08/15/2020
08/15/2020, Nature’s Scientific Reports,
MS#affcf486-cbcb-4658-9c9d-289452a02b30
(Alternative title used in the submissions: Attribution in the presence of
hidden climate drivers—Application shows that El Niño Modoki is predictable at
a lead time of over 10 years)

The 2014-2015 “Monster”/“Super” El Niño failed to be predicted one year
earlier due to the growing importance of a new type of El Niño, El Niño
Modoki, which reportedly has much lower forecast skill with the classical
models. In this study, we show that, so far as of today, this new El Niño
actually can be mostly predicted at a lead time of more than 10 years. This is
achieved through tracing the predictability source with an information flow-
based causality analysis, which is rigorously established from first
principles in the past decade. We show that the information flowing from the
solar activity 45 years ago to the sea surface temperature results in a causal
structure resembling the El Niño Modoki mode. Based on this, a
multidimensional system is constructed out of the sunspot number series with
time delays of 22-50 years. The first 25 principal components are then taken
as the predictors to fulfill the prediction, which through information flow-
based causal deep learning reproduces rather accurately the events 12 years in
advance.

El Niño prediction has become a routine practice in operational centers all
over the world because it helps make decisions in many sectors of our society
such as agriculture, hydrology, health, energy, to name a few. Recently,
however, a number of projections fell off the mark. For example, in 2014, the
initially projected “Monster El Niño” did not show up as expected; in the
meantime, almost no model successfully predicted the 2015 super El Niño at a
1-year lead timeMcPhaden2015 . Now it has been gradually clear that there are
different types of El Niño; particularly, there exists a new type that warms
the Pacific mainly in the center, and tends to be more unpredictable than the
traditional El Niño with the present coupled modelsAshok2009 . This
phenomenon, which has caught much attention recentlyTrenberth2001 Larkin2005
Yu2007 Ashok2007 Kug2009  (c.f. Fu1986  for an earlier account), has been
termed El Niño ModokiAshok2007 , Central-Pacific (CP) type El Niño Yu2007 ,
Date Line El Niño Larkin2005 , Warm Pool El Niño Kug2009 , etc.; see Wang2017
for a review. El Niño Modoki has left climate imprints which are distinctly
different from that caused by the canonical El Niño. For example, during El
Niño Modoki episodes, the western coast region of the United States suffers
from severe drought, whereas during canonical El Niño periods it is usually
wetBehera2018 ; in the western Pacific, opposite impacts have also been
realized for the tropical cyclone activities during the two types of El Niño
periods. El Niño Modoki can reduce more effectively the Indian monsoon
rainfall, exerting influences on the rainfall in Australia and southern China
in a way different from that canonical El Niño does; also notably, during the
2009 El Niño Modoki, a stationary anticyclone is induced outside Antarctica,
causing the melting the iceMcPhaden2015 .

El Niño Modoki prediction is difficult because its generating mechanisms are
still largely unknown. Without a correct and complete attribution, the
dynamical models may not have all the adequate dynamics embedded (e.g.,
Barnston2005 ). Though recently it is reported to correspond to one of the two
most unstable modesXie2018  for the Zebiak-Cane modelZebiak1987 , the
excitation of the mode, if true, remains elusive. It has been argued that the
central Pacific warming is due to ocean advectionKug2009 , wintertime
midlatitude atmospheric variationsYu2010 , wind-induced thermocline
variationsAshok2007 , westerly wind burstsChen2015 , to name a few. But these
proposed mechanisms are yet to be verified. One approach to verification is to
assess the importance of each factor through sensitivity experiments with
numerical models. The difficulty here is, a numerical climate system may
involve many components with parameters yet to be determined empirically, let
alone the highly nonlinear components, say, the fluid earth subsystem, may be
intrinsically uncertain. Besides, for a phenomenon with physics largely
unknown, whether the model setup is dynamically consistent is always a
concern. For El Niño Modoki, it has been reported that the present coupled
models tend to have lower forecast skill for El Niño Modoki than that for
canonical El Niño Ashok2009 , probably due to the fact that, during El Niño
Modoki, the atmosphere and ocean fail to connect. For all the above, without a
correct attribution, it is not surprising to see that the 2014-15 forecasts
fell off the mark–Essentially no model predicted the 2015 El Niño at a one-
year lead time (cf. Tang2018 ).

Rather than based on dynamical models, an alternative approach relies heavily
on observations and gradually becomes popular as data accumulate. This
practice, commonly known as statistical forecast in climate science (e.g.,
vonStorch2001 ) actually appears in a similar fashion in many other
disciplines, and now has evolved into a new science, namely, data science.
Data science is growing rapidly during the past few years, one important topic
being predictability. Recently it has been found that predictability can be
transferred between dynamical events, and the transfer of predictability can
be rigorously derived from first principles in physicsLK2005 Liang2016
Liang2018 . The finding is expected to help unravel the origin(s) of
predictability for El Niño Modoki, and hence guide its prediction. In this
study, we introduce this systematically developed theory, which was in fact
born from atmosphere-ocean science, and, after numerous experiments, show that
El Niño Modoki actually can be mostly predicted at a lead time of 13 years or
over, in contrast to the canonical El Niño forecasts. We emphasize that it is
NOT our intention to do dynamical attributions here. We just present a fact on
accomplished predictions. There is still a long way to go in pursuit of the
dynamical processes underlying the new type of El Niño.

Seeking for predictor(s) for El Niño Modoki with information flow analysis.
The major research methodology is based on a recently rigorously developed
theory of information flow, a fundamental notion in physics which has
applications in a wide variety of scientific disciplines (see Liang2016  and
references therein). Its importance lies beyond the literal meaning in that it
implies causation, transfer of predictability, among others. Recently, it has
just been found to be derivable from first principles (rather than
axiomatically proposed), initially motivated by the predictability study in
atmosphere-ocean scienceLK2005 . Since its birth it has been validated with
many benchmark dynamical systems (cf. Liang2016 ), and has been successfully
applied to different disciplines such as earth system science, neuroscience,
quantitative finance, etc. (cf. Vannitsem2019 , Dionisis2019 , Stips2016 ).
Refer to the Method Section for a brief introduction.

We now apply the information flow analysis to find where the predictability of
El Niño Modoki is from. Based on this it is henceforth possible to select
covariate(s) for the prediction. Covariate selection, also known as sparsity
identification, is performed before parameter estimation. It arises in all
kinds of engineering applications such as helicopter controlJohnson2005  and
other mechanical system modelingBollt2017 , as well as in climate science. We
particularly need to use the formula (7) or its two-dimensional (2D) version
(8) as shown in the Method Section. The data used are referred to the Data
Availability Section.

At the first step we sort of use brute force: compute the information flows
from the index series of some prospect driver to the sea surface temperature
(SST) at each grid point of the tropical Pacific. Because of its quantitative
nature, the resulting distribution of information flow will form a structure
which we will refer to “causal structure” henceforth. While it may be within
expectation that information flow may exist (different from zero) for the SST
at some grid points, it is not coincidental if the flows at all grid points
organize into a certain pattern/structure. (This is similar to how
teleconnection patterns are identifiedWallace1981 .) We then check whether the
resulting causal structure resembles the El Niño Modoki pattern, say, the
pattern in the Fig. 2b of Ashok et al. (2007)Ashok2007 . We have tried many
indices (e.g. those of Indian Ocean Dipole, Pacific Decadal Oscillation, North
Atlantic, Atlantic Multi-decadal Oscillation, to name but a few), and have
found that the series of sunspot numbers (SSN) is the very one. In the
following we present the algorithm and results.

We first use the 2D version of Eq. (7) (cf. (8)) to do a rough estimate for
the information flow from SSN to the SST at each grid point, given the time
series as described in the preceding section. This is the practice how
teleconnection patterns are identified using correlation analysis, but here
the computed is information flow. (Different from the symmetric correlation,
here the information flow in the opposite direction is by computation
insignificant, just as expected.) We then form delayed time series for SSN,
given time delays $n=1,2,...$, and repeat the above step. For convenience,
denote an SSN series with delay $n$ by SSN($-n$). We find that the resulting
information flow is not significant at a 90% confidence level, or has a
structure bearing no resemblance to the desired pattern, until $n$ approaches
45 years; see Supplementary Figure S1, for a number of examples. Here the data
we are using are the SST from 01/1980 through 12/2017 (SST are most reliable
after 1979), and the SSN data are correspondingly from 01/1935 through
12/1972. Note that using the time delayed coordinates we can reconstruct a
dynamical system which is topologically equivalent to the one that originally
generates the time series, thanks to Takens’ embedding theoremTakens1980 .
Here, based on the above rough estimates, we choose to form a dynamical system
with components corresponding to the SST series and the time-delayed series
SSN($-n$) with $n$ ranging from 22-50 years. Note one is advised not to choose
two delays too close. SSN is dominated by low-frequency processes; two close
series thus-formed may not be independent enough to span a subspace, and hence
may lead to singularity. Here we choose SSN series at delays of 22-50 years
every 5 years or 60 steps. We have tried many other sampling intervals and
found this is the best. This makes sense: An interval of 5 years put the two
series approximately at the opposite phases in a solar cycle. From the thus-
generated multivariate series, we then compute the information flow from
SSN($-45$ years) to the Pacific SST, using formula (7).

A note on the embedding coordinates. In general, it is impossible to
reconstruct exactly the original dynamical system from the embedding
coordinates; it is just a topologically equivalent one. That also implies that
the reconstruction is not unique. Other choices with delayed time series may
serve the purpose as well, provided that they are independent of other
coordinates; see Abarbanel1996  for empirical methods for time delay choosing.
Fortunately, this does not make a problem for our causal inference, thanks to
a property of the information flow ab initio, which asserts that

  * The information flow between two coordinates of a system stays invariant upon arbitrary nonlinear transformation of the remaining coordinates.

More details are referred to the Method Section. This remarkable property
implies that the embedding coordinates do not matter in evaluating the
information flow from SSN to SST; what matters is the number of the
coordinates, i.e., the dimension, of the system, which makes a topological
invariant. There exist empirical ways to determine the dimension of a system,
e.g., those as shown in Abarbanel1996 . Here we choose by adding new
coordinates and examining the information flow; if it does not change any
more, the process stops and the dimension is hence determined. For this
problem, as shown above, the SSN series at delays of 22-50 years every 5 years
are chosen. This results in 6 auxiliary coordinates, which together with the
SST and SSN(-45 years) series make a system of dimension 8.

The spatial distribution of the absolute information flow, or “causal
structure” as will be referred to henceforth, from SSN($-45$ years) to the
Pacific SST is shown in Fig. 1.

Figure 1: (a) The absolute information flow from the delayed (by 45 years)
series of sunspot numbers to those of the sea surface temperature in the
Pacific Ocean (in nats/month). The pattern resembles very much the El Niño
Modoki mode as shown in the Fig. 2b of Ashok et al. (2007)Ashok2007 . Here
information flow means the transfer of predictability from one series to
another. It is computed with Eq. (7) in section 2. (b) As (a), but only the
information flow significant at a 90% confidence level is shown.

Remarkably, the causal structure in Fig. 1a is very similar to the El Niño
Modoki pattern, as given by, say, Ashok et al.Ashok2007 . (Please see their
Fig. 2b.) That is to say, the information flow from SSN to the Pacific SST
does form a spatially coherent structure resembling to El Niño Modoki. The
computed information flow is significant at a 90% confidence level, as shown
in Fig. 1b.

Comparing to the rough estimate in the bivariate case (see Supplementary
Figure S1d), Fig. 1a has a horseshoe structure pronounced in the upper part,
plus a weak branch in the Southern Hemisphere. This is just as the El Niño
Modoki structure as obtained by Ashok et al.Ashok2007  (see their Fig. 2b).
This indicates that the bivariate information flow analysis provides a good
initial guess, but the result is yet to be rectified due to the incorrect
dimensionality.

The information flows with other time delays have also been computed. Shown in
Supplementary Figure S2, are a number of examples. As can be seen, except
those around 45 years (about 44-46 years, to be precise), other delays either
yield insignificant information flows or do not result in the desired causal
pattern.

Since we are about to predict, it is desirable to perform the causality
analysis in a forecast fashion, i.e., to pretend that the EMI over the
forecast period be unavailable. We have computed the information flow with
different years of availability. Shown in Fig. S3 is the information flow as
that in Fig. 1, but with SST data available until the year of 2005. Obviously,
the causal pattern is still there and significant, and, remarkably, it appears
even more enhanced, probably due to the frequent occurrences of the event
during that period.

As information flow tells the transfer of predictability, SSN forms a natural
predictor for El Niño Modoki. This will be further confirmed in the next
section.

Projection of El Niño Modoki. The remarkable causal pattern in Fig. 1 implies
that it may be possible to make projections of El Niño Modoki many years ahead
of time. Based on it we hence conjecture that the lagged SSN series be the
desired covariate. Of course, due to the short observation period, the data
needed for parametric estimation are rather limited– This is always a problem
for climate projection.

Let us start with a simple linear regression model for the prediction of El
Niño Modoki Index (EMI) from SSN. Build the model using the EMI data from
1975-2005 ($30\times 12=600$ equations), leaving the remaining years after
2005 for prediction. For each EMI, the model inputs are the SSNs at lead times
of 50 years through 22 years (336 in total), guided by the causality analysis
in the preceding section. For best result, the annual signals and the signals
longer than 11 years are filtered from the SSN series. The filtering is
fulfilled through wavelet analysis with the orthonormal spline wavelets built
in LA2007 . Since it is required that the series length be a power of 2, we
choose a time range for the monthly SSN series from May 1847 to December 2017,
depending on the availability of the data when this research was initialized.
This totals 170 years and 8 months (170.67 years), or $2^{11}=2048$ time
steps. The upper bound scale level for the wavelet analysis is set $8$, which
gives a lower period of $2^{-8}\times 170.67=0.67$ years; the lower bound
scale level is chosen to be 4, resulting an upper period of $2^{-4}\times
170.67=10.67$ years. In doing this the seasonal cycle and the interdecadal
variabilities are effectively removed (Supplementary Figure S4). With the
pretreated SSN series the prediction is launched, and the result is shown in
Fig. 2a.

We would not say that the projection in Fig. 2a with such a simple linear
model is successful, but the 2015/2016 event is clearly seen. The strong 2009
event is not correct; it has a phase error. The general trend between
2005-2018 seems to be fine. Beyond 2018 it becomes off the mark.

Figure 2: Projection of the El Niño Modoki index based solely on the sunspot
numbers 22-50 years ago with a simple linear regression model. (a) The SSN
series is bandpass filtered with the annual signals and the signals longer
than 11 years removed. (b) Only 25 first principal components are used as
inputs. (c) The SSN series is bandpass filtered as in (a), and only the 25
first principal components are used as inputs. The shaded period is for
prediction (note in (c) the period for prediction is 1995-), and the predicted
index is in red.

It would be of some use to see how the simple linear prediction may depend on
some factors, though it cannot be said successful and hence a quantitative
investigation does not make sense. We hope to learn some experience for the
machine learning to be performed soon. First, we check the effect of
filtering. As shown in Supplementary Figure S5, without filtering, the result
is much noisy, the 2015/2016 event is too strong, and the strong 2009/2010
event completely disappears; with low-pass filtering (only the signals below a
year are filtered), the result is similar, only with the noise reduced and the
2015/2016 event better.

Second, we perform an empirical orthogonal function (EOF) analysis for the 336
time-delayed SSN series, which forms a column vector with 336 entries at each
time step: [${\rm SSN}(n-600),{\rm SSN}(n-599),...,{\rm SSN}(n-265)$]T. The
variance for each EOF mode is plotted in Supplementary Figure S6. From it the
first 8, 25, and 50 principal components (PCs) approximately account for 84%,
89%, and 92% of the total variances, respectively. Some examples of the EOF
modes are plotted in Supplementary Figure S7. Now use the first 25 PCs as
inputs and repeat the above process to make the projection. The result is
shown in Fig. 2b. Obviously, during the 12 years of prediction from 2005 to
2017, except for the strong 2009/2010 event, the others are generally fine.
Projections with other numbers of PCs have also been conducted. Shown in
Supplementary Figure S8 are examples with 8 and 50 PCs, respectively. It is
particularly interesting to see that, except for the 2009/2010 event, the case
with 8 PCs already captures the major trend of the EMI evolution. Recall that,
in performing the information flow analysis, we have found that the
reconstructed dynamical system approximately has a dimension of 8. Here the
prediction agrees with the dimension inference.

Because of the encouraging linear model result, it is desirable to achieve a
better projection using more sophisticated tools in deep learning, e.g., the
back propagation (BP) neural network algorithmGoodfellow2016 . (At the
preparation of the original version of this manuscript, we noticed that
recently there have been applications of other machine learning methods to El
Niño forecasts, e.g., Ham2019 .) Three hidden layers are used for the BP
neural network forecast. Again, a major issue here is the short observational
period which may prevent from an appropriate training with many weights. We
hence use only 8, 6, and 1 neurons for the first, second, and third hidden
layers, respectively, as schematized in Fig. 3. This architecture will be
justified soon. To predict an EMI at step (month) $n$, written $y$, a
$336\times 1$ vector ${\bf x}$ is formed with SSN data at steps (months)
$n-600,n-599,...,n-265$, the same as the case with lead times of 50-22 years
for the above simple linear regression model. But even with the modest number
of parameters, the data are still not enough. To maximize the use of the very
limited data, we choose the first 25 principal components (PCs) that account
for 90% of the total variance. We hence train the model with these 25 inputs
(rather than 336 inputs); that is to say, ${\bf x}$ now is a $25\times 1$
vector, which greatly reduces the number of parameters to train. (See below in
Eq. (1): The ${\bf w}_{1}$ is reduced from a $8\times 336$ matrix to a
$8\times 25$ matrix.) Once this is done, ${\bf x}$ is input into the following
equation to arrive at the prediction of $y$:

$\displaystyle y={{\bf w}_{4}}\ {\rm tansig}({{\bf w}_{3}}\ {\rm tansig}({\bf
w}_{2}\ {\rm tansig}({\bf w}_{1}\ {\bf x}+{\bf b}_{1})+{\bf b}_{2})+{\bf
b}_{3})+{\bf b}_{4},$ (1)

where the matrices/vectors of parameters ${\bf w}_{1}$ ($8\times 25$ matrix),
${\bf w}_{2}$ ($6\times 8$ matrix), ${\bf w}_{3}$ ($1\times 6$ row vector),
${\bf w}_{4}$ ($1\times 1$ matrix—a scalar), ${\bf b}_{1}$ ($8\times 1$
vector), ${\bf b}_{2}$ ($6\times 1$ vector), $b_{3}$ (scalar) and $b_{4}$
(scalar) are obtained through training. An illustrative explanation is seen in
Fig. 3. As has been argued (e.g. Heaton2008 ) that the number of hidden
neurons should be kept fewer than 2/3 the size of the input and output layers,
i.e., $\frac{2}{3}\times(25+1)=17$ here. With this architecture we have 15
neurons in total, meeting the requirement.

Figure 3: The neural network architecture used in this study, which is made of
an input layer, an output layer, and three hidden layers. The symbols are the
same as those in (1). The matrices/vectors ${\bf w}_{i}$ and ${\bf b}_{i}$
($i=1,2,3,4$) are obtained through training. To predict the El Niño Modoki
Index (EMI) at month/step $n$, ${\rm EMI}(n)$, find the first 25 principal
components of [${\rm SSN}(n-600),{\rm SSN}(n-599),...,{\rm SSN}(n-265)$]T, a
column vector with 336 entries, and form the input vector ${\bf x}$. The final
output $y$ is the predicted ${\rm EMI}(n)$. Iterating on $n$, we can get a
prediction of EMI for any target time interval.

To predict the EMI at month/step $n$, ${\rm EMI}(n)$, organize ${\rm
SSN}(n-600),{\rm SSN}(n-599),...,{\rm SSN}(n-265)$, that is, the SSNs lagged
by 50 through 22 years, into a $336\times 1$ vector. Project it onto the first
25 EOFs as obtained above, and obtain ${\bf x}$, a vector with 25 entries.
Input ${\bf x}$ into Eq. (1). The final output $y$ is the predicted ${\rm
EMI}(n)$. Iterating on $n$, we can get a prediction of EMI for any target time
interval.

The 5-year data prior to the starting step of prediction are used for
validation. In this study, we start off predicting EMI at January 2008. We
hence take the data over the period from January 2003 through December 2007 to
form the validation set, and the data until December 2002 to train the model.
10,000 runs have been performed, and we pick the one that minimizes the mean
square error (MSE) over the validation set. This is plotted in Fig. 4.
Remarkably, the 12-year long index has been forecast with high skill (corr.
coeff.=0.91). Particularly, the strong 2009/10 event is well predicted, so is
the 2019/20 event. The 2014/16 event, which has made the El Niño forecasts off
the mark in 2014-15, appears a little weaker but also looks fine by trend.

Figure 4: Projection of the El Niño Modoki index based solely on the sunspot
numbers 22-50 years ago with the back propagation neural network algorithm.
The EMI is in blue, while the predicted index is in red. Lightly shaded is the
period over which the validation set is formed, and the period with dark
shading is for prediction.

Out of the 10,000 runs we select the top 10% with lowest MSE over the
validation set (January 2003-December 2007) to form an ensemble of predictions
for the 12 years of EMI (January 2008 till now). The spread, the mean, and the
standard deviation of the ensemble are shown in Fig. 5. Also overlaid is the
observed EMI (blue). As can be seen, the uncertainty is rather limited for a
climate projection.

Figure 5: The 1000 predictions showing the spread of the ensemble. Overlaid
are the observed EMI (blue), the mean of the realizations (cyan), and the mean
plus and minus the standard deviation (black). The light shading marks the
period for validation, while the darker shading marks the prediction period.

We have examined how the projection performance may vary with the choices of
PCs, the number of neurons, and the lags. The PCs of a vector are obtained by
projecting it onto the EOF modes (described above), which possess variances
and structures as plotted in Supplementary Figures S6 and S7, respectively.
Shown in Fig. 6 is the MSE as a function of the number of PCs. Obviously with
25 PCs the MSE is optimized, and that is what we have chosen for the network
architecture. But if one takes a closer look, the MSE does not decrease much
beyond 8. Recall that, in performing the causality analysis, we have seen that
the reconstructed dynamical system approximately has a dimension 8. Fig. 6
hence provides a verification of the dimension inference.

Figure 6: Mean square error for the neural network prediction as a function of
number of principal components.

The EOF modes associated with these 8 PCs are plotted in Supplementary Figure
S7. As can be seen, except for the 11-year process, among others, the 5-year
variability is essential. This is, again, consistent with what we have done in
forming the embedding coordinates: select the delayed series every 5 years.

The dependence of the model performance on the number of neurons has also been
investigated. We choose 1 neuron for layer 3, leaving the numbers of layers 1
and 2 for tuning. The result is contoured in Fig. 7. As can be clearly seen, a
minimum of MSE appears at (8,6), i.e., when the numbers of hidden layer 1 and
layer 2 are, respectively, 8 and 6. This is what we have chosen to launch the
standard prediction.

Figure 7: Mean square error as function of number of neurons.

Also studied is the influence of the lags for the delayed series. In Fig. 8
the MSE is shown as function of the lower and upper bounds of the delays (in
years). From the figure there are a variety of local minima, but the one at
(22, 50) is the smallest one. That is to say, the series with delays from 22
to 50 years make the optimal embedding coordinates, and this is just we are
choosing.

Figure 8: Mean square error as function of lags.

It is known that machine learning suffers from the problem of
reproducibilityHutson2018 . But in this study, as we have shown above, the
prediction guided by the information flow analysis is fairly robust, with
uncertainty acceptably small for climate projection. For all that account, at
least the El Niño Modoki events so far are mostly predictable at a lead time
of more than 10 years.

Concluding remarks. The rigorously developed theory of information flow in
terms of predictability/uncertainty transfer provides a natural way for one to
seek for predictor(s) for a dynamical phenomenon. In this study, it is found
that the delayed causal pattern from SSN to the Pacific SST resembles very
much the El Niño Modoki mode; the former hence can be used to predict the
latter. Indeed, as detailed above, with all the observations we have had so
far, El Niño Modoki can be essentially predicted based solely on SSN at a lead
time of 12 years or over. Particularly, the strong event in 2009/10 and the
elusive event during 2014-16 have been well predicted (Fig. 2b).

We, however, do NOT claim that El Niño Modoki is ultimately driven by solar
activities. It is NOT our intention in this study to investigate on the
dynamical aspects of this climate mode. We just present an observational fact
on fulfilled predictions. There is still a very long way to go in unraveling
the dynamical origin(s) of El Niño Modoki. The success of atmosphere-ocean
coupled models during the past decades confirms that the canonical El Niño is
an intrinsic mode in the climate system; El Niño Modoki may be so also.
Nonetheless, despite the long-standing controversy on the role of solar
activity in climate change (see a review in Bard2006 ), there does exist
evidence on the lagged response as identified here; the North Atlantic climate
response is such an exampleScaife2013 .

As a final remark, it is interesting to see that the maximum information flow
from the SSN lagged by 45 years to the Pacific SST agrees very well with
frequent occurrences of El Niño Modoki during the period of 2000-2010: three
of the four El Niño events are of this type (i.e., the 02/03, 04/05, 09/10
episodes). If we go back to 45 years ago, the sun is most active during the
period of 1955-1965 (Fig. 9) Particularly, in the wavelet spectrum, there is a
distinct high at the scale level of 4 (corresponding to a time scale of
$2^{-4-(-15)}=2048$ days $\approx 5.6$ years)—Recall that this is roughly the
sampling interval we used in choosing the time-delay series to form the
${\nu}$-dimensional system and compute the information flow. After that
period, the sun becomes quiet, correspondingly we do not see much El Niño
Modoki for the decade 2010-20. But two peaks of SSN are seen in the twenty
years since 1975-78. Does this correspondence herald that, after entering
2020, in the following two decades, El Niño Modoki may frequent the equatorial
Pacific again? We don’t know whether this indeed points to a link, but when
this project was initiated early in 2018, the projected El Niño Modoki and La
Niña Modoki already occurred recently. This question, among others, are to be
addressed in future studies. (For reference, the codes used in this study are
available at
http://www.ncoads.org/enso_modoki_data_codes.tar.gz. )

Figure 9: Daily series and orthonormal wavelet spectrum of the sunspot numbers
based on the spline wavelets built in LA2007 .

## METHOD

Estimation of information flow and causality among multivariate time series.
Information flow is a fundamental notion in general physics which has
applications in a wide variety of scientific disciplines (see Liang2016  and
references therein). Its importance lies beyond the literal meaning in that it
implies causation, uncertainty propagation, predictability transfer, etc.
Though with a history of research for more than 30 years, it has just been put
on a rigorous footing and derived from first principles, initially motivated
by the predictability study in atmosphere-ocean scienceLK2005 . Since its
birth it has been validated with many benchmark dynamical systems such as
baker transformation, Hénon map, Rössler system, etc. (cf. Liang2016 Liang2014
), and has been applied with success to different problems in earth system
science, neuroscience and quantitative finance, e.g., Vannitsem2019 ,
Dionisis2019 , Stips2016 , to name a few. Hereafter we just give a very brief
introduction.

Consider a dynamical system

$\displaystyle\frac{d{\bf x}}{dt}={\bf F}(t;{\bf x})+{\bf B}(t;{\bf
x})\dot{{\bf w}},$ (2)

where ${\bf x}$ and ${\bf F}$ are ${\nu}$-dimensional vector, ${\bf B}$ is a
${\nu}$-by-$m$ matrix, and $\dot{{\bf w}}$ an $m$-vector of white noise. Here
we follow the convention in physics and do not distinguish the notation of a
random variable from that of a deterministic variable. Liang (2016)Liang2016
established that the rate of uncertainty (in terms of Shannon entropy)
transferred from $x_{2}$ to $x_{1}$, or rate of information flowing from
$x_{2}$ to $x_{1}$ (in nats per unit time), is:

$\displaystyle T_{2\to
1}=-E\left[\frac{1}{\rho_{1}}\int_{{}^{{\nu}-2}}\frac{\partial(F_{1}\rho_{{\backslash
2}})}{\partial
x_{1}}dx_{3}...dx_{\nu}\right]+\frac{1}{2}E\left[\frac{1}{\rho_{1}}\int_{R^{{\nu}-2}}\frac{\partial^{2}(g_{11}\rho_{{\backslash
2}})}{\partial x_{1}^{2}}dx_{3}...dx_{\nu}\right],$ (3)

where $E$ stands for mathematical expectation, $\rho_{1}=\rho_{1}(x_{1})$ is
the marginal probability density function (pdf) of $x_{1}$, $\rho_{{\backslash
2}}=\int\rho dx_{2}$, and $g_{11}=\sum_{j=1}^{m}b_{1j}b_{1j}$. The units are
in nats per unit time. Ideally if $T_{2\to 1}=0$, then $x_{2}$ is not causal
to $x_{1}$; otherwise it is causal (for either positive or negative
information flow). But in practice significance test is needed.

Generally $T_{2\to 1}$ depends on $(x_{3},...,x_{\nu})$ as well as
$(x_{1},x_{2})$. But it has been established thatLiang2018  it is invariant
upon arbitrary nonlinear transformation of $(x_{3},x_{4},...,x_{\nu})$,
indicating that information flow is an intrinsic physical property. Also
established is the principle of nil causality, which asserts that, if the
evolution of $x_{1}$ is independent of $x_{2}$, then $T_{2\to 1}=0$. This is a
quantitative fact that all causality analyses try to verify in applications,
while in the framework of information flow it is a proven theorem.

In the case of linear systems where ${\bf F}({\bf x},t)={\bf A}{\bf x}$, ${\bf
A}=(a_{ij})$, the formula (3) becomes quite simpleLiang2016 :

$\displaystyle T_{2\to 1}=a_{12}\frac{\sigma_{12}}{\sigma_{11}},$ (4)

where $\sigma_{ij}$ is the population covariance. An immediate corollary is
that causation implies correlation, while correlation does not imply
causation, fixing the long-standing philosophical debate over causation versus
correlation ever since Berkeley (1710)Berkeley1710 .

To arrive at a practically applicable formula, we need to estimate (4), given
${\nu}$ time series. The estimation roughly follows that of Liang
(2014)Liang2014 , and can be found in Liang2021 . The following is just a
brief summary of some relevant pieces as detailed in Liang2021 .

Suppose we have ${\nu}$ time series, $x_{i}$, $i=1,...,{\nu}$, and these
series are equi-spaced, all having $N$ data points $x_{i}(n)$, $n=1,2,...,N$.
Assume a linear model ${\bf F}({\bf x};t)={\bf A}{\bf x}$, ${\bf A}=(a_{ij})$
being a ${\nu}\times{\nu}$ matrix), and ${\bf B}$ being a ${\nu}\times{\nu}$
diagonal matrix. To estimate $T_{2\to 1}$, we first need to estimate $a_{12}$.
As shown before in Liang2014 , when $a_{ij}$ and $b_{i}$ are constant, the
maximal likelihood estimator (mle) is precisely the least square solution of
the following $N$ (overdetermined) algebraic equations

$\displaystyle\sum_{j=1}^{\nu}a_{1j}x_{j}(n)=\dot{x}_{1}(n),\qquad\qquad
n=1,...,N$ (5)

where $\dot{x}_{1}(n)=(x_{1}(n+1)-x_{1}(n))/\Delta t$ is the differencing
approximation of $dx_{1}/dt$ using the Euler forward scheme, and $\Delta t$ is
the time stepsize (not essential; only affect the units). Following the
procedure in Liang2014 , the least square solution of
$(a_{11},...,a_{1{\nu}})$, $(\hat{a}_{11},...,\hat{a}_{1{\nu}})$, satisfies
the algebraic equation

$\displaystyle\left[\begin{array}[]{ccc}C_{1,1}&...&C_{1,{\nu}}\\\
\vdots&\vdots&\vdots\\\
C_{{\nu},1}&...&C_{{\nu},{\nu}}\end{array}\right]\left[\begin{array}[]{c}\hat{a}_{11}\\\
\vdots\\\
\hat{a}_{1{\nu}}\end{array}\right]=\left[\begin{array}[]{c}C_{1,d1}\\\
\vdots\\\ C_{{\nu},d1}\end{array}\right],$

where

$\displaystyle
C_{ij}=\frac{1}{N}\sum_{n=1}^{N}(x_{i}(n)-\bar{x}_{i})(x_{j}(n)-\bar{x}_{j})$
$\displaystyle
C_{i,dj}=\frac{1}{N}\sum_{n=1}^{N}(x_{i}(n)-\bar{x}_{i})(\dot{x}_{j}(n)-{\bar{\dot{x}}}_{j})$

are the sample covariances. Hence $\hat{a}_{12}=\frac{1}{\det{\bf
C}}\cdot\sum_{j=1}^{\nu}\Delta_{2j}C_{j,d1}$ where $\Delta_{ij}$ are the
cofactors. This yields an estimator of the information flow from $x_{2}$ to
$x_{1}$:

$\displaystyle\hat{T}_{2\to 1}=\frac{1}{\det{\bf
C}}\cdot\sum_{j=1}^{\nu}\Delta_{2j}C_{j,d1}\cdot\frac{C_{12}}{C_{11}},$ (7)

i.e., Eq. (7) in section 2. Here $\det{\bf C}$ is the determinant of the
covariance matrix ${\bf C}$, and now $T$ is understood as the MLE of $T$. We
slightly abuse notation for the sake of simplicity.) For two-dimensional (2D)
systems, ${\nu}=2$, the equation reduces to

$\displaystyle\hat{T}_{2\to
1}=\frac{C_{11}C_{12}C_{2,d1}-C_{12}^{2}C_{1,d1}}{C_{11}^{2}C_{22}-C_{11}C_{12}^{2}},$
(8)

recovering the familiar one as obtained in Liang2014  and frequently used in
applications (e.g., Stips2016 , Vannitsem2019 , Dionisis2019 ).

The significance of $\hat{T}_{2\to 1}$ of (7) can be tested following the same
strategy as that used in Liang2014 . By the MLE property (cf. Garthwaite1995 ,
when $N$ is large, $\hat{T}_{2\to 1}$ approximately follows a Gaussian around
its true value with a variance
$\left(\frac{C_{12}}{C_{11}}\right)^{2}{\hat{\sigma}}_{a_{12}}^{2}$. Here
${\hat{\sigma}}_{a_{12}}^{2}$ is the variance of $\hat{a}_{12}$, which is
estimated as follows. Denote by ${\bm{\theta}}$ the vector of parameters to be
estimated; here it is $(a_{11},a_{12},a_{13},...,a_{1{\nu}};b_{1})$. Compute
the Fisher information matrix ${\bf I}=(I_{ij})$

$\displaystyle
I_{ij}=-\frac{1}{N}\sum_{n=1}^{N}\frac{\partial^{2}\log\rho({\bf x}_{n+1}|{\bf
x}_{n};\hat{{\bm{\theta}}})}{\partial\theta_{i}\partial\theta_{j}}$

where $\rho({\bf x}_{n+1}|{\bf x}_{n})$ is a Gaussian for a linear model:

$\displaystyle\rho({\bf x}_{n+1}|{\bf x}_{n})$
$\displaystyle=\frac{1}{[(2\pi)^{2}\det({\bf B}{\bf B}^{T}\Delta t)]^{1/2}}$
$\displaystyle\qquad\times e^{-\frac{1}{2}({\bf x}_{n+1}-{\bf x}_{n}-{\bf
A}{\bf x}_{n}\Delta t)^{T}({\bf B}{\bf B}^{T}\Delta t)^{-1}({\bf x}_{n+1}-{\bf
x}_{n}-{\bf A}{\bf x}_{n}\Delta t)}.$

The computed entries are then evaluated using the estimated parameters, and
hence the Fisher information matrix is obtained. It has been established (cf.
Garthwaite1995 ) that ${\bf I}^{-1}/N$ is the covariance matrix of
$\hat{{\bm{\theta}}}$, from which it is easy to find the entry
$\hat{\sigma}_{a_{12}}^{2}$ [here it is the entry (2,2)]. Given a level, the
significance of an estimated information flow can be tested henceforth.

### Data availability

The data used in this study include the SST and sunspot number (SSN) from
National Oceanic and Atmospheric Administration (NOAA) (SST:
http://www.esrl.noaa.gov/psd/data/gridded/cobe2.html; SSN:
https://wwww.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/Data/sunspot.long.data),
and the El Niño Modoki Index (EMI) from Japan Agency for Marine-Earth Science
and Technology (JAMSTEC)
(http://www.jamstec.go.jp/aplinfo/sintexf/e/elnmodoki/data.html). Among them
SST has a horizontal resolution of $1^{\rm o}\times 1^{\rm o}$. Temporally all
of them are monthly data, and at the time when this study was initiated, SST,
EMI, and SSN respectively have a time coverage of 01/1850-12/2018,
01/1870-11/2018, and 05/1847-12/2017. Daily SSN is also used for spectral
analysis; it covers the period 15/04/1921-31/12/2010.

### Code availability

All the codes (and some of the downloaded data) are available at
http://www.ncoads.org/upload/enso_modoki_data_codes.tar.gz

Acknowledgments. The long discussion with Richard Lindzen on August 17, 2019,
and the email exchange later on in late February 2020 are sincerely
appreciated. The comments from Dmitry Smirnov and four anonymous reviewers
have helped improve the manuscript. We thank NOAA and JAMSTEC for making the
data available, and the NUIST High Performance Computing Center for providing
computational resources. This study is partially funded by National Science
Foundation of China (No. 41975064), and the 2015 Jiangsu Program for
Innovation Research and Entrepreneurship Groups.

Author contributions. XSL: idea, methodology, experiment design, computation
(causality analysis and prediction), writing; FX: experiment design,
computation (causality analysis); YR: experiment design, computation
(prediction). All authors discussed the study results.

Competing interests. The authors declare no competing interests.

## References

  * (1) McPhaden, M.J. (2015), Playing hide and seek with El Niño, Nature Climate Change 5, 791-795.
  * (2) Ashok, K., and T. Yamagata (2009), The El Niño with a difference. Nature 461, 481.
  * (3) Trenberth, K.E., and D.P. Stepaniak (2001), Indices of El Niño evolution. J. Clim. 14, 1697-1701.
  * (4) Larkin, N.K., and D.E. Harrison (2005), Global seasonal temperature and precipitation anomalies during El Niño autumn and winter. Geophys. Res. Lett. 32, L16705.
  * (5) Yu, J.Y., and H.Y. Kao (2007), Decadal changes of ENSO persistence barrier in SST and ocean heat content indices: 1958-2001. J. Geophys. Res. Atmos. D13106b.
  * (6) Lee, T., and M.J. McPhaden (2010), Increasing intensity of El Niño in the central equatorial Pacific. Geophys. Res. Lett. 37, L14603.
  * (7) Yu, J.Y., and S.T. Kim (2010), Three evolution patterns of central-pacific El Niño . Geophys. Res. Lett. 37, L08706.
  * (8) Ashok, K., S.K. Behera, S.A. Rao, H.Y. Weng, and T. Yamagata (2007), El Niño Modoki and its possible teleconnection. J. Geophys. Res. Oceans 112, C11007.
  * (9) Kug, J.S., F.F. Jin, and S.I. An (2009), Two types of El Niño events: cold tongue El Niño and warm pool El Niño . J. Clim. 22, 1499-1515.
  * (10) Fu, C., H.F. Diaz, and J.O. Fletcher (1986), Characteristics of the response of sea surface temperature in the central Pacific associated with warm episodes of the Southern Oscillation. Mon. Weath. Rev. 114, 1716-1738.
  * (11) Wang, C.Z., C. Deser, J.Y. Yu, P. Dinezio, and A. Clement (2017), “El Niño and Southern Oscillation (ENSO): A Review” in Coral Reefs of the Eastern Tropical Pacific, P.W. Glynn et al., Eds., Springer, Netherlands, pp. 85-106.
  * (12) Behera, S., and T. Yamagata (2018), Climate Dynamics of ENSO Modoki Phenomena, Oxford Research Encyclopedias.
  * (13) Mo, K.C. (2010), Interdecadal modulation of the impact of ENSO on precipitation and temperature over the United States. J. Clim. 23, 3639-3656.
  * (14) Barnston, A.G., A. Kumar, L. Goddard, and M.P. Hoerling (2005), Improving seasonal prediction practices through attribution of climate variability. Bull. Amer. Meteorol. Soc. 2005, 59-72.
  * (15) Xie, R., and F.-F. Jin (2018), Two leading ENSO modes and El Niño Types in the Zebiak-Cane Model. J. Clim. 31, 1943-1962.
  * (16) Zebiak, S.E., and M.A. Cane (1987), A model El Niño-Southern Oscillation. Mon. Weath. Rev. 115, 2262-2278.
  * (17) Chen, D., and Coauthors (2015), Strong influence of westerly wind bursts on El Niño diversity. Nature Geoscience DOI:10.1038/NGEO2399.
  * (18) Tang, Y.M., and Coauthous (2018), Progress in ENSO prediction and predictability study. National Science Review 5, 826-839.
  * (19) Von Storch, H. (2001), “Statistics — an indispensable tool in dynamical modeling” in Models in Environmental Research, H. von Storch and G. Flöser, Eds., Springer Verlag.
  * (20) Liang, X.S., and R. Kleeman (2005), Information transfer between dynamical system components. Phys. Rev. Lett. 95, 244101.
  * (21) Liang, X.S. (2014), Unraveling the cause-effect relation between time series. Phys. Rev. E. 90, 052150.
  * (22) Liang, X.S. (2016), Information flow and causality as rigorous notions ab initio. Phys. Rev. E 94, 052201.
  * (23) Liang, X. S. (2018), Causation and information flow with respect to relative entropy. Chaos, 28, 075311.
  * (24) Vannitsem, S., Q. Dalaiden, and H. Goose (2019), Testing for dynamical dependence — Application to the surface mass balance over Antarctica. Geophys. Res. Lett. DOI: 10.1029/2019GL084329.
  * (25) Hristopulos, D.T., A. Babul, S.A. Babul, L.R. Brucar, and N. Virji-Babul (2019), Disrupted information flow in resting-state in adolescents with sports related concussion. Frontiers in Human Neuroscience 13, 419. DOI: 10.3389/fnhum.2019.00419.
  * (26) Stips, A., D. Macias, C. Coughlan, E. Garcia-Gorriz, and X.S. Liang (2016), On the causal structure between CO2 and global temperature. Scientific Reports 6: 21691. DOI: 10.1038/srep21691.
  * (27) Johnson, E.N., and S.K. Kannan (2005), Adaptive trajectory control for autonomous helicopters. Journal of Guidance, Control, and Dynamics 28(3), 524-538.
  * (28) Kim, P., J. Rogers, J. Sun, and E. Bollt (2017), Causation entropy identifies sparsity structure for parameter estimation of dynamic systems. Journal of Computational and Nonlinear Dynamics 12(1).
  * (29) Wallace, J.M., and D.S. Gutzler (1981), Teleconnections in the geopotential height field during the northern hemisphere winter. Mon. Wea. Rev. 109, 784-812.
  * (30) Takens, F. (1981), “Detecting strange attractors in turbulence” in Dynamical Systems and Turbulence, Lecture Notes in Mathematics, D.A. Rand, L.-S. Young, Eds., Springer-Verlag, New York, NY, vol. 898, pp. 366-381.
  * (31) Abarbanel, H.D.I. (1996), Analysis of Observed Chaotic Data, Springer-Verlag, New York.
  * (32) Liang, X.S., and D.G.M. Anderson (2007), Multiscale window transform. SIAM J. Multiscale Model. Simul. 6, No.2, 437-467.
  * (33) Goodfellow, I., Bengio, Y., Courville, A. (2016), Deep Learning. MIT Press.
  * (34) Ham, Y.-G., J.-H. Kim, and J.J. Luo (2019), Deep learning for multi-year ENSO forecasts. Nature. DOI:10.1038/s41586-019-1559-7.
  * (35) Heaton, J. (2008), Introduction to Neural Networks with Java. Heaton Research Inc., 44p pp.
  * (36) Hutson, M. (2018), Artificial intelligence faces reproducibility crisis. Science 359, 725-726.
  * (37) Bard, E., and M. Frank (2006), Climate change and solar variability: What’s new under the sun? Earth Planet. Sci. Lett. 248, 1-14.
  * (38) Scaife, A.A., S. Ineson, J.R. Knight, L. Gray, K. Kodera, and D.M. Smith (2013), A mechanism for lagged North Atlantic climate response to solar variability. Geophys. Res. Lett. 40, 434-439.
  * (39) Berkeley, G. (1710), A Treatise on Principles of Human Knowledge.
  * (40) Liang, X.S., Normalized multivariate time series causality analysis and causal graph reconstruction. The 37th Conference on Uncertainty in Artificial Intelligence, (submitted).
  * (41) Garthwaite, P.H., I.T. Jolliffe, and B. Jones (1995), Statistical Inference. Prentice-Hall, Hertfordshire, UK.

##  El Niño Modoki thus far can be mostly predicted more than 10 years ahead
of time

X. San Liang, Fen Xu, & Yineng Rong

Supplementary Information

Figure S1: The absolute information flow estimated using the bivariate formula
from the delayed series of sunspot numbers (SSN) to those of the sea surface
temperature (SST) in the Pacific Ocean (in nats/month). The SST time series
between January 1980 - December 2017 are used. In each subplot, the lower
panel shows only the value that is significant at a 90% confidence level, and
indicated on it is the delay (in years).

Figure S2:  As Supplementary Figure S1, but the information flows are
estimated using the multivariate formula, with the embedding coordinates
formed with SSN series lagged by 22-50 years every 5 years.  Figure S3:  The
absolute information flow from the SSN series delayed by 45 years to those of
the SST in the Pacific Ocean (in nats/month). It is estimated with the
multivariate formula using embedding coordinates formed with SSN series
delayed by 22-50 years every 5 years. The SST time series between January 1980
- December 2005 are used. In each subplot, the lower panel shows only the
value that is siginificant at a 90% confidence level.  Figure S4:  The
original monthly SSN series (blue), lowpass-filtered series (black), and
bandpass-filtered series (red).

Figure S5:  Projection of the El Niño Modoki Index (EMI) based solely on the
SSN 22-50 years ago with a simple linear regression model. Upper: The SSN
series as input is not filtered. Lower: The SSN series is low-pass filtered.
The observed EMI is in blue, while the projected is in red.  Figure S6:
NormaliEd variance vs. EOF mode. The EOF analysis is performed with the SSN
time series with delays from 22-50 years, which form a vector with 336 members
of series.  Figure S7:  Examples of the EOF modes as described in
Supplementary Figure S6.

Figure S8:  As Figure 2b, but the projection uses, respectively, the first 8
and 50 principal components as inputs.


# Control of spontaneous emission dynamics in microcavities with chiral
exceptional surfaces

Q. Zhong,1 A. Hashemi,1 Ş. K. Özdemir,2,3 and R. El-Ganainy1,4,
<EMAIL_ADDRESS>1Department of Physics, Michigan Technological University,
Houghton, Michigan 49931, USA  2Department of Engineering Science and
Mechanics, The Pennsylvania State University, University Park, Pennsylvania
16802, USA  3Materials Research Institute, The Pennsylvania State University,
University Park, Pennsylvania 16802, USA  4Henes Center for Quantum Phenomena,
Michigan Technological University, Houghton, Michigan 49931, USA

###### Abstract

We investigate spontaneous emission from a quantum emitter located within the
mode volume of a microring resonator that features chiral exceptional points.
We show that this configuration offers enough degrees of freedom to exhibit a
full control to either enhance or suppress the emission process. Particularly,
we demonstrate that the Purcell factor can be enhanced by a factor of two
beyond its value in an identical microring operating at a diabolic point. Our
conclusions, which are derived using a non-Hermitian Hamiltonian formalism,
are confirmed by employing full-wave simulations of realistic photonic
structures and materials. Our results offer a straightforward route to improve
the performance of single photon sources using current photonics technology
without the need for building optical resonators with ultrahigh quality
factors or nanoscale volumes.

##  I Introduction

Quantum engineering seeks to utilize quantum mechanics to build a new
generation of computing machines, encryption schemes, and sensing devices with
unprecedented performance in terms of computational power, security strength,
and sensitivity, among other applications. In the pursuit to achieve these
goals, several material platforms provide complementary solutions to overcome
various practical hurdles. These include trapped atoms Gross and Bloch (2017),
superconducting circuits Blais  _et al._ (2020), and photonics O’Brien (2007).
The latter is particularly interesting due to its mature technology and
natural interface with current optical communication systems. At the heart of
modern quantum optics technology is the ability to control light-matter
interaction at the quantum level for various applications such as building
nonclassical light sources Zeilinger (2017), optical transistors Chang  _et
al._ (2007), and quantum memory Lvovsky  _et al._ (2009). In this regard,
efforts have been recently dedicated to building efficient single photon
sources that can produce individual photons on demand at high repetition rates
Buckley  _et al._ (2012). This progress was enabled by engineering various
optical resonator geometries that can support small modal volumes and large
quality factors to tailor the photonic local density of states (PLDOS)
surrounding quantum emitters (QE), hence controlling their spontaneous
emission (SE) rates as quantified by the Purcell factor (PF) Purcell (1946)
(see Ref. Yokoyama and Ujihara (2019) for detailed discussions). Examples
include planar photonic crystals, vertical Bragg reflectors, microdisks and
plasmonic structures. For comprehensive reviews and performance comparison,
see Ref. Vahala (2003). Despite these promising results, the aforementioned
arrangements are not easy to mass produce or integrate with other photonics
components. An attractive alternative in terms of mass fabrication and large-
scale integration is microring resonators Vahala (2004). On the downside,
however, microring resonators suffer from relatively large mode volumes and
limited quality factors Vahala (2003). It will be thus of interest to devise
new routes for controlling SE in microring resonators and possibly enhance
their PF beyond their current performance.

Motivated by this goal, here, we study the interaction between light and a QE
in a family of microring resonators whose design is tailored to operate in the
vicinity of or at a special type of non-Hermitian singularities known as
chiral exceptional points (EPs). At an EP, two or more eigenmodes of a system
fuse together, leading to a reduction of the dimensionality associated with
the eigenspace of the system (for recent reviews on the physics of non-
Hermitian optical systems and EPs see Refs. El-Ganainy  _et al._ (2018); Feng
_et al._ (2017); Özdemir  _et al._ (2019); Miri and Alù (2019)). In contrast
to recent investigations of SE near isolated EPs Lin  _et al._ (2016); Pick
_et al._ (2017a, b); Hadad and Engheta (2020); Khanbekyan and Wiersig (2020),
here instead, we consider a system that features an exceptional surface Zhong
_et al._ (2019a, b, 2020) which offers robust operation and tunability. As we
will see, this geometry offers several advantages in terms of controlling PF
[by suppressing it completely or enhancing its value by a factor of two
compared with microring resonators operating at diabolic points (DPs)] and
integration with a waveguide to collect the emitted photons from a
predetermined port.

##  II Formalism

The optical platform we consider in this study is shown in Fig. 1. It consists
of a microring resonator coupled to a waveguide terminated by a mirror at one
of its ports. In the absence of the mirror, the optical resonator has two
degenerate modes, i.e., it operates at a DP. On the other hand, the mirror
provides unidirectional coupling between the clockwise (CW) and
counterclockwise (CCW) modes and, hence, induces a chiral EP (in fact an
exceptional surface). This structure was recently proposed as a robust
implementation for EP-based optical sensors Zhong  _et al._ (2019a), optical
amplifiers Zhong  _et al._ (2020), and directional absorbers Zhong  _et al._
(2019b). In these previous papers, the system was externally excited through
the waveguide, and the response was studied by monitoring the transmission or
reflection spectra under different conditions. Here, we study the emission
properties of a QE located within the mode volume of the resonator (Fig. 1) by
monitoring the emitted optical power collected by the waveguide ports. We are
particularly interested in evaluating how the presence of this chiral EP
modifies the PF of the QE, i.e., the decay rate from the excited state to the
ground state in the absence of any external active perturbation (no applied
optical or electrical signal).

Figure 1: A schematic of the proposed geometry. It consists of a microring
resonator evanescently coupled to a waveguide with an end mirror as proposed
in Ref. Zhong  _et al._ (2019a) for implementing exceptional surfaces. Here,
we assume that a quantum emitter (QE) is located inside the microring
resonator. Under resonant conditions (transition frequency between the QE
energy levels matches the eigenfrequencies of the optical modes of the ring
resonator), the QE couples to the clockwise (CW) and counterclockwise (CCW)
modes of the ring. Photons emitted from the QE take different paths and self-
interfere before they arrive at the exit port (P1 in the figure). This can
lead to the enhancement or suppression of Purcell factor (PF).

To simplify the analysis and gain insight into the problem, we will employ a
technique based on a non-Hermitian Hamiltonian
$\hat{H}=\hat{H}_{R}+\hat{H}_{E}+\hat{H}_{I}$ (with the sub-indices $R,E$, and
$I$ referring to radiation, emitter and interaction, respectively) and neglect
the effect of quantum jumps. The system shown in Fig. 1 can be represented by
the Hamiltonian

$\displaystyle\hat{H}_{R}$
$\displaystyle=\hbar(\omega_{o}-i\gamma_{W})(\hat{a}^{\dagger}\hat{a}+\hat{b}^{\dagger}\hat{b})+\hbar\kappa\hat{a}\hat{b}^{\dagger},$
(1) $\displaystyle\hat{H}_{E}$
$\displaystyle=\hbar(\omega_{e}-i\gamma_{e})\ket{e}\bra{e},$
$\displaystyle\hat{H}_{I}$ $\displaystyle=\hbar(J_{\text{CW
}}\hat{a}^{\dagger}+J_{\text{CCW}}\hat{b}^{\dagger})\ket{g}\bra{e}+\text{H.c.}$

In the radiation Hamiltonian $\hat{H}_{R}$, $\hat{a}^{\dagger}$, $\hat{a}$,
and $\hat{b}^{\dagger}$, $\hat{b}$ are the creation and annihilation operators
of the CW and CCW optical modes, which have the same resonant frequency
$\omega_{o}$ and loss rate $\gamma_{W}$ (coupling loss). The unidirectional
coupling from the CW to the CCW mode is denoted by
$\kappa=-2i\gamma_{W}|r|e^{i\phi_{W}}$ Zhong  _et al._ (2019a); Vahala (2004).
Here, $|r|$ is the field reflection amplitude from the mirror, and
$\phi_{W}=2\beta L+\phi_{r}$, where $\beta$ is the propagation constant in the
waveguide, $L$ is the distance between the waveguide-ring junction and the
mirror, and $\phi_{r}$ is the field reflection phase from the mirror (i.e.,
$r=|r|e^{i\phi_{r}}$). Note that $\hat{H}_{R}$ (which can be inferred from Eq.
(1) in Zhong  _et al._ (2019a) by elevating the classical field variables into
operators) is not Hermitian. In general, the Hamiltonian $\hat{H}_{R}$ will
exhibit an EP of order two in the single-photon subspace. In the emitter
Hamiltonian $\hat{H}_{E}$, $\omega_{e}$ denotes the transition frequency
between the ground state $\ket{g}$ and excited state $\ket{e}$, and
$\gamma_{e}$ quantifies the excited state lifetime due to coupling only to
free space continuum modes. We assume that the emitter does not radiate
efficiently into free space directly, which we will justify later when
discussing implementations. In the interaction Hamiltonian $\hat{H}_{I}$,
$J_{\text{CW,CCW}}=\vec{\mu}\cdot\vec{E}_{\text{CW,CCW}}(\vec{r})$ are the
coupling constants between the CW and CCW optical modes and an emitter having
an electric dipole moment $\vec{\mu}$ and located at $\vec{r}$, and
$\vec{E}_{\text{CW,CCW}}(\vec{r})$ are the normalized electric fields of CW
and CCW modes of the bare microring resonator (i.e., without mirror or QE),
which we use as bases for expanding the optical fields at the location of the
emitter $\vec{r}$ Santori  _et al._ (2010). For these CW/CCW modes, the
magnitude of the electric field at a particular transverse position does not
vary along the angular direction. Instead, the field just acquires a phase. If
we chose the ring-waveguide junction as a reference point, we can then write
$J_{\text{CW}}=Je^{-i\phi_{E}}$, and $J_{\text{CCW}}=Je^{i\phi_{E}}$, where
$\phi_{E}=\beta D$ (see Fig. 1), and H.c. stands for Hermitian conjugate.

Within the single excitation subspace, which is relevant to the SE process,
the general wavefunction can be written as
$\ket{\psi(t)}=a(t)\ket{1,0,g}+b(t)\ket{0,1,g}+c(t)\ket{0,0,e}$, where the
coefficients $a(t)$, $b(t)$, and $c(t)$ are the probability amplitudes of
finding the excitation either in the CW mode, CCW mode, or in the QE.
Importantly, it is straightforward to show that $[\hat{N},\hat{H}]=0$, where
$\hat{N}=\hat{a}^{\dagger}\hat{a}+\hat{b}^{\dagger}\hat{b}+\ket{e}\bra{e}$ is
a generalized number operator that accounts for the total excitations in the
bosonic modes and the QE. In other words, the Hamiltonian $\hat{H}$ conserves
the number of excitations. However, as a result of its non-Hermitian
character, it does not conserve the probability of finding the excitation
trapped inside the system, i.e. $p(t)=|a(t)|^{2}+|b(t)|^{2}+|c(t)|^{2}\leq 1$.
By substituting $\psi(t)$ in the Schrödinger’s equations $i\hbar
d\ket{\psi(t)}/dt=\hat{H}\ket{\psi(t)}$ and projecting on the states
$\ket{1,0,g}$, $\ket{0,1,g}$, and $\ket{0,0,e}$, we obtain

$i\frac{d\vec{v}}{dt}=H\vec{v},H=\begin{pmatrix}\omega_{o}-i\gamma_{W}&0&Je^{-i\phi_{E}}\\\
\kappa&\omega_{o}-i\gamma_{W}&Je^{i\phi_{E}}\\\
Je^{i\phi_{E}}&Je^{-i\phi_{E}}&\omega_{e}-i\gamma_{e}\end{pmatrix},$ (2)

where $\vec{v}=\left(a(t),b(t),c(t)\right)^{T}$, with the superscript $T$
indicating matrix transpose. Note that, consistent with the expression for
$\hat{H}$, the effective discrete Hamiltonian $H$ is non-Hermitian and
exhibits an exceptional surface when $J=0$. To study SE using Eq. (2), we
consider the initial condition $a(0)=b(0)=0$ and $c(0)=1$, i.e., the emitter
is initially in the excited state, and the optical modes are in the vacuum
state.

##  III Results

To ensure that the system operates in the vicinity of the EP, we consider the
weak coupling regime $J\ll\gamma_{W}$ with negligible Markovian effects. In
other words, a photon emitted from the QE into the photonic mode will escape
quickly to the waveguide environment before it is able to couple back to the
QE. In microring resonators with large mode volume, this assumption is valid.
In addition, for resonators having quality factors in the order of $Q\sim
5000$, the photon lifetime is $2Q/\omega_{o}\sim 10^{-1}\>\text{ns}$.
Meanwhile, the typical lifetime of the excited state associated with a quantum
dot in free space is $\sim$ 1–10 ns. In this regime of operation, one can
approximate the decay into the two channels as independent processes. This is
mainly justified by the absence of interference between these two processes.
In Appendix A, we discuss the validity of this approximation.

Under these conditions, we can integrate the last equation in Eq. (2) and
obtain $c(t)\sim e^{-i\omega_{e}t}$. By substituting back in the first two
equations (which is equivalent to performing adiabatic elimination), we obtain

$i\frac{d}{dt}\begin{pmatrix}a\\\
b\end{pmatrix}=\begin{pmatrix}\omega_{o}-i\gamma_{W}&0\\\
\kappa&\omega_{o}-i\gamma_{W}\\\ \end{pmatrix}\begin{pmatrix}a\\\
b\end{pmatrix}+J\begin{pmatrix}e^{-i(\omega_{e}t+\phi_{E})}\\\
e^{-i(\omega_{e}t-\phi_{E})}\end{pmatrix}.$ (3)

By seeking a solution of the form
$\left(a(t),b(t)\right)^{T}=\left(A,B\right)^{T}e^{-i\omega_{e}t}$, we find
the steady state solution as

$\displaystyle A=$ $\displaystyle\frac{Je^{-i\phi_{E}}}{\Delta+i\gamma_{W}},$
(4) $\displaystyle B=$
$\displaystyle\frac{Je^{i\phi_{E}}}{\Delta+i\gamma_{W}}+\frac{J\kappa
e^{-i\phi_{E}}}{\left(\Delta+i\gamma_{W}\right)^{2}},$

where $\Delta=\omega_{e}-\omega_{o}$. As expected, the $\kappa$ term, which is
responsible for the existence of the exceptional surface when $J=0$,
introduces a second-order pole in the expression for $B$. From the above
expressions, we observe that, for $J\ll\gamma_{W}$, $|A|\&|B|\ll 1$ (even when
$\Delta=0$), which justifies the weak coupling approximation used to arrive at
these results. From a classical perspective, one can think of $A$ and $B$ as
the steady state field amplitudes of the CW and CCW modes under excitation by
a driven classical dipole antenna. The output signal from ports 1 and 2 will
be thus
$s_{\text{EP}}^{(1)}=-\sqrt{2\gamma_{W}}\left(A|r|e^{i\phi_{W}}+B\right)$ and
$s_{\text{EP}}^{(2)}=-\sqrt{2\gamma_{W}}\left[|t|e^{i(\beta
L+\phi_{t})}A\right]$, where $\phi_{t}$ is the phase associated with the field
transmission coefficient $t$ from the mirror, i.e, $t=|t|e^{i\phi_{t}}$. The
total normalized power emitted from both ports
$P_{\text{EP}}\equiv|s_{\text{EP}}^{(1)}|^{2}+|s_{\text{EP}}^{(2)}|^{2}$ is
given by

${P_{\text{EP}}}=P_{\text{DP}}\frac{|\chi+|r|+e^{i\Delta\phi}|^{2}+|t|^{2}}{2},$
(5)

where the $P_{\text{DP}}\equiv
P_{\text{EP}}|_{r=0}=\frac{4\gamma_{W}J^{2}}{\Delta^{2}+\gamma_{W}^{2}}$ is
the power emitted in the DP case. In Eq. (5), $\Delta\phi\equiv
2\phi_{E}-\phi_{W}$ and $\chi\equiv\frac{\kappa
e^{-i\phi_{W}}}{\Delta+i\gamma_{W}}=\frac{-2i\gamma_{W}|r|}{\Delta+i\gamma_{W}}$.
Note that $P_{\text{EP,DP}}$ are the power emitted only into the waveguide and
do not account for coupling to free space. By taking this latter effect into
consideration as an additive process and denoting the power decay PF
enhancement as $\eta$, we finally arrive at Vuckovic  _et al._ (1999)

$\eta\equiv\frac{PF_{\text{EP}}}{PF_{\text{DP}}}=\frac{2\gamma_{e}+P_{\text{EP}}}{2\gamma_{e}+P_{\text{DP}}}=\frac{\mathcal{F}+\mathcal{R}}{\mathcal{F}+1},$
(6)

where $PF_{\text{EP,DP}}$ are the Purcell enhancement factors in the EP and DP
cases with respect to emission only into free space, i.e.,
$PF_{\text{EP,DP}}=\frac{2\gamma_{e}+P_{\text{EP,DP}}}{2\gamma_{e}}$. The
variables $\mathcal{R}$ and $\mathcal{F}$ in Eq. (6) are defined as
$\mathcal{R}\equiv\frac{P_{\text{EP}}}{P_{\text{DP}}}$ and
$\mathcal{F}\equiv\frac{2\gamma_{e}}{P_{\text{DP}}}=\frac{1}{PF_{\text{DP}}-1}$.
When $r=0$ (and $t=1$), the system operates at a DP, and hence, we have
$\eta=1$, as expected. In the limit of $r=1$ (i.e., perfect mirror), the
system is at the EP, and for $\Delta=0$, we have resonant emission with
$\chi=-2$. Under these conditions, we have $\mathcal{R}=1-\cos\Delta\phi$,
which attains its maximum value $\mathcal{R}_{\text{max}}=2$ at
$\Delta\phi=(2m+1)\pi$ for integer $m$. For $\gamma_{e}\ll P_{\text{DP}}$ or
equivalently $PF_{\text{DP}}\gg 1$ (a valid assumption for typical Purcell
enhancement settings), this translates into $\eta=2$. On the other hand, for
$\Delta\phi=(2m\pm 1/2)\pi$, we obtain $\eta=1$, i.e., equivalent to the DP
case at the same resonant frequency condition. Interestingly, for
$\Delta\phi=2m\pi$, we find that $\mathcal{R}=0$, and $\eta$ attains its
minimum value of $\eta_{\rm
min}=\frac{2\gamma_{e}}{2\gamma_{e}+P_{\text{DP}}}$, i.e., the QE radiates
with a rate $2\gamma_{e}$ as if it was located in free space.

Figure 2: Purcell factor enhancement $\eta$ as a function of the mirror field
reflectivity amplitude $|r|$ and $\Delta\phi$ (characterizing the change in
the position of the mirror or the QE) under the resonant condition $\Delta=0$.
The absorption of the mirror in the above plot is assumed to be zero. The
maximum enhancement $\eta=2$ occurs at $\Delta\phi=\pi$. For more realistic
mirrors with a finite absorption coefficient, the peak enhancement will be
less than its maximum possible value attainable in this ideal case.

These results show that the platform shown in Fig. 1 provides enough degrees
of freedom to tune the SE rate from that of free space emission to an enhanced
SE regime with an enhancement twice that of the DP case. This is also
evidenced in Fig. 2, which depicts $\eta$ as a function of both $|r|$ and
$\Delta\phi$ under the condition $PF_{\text{DP}}\gg 1$. An interesting feature
of Eq. (3) can be observed by writing
$B=\frac{Je^{i\phi_{E}}}{\left(\Delta+i\gamma_{W}\right)^{2}}\times\left(\Delta+i\gamma_{W}-2i\gamma_{W}|r|e^{-i\Delta\phi}\right)$.
When $\cos\Delta\phi=1/(2|r|)$ and $\Delta=\pm\gamma_{W}\sqrt{4|r|^{2}-1}$
(the sign depends on the value of $\Delta\phi$) with $|r|\geq 0.5$, we have
$B=0$, and $\eta=0.5$. In this case, the SE rate is suppressed compared with
that associated with the DP, but the emission becomes chiral with the photon
emitted only in the CW mode. This feature, which arises due to destructive
interference between the CCW mode and the back-reflected wave in the waveguide
inside the ring, has been recently observed experimentally using PT-symmetric
microwave and acoustic setups Chen  _et al._ (2020).

Before we proceed, it is worth mentioning that the above results can be
obtained by using the notion of PLDOS. In this case, the EP associated with
the photonic Hamiltonian gives rise to a double pole in the spectrum of the
Green’s operator. Consequently, the nonvanishing coupling between the QE and
the Jordan vector yields the square Lorentzian term of Eq. (4). Indeed, a
straightforward calculation of the Green’s operator associated with the
Hamiltonian matrix $H$ (after eliminating the third row and column that
account for the QE) yields expressions equivalent to Eq. (4). For a more
detailed mathematical discussion on how the Jordan vector modify the PLDOS and
SE, see Refs. Pick  _et al._ (2017a); Chen  _et al._ (2020); Khanbekyan and
Wiersig (2020).

Figure 3: Plots of (a) the value of Purcell factor (PF) enhancement $\eta$ as
a function of the normalized frequency detuning $\frac{\Delta}{\gamma_{W}}$
for $\Delta\phi=\pi$, and $|r|=0.976$; and (b) $\eta$ as a function of the
mirror (or QE) position as parametrized by $\Delta\phi$, under resonant
condition $\Delta=0$. In producing this plot, we used $|r|=0.976$, as obtained
for our realistic implementation of the mirror. Red solid line represents the
results obtained by using Eq. (6) together with the optical parameters of the
device, obtained as outlined in Appendix B. Black dots represent the results
obtained by full-wave simulations.

Next, we confirm the above predictions by performing two-dimensional full-wave
finite-difference time-domain (FDTD) simulations using realistic structure
dimensions and material systems. In our simulations, the microring resonator
has a refractive index of $n_{r}=3.47$ embedded in a background with
$n_{b}=1.44$. The outer radius of the microring is taken to be
$R=5\>\mu\text{m}$, and its width is $w=0.25\>\mu\text{m}$. The edge-to-edge
separation between the ring and the waveguide is $d=0.2\>\mu\text{m}$, and its
width is identical to that of the ring waveguide. The mirror at end of the
waveguide is made of a 100-nm-thick silver layer. These design parameters lead
to the following optical properties for the TE optical modes: effective
refractive index $n_{\text{eff}}=2.93$, $\omega_{o}=1216\>\text{THz}$ or
equivalently $\lambda_{o}=1549\>\text{nm}$, $\gamma_{W}=124\>\text{GHz}$,
corresponding to $Q=4900$, and $|r|=0.976$. Due to the optical absorption of
silver, we find $|t|=0.008$. The numerical evaluation of the above parameters
is presented in detail in Appendix B. The QE in our simulations is located as
shown in Fig. 1, i.e., at one-quarter of the perimeter as measured from the
ring-waveguide junction in the clockwise direction. Its dipole moment tensor
is assumed to have a component perpendicular to the ring plane and, thus,
couples only to the TE optical mode. In the FDTD simulations, the QE is
modeled by using an oscillating classical electric dipole. Note that these
calculations can provide information about the PF El-Ganainy and John (2013);
Vuckovic  _et al._ (1999) but not about the temporal ‘wavefunction’ of the
photon.

Figure 3(a) plots $\eta$ as a function of $\Delta/\gamma_{W}$ for the optimal
mirror position that maximizes the PF for the device described above. The
theory plot (red solid line) presents the results as obtained from Eq. (6) and
the optical parameters of the device which are extracted from FDTD simulations
(Appendix B for details). On the other hand, the curve denoted as simulation
(black dots) presents the same information as obtained directly by calculating
the optical power at the output ports of the waveguide from the FDTD
simulations (with and without the mirror). In both cases, we applied the
definition of $\eta$ under the condition $PF_{\text{DP}}\gg 1$. In our system,
this approximation can be justified by noting that, in our FDTD simulations,
the power emitted outside the microring resonator to free space is negligible
compared with that emitted to the waveguide ports. Evidently good agreement is
observed especially at the resonant frequency. Due to the finite reflectivity
of the mirror, the maximum enhancement here is $\eta=1.95$. Next, we plot the
values of $\eta$ as a function of $\Delta\phi$. In our simulations, we change
$\Delta\phi$ by scanning the mirror position around its optimal value.
Particularly, we varied the mirror position up to 440 nm with 20 nm
increments. Again, we observe a good agreement between theory and numerical
simulations.

To gain insight into these results, we investigate the electric field
distribution associated with various mirror positions as shown in Fig. 4,
focusing on the cases when $\Delta\phi=\pi$ (optimal mirror position) and
$\Delta\phi=0$ (trapping condition). In our simulations, these values
correspond to $L=5050$ and $L=5180\>\text{nm}$, respectively. From Fig. 4(a),
we observe that the field in the left port of the waveguide is uniform,
indicating an escaping traveling wave as expected. The field inside the ring
forms an imperfect standing wave pattern (the amplitude of the CCW component
three times larger than that of the CW component) having a peak at the QE
location, which explains the enhancement of the PF. On the other hand, the
field distribution in Fig. 4(b) demonstrates a perfect standing wave inside
the ring with a node located at the QE position. As a result, the photonic
mode effectively decouples from the QE, leading to near zero PF. In this
latter scenario, the microring and the right section of the waveguide with its
end mirror form a resonator that traps the excitation, with the lifetime
determined by the mirror reflectivity and the radiation rate from the ring to
free space.

Figure 4: Distribution of the electric field amplitudes in the structure under
semi steady state conditions for the two different extreme scenarios: (a)
$\Delta\phi=\pi$ ($L=5050$ nm) and (b) $\Delta\phi=0$ ($L=5180$ nm). In the
first case, the amplitudes of the clockwise (CW) and counterclockwise (CCW)
waves are different by a factor of three, larger [than the diabolic point (DP)
case] power output from $\text{P}_{1}$, which leads to maximum Purcell factor
(PF) enhancement. In the second case, the CW and CCW waves have identical
amplitudes, forming a perfect standing wave pattern inside the ring, with the
null located at the position of the quantum emitter (QE). This leads to a
decoupling between the QE and the optical modes. As a result, the PF is
suppressed ($\eta=0$) as evidenced by the absence of the optical power in the
output port $\text{P}_{1}$.

##  IV Conclusion

In summary, we have investigated the interaction between light and a QE in a
microring optical resonator exhibiting a chiral EP. Our analysis based on the
non-Hermitian Hamiltonian approach shows that, in the weak coupling regime,
the presence of an EP can enhance the PF by a factor of two. Furthermore,
implementing the chiral EP by side-coupling the resonator to a waveguide
terminated by a mirror can offer enough degrees of freedom to even
significantly suppress the SE process. These conclusions are confirmed by
using full-wave FDTD simulations. Our results open the door for designing more
efficient single photon sources. In future papers, we intend to explore this
possibility by studying the case of continuously driven QE near an EP.

##  Appendix A Evaluation of Purcell effect via eigenvalue analysis

In the main text, we used the method presented in Ref. Vuckovic  _et al._
(1999) to derive the $\eta$ factor. This was done to make a connection with
our computational study. In this section, we derive the same results by
directly solving Eq. (2). The general solution can be expressed as a sum of
eigenvectors, each scaled by an exponential time evolution factor given by the
corresponding eigenvalue of the matrix $H$. We note, however, that the Purcell
regime lies in the deep weak coupling domain where the decay rate of the QE is
only modified without significantly impacting the eigenvectors. Thus, one can
simplify the problem by focusing only on the correction introduced to the free
space decay rate $\gamma_{e}$ by the photonic environment. There are several
equivalent techniques to do so. For instance, one can solve the linear
equations using Laplace transformation and employing the so-called dominant
pole approximation El-Ganainy and John (2013); Teimourpour and El-Ganainy
(2015) or alternatively by using a perturbative expansion. Here, we use this
latter technique. Particularly, by seeking a series solution for the
eigenvalues of $H$ as a function of the perturbation parameter
$\epsilon=J/\alpha$ with $\alpha=\gamma_{W}-\gamma_{e}$, we find that the new
value for the otherwise bare (i.e., obtained with $J=0$ ) eigenvalue
$\omega_{e}-i\gamma_{e}$ is now given by

$\lambda_{e}=\omega_{e}-i\gamma_{e}-\delta+\mathcal{O}\\{\epsilon^{4}\\},$ (7)

where

$\displaystyle\delta=-\frac{J^{2}\left(2\Delta+\kappa
e^{-2i\phi_{E}}+2i\alpha\right)}{\left(\Delta+i\alpha\right){}^{2}}.$ (8)

Note that $\delta$ is in general a complex number. The real part of $\delta$
corresponds to the Lamb shift and is negligible in the Purcell regime. The
imaginary part, on the other hand, corresponds to the modification of the
decay rate of the excited state amplitude. Thus, the total amplitude SE rate
is now given by

$\displaystyle\gamma_{t}=\gamma_{e}+\text{Im}\\{\delta\\},$ (9)

where

$\begin{split}\text{Im}\\{\delta\\}=\frac{2J^{2}}{\left(\Delta^{2}+\alpha^{2}\right){}^{2}}\left[|r|\gamma_{W}\cos(\Delta\phi)(\Delta^{2}-\alpha^{2})\right.\\\
\left.-2\Delta|r|\gamma_{W}\alpha\sin(\Delta\phi)+\alpha(\Delta^{2}+\alpha^{2})\right].\end{split}$
(10)

Equation (9) indeed confirms the validity of our approximation of treating the
two decay processes (into free space and the waveguide) as two additive
independent effects. Moreover, one can quickly confirm the validity of the
expression in Eq. (10) in some simple scenarios. For instance, for the
resonant DP case ($r=0$ and $\Delta=0$), Im{$\delta$}$=2J^{2}/\alpha$, which
corresponds to the expression in the main text when $\gamma_{e}\ll\gamma_{W}$.
On the other hand, for the EP case with $|r|=1$, $\Delta=0$, and
$\Delta\phi=\pi$, we find Im{$\delta$}$=2J^{2}(\alpha+\gamma_{W})/\alpha^{2}$
also consistent with our results in the main text for
$\gamma_{e}\ll\gamma_{W}$. Finally, for $|r|=1$, $\Delta=0$, and
$\Delta\phi=0$, we obtain Im{$\delta$}$=2J^{2}(\alpha-\gamma_{W})/\alpha^{2}$
$\rightarrow$ 0 when $\gamma_{e}\ll\gamma_{W}$, as expected.

Figure 5: Comparison between the spontaneous emission (SE) decay rates as
obtained by an exact numerical solution of Eq. (2) (solid lines) and the
analytic expression of Eq. (9) (dashed lines) when the system operates at a
diabolic point (DP) ($|r|=0$) and when the system operates at an exceptional
point (EP) with maximum $\eta$ enhancement ($|r|=1$ and $\Delta\phi=\pi$). The
simulation parameters were taken to be $\omega_{o}=\omega_{e}$ ,
$\frac{J}{\gamma_{W}}=10^{-1}$, and $\frac{\gamma_{e}}{\gamma_{W}}=10^{-3}$.
For these parameters, using the approximation
$\gamma_{t}\approx\text{Im}\\{\delta\\}$ does not introduce any visible
changes to the plots, which justifies neglecting $\gamma_{e}$ in this regime.

In general, within the Purcell regime, where $\gamma_{e}\ll J\ll\gamma_{W}$
and $\gamma_{e}\ll J^{2}/\gamma_{W}$ Bozhevolnyi and Khurgin (2016), one can
use the approximation $\alpha\approx\gamma_{W}$ in the expression for
$\delta$. Moreover, in the SE enhanced decay regime, we can neglect
$\gamma_{e}$ compared with $\text{Im}\\{\delta\\}$. This is basically the
procedure we followed in solving Eq. (2) in the main text. In this case, the
above formula indeed gives the same result obtained in the main text using the
normalized emitted power, i.e., $2\text{Im}\\{\delta\\}=P_{\text{EP}}$, with
the equality holding also when $|r|=0$, i.e., in the DP case. However, as we
pointed out, in the SE decay suppression regime when
$\text{Im}\\{\delta\\}=0$, the excited QE will relax to its ground state by
emission into free space with rate $\gamma_{e}$.

Finally, we confirm the above results by plotting the function
$e^{-2\gamma_{t}t}$ and the full numerical solution of Eq. (2) on the same
graph for different cases, as shown in Fig. 5. Clearly, a good agreement is
observed for the relevant parameters listed in the figure caption.

##  Appendix B Numerical evaluation of the optical parameters of the structure

Figure 6: (a) An add-drop configuration is used to evaluate the decay rate
from the resonator to the waveguide. (b) Optical power transmission as a
function of frequency as obtained by using finite-difference time-domain
(FDTD). Based on these simulations, we estimate that $4\gamma_{W}=496$ GHz.

In the main text, we have compared the results obtained by using FDTD to those
estimated using the coupled mode theory (CMT) as represented by the
Hamiltonian $\hat{H}$. To use this later approach, however, one must evaluate
the essential optical parameters characterizing the system, such as the
coupling between the microring resonator and the waveguide (together with the
related optical loss from the ring to the waveguide port), the mirror complex
reflectivity, as well as its transmission (and hence its optical absorption).
In our paper, these parameters were computed numerically using FDTD. Here, we
outline the details of these calculations.

Figure 7: Absolute value of the electric field distribution when an optical
mode is launched from the left side of the waveguide onto the mirror (in the
absence of the resonator). Black dots are numerical data obtained from finite-
difference time-domain (FDTD), and the red curve is the fit according to the
expression of Eq. (11). The horizontal axis $x$ is the distance measured from
the position of waveguide-ring junction of the full geometry.

To extract the decay rate from the microring to the waveguide ($\gamma_{W}$),
we consider an add-drop ring resonator filter geometry, as shown in Fig. 6(a).
The amplitude of the power transmission coefficient as a function of the input
wave frequency is then obtained using FDTD [Fig. 6(b)]. By noting that the
full width at half maximum (FWHM) of the transmitted power is given by
$4\gamma_{W}$, we estimate that, in our case, $\gamma_{W}=124\>$GHz, which
together with $\omega_{o}=1216\>\text{THz}$ corresponds to a quality factor of
$Q=4900$ for the all-pass structure.

Next, we calculate the optical parameters of the mirror by removing the ring
resonator altogether and using FDTD simulations to evaluate the field
reflection and transmission coefficient from the mirror. In these simulations,
an incident optical mode having a free space wavelength of
$\lambda_{o}=1549\>$nm is launched into the waveguide from the left. The
absolute values of the reflection and transmission coefficient are easily
obtained by numerically measuring the reflected and transmitted optical
powers, respectively. Doing so gives $|r|=0.976$ and $|t|=0.008$. To evaluate
the phase of the reflection coefficient, we plot the absolute value of the
steady state electric field $|E|$ distributed at the incident wave side, as
shown in Fig. 7, and use curve fitting based on the analytical expression

$|E|=|E_{o}|\sqrt{1+|r|^{2}+2|r|\cos[2\beta(x-L)-\phi_{r}]},$ (11)

where here, $|E_{o}|$ is the incident wave amplitude, $|r|$ is the field
amplitude reflection coefficient, $L=3.95\ \mu\text{m}$ is the distance
between the mirror and the waveguide-ring junction, $\phi_{r}$ is the
reflection phase, and $\beta=\frac{2\pi
n_{\text{eff}}}{\lambda}=11.891\>\mu\text{m}^{-1}$. Based on these values, we
find that $\phi_{r}=3.56$ gives the best fit between Eq. (11) and the data in
Fig. 7. However, to obtain the best match between CMT and FDTD for the full
structure, we used $\phi_{r}=3.61$, which is an accepted discrepancy given the
approximate nature of CMT.

###### Acknowledgements.

R.E. acknowledges support from the Army Research Office (ARO; Grant No.
W911NF-17-1-0481), National Science Foundation (NSF; Grant No. ECCS 1807552),
and the Max Planck Institute for the Physics of Complex Systems. S.K.O.
acknowledges support from ARO (Grant No. W911NF-18-1-0043), NSF (Grant No.
ECCS 1807485), and Air Force Office of Scientific Research (Award No.
FA9550-18-1-0235).

Q.Z. and A.H. contributed equally to this work.

## References

  * Gross and Bloch (2017) C. Gross and I. Bloch, Science 357, 995 (2017).
  * Blais  _et al._ (2020) A. Blais, S. M. Girvin, and W. D. Oliver, Nat. Phys. 16, 247 (2020).
  * O’Brien (2007) J. L. O’Brien, Science 318, 1567 (2007).
  * Zeilinger (2017) A. Zeilinger, Physica Scripta 92, 072501 (2017).
  * Chang  _et al._ (2007) D. E. Chang, A. S. Sørensen, E. A. Demler, and M. D. Lukin, Nat. Phys. 3, 807 (2007).
  * Lvovsky  _et al._ (2009) A. I. Lvovsky, B. C. Sanders, and W. Tittel, Nat. Phys. 3, 706 (2009).
  * Buckley  _et al._ (2012) S. Buckley, K. Rivoire, and J. Vu$\mathrm{\check{c}}$kovi$\mathrm{\acute{c}}$, Rep. Prog. Phys. 75, 126503 (2012).
  * Purcell (1946) E. M. Purcell, Phys. Rev. 69, 681 (1946).
  * Yokoyama and Ujihara (2019) H. Yokoyama and K. Ujihara,  _Spontaneous Emission and Laser Oscillation in Microcavities_ (CRC Press, London, 2019).
  * Vahala (2003) K. J. Vahala, Nature (London) 424, 839 (2003).
  * Vahala (2004) K. Vahala,  _Optical Microcavities_ (World Scientific, Singapore, 2004).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, Z. H. Khajavikhan, S. Rotter, and D. N. Christodoulides, Nat. Phys. 14, 11 (2018).
  * Feng  _et al._ (2017) L. Feng, R. El-Ganainy, and L. Ge, Nat. Photon. 11, 752 (2017).
  * Özdemir  _et al._ (2019) Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nat. Mater. 18, 783 (2019).
  * Miri and Alù (2019) M.-A. Miri and A. Alù, Science 363, 363 (2019).
  * Lin  _et al._ (2016) Z. Lin, A. Pick, M. Lončar, and A. W. Rodriguez, Phys. Rev. Lett. 117, 107402 (2016).
  * Pick  _et al._ (2017a) A. Pick, B. Zhen, O. D. Miller, C. W. Hsu, F. Hernandez, A. W. Rodriguez, M. Soljačić, and S. G. Johnson, Opt. Express 25, 12325 (2017a).
  * Pick  _et al._ (2017b) A. Pick, Z. Lin, W. Jin, and A. W. Rodriguez, Phys. Rev. B 96, 224303 (2017b).
  * Hadad and Engheta (2020) Y. Hadad and N. Engheta, Proc. Natl. Acad. Sci. U.S.A. 117, 5576 (2020).
  * Khanbekyan and Wiersig (2020) M. Khanbekyan and J. Wiersig, Phys. Rev. Research 2, 023375 (2020).
  * Zhong  _et al._ (2019a) Q. Zhong, J. Ren, M. Khajavikhan, D. N. Christodoulides, Ş. K. Özdemir, and R. El-Ganainy, Phys. Rev. Lett. 122, 153902 (2019a).
  * Zhong  _et al._ (2019b) Q. Zhong, S. Nelson, Ş. K. Özdemir, and R. El-Ganainy, Opt. Lett. 44, 5242 (2019b).
  * Zhong  _et al._ (2020) Q. Zhong, S. K. Ozdemir, A. Eisfeld, A. Metelmann, and R. El-Ganainy, Phys. Rev. Applied 13, 014070 (2020).
  * Santori  _et al._ (2010) C. Santori, Y. Yamamoto, and D. Fattal,  _Single-Photon Devices and Applications_ (Wiley-VCH, Weinheim, Germany, 2010).
  * Vuckovic  _et al._ (1999) J. Vuckovic, O. Painter, Y. Xu, A. Yariv, and A. Scherer, IEEE J. Quantum Electron. 35, 1168 (1999).
  * Chen  _et al._ (2020) H.-Z. Chen, T. Liu, H.-Y. Luan, R.-J. Liu, X.-Y. Wang, X.-F. Zhu, Y.-B. Li, Z.-M. Gu, S.-J. Liang, H. Gao, L. Lu, L. Ge, S. Zhang, J. Zhu, and R.-M. Ma, Nat. Physics 16, 571 (2020).
  * El-Ganainy and John (2013) R. El-Ganainy and S. John, New J. Phys. 15, 083033 (2013).
  * Teimourpour and El-Ganainy (2015) M. H. Teimourpour and R. El-Ganainy, Journal of Optics 17, 125203 (2015).
  * Bozhevolnyi and Khurgin (2016) S. I. Bozhevolnyi and J. B. Khurgin, Optica 3, 1418 (2016).