# Analytical Phase Reduction for Weakly Nonlinear Oscillators

Iván León  Hiroya Nakao  Department of Systems and Control Engineering, Tokyo
Institute of Technology, Tokyo 152-8550, Japan

###### Abstract

Phase reduction is a dimensionality reduction scheme to describe the dynamics
of nonlinear oscillators with a single phase variable. While it is crucial in
synchronization analysis of coupled oscillators, analytical results are
limited to few systems. In this work, we analytically perform phase reduction
for a wide class of oscillators by extending the Poincaré-Lindstedt
perturbation theory. We exemplify the utility of our approach by analyzing an
ensemble of Van der Pol oscillators, where the derived phase model provides
analytical predictions of their collective synchronization dynamics.

##  1 Introduction

Populations of coupled oscillators are at the core of many natural and
technological systems, from the beating of the heart to the power transmission
in the electrical grid [1]. Although the relevance of this topic has motivated
extensive research in the last century, the large number of degrees of freedom
and intrinsic complexity of those systems have hampered complete
comprehension.

In the ’60s and ’70s, the seminal works of Winfree [2, 3] and Kuramoto [4, 5]
boosted our understanding of oscillatory dynamics by introducing phase
reduction. Phase reduction is a dimensionality reduction scheme that, assuming
weak coupling, derives a model capturing the dynamics of each oscillator in
terms of a single variable: the oscillator’s phase $\theta$. Derivation of the
Kuramoto model, which provided an analytical explanation of the
synchronization transition [4, 6, 1], clearly demonstrated the advantages of
phase reduction.

Since then, there has been an explosion of research employing phase models to
describe numerous systems, such as ensembles of neurons [7], mechanical
oscillators [8], or flashing of fireflies [9], to cite a few. In addition,
there have also been multiple advances in phase reduction, as generalizations
to include noise [10, 11], non-weak coupling [12, 13], higher-order
corrections [14, 15, 16, 17], or to deal with slow-fast systems [18]. Despite
these advances and the usefulness of phase reduction, analytical results are
limited only to a few types of oscillators for which the limit cycle and
isochrons are easily computable [3, 19].

In this work, we present an analytical perturbative approach to phase
reduction for oscillators whose unperturbed solution is known. The
paradigmatic example is weakly nonlinear oscillators that reduce to the
harmonic oscillator in the absence of nonlinearities. We achieve the phase
reduction by computing the frequency, limit cycle, and phase sensitivity
function (a.k.a. infinitesimal phase response curve) as a perturbative
expansion around the known solution. We illustrate the usefulness of the
approach by analyzing an ensemble of globally coupled Van der Pol oscillators
through the reduced phase model, obtaining analytical predictions of the
collective dynamics.

The method presented in this work broadens our knowledge about analytical
phase reduction and provides a useful tool to better understand the
oscillatory dynamics for a wide variety of systems. It derives simple,
analytical, and analyzable phase models that capture the weak coupling
dynamics, as exemplified by the results obtained for the Van der Pol
oscillator, a well-known model with many experimental applications [1, 20].

##  2 Phase reduction

Before presenting our results, it is convenient to make a brief introduction
to phase reduction [21]. Consider an $m$ dimensional oscillator whose state is
defined by $\boldsymbol{X}=(X_{1},\dots,X_{m})$. The dynamics of the
oscillator under the effect of some small perturbation
$\epsilon\boldsymbol{p}$ is described by the ordinary differential equation
(ODE):

$\dot{\boldsymbol{X}}=\boldsymbol{F}(\boldsymbol{X})+\epsilon\boldsymbol{p},$
(1)

where $\boldsymbol{F}(\boldsymbol{X})$ is the velocity field and overdot is
the derivative with respect to time $t$. In the absence of perturbation, Eq.
(1) displays a linearly stable limit cycle $\boldsymbol{X}_{c}$ with frequency
$\omega$. The perturbation $\boldsymbol{p}$ may depend on time $t$ and
oscillator state $\boldsymbol{X}$.

When the system is unperturbed, we define a phase variable $\theta$ on the
limit cycle that grows constantly from $0$ to $2\pi$ in one period of the
oscillation. This definition can be extended to the whole basin of attraction
of the limit cycle using the concept of isochron [3]. Isochrons foliate the
state space and assign to each point their asymptotic phase,
$\theta=h(\boldsymbol{X})$. Since the asymptotic phase of every point has to
be computed to obtain the isochrons, an analytical expression for the
isochrons is not generally accessible, except for a few types of oscillators
with polar symmetry.

The evolution of the phase in the presence of a perturbation is computed from
Eq. (1) and the definition of phase. Assuming the perturbation to be weak
compared to the second Floquet exponent $\lambda$, which characterizes the
stability of the limit cycle, i.e., $\epsilon\ll\lambda/\omega$, we obtain a
simple expression:

$\dot{\theta}=\omega+\epsilon\boldsymbol{Z}(\theta)\cdot\boldsymbol{p},$ (2)

where terms of order $\epsilon^{2}$ have been neglected [14, 22, 23] and
$\boldsymbol{Z}(\theta)=\nabla_{\boldsymbol{X}}\theta=\nabla_{\boldsymbol{X}}h(\boldsymbol{X}_{c}(\theta))$
is the gradient of the asymptotic phase $h({\boldsymbol{X}})$ evaluated on the
limit cycle $\boldsymbol{X_{c}}$ at phase $\theta$. The function
$\boldsymbol{Z}(\theta)$ is called the phase sensitivity function and
characterizes how much the phase is shifted under the effect of an
infinitesimal perturbation. In order for Eq. (2) to be closed, we
approximately evaluate $\boldsymbol{p}$ on the limit cycle at
$\boldsymbol{X_{c}}(\theta)$. We note the big dimensionality reduction from
the $m$-dimensional system Eq. (1) to the one-dimensional oscillator Eq. (2)
that captures the dynamics for small $\epsilon$.

From this derivation, we learned that to perform phase reduction, we only need
expressions for the frequency $\omega$, the limit cycle
$\boldsymbol{X_{c}}(\theta)$, and the phase sensitivity function
$\boldsymbol{Z}(\theta)$. Although $\boldsymbol{Z}(\theta)$ can be directly
computed from the asymptotic phase, $h(\boldsymbol{X})$, it is easier to
employ the Floquet theory. This theory states that the phase sensitivity
function is the periodic solution to the adjoint equation:

$\omega\frac{d\boldsymbol{Z}(\theta)}{d\theta}=-J^{\top}(\boldsymbol{X}_{c}(\theta))\boldsymbol{Z}(\theta)$
(3)

subject to the normalization condition
$\boldsymbol{Z}(\theta)\cdot\boldsymbol{F}({\boldsymbol{X}}_{c}(\theta))=\omega$
[24, 25, 26, 27], where $J^{\top}(\boldsymbol{X}_{c}(\theta))$ is the
transpose of the Jacobian matrix of $\boldsymbol{F}$ evaluated on the limit
cycle at ${\boldsymbol{X}}_{c}(\theta)$.

##  3 Perturbative method

In this work, we perform analytical phase reduction for a general class of
oscillators whose evolution depends on a parameter $\mu$:

$\dot{\boldsymbol{X}}=\boldsymbol{F}(\boldsymbol{X},\mu),$ (4)

which reduces to an oscillator whose solution $\boldsymbol{X}_{0}$ for $\mu=0$
is known. Weakly nonlinear oscillators belong to this class since they reduce
to a harmonic oscillator when nonlinearities are eliminated.

Provided that the solution for $\mu=0$ is known, we can develop a perturbation
theory around this solution. One has to be careful at this step since ordinary
perturbation theory produces secular terms, and thus, divergent solutions. We
solve this issue by using the Poincaré-Lindstedt method, although other
techniques such as multiple-timescale perturbation theory could also be
applied [28, 29, 30]. We chose the Poincaré-Lindstedt method for its
simplicity.

The first step in the Poincaré-Lindstedt method is to perform a change of
variable from time to phase, $\theta=\omega t$, where $\omega$ is yet unknown.
Then, we expand in powers of $\mu$ the frequency
$\omega=\omega_{0}+\mu\omega_{1}+\dots$ and the limit cycle solution
$\boldsymbol{X}=\boldsymbol{X}_{0}+\mu\boldsymbol{X}_{1}+\dots$, and replace
them in (4).

$\omega_{0}\frac{d\boldsymbol{X}_{0}}{d\theta}+\mu\left(\omega_{1}\frac{d\boldsymbol{X}_{0}}{d\theta}+\omega_{0}\frac{d\boldsymbol{X}_{1}}{d\theta}\right)+\dots=\boldsymbol{F}(\boldsymbol{X}_{0}+\mu\boldsymbol{X}_{1}+\dots,\mu)$
(5)

Expanding the right-hand side and collecting the terms at each order of $\mu$,
an ODE for each order is obtained.

We solve the ODEs recursively subject to an initial condition
${\boldsymbol{X}}(0)={\boldsymbol{A}}_{0}+\mu{\boldsymbol{A}}_{1}+\mu^{2}{\boldsymbol{A}}_{2}\cdots$
where ${\boldsymbol{A}}_{n}$ will be determined later. The solutions to the
ODEs, as in ordinary perturbative theory, contain secular terms. However, now
these terms depend on the undetermined constants $\boldsymbol{A}_{n}$ and
unknown frequencies $\omega_{n}$. We can uniquely choose those quantities so
that secular terms are removed. This method yields the frequency $\omega$ and
the limit cycle solution expressed in terms of the phase $\theta$ to any
desired order $\mu$.

As an example, we apply the method to the Van der Pol oscillator [31, 32]:

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle y,$
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle-x+\mu(1-x^{2})y.$ (6)

