* [23] Yue Yu, Yong-Shi Wu, and Xincheng Xie.  Bulk–edge correspondence, spectral flow and atiyah–patodi–singer theorem for the z2-invariant in topological insulators.  Nuclear Physics B, 916:550–566, 2017.
  * [24] Guillaume Bal.  Topological charge conservation for continuous insulators.  arXiv:2106.08480, 2021.
  * [25] Y. X. Zhao and Shengyuan A. Yang.  Index theorem on chiral landau bands for topological fermions.  Phys. Rev. Lett., 126:046401, Jan 2021.
  * [26] M Nakahara.  Geometry, Topology and Physics.  CRC Press 2003, 2003.
  * [27] Barry Bradlyn, Jennifer Cano, Zhijun Wang, M. G. Vergniory, C. Felser, R. J. Cava, and B. Andrei Bernevig.  Beyond dirac and weyl fermions: Unconventional quasiparticles in conventional crystals.  Science, 353(6299), 2016.
  * [28] F. D. M. Haldane.  Model for a quantum hall effect without landau levels: Condensed-matter realization of the "parity anomaly".  Phys. Rev. Lett., 61:2015–2018, Oct 1988.
  * [29] B. Andrei Bernevig.  Topological Insulators and Topological Superconductors.  Princeton University Press, 2013, 2013.
  * [30] Kai-Yu Yang, Yuan-Ming Lu, and Ying Ran.  Quantum hall effects in a weyl semimetal: Possible application in pyrochlore iridates.  Phys. Rev. B, 84:075129, Aug 2011.
  * [31] D. T. Son and B. Z. Spivak.  Chiral anomaly and classical negative magnetoresistance of weyl metals.  Phys. Rev. B, 88:104412, Sep 2013.
  * [32] L. Balents.  Viewpoint: Weyl electron kiss.  Physics, 4:36, 2011.
  * [33] A. A. Burkov and Leon Balents.  Weyl semimetal in a topological insulator multilayer.  Phys. Rev. Lett., 107:127205, Sep 2011.
  * [34] Claude Cohen-Tannoudji, Bernard Diu, and Franck Laloë.  Quantum Mechanics, volume 1.  Hermann, 1998.
  * [35] K. Landsteiner.  Notes on anomaly induced transport.  Acta Phys. Pol. B, 47:2617, 2016.
  * [36] S. Raghu and F. D. M. Haldane.  Analogs of quantum-hall-effect edge states in photonic crystals.  Phys. Rev. A, 78:033834, Sep 2008.
  * [37] M. F. Atiyah and I. M. Singer.  The index of elliptic operators on compact manifolds.  Bull. Amer. Math. Soc., 69:422, 1963.
  * [38] Giampiero E.  Dirac Operators and Spectral Geometry (Cambridge Lecture Notes in Physics, Series Number 12).  Cambridge University Press, 1998.
  * [39] M. F. Atiyah, V. K. Patodi, and I. M. Singer.  Spectral asymmetry and riemannian geometry. iii.  Mathematical Proceedings of the Cambridge Philosophical Society, 79(1):71–99, 1976.
  * [40] Frédéric Faure and Boris Zhilinskii.  Topological properties of the born–oppenheimer approximation and implications for the exact spectrum.  Letters in Mathematical Physics, 55(3):219–238, 2001.
  * [41] G. B. Folland.  Harmonic Analysis in Phase Space.  Princeton University Press; F First Edition, 1989.
  * [42] M. Zworski.  Semiclassical Analysis, volume 138.  American Mathematical Society, Graduate Studies in Mathematics, 2012.
  * [43] A. Bohm, A. Mostafazadeh, H. Koizumi, Niu Q., and Zwanziger J.  The Geometric Phase in Quantum Systems.  Springer, 2010.
  * [44] M. V. Berry.  Quantal phase factors accompanying adiabatic changes.  Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 392, 1984.
  * [45] Barry Simon.  Holonomy, the quantum adiabatic theorem, and berry’s phase.  Phys. Rev. Lett., 51:2167–2170, Dec 1983.
  * [46] Ganesh Sundaram and Qian Niu.  Wave-packet dynamics in slowly perturbed crystals: Gradient corrections and berry-phase effects.  Phys. Rev. B, 59:14915–14925, Jun 1999.
  * [47] Di Xiao, Ming-Che Chang, and Qian Niu.  Berry phase effects on electronic properties.  Rev. Mod. Phys., 82:1959–2007, Jul 2010.
  * [48] Gregor Jotzu, Michael Messer, Rémi Desbuquois, Martin Lebrat, Thomas Uehlinger, Daniel Greif, and Tilman Esslinger.  Experimental realization of the topological haldane model with ultracold fermions.  Nature, 515(7526):237–240, 2014.
  * [49] Martin Wimmer, Hannah M. Price, Iacopo Carusotto, and Ulf Peschel.  Experimental measurement of the berry curvature from anomalous transport.  Nature Physics, 13(6):545–550, 2017.
  * [50] Tran Trung Luu and Hans Jakob Wörner.  Measurement of the berry curvature of solids using high-harmonic spectroscopy.  Nature Communications, 9(1):916, 2018.
  * [51] Soohyun Cho, Jin-Hong Park, Jisook Hong, Jongkeun Jung, Beom Seo Kim, Garam Han, Wonshik Kyung, Yeongkwan Kim, S.-K. Mo, J. D. Denlinger, Ji Hoon Shim, Jung Hoon Han, Changyoung Kim, and Seung Ryong Park.  Experimental observation of hidden berry curvature in inversion-symmetric bulk $2h\text{$-$}{\mathrm{wse}}_{2}$.  Phys. Rev. Lett., 121:186401, Oct 2018.
  * [52] Yasuhiro Hatsugai.  Chern number and edge states in the integer quantum hall effect.  Phys. Rev. Lett., 71:3697–3700, Nov 1993.
  * [53] Yasuhiro Hatsugai.  Edge states in the integer quantum hall effect and the riemann surface of the bloch function.  Phys. Rev. B, 48:11851–11862, Oct 1993.
  * [54] G.M. Graf and M Porta.  Bulk–edge correspondence for two-dimensional topological insulators.  Commun. Math. Phys, 324:851, 2013.
  * [55] Wang Yao, Shengyuan A. Yang, and Qian Niu.  Edge states in graphene: From gapped flat-band to gapless chiral modes.  Phys. Rev. Lett., 102:096801, Mar 2009.
  * [56] F. Zhang, A. H. MacDonald, and E. J. Mele.  Valley chern numbers and boundary modes in gapped bilayer graphene.  Proc. Natl. Acad. Sci. USA, 110:10546, 2013.
  * [57] Guancong Ma, Meng Xiao, and C. T. Chan.  Topological phases in acoustic and mechanical systems.  Nature Reviews Physics, 1(4):281–294, 2019.
  * [58] Lavi K. Upreti and Pierre Delplace.  Topological chiral interface states beyond insulators.  Phys. Rev. A, 102:023520, Aug 2020.
  * [59] Manolis Perrot, Pierre Delplace, and Antoine Venaille.  Topological transition in stratified fluids.  Nature Physics, 15(8):781–784, 2019.
  * [60] M. Marciani and P. Delplace.  Chiral maxwell waves in continuous media from berry monopoles.  Phys. Rev. A, 101:023827, Feb 2020.
  * [61] Antoine Venaille and Pierre Delplace.  Topology brought to the coast.  arXiv preprint arXiv:2011.03440, 2020.
  * [62] C. Tauber, P. Delplace, and A. Venaille.  Anomalous bulk-edge correspondence in continuous media.  Phys. Rev. Research, 2:013147, Feb 2020.
  * [63] J. B. Pendry, L. Martín-Moreno, and F. J. Garcia-Vidal.  Mimicking surface plasmons with structured surfaces.  Science, 305(5685):847–848, 2004.
  * [64] U. Fano.  Atomic theory of electromagnetic interactions in dense materials.  Phys. Rev., 103:1202–1218, Sep 1956.
  * [65] J. J. Hopfield.  Theory of the contribution of excitons to the complex dielectric constant of crystals.  Phys. Rev., 112:1555–1567, Dec 1958.
  * [66] F Lemoult, N. Kaina, M. Fink, and G Lerosey.  Soda cans metamaterial: A subwavelength-scaled phononic crystal.  Crystals, 6:82, 2016.
  * [67] S. Yves, R. Fleury, F. Lemoult, M. Fink, and G. Lerosey.  Topological acoustic polaritons: robust sound manipulation at the subwavelength scale.  New J. Phys., 19:075003, 2017.
  * [68] Florent Baboux, Eli Levy, Aristide Lemaître, Carmen Gómez, Elisabeth Galopin, Luc Le Gratiet, Isabelle Sagnes, Alberto Amo, Jacqueline Bloch, and Eric Akkermans.  Measuring topological invariants from generalized edge states in polaritonic quasicrystals.  Phys. Rev. B, 95:161114, Apr 2017.
  * [69] Yaroslav V. Kartashov and Dmitry V. Skryabin.  Two-dimensional topological polariton laser.  Phys. Rev. Lett., 122:083902, Feb 2019.
  * [70] F Faure.  Manifestation of the topological index formula in quantum waves and geophysical waves.  arXiv:1901.10592, 2019.
  * [71] J. E. Avron, L. Sadun, J. Segert, and B. Simon.  Chern numbers, quaternions, and berry’s phases in fermi systems.  Commun. Math. Phys., 124:595–627, 1989.
  * [72] B.A. Dubrovin, A.T. Fomenko, and S.P. Novikov.  Modern Geometry— Methods and Applications: Part II: The Geometry and Topology of Manifolds.  Springer Science and Business Media, 1985.
  * [73] Pierre Delplace.  Habilitation thesis : A few aspects of topological waves, 2019.
  * [74] Gian Michele Graf, Hansueli Jud, and Clément Tauber.  Topology in shallow-water waves: A violation of bulk-edge correspondence.  Communications in Mathematical Physics, 383(2):731–761, 2021.
  * [75] Neal H. McCoy.  On the function in quantum mechanics which corresponds to a given function in classical mechanics.  Proceedings of the National Academy of Sciences, 18(11):674–676, 1932.
  * [76] Constantine Callias.  Axial anomalies and index theorems on open spaces.  Communications in Mathematical Physics, 62(3):213–234, 1978.
  * [77] Erick J. Weinberg.  Index calculations for the fermion-vortex system.  Phys. Rev. D, 24:2669–2673, Nov 1981.
  * [78] David H. Adams.  Theoretical foundation for the index theorem on the lattice with staggered fermions.  Phys. Rev. Lett., 104:141602, Apr 2010.
  * [79] Nicole Berline, Ezra Getzler, and Michéle Vergne.  Heat Kernels and Dirac Operators.  Springer-Verlag Berlin Heidelberg NewYork, 1996.
  * [80] Taroh Matsuno.  Quasi-geostrophic motions in the equatorial area.  Journal of the Meteorological Society of Japan. Ser. II, 44(1):25–43, 1966.
  * [81] Jeffrey B. Parker, J. W. Burby, J. B. Marston, and Steven M. Tobias.  Nontrivial topology in the continuous spectrum of a magnetized plasma.  Physical Review Research, 2(3):033425–, 09 2020.
  * [82] Shin-Ming Huang, Su-Yang Xu, Ilya Belopolski, Chi-Cheng Lee, Guoqing Chang, Tay-Rong Chang, BaoKai Wang, Nasser Alidoust, Guang Bian, Madhab Neupane, Daniel Sanchez, Hao Zheng, Horng-Tay Jeng, Arun Bansil, Titus Neupert, Hsin Lin, and M. Zahid Hasan.  New type of weyl semimetal with quadratic double weyl fermions.  Proceedings of the National Academy of Sciences, 113(5):1180–1185, 2016.
  * [83] M. Zahid Hasan, Guoqing Chang, Ilya Belopolski, Guang Bian, Su-Yang Xu, and Jia-Xin Yin.  Weyl, dirac and high-fold chiral fermions in topological quantum matter.  Nature Reviews Materials, 6(9):784–803, 2021.
  * [84] T. Frankel.  The Geometry of Physics.  Cambridge University Press, third edition, 2012.
  * [85] G. K. Vallis.  Atmospheric and Oceanic Fluid Dynamics : Fundamentals and Large-Scale Circulation.  Cambridge University Press, Cambridge, U. K., 2nd edition edition, 2017\.

