# Nonunique fraction of Fock exchange for defects in two-dimensional materials

Wei Chen  Institute of Condensed Matter and Nanoscicence (IMCN), Université
catholique de Louvain, Louvain-la-Neuve 1348, Belgium  Sinéad M. Griffin
Molecular Foundry Division, Lawrence Berkeley National Laboratory, Berkeley,
California 94720, USA  Materials Sciences Division, Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA  Gian-Marco Rignanese  Institute
of Condensed Matter and Nanoscicence (IMCN), Université catholique de Louvain,
Louvain-la-Neuve 1348, Belgium  Geoffroy Hautier  Institute of Condensed
Matter and Nanoscicence (IMCN), Université catholique de Louvain, Louvain-la-
Neuve 1348, Belgium  Thayer School of Engineering, Dartmouth College, Hanover,
New Hampshire 03755, USA

###### Abstract

By investigating the vacancy and substitutional defects in monolayer WS2 with
hybrid functionals, we find that there is no unique amount of Fock exchange
that concurrently satisfies the generalized Koopmans’ condition and reproduces
the band gap and band-edge positions. Fixing the mixing parameter of Fock
exchange based on the band gap can lead to qualitatively incorrect defect
physics in two-dimensional materials. Instead, excellent agreement is achieved
with both experiment and many-body perturbation theory within the $GW$
approximation once the mixing parameters are tuned individually for the defect
species and the band edges. We show the departure from a unique optimized
mixing parameter is inherent to two-dimensional systems as the band edges
experience a reduced screening whilst the localized defects are subject to
bulklike screening.

Two-dimensional (2D) materials are rapidly gaining ground in a range of
technologies owing to both their exotic electronic and optical properties due
to quantum confinement [1] and the increasing ease of manipulation of their
properties on a layer-by-layer basis. Defects have been shown to be viable
active sites for enabling electrocatalysis (e.g., hydrogen evolution [2, 3,
4]) on otherwise inert basal planes of 2D materials. In quantum information
science, 2D materials hold exceptional promise as defect hosts because their
associated point defects typically present high spin states and deep defect
levels [5, 6, 7], and any near-surface defects are easier to manipulate and
characterize. The nature of such quantum defects is often elusive from
experiment, and reliable theoretical modeling is required to understand their
atomistic origins [8, 9]. First-principles defect calculations are now
typically carried out with hybrid density functionals, the accuracy of which
is closely related to the amount of Fock exchange admixed with semilocal
exchange in the framework of generalized Kohn-Sham density-functional theory
(DFT) [10]. The mixing parameter of hybrid functionals is commonly determined
on the basis of the band gap of host systems from experiment or higher levels
of theory. It has been well validated that such band-gap optimized mixing
parameters lead to accurate defect levels for localized point defects in bulk
materials [11].

The success of hybrid functionals can be attributed to the fulfillment of
generalized Koopmans’ condition, which for a localized defect, requires that
the ionization energy and electron affinity to be equal, hence recovering the
exact piecewise linearity of the total energy upon electron occupation [12,
13, 14, 15]. Indeed, one can often find a single mixing parameter that
describes the host band gap and the defect localization equally well in bulk
materials [16]. The uniqueness of optimal mixing parameter promotes the use of
a fixed, band-gap targeted mixing parameter in state-of-the-art defect
calculations [10]. The same approach has been routinely applied to understand
the defect properties in 2D materials [17, 18, 19, 20]. In this Letter, we
show that, in contrast to bulk 3D materials, there exists no unique amount of
Fock exchange that can reproduce the experimental band gap and satisfying
Koopmans’ conditions for the defect level concurrently. Admixing a single
fraction of Fock exchange thus could entail large errors in defect energy
levels for 2D materials. We show that imposing appropriate mixing parameters
individually for defects and host band edges is needed for an accurate
description of defect energy levels, and discuss the failure of the typical
hybrid-functional defect computation scheme.

To illustrate this issue, we investigate various point defects in monolayer
(ML) WS2, including the sulfur monovacancy ($V_{\text{S}}$), the single cobalt
atom substituting sulfur (Co${}_{\text{S}}$) or tungsten (Co${}_{\text{W}}$),
and the single carbon atom substituting sulfur (C${}_{\text{S}}$). Several of
these defects have recently become the subject of increasing theoretical and
experimental efforts following the demonstrated quantum light emitters in 2D
transition-metal dichalcogenides (TMDs) [21, 22, 23, 24, 25, 26], and so are
well-motivated for further consideration here [27, 28, 29, 19].

Figure 1: Single-particle defect levels as a function of mixing parameter
$\alpha$ for various defects in the ML WS2. The eigenvalues are evaluated at
the four representative $\alpha$ values shown by the solid dots and
interpolated in between (solid lines). Linear extrapolations (indicated by
dashed lines) are used when the localized defect state cannot be stabilized
within the band gap. The band edges of ML WS2 are indicated by the blue
shades. All energies are referred to the vacuum level. The orange marker
indicates the mixing parameter $\alpha_{\text{K}}$ at which the generalized
Koopmans’ condition
$\varepsilon^{\text{HO}}$($q$)$=\varepsilon^{\text{LU}}$($q+1$) is fulfilled,
whereas the grey marker indicates the optimal mixing parameter
($\alpha_{\text{G}}=0.22$) for reproducing the band gap of ML WS2. Spin-orbit
coupling (SOC) is not taken into account.

Our analysis of Koopmans’ condition is primarily carried out using a global
hybrid functional in which a single mixing parameter ($\alpha$) controls the
amount of Fock exchange admixed with the semilocal Perdew-Burke-Ernzerhof
(PBE) exchange [30]. At $\alpha=0.25$ this global hybrid functional is
essentially the PBE0 functional [31, 32]. Conventionally we refer to this
family of one-parameter hybrid functionals as PBE0($\alpha$). All hybrid-
functional calculations are performed with the projector-augmented-wave method
[33] as implemented in vasp [34, 35]. As a higher level of theory, many-body
perturbation theory within the $GW$ approximation [36, 37] is used to
establish an accurate reference for the hybrid functional. Specifically, on
top of a PBE starting point, one-shot $G_{0}W_{0}$ calculations are performed
for the ML WS2 and some selected defects. The truncated Coulomb interaction
[38] and the subsampling technique [39] are used to address the slow
convergence of quasiparticle (QP) energies specific to 2D systems [40]. The
$G_{0}W_{0}$ calculations are performed with berkeleygw [41] interfaced to the
DFT code from Quantum-espresso [42]. The details of computational parameters
for the hybrid-functional and $G_{0}W_{0}$ calculations are provided in the
Supplemental Material (SM) 111 See Supplemental Material for computational
details, finite-size corrections, and Koopmans’ condition with monolayer boron
nitride, which includes Refs. [67, 68, 69, 70, 71].

We first survey the electronic structure of the pristine ML WS2. Obtained from
exfoliating 2H-WS2, the ML WS2 is a direct-gap semiconductor with a $K$–$K$
transition. As a result of the trigonal prismatic crystal field splitting, the
valence band maximum (VBM) is of $d_{x^{2}-y^{2}}$ and $d_{xy}$ character
whereas the conduction band minimum (CBM) is of $d_{z^{2}}$ character. The
$G_{0}W_{0}$ band gap of the ML WS2 is 2.9 eV without SOC. To reproduce this
reference $G_{0}W_{0}$ band gap, we find that the mixing parameter (dubbed as
$\alpha_{\text{G}}$) needs to be adjusted to 0.22 for the hybrid functional.
Once the SOC is included, this PBE0($\alpha_{\text{G}}$) functional leads to a
band gap of 2.58 eV, in good agreement with the experimentally reported values
ranging from 2.4 to 2.7 eV [44, 27, 45]. We note that the experimental
determination of the band gap is subject to the dielectric screening effect
arising from the substrate [46, 47].

To identify the mixing parameter fulfilling the generalized Koopmans’
condition for a given localized defect $D$, we calculate the single-particle
eigenvalue $\varepsilon^{\text{HO}}(q)$ associated with the highest occupied
(HO) state at charge state $q$ as well as the eigenvalue
$\varepsilon^{\text{LU}}(q+1)$ of the lowest unoccupied (LU) state at charge
state $q+1$. The optimal mixing parameter $\alpha_{\text{K}}$ is determined
when the equality $\varepsilon^{\text{HO}}(q)=\varepsilon^{\text{LU}}(q+1)$
holds. To calculate this we consider candidate defects embedded in an
orthorhombic supercell of 90 atoms representing the ML WS2. In particular, we
refer the band edges and single-particle defect eigenvalues to the vacuum
level, which is a natural and physical choice for 2D systems [18]. While the
alignment is straightforward for neutral defects, the single-particle levels
of charged defects are subject to finite-size effect arising from the spurious
electrostatic interactions [48, 49, 11, 50]. The alignment is further
complicated by the fact that the electrostatic potential in the vacuum region
experiences an artificial bending due to the dipole moment introduced by the
neutralizing charge background. Here we apply the potential correction scheme
of Chagas da Silva et al. [51], thereby achieving well-defined single-particle
eigenvalues with respect to vacuum for charged defects.