The analytic expressions obtained with the above procedure and their explicit
derivation can be found in A. In Fig. 1, we compare the approximate limit
cycles with the exact ones obtained by direct numerical simulation for
$\mu=0.5$ (a) and $\mu=0.7$ (b). First of all, we remark that the exact limit
cycle (black) is far from the perfect circle of the unperturbed harmonic
oscillator (blue), i.e., the $O(1)$ approximation. The higher-order
corrections contain new harmonics that reduce this discrepancy, as can be seen
in the $O(\mu)$ (orange) and $O(\mu^{2})$ (red) approximations. As expected,
the accuracy improves as higher orders are considered, while decreases as
$\mu$ increases. Strikingly, the order $\mu^{2}$ approximation is accurate
even for non-small values of $\mu$, as can be seen in panel (b).

Figure 1: Limit cycle of the Van der Pol oscillator for $\mu=0.5$ (a) and
$\mu=0.7$ (b). The exact numerical limit cycle and the $O(1)$, $O(\mu)$ and
$O(\mu^{2})$ approximations are depicted in black, blue, orange, and red,
respectively.

The crucial idea in this work is the extension of the Poincaré-Lindstedt
method to compute the phase sensitivity function $\boldsymbol{Z}$ by applying
a similar perturbative scheme to the adjoint equation (3). We expand the phase
sensitivity function in powers of $\mu$:
$\boldsymbol{Z}=\boldsymbol{Z}_{0}+\mu\boldsymbol{Z}_{1}+\dots$. Since
$\boldsymbol{Z}(\theta)$ is a $2\pi$-periodic function, all elements in the
expansion are also $2\pi$-periodic and do not contain secular terms. Plugging
the expansion of $\boldsymbol{Z}$ and the expressions for the frequency and
limit cycle into Eq. (3), we obtain

$\displaystyle\omega_{0}\frac{d\boldsymbol{Z}_{0}}{d\theta}+\mu\left[\omega_{1}\frac{d\boldsymbol{Z}_{0}}{d\theta}+\omega_{0}\frac{d\boldsymbol{Z}_{1}}{d\theta}\right]+\dots$
$\displaystyle=-J^{\top}(\boldsymbol{X}_{0}+\mu\boldsymbol{X}_{1}+\dots)(\boldsymbol{Z}_{0}+\mu\boldsymbol{Z}_{1}+\dots)$
(7)

and the normalization condition becomes

$\omega_{0}+\mu\omega_{1}+\dots=(\boldsymbol{Z}_{0}+\mu\boldsymbol{Z}_{1}+\dots)\cdot\boldsymbol{F}(\boldsymbol{X}_{0}+\mu\boldsymbol{X}_{1}+\dots,\mu).$
(8)

Expanding the right-hand side of both equations and collecting the terms at
each order of $\mu$, we obtain an ODE and normalization condition for each
order. The solutions to the system contain secular terms, which should be
removed by choosing appropriately the integration constants as in the standard
Poincaré-Lindstedt method. This procedure yields the phase sensitivity
function as a power expansion in $\mu$. See A for the derivation and
analytical results on the Van der Pol oscillator.

We now check the validity and accuracy of the method by computing the phase
sensitivity function of the Van der Pol oscillator. We depict the $x$ and $y$
components of $\boldsymbol{Z}$ vs the phase $\theta$ (see A for the
definition) in Figs. 2 (a) and (b), respectively, for $\mu=0.5$. The exact
results in black obtained by numerically solving the adjoint equation (3) are
compared with the $O(1)$, $O(\mu)$, and $O(\mu^{2})$ approximations in blue,
orange, and red, respectively. Figures 2 (c) and (d) depict the same
components of $\boldsymbol{Z}$ for $\mu=0.7$. The conclusions about accuracy
are analogous to those for the limit cycle. Remarkably, the $O(\mu^{2})$
approximation is extremely accurate even for $\mu=0.7$, indicating the
validity of the approximation up to those values.

Figure 2: Phase sensitivity function ${\boldsymbol{Z}}=(Z_{x},Z_{y})$ ($x$ and
$y$ components) of the Van der Pol oscillator for $\mu=0.5$ [(a) and (b)] and
$\mu=0.7$ [(c) and (d)]. The exact numerical $\boldsymbol{Z}$ and the $O(1)$,
$O(\mu)$, and $O(\mu^{2})$ approximations are depicted in black, blue, orange,
and red, respectively.

We note that, although the Poincaré-Lindstedt method is usually applied to
weakly nonlinear oscillators, the above procedure is not limited to them. For
example, it can be applied to oscillators that already possess a limit cycle
at $\mu=0$, as shown for a Stuart-Landau oscillator with a constant bias in B.

Finally, we remark that being able to obtain analytical expressions for the
frequency, limit cycle, and phase sensitivity function implies that we can
analytically perform phase reduction. Moreover, since those quantities are
independent of the perturbation ${\boldsymbol{p}}$, once they have been
computed, they can be used irrespective of the coupling types or topology.
Thus, it provides simple analytical phase models capturing the weak coupling
dynamics.

##  4 Ensemble of globally coupled Van der Pol oscillators

We now illustrate how the reduced phase model obtained by the presented
approach can help understanding the collective synchronization dynamics of
oscillatory systems. We choose, as a simple example, a population of $N$
globally coupled Van der Pol oscillators with identical properties. The $i$th
oscillator ($i=1,...,N$) is coupled to other oscillators through:

$\boldsymbol{p}_{i}(\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{N})=\frac{\epsilon}{N}\sum_{k=1}^{N}\bigg{(}(x_{k}-x_{i})\cos\alpha,(y_{k}-y_{i})\sin\alpha\bigg{)},$
(9)

where $\epsilon$ is the coupling strength and $\alpha$ measures the ratio of
$x$ and $y$ couplings. We remark that the parameter $\alpha$ is not the phase-
lag of the Kuramoto-Sakaguchi model, as will be observed when phase reduction
is performed. Denoting the phase of the $i$th oscillator as $\theta_{i}$ and
replacing (9) evaluated on the limit cycle, i.e., with
${\boldsymbol{X}}_{i}={\boldsymbol{X}}_{c}(\theta_{i})$, and the phase
sensitivity function into (2), we analytically obtain the reduced phase model
up to order $\mu^{2}$. We then perform averaging [5] to remove the fast-
oscillating terms and obtain the phase model:

$\dot{\theta}_{i}=\Omega+\epsilon\bigg{[}\eta_{1s}R_{1}\sin(\Psi_{1}-\theta_{i})+\eta_{1c}R_{1}\cos(\Psi_{1}-\theta_{i})\\\
+\eta_{3}R_{3}\sin(\Psi_{3}-3\theta_{i})\bigg{]},$ (10)

where we have defined the Kuramoto-Daido order parameters
$R_{k}e^{i\Psi_{k}}=\sum_{j}e^{ik\theta_{j}}$ for $k=1,3$ [33]. The new
constants are analytically represented by the original parameters as
$\Omega=-1+\mu^{2}/16-\epsilon\eta_{1c}$,
$\eta_{1c}=\mu(\sin\alpha+9\cos\alpha)/8$,
$\eta_{1s}=\frac{1}{2}\left[\sin\alpha\left(1+\mu^{2}/64\right)+\cos\alpha\left(1-7\mu^{2}/64\right)\right]$,
and $\eta_{3}=\mu^{2}\left(5\cos\alpha-3\sin\alpha\right)/128$. In contrast to
the ensemble of the original Van der Pol oscillators, the reduced phase model
(10) is straightforward to analyze following the standard techniques; see C
for details. In what follows, we consider the thermodynamic limit
$N\rightarrow\infty$.

Before analyzing the model, we point out that (10) displays a wide variety of
stable collective dynamics whose snapshots are shown in Fig. 3. Those dynamics
are triggered by the addition of a higher harmonic to the Kuramoto-Sakaguchi
model [34, 35, 36, 37]. Since $\eta_{3}\propto\mu^{2}$, for the specific
coupling (9), terms of order $\mu^{2}$ are needed to describe those new
dynamics.

Figure 3: Snapshots of FS (a), UIS (b), QPS (c), slow switching (d,e), and
three-cluster state (f) for $\epsilon=0.1$, $\mu=0.5$, and the displayed
$\alpha$. The blue dots are 60 randomly chosen oscillators from a population
of $N=300$ phase oscillators, the orange line is the evolution of the Kuramoto
order parameter $(R_{1}\cos\psi,R_{1}\sin\psi)$ and the red cross indicates
the position of the Kuramoto order parameter at the moment of the snapshot.
Panels (d) and (e) are two consecutive snapshots, where the blue and purple
color indicate which cluster the oscillator belong to.

Similarly to the Kuramoto-Sakaguchi model, Eq. (10) displays full synchrony
(FS), where all oscillators form one point cluster, Fig. 3 (a), and uniform
incoherence state (UIS), where oscillator phases are uniformly distributed,
Fig. 3 (b). It can be proven that FS is stable when $\eta_{1s}+3\eta_{3}>0$ ,
expressed in the original parameters as:

$\left(1-\frac{\mu^{2}}{8}\right)\sin\alpha+\left(1+\frac{\mu^{2}}{8}\right)\cos\alpha>0,$
(11)

while UIS is stable when $\eta_{1s}$ and $\eta_{3}$ are negative:

$\left(1+\frac{\mu^{2}}{64}\right)\sin\alpha+\left(1-\frac{7\mu^{2}}{64}\right)\cos\alpha<0,$
(12a) $-3\sin\alpha+5\cos\alpha<0.$ (12b)

The regions where FS and UIS are stable are depicted in Fig. 4 in blue and
yellow, respectively.

The phase model also allows us to analyze the stability of cluster states, in
which oscillators form $n$ point clusters. The analysis indicates that the
two-cluster states are unstable in the studied parameter region. Instead, we
observe slow switching (Sl-Sw), where the two-cluster states are saddles,
connected in a heteroclinic cycle [36, 37]. As the system approaches the
heteroclinic cycle, it switches between the saddle two-cluster states, Fig 3
(d). During the switches, one of the clusters is disintegrated and formed
again, Fig 3 (e). Slow switching is stable in the region where UIS and FS are
unstable and $\eta_{3}<0$. This region of the parameter space is depicted in
red in Fig. 4.

