.tiffpng.pngconvert #1 .tiff

# Heat conduction in an irregular magnetic field: Part II. Heat transport as a
measure of the effective non-integrable volume

Elizabeth J. Paul1,2<EMAIL_ADDRESS>Stuart R. Hudson2  Per Helander3
1Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA 2Princeton Plasma Physics Laboratory, PO Box 451, Princeton NJ
08543, USA 3Max Planck Institute for Plasma Physics, Greifswald, Germany

###### Abstract

Given the large anisotropy of transport processes in magnetized plasmas, the
magnetic field structure can strongly impact heat diffusion: magnetic surfaces
and cantori form barriers to transport while chaotic layers and island
structures can degrade confinement. When a small but non-zero amount of
perpendicular diffusion is included, the structure of the magnetic field
becomes less important, allowing pressure gradients to be supported across
chaotic regions and island chains. We introduce a metric for the effective
volume over which the local parallel diffusion dominates based on the solution
to the anisotropic heat diffusion equation. To validate this metric, we
consider model fields with a single island chain and a strongly chaotic layer
for which analytic predictions of the relative parallel and perpendicular
transport can be made. We also analyze critically chaotic fields produced from
different sets of perturbations, highlighting the impact of the mode number
spectrum on the heat transport.  Our results indicate that this metric
coincides with the effective volume of non-integrability in the limit
$\kappa_{\perp}\rightarrow 0$. We propose that this metric be used to assess
the impact of non-integrability on the heat transport in stellarator
equilibria.

##  1 Introduction

The heat transport in irregular magnetic fields is critical for understanding
magnetic confinement experiments and astrophysical systems. Tokamak
experiments exploit control coils which produce resonant magnetic
perturbations for control of edge localized modes (ELMs) (Evans et al., 2005).
Such perturbations must induce a small amount of stochasticization near an
edge resonant surface to suppress ELMs without significant degradation of the
core confinement. These coils can impart structure in the field outside the
last-closed flux surface, resulting in homoclinic tangles (Evans et al., 2005)
and stochastic layers (Punjabi et al., 2008) which can impact the heat flux
deposition at the divertor. Reversed Field Pinches (RFPs) provides confinement
through magnetic fields generated largely by the plasma current. The transport
in an RFP is often dominated by parallel diffusion due to stochastic magnetic
fields which form due to tearing mode activity (Sarff et al., 2013). The
presence of a transport barrier (Biewer et al., 2003; Lorenzini et al., 2009)
can enable the support of a temperature gradient. As stellarators provide
confinement with three-dimensional magnetic fields, the impact of non-
integrability is critical for heat transport in both the core and divertor
regions. For example, the stochastic edge of LHD features island remnants,
stochastic layers, and partial barriers. Due to the long connection lengths,
both the parallel and perpendicular transport are important for determining
the heat flux impacting the divertor plates in a stellarator (Feng et al.,
2011). Stochastic magnetic fields are also expected to exist in several
astrophysical environments. For example, the stochasticity of the galactic
magnetic field may impact the transport of cosmic rays (Jokipii & Parker,
1969; Chuvilgin & Ptuskin, 1993).

Given the Hamiltonian nature of magnetic field line flow (Cary & Littlejohn,
1983), there are a number of tools from the field of dynamical systems that
enable the quantification of the transition to chaos. At the most granular
level, the linear stability of a single resonance is quantified through
Greene’s residue (Greene, 1979), $R$, where $R<0$ or $R>1$ indicates unstable
(hyperbolic or hyperbolic-with-reflection) orbits and $R\in[0,1]$ indicates
stable (elliptic) orbits. Thus one technique to reduce the stochasticity and
island size in a stellarator is through minimization of the magnitude of the
residue (Hanson & Cary, 1984; Cary & Hanson, 1986). Similarly, the island
width under the small island approximation can be evaluated through a
linearization about the island center (Cary & Hanson, 1991). A shortcoming of
these approaches is their limitation to a single resonance rather than a
volumetric measure of non-integrability.

One can also evaluate the existence or non-existence of a KAM surface to
determine the transition to chaos in a given region. To estimate the value of
the perturbation required to destroy the last KAM surface between two primary
resonances, one may compute the Chirikov overlap criterion (Chirikov, 1979),
obtained from the overlap of resonances and their second order island chains.
Greene has also conjectured a procedure for determining the existence of a KAM
surface by evaluating the residues of the high-order rational convergents of a
given irrational. By examining the stability of the neighboring periodic
orbits to the most stable KAM surface, one can determine the transition to
stochasticity (Greene, 1979; MacKay, 1992; Falcolini & de La Llave, 1992). In
a similar way, one can evaluate Mather’s $\Delta W$, the difference in action
between the X and O points associated with a given rational rotation number,
$p/q$. If $\Delta W\rightarrow 0$ as $p/q$ approaches a given irrational, then
the corresponding KAM surface survives (Mather, 1986). Furthermore, converse
KAM theory (Mather, 1984; MacKay & Percival, 1985; MacKay, 2018) provides a
sufficient condition for non-existence of a KAM surface by considering the
rotation of nearby trajectories with respect to a given foliation.

While each of these techniques may shed light on the existence or non-
existence of a KAM surface in a given region, another set of techniques are
required to provide insight into the relative impact of non-integrability on
the transport in a dynamical system (Meiss, 2015). For example, Mather’s
$\Delta W$ is related to the flux of trajectories across a given curve in
phase space, providing a bound on the transport (MacKay et al., 1984). One can
then partition phase space into subvolumes bounded by KAM surfaces or cantori
and model transport between states, such as with a Markov model (Meiss & Ott,
1986).

Rather than using nonlinear dynamical surrogates of transport, one could
directly evaluate the transport in a given magnetic field line flow. A logical
tactic for studying the transport is by introducing a tracer, such as the
temperature, $T(\textbf{x})$, which diffuses along the magnetic field subject
to some initial condition $T_{0}(\textbf{x})$,

$\displaystyle\frac{\partial T(\textbf{x},t)}{\partial
t}-\nabla\cdot\left(\kappa_{\|}\hat{\textbf{b}}\hat{\textbf{b}}\cdot\nabla
T\right)$ $\displaystyle=0$ (1a) $\displaystyle T(\textbf{x},0)$
$\displaystyle=T_{0}(\textbf{x}).$ (1b)

The temperature at a given point, $T(\textbf{x},t)$, can be computed by
integrating over all contributions of the initial condition along the magnetic
field line path (del Castillo-Negrete & Chacon, 2012). In the steady-state
$t\rightarrow\infty$ limit, the temperature becomes constant on magnetic field
lines; it is simply the average of the initial condition along the magnetic
field line path,

$\displaystyle
T(\textbf{x},t\rightarrow\infty)=\lim_{L\rightarrow\infty}\frac{\int_{0}^{L}dl\,T_{0}(\textbf{x}(l))}{\int_{0}^{L}dl},$
(2)

where $l$ measures length along a field line. Because the temperature is
constant along field lines, in a chaotic magnetic field this leads to a
fractal-like structure in the temperature profile (Hudson & Breslau, 2008; del
Castillo-Negrete & Chacon, 2012). In this way, local flattening of the
temperature profile can be used as a measure of the local stochasticity.

While the structure of the magnetic field itself can give rise to transport,
in a plasma a non-zero amount of perpendicular transport exists due to the
gyromotion, cross-field drifts, collisions, and fluctuations. Although the
effective perpendicular diffusion can be extremely small in a magnetic
confinement device (experimental measurements on the TEXTOR, DIII-D, and ASDEX
tokamaks indicate $\kappa_{\perp}/\kappa_{\|}\sim 10^{-9}-10^{-7}$ (Meskat et
al., 2001; Hölzl et al., 2009; Bardóczi et al., 2016), where $\kappa_{\perp}$
is the perpendicular diffusion coefficient and $\kappa_{\|}$ is the parallel
diffusion coefficient), non-zero perpendicular diffusion can substantially
alter the heat transport. For example, in the Large Helical Device (LHD) a
$\beta=4.8\%$ discharge demonstrated a 32% reduction in nested magnetic
surfaces in comparison with the vacuum magnetic field. However, significant
flattening of the temperature profile was not observed (Sakakibara et al.,
2008), indicating the relative importance of the perpendicular transport. The
perpendicular diffusion becomes even more critical in the stochastic edge,
where the collisionality is enhanced.

In the steady-state limit, the transport equation (1) can be modified to
include the effect of perpendicular diffusion,

$\displaystyle\nabla\cdot\left(\kappa_{\|}\nabla_{\|}T+\kappa_{\perp}\nabla_{\perp}T\right)$
$\displaystyle=0,$ (3)

where the parallel gradient is
$\nabla_{\|}=\hat{\textbf{b}}\hat{\textbf{b}}\cdot\nabla$ and the
perpendicular gradient is $\nabla_{\perp}=\nabla-\nabla_{\|}$. We will discuss
important properties of this anisotropic diffusion equation (ADE) in §2. The
relative importance of the structure of the magnetic field on transport can
then be quantified by comparing the parallel diffusion,

$\displaystyle
D_{\|}=\frac{1}{2}\int_{V}d^{3}x\,\kappa_{\|}\left|\nabla_{\|}T\right|^{2},$
(4)

with the perpendicular diffusion,

$\displaystyle
D_{\perp}=\frac{1}{2}\int_{V}d^{3}x\,\kappa_{\perp}\left|\nabla_{\perp}T\right|^{2},$
(5)

where $V$ is the volume of the domain of interest. In this way, the non-
integrability of the magnetic field can be quantified by its impact on the
transport in a volume rather than through the analysis of individual
resonances or the existence of a given KAM surface. This allows us to define
an effective volume of perpendicular diffusion in §3,  which we conjecture
converges to the effective volume of non-integrability in the limit
$\kappa_{\perp}\rightarrow 0$. In Section 4 we present the numerical methods
used to solve the diffusion equation and the model field used for our
calculations. In Sections 5 and 6 we perform calculations of the transport in
the presence of a single island chain and a strongly chaotic layer to justify
our metric. We then apply the metric to a critical chaotic layer created from
resonances of different mode numbers in §7. We then conclude in §8.

##  2 Properties of the anisotropic diffusion equation

The anisotropic diffusion equation (ADE) can be expressed as,

$\displaystyle\nabla\cdot\left(\bm{\kappa}\cdot\nabla T\right)$
$\displaystyle=0,$ (6)

with the diffusion tensor,

$\displaystyle\bm{\kappa}=\kappa_{\perp}\textbf{I}+(\kappa_{\|}-\kappa_{\perp})\hat{\textbf{b}}\hat{\textbf{b}},$
(7)

where the perpendicular, $\kappa_{\perp}$, and parallel, $\kappa_{\|}$,
diffusion coefficients are taken to be constants. We compute solutions to (6)
in an annular volume $\Omega$ bounded by two toroidal surfaces, $S_{-}$ and
$S_{+}$ with constant Dirichlet boundary conditions $T_{-}$ and $T_{+}$.

###  2.1 Ellipticity

We note that $\bm{\kappa}$ is positive definite for $\kappa_{\perp}>0$ and
$\kappa_{\|}\geq\kappa_{\perp}$ since,

$\displaystyle\nabla T\cdot\bm{\kappa}\cdot\nabla T=\kappa_{\perp}|\nabla
T|^{2}+(\kappa_{\|}-\kappa_{\perp})\left(\hat{\textbf{b}}\cdot\nabla
T\right)^{2}>0,$ (8)

for all $\nabla T\neq 0$.  (While the assumption
$\kappa_{\|}\geq\kappa_{\perp}$ is not necessary for ellipticity, it is
sufficient according to (8) and not too restrictive since we are interested in
strong anisotropy, $\kappa_{\|}\gg\kappa_{\perp}$.) We can furthermore see
that the anisotropic diffusion operator is elliptic, as it can be expressed
as,

$\displaystyle\sum_{i,j}\partial_{i}\left(\kappa_{ij}\partial_{j}T\right)=\sum_{i,j}\kappa_{ij}\partial_{i}\partial_{j}T+\left(\partial_{i}\kappa_{ij}\right)\partial_{j}T.$
(9)

Ellipticity of the anisotropic diffusion operator (6) follows from the
symmetric positive-definiteness of the tensor which multiplies the second
derivatives, $\kappa_{ij}$. Ellipticity implies the existence of a maximum
principle (Appendix A),

$\displaystyle\max_{\overline{\Omega}}T=\max_{\partial\Omega}T.$ (10)

The corresponding statement also holds for the minimum. In words, the maximum
(minimum) temperature must be obtained on the boundary. Furthermore, no local
minima or maxima can occur within $\Omega$. As a further consequence of
ellipticity, the topology of the isotherms is constrained (Appendix B).
Specifically, an isotherm cannot enclose a net volume, but instead must link
the domain toroidal and poloidally as the boundaries (with the possible
addition of a handle).

###  2.2 Variational principle

The anisotropic diffusion equation can be obtained (Helander et al., 2021)
from a variational principle involving the functional,

$\displaystyle\mathcal{W}[T]=\frac{1}{2}\int_{\Omega}d^{3}x\,\nabla
T\cdot\bm{\kappa}\cdot\nabla T.$ (11)

The first variation of $\mathcal{W}$ with respect to $T$ is computed to be,

