text:
\begin{document}
\maketitle

\begin{abstract}
We prove that with a $(2+1)$-dimensional Toda type system are associated algebraic skeletons which are (compatible assemblings)  of particle-like Lie algebras of dyons and triadons type.  We obtain trix-coaxial and dyx-coaxial Lie algebra structures for the system from algebraic skeletons of some particular choice for compatible associated absolute parallelisms.
In particular, by a first choice of the absolute parallelism, we associate with the $(2+1)$-dimensional Toda type system a trix-coaxial Lie algebra structure
made of two (compatible) base triadons constituting a $2$-catena.  Furthermore, by a second choice of the absolute parallelism, we associate a dyx-coaxial Lie algebra structure
made of two (compatible) base dyons, as well as particle-like Lie algebra structures made of single $3$-dyons. Some explicit examples of applications such as conservation laws related to special solutions, and an inverse spectral problem are worked out.

\end{abstract}

\noindent {\bf Key words}: Particle-like Lie algebra structure, infinitesimal skeleton, tower, Toda system.

\section{Introduction}

Toda type systems are nonlinear models which play a role in a variety of physical and, more in general, natural phenomena.

The problem of integrability of nonlinear models has been recognized to be related to their algebraic properties in discrete and continuous, as well as, classical and quantum formulations.
 Algebraic properties can be interpreted as the counterpart of the concept of integrability given as of having `enough' conservation laws to exhaustively describe the underling field or associated dynamics.
Indeed, from an historical point of view, algebraic-geometric approaches are based on the requirement for the existence of conservation laws which emerge from internal symmetries (given in terms of algebraic structures).

In the Seventies, in fact, Wahlquist and Estabrook \cite{WaEs75,Estabrook} proposed a  technique for systematically deriving, from an integrable  system, what they called a `prolongation structure' in terms of a set of `pseudopotentials' related to the existence of an infinite set of associated conservation laws.
They also conjectured that, as a characterizing feature of the integrability property, the structure was `open' \ie not a set of structure relations of a finite--dimensional Lie group. Since then, `open' Lie algebras have been extensively studied in order to distinguish them from freely generated infinite-dimensional Lie algebras.

Their interest in the study of integrability is in the fact that Lax pairs of the inverse spectral transform containing an isospectral parameter can be obtained by an homomorphism of the infinite-dimensional open Lie algebra in a finite-dimensional `closed' Lie algebra.
In their approach, conservation laws  are written in terms of `prolongation' forms and integrability is intended as a Frobenius
integrability condition for a `prolonged' ideal of differential forms describing intrinsically the given nonlinear model in the sense of \'E. Cartan.

Attempting a {\em description of symmetries in terms of Lie algebras} implies the appearance of an homogeneous space and thus the interpretation of prolongation forms as {\em Cartan--Ehresmann connections}.
It is clear that here the unknowns are both conservation laws and symmetries,
and the main point in this is how to realize the form of the conservation laws
and thus the explicit expression of the prolongation forms.
Different  prolongation ideals give rise to both different algebraic structures (symmetries) and corresponding conservation laws.
By an inverse procedure based on the intrinsic duality between Lie algebras and differential systems \cite{Es82}, open Lie algebraic structures can `generate' whole families of different nonlinear systems bound by the same internal symmetry structure.

In a series of papers \cite{Pa05,Pa13,Pa16_Eco,PaWi02,PaWi10,PaWi11}, we explicated
an algebraic-geometric interpretation of the above mentioned `prolongation' procedure
in terms of {\em towers with infinitesimal algebraic skeletons} (in the sense of  \cite{Mo93}) and we will refer to that framework in this paper.
It is noteworthy that slight modification of the internal symmetry properties generates {\em new} models which can contain possible integrable subcases.
For example, activator-substrate systems have been obtained by performing a slight modification of the internal symmetry algebra of twisted reaction-diffusion equations \cite{Pa16_Eco}.

The structure itself with which the tower forms are postulated can produce open algebraic structures or just Lie algebras.
Our aim in this work is to investigate some common features of them and to show the emerging of particle-like Lie algebras structure as symmetry structures of integrable systems (associated with Poisson structures  the compatibility of which is worthy of study \cite{Fernandes93,Kos10}).
Indeed, in general, infinite dimensional open Lie algebras are the main object of the search in view of the application of the inverse spectral transform to obtain soliton solutions, B\"acklund transformations and so on;  recent examples of applications can be found \eg in \cite{IgMa19,IgMa20,Morozov,Wang,Wang2,Li}. Although these features are not our prominent task in this paper, an inverse problem will be obtained in Section \ref{Lax section}.

Prolongation forms bringing to finite dimensional Lie algebras (without a spectral parameter) are generally discarded when searching for a Lax pair to be used within the inverse spectral transform.

However, integrable systems, admitting infinite-dimensional prolongation Lie algebras can also admit finite-dimensional Lie algebras, which still can be related to some kind of internal symmetries of the systems themselves and to associated conservation laws, or even to B\"acklund transformations.
We refer, in particular, to the paper ``More prolongation structures'' by C. Hoenselaers \cite{Hoenselaers}, which pointed out two important features of the algebraic structures obtained by the method of Wahlquist and Estabrook.

First feature: it can be that the prolongation forms can not always be solved in such a way that one obtains commutators among vector fields depending only on the `pseudopotentials' coordinates. A typical example is, in fact, the most general prolongation problem associated by such a procedure to equation \eqref{1}: in \cite{Palese} the prolongation problem was formally solved by introducing suitable operators of Bessel type, however a prolongation algebra (and then an inverse problem) could not be obtained explicitly.

Furthermore, it is noteworthy that a certain arbitrariness is given by postulating the structure of the tower in the search of a skeleton. Say, to a given equation can be associated different towers with different skeletons. As it was already stressed by Estrabrook himself \cite{Es82} the same algebraic structure `contains' families of equations, associated linear spectral problems and B\"acklund transformations, interrelated by transformations between dependent and independent coordinates.

The motivation of our previous research for the interpretation of prolongation structures as skeletons of towers was the aim of more deeply understanding such a feature: {\em equations, linear problems, B\"acklund transformations, are local coordinates expressions of common intrinsic structures}; this is of help also in practical questions: solutions of systems can be obtained by simpler systems having in common (part of) skeletons. In Section \ref{Conservation laws} we show that each one of the two compatible $4$-triadon constituting the skeleton given by a trix-coaxial Lie algebra structure generates the same conservation law and related special solutions.
This justifies the possible choice of a more restricted (instead of the most general one) form of the tower (then of the algebraic skeleton) still getting `solutions' (with this term meaning analytical solutions as well as particular conservation laws) of the original equation.

Second feature: even if the prolongation algebra is a finite dimensional (even abelian in his example) Lie algebra,  nevertheless there can exist B\"acklund transformations. It is shown that the exterior differential associated to the prolongation structure of the NLS equation being of genus $3$, and it is stressed that we can choose dependent and independent variables in an arbitray way.
We can also lower the genus (in our skeleton formulation this means the choice of different representations $\rho$ or even different vector spaces $V$) so obtaining a reduced ideal where one of the independent coordinates is turned into a dependent coordinate. This turning a global symmetry in a local one provides Miura type transformations between the modified NLS and another system, the prolongation structure of which is finite dimensional and there is only one nontrivial potential entering a B\"acklund transformation acting on the modified NLS equation.

In few words the Wahlquist-Estabrook method not always produces infinite dimensional open Lie algebras, but it could be that by that `procedure' we get  only a part of an algebraic skeleton. Therefore we can not automatically infer that, being the prolongation structure finite dimensional, then the system is not integrable. The results in \cite{Hoenselaers} are a counterexample, which suggest that we could extend a finite dimensional prolongation structure in order to implement the skeleton structure of an integrable family of nonlinear systems. These aspects are related to Olver's symmetry reduction \cite{Olver}.

Moreover, we guess that the skeleton can be further implemented by a finer structure related to particle-like Lie algebra structure and this is the inspiring idea of our investigations.
In this note we show that such (otherwise discarded) symmetries deserve a more careful study. We take a $(2+1)$-dimensional Toda type system as a study case and
show that it posses algebraic properties related to the recently introduced concept of {\em particle-like Lie algebra structures} \cite{Vin17,Vin18}.

Vinogradov developed a completely abstract theory of compatibility of Lie algebra structures starting from the corresponding compatibility theory of Poisson structures.
Although the mathematical aspects of the theory are quite involved the nice point is that simple criteria of compatibility or non compatibility have been obtained which somehow have a certain grade of automatism.

Furthermore, as for the physical side, Vinogradov speculated that this particle-like structures could be related to the ultimate particle structure of the matter: he noted that since
{\quote
`the symmetry algebra $u(2) = so(3)$ of a nucleon can be assembled in one step
from three triadons [...] one might think that this structure of the symmetry reflects the fact that a nucleon is made from three ``quarks'' '.
\quote}
This is of course only a speculation, but it also suggests a quite fascinating new perspective on internal symmetries of integrable systems.

\section{Internal symmetries of Toda type systems in $(2+1)$ dimensions}