In addition, it is possible to show that a three-cluster state characterized
by $R_{1}=0$ and $R_{3}=1$, Fig. 3 (e), is stable when

$\left(1+\frac{19\mu^{2}}{64}\right)\sin\alpha+\left(1-\frac{37\mu^{2}}{64}\right)\cos\alpha<0,$
(13a) and $-3\sin\alpha+5\cos\alpha>0,$ (13b)

which is depicted with a green region in Fig. 4.

The last state present in the phase model is quasi-periodic partial synchrony
(QPS), Fig. 3 (c). In this state, the Kuramoto-Daido order parameters rotate
uniformly, while individual oscillators evolve in a quasi-periodic fashion
[38, 39, 34]. This state arises when UIS becomes unstable and loses its
stability via a subcritical Hopf bifurcation. Although the stability of QPS is
not determined analytically, it can be computed without numerical simulations.
We introduce an appropriate rotating frame of reference such that QPS is a
fixed point and use the Newton-Raphson algorithm to find the fixed point and
frequency of the rotating frame of reference. Then, a linear stability
analysis of the fixed point is used to determine the stability of QPS. The
region with stable QPS is depicted in Fig. 4 with gray hatching.

Figure 4: (a-c) Phase diagram predicted by the phase model (a) and magnified
regions near the bifurcations [(b) and (c)]. In the yellow, blue, green, and
red regions, UIS, FS, three cluster state, and slow switching are stable,
respectively. In the hatched blue-green region, there is bistability between
FS and three cluster state, while in the gray-red hatched region, both QPS and
slow switching are stable. (d-f) Phase diagram obtained by numerical
simulations of the Van der Pol oscillators (d) and magnified regions near the
bifurcations [(e) and (f)]. The same color scheme as in panel (a-c) is used.
The narrow green-white hatched region represents bistability between three-
cluster state and other multi-cluster states.

We sum up the previous discussions in Fig. 4 (a), where we depict the phase
diagram of Eq. (10) for $\epsilon=0.1$. Panels (b) and (c) depict
magnifications of the regions close to bifurcations in Fig. 4 (a). The regions
where FS, UIS, QPS, slow switching, and three-cluster state are stable are
depicted in blue, yellow, gray, red, and green, respectively, and the hatching
indicates bistability between two of those states. We stress that all the
stability boundaries but the one for QPS are analytically determined. See C
for the analytical derivations. We would like to remark that the nonlinearity
of the Van der Pol oscillator is responsible for the complex behavior observed
in the phase diagram, mainly due to the effect of the third- harmonic
component. See also e.g., [40] for the effect of high-harmonic components on
the dynamics of the Van der Pol oscillator.

To check the accuracy of the results obtained from the phase model, we perform
numerical simulations of an ensemble of $N=1000$ Van der Pol oscillators for
$\epsilon=0.1$. In Fig. 4 (d), we depict the phase diagram obtained from
numerical simulations for the ensemble of Van der Pol oscillators, with some
magnified regions in panels (e) and (f), and the same color code as Fig. 4
(a). It is remarkable that almost all the dynamical states predicted
analytically by the phase models are realized in the original system of Van
der Pol oscillators. The similarities between the phase diagrams for the Van
der Pol oscillators and the reduced phase model prove the latter is capturing
the qualitative dynamics with great accuracy.

There is only a small region in the parameter space [white-green hatched in
Fig. 4(d)] where there is bistability between the three-cluster state and a
multi-cluster state not present in the phase model. The fact that this region
shrinks as $\epsilon\rightarrow 0$ indicates that higher-order terms in
$\epsilon$ are needed to capture those dynamics. Additionally, we observe some
discrepancy for small $\mu$. This discrepancy is not surprising; it is due to
the limit cycle being less attractive as $\mu\rightarrow 0$, and thus the
assumption $\epsilon\ll\lambda/\omega$ for the phase reduction does not hold,
given the fixed $\epsilon=0.1$.

##  5 Conclusions

Through this work, we have presented a method to analytically obtain the
frequency, limit cycle, and phase sensitivity function as a perturbative
expansion in $\mu$ for a wide class of oscillators. These results enabled us
to analytically derive reduced phase models for perturbed or coupled
oscillators, facilitating an analytical understanding of their dynamics. We
have exemplified the methodology with the Van der Pol oscillator, showing its
usefulness to describe the dynamics and to obtain analytical phase diagrams.

We point to the importance of these results that broaden analytical phase
reduction to a wider variety of systems. This is especially remarkable for the
Van der Pol oscillator since it has multiple applications and is commonly used
as a test-bed model [1, 41, 20].

A possible extension of this work is to consider heterogeneous oscillators, in
which each unit has a different frequency. The presented methodology is easily
applicable to this case, though we have restricted ourselves to the identical
case for the sake of simplicity. We stress that the present approach can also
yield phase models with higher-order interactions [42] if the system contains
multi-body interactions. Finally, extending the present approach to compute
$\epsilon^{2}$ terms, where higher-order interactions appear naturally [14,
22, 23], is also possible, though it is more challenging and falls beyond the
scope of this work.

All in all, in this paper, we presented a systematic perturbative method to
derive analytical expressions for the phase model, which broadens the
knowledge about phase reduction. We believe the advances made in phase
reduction are an indispensable step to understand the ubiquitous oscillatory
dynamics.

## Acknowledgment

We acknowledge JSPS KAKENHI JP22K11919, JP22H00516, and JST CREST JP-MJCR1913
for financial support. We thank Diego Pazó for the fruitful discussions and
proposing the perturbation of Stuart-Landau oscillator with constant bias.

##  Appendix A Deriving frequency, limit cycle, and phase sensitivity function
of the Van der Pol oscillator

In this section, we explicitly derive the frequency, limit cycle, and phase
sensitivity function of the Van der Pol oscillator [31, 32], following the
method explained in the main text. We consider the evolution of the Van der
Pol oscillator as written in the form of (6) where $\mu$ is a small parameter.

###  A.1 Frequency and limit cycle

In order to obtain the frequency $\omega$ and the limit cycle
${\boldsymbol{X}}_{c}=(x_{c},y_{c})$, we employ the Poincaré-Lindstedt method
[28, 30]. This method only seeks the periodic solution of (6), and any
transient dynamic before the system reaches the limit cycle will be ignored.
The first step of the method is making the change of variable $\theta=\omega
t$, with unknown $\omega$. Then, we expand in powers of $\mu$ the frequency
$\omega=\omega_{0}+\mu\omega_{1}+\dots$, and the limit-cycle solution
$x_{c}=x_{0}+\mu x_{1}+\dots$ and $y_{c}=y_{0}+\mu y_{1}+\dots$. Plugging
those expansions into Eq. (6) and collecting each order of $\mu$, we obtain an
ordinary differential equation (ODE) for each order.

The $O(1)$ ODE yields the harmonic oscillator:

$\displaystyle\omega_{0}\frac{dx_{0}}{d\theta}$ $\displaystyle=$
$\displaystyle y_{0},$ $\displaystyle\omega_{0}\frac{dy_{0}}{d\theta}$
$\displaystyle=$ $\displaystyle-x_{0},$ (14)

where $\omega_{0}=-1$111One can alternatively consider $\omega_{0}=1$, which
is equivalent to defining the phase when $\mu=0$ as
$\theta=\arctan(-y_{0}/x_{0})$. In order to obtain the expressions for this
choice of the sign of $\omega_{0}$, it is only necessary to change
$\omega\rightarrow-\omega$,
$(x_{c}(\theta),y_{c}(\theta))\rightarrow(x_{c}(-\theta),y_{c}(-\theta))$, and
$\boldsymbol{Z}(\theta)\rightarrow-\boldsymbol{Z}(-\theta)$. is determined
from (6) with $\mu=0$ . We consider an initial condition on the limit cycle,
$y=0$ and $x=A_{0}+\mu A_{1}+\mu^{2}A_{2}\dots>0$, where the constants $A_{n}$
have to be determined later. The solution to (A.1) is then:

$x_{0}(\theta)=A_{0}\cos\theta;\quad y_{0}(\theta)=A_{0}\sin\theta.$ (15)

This solution is not yet determined due to the presence of $A_{0}$. We compute
the next order to determine $A_{0}$.

The order $\mu$ ODE is given by:

$\displaystyle\omega_{0}\frac{dx_{1}}{d\theta}$ $\displaystyle=$
$\displaystyle y_{1}-\omega_{1}\frac{dx_{0}}{d\theta},$
$\displaystyle\omega_{0}\frac{dy_{1}}{d\theta}$ $\displaystyle=$
$\displaystyle-x_{1}+(1-x_{0}^{2})y_{0}-\omega_{1}\frac{dy_{0}}{d\theta}.$
(16)

After plugging Eq. (15) into Eq. (A.1), we find the solution to this linear
nonhomogeneous system as the homogeneous solution (harmonic oscillations) plus
a particular solution. The latter is easy to find since all terms in the
r.h.s. are combinations of sines and cosines. This leads us to the solution to
Eq. (A.1):

$x_{1}(\theta)=\left(A_{1}+\theta\frac{A_{0}^{3}-4A_{0}}{8}\right)\cos\theta\\\
+\left(\frac{A_{0}}{2}-\frac{7}{32}A_{0}^{3}-A_{0}\theta\omega_{1}\right)\sin\theta+\frac{1}{32}A_{0}^{3}\sin
3\theta,$ (17a)
$y_{1}(\theta)=\left(A_{1}+\theta\frac{A_{0}^{3}-4A_{0}}{8}\right)\sin\theta\\\
+\frac{1}{32}\left(3A_{0}^{3}+32A_{0}\theta\omega_{1}\right)\cos\theta-\frac{1}{32}3A_{0}^{3}\cos
3\theta.$ (17b)