$\displaystyle\delta\mathcal{W}[\delta T]$
$\displaystyle=\int_{\Omega}d^{3}x\,\nabla\delta T\cdot\bm{\kappa}\cdot\nabla
T=-\int_{\Omega}d^{3}x\,\delta T\nabla\cdot\left(\bm{\kappa}\cdot\nabla
T\right),$ (12)

upon application of the boundary condition $\delta
T\rvert_{\partial\Omega}=0$. Thus stationary points of $\mathcal{W}$ with
respect to $T$ correspond to solutions of (6). We can now compute the second
variations as,

$\displaystyle\delta^{2}\mathcal{W}[\delta T,\delta
T^{\prime}]=\int_{\Omega}d^{3}x\,\nabla\delta
T\cdot\bm{\kappa}\cdot\nabla\delta T^{\prime}.$ (13)

We note that $\delta^{2}\mathcal{W}[\delta T,\delta T]\geq 0$ under the
assumption that $\kappa_{\perp}<\kappa_{\|}$. Therefore, $\mathcal{W}$ is a
convex function, implying that any stationary point of $\mathcal{W}$ is a
global minimum. We will refer to the functional $\mathcal{W}$ that appears in
this variational principle as the diffusion integral,

$\displaystyle
D=\frac{1}{2}\int_{\Omega}d^{3}x\,\left(\kappa_{\|}|\nabla_{\|}T|^{2}+\kappa_{\perp}|\nabla_{\perp}T|^{2}\right),$
(14)

as it is related to the entropy production functional (Hameiri &
Bhattacharjee, 1987). Thus the temperature profile minimizes the entropy
production subject to the boundary conditions.

We note that the total diffusion is related to the heat flux through the
boundaries,

$\displaystyle
Q=-\int_{S_{+}}d^{2}x\,\hat{\textbf{n}}\cdot\textbf{q}=-\int_{S_{-}}d^{2}x\,\hat{\textbf{n}}\cdot\textbf{q},$
(15)

through

$\displaystyle Q=\frac{2D[T]}{T_{+}-T_{-}},$ (16)

where $\textbf{q}=-\kappa_{\|}\nabla_{\|}T-\kappa_{\perp}\nabla_{\perp}T$ and
$\hat{\textbf{n}}$ is the unit normal oriented in the direction from $S_{-}$
to $S_{+}$. This relation follows from,

$\displaystyle 0$ $\displaystyle=\int_{\Omega}d^{3}x\,T\nabla\cdot{\bf
q}=-\int_{\Omega}d^{3}x\,{\bf q}\cdot\nabla
T+\int_{S_{+}}d^{2}x\,T\hat{\textbf{n}}\cdot{\bf
q}-\int_{S_{-}}d^{2}x\,T\hat{\textbf{n}}\cdot{\bf q}$ (17)
$\displaystyle=2D[T]-Q(T_{+}-T_{-}).$

Thus the total diffusion quantifies the heat flux required to support a given
temperature gradient across a volume. In this work we will consider $T_{+}$
and $T_{-}$ to be fixed such that $D$ measures the required heat flux to
support the prescribed temperature gradient. Since the perpendicular transport
is relatively insensitive to the structure of the field, we can assess the
impact of non-integrability on the transport through the local parallel
diffusion,

$\displaystyle d_{\|}=\frac{1}{2}\kappa_{\|}|\nabla_{\|}T|^{2},$ (18)

in comparison with the local perpendicular diffusion,

$\displaystyle d_{\perp}=\frac{1}{2}\kappa_{\perp}|\nabla_{\perp}T|^{2}.$ (19)

##  3 The effective volume of parallel diffusion

Motivated by comparing $d_{\|}$ to $d_{\perp}$ in order to determine the
impact of non-integrability on transport, we propose the following metric:

$\displaystyle\mathcal{V}_{\text{PD}}=\frac{\int_{\Omega}d^{3}x\,\Theta\left(\kappa_{\|}|\nabla_{\|}T|^{2}-\kappa_{\perp}|\nabla_{\perp}T|^{2}\right)}{\int_{\Omega}d^{3}x},$
(20)

where $\Theta$ is the Heaviside step function. In words,
$\mathcal{V}_{\text{PD}}$ measures the fraction of volume over which the local
parallel transport is larger than the perpendicular transport. While one could
imagine many ways in which to compare the parallel and perpendicular
diffusion, the above definition is easy to interpret as it is a positive
scalar in $[0,1]$.

While the definition of $\mathcal{V}_{\text{PD}}$ is not intrinsic to the
field structure as it depends on $\kappa_{\perp}/\kappa_{\|}$, we can infer
general trends in the limit of small $\kappa_{\perp}/\kappa_{\|}$. In the
limit of an integrable magnetic field with sufficiently small
$\kappa_{\perp}$, the isotherms will largely coincide with magnetic surfaces,
$|\hat{\textbf{b}}\cdot\nabla T|$ will be relatively small, and thus the local
parallel transport will be negligible. Thus $\mathcal{V}_{\text{PD}}$ will be
small in the limit of integrability. Since the topology of the isotherms is
constrained, isotherms will not be able to completely align to the structure
of the field in regions of non-integrability, and
$|\hat{\textbf{b}}\cdot\nabla T|$ will be consequently larger, as measured by
$\mathcal{V}_{\text{PD}}$.

We can consider the metric $\mathcal{V}_{\text{PD}}$ to be a “diffusive”
analog of a metric based on converse KAM theory. With a converse KAM approach,
the non-existence of an invariant torus at a given point in phase space is
determined from the rotation of nearby trajectories with respect to a chosen
foliation. This allows one to determine regions of phase space over which
invariant surfaces do not exist with the topology of the chosen foliation,
i.e. a given class of tori (MacKay & Percival, 1985; MacKay, 2018). Similarly,
by comparing the local parallel and perpendicular diffusion, we can assess
whether a given region of physical space is “effectively non-integrable”
(indicated by enhanced parallel diffusion) with respect to its impact on the
transport. This assessment similarly depends on the topology of the isotherms,
which foliate the volume.  Since the isotherms must link the domain as the
boundaries do, the ADE similarly provides an assessment of the non-
integrability with respect to a limited class of tori. For instance, an island
chain will impact the parallel transport if the boundary surfaces are chosen
to have the same topology as the unperturbed Hamiltonian,  given that
isotherms cannot have the same linking as the island. If instead the boundary
surfaces have the same linking as the surfaces within the island chain, the
parallel transport will not be impacted.

To assess this metric, we consider two cases for which analytic predictions
for the parallel transport can be employed: diffusion across a single island
chain (§5) and across a strongly chaotic layer (§6). In §7 we then consider a
model field for which such prediction is more challenging, namely a critically
chaotic field with remnants of island chains and partial barriers.

##  4 Numerical methods

###  4.1 Linear solution

We compute the numerical solution to the anisotropic diffusion equation
(6)-(7) using a Fourier Galerkin discretization in $\theta$ and $\zeta$ and a
fourth-order finite-difference radial discretization. The resulting sparse
linear system is solved with the PETSc library (Balay et al., 2021b, a, 1997)
using the MUMPS package (Amestoy et al., 2001, 2019) to compute the
$LU$-factorization. Throughout we will normalize calculations such that
$T_{-}=0$, $T_{+}=1$, and $\kappa_{\|}=1$.

###  4.2 Model magnetic field

We consider a model Hamiltonian,

$\displaystyle\chi(\psi,\theta,\zeta)=\frac{\iota^{\prime}\psi^{2}}{2}+\sum_{m,n}\epsilon_{m,n}\psi(\psi-\overline{\psi})\cos(m\theta-n\zeta),$
(21)

which yields the magnetic field,

$\displaystyle\textbf{B}=\nabla\psi\times\nabla\theta-\nabla\chi\times\nabla\zeta$
(22)

where $\psi=\overline{\psi}\rho/\overline{\rho}$ and $(\rho,\theta,\zeta)$
forms an orthogonal coordinate system
($\hat{\bm{\rho}}\times\hat{\bm{\theta}}\cdot\hat{\bm{\zeta}}=1$) with
$\rho\in[0,\overline{\rho}]$, $\theta\in[0,2\pi)$, $\zeta\in[0,2\pi)$,
$|\nabla\rho|=1$, $|\nabla\theta|=1/\overline{\rho}$, and
$|\nabla\zeta|=1/L_{\zeta}$. We will use the notation
$\widetilde{\epsilon}_{m,n}(\psi)=\epsilon_{m,n}\psi(\psi-\overline{\psi})$ to
denote the strength of the perturbation modes. In the limit
$\epsilon_{m,n}\rightarrow 0$ we recover an integrable magnetic field lying on
surfaces of constant $\rho$ with rotational transform
$\iota=\iota^{\prime}\psi$ and constant shear
$\iota^{\prime}(\psi)=\iota^{\prime}$. The perturbation is chosen to vanish at
the upper ($\rho=\overline{\rho}$) and lower ($\rho=0$) boundaries such that
there is no parallel flux into the volume. For all of the following numerical
calculations, we will take
$\overline{\psi}=\overline{\rho}=L_{\zeta}=\iota^{\prime}=1$.

##  5 Transport across an island chain

In this Section, we consider the relative perpendicular and parallel transport
in the presence of a single island chain generated by a resonance $\iota=n/m$.
We consider the case in which $\kappa_{\perp}$ remains small but non-zero such
that perpendicular diffusion competes with parallel diffusion.

To analyze the behavior of the transport in such a field, in Appendix C we
perform a perturbation analysis in the smallness of the amplitude
$\epsilon_{m,n}$. Due to the smallness of $\kappa_{\perp}$, outside of a small
boundary layer of width

$\displaystyle\mathcal{W}_{\text{c}}=\kappa_{\perp}^{1/4}\left(1+\frac{n^{2}\overline{\rho}^{2}}{m^{2}L_{\zeta}^{2}}\right)^{1/4}\left(\frac{L_{\zeta}\overline{\rho}}{m\iota^{\prime}\overline{\psi}}\right)^{1/2},$
(23)

the perpendicular diffusion becomes unimportant. The solution is evaluated
both inside and outside the boundary layer. We conclude that for both the
inner and outer solution, the parallel and perpendicular diffusion can only be
comparable within the island separatrices when
$\mathcal{W}_{m,n}/\mathcal{W}_{\text{c}}\gtrsim 1$, where the island half-
width is,

$\displaystyle\mathcal{W}_{m,n}=\frac{2\overline{\rho}}{\overline{\psi}}\sqrt{\frac{\widetilde{\epsilon}_{m,n}(\psi)}{\iota^{\prime}(\psi)}}\bigg{\rvert}_{\psi=n/(m\iota^{\prime})}.$
(24)

In the limit that $\kappa_{\perp}\rightarrow 0$,
$\mathcal{W}_{m,n}/\mathcal{W}_{\text{c}}\gg 1$, and the parallel transport
can compete with perpendicular transport throughout the volume within the
separatrices.  In this sense, an island can be considered unimportant to the
transport if

$\displaystyle\mathcal{W}_{m,n}\ll\mathcal{W}_{\text{c}}.$ (25)

In Figure 1 we present the numerically computed values of the effective volume
(20) for a single island chain created with resonance $\epsilon_{2,1}$ with
varying amplitudes and values of $\kappa_{\perp}$. The critical value,
$\epsilon_{2,1}^{\text{crit}}=\sqrt{\kappa_{\perp}}/2$, predicted for the
parallel transport to compete with the perpendicular transport, is indicated
by the vertical dashed line. We see that, as expected, when the island is
smaller than the critical width $\mathcal{W}_{\text{c}}$, the computed
effective volume vanishes. The point at which $\mathcal{V}_{\text{PD}}$
becomes positive is within an order of magnitude of
$\epsilon_{2,1}^{\text{crit}}$. For large values of $\epsilon_{2,1}$, the
effective volume approaches the expected scaling of the total volume enclosed
by the separatrices, $\mathcal{W}_{2,1}\sim\sqrt{\epsilon_{2,1}}$, shown in
orange. At non-zero values of the perpendicular diffusion, the scaling of
$\mathcal{V}_{\text{PD}}$ is shallower than the scaling of the total island
width since the entire island chain is not dominated by parallel diffusion.

(a) $\kappa_{\perp}=10^{-6}$

(b) $\kappa_{\perp}=10^{-7}$

(c) $\kappa_{\perp}=10^{-8}$

Figure 1: The numerically computed value of the effective volume (20) is shown
(blue) for the model field (22) with the displayed values of $\epsilon_{2,1}$
and $\kappa_{\perp}$. Also displayed is the scaling of the island width with
$\epsilon_{2,1}$ (orange), and the critical value of $\epsilon_{2,1}$ for
which the parallel transport may compete with perpendicular transport (black).
The numerically computed value is predicted to approach the
$\epsilon_{2,1}^{1/2}$ scaling in the limit $\kappa_{\perp}\rightarrow 0$.

In Figure 2 we present the isotherms (red) computed with the displayed values
of $\kappa_{\perp}$ on the Poincaré section of the model magnetic field (22)
with $\epsilon_{2,1}=0.008$. As $\kappa_{\perp}$ decreases, the isotherms
conform more strongly to the structure of the magnetic field. We note the
presence of isotherms which cut across the island chain at smaller values of
$\kappa_{\perp}$ due to the topological restrictions. On the right we display
the Poincaré section with color highlighting the region over which the local
parallel diffusion is dominating. The colorscale indicates the ratio of the
local parallel to perpendicular transport. The parallel transport is strongest
near the X-points, where the misalignment of the isotherms with the local
field is the strongest. The parallel transport is relatively weak at the top
and bottom of the island chain, where the isotherms are more strongly aligned
with the magnetic field and the resonant field is weaker. As $\kappa_{\perp}$
decreases, we see the area over which the parallel transport dominates
increase as well as its magnitude in comparison with the perpendicular
transport.