Figure 2: (a) Defect levels of Co${}_{\text{S}}^{-1}$ and C${}_{\text{S}}^{0}$
obtained with single-particle eigenvalues and total-energy difference with
respect to the vacuum level. For each computational scheme, the PBE0($\alpha$)
defect levels are calculated with $\alpha_{\text{K}}$ and $\alpha_{\text{G}}$
and are benchmarked against the $G_{0}W_{0}$ reference. To exclude the effect
of structural relaxation, we use the PBE equilibrium structures throughout.
The PBE0 band-edge positions are determined by PBE0($\alpha_{\text{G}}$). (b)
Calculated single-particle defect levels with PBE0($\alpha_{\text{K}}$) and
PBE0($\alpha_{\text{G}}$) compared to STS measurements for $V_{\text{S}}^{0}$
and C${}_{\text{S}}^{-}$. The structures are relaxed with the corresponding
$\alpha$ values for the PBE0 calculations. The defect levels are referred to
the VBM and take into account the SOC. The band gap takes the value as
obtained with PBE0($\alpha_{G}$)+SOC throughout. The $G_{0}W_{0}$ results
(denoted by $G_{0}W_{0}^{\ast}$) are taken from Ref. 27.

Figure 1 depicts the evolution of band edges and single-particle defect
eigenvalues with respect to the mixing parameter $\alpha$. The defect
eigenvalues are obtained by linear interpolations or extrapolations based on
the eigenvalues calculated at four $\alpha$ values (0, 0.05, 0.10, and 0.22).
At a given $\alpha$, the defect structure upon which the Koopmans’ condition
is assessed corresponds to the equilibrium structure of the charged state. At
this stage, we do not take into account SOC for geometric relaxations and
hence the determination of $\alpha_{\text{K}}$ due to its computational
complexity. For all defects considered, the values of $\alpha_{\text{K}}$
fulfilling the Koopmans’ condition are found in a narrow range between 0.06
and 0.11, significantly lower than the band-gap enforced mixing parameter
$\alpha_{\text{G}}$ of 0.22. We note that $\alpha_{\text{K}}$ values are
largely unaffected by SOC as the Koopmans’ condition is still fulfilled within
0.1 eV when SOC is included (see Table S1 of the SM). This is our first key
result. Indeed, PBE($\alpha_{\text{G}}$) exhibits strong concavity as the
ionization energy $-\varepsilon^{\text{HO}}(q)$ is consistently larger than
the electron affinity $-\varepsilon^{\text{LU}}(q+1)$ for any given defect
[52]. This is clearly at odds with the existence of a unique mixing parameter
for bulk materials which is able to describe the band gap reasonably well
without compromising the Koopmans’ condition for localized defects [16, 53,
54, 55].

Among all defects considered, the effect of mixing parameter is particularly
notable for Co${}_{\text{S}}$. The defect levels associated with
Co${}_{\text{S}}$ are predominately characterized by the strongly localized
Co-$3d$ states. Specifically, the (0/$-$) and ($+$/0) transitions involve the
$d_{x^{2}-y^{2}}$ and the $d_{z^{2}}$ orbitals of the Co atom, respectively.
Taking the neutral Co${}_{\text{S}}^{0}$ for example, we notice that the
occupied $d_{z^{2}}$ and unoccupied $d_{x^{2}-y^{2}}$ defect levels move
sharply towards the host band edges as $\alpha$ increases and eventually merge
into the valence and conduction band when $\alpha=\alpha_{\text{G}}$ (cf. Fig.
1). Consequently, PBE0($\alpha_{\text{G}}$) would predict that no localized
defect levels exist within the band gap for Co${}_{\text{S}}$, qualitatively
at variance with the deep levels obtained with PBE0($\alpha_{\text{K}}$).

To settle the conflicting results of Co${}_{\text{S}}$ due to the choice of
$\alpha$, we carry out one-shot $G_{0}W_{0}$ calculations for the spin-
unpolarized Co${}_{\text{S}}^{-}$ defect (see the SM for computational
details). We find excellent agreement between the PBE0($\alpha_{\text{K}}$)
and $G_{0}W_{0}$ defect levels (mainly of $d_{x^{2}-y^{2}}$ and $d_{xy}$
characters of Co) [cf. Fig. 2(a)]. Both methods place the single-particle
defect level $-4.1$ eV below the vacuum level, or equivalently about 2.1 eV
above the VBM. By contrast, PBE0($\alpha_{\text{G}}$) overestimates the defect
ionization energy by nearly 1.3 eV, consistent with the strong concavity noted
previously in Fig. 1. The high accuracy of PBE0($\alpha_{\text{K}}$) is
further demonstrated for C${}_{\text{S}}^{0}$ [cf. Fig. 2(a)], while
PBE0($\alpha_{\text{G}}$) still underperforms, albeit to a lesser extent than
for Co${}_{\text{S}}^{-}$.

An arguably more common approach to the determination of defect levels is
through the total-energy difference between different defect charge states
[10]. When referred to the VBM of the pristine host, the (0/$-$) charge
transition level of defect $D$ can be obtained as
$E_{\text{tot}}(D^{-})-E_{\text{tot}}(D^{0})-E_{\text{VBM}}$. Analogous to the
single-particle level, the total energy of charged defects are ill-defined for
periodic systems. Here we apply the finite-size correction scheme of Komsa et
al. [50, 56, *Komsa2018] as implemented in slabcc [58] to the total energy of
charged defect. The computed (vertical) (0/$-$) charge transition levels are
shown in Fig. 2(a) for Co${}_{\text{S}}$ and C${}_{\text{S}}$. In general,
defect levels obtained with the total-energy difference scheme are less
sensitive to the amount of Fock exchange when aligned to a common reference
(e.g., average electrostatic potential or vacuum) [59, 60, 61, 62]. While this
is the case for the (0/$-$) level of C${}_{\text{S}}$, for Co${}_{\text{S}}$
the agreement with the $G_{0}W_{0}$ QP energy found with
PBE0($\alpha_{\text{K}}$) is substantially worse if $\alpha_{\text{G}}$ is
used instead.

Our results so far suggest that fixing the $\alpha$ value, which is typically
chosen based on the host band gap for studying defects in bulk materials, can
be problematic and even qualitatively change the defect physics in 2D
materials. The discernible discrepancy between $\alpha_{\text{G}}$ and
$\alpha_{\text{K}}$ suggests that defect levels need to be treated
individually from band edges, the alignment of which can be easily achieved
through the common vacuum level. To illustrate the validity the proposed
scheme, we refer to to the scanning tunneling spectroscopy (STS) measurements
of the neutral $V_{\text{S}}^{0}$ [27] and the negatively charged
C${}_{\text{S}}^{-}$ [29] in Fig. 2(b). To facilitate the comparison, we
include SOC in the hybrid-functional calculations and focus on the single-
particle eigenvalues with the equilibrium defect structure at the charge state
pertinent to experiment. The STS spectra revealed two prominent peaks at 1.8
eV and 2.0 eV above the VBM for $V_{\text{S}}^{0}$ [27]. The peaks correspond
to the unoccupied $d_{xy}$ and $d_{x^{2}-y^{2}}$ orbitals of the W atoms near
the vacancy, whose degeneracy is lifted by spin-orbit splitting [17]. The
positions of the two defect levels are well reproduced by the single-particle
defect energies of our PBE0($\alpha_{\text{K}}$) calculations in reference to
the band edges obtained with PBE0($\alpha_{\text{G}}$), and the accuracy is
comparable to that of the $G_{0}W_{0}$ calculation of Ref. [27]. The
C${}_{\text{S}}^{-}$ is also characterized by two defect levels deep within
the band gap [29]. The two levels are mainly of the C-$p_{z}$ character, with
one occupied at 1.1 eV and the other unoccupied at 2.0 eV above the VBM. The
importance of enforcing the Koopmans’ condition is again evidenced for
C${}_{\text{S}}^{-}$ as the splitting of the two defect levels would be too
large if the mixing parameter is fixed at $\alpha_{\text{G}}$.

Figure 3: (a) Single-particle defect levels at the Koopmans-fulfilled
$\alpha_{\text{K}}$ in the monolayer (blue) and the bulk (red) WS2. The band-
edge positions of the monolayer and the bulk WS2 are shown accordingly by the
same color scheme. (b) Defect charge density $\rho$($D^{q}$)$-\rho$($D^{0}$)
of Co${}^{-}_{\text{S}}$, C${}^{-}_{\text{S}}$, and $V^{-}_{\text{S}}$ as
obtained with PBE along the in-plane direction in the monolayer (blue) and the
bulk (red) WS2. The horizontal line refers to the delocalized screening charge
density $(1-1/\bar{\epsilon})q/\Omega$ given a localized defect at charge
state $q$ where $\bar{\epsilon}$ stands for the averaged dielectric constant
of the bulk WS2. The defects are situated near the center.