##  Appendix A Operator-Symbol correspondence

In quantum mechanics position and momentum are described by operators
$\hat{x}$ and $\hat{p}$, while in Hamiltonian classical mechanics, those are
scalar conjugate variables in phase space $(x,p)$. A crucial difference
between the two theories being that $\hat{x}$ and $\hat{p}$ do not commute
while $x$ and $p$ do. The mapping from phase space to Hilbert space of
operators is called quantization. This procedure, which was at the root of the
construction of quantum mechanics from classical mechanics during the first
decades of the XXth century, stimulated many mathematical works. Indeed such a
mapping is not unique, and many quantization procedures where developed called
e.g. $x-p$ and $p-x$ quantizations, Wick and anti-Wick quantizations,
geometric quantization, and so on. The inverse procedure, that is the mapping
from Hilbert space to phase space of classical observables is also routinely
used in physics, with for instance the WKB approximation, that consists in
expanding a wave function with respect to a small parameter, and the Wigner
transform that allows the representation of a quantum state in phase space.
Those procedures are based on expansions with respect to a small parameter,
namely the Planck constant $h$. Actually, such procedures are not restricted
to quantum physics, but much more broadly to any field that deals with
differential equations with a small parameter. This field is referred to semi-
classical analysis, or also micro-local analysis in the mathematical context.
When an operator of Hilbert space is mapped to a function on phase space, this
resulting ”classical” function is called the symbol. So in the following, a
symbol will refer to a function $\sigma$ on phase space, and will be
associated to an operator $\text{Op}[\sigma]\equiv\hat{\sigma}$ through a
quantization procedure. Conversely, we shall use the notation
$\sigma[\hat{\sigma}]=\sigma$ to indicate that we take the symbol of an
operator.

In practice, we will be interested in cases where there is a one-to-one
correspondence between a symbol and its operator. For this, we use a standard
quantization procedure, called Weyl quantization. Weyl quantization turns a
function $\sigma(x,p)$ on phase space into an operator $\hat{\sigma}$, through
the expression of the action of this operator on an arbitrary wave function in
position representation $\psi(x)$, called the Weyl transform and that reads

$\displaystyle\hat{\sigma}\psi(x)=\frac{1}{2\pi\epsilon}\int_{\mathbb{R}}{\rm
d}x^{\prime}\int_{\mathbb{R}}{\rm
d}\xi\,\,\sigma\left(\frac{x+x^{\prime}}{2},p\right){\rm e}^{{\rm
i}\xi(x-x^{\prime})/\epsilon}\psi(x^{\prime})\ .$ (155)

This formula is given for 1D systems, but generalizes straightforwardly at any
dimension. Note the existence of a parameter $\epsilon$ in the definition.
Strictly speaking, a quantization procedure is only well-defined with respect
to such a parameter, in order to recover a classical behavior of the dynamics
in the limit $\epsilon\rightarrow 0$. The Weyl transform (155) yields

$\displaystyle\sigma(x,p)=x\quad\rightarrow\quad\hat{\sigma}\psi(x)=x\psi(x)$
(156)
$\displaystyle\sigma(x,p)=p\quad\rightarrow\quad\hat{\sigma}\psi(x)=-{\rm
i}\epsilon\partial_{x}\psi(x)$ (157)

which reproduces the expected quantization rule of position $x$ and momentum
$p$, that is usually written as the substitution

$\displaystyle x\quad\rightarrow\quad\text{Op}[x]=\hat{x}=x$ (158)
$\displaystyle p\quad\rightarrow\quad\text{Op}[p]=\hat{p}=-{\rm
i}\epsilon\partial_{x}$ (159)

in position representation.

Conversely, the symbol $\sigma(x,p)$ of an operator
$\hat{\sigma}(\hat{x},\hat{p})$ can be obtained, at least formally, from the
Wigner transform

$\displaystyle\sigma(x,p)=\int{\rm
d}x^{\prime}\Sigma\left(\frac{x+x^{\prime}}{2},x-\frac{x-x^{\prime}}{2}\right){\rm
e}^{-{\rm i}px^{\prime}/\epsilon}$ (160)

where $\Sigma(x,x^{\prime})$ is the kernel of the operator $\hat{\sigma}$,
that is, it verifies

$\displaystyle\hat{\sigma}\psi(x)=\int{\rm
d}x^{\prime}\Sigma(x,x^{\prime})\psi(x^{\prime})\ .$ (161)

Equations (155), (160) and (161) yield a self-consistent definition of the
Wigner-Weyl calculus. Computing the symbol of an operator via the Wigner
transform necessitates the knowledge of the Kernel of this operator, which, in
general, is a distribution. This is not always the easiest way to compute the
symbol. In practice, one can use the powerful formula of the $\star$ product
between symbols $\sigma_{1}(x,p)$ and $\sigma_{2}(x,p)$

$\displaystyle\sigma_{1}(x,p)\star\sigma_{2}(x,p)\equiv\sum_{m,n=0}^{\infty}\frac{1}{n!}\left(\frac{{\rm
i}\epsilon}{2}\right)^{n}(-1)^{m}\begin{pmatrix}n\\\
m\end{pmatrix}\left(\partial_{x}^{n-m}\partial_{p}^{m}\sigma_{1}\right)\left(\partial_{p}^{n-m}\partial_{x}^{m}\sigma_{2}\right)$
(162)