(a) $\kappa_{\perp}=10^{-6}$

(b) $\kappa_{\perp}=10^{-7}$

(c) $\kappa_{\perp}=10^{-8}$

Figure 2: The Poincaré section of the model field (22) with
$\epsilon_{2,1}=0.008$ is displayed with isotherms (red) on the left computed
with the three values of $\kappa_{\perp}$ displayed. On the right, the
Poincaré section is displayed with the ratio of the parallel (18) to
perpendicular (19) diffusion. The colorscale is set to white where the
parallel diffusion is smaller than the perpendicular diffusion.

##  6 Transport across a strongly chaotic layer

In this Section we consider the relative impact of perpendicular and parallel
transport in the presence of strong chaos. We will assume a quasilinear model
(Chirikov, 1979; Hazeltine & Meiss, 2003; Lichtenberg & Lieberman, 2013) to
obtain a parallel diffusion coefficient. We will assume that the parallel
transport at a given location will be dominated by nearby resonances within an
island width and that the phase of the resonances is decorrelated after one
toroidal transit. We assume strong chaos such that the Chirikov island overlap
parameter is large,

$\displaystyle\frac{\mathcal{W}_{m,n}+\mathcal{W}_{m^{\prime},n^{\prime}}}{(\Delta\rho)_{m,n;m^{\prime},n^{\prime}}}\gg
1,$ (26)

where $\mathcal{W}_{m,n}$ is the island half-width associated with resonances
(24) and $(\Delta\rho)_{m,n;m^{\prime},n^{\prime}}$ is the distance between
the resonances.

Under the assumption that the resonances are strongly overlapping but
$\epsilon_{m,n}$ is sufficiently small that a perturbation analysis is
appropriate, the motion in the field can be treated as a random walk,

$\displaystyle\frac{\partial T}{\partial\zeta}\sim
D_{\|}\frac{\partial^{2}T}{\partial\psi^{2}},$ (27)

with a quasi-linear diffusion coefficient given by,

$\displaystyle
D_{\|}\approx\frac{\pi}{2}\sum_{m}m^{2}\left(\sum_{n}\Theta(\mathcal{W}_{m,n}/2-|\rho-\rho_{\text{r}}|)\widetilde{\epsilon}_{m,n}(\psi_{0})\right)^{2}.$
(28)

See Appendix D for details.  The local parallel transport then scales as,

$\displaystyle|\nabla_{\|}T|^{2}\sim
D_{\|}^{2}|\nabla_{\perp}^{2}T|^{2}\frac{\overline{\rho}^{4}}{\overline{\psi}^{4}L_{\zeta}^{2}}.$
(29)

We can estimate $|\nabla_{\perp}^{2}T|\sim|\nabla_{\perp}T|/L_{\perp}$ for
some perpendicular length scale $L_{\perp}$. Given the analysis in Section 5
and Appendix C, we can note several relevant length scales. Within the
boundary layer, $L_{\perp}\sim\mathcal{W}_{\text{c}}$. Outside a width
$\mathcal{W}_{m,n}$ of a given resonance and outside the boundary layer,
$L_{\perp}\sim\overline{\rho}$. Within a width $\mathcal{W}_{m,n}$ and outside
the boundary layer, $L_{\perp}\sim|\rho-\rho_{\text{r}}|$. Given that the
resonances are strongly overlapping, the perpendicular temperature gradient
will be dominated by nearby resonances. Thus we can estimate
$L_{\perp}\lesssim\mathcal{W}_{m,n}$ and,

$\displaystyle|\nabla_{\|}T|^{2}\gtrsim
D_{\|}^{2}\frac{|\nabla_{\perp}T|^{2}}{\mathcal{W}_{m,n}^{2}}\frac{\overline{\rho}^{4}}{\overline{\psi}^{4}L_{\zeta}^{2}}.$
(30)

We conclude that the parallel transport may compete with the perpendicular
transport when,

$\displaystyle\frac{D_{\|}^{2}\overline{\rho}^{4}}{\mathcal{W}_{m,n}^{2}\overline{\psi}^{4}L_{\zeta}^{2}}\sim\kappa_{\perp}.$
(31)

To test the expected scaling, we consider the model field (22) with $m=12$ and
$n=[2,\dots,10]$. We evaluate three strongly chaotic fields with amplitudes
$\epsilon_{m,n}$ chosen such that the Chirikov parameter (26) is 2,
$2\sqrt{2}$, and 4. We compute solutions to the ADE for a range of
$\kappa_{\perp}$. The effective volume (20) is shown in Figure 3 along with
the critical value of $\kappa_{\perp}$ (vertical dashed) at which the parallel
diffusion may compete (31). We note that the point at which the effective
volume becomes non-zero is within an order of magnitude of that expected from
quasilinear theory for each value of the overlap parameter. As
$\kappa_{\perp}$ decreases, the effective volume computed from the three model
fields approach each other.

Figure 3: The numerically computed value of the effective volume (20) is shown
for the model field (22) with $m=12$ resonances with amplitudes chosen to
provide strong island overlap (values of the overlap parameter (26) are
indicated). Also displayed is the critical value of $\kappa_{\perp}$ for which
the parallel transport is predicted to compete with perpendicular transport
(vertical dashed).

In Figure 4 we display the Poincaré section with isotherms (left) and with the
parallel-transport dominated regions highlighted in color (right). Although
the field is strongly chaotic, we note several isotherms near the boundaries
which adapt to the $m=12$ structure of the residual island chains. For
$\kappa_{\perp}=10^{-4}$, we note only a very small volume over which the
parallel transport dominates. As we decrease the value of $\kappa_{\perp}$, we
see that the parallel transport dominates over much of the volume contained
within the outermost island separatrices and the temperature gradient across
the chaotic layer is significantly reduced. We note an $n=6$ variation in the
structure of the parallel diffusion due to the poloidal dependence of the
radial field (22). Near the edges of the chaotic layer, the isotherms become
strongly shaped in comparison with the relatively flat isotherms in the
chaotic layer. Given that the magnitude of the radial field is smaller than
the toroidal field by the amplitude of the perturbation, this feature enhances
the parallel diffusion near the boundaries of the chaotic layer.

(a) $\kappa_{\perp}=10^{-4}$

(b) $\kappa_{\perp}=10^{-5}$

(c) $\kappa_{\perp}=10^{-6}$

Figure 4: The Poincaré section of the model field (22) with $m=12$ resonances
chosen to provide an overlap parameter (26) of 4 is displayed with isotherms
(red) on the left for the three values of $\kappa_{\perp}$ displayed. On the
right, the Poincaré section is displayed with the ratio of the parallel (18)
to perpendicular (19) diffusion. The colorscale is set to white where the
parallel diffusion is smaller than the perpendicular diffusion.

##  7 Transport across a critically chaotic layer

Although simple analytic scalings can be derived in the case of a single
island chain or strong chaos, in reality a magnetic field may possess a
complicated mixture of KAM surfaces, partial barriers, chaos, and island
chains. In this case, we can use the metric defined in §3 to assess the
relative impact of non-integrability on the transport through the volume of
parallel diffusion.

As a point of comparison, we will consider four marginally chaotic magnetic
fields generated by resonances with different poloidal mode numbers. The
amplitude of the perturbation will be chosen to achieve critical island
overlap with a comparable volume enclosed by the outermost separatrices. For
the first field, we choose $n/m=[1/4,2/4,3/4]$. Since the resonances are
equally spaced with $\Delta\psi=1/4$, we choose the perturbation amplitudes
$\epsilon_{m,n}$ such that the associated half island widths (24) are $1/8$.
We similarly consider fields with $n/m=[2/12,\dots,10/12]$,
$n/m=[3/20,\dots,17/20]$, and $n/m=[5/36,\dots,31/36]$ with critical island
overlap. The separatrices associated with the outermost island chains for each
of these model fields are located at $\rho=1/8$ and $\rho=7/8$ such that they
have a comparable “volume” of chaos.

If the remnants of the island chains are dominating the transport, we
anticipate the parallel transport to dominate when,

$\displaystyle\epsilon_{m,n}^{2}m^{2}\sim\kappa_{\perp},$ (32)

from (25), while if strong chaos is dominating the transport,

$\displaystyle\epsilon_{m,n}^{3}m^{4}\sim\kappa_{\perp},$ (33)

from (31). Since we have chosen $\epsilon_{m,n}\sim 1/m^{2}$ to achieve
critical overlap, we expect the parallel transport to be relatively more
impactful when $m$ is small. In Figure 5 we compare the effective volume
between the model fields. Since the volume contained within the outermost
separatrices differs slightly between the two model fields, we compute the
effective volume only for $\rho\in[0.25,0.75]$. As anticipated, the $m=4$
resonances produce a larger effective volume of parallel diffusion, and the
relative difference between the model fields increases with increasing
$\kappa_{\perp}$. To quantify the effect of the field structure on temperature
flattening, we also compare the total temperature differential between
$\rho=0.75$ and $\rho=0.25$. At large values of $\kappa_{\perp}$, each of the
model fields can support a significant temperature gradient, although it is
reduced for the $m=4$ field.

As can be seen from Figure 6, each of the fields are visually chaotic with
secondary island chains as well as stochastic regions. We furthermore note
that the magnitude of the parallel transport, in addition to the volume over
which it dominates, is considerably larger in the $m=4$ field. For the $m=36$
perturbation, we note that the temperature gradient is not markedly reduced in
the chaotic layer when $\kappa_{\perp}=10^{-6}$, although all of the flux
surfaces are destroyed. Thus while all of the fields would be classified as
very far from integrability by traditional measures such as the Chirikov
overlap criterion or converse KAM theory, the impact of the magnetic field
structure on the overall transport at non-zero $\kappa_{\perp}$ is remarkably
distinct.

(a)

(b)

Figure 5: (a) The effective volume (20) is computed for the model field (22)
with the displayed mode number perturbations with amplitudes $\epsilon_{m,n}$
chosen for critical island overlap. The effective volume is computed over a
subset of the entire volume, $\rho\in[0.25,0.75]$. (b) The total temperature
difference between $\rho=0.75$ and $\rho=0.25$ upon averaging over the angles.

(a) $m=4$

(b) $m=12$

(c) $m=36$

Figure 6: The ADE is solved for the model field (22) with critically
overlapping resonances of $m=4$, $m=12$, and $m=36$ with
$\kappa_{\perp}=10^{-6}$. On the left are the isotherms (red), and on the
right is the ratio of the parallel (18) to perpendicular (19) diffusion. The
colorscale is set to white where the parallel diffusion is smaller than the
perpendicular diffusion.

##  8 Conclusions

We have proposed a metric, the effective volume of parallel diffusion, which
quantifies the fractional volume over which the parallel transport is dominant
over the perpendicular transport. This enables one to determine if a magnetic
field is sufficiently integrable with regards to its impact on the transport
of heat. We validate this metric for two model fields: a single island chain
and strongly chaotic layer. We then consider a critically chaotic magnetic
field generated by resonances of different mode numbers. We find that even if
no KAM surfaces exist, the relative parallel transport can vary widely between
these critically chaotic fields, especially at larger values of
$\kappa_{\perp}$. This highlights the importance of quantifying non-
integrability not only through the stability of a given orbit or the non-
existence of a given KAM surface, but through the resulting impact on the
transport. We anticipate many applications of this analysis, such as for the
optimization of stellarator magnetic fields and analysis of transport in the
stochastic edge.

As discussed in §3, our metric appears to be related to the converse KAM
approach, which enables a quantification of the volume of phase space that
does not contain KAM surfaces of a given foliation. We expect that the
effective volume of parallel diffusion might agree with such a calculation in
the limit of small perpendicular diffusion. We reserve such a comparison for
future publication.

## Acknowledgements

The authors would like to acknowledge fruitful discussions with D. Pfefferle,
P. Constantin, R. Mackay, T. Drivas, D. Ginsberg, and A. Bhattacharjee.

## Funding

This work was supported in part by a grant from the Simons Foundation/SFARI
(560651) and by DOE Contract No. DEAC02–76CH03073.

## Declaration of interests

The authors report no conflict of interest.

## Author ORCID

E. J. Paul, https://orcid.org/0000-0002-9355-5595; S. R. Hudson,
https://orcid.org/0000-0003-1530-2733; P. Helander,
https://orcid.org/0000-0002-0460-590X.

##  Appendix A Maximum principle for ADE

We can now prove a maximum principle for the anisotropic diffusion equation
(Fraenkel, 2000). Let $T$ satisfy the ADE. We define a function $w=T+\epsilon
e^{K\textbf{x}\cdot\textbf{c}}$ for constant scalars $K>0$, $\epsilon>0$, and
constant vector c. Applying the ADE operator to $w$, we obtain,

$\displaystyle\nabla\cdot\left(\bm{\kappa}\cdot\nabla
w\right)=\epsilon\left(\textbf{c}\cdot\bm{\kappa}\cdot\textbf{c}K^{2}+\left(\nabla\cdot\bm{\kappa}\right)\cdot\textbf{c}K\right)e^{K\textbf{x}\cdot\textbf{c}}.$
(34)