To rationalize the departure from the existence of a unique $\alpha$ value as
for bulk materials, we turn to the bulk WS2 and examine the defect levels
therein. The bulk 2H-WS2 comprises repeated layers of ML WS2, and has an
indirect ($\Gamma$–$K$) band gap of 1.3 eV [63, 64]. The markedly smaller band
gap leads to a strongly reduced mixing parameter
$\alpha_{\text{G}}^{\text{bulk}}$ of 0.07. Defects are introduced to an
orthorhombic supercell containing 72 atoms, which is based on the 36-atom
supercell of the ML WS2 comprising two alternating MLs stacked along the out-
of-plane direction. Figure 3(a) shows the single-particle defect levels at the
respective values of $\alpha_{\text{K}}$ for defects in both bulk and ML WS2.
Some defects such as Co${}_{\text{W}}$ no longer exhibit deep levels in the
bulk WS2. Interestingly, irrespective of the dimensionality of the host
material, the various $\alpha_{\text{K}}$ values fall closely in the vicinity
of $\alpha_{\text{G}}^{\text{bulk}}$. The consistency in the
$\alpha_{\text{K}}$ values implies that the localization of defect in the ML
WS2 is comparable to that of the bulk. This is clearly demonstrated in Fig.
3(b) where the defect charge density $\rho(D^{q})-\rho(D^{0})$ is plotted for
the ML and the bulk host. Particularly for Co${}_{\text{S}}^{-}$ and
C${}_{\text{S}}^{-}$, the charge density along the in-plane direction is
nearly intact when going from the bulk to the ML WS2. The charge density
analysis additionally provides an avenue to understand the screening effect.
In the presence of a strictly localized defect, the average charge density
approaches the delocalized screening charge density $(1-1/\epsilon)q/\Omega$
far away from the defect [49], where $\epsilon$ is the dielectric constant and
$\Omega$ the volume of the supercell. Using the (geometrically) averaged
dielectric constant ($\bar{\epsilon}=10.2$) of the bulk WS2, we find that the
delocalized screening charge is recovered for Co${}_{\text{S}}^{-}$ and
C${}_{\text{S}}^{-}$ [cf. Fig. 3(b)]. While this is expected for localized
defects in the bulk WS2, the agreement found in the ML WS2 is a compelling
evidence that the localized defects therein are still subject to the bulklike
screening and can be adequately described by
$\alpha_{\text{G}}^{\text{bulk}}$. The deviation from the ideal screening
charge is more visible for V${}_{\text{S}}$, in accord with its more
delocalized nature and hence the larger deviation of $\alpha_{\text{K}}$ from
$\alpha_{\text{G}}^{\text{bulk}}$ compared to the substitutional defects.

In contrast to $\alpha_{\text{K}}$, the large $\alpha_{\text{G}}$ for the ML
WS2 is inherent to the strong opening of band gaps for 2D materials. We note
in passing that the long-range part of Fock exchange is critical in opening
the band gap of 2D materials given the asymptotic $1/r$ decay in vacuum. In
fact, the short-ranged Heyd–Scuseria–Ernzerhof (HSE) functional [65] struggles
to open up the band gap of the ML WS2 unless an excessively large mixing
parameter is used (0.45 with the PBE structure or 0.55 if the structure is
relaxed self-consistently). However, such a large mixing parameter results in
a highly distorted band structure (e.g., valence bandwidth) and, more
importantly, shifts the whole band edges by $-0.3$ eV compared to $G_{0}W_{0}$
and PBE0($\alpha_{\text{G}}$) (see the SM). Nonetheless, for localized
defects, the HSE defect levels are reasonably aligned to the vacuum level as
the short-range bulklike screening dominates.

In addition to the TMDs, we find that the distinct two-set values
($\alpha_{\text{K}}$ and $\alpha_{\text{G}}$) also apply to the ML hexagonal
boron nitride (h-BN) involving only simple $sp$ elements (see the SM). While
systems characterized by localized $d$ electrons are less amenable to
treatment for fulfilling the Koopmans’ condition [66, 54], we emphasize that
the non-uniqueness of $\alpha$ is a general attribute for hybrid-functional
defect calculations in 2D materials.

Before closing, we note that our study focuses on defect energy levels. It
remains an open question how the nonuniqueness of $\alpha$ affects other
defect properties, such as the defect formation energies, the optical
transitions, and the phonons. Our findings motivate further investigations of
the possible implications for general defect physics in 2D materials.

In conclusion, we have shown that using a fixed amount of Fock exchange in
hybrid functionals could lead to qualitatively incorrect description of
defects in 2D materials, in contrast with the established practice common to
defect calculations in bulk materials. The absence of a unique mixing
parameter stems from the reduced screening for the delocalized band-edge
states whereas the screening is bulklike for localized defects. The distinct
screening behaviors hence require the band edges and the defect states to be
treated separately through their own optimal mixing parameters. We expect
these effects to be present in all 2D materials and they should be carefully
taken into account to ensure accurate hybrid-functional defect computations.

###### Acknowledgements.

This work was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences in Quantum Information Science under Award Number DE-
SC0022289. This research used resources of the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility supported
by the Office of Science of the U.S. Department of Energy under Contract No.
DE-AC02-05CH11231 using NERSC award BES-ERCAP0020966. Computational resources
were also provided by supercomputing facilities of UCLouvain (CISM) and
Consortium des Equipements de Calcul Intensif en Fédération Wallonie-Bruxelles
(CECI). The authors acknowledge Y. Ping, M. C. da Silva, and P. Deák for
fruitful discussions.