Consider the $(2+1)$-dimensional system, a continuous (or long-wave) approximation of a spatially two-dimensional Toda lattice \cite{Toda}:
\bEq\label{1}
u_{xx}+u_{yy}+(e^{u})_{zz}=0\,,
\eEq
where $u=u(x,y,z)$ is a real field, $x,y,z$ are real local coordinates (if we want, $z$ playing the r\^ole of a `time') and the subscripts mean partial derivatives.
It can be seen as the limit for $\gam\to \infty$ of the more general model
$$u_{xx}+u_{yy}+\left[\left(1+u/ \gam \right)^{\gam -1}\right]_{zz}=0$$
covering (for $\gam \neq 0,1$) various continuous approximations of  lattice models, among them the Fermi-Pasta-Ulam ($\gam =3$) \cite{Alfinito2}.
This model is almost {\em ubiquitus},  it appears in differential geometry; in mathematical and theoretical physics (Newman and Penrose);
in the theory of Hamiltonian systems; in general relativity; in the large $n$ limit of the $\textstyle{sl}(n)$ Toda lattice; in extended conformal
symmetries, and theory of gravitational instantons;
in strings theory and statistical mechanics {\em etc.} (see \eg \cite{Boyer,Lebrun,Park,Plebanski}).

It can be seen as the particular case with $d=1$ of so-called  $2d$-dimensional
Toda-type systems \cite{Saveliev} obtained from a `continuum Lie algebra' by means of a zero curvature representation $u_{w\bar{w}}= K(e^{u})$, (in our particular case $w =x+iy$ and $K$ is  the differential operator  given by $K= \frac{\der^{2}}{\der z^{2}}$).
In particular, it has been studied in the context of symmetry reductions \cite{Alfinito1,Grassi} and a $(1+1)$-dimensional version in the context of prolongation structures \cite{Alfinito2}. The $(2+1)$-dimensional system has been associated with a Ka\v c--Moody Lie algebra and related to Saveliev's continuum Lie algebras of particular kind \cite{PaWi10}.

The Toda system  \eqref{1} can be put in the complex form
\beq
\der_\zeta\der_{\bar\zeta} u =  -1/4 \der^{2}_{z} e^{u} \,,
\eeq
by the transformations
$\zeta = g(\eta)$, $\bar{\zeta} = \bar{g}(\bar{\eta})$, $u = \tilde{u}
- \ln(g^{'}\bar{g}^{'})$,
where $\zeta = x+ iy$,  $\bar\zeta = x- iy$, $\der_\zeta =\frac{1}{2}
(\der_x - i \der_y)$, $\der_{\bar\zeta}=\frac{1}{2}
(\der_x + i \der_y)$, $g^{\prime}=g_{\eta}(\eta)$,
$\bar{g}^{'}=g_{\bar\eta}(\bar\eta)$ and $g(\eta)$ is an arbitrary
holomorphic function of $\eta= x^{\prime}+iy^{\prime}$.
A Lax pair for this complex form of the $2D$ Toda equation has been found; see \eg Manakov and Santini \cite{MaSa09} and references therein; original references are  \cite{Saveliev}, as well as \cite{Zakharov}.

\subsection{Skeletons for the (2+1) Toda system}

Let us first recall a few mathematical tools constituting the background for a detailed treatment of which
we refer to \cite{PaWi02,PaWi10,PaWi11} and \cite{Mo93,PRS79}.

From one side global properties of partial differential equations such as internal symmetries and invariance properties having an issue in dynamics can be described by mathematical tools which enable us to deal with global properties at large scales, connecting local data to global ones.
On the other side transformations of configurations of a system can be globally studied by means of the theory of the action of Lie groups on manifolds.
The differential content carried by a Lie group (and its Lie algebra) and by its structure equations provides differential equations.

We observe that two ingredients constitute the  nonlinear phenomena: symmetries on the one side (algebraic content) and changes in time and space on the other side (differential content).
In particular, to keep account of the `interaction' of both aspects, we recognize a refined structure of open Lie algebraic structures associated with them: we introduce
a notion which generalizes the concept of a homogeneous space, \ie that of an {\em algebraic skeleton} $\bE = \mathfrak{g} \oplus \bV$ on a finite-dimensional vector space $\bV$, with $\mathfrak{g} $ a possibly infinite dimensional Lie algebra. The further step is introducing a tower with such a skeleton.

An {\em algebraic skeleton} on a finite-dimensional vector space $\bV$ is a triple $(\bE,\bG,\rho)$, with $\bG$ a (possibly infinite-dimensional) Lie group,
$\bE = \mathfrak{g} \oplus \bV$ is a (possibly infinite-dimensional) vector space {\em not necessarily equipped with a Lie algebra structure},
$\mathfrak{g}$ is the Lie algebra of $\bG$,
and
$\rho$ is a representation of $\mathfrak{g}$ on $\bE$ such that it reduces to the adjoint representation of $\mathfrak{g}$ on itself. The fact that $\bE$ is not a direct sum of Lie algebras, but an open algebraic structure is fundamental in order to be able to generate whole families of nonlinear differential systems, starting from it.

We now consider a suitably constructed differentiable structure which is somewhat modelled on the skeleton above. Let us introduce a differentiable manifold $\bP$ on which a Lie group $\bG$, with Lie algebra $\mathfrak{g}$, acts on the right; $\bP$ is a principal bundle $\bP\to\bZ\simeq\bP/\bG$.
By construction, we have that $\bZ$ is a manifold of type
$\bV$, \ie $\A \bz\in \bZ$,  $T_{\bz} \bZ \simeq \bV$.

Suppose we have a way to define a representation $\rho$ of the Lie algebra $\mathfrak{g}$ on $T_{\bz} \bZ \simeq \bV$, in such a way that it could be possible under certain conditions to find a
homomorphism between the open infinite dimensional Lie algebra, constructed by $\rho$, and a quotient Lie algebra.
Let us call $\mathfrak{k}$ the (possibly infinite dimensional) Lie algebra obtained as the direct sum of  such a quotient Lie algebra  with $\mathfrak{g}$.
From the differentiable side, a {\em tower $\bP(\bZ,\bG)$ on $\bZ$ with skeleton $(\bE, \bG,\rho)$ } is
an {\em  absolute parallelism} $\ome$ on $\bP$ valued in $\bE$,
invariant with respect to $\rho$ and reproducing elements of $\mathfrak{g}$ from the fundamental
vector fields induced on $\bP$, \ie  $R^{*}_{g}\omega = \rho(g)^{-1}\ome$, for
$g\in \bG$; $\ome(\tilde{A}) = A$, for $A\in\mathfrak{g}$;
here $R_{g}$ denotes the right translation and $\tilde{A}$ the fundamental
vector field induced on $\bP$ from $A$.
In general, the absolute parallelism  {\em does not}  define a Lie algebra homomorphism.

Let then $\mathfrak{k}$ be a Lie algebra and $\mathfrak{g}$ a Lie subalgebra of
$\mathfrak{k}$.
Let $\bG$ be a Lie group with Lie algebra $\mathfrak{g}$ and $\bP(\bZ, \bG)$ be a principal fiber bundle with structure group $\bG$ over a manifold $\bZ$ as above.
A {\em Cartan connection} in $\bP$ of type
$(\mathfrak{k}, \bG)$ is
a $1$--form $\ome$ on $\bP$ with values in
$\mathfrak{k}$ such that
$\ome |_{T_{\bp} \bP}: T_{\bp} \bP\to \mathfrak{k}$ is an isomorphism $\forall \bp \in
\bP$,
$R^{*}_{g}\ome=Ad(g)^{-1}\ome$ for $g\in \bG$
and
reproducing elements of $\mathfrak{g}$ from the fundamental
vector fields induced on $\bP$.
It is clear that a Cartan connection
$(\bP, \bZ, \bG, \ome)$ of type $(\mathfrak{k}, \bG)$ is a special case of a tower on $\bZ$.

The vector space $\bV$ is finite dimensional and generated by some of the vector fields in the prolongation structure. It has the property that each  bracket of some of remaining vector fields of the prolongation structure (freely generating an infinite dimensional Lie algebra $\mathfrak{g}$) with its generators is again in $\bV$.
In particular unknown commutators in the freely generated Lie algebra are related in such a way that their assigned relations are elements of $\bV$.

As an example of application of such an abstract formulation to the real world we refer  \eg to \cite{Pa16_Eco}, whereby activator-substrate systems have been obtained by performing a slight modification of the internal symmetry algebra of twisted reaction-diffusion equations: the necessary condition for the generation of stable patterns (related to general integrability properties in the limit of a null normalized diffusion constant) are  formulated in terms of `closeness' properties within the symmetry algebra vector space.

Following \cite{PaWi10}, we recall how to get both some skeletons and towers over them associated with the system \eqref{1}.

On a manifold with local coordinates
$(x,y,z,u,p,q,r)$, we introduce the closed differential ideal defined by the set of $3$--forms:
$
\theta_{1} = du \wedge dx \wedge dy - rdx \wedge dy \wedge
dz$,
$\theta_{2} = du \wedge dy \wedge dz - pdx \wedge dy \wedge
dz$,
$\theta_{3} = du \wedge dx \wedge dz + qdx \wedge dy \wedge
dz$,
$\theta_{4} = dp \wedge dy \wedge dz - dq \wedge dx \wedge
dz$ $ +$ $ e^{u}dr\wedge dx \wedge dy + e^{u}r^{2}dx \wedge dy
\wedge dz$.
It is easy to verify that on every integral submanifold defined by $u = u(x,y,z)$,
$p=u_{x}$, $q=u_{y}$, $r=u_{z}$, with $dx\wedge dy\wedge dz\neq 0$, the above ideal
is equivalent to the Toda system under study.

By an ansatz first introduced in \cite{PaWi02}, we  look for suitable $2$--forms (generating associated conservation laws)
\beq
\Omega^{k}= \tht^k_m\wed\ome^m
\eeq
where
$\tht^k_m = - \hat{A}^{k}_{m} dx-\hat{B}^{k}_{m}dy- \hat{C}^{k}_{m}dz$, with  $\hat{A}^{k}_{m}$, $\hat{B}^{k}_{m}$, $\hat{C}^{k}_{m}$ elements of $N \times N$
constant regular matrices, and the absolute parallelism forms are given by
\bEq\label{omega}
\ome^m= d\hat{\xi}^m+ \hat{F}^m dx  +\hat{G}^m dy  +  \hat{H}^m dz\,,
\eEq
\ie
\bEq \label{Omega}
\Omega^{k} = H^{k}(u, u_{x}, u_{y}, u_{z}; \xi^{m})dx \wedge dy+
F^{k}(u, u_{x},u_{y}, u_{z}; \xi^{m})dx \wedge dz + \\
+ G^{k}(u, u_{x}, u_{y}, u_{z}; \xi^{m})dy \wedge dz+ A^{k}_{m}d \xi^{m} \wedge dx
+ B^{k}_{m}d \xi^{m} \wedge dz + d \xi^{k}
\wedge dy\,, \nonumber
\eEq
where $\xi=\{\xi^{m}\}$, $k,m=1,2,\ldots,
{\rm N}$ (N arbitrary), and $H^{k}$, $F^{k}$ and $G^{k}$ are, respectively, the
pseudopotentials and functions to be determined, while $A^{k}_{m}$ and $B^{k}_{m}$ denote the elements of two $N \times N$
constant regular matrices related to the previous ones and we have rescaled the coordinates $\xi^k$.
In particular note that (see also \cite{PaWi02,Pa05})
\bEq
F^k = \hat{C}^{k}_{m}\hat{F}^m- \hat{A}^{k}_{m}  \hat{H}^m  \label{relations1} \,, \\
G^k = \hat{C}^{k}_{m}\hat{G}^m- \hat{B}^{k}_{m}  \hat{H}^m  \label{relations2} \,, \\
H^k = \hat{B}^{k}_{m}\hat{F}^m- \hat{A}^{k}_{m}  \hat{G}^m \label{relations3} \,, \\
\xi^k = \hat{C}^{k}_{m}\hat{\xi}^m \,.
\eEq

The integrability  condition for the ideal generated by forms
$\theta_{j}$ and $\Omega^{k}$ finally yields
\bEq\label{H}
H^{k} = e^{u}u_{z}L^{k}(\xi^{m}) + P^{k}(u,\xi^{m}) \,,
\eEq
\bEq\label{F}
F^{k} = -\,u_{y}L^{k}(\xi^{m}) + Q^{k}(u, \xi^{m}) \,,
\eEq
\bEq\label{G}
G^{k} = u_{x}L^{k}(\xi^{m}) + M^{k}(u,\xi^{m}) \,,
\eEq
where $L^{k}$, $P^{k}$, $Q^{k}$, $M^{k}$ are functions of integration.

It turns out that $Q^{k}(u, \xi^{m})$ can be written in terms of the others.
Indeed we have (see \eg \cite{Morris,Tondo})
$H^{k}  =A^{k}_{m} G^{m} - B^{k}_{m} F^{m} $
so that $ e^{u}u_{z}L^{k}(\xi^{l}) + P^{k}(u,\xi^{l})=A^{k}_{m} (u_{x}L^{m}(\xi^{l}) + M^{m}(u,\xi^{l})) - B^{k}_{m}(-\,u_{y}L^{m}(\xi^{l}) + Q^{m}(u, \xi^{l}))$, \ie
$B^{k}_{m}Q_u^{m}(u, \xi^{l})= e^{u}u_{z}L^{k}(\xi^{l}) + P_u^{k}(u,\xi^{l})+ A^{k}_{m}  M_u^{m}(u,\xi^{l})$, which can be integrated once the dependence on of $P^{k}(u,\xi^{l})$ and $M^{k}(u,\xi^{l})$ on $u$ is given.
As a consequence, the desired representation $\rho$ for the skeleton is provided by the following equations (we omit the indices for simplicity) \cite{Palese,PaWi10}.
\bEq
P_{u} = e^{u}[L,M]\,,\label{P} \, \quad
 M_{u} = - [L,P]\,,\label{M} \, \quad
 \left[M,P\right] = 0\,.\label{comm}
\eEq
Note that here  $L$ depends only on $\xi^{m}$, while $P$ and $M$ still have a dependence on $u$ determined by the first two differential equations. A tower with $P$ and $M$ given in terms of $L$ has been obtained by suitable {\em operator Bessel coefficients} \cite{Palese}.

Note that formally this tower shall provide the Lax pair of an inverse spectral problem; however, it is a non trivial task to characterize explicitly its algebraic skeleton by means of the representation provided by  the relations $\left[M,P\right] = 0$, \ie to obtain a spectral problem in a manageable form.

Particular choices for the absolute parallelism can provide us explicit representations of the prolongation skeleton; in particular a Ka\v c--Moody Lie algebra has been obtained \cite{PaWi10} (see Proposition \ref{dyonsProposition}, case $1. \, (b)$ below).
In the following we will concentrate on those choices that generate particle-like Lie algebra structures.
We shall see that it is yet possible to obtain a spectral problem with a particular choice of the tower.

\section{Particle-like Lie algebra structure}

Recently, Vinogradov proved that any Lie algebra over an algebraically closed field or over $\R$ can
be assembled in a number of steps from two elementary constituents, that he called dyons
and triadons \cite{Vin17}. He  considered the problems of the construction and classification of
those Lie algebras which can be assembled in one step from base dyons and triadons,
called coaxial Lie algebras. The base dyons and triadons are Lie algebra structures
that have only one non-trivial structure constant in a given basis, while coaxial Lie
algebras are linear combinations of pairwise {\em compatible} base dyons and triadons \cite{Vin18}. Here for the convenience of the reader we recall some basic facts of the theory in the original Vnogradov's notation.

\bDf
Lie algebra structures $\mathfrak{g}_1$  and $\mathfrak{g}_2$ on a vector space $ V$ are called compatible
if $[,]_{\mathfrak{g}_1}+[,]_{\mathfrak{g}_2}$ is also a Lie algebra product.

A Lie algebra $\mathfrak{g}$ is called simply assembled from Lie algebra structures $\mathfrak{g}_1,\dots,\mathfrak{g}_m$ on $|\mathfrak{g}|=V$ if
the Lie algebras  $\mathfrak{g}_i$'s are pairwise compatible and
$[,]_{\mathfrak{g}}=[,]_{\mathfrak{g}_1}+\dots [,]_{\mathfrak{g}_m}$ .
Note that if the Lie algebras $\mathfrak{g}_i$'s are compatible, then any linear combination of compatible Lie algebras commutators is a Lie  algebra commutator (or product) .
\eDf

\bDf
Fix a basis $B=e_1,\dots e _n$ in the representation vector space of a given Lie algebra. Let $i,j,$ and $k$ be integers, $1\leq i,j,k \leq n$, no two of them equal,
and denote by $\{i,j | k\}$ (respectively, $\{i | j\}$) the $n$-triadon (respectively, the $n$-dyon) such that $[e_i,e_j] = - [e_j,e_i] = e_k$ (respectively, $[e_i,e_j] = - [e_j,e_i] = e_j$) are the only non-trivial Lie commutators of basis vectors.
Vinogradov called  them  `base triadon' and `base dyon', respectively
or by the unifying term `base lieon'.

An $n$-dyon is the direct sum of a dyon with an $n-2$-dimensional abelian Lie algebra, $n\geq 2$, (\ie there is only one  non vanishing bracket and it is a dyon).
Analogously an $n$-triadon is the direct sum of a triadon with an $n-3$-dimensional abelian Lie algebra, $n\geq 3$ (\ie there is only one  non vanishing bracket and it is a triadon). They can also be referred generically as $n$-lieons.

A linear combination of pairwise compatible base lieons is called a {\em coaxial Lie algebra structure}.
A Lie algebra structure will be called {\em trix-coaxial} (respectively, {\em dyx-coaxial}) if it consists only of base triadons (respectively,
base dyons).
A coaxial Lie algebra $\mathfrak{g}$ may be presented as a linear combination,
\beq
 \mathfrak{g}=\sum{}{} \alpha_{(i,j|k)} \{i, j|k\} + \sum{}{} \bet_{(m|n)} \{m|n\}
\eeq
of pairwise compatible base lieons.

The vectors $e_i,e_j $, and $e_k$ (respectively, $e_i,e_j$) are called the {\em vertices of the triadon} $\{i,j | k\}$
(respectively, {\em of the dyon }$\{i | j\}$). The vectors $e_i$ and $e_j$ are called  the {\em ends of the triadon} $\{i,j | k\}$, while  $e_k$ is the {\em center of the triadon}.
The {\em origin and the end  of  the dyon $\{i | j\}$} are $e_i$ and $e_j$, respectively.  The base
triadons $\{i,j | k\}$ and $\{j,i | k\}$= - $\{i,j | k\}$ are not distinguished since they have identical compatibility properties.
\eDf
We now recall Proposition $3.1$ of \cite{Vin18} stating some necessary and sufficient conditions for the compatibility or incompatibility  of particle-like Lie algebra structures:
\begin{itemize}
\item Two base triadons are non-trivially compatible if and only if they have a
common center, a common end, or both.
\item Two base dyons are incompatible if and only if the origin of one is the end of the other and
they have no other common vertices.
\item A base dyon is non-trivially compatible with a base triadon if and only if its origin coincides
with one of the ends of the triadon.
\end{itemize}
For further notation and vocabulary we refer the reader to Vinogrados's papers.

\subsection{Trix-coaxial, dyx-coaxial and particle-like Lie algebra structures for the Toda system}

We prove that with a $(2+1)$-dimensional Toda type system are associated algebraic skeletons which are compatible  assemblings  of particle-like Lie algebras of dyons and triadons type.  We obtain trix-coaxial and dyx-coaxial Lie algebra structures for the system from skeletons of some particular choice for compatible associated absolute parallelisms. In particular, we find a trix-coaxial Lie algebra structure
made of two (compatible) base triadons constituting a $2$-catena (see Proposition 3.1, pag 5 \cite{Vin18}).

Let us indeed now look for special skeletons.

\subsection{Trix-coaxial  Lie algebra structures}

\bPr
Associate with the Toda type system  \eqref{1} is a trix-coaxial Lie algebra structure made of two (compatible) base triadons constituting a $2$-catena.
\ePr

\bPf
If we look for operators  $P(u,\xi)= e^{u}\bar{P}(\xi)$, $M(u,\xi)=M(e^u,\xi)$,  we get
$M( e^u;\xi) $ $=$ $ - e^u[L(\xi),\bar{P}(\xi)]+ \bar{M}(\xi)$ and thus
$\bar{P}(\xi)$ $=$ $ - e^u[L(\xi),[L(\xi),\bar{P}(\xi)]]+[L(\xi),\bar{M}(\xi)]$.
There are additional relations determined by the third prolongation equation
$[- e^u[L(\xi),\bar{P}(\xi)]+ \bar{M}(\xi), e^u \bar{P}(\xi)]=0$.

Let us then put $L=X_1,\bar{M}=X_2, \bar{P}=X_3, [X_1 ,X_3 ]= X_4 $. From the above we  have the following prolongation closed Lie algebra
\beq
[X_1 ,X_2 ]=X_3 \,,   [X_1 ,X_3 ]= X_4 \, , [X_1 , X_4 ]= [X_2 ,X_3 ]= [X_2 ,X_4 ]= [X_3 ,X_4 ]= 0 \,.
\eeq
The above is a trix-coaxial Lie algebra structure made of two compatible $4$-triadons.

Indeed,
by taking $X_4=0$, we get
 $ [X_1 ,X_2 ]=X_3 \,,   [X_1 ,X_3 ]=   [X_2 ,X_3 ]= 0$ and  $ [X_1 , X_4 ]=   [X_2 ,X_4 ]= [X_3 ,X_4 ]= 0$ trivially.

On the other hand by taking $X_2=0$, we get
$[X_1 ,X_3 ]= X_4 \,, [X_1 , X_4 ]=   [X_3 ,X_4 ]= 0$ and  $[X_1 ,X_2 ]= [X_2 ,X_3 ]= [X_2 ,X_4 ] =0 $ trivially.

According with \cite{Vin18} the two $4$-triadons above {\em are non trivially compatible having a common end $X_1$}, and they constitute a $2$-catena.
\ePf

\subsubsection{Conservation laws and special solutions associated with a $2$-catena}\label{Conservation laws}

Let us now explicate the tower corresponding to such $4$-triadons.

For the sake of simplicity let us put $A^{k}_{m}=B^{k}_{m}=\delta^{k}_{m}$, were $\delta^{k}_{m}$ is the Kronecker symbol.
By substituting the above commutators into  equations \eqref{H} and \eqref{G} (the expression of \eqref{F} being constrained in this case by the relation $F=G+H$), we get

\bEq\label{H2}
H = e^{u}u_{z}X_1+ e^{u} X_3 \,,
\eEq
\bEq\label{G2}
G = u_{x}X_1 - e^{u} [ X_1, X_3]+ X_2 \,,
\eEq

Now from equation \eqref{Omega}, by sectioning we obtain

\bEq\label{boh}
H^{k}  -  \xi^k_y = -  \xi^k_x \,,
\eEq
\bEq\label{boh2}
G^{k} + \xi^k_y =  \xi^k_z \,,
\eEq
(together with $F^{k} + \xi^k_x =  \xi^k_z $ which depends on the two others).

We note that each one of  $4$-triadons above can be represented in a space of local coordinates $\xi^k$ providing conservation laws related to two compatible Poisson structures.

Indeed let us consider the $4$-triadon given by $X_2=0$. A representation in the coordinates $\{\xi^1,\xi^2,\xi^3\}$ is given by
$X_1= \xi^2 \der/\der \xi^1$, $X_3= \xi^3 \der/\der \xi^2$ and $X_4= -\xi^3 \der/\der \xi^1$.

The tower corresponding to this case gives
\beq\label{H2'}
 e^{u}u_{z}\xi^2 \der/\der \xi^1+ e^{u}  \xi^3 \der/\der \xi^2  -  \xi^k_y \der/\der \xi^k= -  \xi^k_x \der/\der \xi^k  \,,
\eeq
\beq\label{G2'}
 u_{x}\xi^2 \der/\der \xi^1 + e^{u} \xi^3 \der/\der \xi^1+ \xi^k_y\der/\der \xi^k =  \xi^k_z\der/\der \xi^k\,,
\eeq
which gives us the system
\bEq
\bxi_x= \bxi_y +\bM \bxi \, \\
\bxi_z= \bxi_y + \bN \bxi \,
\eEq
where $\bxi =(\xi^1,\xi^2,\xi^3)^T$,  $\bM$ and $\bN$ are $3\times 3$ matrices such that $M_{12}= -e^u u_z$,$ M_{23}= -e^u$, $N_{12}=u_x$, $N_{13}=e^u$, and all the other entries are zeros.
In view of \eqref{relations1}-\eqref{relations3} and \eqref{omega} this system can be interpreted as a conservation law.
It is a  first order system which can be manipulated by resorting to the method of characteristics so that it turns out to be useful to find out special solutions of the Toda system.

Note indeed that this system is equivalent to the following system of coupled equations of Maxwell type
\bEq
& & \xi^1_{yy} - \xi^1_{xx} = (e^u u_z\xi^2)_x +(e^u u_z\xi^2)_y \,, \\
& & \xi^2_{yy} - \xi^2_{xx} = (e^u \xi^3)_x +(e^u \xi^3)_y \,,  \label{second eq} \\
& & \xi^1_{zz} - \xi^1_{yy} =  (e^u \xi^3 +u_x\xi^2)_z + (e^u \xi^3 +u_x\xi^2)_y\,,
\eEq
where $e^u u_z\xi^2$, $e^u \xi^3$ and $e^u \xi^3 +u_x\xi^2$ can be recognized as charge/current densities.
On the other hand, the system can be simplified since $\xi^2_y = \xi^2_z$; we have that \eqref{second eq} can also be written as
\bEq
\xi^2_{zz} - \xi^2_{xx} = (e^u \xi^3)_x +(e^u \xi^3)_y \,,
\eEq
from which we obtain the Maxwell-type equation
\bEq
\xi^2_{yy} = \xi^2_{zz}  \,.
\eEq
Further  manipulations can be made by using $\xi^3_y = \xi^3_y = \xi^3_z$.

We remark that the same conservation law and related outcomes are obtained by the $4$-triadon given by $X_4=0$.
Therefore existence of that tower with a finite dimensional skeleton which is a $2$-catena says us that the two Poisson structures corresponding to each $4$-triadon are compatible also in a sense which is interpretable from a physical point of view: they are structures associated with the same Toda system,  and more precisely with the same conservation law and related special solutions.
Compatibility of Poisson structures is beyond the scope of this paper, however this result suggest interesting links between special solutions and compatible Poisson structures, which will be the object of further investigations (in particular for meron-like configurations or gravitational instantons).

\subsection{Dyx-coaxial and particle-like Lie algebra structures}

In the following we analyze with more detail the case of choice  $P(u,\xi)= \ln u \bar{P}(\xi)$, $M(u,\xi)=M(e^u,\xi)$ studied in \cite{PaWi10} also leading to an infinite dimensional skeleton homomorphic to a Ka\v{c}-Moody Lie algebra. We carefully distinguish the various cases.

This choice of the absolute parallelism associates with the Toda system \eqref{1} dyx-coaxial  and particle-like Lie algebra structures.

First we need a preliminary result (see also \cite{PaWi10}).
\bLm\label{lemma}
Let $P(u,\xi)= \ln u \bar{P}(\xi)$, $M(u,\xi)=M(e^u,\xi)$. We get the following infinitesimal algebraic  skeleton with the structure of an {\em open Lie algebra}:
\bEq\label{open}
[X_1, X_2]=X_4\,, [X_1, X_3]=X_5\,, [X_4, X_5]
=[X_2,X_7] \,, [X_3, X_4] =  [X_2, X_5] \,,
\eEq
\beq
[X_1, X_4]=X_6\,,  [X_1, X_5]=X_7\,, [X_2, X_3] = X_8  \,,
\eeq
\beq
[X_1,X_8]= [X_2, X_4] =  [X_2, X_6] =[X_3, X_7] = 0\,,
\ldots
\eeq
\eLm

\bPf
Put  $L=X_1 (\xi)$.

By derivation
we get
$M(e^u,\xi)=- (\ln u -1)u[X_1 (\xi),X_3(\xi)]+X_2 (\xi)$,
and
$P(u,\xi)= ue^u \ln u[X_1(\xi) , M] $.

For $u \neq 0,1$ (which are trivial solutions of the Toda system), from $[P,M]=0$
we get
\beq
[[X_1, M],M]=0\,;
\eeq
from which we get
\beq
[[X_1,X_2],X_2 ] = 0 \,,
\,
[X_1,[X_1 ,X_3 ]], X_2 ]+ [[X_1 ,X_2 ],[X_1 ,X_3 ]] =0\,,
\eeq
\beq
[[X_1,[X_1 ,X_3 ]],[X_1 ,X_3]] = 0\,.
\eeq
By putting
$[X_1, X_2]=X_4$, $[X_1, X_3]=X_5$, $[X_1, X_4]=X_6$, $[X_1, X_5]=X_7$, $[X_2, X_3] = X_8$, we obtain an infinite dimensional skeleton as follows
\beq
 [X_1, X_8]= [X_2, X_4] =  [X_2, X_6] =[X_3, X_7] = 0\,,
 \eeq
 \bEq\label{commutator relations}
  [X_4, X_5] =[X_2,X_7] \,, [X_3, X_4] =  [X_2, X_5] \,,
   \eEq
 \beq
\ldots
\eeq
\ePf
Here the dots means that we can continue this structure by introducing new generators still obtaining the peculiar relations  of the type \eqref{commutator relations} which distinguish this algebraic structure from a freely generated Lie algebra (see the discussion in \cite{Pa16_Eco}).

\bPr\label{dyonsProposition}
The homomorphism  $X_4=\lam X_2$ and $X_5 =\mu X_3$ associates with the Toda system \eqref{1} dyx-coaxial and particle-like Lie algebra structures   as well as an infinite-dimensional Lie algebra homomorphic with a Ka\v{c}-Moody Lie algebra.
\ePr

\bPf
We essentially distinguish the two cases $X_8\neq 0$ and $X_8 = 0$, together with various different subcases.
\begin{enumerate}
\item
if $[X_2, X_3] = X_8\neq 0$, then $\mu = -\lam$ must old; we can distinguish different subcases
 \begin{enumerate}
 \item in general  the case $X_8\neq 0$ and $\mu= -\lam \neq 0$ can provide infinite-dimensional Lie algebras homomorphic with Ka\v{c}-Moody type Lie algebras.
\item  the particular case $X_8=\nu X_3$ and $\mu = -\lam = 1$ (\ie $X_4=X_2$ and $X_5 = -X_3$) giving an infinite-dimensional Lie algebra homomorphic with a Ka\v{c}-Moody Lie algebra was obtained in \cite{PaWi10}.
 \item
 the particular case $X_8=\nu X_3$ and  $\mu= -\lam = 0$ (\ie $X_4=X_5=0$; see \cite{PaWi10})  gives a particle-like Lie algebra as a base $3$-dyon:
\bEq\label{dyon}
[X_1, X_2]=0\,, [X_1, X_3]=0\,, [X_2, X_3] = \nu X_3\,.
\eEq
\end{enumerate}
\item  if  $[X_2, X_3] = X_8 = 0$,
then $X_6=X_4$, $X_7=X_5$, and we distinguish the following different subcases (the case $\mu= \lam = 0$ giving an abelian Lie algebra):
 \begin{enumerate}
 \item
the case $\mu=0$ and $\lam\neq 0$ provides us with a particle-like Lie algebra as a base $3$-dyon:
\bEq\label{3-dyon}
[X_1, X_2]=\lam X_2\,, [X_1, X_3]=0\,, [X_2, X_3]  =0 \,.
\eEq
\item
the case
$\lam =\mu \neq 0$ provides a  dyx-coaxial Lie algebra structure as an assembling of two compatible base  $3$-dyons
\bEq\label{dyx-family}
[X_1, X_2] = \lam X_2\,, [X_1, X_3] =  \lam X_3\,, [X_2, X_3] = 0 \,.
\eEq
 \item the particular case  $\lam =\mu=1$ (\ie $X_4=X_2$ and $X_5 =X_3$) gives
\beq
[X_1, X_2]=X_2\,, [X_1, X_3]=X_3\,,
[X_2, X_3] = 0 \,,
\eeq
and it  was obtained in  \cite{PaWi10}.
\end{enumerate}
\end{enumerate}
\ePf

\bRm
We can try to check the compatibility of the above particle-like  Lie algebra structures.
The latter Lie algebra \eqref{dyx-family} is constituted of two mutually compatible dyons (dyx-coaxial Lie algebra), whose the first is given by \eqref{3-dyon}, while the Lie algebras \eqref{dyon} and \eqref{3-dyon} are made of a single dyon and they are not compatible.
Indeed we note that the first dyon of \eqref{dyx-family} is not compatible with \eqref{dyon}, while the second dyon of \eqref{dyx-family} is.

The question now is can we still construct a different  dyx-coaxial Lie algebra  from the original skeleton? For example the following would be a dyx-coaxial Lie algebra of compatible dyons
\bEq
[X_1, X_3] =  \lam X_3, [X_2, X_3] = \nu X_3, [X_1, X_2] = 0 \,
\eEq
We ask if we can get it from the prolongation skeleton by a suitable quotienting, \ie if it is somehow compatible with (or derivable from) the skeleton structure.
However, we note that if we put $X_4=0$ from the beginning (which is the case when we assume that  $ [X_1, X_2] = 0$), and  if  $X_8 =\nu X_3$, this would  imply also $[X_1, X_3] = 0$, then we  would get \eqref{dyon} back.

Thus it appears that the case $X_8 =\nu X_3$ corresponds or to a particle-like Lie algebra structure or to a Ka\v{c}-Moody type Lie algebra (it is noteworthy that the latter is anyway  an infinite-dimensional loop Lie algebra of a dyx-coaxial Lie algebra) and these two cases appear to be non compatible.

Let us then investigate from a more general point of view this feature.
We ask whether  we can look for different quotient homomorphisms.

Let now consider the case $X_4 = 0$ from the beginning, and $X_8\neq 0$, and look for a quotient lie algebra given by
$X_8 = - \gamma X_2$, , $X_5=  \mu X_3$, and we obtain the  Lie algebra structure depending on two parameters
\bEq\label{incompatible}
[X_1, X_2] = 0 \,, [X_1, X_3] = \mu X_3 \,, [X_3, X_2] = - \gamma X_2 \,.
\eEq
By appling the Jacobi identity we get $ \mu \gamma X_2 =0$ which, if we require $X_2\neq0$, is verified  either for  $\mu =0$ and $ \gamma \neq 0$ (see below \eqref{other dyon}) or for $\mu \neq0$ and $ \gamma =  0$ (see below \eqref{other dyon 2}), or for $\mu=0$ and $ \gamma= 0$ (trivial case of an abelian Lie algebra).

\bPr \label{new dyon}
The case $X_4 = 0$ from the beginning, and with $X_8\neq 0$,  provides us with two base $3$-dyons.
\ePr

\bPf
\begin{enumerate}
\item the case with  $\mu =0$ and $ \gamma \neq 0$.

By putting   $X_8 = - \gamma X_2$,  and $X_4=X_5= X_6=X_7=0$ we get the $3$-dyon
\bEq\label{other dyon}
 [X_1, X_2] = 0 \,, [X_1, X_3] = 0 \,, [X_3, X_2] = - \gamma X_2 \,.
\eEq
The above dyon is incompatible with \eqref{dyon} while it is compatible with \eqref{3-dyon}.

\item the case $\mu \neq 0$ and $ \gamma =  0$.

We get the $3$-dyon
\bEq\label{other dyon 2}
 [X_1, X_2] = 0 \,, [X_1, X_3] = \mu X_3 \,, [X_3, X_2] = 0 \,.
\eEq

\end{enumerate}
\ePf

\bRm

It appears that the dyx-coaxial Lie algebra \eqref{dyx-family} can be assembled by one step from the case \eqref{3-dyon} and the latter one,  \eqref{other dyon 2}, by putting $\mu=\lam$,.
\eRm
We note in particular  that \eqref{incompatible} can not be seen as a dyx-family of  dyons since the two  dyons $[X_1, X_3] = \mu X_3$, $[X_3, X_2] = - \gamma X_2$ are incompatible and indeed if we apply the Jacobi identity we get particle-like structures made of single base dyons.

It seems therefore that the prolongation skeleton is homomorphic with quotient finite dimensional Lie algebras which have always the structure of a family of compatible dyons or single base $3$-dyons.
We note that the first dyon of \eqref{dyx-family} is compatible with \eqref{other dyon}, while the second dyon of \eqref{dyx-family} is not.

Summing up we were able to associate with the infinitesimal skeleton \eqref{open} a dyx-coaxial Lie algebra structure \eqref{dyx-family} and  particle-like Lie algebra structures made of three base $3$-dyons which are only partially compatible among them, \ie
\begin{itemize}

\item the first dyon of \eqref{dyx-family} is compatible with \eqref{other dyon}, while the second dyon of \eqref{dyx-family} is not.

\item the first dyon of \eqref{dyx-family} is not compatible with \eqref{dyon}, while the second dyon of \eqref{dyx-family} is.

\end{itemize}

Note that \eqref{dyon}, \eqref{3-dyon}, \eqref{other dyon} and \eqref{other dyon 2} are not all compatible among them, even they are not compatible in triples, but they are only  compatible when took in couples.
\eRm

\subsubsection{A Lax pair associated with the dyx-coaxial Lie algebra structure of two compatible base  $3$-dyons}\label{Lax section}

Following a procedure similar to that of Section \ref{Conservation laws} we shall now derive a Lax pair  related to the Lie algebra \eqref{dyx-family}  to which the tower skeleton \eqref{commutator relations} is homomorphic.

We refer again to equations \eqref{H} and \eqref{G}. The tower associated with case $2. (b)$ in Proposition \ref{dyonsProposition} becomes in this specific case
\bEq\label{H3}
H = e^{u}u_{z}X_1- u^2  e^{u} \ln u( \ln u -1 )[ X_1, [X_1,X_3 ] ] +  u  e^{u} \ln u [X_1,X_2] \,,
\eEq
\bEq\label{G3}
G = u_{x}X_1 - u(\ln u -1)[X_1,X_3]+X_2\,,
\eEq
By taking into account the representation $\rho$ given by relations \eqref{dyx-family}
we then get
\bEq\label{H3'}
H = e^{u}u_{z}X_1 -  \lam^2 u^2  e^{u} \ln u( \ln u -1 ) X_3  + \lam u  e^{u} \ln u  X_2 \,,
\eEq
\bEq\label{G3'}
G = u_{x}X_1 - \lam u(\ln u -1)  X_3 +X_2 \,,
\eEq

Now, let us represent the dyx-coaxial Lie algebra above in a space of `{\em pseudopotentials}' $\xi^k$ by
$X_1 = - \xi^1\der/\der\xi^1+\xi^2\der/\der\xi^2$, $X_2 = - \lam \xi^3\der/\der\xi^1$, $X_3 = \lam \xi^2\der/\der\xi^3$.

Again by sectioning the tower, equations \eqref{boh} and \eqref{boh2} provide the following
\beq\label{H2''}
& & e^{u} (u_{z}\xi^1 + \lam^2 u  \ln u \xi^3)   \der/\der\xi^1- e^{u}u_{z} \xi^2\der/\der\xi^2  + \\ & &+  \lam^3 u^2  e^{u} \ln u( \ln u -1 )  \xi^2\der/\der\xi^3    +  \xi^k_y \der/\der \xi^k= \xi^k_x \der/\der \xi^k  \,, \nonumber
\eeq
\beq\label{G2''}
& & - u_{x}( \xi^1 + \lam \xi^3 )\der/\der\xi^1+ u_{x}\xi^2\der/\der\xi^2  - \lam^2 u(\ln u -1)  \xi^2 \der/\der\xi^3 + \\ & &  + \xi^k_y\der/\der \xi^k =  \xi^k_z\der/\der \xi^k\,, \nonumber
\eeq
which gives us the inverse spectral problem
\bEq
\bxi_x= \bxi_y + {\hat \bM} \bxi \, \\
\bxi_z= \bxi_y + {\hat \bN} \bxi \,
\eEq
Where $\bxi =(\xi^1,\xi^2,\xi^3)^T$,  ${\hat \bM}$ and ${\hat \bN}$ are $3\times 3$ matrices such that $\hat{M}_{11}= e^{u} u_{z}$, $\hat{M}_{13}= e^u\lam^2 u  \ln u$, $\hat{M}_{22}= - e^{u}u_{z} $, $\hat{M}_{32} = \lam^3 u^2  e^{u} \ln u( \ln u -1 )$
$\hat{N}_{11} = -u_x$, $ \hat{N}_{13}=  - \lam u_x $, $ \hat{N}_{22} = u_x$, $ \hat{N}_{32} = - \lam^2 u(\ln u -1)$ and all the other entries are zeros.
Here $\lam$ plays the role of a spectral parameter and, in view of \eqref{relations1}-\eqref{relations3} and \eqref{omega}, ${\hat \bM}$ and ${\hat \bN}$ can be considered a Lax pair related to the multidimensional Toda system  \eqref{1} (for spectral problems related to multidimensional nonlinear system see, \eg   \cite{Morris,Tondo}). This Lax pair should be compared with \cite{MaSa09}.

Compatibility of Lie algebraic structures being expressions of  compatibility of the corresponding Poisson structures, we note  here that Fernandes \cite{Fernandes93} studied the relationship between the master symmetries and bi-Hamiltonian structure of the Toda lattice.
We stress that {\em dyons provide indeed particular examples of master symmetries} of related ordinary differential equations.

\subsection{Concluding remarks}

The structure of trix-coaxial and dyx-coaxial Lie algebras assembled in one step from couples of particle-like Lie algebra structures appears as an intrinsic feature of the Toda system \eqref{1}, at least associated with the chosen absolute parallelisms.
Indeed the similitude transformations seems to be the fundamental internal symmetries of the system (see \eg \cite{Alfinito1}).

As final remark, since \eqref{other dyon} is compatible with \eqref{3-dyon},  and since \eqref{dyon} is compatible with \eqref{other dyon 2}, we could construct the following dyx-coaxial Lie algebras:
\bEq\label{3-dyon 2}
[X_1, X_2]=\lam X_2\,, [X_1, X_3] = 0\,, [X_3, X_2]  = - \gamma X_2 \,,
\eEq
and
\bEq\label{3-dyon 3}
[X_1, X_2]=0\,, [X_1, X_3] = \mu X_3\,, [X_2, X_3] = \nu X_3\,.
\eEq
However, it is important to realize that they could {\em not} be obtained from the skeleton \eqref{open} by the choice of an homomorphism, and therefore they are not identified as internal symmetries of the Toda system by the choice of the absolute parallelism given by Lemma \ref{lemma}.
The question if the choice of other forms of the absolute parallelism could identify them is open and will be the object of future investigations.

\section*{Acknowledgements}
Research partially supported by Department of Mathematics - University of Torino through the  projects PALM$\_$RILO$\_16\_ 01$ and FERM$\_$RILO$\_17\_ 01$ (MP)
 and written under the auspices of GNSAGA-INdAM.
The first author (MP) would like to acknowledge the contribution of the COST Action CA17139. The authors would like to thank M.V. Pavlov for bringing to their attention the reference \cite{Zakharov}.

config: arxiv


TITLE: Taylor series with point of expansion.
QUESTION [0 upvotes]: Write the Taylor series for the function $f(x)= x^2-3x+1$ using $x=2$ as the point of expansion , i.e. write a formula for $f(2+h)$
Even though I know what Taylor Expansion is I am really confused what the point of expansion is and how to find $f(2+h)$

REPLY [3 votes]: You're probably accustomed to Taylor series about $x=0$.  In general, the Taylor series about $x=a$ is $$f(a)+\frac{f'(a)}{1!}(x-a)+\frac{f''(a)}{2!}(x-a)^2+\cdots.$$
So, we compute the derivatives $f'(x)$, $f''(x)$, and so on (since this is a polynomial, it'll have a finite number of non-zero derivatives).  We use them to find the values $f'(a)$, $f''(a)$, and so on, then substitute them into the above equation.  In this case $a=2$.


