1 Introduction
In plasma physics, a vast number of particles, interacting through their own fields, generate collective behaviors of tremendous complexity[1]. This complexity provides new uncertainties and unexpected phenomena. To study the complex behavior of a plasma system, pure theory and analytic solutions are not sufficient. Therefore computer simulations have become an essential part of plasma physics research. Not only these simulations must be large in terms of simulation domains but also in terms of simulation time. Even the simulation of an entire plasma discharge, of say 200 seconds, is an eternity on the time scale of small scale turbulence.
To be able to do long and large domain simulations is quite a challenge. Minimization of local errors does not limit global error growth and accumulating errors at each time step, may eventually lead to unphysical solutions. For example in turbulence phenomena, which involve transfer processes between different scales, it might be difficult to distinguish what is physics and what is numerical nonsense. Therefore, because in many cases is impossible to have an aposteriori evaluation of the reliability of the results, it is essential to have absolute confidence on the stability of the numerical algorithm. Towards this goal there are two trends of thought:
1  In the first approach one tries to obtain algorithms with contracting (or at least nonexpanding) evolution operators. If the finite difference explicit operator is uniformly expanding then stability is obtained by making the algorithm implicit. If however the evolution operator has both contracting and expanding sectors the situation is more delicate, a possible approach being the semiimplicit scheme discussed in section 2. One should also bear in mind that in nonlinear systems the spectral nature of the evolution operator changes in time and the algorithm must adapt accordingly.
2  The second approach implements in the numerical scheme exact conservation laws which constrain the evolution. The evolution operators might even be unstable but if enough conservation laws are implemented one is sure that the solution oscillates in a domain between the level hypersurfaces of the conservation laws. Rather than discretizing the continuum equations of motion, the dynamics is directly formulated in terms of a discrete Lagrangian or Hamiltonian, from which the evolution equations are derived. By construction the symmetries and invariants of the Lagrangian are then preserved by the dynamics. This technique has been developed under the name of variational integration and we refer to [2] [3] [4] for further details. However this approach is only applicable if enough conservation laws are identified in the system.
In this paper, in line with trend 1, we develop a semiimplicit scheme that guarantees stability (or marginal stability) for arbitrary dynamical systems. As an illustration we apply it to a twofield model of the plasma in the scrapeoff layer (SOL). This model, not only is a good testing ground for the numerical method as it is important in its own right to model an important region of fusion devices. For more details on the model and the physics underlying it we refer to [6] and references therein.
2 A semiimplicit integration scheme
2.1 Explicit vs. implicit integration schemes
Explicit methods for differential equations compute the solution at a later time from the values at the present and previous times, whereas implicit methods find the solution by solving an equation that involves both the current value and the one at some later time. The extra effort involved in solving the equation is compensated by better numerical stability as well as by the use of larger time steps.
Consider, for example, the numerical solution of a diffusion problem
(1) 
An explicit discretization leads to
(2) 
where is the operator
(3) 
On a lattice of points the
operator has positive eigenvalues
(4) 
The exact time evolution of each one of the eigenvectors
of the operator would be(5) 
a stable evolution because .
However one sees from (2) that not only , but also may be larger than leading to runaway behavior. In particular, this limits the time step to be
In contrast with this situation, an implicit scheme
(6) 
solves as
(7) 
Now is always smaller than one and larger time steps may be used. Now and coincide to first order. An even better approximation is obtained by applying the operator to the average of and which leads to
(8) 
This is a better approximation because
(9) 
coinciding up to order with the Taylor expansion of the exponential. On the other hand is still smaller than one.
All the benefits of the implicit scheme in this case follow from the fact that all eigenvalues are positive. Stability would be lost if some of the eigenvalues of the differential operator were negative. This suggests a splitting of the function space into components associated respectively to the positive and the negative spectrum of the evolution operator.
2.2 A splitting scheme
The example in the previous subsection shows the stabilizing effect of the operator (or more generally of the Laplacian for higher dimensions) for the implicit scheme. This effect persists even with some additional differential terms in the equation. For example for an equation
(10) 
with and arbitrary, the explicit scheme is
(11) 
and the implicit one
(12) 
the operator being
(13)  
In a lattice of points the eigenvalues of this operator are
(14) 
therefore and the implicit scheme is stable.
In contrast, for many other equations of practical interest, in particular when the spectrum has both positive and negative real parts one possible solution would be to split the evolution matrix in
(15) 
into
(16) 
with having eigenvalues smaller or equal to one and having eigenvalues greater or equal to one. Then a semiimplicit scheme would be
(17) 
For this scheme not to be unstable it is necessary that the operators
be nonexpansive. One discusses this question in the context of finitedimensional matrices. Given a norm in , for example the Euclidean norm,
(18) 
is the matrix norm induced by . A matrix is said to be nonexpansive if
(19) 
A matrix is nonexpansive if for the induced norm [5].
Lemma: An arbitrary matrix may always be decomposed into in such a way that both and are nonexpansive.
Proof:
(i) Suppose first that the matrix can be diagonalized with eigenvalues . Then it is always possible to decompose it
in such a way that the eigenvalues and of and satisfy and .
Consider the eigenvalue decomposition
and decompose the diagonal matrix into .
Now, given an eigenvalue of if this eigenvalue is assigned to because .
If one has two possibilities. If it is assigned to . Otherwise decompose
being a real positive number between and to make . This is always possible. Let with and arbitrary. Then
which is if . This last condition may always hold with an in the interval . Therefore the part is assigned to and to because . Finally and .
(ii) If the matrix
is not diagonalizable, for example if it is a singular matrix as it very frequently happens in practical situations, one uses a singular value decomposition
(20) 
where and are orthogonal matrices and is a diagonal matrix with diagonal elements . Notice that the ’s are the eigenvalues of or and, denoting by and the normalized columns of and ,
(21) 
This lemma insures the feasibility of the semiimplicit scheme. From the proof it is also clear that the decomposition is not necessarily unique and may, to some extent, still be adjusted to other numerical conveniences. Notice also that and would be contractive if .
As an illustration, in the following section the semiimplicit scheme will be used to develop a numerical code for a model of the scrapeoff layer of tokamak plasmas.
3 Example: A twofield model of the SOL
A numerical code was developed for the conservative part of a twofield model of the SOL [6], namely
(22) 
The two fields and (logdensity and potential) are defined in a bounded two dimensional space and . The full model contains also loss, diffusion and source terms, but it is the conservative part that is responsible for the main dynamical features of the model [6]. The system (22) is conservative in the sense that it is a noncanonical Hamiltonian system with Hamiltonian
(23) 
and Poisson structure
(24)  
the equations (22) being obtained from by
(25) 
However, no invariants other than the Hamiltonian are explicitly known, which hinders the applicability of variational techniques to the numerical simulation. Nevertheless, we also impose energy conservation on our semiimplicit code.
For the numerical code we choose a grid of points. We have chosen periodic boundary conditions on both and but other boundary conditions may as easily be chosen. The quantities and are
dimensional vectors which evolve according to
(26) 
and being the dependent operators
(27) 
The matrices that must be split, if they are expansive, are
(28) 
Computed for a typical time step along the evolution, Fig.1 shows the singular values obtained on the singular value decomposition of, respectively, , , and . One sees that since both and have contractive and expansive singular values, neither an explicit nor an implicit scheme would be stable. In contrast all the evolution operators involved in the semiimplicit method are nonexpansive. They are not strictly contractive because has some zero singular values. The iteration then becomes
(29) 
the potential being obtained from by a standard Poisson code with the physically irrelevant constraint
(30) 
The splitting guarantees that all evolution directions are either stable or marginally stable. We also correct for some numerical drift on energy by a small change on the magnitude of and at each evolution step. This is equivalent (but simpler) to the general scheme proposed by Shadwick, Bowan and Morrison [7] [8] who add to the dynamical equations extra terms to compensate for the numerical drift of the constants of motion.
We have run the code for a number of different initial conditions and the algorithm shows great stability for large times. Notice that the size of the step affects precision but not stability because it is included in the matrix that is split. An important problem when analyzing solutions of models for physical phenomena is the stability of the solutions under small random perturbations, because that is what one expects to be the real world situation. Of course to judge the effect of the perturbations on the solutions it is important to have a stable simulation algorithm without which it would be difficult to distinguish between the effect of the perturbations and numerical instabilities. Our stable algorithm is particularly adequate to deal with this problem and, in particular, to test whether in conservative complex differential systems there are preferential large scale structures which remain stable under stable perturbations. The role of selection of large scale structures by invariance under random perturbation of invariant measures and solutions has been explored in the Euler context elsewhere [9]. Here we explore the possible emergence of this phenomenon for the twofield model by starting from a periodic exact solution of (22), namely
(31) 
with and , .
From this as an initial condition (at ) we make it evolve with, at each time step, the addition of a symmetrical random perturbation. The result of one such simulation is shown in Fig.2 which shows snapshots of the evolution of the logdensity taken at 1, 50, 100, 200, 300 and 400 evolution steps with .
Without random perturbations the pattern shown in the upper left corner simply translates in time. This is what follows from the exact solution (31), also checked by the application of our algorithm. With the random perturbation the pattern eventually becomes unstable and converges to a concentrated pattern shown in the last plot, which is then very stable for a long time. Inspired by the theory developed in [9], we might conjecture that this type of approximate stochastic stable patterns are those compatible with a unique stochastically stable invariant measure which depends on the boundary conditions. The emergence of these quasistable large scale structures associated to large density concentrations are of great relevance for plasma physics [6].
References
 [1] J. F. Freidberg; Plasma Physics and Fusion Energy, Cambridge University Press, 2010.
 [2] J. E. Marsden and M. West; Discrete mechanics and variational integrators, Acta Numerica 10 (2001) 357514.
 [3] J. E. Marsden, G. W. Patrick and S. Shkoller; Multisymplectic geometry, variational integrators and nonlinear PDE’s, Comm. Math. Phys. 199 (1998) 351395.
 [4] M. Kraus; Variational integrators in plasma physics, Ph. D. Thesis T. U. München, arXiv:13.5665.
 [5] T.C. Lim; Nonexpansive Matrices with Applications to Solutions of Linear Systems by Fixed Point Iterations, Fixed Point Theory and Applications, vol. 2010, Article ID 821928.

[6]
R. Vilela Mendes and J. P. Bizarro;
Analytical study of growth estimates, control and structures in a twofield model of the scrapeoff layer
, Phys. of Plasmas 24 (2017) 012303.  [7] B. A. Shadwick, J. C. Bowan and P. J. Morrison; Exactly conservative integrators, SIAM J. Appl. Math. 59 (1999) 1112  1133..
 [8] P. J. Morrison; Structure and structurepreserving algorithms for plasma physics, Phys. of Plasmas 24 (2017) 055502.
 [9] F. Cipriano, H. Ouerdiane, R. Vilela Mendes; Stochastic stability of invariant measures: The 2D Euler equation, arXiv:1802.06422.