thats is related to the product of their operators as

$\displaystyle\hat{\sigma}_{1}\hat{\sigma}_{2}=\text{Op}[\sigma_{1}\star\sigma_{2}]\
.$ (163)

Using the fact that identity verifies $\text{Op}\mathbb{I}=\mathbb{I}$, the
symbol of the operators $\hat{x}$ and $\hat{p}$ are then easily obtained by
using (162) and (163) and one gets

$\displaystyle x\quad\leftarrow\quad\text{Op}[x]=\hat{x}=x$ (164)
$\displaystyle p\quad\leftarrow\quad\text{Op}[p]=\hat{p}=-{\rm
i}\epsilon\partial_{x}$ (165)

as expected. When $\epsilon<<1$, this formula can be used to approximate a
product of operators as the operator of a product of their symbols in powers
of $\epsilon$. In particular, it yields that the first order correction in
$\epsilon$ to
$\hat{\sigma}_{1}\hat{\sigma}_{2}=\text{Op}[\sigma_{1}\sigma_{2}]$ is given by
the operator of the Poisson bracket $\\{\sigma_{1},\sigma_{2}\\}$, that is
precisely the commutator $[\hat{\sigma}_{1},\hat{\sigma}_{2}]$. For our
purpose however, the operators we have at hand are polynomials in $\hat{x}$
and $\hat{p}$, so that the sum in the $\star$ product (162) yields an exact
result with a finite number of terms. We shall consider all those terms, and
thus take $\epsilon=1$ from now. This yields the correspondence

$\displaystyle x+{\rm i}p\quad\leftrightarrow\quad{\hat{a}}$ (166)
$\displaystyle x-{\rm i}p\quad\leftrightarrow\quad{\hat{a}^{\dagger}}$ (167)

where $x$ and $p$ designate two canonical conjugate observables in phase
space, not necessarily position and impulsion. Note that the relation (163) is
particularly useful to compute the symbol of a product of operators. By taking
formally the symbol of each members of this equation, one has

$\displaystyle\sigma[\hat{\sigma}_{1}\hat{\sigma}_{2}]=\sigma_{1}\star\sigma_{2}\
.$ (168)

We can use the $\star$ product to establish the correspondences (86) and (87)
between ${H_{S}}$ and ${\hat{\mathcal{H}}}_{\text{op}}$ for the $h.\hat{S}$
model. To demonstrate this relation, let us proceed by induction. First, we
know that

$\displaystyle\text{Op}[\alpha x+\beta p]=\alpha\hat{x}+\beta\hat{p}$ (169)

where $\alpha$ and $\beta$ are arbitrary complex numbers. Then we assume the
relation

$\displaystyle\text{Op}[(\alpha x+\beta
p)^{m-1}]=(\alpha\hat{x}+\beta\hat{p})^{m-1}$ (170)

for a given $m-1$. Let us now prove that this relation remains true when
$m-1\rightarrow m$. We have

$\displaystyle(\alpha\hat{x}+\beta\hat{p})^{m}=\ $
$\displaystyle(\alpha\hat{x}+\beta\hat{p})^{m-1}(\alpha\hat{x}+\beta\hat{p})$
(171) $\displaystyle=\ $ $\displaystyle\text{Op}[(\alpha x+\beta
p)^{m-1}]\text{Op}[\alpha x+\beta p]$ (172) $\displaystyle=\ $
$\displaystyle\text{Op}[(\alpha x+\beta p)^{m-1}\star(\alpha x+\beta p)]$
(173) $\displaystyle=\ $
$\displaystyle\text{Op}\Big{[}\sum_{m,n=0}^{\infty}\frac{1}{n!}\left(\frac{{\rm
i}\epsilon}{2}\right)^{n}(-1)^{m}\begin{pmatrix}n\\\ m\end{pmatrix}$
$\displaystyle\times$ $\displaystyle\
\left(\partial_{x}^{n-m}\partial_{p}^{m}(\alpha x+\beta
p)^{m-1}\right)\left(\partial_{p}^{n-m}\partial_{x}^{m}(\alpha x+\beta
p)\right)\Big{]}\ .$ (174)

The last term in (174) is nonzero only when $(n,m)$ takes the values $(0,0)$,
$(1,0)$ and $(1,1)$. The contribution $(0,0)$ yields the term $(\alpha x+\beta
p)$, while the two other contributions $(1,0)$ and $(1,1)$ give two terms
proportional to $\epsilon$, such that

$\begin{split}(\alpha\hat{x}+\beta\hat{p})^{m}=\ \text{Op}[&(\alpha x+\beta
p)^{m-1}(\alpha x+\beta p)\\\ &+\frac{{\rm i}\epsilon}{2}(m-1)\alpha(\alpha
x+\beta p)^{m-2}\beta-\alpha(m-1)\beta(\alpha x+\beta p)^{m-2}]\ .\end{split}$
(175)

The two terms proportional to $\epsilon$ compensate each other, and we end up
with

$\displaystyle\text{Op}[(\alpha x+\beta
p)^{m}]=(\alpha\hat{x}+\beta\hat{p})^{m}$ (176)

which completes the proof. This formula is a famous result of Weyl calculus,
that was first obtained long ago by McCoy before the $\star$ product was
introduced [75]. Conversely, we also get that

$\displaystyle\sigma[(\alpha\hat{x}+\beta\hat{p})^{m}]=(\alpha x+\beta p)^{m}$
(177)

which leads to the relations (86) and (87) when $\alpha=1$ and $\beta=\pm{\rm
i}$.

##  Appendix B Toolbox of differential calculus: Stokes theorem and Brouwer
degree formula

There are many good textbooks that give detailed and consistent introductions
to differential forms, such as [84, 43, 26] that are dedicated to physicists.
Here we present a digest summary that sketches a few basic definitions and
results that are useful for the following.

Differential forms constitute a generalization of functions. A usual function
${\boldsymbol{\lambda}}:\Lambda\rightarrow
f({\boldsymbol{\lambda}})\equiv\omega_{0}$ is then a 0-form, while its
differential (provided it is differentiable) ${\rm d}f$ is an example of a
1-form. More generally, a 1-form is a kind of vector, called covector as it
transforms as a covariant vector [84], and decomposes as
$\omega_{1}=\omega^{i}_{1}{\rm d}\lambda_{i}$ (where we use the implicit sum
convention). If the dimension of parameter space $\Lambda$, that is our base
space, is $N$, then the set of 1-forms $\Omega^{1}(\Lambda)$ is also of
dimension $N$. A 2-form $\omega_{2}$ is a anti-symmetric tensor that
decomposes as $\omega_{2}=\omega_{jk}{\rm d}\lambda_{j}\wedge{\rm
d}\lambda_{k}$ (with $j<k$), where the wedge product $\wedge$ generalizes the
vector product defined in $\mathbb{R}^{3}$ to any dimension, by satisfying

$\displaystyle{\rm d}\lambda_{j}\wedge{\rm d}\lambda_{k}=-{\rm
d}\lambda_{k}\wedge{\rm d}\lambda_{j}$ (178)

and in particular ${\rm d}\lambda_{j}\wedge{\rm d}\lambda_{j}=0$. By extension
an r-form is a totally anti-symmetric tensor that decomposes as
$\omega_{r}=\omega_{j_{1},j_{2},\cdots,j_{r}}^{r}{\rm
d}\lambda_{j_{1}}\wedge{\rm d}\lambda_{j_{2}}\wedge\cdots\wedge{\rm
d}\lambda_{j_{r}}$, and in particular the set $\Omega^{N}(\Lambda)$ of N-forms
contains a single element $\omega_{N}=\omega_{1,2,\cdots,N}^{N}{\rm
d}\lambda_{1}\wedge{\rm d}\lambda_{2}\wedge\cdots\wedge{\rm d}\lambda_{N}$.

As recalled at the beginning of this section, if $f$ is a 0-form (a function)
then ${\rm d}f$ is a 1-form. This familiar result generalizes to any
differential forms, with the exterior derivative on forms ${\rm d}$ that
constitutes a map $\Omega^{r}(\Lambda)\overset{{\rm
d}}{\rightarrow}\Omega^{r+1}(\Lambda)$; namely, the (exterior) derivative of
an r-form is an (r+1)-form. In practice it is obtained as

$\displaystyle{\rm
d}\omega_{r}\equiv{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\frac{\partial}{\partial\lambda_{r+1}}}\left(\omega_{j_{1},j_{2},\cdots,j_{r}}^{r}\right){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
d}\lambda_{r+1}\wedge}{\rm d}\lambda_{j_{1}}\wedge{\rm
d}\lambda_{j_{2}}\wedge\cdots\wedge{\rm d}\lambda_{j_{r}}\ .$ (179)