From the positive definiteness of $\bm{\kappa}$, we note that
$\textbf{c}\cdot\bm{\kappa}\cdot\textbf{c}\geq\lambda_{0}|\textbf{c}|^{2}$
where $\lambda_{0}>0$ is the smallest eigenvalue of $\bm{\kappa}$. This yields
the inequality,

$\displaystyle\nabla\cdot\left(\bm{\kappa}\cdot\nabla w\right)$
$\displaystyle\geq\epsilon\left(\lambda_{0}|\textbf{c}|^{2}K^{2}-\max_{\Omega}\left|\left(\nabla\cdot\bm{\kappa}\right)\cdot\textbf{c}\right|K\right)e^{K\textbf{x}\cdot\textbf{c}}.$
(35)

We choose
$K>\max_{\Omega}\left|\left(\nabla\cdot\bm{\kappa}\right)\cdot\textbf{c}\right|/(\lambda_{0}|\textbf{c}|^{2})$
so that $\nabla\cdot\left(\bm{\kappa}\cdot\nabla w\right)>0$.

Suppose now that the maximum value of $w$ is obtained within the domain,
$\Omega$, at a point $\textbf{x}_{0}$. This implies that $\nabla
w\rvert_{\textbf{x}_{0}}=0$ and the Hessian matrix $\nabla\nabla
w\rvert_{\textbf{x}_{0}}$ is negative semi-definite. This implies that,

$\displaystyle\nabla\cdot\left(\bm{\kappa}\cdot\nabla
w\right)\rvert_{\textbf{x}_{0}}=\sum_{i,j}\kappa_{ij}(\textbf{x}_{0})\partial_{i}\partial_{j}w(\textbf{x}_{0})=\mathrm{tr}\left(\bm{\kappa}\nabla\nabla
w\right)_{\textbf{x}_{0}}\leq 0.$ (36)

The final inequality follows from expressing $\bm{\kappa}$ in terms of the
eigenvalues, $\mu_{i}$, and eigenvectors, $\overrightarrow{\textbf{x}}_{i}$,
of $\nabla\nabla w$ with $\mu_{i}\leq 0$,

$\displaystyle\bm{\kappa}=\sum_{i}\omega_{i}\overrightarrow{\textbf{x}}_{i}\overrightarrow{\textbf{x}}_{i}^{T}.$
(37)

From the positive-definiteness of $\bm{\kappa}$, we find that $\omega_{i}>0$.
We can then note that,

$\displaystyle\text{tr}\left(\bm{\kappa}\nabla\nabla
w\right)_{\textbf{x}_{0}}=\text{tr}\left(\sum_{i}\omega_{i}\mu_{i}\overrightarrow{\textbf{x}}_{i}\overrightarrow{\textbf{x}}_{i}^{T}\right)\leq
0.$ (38)

As this is in contradiction with the choice we have made for $K$, we conclude
that the maximum of $w$ must occur on $\partial\Omega$ and,

$\displaystyle
T<w\leq\max_{\overline{\Omega}}(w)=\max_{\partial\Omega}w\leq\max_{\partial\Omega}T+\epsilon\max_{\partial\Omega}e^{K\textbf{x}\cdot\textbf{c}}.$
(39)

This implies that $\max_{\overline{\Omega}}T=\max_{\partial\Omega}T$.
Otherwise, we would have $\max_{\Omega}T=\max_{\partial\Omega}T+\delta$ for
some $\delta>0$. As we are free to choose
$\epsilon<\delta/\max_{\partial\Omega}e^{K\textbf{x}\cdot\textbf{c}}$, this
would lead to a contradiction. A similar argument can be constructed to show
that the minimum value of $T$ must occur on the boundary.

##  Appendix B Topology of isotherms

In this Appendix, we discuss the possible topology of the isotherms given
properties of the ADE.

To begin, we consider the region near the boundaries. We first note that the
normal derivative of the temperature must be sign-definite on the boundaries
unless $T$ is constant in $\Omega$. In other words,
$\hat{\textbf{n}}\cdot\nabla T>0$ on $S_{-}$ and $S_{+}$ assuming
$T_{+}>T_{-}$ and $\hat{\textbf{n}}$ is oriented in the direction of
increasing $\rho$. This result is known as the Hopf lemma (Fraenkel, 2000).
Intuitively, the flux through the boundary cannot change sign since $S_{-}$
and $S_{+}$ are isotherms. As a result, in a small neighborhood of the
boundary $\partial T/\partial\rho\neq 0$, and isotherms close with the same
topology as the boundary. We conclude that at least within a small
neighborhood of $S_{-}$ and $S_{+}$, isotherms must be deformable to the
boundaries.

We now argue that the temperature gradient cannot vanish within a subvolume of
$\Omega$ if the specified temperature differential $T_{+}-T_{-}$ is non-zero.
To see this, suppose that the temperature is constant within a subvolume,
$\mathcal{V}$. We remark that under the assumption that the coefficients
appearing in the elliptic problem (9), $\bm{\kappa}$ and
$\nabla\cdot\bm{\kappa}$, are analytic, $T$ must also be analytic due to the
ellipticity of the ADE operator. If $T$ is analytic and is constant within
$\mathcal{V}$, this implies that $T$ is constant in $\Omega$. Since such a
solution is inconsistent with the boundary conditions, we conclude that the
temperature gradient can only vanish on isolated lines, surfaces, or points.
These must be saddle points given the maximum principle.

We now consider the possible linking of isothermal surfaces (Figure 7).
Suppose that there exists an isotherm, $\mathcal{S}$, with the linking of a
sphere (i.e., does not close toroidally or poloidally) or the linking of an
“island” (i.e., closes toroidally or poloidally but not both). In that case,
we can integrate the ADE over a volume $\mathcal{V}$ bounded by one of these
isotherms, $\mathcal{S}=\partial\mathcal{V}$,

$\displaystyle\int_{\mathcal{V}}d^{3}x\,\nabla\cdot\left(\bm{\kappa}\cdot\nabla
T\right)=\int_{\partial\mathcal{V}}d^{2}x\,\hat{\textbf{n}}\cdot\bm{\kappa}\cdot\nabla
T=0.$ (40)

We note that $\nabla T\propto\hat{\textbf{n}}$ and
$\hat{\textbf{n}}\cdot\nabla T$ is sign-definite on $\partial\mathcal{V}$.
Because $\bm{\kappa}$ is positive-definite,
$\hat{\textbf{n}}\cdot\bm{\kappa}\cdot\nabla T$ will also be sign-definite on
$\partial\mathcal{V}$. This leads to a contradiction unless $\nabla T=0$.
Given that the temperature gradient cannot vanish within a finite volume, we
conclude that the isotherms must link as the boundaries such that they do not
enclose a net volume.

(a)

(b)

(c)

(d)

(e)

Figure 7: Possible linking of closed isotherms (blue) which are periodic and
do not intersect the boundaries (red): (a) “spherical” (not closed toroidally
or poloidally), (b) poloidally closed “island” structure, (c) toroidally
closed “island” structure, (d) toroidally and poloidally closed (deformable to
boundary), and (e) a toroidally and poloidally closed surface with the
addition of a handle.

Now that we have argued that the isotherms must link toroidally and poloidally
as the boundaries, we now discuss their topology. We consider an isotherm
which links as the boundaries but with the addition of a small handle (Figure
7(e)). Such an isotherm is consistent with the integral condition (40) as it
does not enclose a volume. We conclude that the isotherms must link toroidally
and poloidally, with the possible addition of a handle. However, isotherms
cannot enclose a net volume, as an “island” or “sphere” structure.

##  Appendix C Single island chain analysis

In this Appendix we analyze the ADE for the case of a single small island
chain, as discussed in Section 5. We perform a perturbation analysis in the
small parameter $\epsilon=|\chi_{1}|/|\chi_{0}|\ll 1$ where
$\chi_{0}=\iota^{\prime}\psi^{2}/2$ and $\chi_{1}=\chi-\chi_{0}$. In the limit
$\chi=\chi_{0}$, the solution of the ADE is given by,

$\displaystyle T_{0}(\rho)=T_{-}+\frac{\rho(T_{+}-T_{-})}{\overline{\rho}}.$
(41)

We will call this the unperturbed temperature. In computing the first order
correction, the perturbed temperature, we apply the ansatz
$T_{1}(\rho,\theta,\zeta)=\hat{T}(\rho)\cos(m\theta-n\zeta)$. The first order
perturbation to the ADE reads,

$\displaystyle\kappa_{\perp}\nabla^{2}T_{1}+(1-\kappa_{\perp})\hat{\textbf{b}}_{0}\cdot\nabla\left(\hat{\textbf{b}}_{0}\cdot\nabla
T_{1}+\hat{\textbf{b}}_{1}\cdot\nabla T_{0}\right)=0,$ (42)

where $\hat{\textbf{b}}_{0}$ and $\hat{\textbf{b}}_{1}$ are the unit vectors
in the direction of the magnetic field at $\mathcal{O}(\epsilon^{0})$ and
$\mathcal{O}(\epsilon^{1})$. Applying the ansatz for $T_{1}$ results in the
ODE,

$\kappa_{\perp}\left(\hat{T}^{\prime\prime}(\rho)-\hat{T}(\rho)\left(\frac{m^{2}}{\overline{\rho}^{2}}+\frac{n^{2}}{L_{\zeta}^{2}}\right)\right)\\\
-(1-\kappa_{\perp})\frac{m^{2}\iota^{\prime
2}\overline{\psi}^{2}}{L_{\zeta}^{2}\overline{\rho}^{2}\left(1+\frac{n^{2}\overline{\rho}^{2}}{m^{2}L_{\zeta}^{2}}\right)}\left(\hat{T}(\rho)\left(\rho-\rho_{\text{r}}\right)^{2}-\frac{\mathcal{W}_{m,n}^{2}(\rho-\rho_{\text{r}})}{4\overline{\rho}}(T_{+}-T_{-})\right)=0,$
(43)

where $\mathcal{W}_{m,n}$ is the island half-width (24). In obtaining the
above expression, we focus on a region near the $n/m$ rational surface such
that the unperturbed field strength can be approximated as
$|\textbf{B}_{0}|=\left(\overline{\psi}/\overline{\rho}^{2}\right)\left(1+\left(\iota^{\prime}\overline{\psi}\rho/L_{\zeta}\right)^{2}\right)^{1/2}\approx\left(\overline{\psi}/\overline{\rho}^{2}\right)\left(1+\left(n\overline{\rho}/(mL_{\zeta})\right)^{2}\right)^{1/2}$.
Due to the smallness of $\kappa_{\perp}$, outside of a small boundary layer of
width $\mathcal{W}_{\text{c}}$, the perpendicular diffusion becomes
unimportant. We can express the first order perturbation to the ADE in terms
of the layer variable $X=(\rho-\rho_{\text{r}})/\mathcal{W}_{\text{c}}$, where
$\rho_{\text{r}}=\overline{\rho}n/(m\iota^{\prime}\overline{\psi})$ is the
location of the resonant surface,

$\kappa_{\perp}\left(\frac{\hat{T}^{\prime\prime}(X)}{\mathcal{W}_{\text{c}}^{2}}-\hat{T}(X)\left(\frac{m^{2}}{\overline{\rho}^{2}}+\frac{n^{2}}{L_{\zeta}^{2}}\right)\right)\\\
-(1-\kappa_{\perp})\frac{m^{2}\iota^{\prime
2}\overline{\psi}^{2}}{L_{\zeta}^{2}\overline{\rho}^{2}\left(1+\frac{n^{2}\overline{\rho}^{2}}{m^{2}L_{\zeta}^{2}}\right)}\left(\hat{T}(X)\mathcal{W}_{\text{c}}^{2}X^{2}-\frac{\mathcal{W}_{m,n}^{2}\mathcal{W}_{\text{c}}X}{4\overline{\rho}}(T_{+}-T_{-})\right)=0.$
(44)

We will analyze the solution both inside and outside of the layer, similar to
the analysis in Fitzpatrick (1995).

###  C.1 Inner solution

Inside the boundary layer, the second term in (44) is unimportant as
$\mathcal{W}_{\text{c}}/\overline{\rho}\ll 1$ by assumption. Thus the first
and third terms must balance each, allowing us to obtain the scaling for the
boundary layer width as (23). Balance between the first and last terms yields
a scaling for the temperature perturbation as,

$\displaystyle\hat{T}(X)\sim\frac{\mathcal{W}_{m,n}^{2}}{\overline{\rho}\mathcal{W}_{\text{c}}}(T_{+}-T_{-}).$
(45)

This matches the scaling in (32) of Fitzpatrick (1995). This implies that
$|\nabla_{\perp}T_{1}|\sim\left(\mathcal{W}_{m,n}/\mathcal{W}_{\text{c}}\right)^{2}|\nabla_{\perp}T_{0}|$.
When $\mathcal{W}_{m,n}\ll\mathcal{W}_{\text{c}}$,
$|\nabla_{\perp}T_{0}|\gg|\nabla_{\perp}T_{1}|$, and the ratio of parallel to
perpendicular diffusion scales as,