We observe that all the secular terms, for example
$\theta\frac{A_{0}^{3}-4A_{0}}{8}\cos\theta$ or
$\theta\omega_{1}A_{0}\sin\theta$, are removed by choosing $A_{0}=2$ and
$\omega_{1}=0$. Thus, the $O(1)$ solution, Eq. (15), is completely determined.
The $O(\mu)$ solution, Eq. (17), is now free of the secular terms, but it is
still not completely determined due to the presence of the constant $A_{1}$.
We can proceed to the following orders to determine $A_{1}$ and to obtain the
frequency and limit cycle as a power expansion in $\mu$ to any desired order.
Following this procedure, we compute the frequency and limit-cycle solution up
to order $\mu^{2}$,

$\omega=-1+\frac{\mu^{2}}{16}+O(\mu^{4}),$ (18)
$x_{c}(\theta)=2\cos\theta-\frac{\mu}{4}\left(3\sin\theta-\sin
3\theta\right)\\\ -\frac{\mu^{2}}{8}\left(\cos\theta-\frac{3}{2}\cos
3\theta+\frac{5}{12}\cos 5\theta\right)+O(\mu^{3}),$ (19)
$y_{c}(\theta)=2\sin\theta+\frac{3\mu}{4}\left(\cos\theta-\cos
3\theta\right)\\\ -\frac{\mu^{2}}{8}\left(2\sin\theta-\frac{9}{2}\sin
3\theta+\frac{25}{12}x\sin 5\theta\right)+O(\mu^{3}),$ (20)

which represents the limit cycle ${\boldsymbol{X}}_{c}$ of Eq. (6). Two
details that worth mentioning: first, the frequency $\omega$ only contains
even powers of $\mu$; and second, each new order in the limit-cycle solution
adds a new harmonic, and correct the coefficients of the previous orders.

###  A.2 Phase sensitivity function

Now that the limit cycle ${\boldsymbol{X}}_{c}$ and frequency $\omega$ have
been computed, we can use those expressions to obtain the phase sensitivity
function ${\boldsymbol{Z}}$ as a power expansion in $\mu$. We extend the
Poincaré-Lindstedt method to the adjoint equation for ${\boldsymbol{Z}}$:

$\omega\frac{d\boldsymbol{Z}(\theta)}{d\theta}=-J^{\top}(\boldsymbol{X}_{c}(\theta))\boldsymbol{Z}(\theta)$
(21)

with the normalization condition

$\boldsymbol{Z}(\theta)\cdot\boldsymbol{F}(\theta)=\omega.$ (22)

We expand the phase sensitivity function, $\boldsymbol{Z}=(Z_{x},Z_{y})$, in
powers of $\mu$:

$Z_{x}=Z_{x0}+\mu Z_{x1}+\mu^{2}Z_{x2}+\dots;\quad Z_{y}=Z_{y0}+\mu
Z_{y1}+\mu^{2}Z_{y2}+\dots$ (23)

It is important to have in mind that since $\boldsymbol{Z}$ is
$2\pi$-periodic, all terms in the expansion are also $2\pi$-periodic, and
thus, they do not contain secular terms. In other words, we seek for the
solution to the adjoint equation in the space of $2\pi$-periodic functions. We
plug this expansion and the expressions for the frequency and the limit cycle
into the adjoint equation, Eq. (21), and normalization condition, Eq. (22),
and collect the terms at each order of $\mu$. This provides an ODE and a
normalization condition for each order.

The $O(1)$ yields:

$\displaystyle\omega_{0}\frac{dZ_{x0}}{d\theta}$ $\displaystyle=$
$\displaystyle Z_{y0},$ $\displaystyle\omega_{0}\frac{dZ_{y0}}{d\theta}$
$\displaystyle=$ $\displaystyle-Z_{x0},$ (24)

with the normalization condition

$y_{0}Z_{x0}-x_{0}Z_{y0}=\omega_{0}.$ (25)

The solution to this system is

$Z_{x0}(\theta)=B_{0}\cos\theta-\frac{1}{2}\sin\theta,$ (26a)
$Z_{y0}(\theta)=B_{0}\sin\theta+\frac{1}{2}\cos\theta,$ (26b)

where $B_{0}$ is an integration constant, which is not determined yet. We
compute the next order to determine $B_{0}$.

The ODE and normalization condition for order $\mu$ are:

$\displaystyle\omega_{0}\frac{dZ_{x1}}{d\theta}$ $\displaystyle=$
$\displaystyle Z_{y0}+2x_{0}y_{0}Z_{y0}-\omega_{1}\frac{dZ_{x0}}{d\theta},$
$\displaystyle\omega_{0}\frac{dZ_{y1}}{d\theta}$ $\displaystyle=$
$\displaystyle-Z_{x0}-(1-x_{0}^{2})Z_{y0}-\omega_{1}\frac{dZ_{y0}}{d\theta},$
(27)

and

$y_{0}Z_{x1}-x_{0}Z_{y1}+y_{1}Z_{x0}-x_{1}Z_{y0}+(1-x_{0}^{2})y_{0}Z_{y0}=\omega_{1},$
(28)

whose solution is computed as a sum of the homogeneous solution and a
particular solution as

$Z_{x1}(\theta)=\left(-B_{0}\theta+B_{1}+\frac{9}{8}\right)\cos\theta\\\
-\frac{7}{8}B_{0}\sin\theta+\frac{5}{8}B_{0}\sin 3\theta+\frac{5}{16}\cos
3\theta,$ (29a)
$Z_{y1}(\theta)=\left(-B_{0}\theta+B_{1}+\frac{1}{8}\right)\sin\theta-\frac{1}{8}B_{0}\cos\theta\\\
+\frac{1}{8}B_{0}\cos 3\theta-\frac{1}{16}\sin 3\theta,$ (29b)

where $B_{1}$ is an undetermined integration constant. Since
${\boldsymbol{Z}}$ should be $2\pi$-periodic, we need to eliminate the secular
terms from these equations. It is easy to spot that, if $B_{0}=0$, all the
secular terms are removed. Proceeding forward, we can determine the phase
sensitivity function to any desired order. The result up to order $\mu^{2}$
is:

$Z_{x}(\theta)=-\frac{1}{2}\sin\theta+\frac{\mu}{16}\left(15\cos\theta+5\cos
3\theta\right)\\\ -\frac{\mu^{2}}{64}\left(21\sin\theta+\sin
3\theta-\frac{29}{6}\sin 5\theta\right)+O(\mu^{3}),$ (30)
$Z_{y}(\theta)=\frac{1}{2}\cos\theta-\frac{\mu}{16}\left(\sin\theta+\sin
3\theta\right)\\\ +\frac{\mu^{2}}{64}\left(3\cos\theta-5\cos
3\theta-\frac{1}{6}\cos 5\theta\right)+O(\mu^{3}).$ (31)

As in the expression for the limit cycle, higher-order corrections add new
harmonics and correct the coefficients in the previous orders.

##  Appendix B Stuart-Landau oscillator with a constant bias

In this section, to demonstrate the generality of our approach, we consider as
an example a Stuart-Landau oscillator with a constant bias,
$\dot{A}=(1+i\nu)A-|A|^{2}A+\mu$, written in polar coordinates $A=re^{i\phi}$
as

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle
r(1-r^{2})+\mu\cos\phi,$ (32) $\displaystyle\dot{\phi}$ $\displaystyle=$
$\displaystyle\nu-\frac{\mu}{r}\sin\phi,$ (33)

where $\nu$ is a free parameter. Following the procedure explained above, we
obtain the expressions for the frequency, limit cycle, and phase sensitivity
function up to order $\mu^{2}$:

$\omega=\nu-\mu^{2}\frac{2\mu^{2}}{\nu(4+\nu^{2})}+O(\mu^{4}),$ (34)
$r=1+\mu\frac{2\cos\theta+\nu\sin\theta}{\nu^{2}+4}-\frac{\mu^{2}}{\nu^{2}+4}\bigg{[}\frac{3}{4}+\cos\theta-\frac{2\sin\theta}{\nu}\\\
-\frac{\left(\nu^{4}+23\nu^{2}+4\right)\cos
2\theta}{4\left(\nu^{2}+1\right)\left(\nu^{2}+4\right)}+\frac{\left(-\nu^{4}+16\nu^{2}+8\right)\sin
2\theta}{2\nu\left(\nu^{2}+4\right)\left(\nu^{2}+1\right)}\bigg{]},$ (35a)
$\phi=\theta+\mu\frac{\cos\theta-1}{\nu}+\mu^{2}\bigg{[}\frac{1-\cos
2\theta}{2\nu\left(\nu^{2}+4\right)}+\frac{\sin\theta}{\nu^{2}}\\\
-\frac{\left(2+\nu^{2}\right)\sin
2\theta}{2\nu^{2}\left(\nu^{2}+4\right)}\bigg{]},$ (35b)
$Z_{r}=\mu\frac{\nu\cos\theta+2\sin\theta}{\nu^{2}+4}+\frac{\mu^{2}}{\nu^{2}+4}\bigg{[}\frac{8-4\nu^{2}}{\nu^{3}+4\nu}-\frac{2\cos\theta}{\nu}\\\
-\frac{2\nu\cos 2\theta}{1+\nu^{2}}+\sin\theta-\frac{(5+\nu^{2})\sin
2\theta}{2(\nu+\nu^{2})}\bigg{]},$ (36a)
$Z_{\phi}=1+\mu\frac{\sin\theta}{\nu}+\frac{\mu^{2}}{(\nu^{2}+4)}\bigg{[}\frac{8+6\nu^{2}}{\nu^{2}\left(\nu^{2}+4\right)}-\frac{2}{\nu^{2}}\cos\theta-\frac{\sin
2\theta}{\nu}\bigg{]}.$ (36b)