Such an (r+1)-form, $\omega^{r+1}={\rm d}\omega^{r}$ is called an exact form,
as it is derived from an r-form (such as ${\rm d}f$). It may also happen that
a form $\omega$ satisfies ${\rm d}\omega=0$. This is called a closed form.
Importantly, the anti-symmetric relation (178) imposes that ${\rm d}({\rm
d}\omega_{r})=0$, implying that any exact form is closed. The reciprocal is
not true.

More generally, differential forms, such as the 2-form Berry curvature, are
objects that one can integrate over manifolds. In particular, if the dimension
of a manifold $M$ is $r$, then $\int_{M}\omega_{r}$ is a number. Of particular
importance is the volume form $\Omega_{S^{n}}$, that is an n-form whose
integral over a sphere $S^{n}$ embedded in $\mathbb{R}^{n+1}$ is $1$. In
cartesian coordinates, it reads

$\displaystyle\Omega_{S^{n}}\equiv\sum_{i=1}^{n+1}(-1)^{n+1}\frac{x_{i}{\rm
d}x_{1}\wedge\cdots\wedge{\rm d}x_{i-1}\wedge{\rm d}x_{i+1}\wedge\cdots{\rm
d}x_{n+1}}{\gamma_{n}\left(x_{1}^{2}+x_{2}^{2}\cdots+x_{n+1}^{2}\right)^{(n+1)/2}}$
(180)

where $\gamma_{n}=2\pi^{\frac{n+1}{2}}/\Gamma(\frac{n+1}{2})$ with $\Gamma(x)$
the Euler function, is the surface of the unit sphere $S^{n}$. Of particular
importance for the following are the two examples in $\mathbb{R}^{2}$ and
$\mathbb{R}^{3}$

$\displaystyle\Omega_{S^{1}}=\frac{x{\rm d}y-y{\rm
d}x}{2\pi\left(x^{2}+y^{2}\right)}\quad\text{and}\quad\Omega_{S^{2}}=\frac{x{\rm
d}y\wedge{\rm d}z+y{\rm d}z\wedge{\rm d}y+z{\rm d}x\wedge{\rm
d}y}{4\pi\left(x^{2}+y^{2}+z^{2}\right)^{3/2}}\ .$ (181)

A key result in differential calculus is Stokes theorem, that relates the
integral of a differential form $\omega$ with that of its derivative ${\rm
d}\omega$, when it exists everywhere over $M$, as

$\displaystyle\text{Stokes theorem}~{}:\qquad\int_{M}{\rm
d}\omega=\int_{\partial M}\omega$ (182)

where $\partial M$ denotes the boundary of the manifold $M$.

The expression (40) of the Chern number is very general. In particular, it
does not depend on the model at hand. In order to derive an explicit
expression for the Berry-Chern monopole model, it is useful to introduce the
Brouwer theorem, that relates the integrals of a differential form over two
different manifolds
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A}$ and
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}B}$ of same
dimension $n$. Consider a differential form $\omega$ defined over the manifold
$B$ of local coordinates $\mathbf{y}=(y_{1},\cdots,y_{n})$, i.e.
$\omega=\omega_{j_{1},\cdots,j_{n}}(\mathbf{y}){\rm
d}y_{j_{1}}\wedge\cdots\wedge{\rm d}y_{j_{n}}$, and a map $f$ between the two
manifolds
$\mathbf{x}:{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A}\rightarrow\mathbf{y}=f(\mathbf{x})\in{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}B}$.
Then one can define a form on
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A}$ from
$\omega$, through the action of $f$. It is called the pull-back of $\omega$
and reads

$\displaystyle
f^{\star}\omega=\omega_{j_{1},\cdots,j_{n}}(f(\mathbf{y}))\operatorname{det}\left(\frac{\partial
y^{\alpha}}{\partial x^{\beta}}\right){\rm d}x_{j_{1}}\wedge\cdots\wedge{\rm
d}x_{j_{n}}$ (183)

where the Jacobian matrix accounts for the change of local coordinates.

Then one has the following

$\displaystyle\text{Brouwer
theorem:}\qquad\int_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A}}f^{\star}\omega=\deg
f\int_{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}B}}\omega$
(184)

where the degree of $f$ was defined in (79). An important application of the
Brouwer theorem is that it gives an integral formulation of the degree of a
map, which one obtained when applying the formula (184) when
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}B=S^{n}}$
for the volume form $\Omega_{S^{n}}$.

For instance, in the case of the punctured plane
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}B=\mathbb{R}^{2}\backslash\\{0\\}\cong
S^{1}}$, one gets

$\displaystyle\deg h$
$\displaystyle=\int_{{\boldsymbol{\lambda}}\in{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{1}}}h^{\star}\Omega_{S^{1}}$
(185)
$\displaystyle=\int_{{\boldsymbol{\lambda}}\in{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{1}}}\frac{1}{2\pi
h^{2}}\left(h_{x}{\rm d}h_{y}-h_{y}{\rm d}h_{x}\right)\ .$ (186)

Interpreting $\mathbf{h}=(h_{x},h_{x})=(x(t),y(t))$ as a vector position that
depends on a parameter
$t\in{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{1}}$,
one finds the standard expression of the winding number.

The degree is of main interest to characterize the topology of Berry monopole
Hamiltonians i.e. that read
$H=\boldsymbol{h}({\boldsymbol{\lambda}}).{\hat{\mathbf{S}}}$, since the
vector $\boldsymbol{h}$ defines a map
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}\overset{\boldsymbol{h}/|\boldsymbol{h}|}{\rightarrow}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}S^{2}}$
whose degree is an integer that reads

$\displaystyle\deg h$
$\displaystyle=\int_{{\boldsymbol{\lambda}}\in{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}}h^{\star}\Omega_{S^{2}}$
(187)
$\displaystyle=\int_{{\boldsymbol{\lambda}}\in{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}}\frac{1}{4\pi
h^{3}}\left(h_{x}{\rm d}h_{y}\wedge{\rm d}h_{z}+h_{y}{\rm d}h_{z}\wedge{\rm
d}h_{x}+h_{z}{\rm d}h_{y}\wedge{\rm d}h_{x}\right)$ (188)

where the $h_{i}$’s are functions of $\lambda_{\alpha}$ so that ${\rm
d}h_{i}=\frac{\partial h_{i}}{\partial\lambda_{\alpha}}{\rm
d}\lambda_{\alpha}$, and one gets

$\displaystyle\deg
h=\frac{1}{4\pi}\int_{{\boldsymbol{\lambda}}\in{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}}\frac{\epsilon^{ijk}}{h^{3}}h_{i}\frac{\partial
h_{j}}{\partial\lambda_{\alpha}}\frac{\partial
h_{k}}{\partial\lambda_{\beta}}\,{\rm d}\lambda_{\alpha}\wedge{\rm
d}\lambda_{\beta}$ (189)

Finally, using the usual vectorial notation in $\mathbb{R}^{3}$, this
expression takes the form

$\displaystyle\deg
h=\frac{1}{4\pi}\sum_{\alpha,\beta}\int_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}}\frac{\mathbf{h}}{h^{3}}\cdot\left(\frac{\partial\mathbf{h}}{\partial\lambda_{\alpha}}\times\frac{\partial\mathbf{h}}{\partial\lambda_{\beta}}\right)\,{\rm
d}\lambda_{\alpha}\wedge{\rm d}\lambda_{\beta}$ (190)

This formula has a direct meaningful geometrical interpretation : it tells
that the degree of the maps
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}\overset{h}{\rightarrow}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}S^{2}}$
is a wrapping number, meaning that it counts the number of times the vector
$\mathbf{h}/h=\mathbf{n}$ wraps the target sphere
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}S^{2}}$
when ${\boldsymbol{\lambda}}$ spans the entire base space
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}S^{2}}$.

##  Appendix C The equatorial shallow water model and its analytical indices
from the $h.\hat{S}$ model

Here we come back on the celebrated shallow water model that describes
equatorial fluids waves in the Earth’s atmosphere and oceans. Close to the
equator, the modes of frequency $\tilde{\omega}$ and azimutal wave number
$k_{x}$ (along the equator), around a state of rest, are given by the
following equation [85, 80]

$\displaystyle\tilde{\omega}\Psi=\mathcal{H}_{SW}\Psi$ (191)

with

$\displaystyle\mathcal{H}_{SW}=\begin{pmatrix}0&-{\rm i}f(y)&ck_{x}\\\ {\rm
i}f(y)&0&c{\rm i}\partial_{y}\\\ ck_{x}&c{\rm i}\partial_{y}&0\end{pmatrix}$
(192)