$\displaystyle\frac{|\nabla_{\|}T|^{2}}{\kappa_{\perp}|\nabla_{\perp}T|^{2}}\sim\frac{\mathcal{W}_{m,n}^{4}}{\mathcal{W}_{\text{c}}^{4}}\ll
1,$ (46)

since
$|\nabla_{\|}T_{1}|\sim\mathcal{W}_{\text{c}}^{-1}\kappa_{\perp}^{1/2}T_{1}$.
Otherwise, if $\mathcal{W}_{m,n}/\mathcal{W}_{\text{c}}\gtrsim 1$, the ratio
can be $\sim 1$.

###  C.2 Outer solution

Outside of the boundary layer, the perpendicular diffusion can be ignored to a
good approximation, and the ADE reduces to,

$\displaystyle\nabla\cdot\left(\hat{\textbf{b}}\hat{\textbf{b}}\cdot\nabla
T\right)=0.$ (47)

Outside of the island separatrix, we apply a perturbation analysis to obtain,

$\displaystyle\hat{\textbf{b}}_{0}\cdot\nabla\left(\hat{\textbf{b}}_{0}\cdot\nabla
T_{1}+\hat{\textbf{b}}_{1}\cdot\nabla T_{0}\right)=0.$ (48)

This yields the solution,

$\displaystyle\hat{T}^{(0)}(\rho)=\mathcal{W}_{m,n}^{2}\frac{T_{+}-T_{-}}{4\overline{\rho}(\rho-\rho_{\text{r}})},$
(49)

where the superscript $(0)$ indicates the solution which neglects
perpendicular diffusion. This matches the scaling of (38) in Fitzpatrick
(1995). Given (48), the parallel temperature gradient is negligible in this
region if perpendicular diffusion is ignored. Now retaining the perpendicular
diffusion terms from (44), the correction to $\hat{T}(\rho)$ scales as
$\hat{T}^{(1)}(\rho)\sim\hat{T}(\rho)\mathcal{W}_{\text{c}}^{4}/\left(\overline{\rho}^{2}(\rho-\rho_{\text{r}})^{2}\right)$,
and the parallel diffusion will scale as
$|\nabla_{\|}T|\sim\kappa_{\perp}^{1/2}\mathcal{W}_{\text{c}}^{-2}(\rho-\rho_{\text{r}})\hat{T}^{(1)}(\rho)$.
Given that $|\nabla_{\perp}T_{0}|\gg|\nabla_{\perp}T_{1}|$ in this region, the
ratio of the parallel to perpendicular diffusion will have the scaling,

$\displaystyle\frac{|\nabla_{\|}T|^{2}}{\kappa_{\perp}|\nabla_{\perp}T|^{2}}\sim\frac{\mathcal{W}_{m,n}^{4}\mathcal{W}_{\text{c}}^{4}}{\overline{\rho}^{4}(\rho-\rho_{\text{r}})^{4}}\ll
1.$ (50)

Within the island separatrix, (47) implies that the temperature must be a
function of the island flux surfaces. This is not possible unless the
temperature is flattened across the island chain given the maximum principle.
At non-zero $\kappa_{\perp}$, a finite temperature gradient is allowed. Given
that the temperature is flattened within the separatrices if perpendicular
diffusion is not included, both the parallel and perpendicular temperature
gradients will be driven by $\hat{T}^{(1)}(\rho)$ with scaling given by
$|\nabla_{\|}T|\sim\kappa_{\perp}^{1/2}\mathcal{W}_{\text{c}}^{-2}(\rho-\rho_{\text{r}})\hat{T}^{(1)}(\rho)$
and $|\nabla_{\perp}T|\sim\hat{T}^{(1)}(\rho)/(\rho-\rho_{\text{r}})$. This
yields the scaling for the ratio of the parallel to perpendicular diffusion
as,

$\displaystyle\frac{|\nabla_{\|}T|^{2}}{\kappa_{\perp}|\nabla_{\perp}T|^{2}}\sim\frac{(\rho-\rho_{\text{r}})^{4}}{\mathcal{W}_{\text{c}}^{4}}\gg
1.$ (51)

We conclude that for the outer solution, the ratio
$|\nabla_{\|}T|^{2}/\left(\kappa_{\perp}|\nabla_{\perp}T|^{2}\right)$ can only
be $\gtrsim 1$ if $|\rho-\rho_{\text{r}}|<\mathcal{W}_{m,n}$ and
$\mathcal{W}_{\text{c}}\ll\mathcal{W}_{m,n}$.

##  Appendix D Strong chaos analysis

In this Appendix we analyze the ADE for the case of strong chaos, as discussed
in Section 6.  If we assume that $\epsilon_{m,n}$ is sufficiently small, we
can linearize the trajectories with respect to the small parameter
$\epsilon=|\chi_{1}|/|\chi_{0}|\ll 1$ where
$\chi_{0}=\iota^{\prime}\psi^{2}/2$ and $\chi_{1}=\chi-\chi_{0}$,

$\displaystyle\psi$
$\displaystyle=\psi_{0}+\psi_{1}(\zeta)+\mathcal{O}(\epsilon^{2})$ (52a)
$\displaystyle\theta$
$\displaystyle=\theta_{0}+\iota(\psi_{0})\zeta+\theta_{1}(\zeta)+\mathcal{O}(\epsilon^{2}),$
(52b)

where $\iota(\psi_{0})=\iota^{\prime}\psi_{0}$ is the rotational transform
associated with the unperturbed Hamiltonian (22). The linearized trajectories
satisfy,

$\displaystyle\frac{d\psi_{1}}{d\zeta}$
$\displaystyle=\sum_{m,n}m\widetilde{\epsilon}_{m,n}(\psi_{0})\sin(m(\theta_{0}+\iota(\psi_{0})\zeta)-n\zeta)$
(53a) $\displaystyle\frac{d\theta_{1}}{d\zeta}$
$\displaystyle=\iota^{\prime}\psi_{1}+\sum_{m,n}\widetilde{\epsilon}_{m,n}^{\prime}(\psi_{0})\cos(m(\theta_{0}+\iota(\psi_{0})\zeta)-n\zeta),$
(53b)

and are computed to be,

$\displaystyle\psi_{1}(\zeta)$
$\displaystyle=-\sum_{m,n}m\widetilde{\epsilon}_{m,n}(\psi_{0})\frac{\cos(m(\theta_{0}+\iota(\psi_{0})\zeta)-n\zeta)}{m\iota(\psi_{0})-n}$
(54a) $\displaystyle\theta_{1}(\zeta)$
$\displaystyle=\sum_{m,n}\left[\frac{-m\iota^{\prime}\widetilde{\epsilon}_{m,n}(\psi_{0})}{(m\iota(\psi_{0})-n)^{2}}+\frac{\widetilde{\epsilon}_{m,n}^{\prime}(\psi_{0})}{m\iota(\psi_{0})-n}\right]\sin(m(\theta_{0}+\iota(\psi_{0})\zeta)-n\zeta).$
(54b)

Thus we can see that the resonant terms in $\theta_{1}(\zeta)$ within a width
$\mathcal{W}_{m,n}/2$ (24) of a given rational surface will be
$\mathcal{O}(1/m)$ and will thus dominate the evolution. We can then
approximate the trajectories as,

$\displaystyle\frac{d\psi_{1}}{d\zeta}$
$\displaystyle\approx\sum_{m,n}\Theta(\mathcal{W}_{m,n}/2-|\rho-\rho_{\text{r}}|)m\widetilde{\epsilon}_{m,n}(\psi_{0})\sin(m\theta_{0})$
(55a) $\displaystyle\frac{d\theta_{1}}{d\zeta}$
$\displaystyle\approx\iota^{\prime}\psi_{1}+\sum_{m,n}\Theta(\mathcal{W}_{m,n}/2-|\rho-\rho_{\text{r}}|)\widetilde{\epsilon}_{m,n}^{\prime}(\psi_{0})\cos(m\theta_{0}),$
(55b)

where $\Theta$ is the Heaviside step function. In one toroidal period, the
flux changes by an amount given by  $\Delta\psi\approx 2\pi d\psi_{1}/d\zeta$
and the poloidal angle changes by an amount given by $\Delta\theta\approx 2\pi
d\theta_{1}/d\zeta$. We note that the uncertainty in the poloidal angle will
be $\mathcal{O}(1)$ if the resonances are strongly overlapping. In this way,
with every toroidal period, the poloidal angle will be kicked by an amplitude
that depends on the surrounding resonances. We assume that $\theta_{0}$ can,
therefore, be treated as a random variable, and we approximate the motion in
such a field as a random walk,

$\displaystyle\frac{\partial T}{\partial\zeta}\sim
D_{\|}\frac{\partial^{2}T}{\partial\psi^{2}},$ (56)

with a diffusion coefficient given by,

$\displaystyle D_{\|}\sim\frac{\langle(\Delta\psi)^{2}\rangle}{2\Delta\zeta},$
(57)

where $\langle...\rangle$ indicates an average over the possible initial
conditions, $\theta_{0}$.  The full expression for $D_{\|}$ is given in (28).