text:
\section{Numerical Experiments}\label{sec:numexp}
In this section we compare the mRKC scheme from \cref{sec:fdrkcsq} against the classical RKC method of \cref{sec:rkc} through a series of experiments. First, we apply mRKC to a stiff nonlinear dynamical system to verify convergence in the standard ``ODE sense'' and underpin its efficiency. Then, we apply mRKC to the heat equation and verify convergence in the ``PDE sense'', i.e. when both the mesh size $H$ and the time step $\tau$ decrease simultaneously. In the third experiment, we compare the performance and efficiency of the mRKC and RKC schemes when applied to a linear diffusion problem in complex geometry; here, we also compare mRKC to a second-order accurate RKC scheme (RKC2) \cite{SSV98,VHS90} and the implicit Euler method. In the fourth experiment, we apply the mRKC scheme with the RKC, RKC2 and the implicit Euler method to a nonlinear integro-differential problem.
Finally, we study numerically the stability of mRKC when it is applied to various advection-diffusion-reaction problems.

Both the RKC and mRKC methods need bounds on the spectral radii of the Jacobians of $\ff$ and $\fs$
 to determine the number of stages $s,m$ needed for stability. In our experiments, we estimate them with a cheap nonlinear power method \cite{Lin72,Ver80}. The numerical experiments in \cref{sec:exp_conv2,sec:exp_eff2,sec:exp_stab} were performed using the C++ library \texttt{libMesh} \cite{KPS06}, while for the experiments of \cref{sec:exp_conv0,sec:intdiff} we used the \texttt{Eigen} library \cite{GuB10}.