in the basis of the perturbed velocity fields in $x$ and $y$ direction and the
perturbed height elevation $(\delta u,\delta v,\delta h)$, where $c$ is the
celerity of gravity waves and $f(y)$ is called the Coriolis parameter: it
accounts for the Coriolis force that depends on the latitude. It has the
dimension of a frequency (it is proportional to Earth angular rotation
$\Omega$) and changes sign at the equator. It is convenient to work in the
tangent plane picture, where the latitude dependence of the Coriolis parameter
reads in terms of the $y$ coordinate that points toward the North pole in a
tangent plane to Earth. In the vicinity of the equator, the Coriolis parameter
can be approximated by $f(y)\sim 2\Omega y/R\equiv\beta y$ where $R$ is the
Earth radius. This simplification is called the $\beta$-plane approximation.
The operator Hamiltonian $\mathcal{H}_{SW}$ is known to display a spectral
flow $\Delta\mathcal{N}_{\pm}=\pm 2$ and $\Delta\mathcal{N}_{0}=0$ for the
three branches $(-1,0,1)$, and the Chern numbers of its symbol symbol are
accordingly $\mathcal{C}_{\pm}=\mp 2$ and $\mathcal{C}_{0}=0$ [15]. Here we
make explicit its form as a particular case of the $h.\hat{S}$ model for $d=1$
and $S=1$ and thus gets its analytical indices.

For that purpose, let us notice that $\mathcal{H}_{SW}$ decomposes as

$\displaystyle\mathcal{H}_{SW}=c{\rm
i}\partial_{y}\hat{S}_{1}+f(y)\hat{S}_{2}+ck_{x}\hat{S}_{3}$ (193)

where the $SU(2)$ commutation relations $[\hat{S}_{1},\hat{S}_{2}]={\rm
i}\hat{S}_{3}$, $[\hat{S}_{2},\hat{S}_{3}]={\rm i}\hat{S}_{1}$ and
$[\hat{S}_{3},\hat{S}_{1}]={\rm i}\hat{S}_{2}$. From this point, it is clear
that the shallow water model falls in the category of $h.\hat{S}$ models with
$S=1$ and $d=1$, and thus its symbol Hamiltonian has the Chern numbers
$\mathcal{C}_{m}=-2m$ with $m=-1,0,1$. To construct the index of
$\mathcal{H}_{SW}$, one needs to transform it into
${\hat{\mathcal{H}}}_{\text{op}}$,which can be done by a unitary
transformation

$\displaystyle U\mathcal{H}_{SW}U^{\dagger}=c{\rm
i}\partial_{y}{\hat{S}_{x}}+f(y){\hat{S}_{y}}+ck_{x}{\hat{S}_{z}}$ (194)

where ${\hat{S}_{z}}$ is diagonal. We find such a unitary operator to be

$\displaystyle U=\frac{1}{\sqrt{2}}\begin{pmatrix}1&0&1\\\ 0&\sqrt{2}&0\\\
-1&0&1\end{pmatrix}$ (195)

and one gets

$\displaystyle U\mathcal{H}_{SW}U^{\dagger}=\begin{pmatrix}k_{x}c&-{\rm
i}\frac{1}{\sqrt{2}}\frac{2\Omega}{R}y-c\partial_{y}&0\\\ {\rm
i}\frac{1}{\sqrt{2}}\frac{2\Omega}{R}y-c\partial_{y}&0&-{\rm
i}\frac{1}{\sqrt{2}}\frac{2\Omega}{R}y-c\partial_{y}\\\ 0&{\rm
i}\frac{1}{\sqrt{2}}\frac{2\Omega}{R}y-c\partial_{y}&-k_{x}c\end{pmatrix}\ .$
(196)

This rotated shallow water Hamiltonian has the same spectrum as the original
one, owing to the unitarity of $U$.

There are two natural length scales in this problem : the Earth radius $R$ and
the Rossby deformation radius $L\equiv\frac{c}{2\Omega}$ which gives the scale
above which rotation is relevant. From this two lengths, one can construct the
parameter $\alpha\equiv=L/R$, which is small in our approach, and the
characteristic length $\ell\equiv\sqrt{RL}$ that allows us to adimensionalize
the shallow water Hamiltonian by introducing
${\hat{a}}=\frac{1}{\sqrt{2}}(\frac{y}{\ell}+\ell\partial_{y})$ and
${\hat{a}^{\dagger}}=\frac{1}{\sqrt{2}}(\frac{y}{\ell}-\ell\partial_{y})$ so
that

$\displaystyle
U\mathcal{H}_{SW}U^{\dagger}=2\Omega\begin{pmatrix}\lambda&-{\rm
i}\alpha{\hat{a}^{\dagger}}&0\\\ {\rm i}\alpha{\hat{a}}&0&-{\rm
i}\alpha{\hat{a}^{\dagger}}\\\ 0&{\rm i}\alpha{\hat{a}}&-\lambda\end{pmatrix}$
(197)

with $\lambda=\frac{k_{x}c}{2\Omega}$, which also reads

$\displaystyle\mathcal{H}_{SW}=2\Omega U^{\dagger}\left(-{\rm
i}\alpha{\hat{a}^{\dagger}}{\hat{S}_{+}}+{\rm
i}\alpha{\hat{a}}{\hat{S}_{-}}+\lambda{\hat{S}_{z}}\right)U\ .$ (198)

From here, we can introduce the indices
$\text{ind}(\mathcal{D}_{SW}^{2m_{S}})$ of the shallow water model with
$\mathcal{D}_{SW}={\rm i}\alpha U^{\dagger}{\hat{a}}{\hat{S}_{-}}U$ Since the
index is invariant by homotopy, one gets

$\displaystyle\text{ind}(\mathcal{D}_{SW}^{2m_{S}})=\text{ind}(U^{\dagger}{\hat{a}}{\hat{S}_{-}}U)^{2m_{S}}$
(199)

and since it is also invariant by a unitary transformation, one finally finds

$\displaystyle\text{ind}(\mathcal{D}_{SW}^{2m_{S}})=\text{ind}(\mathcal{D}^{2m_{S}})=2m_{S}\
.$ (200)

with $m_{S}=0,\pm 1$. This is precisely the index of the $h.\hat{S}$ model for
$d=1$ and $S=1$, which yields the spectral flows

$\displaystyle\Delta\mathcal{N}_{\pm}=\pm
2\quad\text{and}\quad\Delta\mathcal{N}_{0}=0$ (201)

as expected.

We now express the corresponding eigenmodes of the spectral flow. Following
the procedure detailed in the section 3.4, one can construct the two
eigenstates for the spectral flow in the case $S=1$ and $d=1$:
$\ket{\Psi_{1}}$ with $1$ zero component and $\ket{\Psi_{2}}$ with $2$ zero
components in the $\ket{m_{S}}$ basis, as

$\displaystyle\ket{\Psi_{2}}=\begin{pmatrix}\ket{1}\\\ 0\\\
0\end{pmatrix}\qquad\ket{\Psi_{1}}=\begin{pmatrix}\psi_{1}\ket{1}\\\
\ket{0}\\\ 0\end{pmatrix}$ (202)

where $\psi_{1}$ is a coefficient. The profiles in the $y$ direction of the
eigenstates for the spectral flow of the original shallow water model in the
$\beta$-plane have therefore the following form

$\displaystyle\braket{y}{U^{\dagger}\Psi_{2}}=\begin{pmatrix}1\\\ 0\\\
-1\end{pmatrix}{\rm e}^{-\frac{1}{2}(\frac{y}{\ell})^{2}}\equiv\Psi_{K}(y)$
(203) $\displaystyle\braket{y}{U^{\dagger}\Psi_{1}}=\begin{pmatrix}y\\\
\psi_{1}\\\ -y\end{pmatrix}{\rm
e}^{-\frac{1}{2}(\frac{y}{\ell})^{2}}\equiv\Psi_{Y}(y)\ .$ (204)

The mode $\Psi_{K}$ corresponds to a longitudinal velocity field $\delta u$
and a height elevation $\delta h$ which are symmetric with respect to the
latitude $y$, while the velocity perpendicular to the equator $\delta v$ is
zero. It corresponds to the well-known equatorial Kelvin mode. The second
mode, $\Psi_{Y}$, has a longitudinal velocity field $\delta u$ and a height
elevation $\delta h$ which are antisymmetric with respect to the latitude $y$,
and the velocity perpendicular to the equator $\delta v$ is nonzero and
symmetric. It is known as the Yanai or mixed planetary-gravity mode [85, 80].


# TeleViT: Teleconnection-driven Transformers Improve Subseasonal to Seasonal
Wildfire Forecasting

Ioannis Prapas1, 3 Nikolaos Ioannis Bountos1,2 Spyros Kondylatos1, 3
Dimitrios Michail2 Gustau Camps-Valls3 Ioannis Papoutsis1
1Orion Lab, Institute for Astronomy, Astrophysics, Space Applications,
and Remote Sensing (IAASARS), National Observatory of Athens
2Department of Informatics and Telematics, Harokopio University of Athens
3Image Processing Laboratory (IPL), Universitat de València