## References

  * Amestoy et al. (2019) Amestoy, P. R., Buttari, A., L’Excellent, J.-Y. & Mary, T. 2019 Performance and scalability of the block low-rank multifrontal factorization on multicore architectures. ACM Transactions on Mathematical Software 45, 2:1–2:26.
  * Amestoy et al. (2001) Amestoy, P. R., Duff, I. S., Koster, J. & L’Excellent, J.-Y. 2001 A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis and Applications 23 (1), 15–41.
  * Balay et al. (2021a) Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D., Knepley, M. G., May, D. A., McInnes, L. C., Mills, R. Tran, Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S. & Zhang, H. 2021a PETSc users manual. Tech. Rep. ANL-95/11 - Revision 3.15. Argonne National Laboratory.
  * Balay et al. (2021b) Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D., Knepley, M. G., May, D. A., McInnes, L. Curfman, Mills, R. Tran, Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S. & Zhang, H. 2021b PETSc Web page.
  * Balay et al. (1997) Balay, S., Gropp, W. D., McInnes, L. Curfman & Smith, B. F. 1997 Efficient management of parallelism in object oriented numerical software libraries. In Modern Software Tools in Scientific Computing (ed. E. Arge, A. M. Bruaset & H. P. Langtangen), pp. 163–202. Birkhäuser Press.
  * Bardóczi et al. (2016) Bardóczi, L., Rhodes, T. L., Carter, T. A., Crocker, N. A., Peebles, W. A. & Grierson, B. A. 2016 Non-perturbative measurement of cross-field thermal diffusivity reduction at the O-point of 2/1 neoclassical tearing mode islands in the DIII-D tokamak. Physics of Plasmas 23 (5), 052507.
  * Biewer et al. (2003) Biewer, T. M., Forest, C. B., Anderson, J. K., Fiksel, G., Hudson, B., Prager, S. C., Sarff, J. S., Wright, J. C., Brower, D. L., Ding, W. X. & others 2003 Electron heat transport measured in a stochastic magnetic field. Physical Review Letters 91 (4), 045004.
  * Cary & Hanson (1986) Cary, J. R. & Hanson, J. D. 1986 Stochasticity reduction. The Physics of Fluids 29 (8), 2464–2473.
  * Cary & Hanson (1991) Cary, J. R. & Hanson, J. D. 1991 Simple method for calculating island widths. Physics of Fluids B: Plasma Physics 3 (4), 1006–1014.
  * Cary & Littlejohn (1983) Cary, J. R. & Littlejohn, R. G. 1983 Noncanonical Hamiltonian mechanics and its application to magnetic field line flow. Annals of Physics 151 (1), 1–34.
  * del Castillo-Negrete & Chacon (2012) del Castillo-Negrete, D. & Chacon, L. 2012 Parallel heat transport in integrable and chaotic magnetic fields. Physics of Plasmas 19 (5), 056112.
  * Chirikov (1979) Chirikov, B. V. 1979 A universal instability of many-dimensional oscillator systems. Physics Reports 52 (5), 263–379.
  * Chuvilgin & Ptuskin (1993) Chuvilgin, L. G. & Ptuskin, V. S. 1993 Anomalous diffusion of cosmic rays across the magnetic field. Astronomy and Astrophysics 279, 278–297.
  * Evans et al. (2005) Evans, T. E., Moyer, R. A., Watkins, J. G., Thomas, P. R., Osborne, T. H., Boedo, J. A., Fenstermacher, M. E., Finken, K. H., Groebner, R. J., Groth, M. & others 2005 Suppression of large edge localized modes in high confinement DIII-D plasmas with a stochastic magnetic boundary. Journal of Nuclear Materials 337, 691–696.
  * Falcolini & de La Llave (1992) Falcolini, C. & de La Llave, R. 1992 A rigorous partial justification of Greene’s criterion. Journal of Statistical Physics 67 (3), 609.
  * Feng et al. (2011) Feng, Y., Kobayashi, M., Lunt, T. & Reiter, D. 2011 Comparison between stellarator and tokamak divertor transport. Plasma Physics and Controlled Fusion 53 (2), 024009\.
  * Fitzpatrick (1995) Fitzpatrick, R. 1995 Helical temperature perturbations associated with tearing modes in tokamak plasmas. Physics of Plasmas 2 (3), 825–838.
  * Fraenkel (2000) Fraenkel, L. E. 2000 An Introduction to Maximum Principles and Symmetry in Elliptic Problems. Cambridge University Press.
  * Greene (1979) Greene, J. M. 1979 A method for determining a stochastic transition. Journal of Mathematical Physics 20 (6), 1183–1201.
  * Hameiri & Bhattacharjee (1987) Hameiri, E. & Bhattacharjee, A. 1987 Entropy production and plasma relaxation. Physical Review A 35 (2), 768\.
  * Hanson & Cary (1984) Hanson, J. D. & Cary, J. R. 1984 Elimination of stochasticity in stellarators. The Physics of Fluids 27 (4), 767–769.
  * Hazeltine & Meiss (2003) Hazeltine, R. D. & Meiss, J. D. 2003 Plasma Confinement. Courier Corporation.
  * Helander et al. (2021) Helander, P., Hudson, S. R. & Paul, E. J. 2021 On heat conduction in an irregular magneticfield: Part I. submitted to Journal of Plasma Physics.
  * Hölzl et al. (2009) Hölzl, M., Günter, S., Classen, I. G. J., Yu, Q., Delabie, E., Team, TEXTOR & others 2009 Determination of the heat diffusion anisotropy by comparing measured and simulated electron temperature profiles across magnetic islands. Nuclear Fusion 49 (11), 115009.
  * Hudson & Breslau (2008) Hudson, S. R. & Breslau, J. 2008 Temperature contours and ghost surfaces for chaotic magnetic fields. Physical Review Letters 100 (9), 095001.
  * Jokipii & Parker (1969) Jokipii, J. R. & Parker, E. N. 1969 Cosmic-ray life and the stochastic nature of the galactic magnetic field. The Astrophysical Journal 155, 799.
  * Lichtenberg & Lieberman (2013) Lichtenberg, A. J. & Lieberman, M. A. 2013 Regular and Stochastic Motion, Applied Mathematical Sciences, vol. 38. Springer Science & Business Media.
  * Lorenzini et al. (2009) Lorenzini, R., Martines, E., Piovesan, P., Terranova, D., Zanca, P., Zuin, M., Alfier, A., Bonfiglio, D., Bonomo, F., Canton, A. & others 2009 Self-organized helical equilibria as a new paradigm for ohmically heated fusion plasmas. Nature Physics 5 (8), 570–574.
  * MacKay (1992) MacKay, R. S. 1992 Greene’s residue criterion. Nonlinearity 5 (1), 161.
  * MacKay (2018) MacKay, R. S. 2018 Finding the complement of the invariant manifolds transverse to a given foliation for a 3D flow. Regular and Chaotic Dynamics 23 (6), 797–802.
  * MacKay et al. (1984) MacKay, R. S., Meiss, J. D. & Percival, I. C. 1984 Transport in Hamiltonian systems. Physica D: Nonlinear Phenomena 13 (1-2), 55–81.
  * MacKay & Percival (1985) MacKay, R. S. & Percival, I. C. 1985 Converse KAM: theory and practice. Communications in Mathematical Physics 98 (4), 469–512.
  * Mather (1984) Mather, J. N. 1984 Non-existence of invariant circles. Ergodic theory and dynamical systems 4 (2), 301–309.
  * Mather (1986) Mather, J. N. 1986 A criterion for the non-existence of invariant circles. Publications Mathématiques de l’IHÉS 63, 153–204.
  * Meiss (2015) Meiss, J. D. 2015 Thirty years of turnstiles and transport. Chaos: An Interdisciplinary Journal of Nonlinear Science 25 (9), 097602.
  * Meiss & Ott (1986) Meiss, J. D. & Ott, E. 1986 Markov tree model of transport in area-preserving maps. Physica D: Nonlinear Phenomena 20 (2-3), 387–402.
  * Meskat et al. (2001) Meskat, J. P., Zohm, H., Gantenbein, G., Günter, S., Maraschek, M., Suttrop, W., Yu, Q. & ASDEX Upgrade Team 2001 Analysis of the structure of neoclassical tearing modes in ASDEX Upgrade. Plasma Physics and Controlled Fusion 43 (10), 1325.
  * Punjabi et al. (2008) Punjabi, A., Ali, H., Evans, T. & Boozer, A. 2008 Modeling of stochastic broadening in a poloidally diverted discharge with piecewise analytic symplectic mapping flux functions. Physics of Plasmas 15 (8), 082507.
  * Sakakibara et al. (2008) Sakakibara, S., Watanabe, K. Y., Suzuki, Y., Narushima, Y., Ohdachi, S., Nakajima, N., Watanabe, F., Garcia, L., Weller, A., Toi, K. & others 2008 MHD study of the reactor-relevant high-beta regime in the Large Helical Device. Plasma Physics and Controlled Fusion 50 (12), 124014.
  * Sarff et al. (2013) Sarff, J. S., Almagri, A. F., Anderson, J. K., Borchardt, M., Carmody, D., Caspary, K., Chapman, B. E., Den Hartog, D. J., Duff, J., Eilerman, S. & others 2013 Overview of results from the MST reversed field pinch experiment. Nuclear Fusion 53 (10), 104017.


††thanks: Present affiliation: Department of Physics, Payame Noor University,
P.O. Box 19395-3697 Tehran, Iran††thanks<EMAIL_ADDRESS>
# Holes’ character and bond versus charge disproportionation in $s\textendash
p$ $ABX_{3}$ perovskites

Mohammad Reza Benam  Department of Physics $\&$ Astronomy, University of
British Columbia, Vancouver, British Columbia, Canada V6T 1Z1  Stewart Blusson
Quantum Matter Institute, University of British Columbia, Vancouver, British
Columbia, Canada V6T 1Z4  Kateryna Foyevtsova  Arash Khazraie  Ilya Elfimov
George A. Sawatzky  Department of Physics $\&$ Astronomy, University of
British Columbia, Vancouver, British Columbia, Canada V6T 1Z1  Stewart Blusson
Quantum Matter Institute, University of British Columbia, Vancouver, British
Columbia, Canada V6T 1Z4

###### Abstract

We use density functional theory methods to study the electronic structures of
a series of $s\textendash p$ cubic perovskites $ABX_{3}$: the experimentally
available SrBiO3, BaBiO3, BaSbO3, CsTlF3, and CsTlCl3, as well as the
hypothetical MgPO3, CaAsO3, SrSbO3, and RaMcO3. We use tight-binding modeling
to calculate the interatomic hopping integrals $t_{sp\sigma}$ between the
$B\textendash s$ and $X\textendash p$ atomic orbitals and charge-transfer
energies $\Delta$, which are the two most important parameters that determine
the low-energy electron and hole states of these systems. Our calculations
elucidate several trends in $t_{sp\sigma}$ and $\Delta$ as one moves across
the periodic table, such as the relativistic energy lowering of the
$B\textendash s$ orbital in heavy $B$ cations leading to strongly negative
$\Delta$ values. Our results are discussed in connection with the general
phase diagram for $s\textendash p$ cubic perovskites proposed in Ref.
Khazraie2, , where the parent superconductors SrBiO3 and BaBiO3 are found to
be in the regime of negative $\Delta$ and large $t_{sp\sigma}$. Here, we
explore this further and search for new materials with similar parameters,
which could lead to the discovery of new superconductors. Also, some
considerations are offered regarding a possible relation between the physical
properties of a given $s\textendash p$ compound (such as its tendency to bond
disproportionate and the maximal achievable superconducting transition
temperature) and its electronic structure.

##  I Introduction

Materials with a cubic perovskite structure $ABX_{3}$, where the anion $X$ can
be an oxygen or a halogen and the possible cations $A$ and $B$ include a broad
variety of elements or even molecules, have attracted considerable attention
due to their rich physics. Indeed, among their intriguing properties are
metal-insulator transitionsMedarde ; Imada , high transition temperature
($T_{\text{c}}$) superconductivityBendorz ; Sleight ; Cava , ferroelectricity,
ferromagnetism, applicability in photovoltaicsGreen , colossal
magnetoresistance Ramirez , magnetoelectricity Kimura , and a topological
insulating stateJin ; Jin2 ; Yan . The crystal structure of the $ABX_{3}$
cubic perovskites consists of a three-dimensional network of corner-sharing
$BX_{6}$ octahedra intercalated with $A$ cations at the twelve-fold anion-
coordinated sites. One of the well-known and widely studied $ABX_{3}$
compounds is $A$BiO3, with $A$ = Ba or Sr. Upon hole doping, achieved via
chemical substitutions, these systems become superconducting with a
surprisingly high maximal $T_{\text{c}}$ of 30 KSleight ; Cava ; Mattheiss ;
Kazakov . As stoichiometry is approached in the pure parent compound, however,
the superconductivity gives way to an insulating state featuring a so-called
breathing structural distortion, where the BiO6 octahedra disproportionate
into small and large ones in a rock-salt-like patternCox ; Cox2 ; Zhou ;
Sleight2 .

Although in the early years following the discovery of $A$BiO3, their
breathing distortion was viewed as a result of charge disproportionation of
the nominally tetravalent Bi4+ ions into Bi3+ and Bi5+Varma ; Taraphder ; Hase
, recent theoreticalMattheiss2 ; Harrison ; Foyevtsova ; Khazraie ; Dalpian ;
Khazraie2  as well as experimentalHair ; Orchard ; Wertheim ; Salem ; Plumb ;
Balandeh  studies have seriously challenged this idea. In a more realistic
microscopic picture, developed by some of us in Refs. Foyevtsova, ; Khazraie,
; Khazraie2, , one starts by recognizing the negative charge-transfer nature
of the $A$BiO3 electronic states, i. e., that the O$\textendash 2p$ states are
in fact higher in energy than the semi-core Bi$\textendash 6s$ states (by
amount $\Delta$), as depicted in the top panel of Fig. 1 (a). However, the
most important parameter of all shaping the $A$BiO3 electronic structure is
the strong hybridization between the Bi$\textendash 6s$ atomic orbital and the
$a_{1g}$ molecular orbital (MO) formed by the O$\textendash 2p_{\sigma}$
orbitals of the oxygen octahedral cage [see Fig. 1 (b)]. It produces a huge
splitting between the bonding and anti-bonding bands, much larger than the
charge transfer energy, with the latter band landing at the Fermi energy above
the O non-bonding states, as depicted in the middle panel of Fig. 1 (a). Since
the character of the conductance anti-bonding band is predominantly that of
the O$\textendash a_{1g}$ MO, the average Bi oxidation state approaches 3+,
leaving two self-doped ligand holes, $\underline{L}$, per oxygen octahedron as
$2{\text{Bi}}^{4+}\rightarrow 2{\text{Bi}}^{3+}\underline{L}^{2}$, in what
Alex Zunger and co-workers called a “self-regulating response”Dalpian ;
Raebiger . Upon the breathing distortion, resulting from the strong electron -
breathing phonon interaction, which increases the short-bond length
Bi$\textendash 6s-$O$\textendash 2p$ hopping integrals and therefore
stabilizes further the bonding state, the ligand holes condense pairwise onto
the small octahedra as
$2{\text{Bi}}^{3+}\underline{L}^{2}\rightarrow[{\text{Bi}}^{3+}]_{\text{large}}+[{\text{Bi}}^{3+}\underline{L}^{2}]_{\text{small}}$,
resulting in nearly the same valence states for the two inequivalent Bi
ionsFoyevtsova , a situation to be called bond, rather than charge,
disproportionation. As shown in the bottom panel of Fig. 1 (a), this process
is associated with opening of a charge gap at the center of the anti-bonding
band of the cubic structure.

Figure 1: (a) A schematic diagram of the Bi$\textendash 6s$ and O$\textendash
2p$ energy levels in $A$BiO3 before (top panel) and after (middle panel)
hybridization. ”A“ and ”B“ denote an anti-bonding and a bonding band,
respectively. The bottom panel demonstrates the effect of Bi-O bond
disproportionation, whereby the bonding and the anti-bonding bands are each
split into two subbands associated with large and small BiO6 octahedra,
denoted as ”B(L)“, ”B(S)“, ”A(L)“, and ”A(S)“, respectively, and a charge gap
is opened as a result. (b) The six O$\textendash 2p_{\sigma}$ orbitals of an
O6 octahedron and their molecular orbital combinations, $a_{1g}$, $t_{1u}$,
and $e_{g}$, with their corresponding energies in units of nearest-neighbor
$pp\sigma$ hopping integral $-t_{pp}$.

With this picture of the $A$BiO3 electronic structure in mind, a further step
was taken in Ref. Khazraie2,  and a general phase diagram was proposed to
describe a crossover from a bond- to a charge-disproportionated regime in
$s\textendash p$ cubic perovskites similar to $A$BiO3. It was shown that there
are two main electronic parameters that determine the regime a given system
will end up in: the charge-transfer energy $\Delta$ and the hybridization
between the $B$-cation $s$ orbital and the oxygen $a_{1g}$ MO characterized by
the hopping integral $T_{sp\sigma}$. However, even though a number of
$s\textendash p$ cubic perovskites other than $A$BiO3 are known, such as the
recently synthesized CsTlCl3 and CsTlF3Retuerto , and some of them even
superconduct (BaPb1-xSbxO3Singh ; Cava2 ), no real examples were discussed in
Ref. Khazraie2,  in relation with the proposed phase diagram. In order to fill
this gap, in the present paper we use ab initio theoretical methods to study
the electronic structures of the above mentioned existing $s\textendash p$
cubic perovskites and also of a systematic series of hypothetical $AB$O3
systems, with $A$ and $B$ cations being the group IIa and group Va elements,
respectively. We hope that this study, conducted in the light of the notion of
charge versus bond disproportionation, will add to our understanding of
superconductivity in the $s\textendash p$ cubic perovskites and also guide the
discovery of new superconductors.