In figure 5 (a), we compare the approximate limit cycles of the Stuart-Landau
with a constant bias to the exact one obtained by direct numerical simulation
for $\mu=0.5$ and $\nu=2$. Although the exact limit cycle (black) is close to
be circular, it is far from the unperturbed Stuart-Landau limit cycle (blue),
i.e., the $O(1)$ approximation. The higher-order corrections contain new
harmonics that reduce this discrepancy, as can be seen in the $O(\mu)$
(orange) and $O(\mu^{2})$ (red) approximations. In panel (b) we depict the
angular component of the phase sensitivity function for the exact numerical
simulations in black and the $O(1)$, $O(\mu)$ and $O(\mu^{2})$ approximations
in blue, orange, and red, respectively. The conclusions about accuracy are
equivalent to the ones obtained for the Van der Pol oscillator.

Figure 5: Limit cycle of the Stuart-Landau oscillator with a constant bias for
$\mu=0.5$ and $\nu=2$ (a) and the angular component of its phase sensitivity
function. The exact numerical limit cycle and the $O(1)$, $O(\mu)$ and
$O(\mu^{2})$ approximations are depicted in black, blue, orange, and red,
respectively.

##  Appendix C Analyzing the phase model Eq. (10)

In this section, we derive the analytical stability conditions of the
dynamical regimes of the reduced phase model for the ensemble of globally
coupled Van der Pol oscillators, Eq. (10). Without loss of generality, we
bring the phase model to the form
$\dot{\theta}_{j}=\sum_{k}\Gamma(\theta_{k}-\theta_{j})/N$, with

$\Gamma(x)=\epsilon\bigg{[}\eta_{1s}\sin x+\eta_{1c}\cos
x+\eta_{3}\sin(3x)\bigg{]},$ (37)

by choosing a rotating reference frame with the frequency $\Omega$ of the
collective oscillation.

Although the stability analysis of the states can be performed for a finite
number of oscillators, we consider the thermodynamic limit
$N\rightarrow\infty$ for simplicity. In this limit, we can define a
probability density function of the oscillators $\rho(\theta,t)$, where
$\rho(\theta,t)d\theta$ measures the fraction of oscillators between $\theta$
and $\theta+d\theta$ at time $t$. This probability density evolves following
the continuity equation
$\partial_{t}\rho=-\partial_{\theta}(\rho\dot{\theta})$. We remark that the
Kuramoto-Daido order parameters are defined in this limit as
$R_{n}e^{i\psi_{n}}=\int_{0}^{2\pi}\rho e^{in\theta}d\theta$ for
$n=1,2,\ldots$.

###  C.1 Full synchronization

We begin our analysis with the study of full synchronization (FS). Full
synchronization is the state in which all oscillators form one point cluster,
and thus $R_{1}=R_{3}=1$. It is known that for coupling functions of the form
(37), the linear stability analysis of FS yields a zero exponent, associated
with a global phase shift, plus the exponents associated with variations of
individual oscillators. These perturbations grow with an exponent

$\lambda_{FS}=-\eta_{1s}-3\eta_{3}.$ (38)

Thus, when $\eta_{1s}+3\eta_{3}>0$, FS is stable as stated in the main text.

###  C.2 Uniform incoherent state

The uniform incoherent state (UIS) is characterized by vanishing order
parameters $R_{1}=R_{3}=0$, and thus the oscillator phases are uniformly
distributed. It is much easier to compute the stability of UIS if we decompose
the probability density of the oscillators, $\rho$, into its Fourier modes:

$\rho(\theta,t)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\rho_{k}(t)e^{ik\theta},$
(39)

where $\rho_{-k}=\rho_{k}^{*}$ and $\rho_{0}=1$. Notice that the coefficients
$\rho_{k}$ correspond to Kuramoto-Daido order parameters
$R_{k}e^{i\psi_{k}}=\int{{}_{0}^{2\pi}}\rho e^{ik\theta}{d\theta}=\rho_{-k}$.
Plugging this decomposition into the continuity equation, we obtain the
evolution of each mode $\rho_{k}$:

$\dot{\rho}_{k}=\frac{\epsilon
k}{2}\bigg{[}(\eta_{1s}-i\eta_{1c})\rho_{1}\rho_{k-1}-(\eta_{1s}+i\eta_{1c})\rho_{-1}\rho_{k+1}\\\
+\eta_{3}\big{(}\rho_{3}\rho_{k-3}-\rho_{-3}\rho_{k+3}\big{)}\bigg{]}.$ (40)

In this representation, the UIS is characterized by $\rho_{n\not=0}=0$. The
linear stability analysis of UIS is now straightforward to perform. It
indicates that UIS is stable when $\eta_{1s}$ and $\eta_{3}$ are negative, as
stated in the main text.

###  C.3 Quasiperiodic partial synchrony

In this section, we indicate how the analysis of QPS is performed. As in the
analysis of UIS, it is convenient to define a probability density $\rho$ of
the oscillators and work in the Fourier modes. Then, QPS is characterized by
the density of the oscillators rotating uniformly with frequency
$\tilde{\Omega}$, which means that the Fourier mode $\rho_{k}$ of $\rho$
evolves as

$\rho_{k}=e^{ik\tilde{\Omega}t}\tilde{\rho}_{k},$ (41)

where $\tilde{\rho}_{k}$ is a constant. In the rotating frame of reference
with frequency $\tilde{\Omega}$, QPS corresponds to a continuous family of
fixed points of equation (40), where all the fixed points are related through
a rotation. We can select one of the fixed points, for example, by imposing
$\operatorname{Re}(\tilde{\rho}_{1})=0$, and determine the values of
$\tilde{\rho}_{k}$ and frequency $\tilde{\Omega}$ using the Newton-Raphson
algorithm. Once those values are determined, we perform linear stability
analysis of the state. The spectrum will contain one zero exponent related to
the rotational symmetry of QPS and the other exponents determine the stability
of QPS. Because (40) is an infinite set of ODEs, it is is necessary to
truncate the system numerically, keeping only a finite number of modes.

###  C.4 Two-cluster states

The two-cluster state is characterized by a fraction $p$ of the oscillators
being in a point cluster $A$ with phase $\theta_{A}$, while the other $(1-p)$
fraction of the oscillators form cluster $B$ with phase $\theta_{B}$. Since
all oscillators are identical, to completely specify a two-cluster state, it
is enough to know the phase difference $\Delta=\theta_{A}-\theta_{B}$ and the
fraction $p$ of the oscillators in cluster $A$.

The evolution of the phase of each cluster is given by:

$\displaystyle\dot{\theta}_{A}$ $\displaystyle=$ $\displaystyle
p\Gamma(0)+(1-p)\Gamma(\theta_{B}-\theta_{A}),$ (42)
$\displaystyle\dot{\theta}_{B}$ $\displaystyle=$
$\displaystyle(1-p)\Gamma(0)+p\Gamma(\theta_{A}-\theta_{B}),$ (43)

and thus the phase difference $\Delta$ evolves according to:

$\dot{\Delta}=(2p-1)\Gamma(0)+(1-p)\Gamma(-\Delta)-p\Gamma(\Delta).$ (44)

The fixed points of this equation are the possible two-cluster states, and
thus the relation between $p$ and $\Delta$ is:

$p(\Delta)=\frac{\Gamma(0)-\Gamma(-\Delta)}{2\Gamma(0)-\Gamma(\Delta)-\Gamma(-\Delta)}.$
(45)

Any $\Delta$ such that $p(\Delta)\in(0,1)$ gives a possible two-cluster state,
nevertheless, it is not necessarily stable. To analyze the stability, we use
the fact that any perturbation given to the state can be decomposed into three
orthogonal modes [43, 37]. One mode corresponds to the locking of the two
clusters and the other two represent the disintegration of the individual
clusters. We denote their corresponding eigenvalues as $\lambda_{L}$,
$\lambda_{A}$, and $\lambda_{B}$.

We first analyze the locking of the two clusters by studying a perturbation
given to $\Delta$. Linearizing (44) around the fixed point $\Delta^{*}$, we
observe that the perturbation grows with the associated exponent

$\lambda_{L}=-(1-p)\Gamma^{\prime}(-\Delta^{*})-p\Gamma^{\prime}(\Delta^{*}).$
(46)

The disintegration of cluster $A$ is analyzed by computing the effect of
perturbing one oscillator from cluster $A$. This perturbation will grow with
the exponent

$\lambda_{A}=-p\Gamma^{\prime}(0)-(1-p)\Gamma^{\prime}(-\Delta^{*}).$ (47)

The exponent associated with the disintegration of cluster $B$ is obtained by
changing $p\rightarrow(1-p)$ and $\Delta^{*}\rightarrow-\Delta^{*}$ in the
previous equation. For the coupling function Eq. (37), at least one of the
eigenvalues is positive in the studied parameter space, hence two-cluster
states are always unstable in the studied parameter region.

###  C.5 Slow switching

A dynamical regime closely related to the two-cluster states is the slow
switching, a state in which two unstable two-cluster states are connected
forming an attractive heteroclinic cycle. It has been shown that for this
state to be stable, the following conditions have to be satisfied [37]:

  * 1.

A value of $p$ exists such that there are three unstable two-cluster states
with phase differences $0<\Delta_{1}<\Delta_{2}<\Delta_{3}<2\pi$. We denote
their eigenvalues by $\lambda_{L,A,B}^{n}$.

  * 2.

Full synchrony is unstable.

  * 3.

$\lambda_{L}^{2}>0$ while $\lambda_{L}^{1}<0$ and $\lambda_{L}^{3}<0$.

  * 4.

$\lambda_{A}^{1}>0$ and $\lambda_{B}^{1}<0$ while $\lambda_{A}^{3}<0$ and
$\lambda_{B}^{3}>0$.

The coupling function Eq. (37) verify those conditions in the parameter space
where UIS and FS are unstable and $\eta_{3}<0$, as stated in the main text.