###### Abstract

Wildfires are increasingly exacerbated as a result of climate change,
necessitating advanced proactive measures for effective mitigation. It is
important to forecast wildfires weeks and months in advance to plan forest
fuel management, resource procurement and allocation. To achieve such accurate
long-term forecasts at a global scale, it is crucial to employ models that
account for the Earth system’s inherent spatio-temporal interactions, such as
memory effects and teleconnections. We propose a teleconnection-driven vision
transformer (TeleViT), capable of treating the Earth as one interconnected
system, integrating fine-grained local-scale inputs with global-scale inputs,
such as climate indices and coarse-grained global variables. Through
comprehensive experimentation, we demonstrate the superiority of TeleViT in
accurately predicting global burned area patterns for various forecasting
windows, up to four months in advance. The gain is especially pronounced in
larger forecasting windows, demonstrating the improved ability of deep
learning models that exploit teleconnections to capture Earth system dynamics.
Code available at github.com/Orion-Ai-Lab/TeleViT.

##  1 Introduction

Global warming increases the frequency and intensity of fire weather,
amplifying the likelihood of the conditions that lead to extreme wildfire
events [10]. In that context, it is important to improve our understanding of
wildfires, anticipating wildfire patterns weeks and months in advance.
Operationally, subseasonal to seasonal wildfire forecasts are typically
treated as anticipated anomalies in temperature and precipitation, derived
from process-based climate models [20], while disregarding crucial fire-
related factors like soil moisture and vegetation dynamics. Deep Learning (DL)
methods, able to learn from data, offer a promising avenue to model Earth
system processes such as wildfires [17], arising from the dynamic interactions
of all the different fire drivers, namely climate, vegetation and human
activity [8].

### Previous work

In fact, several studies have successfully applied DL to wildfire forecasting
tasks [9, 7]. When it comes to global predictions on subseasonal to seasonal
scales, existing work relies on traditional Machine Learning (ML) approaches
[16] that do not effectively capture Earth system dynamics that are important
for long-term forecasting, even if sometimes they use teleconnection indices
as input features [14, 21]. In this work, we argue that for predictions on
long temporal scales, it is crucial to train models that consider the Earth as
one interconnected system, accounting for spatio-temporal interactions such as
memory effects and teleconnections. There is substantial evidence that
teleconnections modulate global wildfires [11, 12, 5, 2]. For example, extreme
wildfires in Siberia have been linked to preceding arctic oscillation [12] and
previous-year soil moisture anomalies [5]. Despite the evidence, there is very
limited ML work to predict wildfires that combines local information from the
fire drivers with Earth system variables. Chen _et al_. [3] use simple
autoregressive statistical methods to combine volume pressure deficit values
with oceanic indicators. Yu _et al_. [21] use a statistical pre-processing to
identify the most prominent oceanic indicators modulating burned area in
Africa and then use the findings to select input features for tree-based
machine learning models. AttentionFire [14] uses an attention-enhanced
recurrent neural network architecture that considers temporal context and
shows slight improvements at longer horizons when adding information from
oceanic indices. More sophisticated models that handle the Earth as a system
are arising naturally for weather and climate prediction. Earthformer [6]
combines spatio-temporal attention with learnable global vectors that are
meant to summarize the dynamics of a system. GraphCast [13] uses graph neural
networks that capture long-range interactions by defining graphs at several
resolutions. Climax [15] is proposed as a transformer-based foundational model
for weather and climate at coarse resolutions, working with global input or on
a very large spatial domain.

### Application context

Our study has direct implications for humanitarian assistance and disaster
response (HADR) operations, particularly in addressing the challenges of
intensified wildfires due to climate change. Forecasting wildfires weeks or
months in advance can inform the deployment of HADR-related anticipatory
actions. This includes measures for vegetation and forest management (_e.g_.
controlled burns and targeted firebreaks), coordination with local authorities
to form evacuation plans, procurement of resources such as leasing
firefighting equipment, and timely mobilisation of international and cross-
border cooperation and aid.

Our transformer-based models combine local and non-local inputs like
teleconnection indices and coarsened global variables, and thus offer a novel
approach to model Earth system processes and improve our capabilities to
forecast wildfires at subseasonal to seasonal scales. Outside of wildfire
forecasting, these methods hold promise for diverse HADR applications,
enabling foresight across extended timeframes, which can empower proactive
planning and optimized resource allocation before disastrous events.

### Contributions

In this work, we propose to explicitly model short (local) and long (global)
range interactions of the Earth system for long-term global wildfire pattern
forecasting. For this, we develop a Teleconnection-driven Vision Transformer
(TeleViT), that expands ViT [4] with an asymmetric tokenisation procedure in
order to seamlessly combine local and global scale inputs, _i.e_. climate
indices and coarsened global variables. We thoroughly examine the performance
of the proposed model, demonstrating superior performance in various
forecasting windows compared to models that do not leverage teleconnections.

##  2 Teleconnection Vision Transformer (TeleViT)

To produce reliable long-range forecasts (weeks or even months in advance), it
is imperative to treat the Earth as a connected system. The Earth system is
characterised by continuous interactions of processes that spread over large
spatio-temporal windows, manifesting as memory effects and teleconnections.
Memory effects are the realisation of the persistence or influence of past
events on current and future states. They represent the temporal aspect of
system behaviour, where the history of a system influences its present or
future behaviour. In the context of wildfire prediction, memory effects
capture how past events such as fuel accumulation, drought conditions, and
weather patterns can impact the likelihood, behaviour, and extent of
wildfires. Teleconnections constitute long-distance interactions between
different regions in the Earth system. They describe how changes in one region
can influence atmospheric or oceanic conditions in another, often through
large-scale atmospheric circulation patterns or oceanic phenomena [19].
Teleconnections are often described in terms of teleconnection indices, _i.e_.
Oceanic and Climatic Indices (OCIs), calculated as large-scale anomalies of
specific parameters, such as temperature, pressure or sea surface temperature.
By teleconnections, however, we refer to the long-range spatio-temporal
interactions and not to the indices per se, which we consider a mere proxy to
the state of the Earth system.

Figure 1: Full pipeline of the TeleViT architecture. The different multi-scale
inputs _i.e_. local, global and teleconnection indices, are tokenized at
different resolutions and fed to a Transformer encoder along with a prepended
classification token. The linear decoder is based on the output of the
classification token.

We propose to capture such distant interactions that are omnipresent in the
Earth system, using data that can inform on the state and dynamics of the
system. Particularly, we propose to enhance fine-resolution local data, _i.e_.
information for a small area on Earth, that are commonly used in isolation,
with information from i) OCIs and ii) coarse-resolution global data.
Effectively combining these inputs is particularly challenging due to high
dimensionality and the discrepancy between the spatial resolution of the
different data sources.

In that direction, we propose to utilise the versatility of the Transformer
architecture[18]. We build upon the Vision Transformer (ViT) [4], which adapts
the architecture to computer vision problems by splitting an image into non-
overlapping chips (tokens). The resulting token sequence serves as input to a
standard Transformer. Extending ViT, we introduce the Teleconnection-driven
Vision Transformer (TeleViT). TeleViT relies on an asymmetric tokenization
method which increments the sequence of input tokens with tokens stemming from
different data sources, with potentially varying spatial and temporal
resolutions. To produce the token sequence, each input source is tokenised
independently, taking into account its inherent characteristics along with the
scale it operates on the Earth system. Local information and time series of
climatic indices should be represented in high detail with smaller token
sizes, while global information, operating at a greater scale, can be less
detailed benefitting from larger token sizes.

TeleViT can be adapted to consume any type of input, given a tokenisation
strategy. For our setting, we assume local input $x_{l}\in R^{C_{l}\times
H_{l}\times W_{l}}$, global input $x_{g}\in R^{C_{g}\times H_{g}\times W_{g}}$
and indices input $x_{i}\in R^{C_{i}\times T}$, with C the number of features,
H the height, W the width and T the time-series length. According to our
asymmetric tokenisation, each of the inputs is split into tokens individually,
resulting in $N_{l}$, $N_{g}$ and $N_{i}$ tokens with dimensions $P_{l}$,
$P_{g}$ and $P_{i}$ that depend on the selected tokenisation strategy. Tokens
are mapped to the embedding dimension $D$ of the Transformer with a different
trainable linear projection function for each input $f_{k}:R^{N_{k}\times
P_{k}}\rightarrow R^{N_{k}\times D},k\in{l,g,i}$. Similar to ViT, we add
learnable positional encoding to the sequence of $N_{l}+N_{g}+N_{i}$ number of
tokens, which is fed to a standard Transformer with depth $K$ and $A$
attention heads for each layer. Pixel-level prediction is performed by a
simple decoder attached to the prepended classification token (cls). We choose
a trainable linear layer that maps to the input resolution. The architecture
is demonstrated in Figure 1.