\subsection{Robertson's stiff test problem}\label{sec:exp_conv0}
First, we study the convergence of the mRKC scheme on a popular stiff test problem, Robertson's nonlinear chemical reaction model \cite{EHL75,HaW02}:
\begin{equation}\label{eq:robertson}
\begin{aligned}
y_1' =& -0.04\, y_1+10^4\, y_2\, y_3, & y_1(0)=&1, \\
y_2' =& \,0.04\, y_1-10^4\, y_2\, y_3-3\cdot 10^7\, y_2^2,\qquad \qquad & y_2(0)=&2\cdot 10^{-5},\\
y_3' =& \,3\cdot 10^7\, y_2^2 , & y_3(0)=& 10^{-1},
\end{aligned}
\end{equation}
where $t\in [0,100]$. With this set of parameters and initial conditions, the only term inducing severe stiffness is $-10^4\,y_2\,y_3$. Thus, we let
\begin{align}
\ff(y)=& \begin{pmatrix}
0 \\ -10^4\, y_2\, y_3\\ 0
\end{pmatrix}, &
\fs(y) =& \begin{pmatrix}
-0.04\, y_1+10^4\, y_2\, y_3\\ 0.04\, y_1-3\cdot 10^7\, y_2^2 \\ 3\cdot 10^7\, y_2^2
\end{pmatrix}, &
f(y)=&\ff(y)+\fs(y).
\end{align}
Now, we solve \eqref{eq:robertson} either with the RKC or the mRKC scheme using step sizes $\tau = 1/2^k$, $k=0,\ldots,7$. For comparison, we use a reference solution obtained with the standard fourth-order Runge--Kutta scheme using $\tau=10^{-4}$. In \cref{fig:conv0_err}, we observe that both the RKC and the mRKC method achieve first-order convergence. In fact, both errors are hardly distinguishable, indicating that the error introduced by the approximation of $f$ by $\fe$ is negligible. We observe in \cref{fig:conv0_eta} that the mean value of $\eta$ during integration is indeed considerably smaller than $\tau$.
\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,1)},anchor=north east},log basis y={10},log basis x={2},legend cell align={left},x dir=reverse,
			xlabel={$\tau$}, ylabel={$\ell_\infty$ error},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=dt,y=err,col sep=comma]
			{data/robertson/res_mRKC.csv};\addlegendentry{mRKC}
			\addplot[color=\classicalcolor,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=dt,y=err,col sep=comma]
			{data/robertson/res_RKC.csv};\addlegendentry{RKC}
			\addplot[black,dashed,line width=\plotdashedlinewidth pt,domain=0.0078125:1] (x,0.00025*x);\addlegendentry{$\bigo{\tau}$}
			\end{loglogaxis}
			\end{tikzpicture}
			\caption{$\ell_\infty$-error vs. the step size $\tau$.}
			\label{fig:conv0_err}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,1)},anchor=north east},log basis y={10},log basis x={2},legend cell align={left},x dir=reverse,
			xlabel={$\tau$}, ylabel={$[t]$},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=\markone,mark size=\plotmarksizeu pt] table [x=dt,y=eta,col sep=comma]
			{data/robertson/res_mRKC.csv};\addlegendentry{Mean of $\eta$}
			\addplot[color=\colortwo,line width=\plotlinewidth pt,mark=\marktwo,mark size=\plotmarksized pt] table [x=dt,y=dt,col sep=comma]
			{data/robertson/res_mRKC.csv};\addlegendentry{$\tau$}
			\end{loglogaxis}
			\end{tikzpicture}
			\caption{Comparing the size of $\eta$ and $\tau$.}
			\label{fig:conv0_eta}
		\end{subfigure}
	\caption{Robertson's stiff test problem. Convergence and comparison of $\eta$ against $\tau$.}
	\label{fig:conv0}
\end{figure}

Next, we compare the two schemes for a fixed step size $\tau=1$. In \cref{fig:conv0_stages}, we display the number of stages taken by the mRKC and the RKC method at each time step with respect to $t\in [0,100]$. Moreover, \cref{fig:conv0_rho} depicts the evolution of the spectral radii $\rho,\rhof,\rhos$ of the Jacobians of $f,\ff,\fs$, respectively. We observe that $\rhos$ decreases with time and consequently the mRKC scheme decreases the number $s$ of expensive function evaluations $\fs$ per step.
In contrast, $\rho$ increases and thus the RKC scheme must increase the number $s$ of $\fs$ function evaluations, although this term does not introduce any stiffness; indeed, $\rho$ increases only because of the term contained in $\ff$.
Finally, we notice in \cref{fig:conv0_stages} that the mRKC scheme increases the number $m$ of (cheap) function evaluations $\ff$ because of the increase in $\rhof$ and $\eta$; indeed, $\eta$ also increases due to the decrease in $s$ and \eqref{eq:defsmeta}. This added cost, however, is much smaller than that from the many additional (expensive) evaluations of $\fs$ required by the RKC method.
\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,0)},anchor=south east},log basis y={10},log basis x={2},legend cell align={left},ymin=0.2,ymax=100,
			xlabel={$t$}, ylabel={Stages},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}},grid=none]
			\addplot[color=\localcolorF,solid,line width=\plotlinewidth pt,mark=\localmarkF,mark size=\plotmarksizeu pt,mark repeat=24,mark phase=0] table [x=t,y=m,col sep=comma]
			{data/robertson/rho_sm_mRKC.csv};\addlegendentry{mRKC, $m$.}
			\addplot[color=\localcolorS,solid,line width=\plotlinewidth pt,mark=\localmarkS,mark size=\plotmarksizeu pt,mark repeat=24,mark phase=8] table [x=t,y=s,col sep=comma]
			{data/robertson/rho_sm_mRKC.csv};\addlegendentry{mRKC, $s$.}
			\addplot[color=\classicalcolor,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksizeu pt,mark repeat=24,mark phase=16] table [x=t,y=s,col sep=comma]
			{data/robertson/rho_s_RKC.csv};\addlegendentry{RKC, $s$.}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{Stages needed by RKC and mRKC vs. time $t$, for a fixed step size $\tau=1$.}
			\label{fig:conv0_stages}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,0)},anchor=south west},log basis y={10},legend cell align={left},ymin=100,ymax=10000,
			xlabel={$t$}, ylabel={Spectral radii},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}},grid=none]
			\addplot[color=\localcolorF,solid,line width=\plotlinewidth pt,mark=\localmarkF,mark size=\plotmarksizeu pt,mark repeat=24,mark phase=0] table [x=t,y=rhoF,col sep=comma]
			{data/robertson/rho_sm_mRKC.csv};\addlegendentry{$\rhof$}
			\addplot[color=\localcolorS,solid,line width=\plotlinewidth pt,mark=\localmarkS,mark size=\plotmarksizeu pt,mark repeat=24,mark phase=8] table [x=t,y=rhoS,col sep=comma]
			{data/robertson/rho_sm_mRKC.csv};\addlegendentry{$\rhos$}
			\addplot[color=\classicalcolor,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksizeu pt,mark repeat=24,mark phase=16] table [x=t,y=rho,col sep=comma]
			{data/robertson/rho_s_RKC.csv};\addlegendentry{$\rho$}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{Evolution of the spectral radii $\rho$, $\rhof$, $\rhos$ vs. time $t$.}
			\label{fig:conv0_rho}
		\end{subfigure}
	\caption{Robertson's stiff test problem. Comparison of spectral radii and number of stages taken by the mRKC and the RKC scheme.}
	\label{fig:conv0_rho_stages}
\end{figure}

\subsection{Heat equation in the unit square}\label{sec:exp_conv2}
Next, we verify the space-time convergence properties of the mRKC method. To do so, we consider the heat equation in the unit square $\Omega=[0,1]\times [0,1]$,
\begin{equation}\label{eq:par1}
\begin{aligned}
\dt u-\Delta u &=g \qquad &&\text{in }\Omega\times [0,T],\\
u&=0 &&\text{in } \partial\Omega \times [0,T],\\
u&=0 &&\text{in } \Omega\times\{0\},
\end{aligned}
\end{equation}
where $T=1/2$ and $g$ is chosen such that $u(\bx,t)=\sin(\pi x_1)^2\sin(\pi x_2)^2\sin(\pi t)^2$ is the exact solution.

Starting from a mesh of $2^j\times 2^j$ simplicial elements with $j=2,\ldots,5$, we locally refine twice all the elements inside the square $\OF=(1/4,3/4)\times (1/4,3/4)$. Each refinement step is performed by splitting all edges of any simplex, i.e. every triangle is split into four self-similar children.
Let $\mathcal{M}$ be the set of elements in the mesh and $\mathcal{M}_F=\{T\in\mathcal{M}\,:\, \overline{T}\cap \overline{\Omega}_F\neq \emptyset\}$ the set of refined elements or their direct neighbors. Then $h=H/4$ is the diameter of the elements inside of $\OF$, with $H$ the diameter of the elements outside of $\OF$.

Next, we discretize \eqref{eq:par1} in space with first-order DG-FE \cite{PiE12} on the mesh $\mathcal{M}$. After inverting the block-diagonal mass matrix, the resulting system is
\begin{equation}
y'=A\, y+G,  \qquad\qquad y(0)=y_0,
\end{equation}
where $A\in\Rb^{N\times N}$ and $G\in C([0,T],\Rb^N)$ corresponds to the spatial discretization of $g(\cdot,t)$. Let $D\in\Rb^{N\times N}$ be a diagonal matrix with $D_{ii}=1$ if the $i$th degree of freedom belongs to an element in $\mathcal{M}_F$ and $D_{ii}=0$ otherwise. We also introduce
\begin{align}\label{eq:splitpar}
A_F=&DA, & A_S=&(I-D)A &\mbox{and} &&\ff(y)=&A_F\, y, & \fs(t,y)=&A_S\, y+G(t),
\end{align}
with $I$ the identity. It is well-known that the spectral radii $\rhos$ and $\rhof$ of $A_S$ and $A_F$ behave as $\bigo{1/H^2}$ and $\bigo{1/h^2}=\bigo{16/H^2}$, respectively.

We now consider a sequence of meshes with $j=2,\ldots,5$ and solve \eqref{eq:par1} either with
the mRKC or the RKC scheme using the same step size $\tau=1/2^j$. The parameters $s$ and $m$ for mRKC are chosen according to \eqref{eq:defsmetaweak}. In \cref{fig:conv2}, we display the $H^1(\Omega)$ errors at final time for mRKC and RKC. Both methods yield space-time first-order convergence and result in similar errors. In \cref{fig:stages}, we show the number of stages needed by RKC and mRKC. For both schemes, $s$ increases as the mesh size $H$ decreases, but for mRKC, $s$ is much smaller, since it only depends on the coarse elements, while $m$ remains constant due to the constant ratio between $\rhof$ and $\rhos$.
\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,1)},anchor=north east},log basis y={10},ymax=1.5,legend cell align={left},
			xlabel={$j$}, ylabel={$H^1(\Omega)$ error},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=k,y=errh1,col sep=comma]
			{data/heat_eq/sol_mrkc.csv};\addlegendentry{mRKC}
			\addplot[color=\classicalcolor,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=k,y=errh1,col sep=comma]
			{data/heat_eq/sol_rkc.csv};\addlegendentry{RKC}
			\addplot[black,dashed,line width=\plotdashedlinewidth pt,domain=2:5] (x,2/2^x);\addlegendentry{$\bigo{\tau+H}$}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{Convergence of RKC and mRKC.}
			\label{fig:conv2}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis y={10},ymax=3200,legend cell align={left},
			xlabel={$j$}, ylabel={Stages},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}},grid=none]
			\addplot[color=\localcolorF,solid,line width=\plotlinewidth pt,mark=\localmarkF,mark size=\plotmarksizeu pt] table [x=k,y=m,col sep=comma]
			{data/heat_eq/sol_mrkc.csv};\addlegendentry{mRKC, $m$.}
			\addplot[color=\localcolorS,solid,line width=\plotlinewidth pt,mark=\localmarkS,mark size=\plotmarksizeu pt] table [x=k,y=s,col sep=comma]
			{data/heat_eq/sol_mrkc.csv};\addlegendentry{mRKC, $s$.}
			\addplot[color=\classicalcolor,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksizeu pt] table [x=k,y=s,col sep=comma]
			{data/heat_eq/sol_rkc.csv};\addlegendentry{RKC, $s$.}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{Number of stages needed by RKC and mRKC vs. refinement level $j$.}
			\label{fig:stages}
		\end{subfigure}
	\caption{Heat equation in the unit square. Space-time convergence and number of stages.}
	\label{fig:convstages}