###  C.6 Three-cluster states

The last state we analyze is the three-cluster state characterized by
$R_{1}=0$ and $R_{3}=1$. The values of the Kuramoto-Daido order parameters
imply that each cluster is formed by one third of the oscillators and the
phase difference between consecutive clusters is $2\pi/3$. As was done for the
two-cluster state, we can decompose the perturbation into different modes,
i.e., the locking of the clusters and their disintegration.

The locking of the clusters is analyzed by computing the evolution of the
phase differences between the clusters, $\Delta_{1}=\theta_{A}-\theta_{B}$ and
$\Delta_{2}=\theta_{B}-\theta_{C}$. The phase differences evolve as:

$\displaystyle 3\dot{\Delta}_{1}$ $\displaystyle=$
$\displaystyle\Gamma(-\Delta_{1})-\Gamma(\Delta_{1})+\Gamma(-\Delta_{1}-\Delta_{2})-\Gamma(-\Delta_{2}),$
(48) $\displaystyle 3\dot{\Delta}_{2}$ $\displaystyle=$
$\displaystyle\Gamma(\Delta_{1})-\Gamma(\Delta_{2})-\Gamma(\Delta_{1}+\Delta_{2})+\Gamma(-\Delta_{2}).$
(49)

The linear stability analysis of the above equations indicates that the
perturbation applied to the fixed point $\Delta_{1}=\Delta_{2}=2\pi/3$ grows
with exponents:

$\lambda_{1}^{\pm}=\frac{3}{2}(\eta_{1s}-6\eta_{3}\pm i\eta_{1c}),$ (50)

which are mutually complex conjugate.

To analyze the disintegration of one cluster, we follow the same procedure as
in full synchrony but with $R_{1}=0$. This yields the exponent

$\lambda_{2}=-\eta_{3}.$ (51)

Thus, the three-cluster state is stable when $\eta_{1s}-6\eta_{3}<0$ and
$\eta_{3}>0$, yielding the expressions Eq. (13) in the main text.

## References

  * [1] A. S. Pikovsky, M. G. Rosenblum, J. Kurths, Synchronization, a Universal Concept in Nonlinear Sciences, Cambridge University Press, Cambridge, 2001.
  * [2] A. T. Winfree, Biological rhythms and the behavior of populations of coupled oscillators., J. Theor. Biol. 16 (1967) 15–42.
  * [3] A. T. Winfree, The Geometry of Biological Time, Springer, New York, 1980.
  * [4] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in: H. Araki (Ed.), International Symposium on Mathematical Problems in Theoretical Physics, Vol. 39 of Lecture Notes in Physics, Springer, Berlin, 1975, pp. 420–422.
  * [5] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer-Verlag, Berlin, 1984.
  * [6] S. H. Strogatz, From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators, Physica D 143 (2000) 1–20.
  * [7] C. Bick, M. Goodfellow, C. R. Laing, E. A. Martens, Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review, The Journal of Mathematical Neuroscience 10 (1) (2020) 1–43.
  * [8] D. Mertens, R. Weaver, Synchronization and stimulated emission in an array of mechanical phase oscillators on a resonant support, Phys. Rev. E 83 (2011) 046221\.
  * [9] B. Ermentrout, N. Kopell, Parabolic bursting in an excitable system coupled with a slow oscillation, SIAM J. Appl. Math. 46 (1986) 233–253.
  * [10] D. S. Goldobin, J.-n. Teramae, H. Nakao, G. B. Ermentrout, Dynamics of limit-cycle oscillators subject to general noise, Phys. Rev. Lett. 105 (2010) 154101.
  * [11] J.-n. Teramae, H. Nakao, G. B. Ermentrout, Stochastic phase reduction for a general class of noisy limit cycle oscillators, Phys. Rev. Lett. 102 (2009) 194102.
  * [12] I. León, D. Pazó, Quasi phase reduction of all-to-all strongly coupled $\lambda-\omega$ oscillators near incoherent states, Phys. Rev. E 102 (2020) 042203.
  * [13] W. Kurebayashi, T. Yamamoto, S. Shirasaka, H. Nakao, Phase reduction of strongly coupled limit-cycle oscillators, Phys. Rev. Res. 4 (2022) 043176.
  * [14] I. León, D. Pazó, Phase reduction beyond the first order: The case of the mean-field complex Ginzburg-Landau equation, Phys. Rev. E 100 (2019) 012211.
  * [15] D. Wilson, B. Ermentrout, Phase models beyond weak coupling, Phys. Rev. Lett. 123 (2019) 164101.
  * [16] Y. Kuramoto, H. Nakao, On the concept of dynamical reduction: the case of coupled oscillators, Phil. Trans. R. Soc. A 377 (2160) (2019) 20190041.
  * [17] M. Rosenblum, A. Pikovsky, Numerical phase reduction beyond the first order approximation, Chaos 29 (1) (2019) 011105\.
  * [18] J. Zhu, Y. Kato, H. Nakao, Phase dynamics of noise-induced coherent oscillations in excitable systems, Phys. Rev. Res. 4 (2022) L022041.
  * [19] A. T. Winfree, Twenty-four hard problems about the mathematics of 24-hour rhythms, in: Lectures in Applied Mathematics, Nonlinear Oscillations in Biology, Vol. 17, 1979, pp. 93–126.
  * [20] L. Glass, M. C. Mackey, From Clocks to Chaos: The Rhythms of Life, Princeton University Press, 1988.
  * [21] H. Nakao, Phase reduction approach to synchronisation of nonlinear oscillators, Contemporary Physics 57 (2) (2016) 188–214.
  * [22] E. Gengel, E. Teichmann, M. Rosenblum, A. Pikovsky, High-order phase reduction for coupled oscillators, J. Phys.: Complexity 2 (1) (2021) 015005.
  * [23] I. León, D. Pazó, Enlarged Kuramoto model: Secondary instability and transition to collective chaos, Phys. Rev. E 105 (2022) L042201.
  * [24] F. C. Hoppensteadt, E. M. Izhikevich, Weakly connected neural networks., Spinger Verlag, N.Y., 1997.
  * [25] E. M. Izhikevich, Dynamical Systems in Neuroscience, The MIT Press, Cambridge, Massachusetts, 2007.
  * [26] B. Ermentrout, Type i membranes, phase resetting curves, and synchrony, Neural Comp. 8 (1996) 979–1001.
  * [27] E. Brown, J. Moehlis, P. Holmes, On the Phase Reduction and Response Dynamics of Neural Oscillator Populations, Neural Computation 16 (4) (2004) 673–715.
  * [28] A. H. Nayfeh, Introduction to Perturbation Techniques, Wiley, London, 1993.
  * [29] S. H. Strogatz, Nonlinear Dynamics and Chaos, Addison-Wesley, Reading, 1995.
  * [30] P. Amore, J. P. Boyd, F. M. Fernández, High order analysis of the limit cycle of the van der pol oscillator, Journal of Mathematical Physics 59 (1) (2018) 012702.
  * [31] B. van der Pol, Vii. forced oscillations in a circuit with non-linear resistance. (reception with reactive triode), The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 3 (13) (1927) 65–80.
  * [32] B. Van der Pol, J. Van der Mark, Frequency demultiplication, Nature 120 (1927) 363–364.
  * [33] H. Daido, Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all interactions: bifurcation of the order function, Physica D 91 (1996) 24–66.
  * [34] P. Clusella, A. Politi, M. Rosenblum, A minimal model of self-consistent partial synchrony, New Journal of Physics 18 (9) (2016) 093037\.
  * [35] P. Clusella, A. Politi, Irregular collective dynamics in a kuramoto–daido system, Journal of Physics: Complexity 2 (1) (2020) 014002.
  * [36] D. Hansel, G. Mato, C. Meunier, Clustering and slow switching in globally coupled phase oscillators, Phys. Rev. E 48 (1993) 3470–3477.
  * [37] H. Kori, Y. Kuramoto, Slow switching in globally coupled oscillators: robustness and occurrence through delayed coupling, Phys. Rev. E 63 (2001) 046214.
  * [38] M. Rosenblum, A. Pikovsky, Two types of quasiperiodic partial synchrony in oscillator ensembles, Phys. Rev. E 92 (2015) 012919.
  * [39] A. Politi, M. Rosenblum, Equivalence of phase-oscillator and integrate-and-fire models, Phys. Rev. E 91 (2015) 042916.
  * [40] D. S. Goldobin, Anharmonic resonances with recursive delay feedback, Physics Letters A 375 (39) (2011) 3410–3414.
  * [41] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer-Verlag, New York, 1983.
  * [42] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, G. Petri, Networks beyond pairwise interactions: Structure and dynamics, Phys. Rep. 874 (2020) 1 – 92.
  * [43] D. Golomb, D. Hansel, B. Shraiman, H. Sompolinsky, Clustering in globally coupled phase oscillators, Phys. Rev. A 45 (1992) 3516–3530.


# High-throughput optical absorption spectra for inorganic semiconductors

Ruo Xi Yang<EMAIL_ADDRESS>Materials Science Division, Lawrence Berkeley
National Lab  Matthew K. Horton  Materials Science Division, Lawrence Berkeley
National Lab  Jason Munro  Materials Science Division, Lawrence Berkeley
National Lab  Kristin A. Persson  Molecular Foundry, Lawrence Berkeley
National Lab  Department of Materials, University of California, Berkeley

###### Abstract

An optical absorption spectrum constitutes one of the most fundamental
material characteristics, with relevant applications ranging from material
identification to energy harvesting and optoelectronics. However, the database
of both experimental and computational spectra are currently lacking. In this
study, we designed a computational workflow for optical absorption spectrum
and integrated the simulated spectra into the Materials Project. Using
density-functional theory, we computed the frequency dependent dielectric
function and the corresponding absorption coefficient for more than 1000 solid
compounds of varying crystal structure and chemistry. The computed spectra
show excellent agreement, as quantified by a high value of the Pearson
correlation, with experimental results when applying the band gap correction
from the HSE functional. The demonstrated calculated accuracy in the spectra
suggests that the workflow can be applied in screening studies for materials
with specific optical properties.