As the Transformer architecture remains intact, most of the inductive biases
come from the tokenisation procedure. By keeping the inductive bias to a
minimum we remove any restrictions to known concepts, enabling unprecedented
information combination of different input datasets, representing local and
global inputs. For example, attention allows i) inter-dataset interactions,
_i.e_. information flows from distant large-scale regions of coarse resolution
to high-resolution localised windows, and ii) intra-dataset interactions,
where several Earth system processes and interactions are modelled at a global
scale while identifying location-specific characteristics for each region
independently (depicted in Figure 4, Appendix D). The simplicity of the method
makes it easily extensible to other data sources, as well as to the inclusion
of a temporal dimension that we do not address in this study.

##  3 Experiments

We conduct our experiments on the SeasFire cube [1], a spatio-temporal dataset
for subseasonal to seasonal wildfire forecasting. It contains 21 years of data
(2001-2021) at a global scale, in an 8-days temporal and $0.25^{\circ}$
spatial resolution. The cube includes a diverse range of seasonal fire
drivers, combining atmospheric, vegetation, and anthropogenic variables along
with climate indices, in addition to target variables related to wildfires
such as burned areas, fire radiative power, and wildfire-related emissions.

Building on the SeasFire cube, Prapas _et al_. [16] defined burned area
pattern forecasting as a segmentation task. As input, they use local patches
that contain different channels of the fire driver variables and train a
U-Net++[22] to predict the presence of burned areas at a future time step.
Their U-Net++ demonstrated a predictive skill greater than the burned area
mean seasonal cycle for a lead forecasting time of up to $2$ months. We follow
a similar setup, where given a snapshot of the fire driver variables at
timestep $t$ we want to predict the presence of burned areas at a future
timestep $t+h$. For the local input, we extract patches of size $80\times 80$,
and as such, the world represented at $0.25^{\circ}$ ($1440\times 720$ cells),
is split into $18\times 9=162$ local input patches. To extract the global
input, we coarsen the cube to $1^{\circ}$ (see Appendix C), reducing its size
by a factor of 16, making the global input size 360$\times$180\. Along with 10
fire driver variables extracted from the cube, the same for both local and
global inputs, we calculate a global positional encoding, _i.e_. sine and
cosine of the longitude and latitude. For each sample, 10 OCIs are extracted
for the 10 months preceding $t$. The target is the presence of burned area at
time-step $t+h$ for the region of the local input, where h is the lead time
forecasting horizon. As such, a sample is comprised of four vectors; i) a
local input $x_{l}$ of size $(14,80,80)$, ii) a global input $x_{g}$ of size
$(14,360,180)$, iii) an OCI input $x_{i}$ of size $(10,10)$ and iv) a target
vector of size $(1,80,80)$. The variables used are shown in Appendix A. The
model’s performance is evaluated using the Area Under the Precision-Recall
Curve (AUPRC). The train, validation, and test split is time-based, using
years 2002 - 2017 for training, 2018 for validation and 2019 for testing. Only
samples that contain burned areas are considered.

We assess the following models: i) a U-Net++, which uses only local input
$x_{l}$ [16], ii) a simple ViT, which uses only local input $x_{l}$, iii)
$\mathrm{TeleViT}_{i}$, which uses OCIs $x_{i}$ along with local input
$x_{l}$, iv) $\mathrm{TeleViT}_{g}$, which uses only global input $x_{g}$
along with local input $x_{l}$, and v) $\mathrm{TeleViT}_{i,g}$, which uses
both OCIs $x_{i}$ and global input $x_{g}$ along with the local input $x_{l}$.
The hyper-parameters have been tuned for simple ViT and applied to TeleViT
models that use the same core architecture. We provide detailed information on
architectural and training choices in Appendix B.

The performance of the five models is examined in several forecasting horizons
$h\in\\{1,2,4,8,16\\}$. A model that predicts at the maximum forecasting
horizon, _i.e_. 16$\times$8-days in advance, learns to predict the burned area
pattern of a particular 8-day period approximately four months in advance. A
different model is trained for each $h$, for a total of $25$ experiments.

##  4 Results and Discussion

Figure 2 summarises the results of the experiments. In general, all the models
exhibit a declining trend in performance as the forecasting window increases.
The decline, however, is much less steep for teleconnection-driven models,
than for the baselines (ViT and U-Net++). ViT proves to be a stronger
baseline, outperforming U-Net++ with the exception of the 16$\times$8-days,
where U-Net++ overtakes it by a small margin. $\mathrm{TeleViT}_{g}$ achieves
a comparable performance to ViT for short forecasting windows of up to
4$\times$8-days, while it shows greater robustness to the increase of the
forecasting window. $\mathrm{TeleViT}_{i}$ and $\mathrm{TeleViT}_{i,g}$
consistently surpass the baselines, with larger performance gaps as the lead
time increases. Notably, $\mathrm{TeleViT}_{i,g}$ has the dominant
performance, which indicates the benefits of a synergistic effect between
teleconnection indices and global-scale representations. Interestingly, it
achieves high gains even in short forecasting windows, which suggests that
contextual information brought by this combination is helpful beyond
teleconnections that operate on larger temporal scales. Figure 3 shows how the
predictions compare to the target for a specific sample date, demonstrating
high agreement.

Figure 2: AUPRC performance of the different models for forecasting windows of
1, 2, 4, 8 and 16$\times$8-days in advance.

In general, the results demonstrate TeleViT’s ability to successfully fuse
local and non-local Earth system information and pave the road for exciting
future research. Both global-scale inputs and OCIs can bring performance
gains, but it is not clear if this is due to the same reasons. Further work is
needed to investigate the relative contribution of each input. This urges us
to comprehend the models, elucidating both familiar and undiscovered
interactions between teleconnections and their influence on burned area
patterns. A thorough examination of the attention maps may offer valuable
insights in this regard, potentially shedding light on the underlying
mechanisms. There is also much potential for further investigation of
tokenisation schemes. The performance improvement induced by the introduction
of coarsened global views may suggest that OCIs are a simplified proxy of the
Earth system state that could be enhanced by global state representations.
Further investigation will reveal the contexts where this enhancement is most
pronounced. Future work can exploit time series for both local and global
inputs. In fact, it can be tested if global inputs with a temporal component
can replace the indices. Finally, using models that treat the Earth as a
system holds immense promise for scientific knowledge discovery, especially
when incorporating knowledge of the physical systems.

Figure 3: A sample prediction versus the target for the best model, predicting
at 4$\times$8-days lead forecasting time. Sea and values lower than 0.05 are
masked out. Confidence is determined as the softmax score of the positive
prediction.

##  5 Conclusions

In this work, we propose TeleViT, a transformer-based architecture that can
leverage local and non-local information to model global Earth system
processes. Using the SeasFire cube, we showcase our model’s ability to improve
burned area pattern forecasting using climatic indices and coarsened global
views of fire driver variables. This improvement is intensified for longer
temporal forecasting horizons, where teleconnections are expected to have a
more significant impact, strengthening our hypothesis that TeleViT can
discover such long spatio-temporal Earth system interactions.

The combination of local information, climatic indices, and global views of
the Earth holds significant potential for various HADR applications, promising
enhanced anticipation capabilities at longer temporal horizons. This includes
forecasting extreme events such as floods or droughts, the intensity of
tropical cyclones, impacts on food security or human displacement. These
applications can benefit from the enhanced prognosis, potentially improving
proactive strategies and resource allocation in the face of disasters.

## Acknowledgements

This work is part of the SeasFire project, which deals with ”Earth System Deep
Learning for Seasonal Fire Forecasting” and is funded by the European Space
Agency (ESA) in the context of the ESA Future EO-1 Science for Society Call.