##  II Method

Our electronic structure calculations are performed within density functional
theory (DFT)Kohn  using the full-potential linearized augmented plane-wave
method as implemented in the WIEN2k packageBlaha . We employ the generalized
gradient approximation (GGA)Perdew  for the exchange-correlation potential.
For all our $s\textendash p$ systems, a simplified cubic $Pm\bar{3}m$ crystal
structure is assumed, with both the tilting and breathing distortions
neglected, and the volume is fully relaxed within GGA. The basis set size is
fixed by setting $R_{\text{MT}}K_{\text{max}}=7$, where $R_{\text{MT}}$ is the
smallest muffin-tin sphere radius and $K_{\text{max}}$ is the cut-off wave
vector. A $12\times 12\times 12$ grid of $\boldsymbol{k}$-points is used for
integrating over the first Brillouin zone. Atomic and molecular orbital
projections are done within muffin-tin spheres. Projections onto molecular
orbitals are done with a modified version of WIEN2k, as discussed in Ref.
Foyevtsova2, . Tight-binding (TB) parameters are obtained by using the
maximally localized Wannier functions (MLWF) method as implemented in the
wannier90 codeMostofi .

Figure 2: The electronic band-structures of (a) MgPO3, (b) CaAsO3, (c) SrSbO3, (d) BaBiO3, and (e) RaMcO3 calculated in GGA and plotted with black solid lines. The Fermi energy is set to zero and marked by a horizontal black dashed line. The red and yellow circles indicate the presence of, respectively, the O$\textendash a_{1g}$ and cation $B\textendash s$ orbital characters in a given Bloch eigenstate, with the amount of their contribution being proportional to the circles’ radii. The dashed blue lines represent the eigenstates of the MLWF-based TB model. Panels (f) and (g) use a two-dimensional (2D) analogue of the cubic perovskite structure to explain the absence of hybridization between the O$\textendash 2p_{\sigma}$ and $B\textendash s$ Bloch functions at the $\Gamma$ point [$\boldsymbol{k}=(0,0,0)$] and its maximum strength at the $R$ point [$\boldsymbol{k}=(\pi,\pi,\pi)$], respectively. The dashed lines mark boundaries between 2D unit cells with the lattice constant $a$. (h) The amount of the O$\textendash a_{1g}$ molecular orbital character in the anti-bonding band of the $AB$O3 TB models, $|\langle a_{1g}(\boldsymbol{k})|\psi^{\text{TB}}(\boldsymbol{k})\rangle|^{2}$, at the $R$ point as a function of chemical composition. Figure 3: (a), (b) The electronic band-structures of CsTlF3 and CsTlCl3. Notations are similar to those in Fig. 2. (c), (d) Partial densities of states of CsTlF3 and CsTlCl3, projected onto the Tl$\textendash 6s$ atomic orbital and the molecular combinations of the F$\textendash 2p_{\sigma}$ or Cl$\textendash 3p_{\sigma}$ orbitals. The zero of energy is at the Fermi energy.  Figure 4: The partial densities of state projected onto the $B\textendash s$ atomic orbital and the molecular combinations of O$\textendash 2p_{\sigma}$ orbitals of (a) MgPO3, (b) CaAsO3, (c) SrSbO3, (d) BaBiO3 and (e) RaMcO3. The Fermi energy is set to zero and marked by a vertical black solid line. Table 1: GGA equilibrium lattice constants $a$, $B\textendash s-X\textendash p$ hopping integrals $t_{sp\sigma}$, and charge-transfer energies $\Delta$ of the studied $ABX_{3}$ cubic perovskites. | Has been synthesized? |  $a$ (Å) |  $t_{sp\sigma}$ (eV) |  $\Delta$ (eV)
---|---|---|---|---
MgPO3 | No | 3.667 | 2.65 | 1.83
CaAsO3 | No | 3.919 | 2.30 | -1.44
SrSbO3 | No | 4.233 | 2.08 | -1.08
BaSbO3 | YesCava2  | 4.280 | 2.01 | -1.53
SrBiO3 | YesKazakov  | 4.372 | 1.88 | -3.68
BaBiO3 | YesSleight  | 4.417 | 1.81 | -3.96
RaMcO3 | No | 4.676 | 1.46 | -8.86
CsTlF3 | YesRetuerto  | 4.799 | 1.24 | 0.65
CsTlCl3 | YesRetuerto  | 5.604 | 1.16 | -0.88

##  III Results and discussion

Let us first discuss the systematic series of $AB$O3 $s\textendash p$ cubic
perovskites, where the $A$ and $B$ cations are varied down the periodic table
as MgPO3, CaAsO3, SrSbO3, BaSbO3, SrBiO3, BaBiO3, and RaMcO3. Among them, only
SrBiO3, BaSbO3, and BaBiO3 exist in nature (see Table 1), but our prime
interest is to identify general trends in the electronic structure of such
$s\textendash p$ systems.

Figures 2 (a)-(e) and 4 (a)-(e) show the band-structures and the projected
densities of states (DOS) of the studied $AB$O3 series, respectively. Note
that no results for SrBiO3 and BaSbO3 are shown as they are very similar to
those for, respectively, BaBiO3 and SrSbO3. All the systems demonstrate the
same strong bonding–anti-bonding splitting between the $B\textendash s$ atomic
and O$\textendash a_{1g}$ molecular orbitals, with the anti-bonding band
landing at the Fermi level and becoming half occupied. Only in MgPO3, there is
a second band crossing the Fermi level, which is mainly of the Mg$\textendash
3s$ orbital character. As a general trend, the overall band width of the
$B\textendash s-{\text{O}}\textendash 2p$ states decreases as we go down from
MgPO3 to RaMcO3. Indeed, as the lattice constant increases due to the large
ionic radii of the $A$ and $B$ cations (see Table 1), the direct hopping
between oxygen orbitals is decreasing, which results in reduction of the
oxygen orbitals’ band width, and so does the hopping between oxygen orbitals
and the $B$ cation $s$ orbitals, which results in reduction of the band width
of the $B\textendash s$ and O$\textendash a_{1g}$ hybrid. On one hand, the
narrowing of the oxygen band makes it easier to push the $a_{1g}$ states up
and out of the top of the oxygen band. On the other hand, the reduced
$B\textendash s-{\text{O}}\textendash 2p_{\sigma}$ hopping integral
$t_{sp\sigma}$ leads to a flattening of the anti-bonding conduction band in
BaBiO3 and, especially, RaMcO3. Correspondingly, the density of states at the
Fermi level is strongly increased in these two end members of the series,
which makes them more strongly driven towards bond disproportionation and
other types of structural distortions. Another important trend in the band
structures of the $AB$O3 series is the gradual change of the dominating
character of the anti-bonding conduction band from one of more $B\textendash
s$ atomic orbital to one of more O$\textendash a_{1g}$ MO character. This is
illustrated in Fig. 2 (h) showing the amount of the O$\textendash a_{1g}$ MO
character in the anti-bonding band of the $AB$O3 TB models (which will be
discussed in more detail shortly), $|\langle
a_{1g}(\boldsymbol{k})|\psi^{\text{TB}}(\boldsymbol{k})\rangle|^{2}$, at the
$R$ point as a function of chemical composition.

It is important to note that the bonding–anti-bonding splitting is strongly
$\boldsymbol{k}$-vector dependent, due to the changing symmetry of the Bloch
functions involved. The splitting vanishes at $\Gamma$
[$\boldsymbol{k}=(0,0,0)$], because at this point the Bloch wave function has
no $a_{1g}$ molecular orbital component [Fig. 2 (f)]. In contrast, at the $R$
point [$\boldsymbol{k}=(\pi,\pi,\pi)$], the Bloch wave function is of a pure
$a_{1g}$ molecular orbital character in the oxygen $p_{\sigma}$ orbitals’
domain [Fig. 2 (g)], and the splitting reaches its maximum.

In accordance with this logic, the band structure plots of Fig. 2 feature no
$a_{1g}$ character in any of the Bloch states at the $\Gamma$ point. Also,
there is only one state with dominating $B\textendash s$ character at this
$\boldsymbol{k}$-vector, which moves to lower energies towards the end of the
$AB$O3 series, from around -5 eV in MgPO3 to -10 eV in RaMcO3. Along the
$\Gamma-R$ path (and, actually, at any point other than $\Gamma$), the finite
$B\textendash s-{\text{O}}\textendash 2p_{\sigma}$ hybridization is expected
to make the bonding band disperse downwards and reach a minimum at $R$. While
this is what we indeed observe in the band structures of BaBiO3 and RaMcO3,
the behavior of the bonding band in the band structures of MgPO3, CaAsO3 and
SrSbO3 appears to be more complex. There, the $B\textendash s$ state at
$\Gamma$ happens to be energetically above the triply degenerate oxygen
molecular orbital $t_{1u}$ states (positioned at around -11 eV in MgPO3 and -8
eV in CaAsO3 and SrSbO3). Away from $\Gamma$, these oxygen bands disperse
upwards and get entangled with the bonding band. This is the reason why in
MgPO3, CaAsO3 and SrSbO3 the character of the bonding $a_{1g}$ and
$B\textendash s$ combination is not continuous.

A similar situation occurs for the anti-bonding band close to the $R$ point.
At this point, in the same three systems, the interatomic orbital
hybridizations are strong enough to push the anti-bonding state above the
triply degenerate $B\textendash p$ states and position it at around 4 to 6 eV.
Away from $R$, the $B\textendash p$ bands and the anti-bonding band mix with
each other, and this is the reason why in these three systems the character of
the anti-bonding $a_{1g}$ and $B\textendash s$ combination is not continuous
either.

Let us now compare our findings about the $AB$O3 series with the calculated
electronic structures of CsTlF3 and CsTlCl3, shown in Fig. 3. These recently
synthesized halidesRetuerto  also demonstrate the bond disproportionation, and
so qualitatively we expect their electronic structure near $E_{\text{F}}$ to
be very similar to the oxides, with O replaced by halogen and the divalent
cation with monovalent Cs, also resulting in divalent Tl formally with one
electron in a $6s$ orbital, similar to the tetravalent Bi-based problem. It
remains unknown, however, whether hole doping can make the thallium halides
superconductRetuerto2 . While we note an overall similarity with the
electronic structure of the previously discussed $AB$O3 compounds, the
$\text{Tl}\textendash 6s-\text{halogen}\textendash a_{1g}$ band splitting is
considerably smaller. Also, even though the difference between the CsTlF3 and
CsTlCl3 lattice constants is quite significant (Table 1), we see only a slight
change in the band-widths of their $\text{Tl}\textendash
6s-\text{halogen}\textendash a_{1g}$ band manifolds.

We proceed now to quantifying the observed differences in the electronic
structures of the $ABX_{3}$ compounds in terms of the $B\textendash
s-X\textendash p$ hybridization strength $t_{sp\sigma}$ and the charge-
transfer energy $\Delta$. As was shown in Ref. Khazraie2, , these are the two
main parameters that determine the character of the empty states crossing the
Fermi level. We should emphasize that in our convention $\Delta$ is defined as
the difference between the on-site energies of the $B\textendash s$ atomic
orbital and of the $X\textendash a_{1g}$ MO: $\Delta=\epsilon(B\textendash
s)-\epsilon(X\textendash a_{1g})$Khazraie2 . $t_{sp\sigma}$ and $\Delta$ are
obtained by calculating MLWF-based TB models for our systems. Similarly to
Ref. Khazraie, , here we consider 10 orbitals per formula unit: one
$B\textendash s$ and three $X\textendash p$ orbitals per each of the three
anions in a simple cubic unit cell. As one can see in Figs. 2 (a)-(e) and 3
(a)-(b), the resulting TB models have band dispersions that agree well with
the GGA band structures. However, small deviation exists at the $R$ point of
MgPO3, CaAsO3, and SrSbO3, due to hybridization with higher energy states, as
discussed earlier in the paper. We have chosen not to consider the higher
energy $B\textendash p$ orbitals in the TB models because our main focus is on
the low-energy scale electron removal and addition states relevant for the
physical properties at relatively low temperatures.

The hopping integrals $t_{sp\sigma}$ are given in Table 1. While for the
$AB$O3 series the value of $t_{sp\sigma}$ decreases rather monotonically from
MgPO3 to RaMcO3 and is inversely related to the unit cell’s size, the
$t_{sp\sigma}$ values of the two halides are surprisingly similar given the
big difference between their unit cells’ sizes. In order to obtain the charge-
transfer energy $\Delta=\epsilon(B\textendash s)-\epsilon(X\textendash
a_{1g})$, we have applied a basis set transformation from oxygen atomic to
oxygen molecular orbitals using the table in Fig. 1 (b). The resulting charge-
transfer energy values $\Delta$ are listed in Table 1. In the $AB$O3 series,
$\Delta$ varies widely from a positive value of 1.83 eV in MgPO3 to a negative
value of -8.86 eV in RaMcO3, thus marking a difference between MgPO3 and
CsTlF3, as being _positive charge-transfer energy_ compounds, and the rest of
the $AB$O3 systems and CsTlCl3, as being _negative charge-transfer energy_
compounds. We note again that in all cases, except perhaps RaMcO3, the total
bonding–anti-bonding splitting strongly dominates over the charge transfer
energy, making the latter less important than the hopping integrals.