††preprint: APS/123-QED

##  I Introduction

Optical absorption spectroscopy is an important material characterization for
identifying material species, atomic structures, reactions, and for
understanding the electronic structure and light-matter interaction.[1, 2]
Optical light wavelengths, spanning UV, visible light and near-IR (350 - 750
nm), are used to probe materials including crystalline semiconductors,
molecules, metals and amorphous materials. Due to the direct relation between
structure and physical properties, in situ UV-vis spectroscopy is commonly
used to monitor the crystallization, composition, and phase transition at
relevant time and length scales by tracking the evolving optical
properties.[3, 4, 5] Its non-destructive nature and inexpensive access makes
it a popular tool to probe dynamic excitation processes in materials.

The absorption process entails the interaction between incident
electromagnetic waves and the material electronic structure. Specifically in
semiconductors, the electrons are excited from the valence band to the
conduction band by absorbing the photon energy. Therefore, optical
spectroscopy is often used to determine the optical bandgap
($E_{\mathrm{g,opt}}$) by considering the onset of the absorption edge.[6]
Intrinsically, the optical bandgap is largely determined by the electronic
bandgap $E_{\mathrm{g,e}}$ between the band edges, considering only vertical
transitions due to lack of momentum transfer from the photons ($q\rightarrow
0$). However, in some materials where direct forbidden transitions are
present, $E_{\mathrm{g,opt}}$ can be larger than $E_{\mathrm{g,e}}$.
Additionally, when electron-hole Coulomb interaction is strong, resonant peaks
of exitons form below the band-to-band absorption edge, leading to a below-
bandgap $E_{\mathrm{g,opt}}$. [7, 8, 9, 10]. Moreover, charged defects with
energy levels close to the band edges can cause the absorption edge to broaden
and extend into lower energy range. [11] Extrinsically, the measured
$E_{\mathrm{g,opt}}$ is influenced by many factors, including experimental
setup, sample thickness etc. Various mathematical ways of extrapolating a
linear fitting of the spectral tailing often result in inconsistent
$E_{\mathrm{g,opt}}$ even in the same sample. [12, 13]

These factors interplay and complicate the determination of Eg,opt. To this
end, a database of internally consistent, computed optical absorption spectra
using a standardized approach will be of benefit in providing a common
reference point for experiments. By comparing with the absorption spectra of
crystalline, defect-free materials, experiments can infer whether other
influences, beyond band-to-band transitions, such as experimental setup and
interpretation, defect formation, excitons, electron-phonon coupling etc., are
present. From materials design perspective, these data provide an important
yet currently not widely available material descriptor, the optical band gap
$E_{\mathrm{g,opt}}$, for optoelectronic materials design.

In this study, we developed a high-throughput computational workflow for
generating optical absorption spectrum using density functional theory (DFT)
as implemented in the VASP package,[14, 15] with projector augmented wave core
potentials (reciprocal space projection). [16] We used the independent-
particle approximation (IPA) to compute the frequency dependent dielectrics
for more than 1000 solid materials chosen from the Materials Project.[17]
Benchmarks against the random-phase approximations (RPA) and with experiments,
using the Pearson correlation coefficients as metric, show excellent agreement
when incorporating a scissor shift of the difference between PBE band gap and
HSE band gap. We anticipate that the workflow and growing dataset will be
useful to the general materials community, and in particular for designing
optoelectronic materials.

##  II Methods

For each material, the structure was relaxed using the PBEsol functional,[18,
19] followed by a more rigorous relaxation with the r2SCAN functional,[20] as
implemented in VASP. [14, 15]. The structures were relaxed until the forces
were converged to 1 meV/Å, and the energy to 10-6 eV.

The static DFT calculation was performed using the tetrahedron smearing
method, with the final wavefunctions then used for the following optical
calculation.

The optical calculation was performed within the independent-particle
approximation (IPA), with two times the number of bands initially considered
in the static calculation to ensure the absorption spectrum covers a
sufficient energy range. The k-point density is chosen to be high (400/atom)
to ensure convergence against the density of states.

Here, the imaginary part of the frequency dependent microscopic dielectric
function is calculated :

$\begin{split}\epsilon_{2}=\frac{4\pi e^{2}}{\Omega}\lim_{q\rightarrow
0}\frac{1}{q^{2}}&\sum_{v,c,\mathrm{k}}2w_{\mathrm{k}}\delta(E_{ck}-E_{vk}-h\omega)\\\
&\times\innerproduct{u_{c,k+qe_{\alpha}}}{u_{vk}}\innerproduct{u_{vk}}{u_{c,k+qe_{\beta}}}\end{split}$
(1)

where $\Omega$ is the unit cell volume, $e$ is the electron charge, $q$ is the
photon momentum, $\omega$ is the photon frequency, and $c$, $v$ and $k$ denote
the conduction band, valence band, and the k-point, respectively. $u$ is the
electron wavefunction and the subscript denotes a k-point and band. Here,
$\ket{u_{nk+qe_{\alpha}}}$ can be obtained using first-order perturbation
theory.

The real part of the dielectric function is calculated via the Kramers-Kronig
relation:

$\epsilon_{1}(\omega)=1+\frac{2}{\pi}\int_{0}^{\infty}\frac{\epsilon_{2}(\omega^{\prime})\omega^{\prime}}{\omega^{\prime
2}-\omega^{2}}d\omega^{\prime}$

The optical absorption coefficient relies on the extinction coefficient
$\tilde{k}$ and the refractive index $\tilde{n}$, where

$\tilde{k}^{2}=\tilde{n}^{2}-\epsilon_{1}$

and $\tilde{n}$ follows:

$\tilde{n}=\frac{\epsilon_{2}}{2\tilde{k}}$

The frequency-dependent absorption coefficient depends on the extinction
coefficient and the incident photon frequency:

$\alpha=\frac{2\omega}{c}\tilde{k}$

Combining the above equations, we arrive at the absorption coefficient
$\alpha$ as a function of both the real and imaginary parts of the dielectric
function:

$\alpha=\frac{2\pi
E}{hc}\sqrt{2}\times\sqrt{\sqrt{\epsilon_{1}^{2}+\epsilon_{2}^{2}}-\epsilon_{1}}$

where $E$ is the incident photon energy and $c$ is the speed of light in
vacuum.

At the long-wavelength limit ($q\rightarrow 0$) the dielectric matrix
determines the optical properties accessible to optical probes. In IPA, the
macroscopic dielectric is approximated to be the same as the microscopic
dielectric by neglecting the local field effect, hence
$\epsilon_{\mathrm{mac}}\approx\mathrm{lim}_{\mathbf{q}\rightarrow
0}\epsilon_{0,0}(\mathbf{q},\omega)$. In order to include the local field
effect, the random-phase-approximation (RPA) calculation is needed to obtain
the response function $\epsilon_{\mathrm{mac}}=(\lim_{q\to
0}\epsilon_{0,0}(\mathbf{q},\omega))^{-1}$. The RPA calculation uses the
wavefunctions $u$ and the corresponding k-space derivatives $\partial
u/\partial k$ generated by the previous IPA calculation and computes the
response function $\chi$.

##  III High-throughput workflow

###  III.1 Screening

As a first step towards building an absorption coefficient database, we
consider solids suitable for photolvoltaic applications, specifically
semiconductors with an overlapping absorption spectrum with visible light.
Hence we screened compounds calculated by the Materials Project with three
criteria. First, an electronic band gap between 0.3 - 3 eV, corresponding to
4100 nm (infra-red) - 413 nm (violet) in photon wavelength; Second, to ensure
structural stability, the energy above hull was selected to be less than 0.02
eV/atom. This is comparable to the thermal energy $k_{B}T$ of 25 meV at room
temperature which can lead to entropic term that reduce the free energy.
Finally, we limit the number of sites in the unit cell to be less than 10 to
ensure a manageable memory footprint for the optical calculations. (Fig.1)
This set of filtering results in 1112 candidates, for which the IPA absorption
spectra were calculated and integrated into the Materials Project.

###  III.2 Workflow

We carried out a high-throughput computational workflow on the candidates
resulted from the screening step. The workflow consists of three parts: 1)
structure relaxation 2) static calculation and 3) optic calculation as
described in Fig.1. The structural relaxations are carried out using PBEsol
functional and an additional relaxation by r2SCAN functional to improve its
accuracy. The static calculation provides a more accurate total energyand the
ground state wavefunction. The following optical calculation at the IPA level
then calculates the dielectric functions, by which the absorption coefficient
is derived. To examine the importance of including the local field effects, we
benchmarked the RPA calculation against IPA for a subset of the materials
where experimental results are available. The validation is shown in the
following section.

Figure 1: The filtering process and the workflow steps used in this study.

##  IV Technical Validation

###  IV.1 IPA vs. RPA

The calculated spectra was benchmarked against experimental results, some of
which are shown in Fig 2. Even with the lower-level IPA calculations, the
experimental peak shapes and intensities are captured reasonably well. Small
spectral difference are observed between the IPA and RPA at higher energies,
where the local field effects are more prominent. [21, 22]

In order to numerically quantify spectral matching, we calculated the Pearson
correlation coefficient between the theoretical spectra and the experimental
one, using the equation:

$P=\frac{\Sigma(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\Sigma(x_{i}-\bar{x})^{2}\Sigma(y_{i}-\bar{y})^{2}}}$
(2)