\end{figure}

\subsection{Diffusion across a narrow channel}\label{sec:exp_eff2}
To illustrate the efficiency of the mRKC method in a situation where geometry constraints require local mesh refinement, we consider the heat equation
\begin{equation}\label{eq:par3}
\begin{aligned}
\dt u-\Delta u &=g \qquad &&\text{in }\Omega_\delta\times [0,T],\\
\nabla u\cdot \bm{n}&=0 &&\text{in } \partial\Omega_\delta \times [0,T],\\
u&=0 &&\text{in } \Omega_\delta\times\{0\},
\end{aligned}
\end{equation}
with $T=0.1$ inside $\Omega_\delta$, which consists of two $10\times 5$ rectangles linked by a narrow $\delta\times 0.05$ channel of width $\delta>0$, see \cref{fig:eff2_sol}. The right-hand side $g(\bx,t)=\sin(10\pi t)^2 e^{-5\Vert \bx-\bm{c}\Vert^2}$ corresponds to a smoothed Gaussian point source centered at $\bm{c}$ in the middle of the upper rectangle.
\begin{figure}
		\hfil
		\begin{subfigure}[t]{\subfigsizemed\textwidth}
			\begin{tikzpicture}[spy using outlines= {circle, magnification=6, connect spies}]
			\coordinate (a) at (-2.19,0);
			\coordinate (b) at (1,0);
			\coordinate (spypoint) at (b);
			\coordinate (magnifyglass) at (2.2,-1.2);
			\node at (a) {\includegraphics[width=0.15\textwidth]{images/narrow_channel/bar.png}};
			\node at (b) {\includegraphics[width=0.64\textwidth]{images/narrow_channel/sol_k_3_T_1.png}};
			\spy [WildStrawberry, size=2cm,magnification=6] on (spypoint) in node[fill=white] at (magnifyglass);
			\end{tikzpicture}
			\caption{Solution for $\delta=1/2^2$.}
			\label{fig:eff2_sola}
		\end{subfigure}\hfil
		\begin{subfigure}[t]{\subfigsizemed\textwidth}
			\begin{tikzpicture}[spy using outlines= {circle, magnification=6, connect spies}]
			\coordinate (a) at (-2.19,0);
			\coordinate (b) at (1,0);
			\coordinate (spypoint) at (b);
			\coordinate (magnifyglass) at (2.2,-1.2);
			\node at (a) {\includegraphics[width=0.15\textwidth]{images/narrow_channel/bar.png}};
			\node at (b) {\includegraphics[width=0.64\textwidth]{images/narrow_channel/sol_k_8_T_1.png}};
			\spy [WildStrawberry, size=2cm,magnification=6] on (spypoint) in node[fill=white] at (magnifyglass);
			\end{tikzpicture}
			\caption{Solution for $\delta=1/2^7$.}
			\label{fig:eff2_solb}
		\end{subfigure}
	\hfil
	\caption{Narrow channel. Numerical solutions of \cref{eq:par3} at $t=1$ using mRKC for a channel width $\delta=1/2^2$ or $\delta=1/2^7$.}
	\label{fig:eff2_sol}