## References

  * Duan  _et al._ [2015] X. Duan, C. Wang, A. Pan, R. Yu, and X. Duan, Two-dimensional transition metal dichalcogenides as atomically thin semiconductors: opportunities and challenges, Chem. Soc. Rev. 44, 8859 (2015).
  * Ye  _et al._ [2016] G. Ye, Y. Gong, J. Lin, B. Li, Y. He, S. T. Pantelides, W. Zhou, R. Vajtai, and P. M. Ajayan, Defects engineered monolayer MoS2 for improved hydrogen evolution reaction, Nano Lett. 16, 1097 (2016).
  * Ouyang  _et al._ [2016] Y. Ouyang, C. Ling, Q. Chen, Z. Wang, L. Shi, and J. Wang, Activating inert basal planes of MoS2 for hydrogen evolution reaction through the formation of different intrinsic defects, Chem. Mater. 28, 4390 (2016).
  * Tsai  _et al._ [2017] C. Tsai, H. Li, S. Park, J. Park, H. S. Han, J. K. Nørskov, X. Zheng, and F. Abild-Pedersen, Electrochemical generation of sulfur vacancies in the basal plane of MoS2 for hydrogen evolution, Nat. Commun. 8, 15113 (2017).
  * Tran  _et al._ [2016] T. T. Tran, K. Bray, M. J. Ford, M. Toth, and I. Aharonovich, Quantum emission from hexagonal boron nitride monolayers, Nat. Nanotechnology 11, 37 (2016).
  * Caldwell  _et al._ [2019] J. D. Caldwell, I. Aharonovich, G. Cassabois, J. H. Edgar, B. Gil, and D. N. Basov, Photonics with hexagonal boron nitride, Nat. Rev. Mater. 4, 552 (2019).
  * Mendelson  _et al._ [2021] N. Mendelson, D. Chugh, J. R. Reimers, T. S. Cheng, A. Gottscholl, H. Long, C. J. Mellor, A. Zettl, V. Dyakonov, P. H. Beton, S. V. Novikov, C. Jagadish, H. H. Tan, M. J. Ford, M. Toth, C. Bradac, and I. Aharonovich, Identifying carbon as the source of visible single-photon emission from hexagonal boron nitride, Nat. Mater. 20, 321 (2021).
  * Turiansky  _et al._ [2019] M. E. Turiansky, A. Alkauskas, L. C. Bassett, and C. G. Van de Walle, Dangling bonds in hexagonal boron nitride as single-photon emitters, Phys. Rev. Lett. 123, 127401 (2019).
  * Ping and Smart [2021] Y. Ping and T. J. Smart, Computational design of quantum defects in two-dimensional materials, Nat. Comput. Sci. 1, 646 (2021).
  * Freysoldt  _et al._ [2014] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, First-principles calculations for point defects in solids, Rev. Mod. Phys. 86, 253 (2014).
  * Chen and Pasquarello [2013] W. Chen and A. Pasquarello, Correspondence of defect energy levels in hybrid density functional theory and many-body perturbation theory, Phys. Rev. B 88, 115104 (2013).
  * Perdew  _et al._ [1982] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Density-functional theory for fractional particle number: Derivative discontinuities of the energy, Phys. Rev. Lett. 49, 1691 (1982).
  * Lany and Zunger [2009] S. Lany and A. Zunger, Polaronic hole localization and multiple hole binding of acceptors in oxide wide-gap semiconductors, Phys. Rev. B 80, 085202 (2009).
  * Stein  _et al._ [2012] T. Stein, J. Autschbach, N. Govind, L. Kronik, and R. Baer, Curvature and frontier orbital energies in density functional theory, J. Phys. Chem. Lett. 3, 3740 (2012).
  * Perdew  _et al._ [2017] J. P. Perdew, W. Yang, K. Burke, Z. Yang, E. K. U. Gross, M. Scheffler, G. E. Scuseria, T. M. Henderson, I. Y. Zhang, A. Ruzsinszky, H. Peng, J. Sun, E. Trushin, and A. Görling, Understanding band gaps of solids in generalized Kohn-Sham theory, Proc. Natl. Acad. Sci. 114, 2801 (2017).
  * Miceli  _et al._ [2018] G. Miceli, W. Chen, I. Reshetnyak, and A. Pasquarello, Nonempirical hybrid functionals for band gaps and polaronic distortions in solids, Phys. Rev. B 97, 121112(R) (2018).
  * Li  _et al._ [2016] W.-F. Li, C. Fang, and M. A. van Huis, Strong spin-orbit splitting and magnetism of point defect states in monolayer WS2, Phys. Rev. B 94, 195425 (2016).
  * Smart  _et al._ [2018] T. J. Smart, F. Wu, M. Govoni, and Y. Ping, Fundamental principles for calculating charged defect ionization energies in ultrathin two-dimensional materials, Phys. Rev. Materials 2, 124002 (2018).
  * Li  _et al._ [2022a] S. Li, G. Thiering, P. Udvarhelyi, V. Ivády, and A. Gali, Carbon defect qubit in two-dimensional WS2, Nat. Commun. 13, 1210 (2022a).
  * Li  _et al._ [2022b] K. Li, T. J. Smart, and Y. Ping, Carbon trimer as a 2 eV single-photon emitter candidate in hexagonal boron nitride: A first-principles study, Phys. Rev. Materials 6, L042201 (2022b).
  * He  _et al._ [2015] Y.-M. He, G. Clark, J. R. Schaibley, Y. He, M.-C. Chen, Y.-J. Wei, X. Ding, Q. Zhang, W. Yao, X. Xu, C.-Y. Lu, and J.-W. Pan, Single quantum emitters in monolayer semiconductors, Nat. Nanotechnology 10, 497 (2015).
  * Srivastava  _et al._ [2015] A. Srivastava, M. Sidler, A. V. Allain, D. S. Lembke, A. Kis, and A. Imamoğlu, Optically active quantum dots in monolayer WSe2, Nat. Nanotechnology 10, 491 (2015).
  * Koperski  _et al._ [2015] M. Koperski, K. Nogajewski, A. Arora, V. Cherkez, P. Mallet, J.-Y. Veuillen, J. Marcus, P. Kossacki, and M. Potemski, Single photon emitters in exfoliated WSe2 structures, Nat. Nanotechnology 10, 503 (2015).
  * Chakraborty  _et al._ [2015] C. Chakraborty, L. Kinnischtzke, K. M. Goodfellow, R. Beams, and A. N. Vamivakas, Voltage-controlled quantum light from an atomically thin semiconductor, Nat. Nanotechnology 10, 507 (2015).
  * Palacios-Berraquero  _et al._ [2017] C. Palacios-Berraquero, D. M. Kara, A. R.-P. Montblanch, M. Barbone, P. Latawiec, D. Yoon, A. K. Ott, M. Loncar, A. C. Ferrari, and M. Atatüre, Large-scale quantum-emitter arrays in atomically thin semiconductors, Nat. Commun. 8, 15093 (2017).
  * Schuler  _et al._ [2020] B. Schuler, K. A. Cochrane, C. Kastl, E. S. Barnard, E. Wong, N. J. Borys, A. M. Schwartzberg, D. F. Ogletree, F. J. G. de Abajo, and A. Weber-Bargioni, Electrically driven photon emission from individual atomic defects in monolayer WS2, Sci. Adv. 6, eabb5988 (2020).
  * Schuler  _et al._ [2019] B. Schuler, D. Y. Qiu, S. Refaely-Abramson, C. Kastl, C. T. Chen, S. Barja, R. J. Koch, D. F. Ogletree, S. Aloni, A. M. Schwartzberg, J. B. Neaton, S. G. Louie, and A. Weber-Bargioni, Large spin-orbit splitting of deep in-gap defect states of engineered sulfur vacancies in monolayer WS2, Phys. Rev. Lett. 123, 076801 (2019).
  * Zhang  _et al._ [2019] F. Zhang, Y. Lu, D. S. Schulman, T. Zhang, K. Fujisawa, Z. Lin, Y. Lei, A. L. Elias, S. Das, S. B. Sinnott, and M. Terrones, Carbon doping of WS2 monolayers: Bandgap reduction and p-type doping transport, Sci. Adv. 5, eaav5003 (2019).
  * Cochrane  _et al._ [2021] K. A. Cochrane, J.-H. Lee, C. Kastl, J. B. Haber, T. Zhang, A. Kozhakhmetov, J. A. Robinson, M. Terrones, J. Repp, J. B. Neaton, A. Weber-Bargioni, and B. Schuler, Spin-dependent vibronic response of a carbon radical ion in two-dimensional WS2, Nat. Commun. 12, 7287 (2021).
  * Perdew  _et al._ [1996a] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996a).
  * Perdew  _et al._ [1996b] J. P. Perdew, M. Ernzerhof, and K. Burke, Rationale for mixing exact exchange with density functional approximations, J. Chem. Phys. 105, 9982 (1996b).
  * Adamo and Barone [1999] C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, J. Chem. Phys. 110, 6158 (1999).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996b).
  * Hedin [1965] L. Hedin, New method for calculating the one-particle green’s function with application to the electron-gas problem, Phys. Rev. 139, A796 (1965).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34, 5390 (1986).
  * Ismail-Beigi [2006] S. Ismail-Beigi, Truncation of periodic image interactions for confined systems, Phys. Rev. B 73, 233103 (2006).
  * da Jornada  _et al._ [2017] F. H. da Jornada, D. Y. Qiu, and S. G. Louie, Nonuniform sampling schemes of the Brillouin zone for many-electron perturbation-theory calculations in reduced dimensionality, Phys. Rev. B 95, 035109 (2017).
  * Qiu  _et al._ [2016] D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Screening and many-body effects in two-dimensional crystals: Monolayer MoS2, Phys. Rev. B 93, 235435 (2016).
  * Deslippe  _et al._ [2012] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G. Louie, BerkeleyGW: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures, Comput. Phys. Commun. 183, 1269 (2012).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with Quantum espresso, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Note [1] See Supplemental Material for computational details, finite-size corrections, and Koopmans’ condition with monolayer boron nitride, which includes Refs.[67, 68, 69, 70, 71].
  * Chernikov  _et al._ [2014] A. Chernikov, T. C. Berkelbach, H. M. Hill, A. Rigosi, Y. Li, O. B. Aslan, D. R. Reichman, M. S. Hybertsen, and T. F. Heinz, Exciton binding energy and nonhydrogenic rydberg series in monolayer WS2, Phys. Rev. Lett. 113, 076802 (2014).
  * Zhu  _et al._ [2015] B. Zhu, X. Chen, and X. Cui, Exciton binding energy of monolayer WS2, Sci. Rep. 5, 9218 (2015).
  * Naik and Jain [2018] M. H. Naik and M. Jain, Substrate screening effects on the quasiparticle band gap and defect charge transition levels in MoS2, Phys. Rev. Mater. 2, 084002 (2018).
  * Zibouche  _et al._ [2021] N. Zibouche, M. Schlipf, and F. Giustino, $GW$ band structure of monolayer MoS2 using the sternheimergw method and effect of dielectric environment, Phys. Rev. B 103, 125401 (2021).
  * Freysoldt  _et al._ [2009] C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Fully ab initio finite-size corrections for charged-defect supercell calculations, Phys. Rev. Lett. 102, 016402 (2009).
  * Komsa  _et al._ [2012] H.-P. Komsa, T. T. Rantala, and A. Pasquarello, Finite-size supercell correction schemes for charged defect calculations, Phys. Rev. B 86, 045112 (2012).
  * Komsa and Pasquarello [2013] H.-P. Komsa and A. Pasquarello, Finite-size supercell correction for charged defects at surfaces and interfaces, Phys. Rev. Lett. 110, 095505 (2013).
  * Chagas da Silva  _et al._ [2021] M. Chagas da Silva, M. Lorke, B. Aradi, M. Farzalipour Tabriz, T. Frauenheim, A. Rubio, D. Rocca, and P. Deák, Self-consistent potential correction for charged periodic systems, Phys. Rev. Lett. 126, 076401 (2021).
  * Bruneval [2009] F. Bruneval, $GW$ approximation of the many-body problem and changes in the particle number, Phys. Rev. Lett. 103, 176403 (2009).
  * Bischoff  _et al._ [2019] T. Bischoff, I. Reshetnyak, and A. Pasquarello, Adjustable potential probes for band-gap predictions of extended systems through nonempirical hybrid functionals, Phys. Rev. B 99, 201114(R) (2019).
  * Yang  _et al._ [2022] J. Yang, S. Falletta, and A. Pasquarello, One-shot approach for enforcing piecewise linearity on hybrid functionals: Application to band gap predictions, J. Phys. Chem. Lett. 13, 3066 (2022).
  * Deák  _et al._ [2019] P. Deák, M. Lorke, B. Aradi, and T. Frauenheim, Optimized hybrid functionals for defect calculations in semiconductors, J. Appl. Phys. 126, 130901 (2019).
  * Komsa  _et al._ [2014] H.-P. Komsa, N. Berseneva, A. V. Krasheninnikov, and R. M. Nieminen, Charged point defects in the flatland: Accurate formation energy calculations in two-dimensional materials, Phys. Rev. X 4, 031044 (2014).
  * Komsa  _et al._ [2018] H.-P. Komsa, N. Berseneva, A. V. Krasheninnikov, and R. M. Nieminen, Erratum: Charged point defects in the flatland: Accurate formation energy calculations in two-dimensional materials [Phys. Rev. X 4, 031044 (2014)], Phys. Rev. X 8, 039902(E) (2018).
  * Farzalipour Tabriz  _et al._ [2019] M. Farzalipour Tabriz, B. Aradi, T. Frauenheim, and P. Deák, slabcc: Total energy correction code for charged periodic slab models, Comput. Phys. Commun. 240, 101 (2019).
  * Alkauskas  _et al._ [2008] A. Alkauskas, P. Broqvist, and A. Pasquarello, Defect energy levels in density functional calculations: Alignment and band gap problem, Phys. Rev. Lett. 101, 046405 (2008).
  * Komsa  _et al._ [2010] H.-P. Komsa, P. Broqvist, and A. Pasquarello, Alignment of defect levels and band edges through hybrid functionals: Effect of screening in the exchange term, Phys. Rev. B 81, 205118 (2010).
  * Freysoldt  _et al._ [2016] C. Freysoldt, B. Lange, J. Neugebauer, Q. Yan, J. L. Lyons, A. Janotti, and C. G. Van de Walle, Electron and chemical reservoir corrections for point-defect formation energies, Phys. Rev. B 93, 165206 (2016).
  * Lyons and Van de Walle [2017] J. L. Lyons and C. G. Van de Walle, Computationally predicted energies and properties of defects in GaN, npj Comput. Mater. 3, 12 (2017).
  * Kam and Parkinson [1982] K. K. Kam and B. A. Parkinson, Detailed photocurrent spectroscopy of the semiconducting group VIB transition metal dichalcogenides, J. Phys. Chem. 86, 463 (1982).
  * Gusakova  _et al._ [2017] J. Gusakova, X. Wang, L. L. Shiau, A. Krivosheeva, V. Shaposhnikov, V. Borisenko, V. Gusakov, and B. K. Tay, Electronic properties of bulk and monolayer TMDs: Theoretical study within dft framework (GVJ-2e method), Phys. Status Solidi A 214, 1700218 (2017).
  * Heyd  _et al._ [2003] J. Heyd, G. E. Scuseria, and M. Ernzerhof, Hybrid functionals based on a screened coulomb potential, J. Chem. Phys. 118, 8207 (2003).
  * Ivády  _et al._ [2013] V. Ivády, I. A. Abrikosov, E. Janzén, and A. Gali, Role of screening in the density functional applied to transition-metal defects in semiconductors, Phys. Rev. B 87, 205201 (2013).
  * Schutte  _et al._ [1987] W. Schutte, J. de Boer, and F. Jellinek, Crystal structures of tungsten disulfide and diselenide, J. Solid State Chem. 70, 207 (1987).
  * Lin [2016] L. Lin, Adaptively compressed exchange operator, J. Chem. Theory Comput. 12, 2242 (2016).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * van Setten  _et al._ [2018] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comput. Phys. Commun. 226, 39 (2018).
  * Godby and Needs [1989] R. W. Godby and R. J. Needs, Metal-insulator transition in Kohn-Sham theory and quasiparticle theory, Phys. Rev. Lett. 62, 1169 (1989).