where $x_{i}$ and $y_{i}$ are the calculated and experimental absorption
coefficients respectively, and $\bar{x}$ and $\bar{y}$ are the average value
of the $x_{i}$ and $y_{i}$, respectively. The Pearson coefficients range
between zero and one, and a value closer to 1.0 indicates highly correlated
spectra. The coefficients are plotted as color map in Fig.3. For the set of
materials with experimental data available, the RPA calculation shows marginal
improvement in $P$ in some materials. In most cases, the differences are
negligible, suggesting that the IPA is sufficient to reproduce the spectra.
Considering the significantly higher cost of the RPA calculation, only the IPA
spectrum was calculated for the rest of other compounds ($>$1000).

We note that for a few materials such as MgF2 and WS2, a sharp peak before the
main onset is measured but is absent in our calculations, indicating possible
formation of excitons which can only be captured by many-body theories, for
example as captured by the Bethe-Salpeter-Equation.

Figure 2: Example comparison between absorption spectrum by IPA, RPA and
experiments. Full data sets are included in SI. Figure 3: Pearson correlation
between simulated spectra and experiments

###  IV.2 HSE band gap correction

For some materials (e.g Si, GaAs, CdS, ZnSe etc.), the onsets of the
calculated IPA/RPA absorption are lower than experiments. This is an expected
outcome due to the underestimation of bandgaps by the PBE functional. To
correct the spurious self-interaction in GGA, we used the HSE functional to
obtain a more accurate bandgap,[23, 24] which includes short-range Hartree-
Fock exchange and is known to decrease the error and improve the agreement
between the calculated $E_{g}$ and corresponding experimental values. [25, 26]
The comparisons in $E_{\mathrm{g}}$ between HSE and PBE band gaps are shown in
Fig.4. As expected, most of the HSE band gaps are above or close ($<$ 0.2 eV)
to the PBE band gaps. A few outliers include: B ($\Delta
E_{g,\mathrm{HSE}}=-1.35$ eV), Al2O3 ($\Delta E_{g,\mathrm{HSE}}=-0.64$), H2O
($\Delta E_{g,\mathrm{HSE}}=-0.41$ eV), CaSO4 ($\Delta
E_{g,\mathrm{HSE}}=-4.09$ eV), Se ($\Delta E_{g,\mathrm{HSE}}=-1.18$ eV)

Using the HSE results, the spectra are shifted by $\Delta E_{g,\mathrm{HSE}}$,
the energy difference between the PBE band gap and HSE band gap. This results
in a rigid blue shift for most of the spectra, and in many cases better match
the experimental data (see Figure 2). A numerical description is reflected in
the Pearson correlation for the shifted spectra.Fig.3, where for materials
(e.g. CsI, ZnO, Te) with previously lower correlation, this scissor correction
improved $P$ significantly. This indicates that, to simulate a reasonably
accurate absorption spectra, carrying out a PBE-level frequency-dependent
dielectric calculation with correction from an HSE band gap is sufficient.

Figure 4: HSE and PBE band gaps for the materials used for benchmark. The
color bar represents $\Delta E_{g,\mathrm{,HSE}}$. The gray dashed line
corresponds to 1:1 ratio of $E_{g,\mathrm{,HSE}}$/$E_{g,\mathrm{,PBE}}$, and
the shaded region represent a $\pm 0.2$ eV window.

##  V Data Records

The calculated spectra are available through the Materials Project website
(www.materialsproject.org) and can be downloaded using the REST API. The
results are also available in the form of a JSON file that can be downloaded
directly from Materials Project.

#### Data representations:

The data for each of the calculated compounds are stored in a list and are
provided as a JSON file. For each compound, there are key values, such as
“energy” and “density”, that point to the appropriate property (Table 2).
Other keys include electronic band gap, imaginary and real part of the
dielectric functions, as well as the number of k-points.

json fields | units | physical property
---|---|---
energies | eV | Absorption energy corresponding to optic wavelength
energy_max | eV | The highest energy of the spectrum for this calculation
optical_absorption_co | cm-1 | Absorption coefficient
average_real_dielectric | none | Real part of the frequency-dependent dielectrics, averaged for each tensor
average_imag_dielectric | none | Imaginary part of the frequency-dependent dielectrics, averaged for each tensor
bandgap | eV | The electronic band gap
nkpoints | none | The number of kpoints used in the optic calculation

Table 1: The data fields of a downloadable absorption spectrum from MP.

##  VI Usage Note

We present a database of calculated optical absorption spectra for 1,116
compounds. The data should be of broad interest to materials applications
relating to electronic structure, particularly optoelectronic and other
electronic properties. For example, we expect this database to be used in the
understanding of optical absorption and in the search for new solar absorbers
with unique and tailored properties. The above use cases are facilitated by
the Materials Project website interface which allows users to search for
materials with associated metrics, such as stability, band gap and/or density.
With this new dataset and underlying data and software infrastructure, users
will now be able to request calculated optical absorption coefficients and
frequency-dependent dielectric functions. Furthermore, these capabilities open
up opportunities in data-driven endeavors, such as the application of machine
learning techniques to identify structural and chemical features that are key
to highly absorbing materials, hence extrapolating the data beyond its current
coverage and accelerating the discovery of novel optoelectronic materials.

##  VII Acknowledgement

This research used resources of the National Energy Research Scientific
Computing Center (NERSC), a U.S. Department of Energy Office of Science User
Facility located at Lawrence Berkeley National Laboratory, operated under
Contract No. DE-AC02-05CH11231. We thank Rachel Woods-Robinson for helpful
discussions. The work is supported by the US Department of Energy, Office of
Science, Office of Basic Energy Sciences, Materials Sciences and Engineering
Division under contract no. DE-AC02-05-CH11231 (Materials Project program
KC23MP).

## References

  * Tissue [2002] B. M. Tissue, Ultraviolet and Visible Absorption Spectroscopy, in  _Characterization of Materials_, edited by E. N. Kaufmann (John Wiley & Sons, Inc., Hoboken, NJ, USA, 2002) p. com059.
  * Zhang  _et al._ [2019a] L. Zhang, J. Ran, S.-Z. Qiao, and M. Jaroniec, Characterization of semiconductor photocatalysts, Chemical Society Reviews 48, 5184 (2019a).
  * Zhang  _et al._ [2019b] L. Zhang, T. Qian, X. Zhu, Z. Hu, M. Wang, L. Zhang, T. Jiang, J.-H. Tian, and C. Yan,  _In Situ_ optical spectroscopy characterization for optimal design of lithium–sulfur batteries, Chemical Society Reviews 48, 5432 (2019b).
  * Buckner and Perry [2019] H. B. Buckner and N. H. Perry, In Situ Optical Absorption Studies of Point Defect Kinetics and Thermodynamics in Oxide Thin Films, Advanced Materials Interfaces 6, 1900496 (2019).
  * Babbe and Sutter-Fella [2020] F. Babbe and C. M. Sutter-Fella, Optical Absorption-Based In Situ Characterization of Halide Perovskites, Advanced Energy Materials 10, 1903587 (2020).
  * Dolgonos  _et al._ [2016] A. Dolgonos, T. O. Mason, and K. R. Poeppelmeier, Direct optical band gap measurement in polycrystalline semiconductors: A critical look at the Tauc method, Journal of Solid State Chemistry 240, 43 (2016).
  * Koch  _et al._ [2006] S. W. Koch, M. Kira, G. Khitrova, and H. M. Gibbs, Semiconductor excitons in new light, Nature Materials 5, 523 (2006).
  * Rohlfing and Louie [1998] M. Rohlfing and S. G. Louie, Excitonic Effects and the Optical Absorption Spectrum of Hydrogenated Si Clusters, Physical Review Letters 80, 3320 (1998).
  * Elliott [1957] R. J. Elliott, Intensity of Optical Absorption by Excitons, Physical Review 108, 1384 (1957).
  * Ogawa and Takagahara [1991] T. Ogawa and T. Takagahara, Optical absorption and Sommerfeld factors of one-dimensional semiconductors: An exact treatment of excitonic effects, Physical Review B 44, 8138 (1991).
  * Redfield [1963] D. Redfield, Effect of Defect Fields on the Optical Absorption Edge, Physical Review 130, 916 (1963).
  * Raciti  _et al._ [2017] R. Raciti, R. Bahariqushchi, C. Summonte, A. Aydinli, A. Terrasi, and S. Mirabella, Optical bandgap of semiconductor nanostructures: Methods for experimental data analysis, Journal of Applied Physics 121, 234304 (2017).
  * Zanatta [2019] A. R. Zanatta, Revisiting the optical bandgap of semiconductors and the proposal of a unified methodology to its determination, Scientific Reports 9, 11225 (2019).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficient iterative schemes for _ab initio_ total-energy calculations using a plane-wave basis set, Physical Review B 54, 11169 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996b).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Physical Review Letters 100, 136406 (2008).
  * Perdew  _et al._ [2009] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Erratum: Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces [Phys. Rev. Lett. 100 , 136406 (2008)], Physical Review Letters 102, 039902 (2009).
  * Furness  _et al._ [2020] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, Accurate and Numerically Efficient r 2 SCAN Meta-Generalized Gradient Approximation, The Journal of Physical Chemistry Letters 11, 8208 (2020).
  * Wiser [1963] N. Wiser, Dielectric Constant with Local Field Effects Included, Physical Review 129, 62 (1963).
  * Louie  _et al._ [1975] S. G. Louie, J. R. Chelikowsky, and M. L. Cohen, Local-Field Effects in the Optical Spectrum of Silicon, Physical Review Letters 34, 155 (1975).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential, The Journal of Chemical Physics 118, 8207 (2003).
  * Heyd  _et al._ [2006] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)], The Journal of Chemical Physics 124, 219906 (2006).
  * Crowley  _et al._ [2016] J. M. Crowley, J. Tahir-Kheli, and W. A. Goddard, Resolution of the Band Gap Prediction Problem for Materials Design, The Journal of Physical Chemistry Letters 7, 1198 (2016).
  * Garza and Scuseria [2016] A. J. Garza and G. E. Scuseria, Predicting Band Gaps with Hybrid Density Functionals, The Journal of Physical Chemistry Letters 7, 4165 (2016).