\end{figure}

Inside $\Omega_\delta$, we use a Delaunay triangulation with maximal element size $H\approx 0.015$. As $\delta$ approaches zero, the elements inside the channel become increasingly smaller and the system stiffer. For each $\delta>0$, we define a neighborhood $\Omega_{F,\delta}\subset\Omega_\delta$ of the channel and $\mathcal{M}_F$, $A$, $A_F$, $A_S$, $\ff$, $\fs$ as in \cref{sec:exp_conv2}. Here, $\Omega_{F,\delta}$ is chosen such that the spectral radius of $A_S$ is almost independent of $\delta$ and only that of $A_F$ increases with decreasing $\delta$. Hence, $\Omega_{F,\delta}$ contains the channel together with all neighboring elements of mesh size smaller than $H$, see \cref{fig:eff2_grida,fig:eff2_gridb}.
\begin{figure}
		\hfil
		\begin{subfigure}[t]{\subfigsizemed\textwidth}
			\centering
			\includegraphics[trim=0cm 0cm 0cm 0cm, clip, width=0.64\textwidth]{images/narrow_channel/grid_k_3_lw_2.png}
			\caption{Zoom over the channel with $\delta=1/2^2$.}
			\label{fig:eff2_grida}
		\end{subfigure}\hfil
		\begin{subfigure}[t]{\subfigsizemed\textwidth}
			\centering
			\includegraphics[trim=0cm 0cm 0cm 0cm, clip, width=0.64\textwidth]{images/narrow_channel/grid_k_8_lw_2.png}
			\caption{Zoom over the channel with $\delta=1/2^7$.}
			\label{fig:eff2_gridb}
		\end{subfigure}
	\hfil
	\caption{Narrow channel. Zoom of the FE mesh for a channel width $\delta=1/2^2$ or $1/2^7$, with the subdomain $\Omega_{F,\delta}$ (in blue).}
	\label{fig:eff2_grid}
\end{figure}

For varying channel width $\delta=1/2^{k}$, $k=0,\ldots,15$, we now solve \eqref{eq:par3} with the RKC and mRKC method using the choice of parameters \eqref{eq:defsmetaweak} with $\tau=0.01$. In \cref{fig:eff2_speed}, the relative speed-up defined as the ratio between the computational times of RKC and mRKC always exceeds one and reaches a value as high as 40.
Note that the relative error between the two solutions in $L^2(\Omega_\delta)$ or $H^1(\Omega_\delta)$ norm is at most $3\cdot 10^{-4}$, as shown in \cref{fig:eff2_relerr}.

In \cref{fig:eff2_rho}, we display for varying $\delta$ also the spectral radii $\rho,\rhof,\rhos$ of $A,A_F,A_S$, respectively; note that $\rhof$ and $\rho$ essentially coincide. For large $\delta$, we also have $\rhof\approx\rhos$ since the typical element size is sufficiently small to resolve the channel (\cref{fig:eff2_grida}). For $\delta$ small, we observe that $\rho$, $\rhof$ increase as $1/\delta^2$ while $\rhos$ remains almost constant. \Cref{fig:eff2_stages} shows that the number $s$ of stages in the mRKC scheme remains constant, as does $\rhos$ in \cref{fig:eff2_rho}, while $m$ increases (as $\rhof$).
For large $\delta$, we have $\rhof\approx\rhos$ and thus $m=1$; then, the RKC and mRKC schemes coincide. Indeed, as is shown in \cref{fig:eff2_relerr}, for $m=1$ then the relative error between the RKC and mRKC solutions is of the order of machine precision.

In \cref{fig:eff2_comptime}, we observe that for $\delta$ large the CPU times of the two methods are similar; thus, despite $\rhof\approx\rhos$, there is no loss in efficiency and the speed-up is at least one (\cref{fig:eff2_speed}). For moderate values of $\delta$, the cost of RKC increases proportionally to $1/\delta$, while the cost of mRKC is hardly affected. For even smaller $\delta$, the number of evaluations of $\ff$ increases and so does its cost with respect to $\fs$ (see \cref{fig:eff2_cf}), since the number of elements in $\mathcal{M}_F$ increases (\cref{fig:eff2_grid}). In this regime, evaluation of $\ff$ dominates the computational cost of mRKC, which increases linearly in $1/\delta$, too. Still, the mRKC method remains about forty times faster than the classical RKC method for this particular discretization inside $\Omega_\delta$, see \cref{fig:eff2_speed}.