# Variational methods for finding periodic orbits in the incompressible
Navier-Stokes equations

J. P. Parker1  T. M. Schneider1  1 Emergent Complexity in Physical Systems
Laboratory (ECPS), École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne,
Switzerland

###### Abstract

Unstable periodic orbits are believed to underpin the dynamics of chaotic
dynamical systems, but by their nature are hard to find computationally. We
present a family of methods to converge such unstable periodic orbits for the
incompressible Navier-Stokes equations, based on variations of an integral
objective functional, and using traditional gradient-based optimisation
strategies. Different approaches for handling the incompressibility condition
are considered. The variational methods are applied to the specific case of
periodic, two-dimensional Kolmogorov flow and compared against existing Newton
iteration-based shooting methods. While computationally slow, our methods
converge from very inaccurate initial guesses.

##  1 Introduction

Exact coherent structures, and particularly unstable periodic orbits, are
believed to underpin the dynamics of chaotic dynamical systems, by acting as a
‘backbone’ of the chaos-supporting set in state space (Kawahara & Kida, 2001;
Kawahara et al., 2012; Cvitanović et al., 2016). The most successful method
for finding periodic orbits has been the shooting method based on Newton
iteration, which though fast and easy to implement on top of existing time-
marching codes, is notable for a relatively small radius of convergence,
requiring extremely accurate initial guesses. Extensions to Newton’s method,
employing a hookstep procedure, have been successful in enlarging this radius
of convergence (Viswanath, 2007; Schneider et al., 2010; Chandler & Kerswell,
2013; Dijkstra et al., 2014; van Veen et al., 2019; Tuckerman, 2020), but the
issue of accurate guesses remains. Alternative methods for finding periodic
orbits have been proposed, for example Lucas (2020) applied the idea of time-
delayed feedback to stabilise otherwise unstable periodic orbits, if a good
guess of the period is already known.

Lan & Cvitanović (2004) developed an alternative to the shooting method which
continuously deforms temporal loops in state space into periodic orbits. The
method was successfully applied to the Kuramoto-Sivshinsky equation (Lan &
Cvitanović, 2008), but it does not scale well to larger systems because of the
need for explicit Jacobian matrices. Boghosian et al. (2011) applied this
method to a Lattice-Boltzman formulation of Navier-Stokes, but noted the
prohibitively large computational requirements. Azimi et al. (2020) used a
matrix-free method inspired by Farazmand (2016) to find periodic orbits in the
Kuramoto-Sivashinsky equation. Though they successfully demonstrated the
method, challenges remain towards the ultimate goal of applying this method to
find unstable periodic orbits in fluid dynamics, as the studied one
dimensional equation differs from the two- and three-dimensional
incompressible Navier-Stokes equations in important aspects. In particular,
the divergence-free condition and the presence of the pressure field,
satisfying a Poisson equation with velocity-dependent source-term, render the
problem non-local. These non-local effects due to incompressibility are absent
in simpler, local, partial differential equations but require careful
consideration when formulating variational methods for computing periodic
orbits of Navier-Stokes.

As a stepping-stone towards full three dimensional turbulence, several
previous authors (Chandler & Kerswell, 2013; Lucas & Kerswell, 2015; Lucas,
2020) have studied the particular case of periodic, two-dimensional Kolmogorov
flow, at relatively low Reynolds number. This has several attractions: in two
dimensions the number of degrees of freedom in the system is greatly reduced,
yet this flow supports turbulence-like chaos; the lack of physical boundaries
simplifies computations, allowing the use of a Fourier (pseudo-)spectral
discretisation; and non-trivial invariant solutions, which have been studied
in detail for decades (Meshalkin & Sinai, 1961; Obukhov, 1983; Platt et al.,
1991; Fylladitakis et al., 2018) appear at relatively low Reynolds number,
allowing the use of low spatial resolutions. We also concentrate on the case
of 2D Kolmogorov flow, for the reasons listed above, as well as for comparison
with this previous work.

In this paper, we study methods which follow ideas of Azimi et al. (2020), but
which are formulated for 2D Kolmogorov flow and explicitly address the
incompressibility constraints of Navier-Stokes. In section 2 we present the
governing equations of the flow, formulate objective functionals that are
minimised at periodic orbits, and derive the necessary gradients of these.
Formulations are presented which tackle incompressibility in different ways:
either it is implicitly enforced by the governing equations; it is included in
the objective functional; or we explicitly project onto and restrict to the
incompressible subspace. In section 3, we discuss possible optimisation
strategies for these objective functionals. Section 4 compares three
variational methods for two test cases and then applies the most successful to
a larger number of possible periodic orbits, comparing against the traditional
shooting method. Concluding remarks are given in section 5.

##  2 Objective functionals for unstable periodic orbits

The well-known incompressible Navier-Stokes equations, in non-dimensional
form, describe the evolution of the 2D flow velocity $(u,v)$ with time $t$

$\displaystyle\partial_{t}u+u\partial_{x}u+v\partial_{y}u$
$\displaystyle=-\partial_{x}p+\frac{1}{Re}\Delta u+f_{x},$ (1)
$\displaystyle\partial_{t}v+u\partial_{x}v+v\partial_{y}v$
$\displaystyle=-\partial_{y}p+\frac{1}{Re}\Delta v+f_{y},$ (2)

subject to the incompressibility constraint

$\partial_{x}u+\partial_{y}v=0,$ (3)

where the pressure $p$ acts as a Lagrange multiplier to enforce the
incompressibility. Here $Re$ is the usual non-dimensional Reynolds number, and
$f_{x}$ and $f_{y}$ are the x- and y-components of a body force. Henceforth,
following previous authors (Chandler & Kerswell, 2013; Lucas, 2020), we will
take $f_{x}=\sin{4y}$, $f_{y}=0$, and consider a doubly periodic box of size
$2\pi\times 2\pi$. Throughout, we take $Re=40$, which is sufficiently large to
allow extensive chaotic behaviour, but permits the use of relatively low
resolutions to expedite computations. The system exhibits continuous symmetry
when shifting in the $x$-direction, and two discrete symmetries which have
been discussed extensively by previous authors. By differentiating (1) and
(2), we can derive a governing equation for the vorticity
$\omega=\partial_{x}v-\partial_{y}u$,