## References

  * [1] Lazaro Alonso, Fabian Gans, Ilektra Karasante, Akanksha Ahuja, Ioannis Prapas, Spyros Kondylatos, Ioannis Papoutsis, Eleannna Panagiotou, Dimitrios Mihail, Felix Cremer, Ulrich Weber, and Nuno Carvalhais.  SeasFire cube: A global dataset for seasonal fire modeling in the earth system.
  * [2] Adrián Cardil, Marcos Rodrigues, Mario Tapia, Renaud Barbero, Joaquin Ramírez, Cathelijne R. Stoof, Carlos Alberto Silva, Midhun Mohan, and Sergio de Miguel.  Climate teleconnections modulate global burned area.  14(1):427.  Number: 1 Publisher: Nature Publishing Group.
  * [3] Yang Chen, James T. Randerson, Shane R. Coffield, Efi Foufoula‐Georgiou, Padhraic Smyth, Casey A. Graff, Douglas C. Morton, Niels Andela, Guido R. Werf, Louis Giglio, and Lesley E. Ott.  Forecasting global fire emissions on subseasonal to seasonal (s2s) time scales.  12(9).
  * [4] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.
  * [5] Matthias Forkel, Kirsten Thonicke, Christian Beer, Wolfgang Cramer, Sergey Bartalev, and Christiane Schmullius.  Extreme fire events are related to previous-year surface moisture conditions in permafrost-underlain larch forests of siberia.  7(4):044021.  Publisher: IOP Publishing.
  * [6] Zhihan Gao, Xingjian Shi, Hao Wang, Yi Zhu, Yuyang Wang, Mu Li, and Dit-Yan Yeung.  Earthformer: Exploring space-time transformers for earth system forecasting.
  * [7] Rafik Ghali and Moulay A. Akhloufi.  Deep learning approaches for wildland fires using satellite remote sensing data: Detection, mapping, and prediction.  6(5):192.  Number: 5 Publisher: Multidisciplinary Digital Publishing Institute.
  * [8] Stijn Hantson, Almut Arneth, Sandy P. Harrison, Doug I. Kelley, I. Colin Prentice, Sam S. Rabin, Sally Archibald, Florent Mouillot, Steve R. Arnold, and Paulo Artaxo.  The status and challenge of global fire modelling.  13(11):3359–3375.  Publisher: Copernicus Publications.
  * [9] Piyush Jain, Sean CP Coogan, Sriram Ganapathi Subramanian, Mark Crowley, Steve Taylor, and Mike D. Flannigan.  A review of machine learning applications in wildfire science and management.
  * [10] Matthew W. Jones, John T. Abatzoglou, Sander Veraverbeke, Niels Andela, Gitta Lasslop, Matthias Forkel, Adam J. P. Smith, Chantelle Burton, Richard A. Betts, Guido R. van der Werf, Stephen Sitch, Josep G. Canadell, Cristina Santín, Crystal Kolden, Stefan H. Doerr, and Corinne Le Quéré.  Global and regional trends and drivers of fire under climate change.  60(3):e2020RG000726.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2020RG000726.
  * [11] Flavio Justino, David H. Bromwich, Vanucia Schumacher, Alex daSilva, and Sheng-Hung Wang.  Arctic oscillation and pacific-north american pattern dominated-modulation of fire danger and wildfire occurrence.  5(1):1–13.  Number: 1 Publisher: Nature Publishing Group.
  * [12] Jin-Soo Kim, Jong-Seong Kug, Su-Jong Jeong, Hotaek Park, and Gabriela Schaepman-Strub.  Extensive fires in southeastern siberian permafrost linked to preceding arctic oscillation.  6(2):eaax3308.  Publisher: American Association for the Advancement of Science.
  * [13] Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Alexander Pritzel, Suman Ravuri, Timo Ewalds, Ferran Alet, Zach Eaton-Rosen, Weihua Hu, Alexander Merose, Stephan Hoyer, George Holland, Jacklynn Stott, Oriol Vinyals, Shakir Mohamed, and Peter Battaglia.  GraphCast: Learning skillful medium-range global weather forecasting.
  * [14] Fa Li, Qing Zhu, William J. Riley, Lei Zhao, Li Xu, Kunxiaojia Yuan, Min Chen, Huayi Wu, Zhipeng Gui, Jianya Gong, and James T. Randerson.  AttentionFire_v1.0: interpretable machine learning fire model for burned-area predictions over tropics.  16(3):869–884.
  * [15] Tung Nguyen, Johannes Brandstetter, Ashish Kapoor, Jayesh K. Gupta, and Aditya Grover.  ClimaX: A foundation model for weather and climate.
  * [16] Ioannis Prapas, Akanksha Ahuja, Spyros Kondylatos, Ilektra Karasante, Eleanna Panagiotou, Lazaro Alonso, Charalampos Davalas, Dimitrios Michail, Nuno Carvalhais, and Ioannis Papoutsis.  Deep learning for global wildfire forecasting.
  * [17] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, and Nuno Carvalhais.  Deep learning and process understanding for data-driven earth system science.  566(7743):195–204.  Publisher: Nature Publishing Group.
  * [18] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, {\textbackslash}Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Advances in neural information processing systems, pages 5998–6008.
  * [19] John M. Wallace and David S. Gutzler.  Teleconnections in the geopotential height field during the northern hemisphere winter.  109(4):784–812.  Publisher: American Meteorological Society Section: Monthly Weather Review.
  * [20] Wikipedia contributors.  Effis - seasonal forecast.  https://effis.jrc.ec.europa.eu/about-effis/technical-background/seasonal-forecast-explained, 2023\.  [Online; accessed 16-June-2023].
  * [21] Yan Yu, Jiafu Mao, Peter E. Thornton, Michael Notaro, Stan D. Wullschleger, Xiaoying Shi, Forrest M. Hoffman, and Yaoping Wang.  Quantifying the drivers and predictability of seasonal changes in african fire.  11(1):2893.  Number: 1 Publisher: Nature Publishing Group.
  * [22] Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang.  UNet++: Redesigning skip connections to exploit multiscale features in image segmentation.

##  Appendix A Input Data from the SeasFire cube

Table 1 shows which variables are used from the SeasFire cube[1], along with
some pre-processing made for each variable. Total precipitation and population
are log-transformed with $log(1+x)$ to follow a less skewed distribution. For
more details on the variables, the reader is referred to the cited dataset.

Full name | Pre-processing
---|---
Local/Global Variables |
Mean sea level pressure |
Total precipitation | Log-transformed
Vapour Pressure Deficit |
Sea Surface Temperature |
Mean Temperature at 2 meters |
Surface solar radiation downwards |
Volumetric soil water level 1 |
Land Surface Temperature at day |
Normalized Difference Vegetation Index |
Population density | Log-transformed
Cosine of longitude | Calculated
Sine of longitude | Calculated
Cosine of latitude | Calculated
Sine of latitude | Calculated
Climatic Indices |
Western Pacific Index |
Pacific North American Index |
North Atlantic Oscillation |
Southern Oscillation Index |
Global Mean Temperature |
Pacific Decadal Oscillation |
Eastern Asia/Western Russia |
East Pacific/North Pacific Oscillation |
Niño 3.4 Anomaly |
Bivariate ENSO Timeseries |
Target Variable |
Burned Areas from GWIS | Made binary
Table 1: Input and target variables used from the SeasFire cube for all
settings. The same variables are used for both local and global views.

##  Appendix B Model Details and Hyperparameters

Models are trained for 50 epochs. We use the cross-entropy loss and the Adam
optimizer to train the models. For the U-Net++ model, the initial learning
rate is set to 0.001, while for the Transformer models, it is set to 0.0001.
The learning rate is reduced on the plateau and the weight decay is set to
0.000001. The model with the lowest validation loss is used for testing.
Before entering the models, local and global inputs are normalized. OCIs,
which are anomalies are divided by their standard deviation.

The encoders of both ViT and TeleViT consist of $K=8$ layers, with $A=12$
attention heads each, and an embedding dimension $D=768$. For the asymmetric
tokenisation, we set $P_{l}=(1,16,16)$, $P_{g}=(1,30,30)$ and $P_{i}=1$. This
means the following:

  * •

The local input $x_{l}$ of size $(14,80,80)$ is tokenised spatially in
$5\times 5$ number of tokens with size $16\times 16$.

  * •

The global input $x_{g}$ of size $(14,360,180)$ is tokenised spatially in
$12\times 6$ number of tokens with size $30\times 30$.

  * •

The OCI input $x_{i}$ of size $(10,10)$ is tokenised in $10\times 10$ number
of tokens.

##  Appendix C Coarsening the SeasFire cube

The SeasFire cube is provided as a xarray-compatible file that is in Zarr
format, which makes it easy to coarsen with xarray. All that is needed is to
provide an aggregation function for each of the variables and a coarsening
factor for each dimension. We use a coarsening factor of 4 along the longitude
and latitude dimensions to convert the cube from $0.25^{\circ}$ to $1^{\circ}$
spatial resolution and all the input variables are mean-aggregated.

##  Appendix D Interactions between the different input datasets

Figure 4 shows different interactions between the different input datasets
captured by attention.

$a_{11}$$a_{12}$$a_{13}$$a_{14}$$a_{15}$$a_{16}$$a_{21}$$a_{22}$$a_{23}$$a_{24}$$a_{25}$$a_{26}$$a_{31}$$a_{32}$$a_{33}$$a_{34}$$a_{35}$$a_{36}$$a_{41}$$a_{42}$$a_{43}$$a_{44}$$a_{45}$$a_{46}$$a_{51}$$a_{52}$$a_{53}$$a_{54}$$a_{55}$$a_{56}$$a_{61}$$a_{62}$$a_{63}$$a_{64}$$a_{65}$$a_{66}$LocalIndicesGlobalLocalIndicesGlobal
Figure 4: Depiction of intra-dataset (blue) and inter-dataset interactions
(orange) in an attention matrix for a sequence of six tokens with an equal
number of tokens for each data source.