\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogxaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis x={2},x dir=reverse,ytick={1,10,20,30,40,50},ymax=50,
			xlabel={$\delta$}, ylabel={Speed-up},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=delta,y=speedup,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{Speed-up}
			\addplot[black,dashed,line width=\plotdashedlinewidth pt,domain=3.0517578125e-05:1] (x,1);\addlegendentry{$1$}
			\end{semilogxaxis}
			\end{tikzpicture}
			\caption{Relative speed-up of mRKC over RKC.}
			\label{fig:eff2_speed}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis x={2},log basis y={10},x dir=reverse,ylabel style={yshift=-2pt},
			xlabel={$\delta$}, ylabel={CPU time [sec.]},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=delta,y=mRKC_cost,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{Cost mRKC}
			\addplot[color=\classicalcolor,solid,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=delta,y=RKC_cost,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{Cost RKC}
			\addplot[black,dashed,line width=\plotdashedlinewidth pt,domain=3.0517578125e-05:0.0078] (x,1/x/64);\addlegendentry{$\bigo{1/\delta}$}
			\end{loglogaxis}
			\end{tikzpicture}
			\caption{Total CPU time w.r.t. $\delta$.}
			\label{fig:eff2_comptime}
		\end{subfigure}\\ \vspace{0.5cm}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,0)},anchor=south east},log basis x={2},log basis y={10},x dir=reverse, ytick={1e-14,1e-8,1e-4},ymax=1,ylabel style={yshift=-8pt},yticklabel style={xshift=3pt},
			xlabel={$\delta$}, ylabel={Relative error},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=\markone,mark size=\plotmarksizeu pt] table [x=delta,y=relerrl2,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{$L^2(\Omega_\delta)$ norm}
			\addplot[color=\colortwo,solid,line width=\plotlinewidth pt,mark=\marktwo,mark size=\plotmarksizeu pt] table [x=delta,y=relerrh1,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{$H^1(\Omega_\delta)$ norm }
			\end{loglogaxis}
			\end{tikzpicture}
			\caption{RKC and mRKC solutions' relative error $\Vert u^{\mRKCop}-u^{\RKCop}\Vert/\Vert u^{\RKCop}\Vert$.}
			\label{fig:eff2_relerr}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis x={2},log basis y={10},x dir=reverse,legend cell align={left},
			xlabel={$\delta$}, ylabel={Spectral radii},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\localcolorF,solid,line width=\plotlinewidth pt,mark=\localmarkF,mark size=\plotmarksizeu pt] table [x=delta,y=mRKC_rhoF,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{$\rhof$}
			\addplot[color=\localcolorS,solid,line width=\plotlinewidth pt,mark=\localmarkS,mark size=\plotmarksizeu pt] table [x=delta,y=mRKC_rhoS,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{$\rhos$}
			\addplot[color=\classicalcolor,solid,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=delta,y=RKC_rho,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{$\rho$}
			\addplot[black,dashed,line width=\plotdashedlinewidth pt,domain=3.0517578125e-05:0.0078] (x,8/x/x);\addlegendentry{$\bigo{1/\delta^2}$}
			\end{loglogaxis}
			\end{tikzpicture}
			\caption{Spectral radii w.r.t. $\delta$.}
			\label{fig:eff2_rho}
		\end{subfigure}\\ \vspace{0.5cm}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis x={2},log basis y={10},x dir=reverse,legend cell align={left},legend cell align={left},
			xlabel={$\delta$}, ylabel={Stages},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\localcolorF,solid,line width=\plotlinewidth pt,mark=\localmarkF,mark size=\plotmarksizeu pt] table [x=delta,y=mRKC_m,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{mRKC, $m$}
			\addplot[color=\localcolorS,solid,line width=\plotlinewidth pt,mark=\localmarkS,mark size=\plotmarksizeu pt] table [x=delta,y=mRKC_s,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{mRKC, $s$}
			\addplot[color=\classicalcolor,solid,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=delta,y=RKC_s,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{RKC, $s$}
			\end{loglogaxis}
			\end{tikzpicture}
			\caption{Number of stages needed by RKC and mRKC.}
			\label{fig:eff2_stages}
		\end{subfigure} \hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogxaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis x={2},x dir=reverse,y tick label style={/pgf/number format/.cd,fixed,precision=2,/tikz/.cd},
			xlabel={$\delta$}, ylabel={$\costf$},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=delta,y=mRKC_cF,col sep=comma]
			{data/narrow_channel/results.csv};\addlegendentry{Relative cost $\costf$ of $\ff$}
			\end{semilogxaxis}
			\end{tikzpicture}
			\caption{Relative evaluation cost of $\ff(y)$ w.r.t. $\ff+\fs$ as a function of $\delta$.}
			\label{fig:eff2_cf}
		\end{subfigure}
	\caption{Narrow channel. Speed-up, error, spectral radii and stages number w.r.t. channel width $\delta$.}
	\label{fig:eff2_results}
\end{figure}
Finally, we compare mRKC against the second-order accurate version of the RKC scheme (RKC2) from \cite{SSV98,VHS90} and the implicit Euler (IE) method. To do so, we consider the two channel widths $\delta=1/2^6,1/2^{12}$ and solve \cref{eq:par3} with the RKC, RKC2, mRKC and IE schemes. For both values of $\delta$ we display in \cref{fig:WvsP} the computational times against the final error, with varying step size $\tau=T/2^j$, $j=0,\ldots,16$ and $T=0.1$. Since the exact solution is unknown, the final error is computed against a reference solution obtained from the second-order RKC2 scheme with step size $\tau=10^{-6}$. For the wider channel with $\delta=1/2^6$, we observe in \cref{fig:WvsPa} that the mRKC scheme is more efficient than RKC2 for most of the step sizes, but becomes less efficient at higher accuracy. The IE and the mRKC method are about equally efficient. In contrast, for the narrow channel with $\delta=1/2^{12}$, that is, in a situation of even more severe stiffness, we observe in \cref{fig:WvsPb} that mRKC is always much faster than RKC2: the speed-up ranges from ten to seventy times faster, depending on the step size imposed by the desired accuracy. In this case of extreme stiffness, the IE method always remains slightly faster than mRKC; clearly, the IE method is also particularly efficient here thanks to the symmetry and linearity of the discrete Laplacian.
\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
				\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,0)},anchor=south east},log basis x={10},log basis y={10},x dir=reverse,
					ylabel style={yshift=0pt},yticklabel style={xshift=2pt},
					xlabel={Relative error in $L^2(\Omega_\delta)$ norm}, ylabel={CPU time [sec.]},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
					\addplot[color=\classicalcolor,solid,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=relerrl2,y=rkc_step_inc_sub,col sep=comma]
					{data/narrow_channel/RKC1_k_7.csv};\addlegendentry{RKC}
					\addplot[color=\colortwo,solid,line width=\plotlinewidth pt,mark=\marktwo,mark size=\plotmarksizeu pt] table [x=relerrl2,y=rkc_step_inc_sub,col sep=comma]
					{data/narrow_channel/RKC2_k_7.csv};\addlegendentry{RKC2}
					\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=relerrl2,y=rkc_step_inc_sub,col sep=comma]
					{data/narrow_channel/mRKC_k_7.csv};\addlegendentry{mRKC}
					\addplot[color=\colorfour,solid,line width=\plotlinewidth pt,mark=\markfour,mark size=\plotmarksizeu pt] table [x=relerrl2,y=linsystemtot,col sep=comma]
					{data/narrow_channel/IE_k_7.csv};\addlegendentry{IE}
				\end{loglogaxis}
			\end{tikzpicture}
			\caption{$\delta=1/2^6$}
			\label{fig:WvsPa}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
				\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,0)},anchor=south east},log basis x={10},log basis y={10},x dir=reverse,
					ylabel style={yshift=0pt},yticklabel style={xshift=2pt},
					xlabel={Relative error in $L^2(\Omega_\delta)$ norm}, ylabel={CPU time [sec.]},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
					\addplot[color=\classicalcolor,solid,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=relerrl2,y=rkc_step_inc_sub,col sep=comma]
					{data/narrow_channel/RKC1_k_13.csv};\addlegendentry{RKC}
					\addplot[color=\colortwo,solid,line width=\plotlinewidth pt,mark=\marktwo,mark size=\plotmarksizeu pt] table [x=relerrl2,y=rkc_step_inc_sub,col sep=comma]
					{data/narrow_channel/RKC2_k_13.csv};\addlegendentry{RKC2}
					\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=relerrl2,y=rkc_step_inc_sub,col sep=comma]
					{data/narrow_channel/mRKC_k_13.csv};\addlegendentry{mRKC}
					\addplot[color=\colorfour,solid,line width=\plotlinewidth pt,mark=\markfour,mark size=\plotmarksizeu pt] table [x=relerrl2,y=linsystemtot,col sep=comma]
					{data/narrow_channel/IE_k_13.csv};\addlegendentry{IE}
				\end{loglogaxis}
			\end{tikzpicture}
			\caption{$\delta=1/2^{12}$}
			\label{fig:WvsPb}
		\end{subfigure}
	\caption{Narrow channel. Work vs. accuracy diagram of the RKC, RKC2, mRKC and implicit Euler (IE) scheme for two channel widths $\delta$. Symbols represent step sizes $\tau=T/2^j$, for $j=0,\ldots,16$.}
	\label{fig:WvsP}
\end{figure}

\subsection{Integro-differential equation}\label{sec:intdiff}
To compare the mRKC and implicit Euler scheme on a nonsymmetric and truly nonlinear problem, we now consider the one-dimensional integro-differential problem from \cite{VaV92}, also considered in \cite{AbV13,Zbi11}:
\begin{equation}\label{eq:intdiff}
\begin{aligned}
	\dt u(t,x)  &=\Delta u(t,x)-\sigma  \int_0^1 \frac{u(t,s)^4}{(1+|x-s|)^2}\dif s \qquad && (t,x)\in [0,1]\times [0,1],\\
	u(0,x)&= \cos(x\pi/2)^2 && x\in [0,1],\\
	u(t,0) &= 1-\sqrt{t}/2 && t\in[0,1],\\
	\partial_x u(t,1) &= 0 && t\in[0,1],
\end{aligned}
\end{equation}
with $\sigma=0.01$. Problem \cref{eq:intdiff} models an idealized temperature profile of air near the ground. We discretize \cref{eq:intdiff} in space on a uniform grid of $N$ cells using central finite differences  for the Laplacian and the composite trapezoidal rule for the integral term. For the mRKC scheme, we assign the Laplacian to $\ff$ and the integral term to $\fs$. For the implicit Euler method, the Jacobian (evaluated analytically) and its LU decomposition are computed only once per time step following \cite[IV.8]{HaW02}.

For $N=100, 3200$ (varying problem size and degree of stiffness), we apply the RKC, RKC2, mRKC and IE schemes for $\tau=1/2^{j}$ with $j=2,\ldots,12$. For each run, we monitor the computational time and the $L^2([0,1])$ error against a reference solution: the resulting efficiency graphs are shown in \cref{fig:eff_intdiff}.  For $N=100$, IE method is always faster than RKC, which remains the most expensive method. For large $\tau$, IE is comparable to mRKC, but as $\tau$ decreases, the mRKC scheme becomes significantly more efficient. The second-order scheme RKC2 becomes faster than mRKC only at high accuracy. For the stiffer case with $N=3200$, the mRKC scheme is the fastest method, as its efficiency is only marginally affected by the increased stiffness, while the cost of the direct solver in the Newton iteration clearly starts to dominate the overall cost of IE. For those parameter settings, the RKC and RKC2 methods were overly expensive and could not be run to completion.
\begin{figure}
	\begin{subfigure}[t]{\subfigsize\textwidth}
		\centering
		\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,0)},anchor=south east},log basis x={10},log basis y={10},x dir=reverse,
				ylabel style={yshift=0pt},yticklabel style={xshift=2pt},
				xlabel={Error in $L^2([0,1])$ norm}, ylabel={CPU time [sec.]},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
				\addplot[color=\classicalcolor,solid,line width=\plotlinewidth pt,mark=\classicalmark,mark size=\plotmarksized pt] table [x=err,y=cpu,col sep=comma]
				{data/integrodifferential/RKC1_N_100.csv};\addlegendentry{RKC}
				\addplot[color=\colortwo,solid,line width=\plotlinewidth pt,mark=\marktwo,mark size=\plotmarksizeu pt] table [x=err,y=cpu,col sep=comma]
				{data/integrodifferential/RKC2_N_100.csv};\addlegendentry{RKC2}
				\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=err,y=cpu,col sep=comma]
				{data/integrodifferential/mRKC_PAR_N_100.csv};\addlegendentry{mRKC}
				\addplot[color=\colorfour,solid,line width=\plotlinewidth pt,mark=\markfour,mark size=\plotmarksizeu pt] table [x=err,y=cpu,col sep=comma]
				{data/integrodifferential/IE_LU_N_100.csv};\addlegendentry{IE}
			\end{loglogaxis}
		\end{tikzpicture}
		\caption{$N=100$}
		\label{fig:eff_intdiffa}
	\end{subfigure}\hfill
	\begin{subfigure}[t]{\subfigsize\textwidth}
		\centering
		\begin{tikzpicture}[scale=\plotimscale]
			\begin{loglogaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend cell align={left},legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(1,0)},anchor=south east},log basis x={10},log basis y={10},x dir=reverse,
				ylabel style={yshift=0pt},yticklabel style={xshift=2pt},
				xlabel={Error in $L^2([0,1])$ norm}, ylabel={CPU time [sec.]},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
				\addplot[color=\localcolor,solid,line width=\plotlinewidth pt,mark=\localmark,mark size=\plotmarksizeu pt] table [x=err,y=cpu,col sep=comma] {data/integrodifferential/mRKC_PAR_N_3200.csv};\addlegendentry{mRKC}
				\addplot[color=\colorfour,solid,line width=\plotlinewidth pt,mark=\markfour,mark size=\plotmarksizeu pt] table [x=err,y=cpu,col sep=comma] {data/integrodifferential/IE_LU_N_3200.csv};\addlegendentry{IE}
			\end{loglogaxis}
		\end{tikzpicture}
		\caption{$N=3200$}
		\label{fig:eff_intdiffb}
	\end{subfigure}
	\caption{Integro-differential problem. Work vs. accuracy diagram of the RKC, RKC2, mRKC and implicit Euler (IE) scheme for two mesh sizes $1/N$. Symbols represent step sizes $\tau=1/2^j$, for $j=2,\ldots,12$.}
	\label{fig:eff_intdiff}