$\partial_{t}\omega+\partial_{y}\psi\partial_{x}\omega-\partial_{x}\psi\partial_{y}\omega=\frac{1}{Re}\Delta\omega+4\cos{4y},$
(4)

where the streamfunction $\psi$ is defined such that $u=\partial_{y}\psi$ and
$v=-\partial_{x}\psi$ and so $\omega=-\Delta\psi.$

These two formulations, called respectively primitive variables (PV) and
streamfunction-vorticity (SV), are equivalent in two dimensions. The latter
uses only two variables rather than three, since incompressibility is directly
built into the formulation, but has the major disadvantage that is does not
extend easily to three dimensions. It is further possible to eliminate
$\omega$ to give one equation solely in terms of $\psi$. As variational
methods based on the resulting fourth order equation were found to be very
inefficient, we do not discuss this formulation further.

Care must be taken when comparing the boundary conditions between the two
formulations. Assuming a doubly periodic streamfunction in the streamfunction-
vorticity formulation implies that $u$ and $v$ have no mean component over the
domain, i.e. no net flow, which is not a priori imposed by the primitive
variables formulation with periodic boundary conditions on $u$ and $v$.
However, it is straightforward to show that (1-3) conserve the net flow, which
is implicitly specified by the choice of initial conditions. Therefore, the
system modelled by the streamfunction-vorticity formulation represents an
invariant subspace of the primitive variables system. For initial conditions
with no net flow both formulations are equivalent, and we henceforth make this
choice.

Let us consider ‘loops’ in state space defined as for a field component

$u(x,y,s):[0,2\pi)\times[0,2\pi)\times[0,2\pi)\to\mathbb{R},$

periodic in all three dimensions, and assumed to be sufficiently smooth.
Similarly we can define loops of $v$, $p$, $\psi$ and $\omega$. Introducing a
period $T>0$ and a phase speed $c$, such loops are exact relative periodic
orbit solutions of the system if they satisfy the governing equations with
$\partial_{t}=\frac{2\pi}{T}\partial_{s}-c\partial_{x}$. In this perspective,
$s$, which parameterises the loop in phase space, is a scaled time variable.
From the governing equations of the two different formulations, we can define
objective functionals based on the PV formulation

$J_{PV}[u,v,p,c,T]=\frac{1}{2}\int\Bigg{\\{}\left(\frac{2\pi}{T}\partial_{s}u-c\partial_{x}u+u\partial_{x}u+v\partial_{y}u+\partial_{x}p-\frac{1}{Re}\Delta
u-\sin{4y}\right)^{2}\\\
+\left(\frac{2\pi}{T}\partial_{s}v-c\partial_{x}v+u\partial_{x}v+v\partial_{y}v+\partial_{y}p-\frac{1}{Re}\Delta
v\right)^{2}+\left(\partial_{x}u+\partial_{y}v\right)^{2}\Bigg{\\}}\mathrm{d}V,$
(5)

and for the SV formulation

$J_{SV}[\psi,\omega,c,T]=\frac{1}{2}\int\Bigg{\\{}\left(\frac{2\pi}{T}\partial_{s}\omega-c\partial_{x}\omega+\partial_{y}\psi\partial_{x}\omega-\partial_{x}\psi\partial_{y}\omega-\frac{1}{Re}\Delta\omega-4\cos{4y}\right)^{2}\\\
+\left(\omega+\Delta\psi\right)^{2}\Bigg{\\}}\mathrm{d}V,$ (6)

where the integrals are taken over $(x,y,s)\in[0,2\pi)^{3}$. Both functionals
are necessarily greater than or equal to zero, with equality if and only if
the governing equations are satisfied everywhere on the loop. Therefore,
finding zeros of these functionals gives us, in general, relative periodic
orbits of the system. Special cases include periodic orbits for which $c=0$,
travelling waves for which $\partial_{s}=0$, and equilibria, for which
$\partial_{s}=0$ and $c=0$. In these latter two cases, any value of $T$ will
give a zero of the functionals. While there is an exact correspondence between
zeros, i.e. global minima, of both objective functionals, local minima do not
generally translate from one formulation to the other.

To minimise these objective functionals, we derive gradients. For the scalar
variables $T$ and $c$, these are simple partial derivatives, for example

$\frac{\partial J_{SV}}{\partial
T}=\int-\frac{2\pi}{T^{2}}\partial_{s}\omega\left(\frac{2\pi}{T}\partial_{s}\omega-c\partial_{x}\omega+\partial_{y}\psi\partial_{x}\omega-\partial_{x}\psi\partial_{y}\omega-\frac{1}{Re}\Delta\omega-4\cos{4y}\right)\mathrm{d}V.$

For the loop variables $u$, $v$, $p$, $\psi$ and $\omega$, we must instead
derive variational derivatives which are defined at all
$(x,y,s)\in[0,2\pi)^{3}$, such as

$\frac{\delta
J_{SV}}{\delta\omega}=-\frac{2\pi}{T}\partial_{s}I+c\partial_{x}I-\partial_{x}(I\partial_{y}\psi)+\partial_{y}(I\partial_{x}\psi)-\frac{1}{Re}\Delta
I+(\omega+\Delta\psi),$

where for clarity we have defined

$I=\frac{2\pi}{T}\partial_{s}\omega-c\partial_{x}\omega+\partial_{y}\psi\partial_{x}\omega-\partial_{x}\psi\partial_{y}\omega-\frac{1}{Re}\Delta\omega-4\cos{4y}.$

Since we have no physical boundaries in this problem, such variational
derivatives are relatively straightforward to derive. All the gradients used
are given in the supplementary materials. It is important to note that the
gradients for the PV formulation have no mean component, assuming that $u$ and
$v$ do not, and so the (zero) net flow is conserved along gradients.

##  3 Optimisation methods for minimizing the objective functionals

Minimising the value of the objective functionals we have defined, which are
single valued functions of (at least in a discretised sense) a very large
number of unknown variables with known derivatives, is best performed through
gradient-based methods. This contrasts with methods by which we would directly
find zeros of the integrands through Newton iteration and related methods,
which require the same number of equations as unknowns. Descent methods find
only local minima of the function, but to find exact solutions of the
governing equations, we need to find zeros, i.e. global minima, of these non-
negative functionals. The methods we present give monotonically decreasing
and/or globally convergent results, but this does not mean we are guaranteed
to converge on exact solutions to the governing equations.

The simplest method is gradient descent with a fixed step size $\epsilon$ at
each iteration. Previous authors (Farazmand, 2016; Azimi et al., 2020) have
considered a dynamical system by the introduction of a fictitious time $\tau$
such that $\frac{du}{d\tau}=-\frac{\delta J}{\delta u}$. Simple gradient
descent is equivalent to discretising this system with a forward-Euler scheme,
and some improvements may be found by considering terms implicitly or using
higher order schemes, or adaptive timesteppers. However, we are not interested
in solving the fictitious system accurately, but merely in finding its stable
fixed points.

Faster convergence is found using the nonlinear conjugate-gradient method,
which requires only the gradients. We implemented this with the Fletcher &
Reeves (1964) choice of $\beta$, using the strong Wolfe (1969) conditions for
an inexact line search. Through trial and error, we found acceptable values
for the tunable parameters of the line search to be $c_{1}=10^{-5}$ and
$c_{2}=0.999$, though we make no claims that these are optimal. For $J_{SV}$
we used an initial step size $\alpha=10^{-7}$ and $\alpha=10^{-5}$ for
$J_{PV}$.

Since $x$, $y$ and $s$ are all periodic, the fields can be easily expressed as
Fourier series, using a pseudo-spectral approach – with 2/3 dealiasing – for
the nonlinear terms. A new code was developed for this work in C++, with
OpenMP parallelisation. Following Chandler & Kerswell (2013), we perform time
integrations using Heun’s method, with Crank-Nicolson on the viscous terms at
every substep. We concentrate on the specific case $Re=40$, which is
sufficiently high to allow fully chaotic behaviour, but low enough so that a
spatial resolution of $64\times 64$ gridpoints gives accurate results when
compared with previous work. A resolution of $64$ gridpoints was also used in
the temporal dimension, which was found to be more than sufficient.

##  4 Comparison of methods

In section 2 we gave objective functionals for two different formulations of
the system. The PV formulation does not assume the velocity field is
divergence-free, and instead includes the incompressibility condition within
the objective functional, whereas the alternative SV formulation automatically
enforces incompressibility. Since we know that any converged solution must be
divergence free, it may also be desirable to enforce this during the
optimisation within the PV formulation. However, care must be taken to ensure
that the resulting method still ensures the monotonic decrease of the
objective functional. For this, we recall the two-dimensional Leray projection
operator (Temam, 2001)

$\mathbf{P}:(u,v)\mapsto(u,v)-(\partial_{x},\partial_{y})\Delta^{-1}(\partial_{x}u+\partial_{y}v),$

where $\Delta^{-1}$ is the inverse Laplacian, here taken with periodic
boundary conditions. For any $(u,v)$, observe that $\mathbf{P}(u,v)$ is
divergence-free. Then for the usual inner product,