The dramatic decrease of $\Delta$ in the $AB$O3 series with $A$ and $B$ moving
down the periodic table can be understood in terms of relativistic lowering of
the $6s$ and $7s$ orbital energies in the heavy elements Bi and McPyykko . In
RaMcO3, it takes an extreme form, which, combined with the reduced
$t_{sp\sigma}$ hybridization, results in an almost ionic character of the Mc3+
ionization state, in agreement with earlier studies Keller . This can be
clearly seen from the RaMcO3 projected DOS shown in Fig. 4 (e). Similarly to
Ba(Sr)BiO3, RaMcO3 is therefore expected to bond-disproportionate in its
ground state as
$2\text{Mc}^{3+}\underline{L}^{2}\rightarrow[\text{Mc}^{3+}]_{\text{large}}+[\text{Mc}^{3+}\underline{L}^{2}]_{\text{small}}$,
with most of the action happening on oxygens. On the other hand, also
interesting is the fact that BaSbO3 can become superconducting upon
substituting Sb with Pb but not with SnCava2 . This again was explained in
terms of the relativistic energy lowering of Pb$\textendash 6s$ states with
respect to Sn$\textendash 5s$ statesSingh . The maximal $T_{\text{c}}$
observed in BaPbxSb1-xO3 is 3.5 K, which is unexpectedly low compared with
that of the bismuthates. At this point, we can only speculate that if there is
a relation between $T_{\text{c}}$ and how negative $\Delta$ is, then optimally
hole-doped RaMcO3 might have a very high $T_{\text{c}}$. Such a relation can
be rooted in electron-phonon coupling which might become most efficient in the
presence of ligand holes, but of course more theoretical and experimental
studies are required to give some validity to this idea. We should also note
that the density of states at $E_{\text{F}}$ is very high in RaMcO3 because of
the narrowing of the anti-bonding band and this could also support a higher
$T_{\text{c}}$.

Figure 5:  The phase diagram of Ref. Khazraie2,  representing the dominant
character of holes as a function of the charge-transfer energy $\Delta$ and
hybridization $t_{sp\sigma}$. Symbols $+$ and $\times$ mark the parameters
relevant to the hypothetical and existing cubic $s\textendash p$ $ABX_{3}$
perovskites, respectively. Green, yellow, and red colors represent the amount
of the $X\textendash e_{g}$, $X\textendash a_{1g}$, and $B$ cation $s$ orbital
contributions to the holes’ character, respectively.

Finally, with the hopping integrals $t_{sp\sigma}$ and charge-transfer
energies $\Delta$ of our studied $ABX_{3}$ compounds at hand, we can mark
their positions on the phase diagram proposed in Ref. Khazraie2, . There will
be some degree of approximation involved because this phase diagram was
obtained for BaBiO3 in a bond-disproportionated state using the value of the
oxygen band-width $W$ specific to BaBiO3, but this is not going to obscure
observation of general trends that we are most interested in. As explained in
Ref. Khazraie2, , the colors in this phase diagram represent the dominant
character of the empty (hole) states above the Fermi level as a function of
hybridization and charge-transfer energy. While there is a sharp boundary
around the green region where holes reside on $X\textendash e_{g}$ orbitals,
the $B\textendash s$ (yellow) and $X\textendash a_{1g}$ (red) orbitals are
always mixed by hybridization and therefore there is a gradual crossover
between the yellow and red regions. The black dashed line marks equal
contributions from the $B\textendash s$ and $X\textendash a_{1g}$ orbitals to
the holes’ character. Figure 5 is showing now that SrBiO3 and BaBiO3 are
relatively deep in the $X\textendash a_{1g}$ region, while MgPO3, CaAsO3,
SrSbO3, BaSbO3, and CsTlF3 are close to having equal $B\textendash s$ and
$X\textendash a_{1g}$ orbital contributions. As for RaMcO3 and CsTlCl3, they
land in the O$\textendash e_{g}$ region but from their projected DOS [Fig. 4
(e) and Fig. 3 (d)] we know that their hole character is strongly
O$\textendash a_{1g}$. This discrepancy is due to the mentioned
approximations, but it is obvious that RaMcO3 must in any case be located very
far inside the O$\textendash a_{1g}$ region close to the O$\textendash e_{g}$
border.

##  IV Conclusions

In this paper, we have used ab initio methods to study the electronic
structures of the following $s\textendash p$ cubic perovskites $ABX_{3}$: the
experimentally available BaSbO3, SrBiO3, BaBiO3, CsTlF3, and CsTlCl3, as well
as the hypothetical MgPO3, CaAsO3, SrSbO3, and RaMcO3. We have used Wannier
functions based tight-binding modeling to calculate their hybridization
strengths $t_{sp\sigma}$ between the $B\textendash s$ and $X\textendash p$
atomic orbitals and charge-transfer energies $\Delta$, which are the two most
important parameters that determine the nature of the systems’ holes. These
calculations have elucidated several trends in $t_{sp\sigma}$ and $\Delta$ as
one moves across the periodic table, such as the relativistic energy lowering
of the $B\textendash s$ orbital in heavy $B$ cations leading to strongly
negative $\Delta$ values. Our results have been discussed in connection with
the general phase diagram for $s\textendash p$ cubic perovskites proposed in
Ref. Khazraie2, .

Also, some considerations have been offered regarding a possible relation
between the highest achievable superconducting transition temperatures and
certain features of the systems’ electronic structures, such as the charge-
transfer energy $\Delta$ and the interatomic orbital hybridization
$B\textendash s-{\text{O}}\textendash 2p_{\sigma}$ $t_{sp\sigma}$, with the
latter primarily controlling the conduction band width and the DOS at the
Fermi level. In particular, we observed that the more negative the $\Delta$
value is the more of the $X\textendash a_{1g}$ character will be there in the
valence and conduction bands provided that the hybridization strength is large
enough to still push out the $X\textendash a_{1g}$ bound state. In the case of
RaMcO3 and also Ba(Sr)BiO3, $\Delta$ is very negative and so the amount of the
$X\textendash a_{1g}$ character in the anti-bonding state will be the largest.
As $\Delta$ becomes even more negative there will be a point where the
$X\textendash a_{1g}$ state no longer is pushed above the $X\textendash e_{g}$
state, as illustrated in the phase diagram of Fig. 5. Given the significantly
larger maximal $T_{\text{c}}$ in the hole-doped bismuthates Ba(Sr)BiO3 than in
the hole-doped antimonates BaSbO3, superconductivity is apparently enhanced if
the charge character of the hole states is mostly $X\textendash a_{1g}$, which
can possibly be traced down to the enhanced coupling with the breathing mode
phonon. However, an even more important feature of the electronic structure
that can be responsible for an enhanced $T_{\text{c}}$ is the increase of the
DOS at the Fermi level in BaBiO3, due to its reduced $t_{sp\sigma}$ value. For
these reasons, we expect that optimally hole-doped RaMcO3 might have a very
high $T_{\text{c}}$.

## Acknowledgments

This work was supported by Natural Sciences and Engineering Research Council
(NSERC) for Canada, CIfAR, and the Max Planck - UBC Stewart Blusson Quantum
Matter Institute.

## References

  * (1) Imada et al. Rev. Mod. Phys. 70, 1039(1998).
  * (2) M. L. Medarde, J. Phys. Condens. Matter 129, 1679 (1997).
  * (3) J. G. Bednorz and K. A. Muller, Zeitschrift fur Physik B Condensed Matter 64, 189 (1986).
  * (4) A. W. Sleight, J. L. Gillson, and P. E. Bierstedt, Solid State Commun. 17, 27 (1975).
  * (5) R. J. Cava, B. Batlogg, J. J. Krajewski, R. Farrow, L. W. Rupp, Jr., A. E. White, K. Short, W. F. Peck, and T. Kometani, Nature 332, 814 (1988).
  * (6) A. M. Green; A. Ho-Baillie; H. J. Snaith; Nature Photonics, 8(7), 506-514 (2014).
  * (7) A. P. Ramirez, J. Phys.: Condens. Matter 9, 8171(1997).
  * (8) T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, Y. Tokura, Nature 426, 58(2003).
  * (9) H. Jin, S. H. Rhim, J. Im and A. J. Freeman, Scientific Reports 3:1651, 1(2012).
  * (10) H. Jin, J. Im, and A. J. Freeman, Phys. Rev. B 86, 121102(R) (2012).
  * (11) B. Yan, M. Jansen and C. Felser, Nature Phys. 9, 709(2013).
  * (12) L. R. Mattheiss, E. M. Gyorgy, and D. W. Johnson, Phys. Rev. B37, 3746 (1988).
  * (13) S. M. Kazakov, C. Chaillout, P. Bordet, J. J. Capponi, M. Nunez-Regueiro, A. Rysak, J. L. Tholence, P. G. Radaelli, S. N. Putilin, and E. V. Antipov, Nature 390, 148 (1997).
  * (14) D. Cox and A. W. Sleight, Solid State Commun. 19, 969 (1976).
  * (15) D. E. Cox and A. W. Sleight, Acta Crystallographica Section B 35, 1 (1979).
  * (16) Q. Zhou, B. J. Kennedy Solid State Commun. 132, 392(2004).
  * (17) A. Sleight, Physica C 514, 152 (2015).
  * (18) C. M. Varma, Phys. Rev. Lett. 61, 2713 (1988).
  * (19) A. Taraphder, H. R. Krishnamurthy, R. Pandit, and T.V. Ramakrishnan, Negative-U Extended Hubbard Model for Doped Barium Bismuthates, Phys. Rev. B 52, 1368 (1995).
  * (20) I. Hase and T. Yanagisawa, Phys. Rev. B 76, 174103 (2007).
  * (21) L. F. Mattheiss and D. R. Hamann, Phys. Rev. B 28, 4227 (1983).
  * (22) W. A. Harrison, Phys. Rev. B 74, 245128 (2006).
  * (23) K. Foyevtsova, A. Khazraie, I. Elfimov and G. Sawatzky, Phys. Rev. B 91, 121114(R) (2015).
  * (24) A. Khazraie, K. Foyevtsova, I. Elfimov, and G. A. Sawatzky, Phys. Rev. B97, 075103 (2018).
  * (25) G. M. Dalpian, Q. Liu, J. Varignon, M. Bibes, and Alex Zunger, Phys. Rev. B 98, 075135 (2018).
  * (26) A. Khazraie, K. Foyevtsova, I. Elfimov, and G. A. Sawatzky, Phys. Rev. B 98, 205104 (2018).
  * (27) J. de Hair and G. Blasse, Solid State Commun. 12, 727 (1973).
  * (28) A. F. Orchard and G. Thornton, J. Chem. Soc. Dalton Trans., 1238 (1977).
  * (29) G. K. Wertheim, J. P. Remeika, and D. N. E. Buchanan, Phys. Rev. B 26, 2120 (1982).
  * (30) S. Salem-Sugui, Jr., E. E. Alp, S. M. Mini, M. Ramanathan, J. C. Campuzano, G. Jennings, M. Faiz, S. Pei, B. Dabrowski, Y. Zheng, D. R. Richards, and D. G. Hinks, Phys. Rev. B43, 5511-5515(1991).
  * (31) N. C. Plumb, D. J. Gawryluk, Y. Wang, Z. Ristic, J. Park, B. Q. Lv, Z. Wang, C. E. Matt, N. Xu, T. Shang, K. Conder, J. Mesot, S. Johnston, M. Shi, and M. Radovic, Phys. Rev. Lett. 117, 037002 (2016).
  * (32) S. Balandeh, R. J. Green, K. Foyevtsova, S. Chi, O. Foyevtsov, F. Li, and G. A. Sawatzky, Phys. Rev. B 96, 165127 (2017).
  * (33) H. Raebiger, S. Lany, and A. Zunger, Charge self-regulation upon changing the oxidation state of transition metals in insulators. Nature 453, 763?766 (2008).
  * (34) M. Retuerto et al., Chem. Mater. 25, 4079 (2013).
  * (35) D. J. Singh, D. A. Papaconstantopoulos, J. P. Julien, F. Cyrot-Lackmann, Phys. Rev. B44, 9519(1991).
  * (36) R. J. Cava et al., Nature 339, 293 (1989).
  * (37) W. Kohn, L.J. Sham, Phys. Rev. 140, A1133(1965).
  * (38) P. Blaha, K. Schwarz, G. K. H. Madsen, D.Kvasnicka, and J. Luitz, WIEN2K, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties (Karlheinz Schwarz, Techn. Universit?at Wien, Austria, 2001).
  * (39) J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77, 3865(1996).
  * (40) K. Foyevtsova and G. A. Sawatzky, Journal of Modern Physics 10, 953?965 (2019).
  * (41) A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza,D. Vanderbilt, and N. Marzari, Computer Physics Commun. 178, 685 (2008).
  * (42) M. Retuerto et al., Inorg. Chem. 54, 1066?1075 (2015).
  * (43) P. Pyykko, Chem. Rev. 88, 563 (1988).
  * (44) O. L. Keller Jr., C. W. Nestor Jr., and Burkhard Fricke. J. Phys. Chem., 78(19), 1949 (1974).