\end{figure}

\subsection{Reaction-convection-diffusion problem}\label{sec:exp_stab}
In \cref{sec:stabanalysis} we proved that the stability conditions of the mRKC method are the same for the $2\times 2$ model problem \eqref{eq:twodim} and for the scalar multirate test equation \eqref{eq:mtesteq}. The splitting of the discrete Laplace operator in \eqref{eq:splitpar} in fact is similar to that in \eqref{eq:tauAfAc1} for the $2\times 2$ model problem. Thus, one could expect that the stability conditions \eqref{eq:defsmeta} are also necessary for more general parabolic problems. However, spatial discretizations of parabolic problems are much more complex than \eqref{eq:twodim}. Here we shall demonstrate via numerical experiment that the weaker stability conditions \eqref{eq:defsmetaweak} in fact are also necessary and sufficient for general parabolic reaction-convection-diffusion problems, such as
\begin{equation}\label{eq:parmod}
\begin{aligned}
\dt u-\nabla\cdot(K\nabla u)+\bm{\beta}\cdot\nabla u +\mu u &=g \qquad &&\text{in }\Omega\times [0,T],\\
u&=0 &&\text{in } \partial\Omega \times [0,T],\\
u&=u_0 &&\text{in } \Omega\times\{0\}.
\end{aligned}
\end{equation}
These experiments also illustrate that the mRKC method indeed requires no scale separation.

We now consider three distinct parameter regimes. First, we let $\Omega=[0,2]\times [0,1]$, $K=I_{2\times 2}$, $\bm{\beta}=\bm{0}$ and $\mu=0$. Inside $\Omega$, we build a $16\times 8$ uniform mesh and refine twice the elements inside of $\OF=(1,2)\times (0,1)$ (see \cref{fig:stab_test_dom}). Again, we use DG-FE
for the spatial discretization, which yields the two matrices $A_F$ and $A_S$, as described in \cref{sec:exp_conv2}. Next, we set $\tau=1$, $s,m,\eta$ as in \eqref{eq:defsmeta} and $A_\eta$ as in \eqref{eq:defAeta}. One step of the mRKC scheme is given by $y_1=\Pe_s(\tau A_\eta)y_0$. We recall that a necessary condition for stability of the scheme (at least for linear problems) is $\tau \rho_\eta\leq \beta s^2$, where $\rhoe$ is the spectral radius of $A_\eta$.

Let $\overline\beta$ be as in \eqref{eq:defsmetaweak}, $\overline\eta\in [0,\eta]$, $\overline m$ such that $\overline\eta\rhof\leq\overline\beta\overline m^2$,
\begin{equation}
\overline A_\eta=\Ps_{\overline m}(\overline \eta A_F)A
\end{equation}
and $\overline \rho_\eta$ be the spectral radius of $\overline A_\eta$. We wish to study for which $\overline\eta$ it holds $\tau\overline{\rho}_\eta\leq\beta s^2$. In \cref{fig:rho0}, we display $\tau \overline\rho_\eta$ for $\overline\eta\in (0,\eta)$ with respect to $w(\overline{\eta})=-\overline\eta\beta s^2/\tau$: for $|w(\overline \eta)|\geq 2$, it holds $\tau\overline\rho_\eta\leq\beta s^2$ and thus the scheme is stable. Observe that $|w(\overline \eta)|\geq 2$ is equivalent to $\overline \eta \geq 2\tau/(\beta s^2)$, as in \eqref{eq:defsmetaweak}.
\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering

			\begin{tikzpicture}
			\node at (0,0) {\includegraphics[width=0.96\textwidth]{images/exp_stab/domain_white.png}};
			\draw[white,fill=white,opacity=0.8] (1.76,0.) circle (0.4cm);
			\draw[white,fill=white,opacity=0.8] (-1.76,0.) circle (0.4cm);
			\node at (1.79,0.) {{\Large $\Omega_F$}};
			\node at (-1.73,0.) {{\Large $\Omega_S$}};
			\draw[white,fill=white,opacity=0.0] (0,-2.4) circle (0.4cm);
			\end{tikzpicture}
			\caption{Domains $\Omega_S$ and $\Omega_F$}
			\label{fig:stab_test_dom}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis y={2},scaled x ticks=false,
			xtick={-6,-4,-2,0}, xticklabels={$w(\eta)$, $-4$, $-2$, $0$}, ymax=4850,
			xlabel={$w(\overline\eta)$},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=weta,y=rho,col sep=comma]
			{data/num_stab/rho_0.csv};\addlegendentry{$\tau \overline \rho_\eta$}
			\addplot[color=black,dashed,line width=\plotdashedlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=weta,y=bssq, col sep=comma]
			{data/num_stab/rho_0.csv};\addlegendentry{$\beta s^2$}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{$\tau\overline\rho_\eta$ vs. $w(\overline\eta)$, first problem.}
			\label{fig:rho0}
		\end{subfigure}
	\caption{Reaction-convection-diffusion problem. Mesh and spectral radius $\tau \overline{\rho}_\eta$ of the first problem setting.}
	\label{fig:stab_test}
\end{figure}
Since the smallest (in magnitude) nonzero eigenvalues of the discrete Laplacians, $A_F$ and $A_S$,
do not depend on the mesh size, but only depends on the size of the domain,
they essentially coincide; hence, this problem exhibits no scale separation assumption. Nevertheless
the mRKC scheme remains stable, as expected from theory.

Finally, we consider two additional cases that further corroborate the previous findings. First, we set $\Omega=[0,1]\times [0,1]$, $K=I_{2\times 2}$, $\bm{\beta}=(1,1)^\top$ and $\mu=1$.  In  $\Omega$,
we build a $8\times 8$ uniform mesh and refine three times the elements inside the small inner square $\OF=(1/4,3/4)\times (1/4,3/4)$. In \cref{fig:rho1}, we show again $\tau \overline\rho_\eta$ for $\overline\eta\in (0,\eta)$ with respect to $w(\overline\eta)$: for $|w(\overline \eta)|\geq 2$, $\tau\overline\rho_\eta\leq\beta s^2$ holds. Next, we use a uniform $32\times 32$ mesh in $\Omega=[0,1]\times [0,1]$ which is refined twice in the lower left corner $\OF=(0,1/32)\times (0,1/32)$. We also set $\bm{\beta}=0$, $\mu=0$, $K(\bx)=1$ for $x_1\geq x_2$ and $K(\bx)=0.1$ elsewhere. The results, shown in \cref{fig:rho2}, again confirm the stability of the mRKC with parameters chosen according to \eqref{eq:defsmetaweak}.
\begin{figure}
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis y={2},ymax=32768,scaled x ticks=false,
			xtick={-6,-4,-2,0}, xticklabels={$w(\eta)$, $-4$, $-2$, $0$},
			xlabel={$w(\overline\eta)$},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=weta,y=rho,col sep=comma]
			{data/num_stab/rho_1.csv};\addlegendentry{$\tau \overline \rho_\eta$}
			\addplot[color=black,dashed,line width=\plotdashedlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=weta,y=bssq, col sep=comma]
			{data/num_stab/rho_1.csv};\addlegendentry{$\beta s^2$}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{$\tau\overline\rho_\eta$ vs. $w(\overline\eta)$, second problem.}
			\label{fig:rho1}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\begin{semilogyaxis}[height=\aspectratio*\plotimsized\textwidth,width=\plotimsized\textwidth,legend columns=1,legend style={draw=\legendboxdraw,fill=\legendboxfill,at={(0,1)},anchor=north west},log basis y={2},ymax=524288,scaled x ticks=false,
			xtick={-6,-4,-2,0}, xticklabels={$w(\eta)$, $-4$, $-2$, $0$},
			xlabel={$w(\overline\eta)$},label style={font=\normalsize},tick label style={font=\normalsize},legend image post style={scale=\legendmarkscale},legend style={nodes={scale=\legendfontscale, transform shape}}]
			\addplot[color=\colorone,solid,line width=\plotlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=weta,y=rho,col sep=comma]
			{data/num_stab/rho_2.csv};\addlegendentry{$\tau \overline \rho_\eta$}
			\addplot[color=black,dashed,line width=\plotdashedlinewidth pt,mark=none,mark size=\plotmarksizeu pt] table [x=weta,y=bssq, col sep=comma]
			{data/num_stab/rho_2.csv};\addlegendentry{$\beta s^2$}
			\end{semilogyaxis}
			\end{tikzpicture}
			\caption{$\tau\overline\rho_\eta$ vs. $w(\overline\eta)$, third problem.}
			\label{fig:rho2}
		\end{subfigure}
	\caption{Stability experiment. Illustration of $\tau\overline\rho_\eta$ versus $w(\overline\eta)$, second and third problem setting.}
	\label{fig:num_stab}
\end{figure}


TITLE: Interpretation of 21 cm intensity mapping
QUESTION [2 upvotes]: I'm currently trying to understand the following plot from Astrophysics for physicists by A.R. Choudhuri.

I think I understand the basic concept of measuring the line-of-sight velocities with the redshift and the intensities for specific frequencies to generate the plot (mostly due to the help of another SE post (see here)).
What I still don't understand is how exactly it is possible to have positive and negative $v_R$ on the same side (left or right from the $l=0$ line). No matter how we chose the galactic coordinate system exactly (to determine the sign of $l$), in either case the mapped gas will move away from or move towards me, but not both at the same time, right? So how can we measure negative and positive radial velocities at the same time.

REPLY [4 votes]: I suspect what you are forgetting is that the Sun is also in motion around the Galactic centre, not just the gas that is being observed.
The basic equation for determining the line of sight velocity $V$ of a parcel of gas that is moving tangentially (i.e. assuming circular orbits) to the Galactic centre with speed $V_g$ and which is at a distance $r_g$ from the Galactic centre is:
$$ V = \frac{V_g}{r_g} r_0 \sin l - V_0 \sin l,$$
(derivation can be found here, but I note that in your previous question you present the same equation from Choudhuri's book, so presumably it is derived there too) where $r_0$ is the distance from the Sun to the Galactic centre, $V_0$ is the (assumed tangential) velocity of the Sun with respect to the Galactic centre and $l$ is the Galactic longitude ($l=0$ is towards the Galactic centre).
Clearly $V$ can be either positive or negative. For $0<l<90^{\circ}$, $V>0$ if $V_g > V_0 r_g/r_0$ and $V<0$ if $V_g < V_0 r_g/r_0$ (and vice versa for $270<l<360$ - note that usually $l$ is defined to run from 0 to 360).
So, whether $V$ is positive or negative depends on $V_g(r)$. To a zeroth order approximation, the tangential velocity as a function of radius is roughly constant at Galactocentric radii within a few kpc of the Sun's radius - i.e. $V_g \simeq V_0$. In that case we have a positive velocity, for $0<l<90^{\circ}$, if $r_g<r_0$ and negative if $r_g>r_0$ (and vice-versa for $270<l<360^{\circ}$).
When we look in a direction that is within about 45 degrees of the Galactic centre, it is possible to see gas that is inside the solar circle (which is on "our side" of the Galaxy and has $r_g < r_0$), but also to see gas which is beyond the solar circle on the other side of the Galaxy with $r_g>r_0$. Thus it is possible to see gas with both positive and negative line of sight velocities.