$\displaystyle\left<\mathbf{P}(u,v),(u,v)\right>$
$\displaystyle=\left<\mathbf{P}(u,v),\mathbf{P}(u,v)\right>+\left<\mathbf{P}(u,v),\nabla\Delta^{-1}(\partial_{x}y+\partial_{y}v)\right>$
$\displaystyle=\left\|\mathbf{P}(u,v)\right\|^{2}-\left<\nabla\cdot\mathbf{P}(u,v),\Delta^{-1}(\partial_{x}u+\partial_{y}v)\right>$
$\displaystyle=\left\|\mathbf{P}(u,v)\right\|^{2}\geq 0,$

where we have used the fact that $\nabla$ is adjoint to $-\nabla\cdot$ with
respect to the inner product. Therefore, using an update direction of
$-\mathbf{P}\left(\frac{\delta J_{PV}}{\delta u},\frac{\delta J_{PV}}{\delta
v}\right)$ instead of $-\left(\frac{\delta J_{PV}}{\delta u},\frac{\delta
J_{PV}}{\delta v}\right)$ will still converge to a minimum of the objective
functional, since the projected gradient is never directed opposite to the
unconstrained gradient.

This leaves us with three potential variational methods: optimising $J_{SV}$
(SV), optimising $J_{PV}$ (PV), or optimising $J_{PV}$ with Leray projection
of the gradient at each step to enforce the divergence-free condition (PV-LP).
To compare the performance of these methods, we discuss the convergence
starting from two initial candidates, which were found through a recurrent
flow analysis (see section 4.1). Figure 1 depicts a state-space projection of
these two candidates as well as the converged solutions, and figure 2 shows
the flow structure at single snapshots in time. The first candidate converged
to a simple periodic orbit of period $T=5.38$, given as ‘P1’ in Chandler &
Kerswell (2013), and the second to the relative periodic orbit ‘R19’ with
$T=12.2$, although the algorithm has converged to two periods of this orbit,
to give $T=24.4$. Note how vastly more complicated the initial guess is than
the converged solution in this latter case, at least in this simple
projection, hinting at the ability of the variational methods to converge to
solutions from very distant initial guesses.

Figure 3 shows the convergence of the relevant residual against the number of
conjugate-gradient iterations, for the first of these two candidates. The
algorithms were each run for three days of wall-clock compute time on a 28
core CPU. Due to the line search requiring more evaluations, each SV iteration
is on average around 40 percent slower than the two PV methods, which leads to
fewer iterations in 72h. More importantly, SV shows a significantly slower
convergence rate. All three methods show an initial fast improvement of the
objective function, followed by a much slower period, consistent with Azimi et
al. (2020); the convergence rate in this latter region is much faster with the
PV methods. The precise reason for this is unclear, and we make no claim that
our conjugate-gradient algorithm uses optimal parameters, but this result was
robust after significant trial-and-error tweaking of them. Whether Leray
projection is used has negligible impact on the rate of convergence, but the
path taken towards the solution is certainly different as PV allows for
compressible intermediates. Consequently, the methods may converge to
different solutions.

Figure 1: Projections of the guess (thin) and converged (thick) loops for the
two examples discussed. The projections $P_{1}$ and $P_{2}$ are defined to be,
respectively, the imaginary parts of the $1,1$ the $2,2$ spatial Fourier modes
of vorticity. The labels correspond to the snapshots shown in figure 2.

\begin{overpic}[width=85.35826pt,angle={90},origin=c]{data/c5guess.png}
\put(10.0,80.0){(a)} \end{overpic}
\begin{overpic}[width=85.35826pt,angle={90},origin=c]{data/c5exact.png}
\put(10.0,80.0){(b)} \end{overpic}
\begin{overpic}[width=85.35826pt,angle={90},origin=c]{data/c68guess.png}
\put(10.0,80.0){(c)} \end{overpic}
\begin{overpic}[width=85.35826pt,angle={90},origin=c]{data/c68exact.png}
\put(10.0,80.0){(d)} \end{overpic}

Figure 2: Slices at $s=0$ of the vorticity field for the guess (a,c) and
converged (b,d) loops discussed in the text, at the snapshots labelled in
figure 1. Movies are available in the supplementary materials. Figure 3:
Convergence of the respective residual $J_{PV}$ or $J_{SV}$ for the primitive
variables formulation, with (PV-LP) and without (PV) projection and the
streamfunction-voriticity (SV) methods applied to the first guess in figure 1.
In the SV case, $J_{SV}$ starts at a significantly larger value, but also
converges much more slowly.

To investigate the range of convergence, we systematically vary the initial
guess by linearly interpolating between the converged solution and the
candidate extracted from flow recurrences,
$u_{\mathrm{initial}}=(1-\alpha)u_{\mathrm{solution}}+\alpha
u_{\mathrm{candidate}}$, with $\alpha=0$ being the exact (relative) periodic
orbit (P1 and R19) and $\alpha=1$ being the candidate. As shown on the top-
left in figure 4, for the first candidate, the PV formulation, both with and
without Leray projection, shows good convergence over the full range of
$\alpha$, with $J_{PV}$ consistently reaching less than $10^{-13}$. However,
extrapolating to $\alpha>1$, where the initial guess is very different from
the periodic orbit P1, the converged period $T$ is shown to rapidly vary. This
is because, in this case, the algorithm is no longer converging to P1 but onto
a travelling wave solution, for which the period is indeterminate. Note that
around $\alpha=1.2$, the Leray-projected method does converge onto P1 whereas
without it does not, which confirms that Leray projection does indeed have a
significant role, and that without it, the PV formulation loses the
divergence-free property of the guess in the course of converging before again
becoming divergence-free as the solution converged, a fact which was confirmed
when the divergence was examined. In contrast with the primitive-variables
formulations, the SV formulation performed poorly, with the residual not
reaching an acceptably low level within the allotted time for all but the
smallest values of $\alpha$, and the final value $T$ being noticeably wrong
for $\alpha>0.3$, without the rapid fluctuations indicative of convergence
onto a different solution.

For the second, longer period candidate, the results are similar, with
convergence good for the primitive variables formulations but much worse for
the streamfunction-vorticity formulation. However, in this case, for
$\alpha\gtrsim 1.3$, the primitive variables formulations appear to converge
to a consistent value of $T$ and a consistent value of $J_{PV}$, which is
considerably greater than zero. This, therefore, is a local minimum of the
objective functional. This demonstrates the importance of ensuring that the
results of the algorithms we find are genuine solutions, rather than just
local minima.

Figure 4: Comparison of methods applied to an initial guess linearly
interpolated/extrapolated between the converged periodic orbit ($\alpha=0$)
and the recurrent flow analysis candidate ($\alpha=1$). The top plots show the
final value of the residual $J_{PV}$, where a value of less than $10^{-10}$ is
well converged. The lower plots show final value of the period $T$. Left and
right as per the candidates in figure 1.

###  4.1 Recurrent flow analysis and comparison with Newton-based shooting

As the streamfunction-vorticity formulation is signficiantly slower and has
worse convergence properties than the primitive variables formulation, in this
section we concentrate on the primitive variables formulations, with (PV-LP)
and without (PV) Leray projection and compare its performance to a traditional
shooting method.

To test the methods, we implemented a recurrent flow analysis to generate
candidates for exact coherent solutions. Our implementation follows that of
Chandler & Kerswell (2013), though the precise details, such as how many
timesteps and phase shifts to minimise over, may differ. Subsequent authors
(Lucas & Kerswell, 2015; Page & Kerswell, 2020) have improved upon this method
to find more candidates, but the focus of the present work is merely to
investigate the properties of the variational methods. We compared against a
Newton shooting (SN) method. We used the Newton solver from Channelflow
(channelflow.ch) (Gibson et al., 2021), with the maximum number of Newton
steps set to 100. We made no attempt to search for relative periodic orbits
which exhibit a shift in the $y$ direction, exploiting the discrete symmetry
of the equations. There is no reason to suspect that including this would
affect the methods we study.

A time series of length $200$ was searched for near-recurrences up to period
$T=50$, giving 106 candidate solutions. Table 1 summarises the results from
trying to converge each of these candidates via PV, PV-LP and SN. We counted
as converged any solution for which the residual $J_{PV}$ was less than
$10^{-8}$, as the known local minima all had $J_{PV}>10^{-6}$. For all cases
except those with particularly long periods, in practice the final value of
$J_{PV}$ was significantly smaller than $10^{-8}$. All three methods were
successful at converging solutions, including several that had been found
previously and a number of new solutions. Even with this relatively small
sample, it is clear that our new variational methods converge to a greater
variety of solutions within the chaotic attractor than the traditional
shooting method. Instances were found where the Newton method converged for a
given candidate and the variational methods did not, and vice-versa.

Both our variational methods and the shooting method tend to struggle more in
converging orbits with longer periods compared with short ones, though the
reasons for this differ. As the shooting method time-marches the orbit,
exponential error amplification leads to sensitive dependence of the
recurrence condition on initial conditions and thus an ill-conditioned root
search. It is for this reason that multi-shooting methods have been developed
(Sánchez & Net, 2010). In our variational methods, the loop structure means
that recurrence is ‘built in’. Nevertheless we observe that the algorithm
converges more slowly for longer orbits. Indeed, a large number of candidates
appeared to be converging to long-period orbits when the algorithms were
terminated (after 3 days of wall-clock time) and may well have converged had
they been allowed to continue. In other cases, definite local minima were
found by the variational methods, at which point convergence is no longer
possible. However, switching to a different variational method (e.g. from PV
to PV-LP) at this point was found to subsequently lead to convergence to a
solution.

$c$ | $T$ | Name | Number converged
---|---|---|---
|  |  | PV | PV-LP | SN
0 | - | E1 | 7 | 7 | 7
0.0198 | - | T1 | 26 | 24 | 22
0 | 5.38 | P1 | 9 | 11 | 3
0 | 2.83 | P2 | 2 | 4 | 1
0 | 2.92 | P3 | 2 | 0 | 0
0.0352 | 12.2 | R19 | 2 | 2 | 0
0.0173 | 36.8 | R47 | 2 | 3 | 0
0.00446 | 7.16 | - | 0 | 0 | 1
0.0213 | 18.1 | - | 0 | 0 | 2

$c$ | $T$ | Name | Number converged
---|---|---|---
|  |  | PV | PV-LP | SN
0.00243 | 3.78 | - | 1 | 0 | 0
0.0472 | 8.46 | - | 1 | 0 | 0
0.00780 | 8.64 | - | 1 | 0 | 0
0.00352 | 9.62 | - | 1 | 2 | 0
0.0106 | 12.5 | - | 1 | 1 | 0
0.00993 | 15.7 | - | 1 | 1 | 0
0.0208 | 18.6 | - | 1 | 0 | 0
Total | 57 | 55 | 36
Unique | 14 | 9 | 6

Table 1: Results of different methods applied to candidates from the recurrent
flow analysis. The names are from Chandler & Kerswell (2013), for those
solutions which are not new. All values are given to three significant
figures.

##  5 Conclusions

In this paper, we have presented the first use of variational methods to find
unstable periodic orbits in the incompressible Navier-Stokes equations. Though
the algorithms were slow to converge compared with existing techniques, we
have demonstrated their superior radius of convergence, raising the
possibility of converging a larger number and a greater variety of exact
solutions in turbulent flows.

We showed that the primitive variables formulation performs notably better
than the streamfunction-vorticity formulation, despite the fact that the
latter automatically enforces incompressibility and eliminates the pressure
variable. Projecting onto the incompressible subspace or including the
divergence-free condition in the objective functional makes little practical
difference to the algorithm, suggesting that including incompressibility in
the cost function but allowing for compressible intermediate fields during
minimization is a viable approach.

The natural extension to this work, and an important step to understanding
turbulence in practical flows, is to consider three-dimensional flow models
with boundaries, such as plane Couette flow or plane Poiseuille flow. This
presents a number of issues not considered here. Firstly, with a third
dimension at moderately high Reynolds number, the number of grid points
required would be several orders of magnitude higher. However, all
computations in this work were performed on a single CPU with OpenMP
paralellisation, and with MPI paralellisation over many CPUs, or indeed GPU
parallelisation, it should be possible to apply such methods to much higher
dimensional systems, since all operations except the Fourier transforms are
local in either spectral or physical space. Secondly, the presence of
boundaries significantly complicates the derivation of gradients of the
objective functionals, and care would need to be taken in particular for
boundary conditions of the pressure field, which can be a complicated issue
(Gresho & Sani, 1987). Thirdly, the streamfunction-vorticity formulation does
not have a simple analogue in three dimensional flow, though since it
performed poorly here, this is a minor consideration.

In section 4.1 we showed that our method is able to converge more solutions
than a shooting method, but in some cases the shooting method converged when
ours did not. This suggests that the best strategy, if the aim is to converge
as many different solutions as possible, would be to run both methods in
parallel. Likewise, a hybrid approach combining both methods may alleviate the
slow convergence of our method by switching to a Newton-based shooting method
as soon as consistent convergence of the variational method is observed.

## References

  * Azimi et al. (2020) Azimi, S, Ashtari, O & Schneider, TM 2020 Adjoint-based variational method for constructing periodic orbits of high-dimensional chaotic systems. arXiv:2007.06427 .
  * Boghosian et al. (2011) Boghosian, Bruce M, Fazendeiro, Luis M, Laett, Jonas, Tang, Hui & Coveney, Peter V 2011 New variational principles for locating periodic orbits of differential equations. Philos. Trans. R. Soc. A 369 (1944), 2211–2218.
  * Chandler & Kerswell (2013) Chandler, GJ & Kerswell, RR 2013 Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow. J. Fluid Mech. 722, 554–595.
  * Cvitanović et al. (2016) Cvitanović, P., Artuso, R., Mainieri, R., Tanner, G. & Vattay, G. 2016 Chaos: Classical and Quantum. Copenhagen: Niels Bohr Inst.
  * Dijkstra et al. (2014) Dijkstra, HA, Wubs, FW, Cliffe, AK, Doedel, E, Dragomirescu, IF, Eckhardt, B, Gelfgat, AY, Hazel, AL, Lucarini, V, Salinger, AG & others 2014 Numerical bifurcation methods and their application to fluid dynamics: analysis beyond simulation. Commun. Comput. Phys. 15 (1), 1–45.
  * Farazmand (2016) Farazmand, M 2016 An adjoint-based approach for finding invariant solutions of Navier–Stokes equations. J. Fluid Mech. 795, 278–312.
  * Fletcher & Reeves (1964) Fletcher, R & Reeves, CM 1964 Function minimization by conjugate gradients. The computer journal 7 (2), 149–154.
  * Fylladitakis et al. (2018) Fylladitakis, ED & others 2018 Kolmogorov flow: Seven decades of history. J. Appl. Math. Phys. 6 (11), 2227\.
  * Gibson et al. (2021) Gibson, JF, Reetz, F, Azimi, S, Ferraro, A, Kreilos, T, Schrobsdorff, H, Farano, M, Yesil, AF, Schütz, SS, Culpo, M & Schneider, TM 2021 Channelflow 2.0. in preparation .
  * Gresho & Sani (1987) Gresho, PM & Sani, RL 1987 On pressure boundary conditions for the incompressible journalStokes equations. Int. J. Numer. Methods Fluids 7 (10), 1111–1145.
  * Kawahara & Kida (2001) Kawahara, G & Kida, S 2001 Periodic motion embedded in plane Couette turbulence: regeneration cycle and burst. J. Fluid Mech. 449, 291–300.
  * Kawahara et al. (2012) Kawahara, G, Uhlmann, M & Van Veen, L 2012 The significance of simple invariant solutions in turbulent flows. Annu. Rev. Fluid Mech. 44, 203–225.
  * Lan & Cvitanović (2004) Lan, Y & Cvitanović, P 2004 Variational method for finding periodic orbits in a general flow. Phys. Rev. E 69 (1), 016217.
  * Lan & Cvitanović (2008) Lan, Y & Cvitanović, P 2008 Unstable recurrent patterns in Kuramoto-Sivashinsky dynamics. Phys. Rev. E 78 (2), 026208.
  * Lucas (2020) Lucas, D 2020 Stabilisation of exact coherent structures in two-dimensional turbulence using time-delayed feedback. arXiv:2008.08388 .
  * Lucas & Kerswell (2015) Lucas, D & Kerswell, RR 2015 Recurrent flow analysis in spatiotemporally chaotic 2-dimensional Kolmogorov flow. Phys. Fluids 27 (4), 045106.
  * Meshalkin & Sinai (1961) Meshalkin, LD & Sinai, Ia G 1961 Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid. J. Appl. Math. Phys. 25 (6), 1700–1705.
  * Obukhov (1983) Obukhov, AM 1983 Kolmogorov flow and laboratory simulation of it. Russ. Math. Surv 38 (4), 113–126.
  * Page & Kerswell (2020) Page, J & Kerswell, RR 2020 Searching turbulence for periodic orbits with dynamic mode decomposition. J. Fluid Mech. 886.
  * Platt et al. (1991) Platt, Nathan, Sirovich, L & Fitzmaurice, N 1991 An investigation of chaotic Kolmogorov flows. Phys. Fluids A 3 (4), 681–696.
  * Sánchez & Net (2010) Sánchez, J & Net, M 2010 On the multiple shooting continuation of periodic orbits by Newton–Krylov methods. Int. J. Bifurc. Chaos. 20 (01), 43–61.
  * Schneider et al. (2010) Schneider, Tobias M, Gibson, John F & Burke, John 2010 Snakes and ladders: localized solutions of plane Couette flow. Phys. Rev. Lett. 104 (10), 104501.
  * Temam (2001) Temam, Roger 2001 Navier-Stokes equations: theory and numerical analysis. American Mathematical Soc.
  * Tuckerman (2020) Tuckerman, Laurette S 2020 Computational challenges of nonlinear systems. In Emerging Frontiers in Nonlinear Science, pp. 249–277. Springer.
  * van Veen et al. (2019) van Veen, L, Vela-Martín, A, Kawahara, G & Yasuda, T 2019 Periodic orbits in large eddy simulation of box turbulence. Fluid Dynamics Research 51 (1), 011411.
  * Viswanath (2007) Viswanath, D 2007 Recurrent motions within plane Couette turbulence. J. Fluid Mech. 580, 339–358.
  * Wolfe (1969) Wolfe, P 1969 Convergence conditions for ascent methods. SIAM review 11 (2), 226–235.