# Explicit high-order noncanonical symplectic algorithms for ideal two-fluid systems

## Abstract

An explicit high-order noncanonical symplectic algorithm for ideal
two-fluid systems is developed. The fluid is discretized as particles
in the Lagrangian description, while the electromagnetic fields and
internal energy are treated as discrete differential form fields on
a fixed mesh. With the assistance of Whitney interpolating forms Whitney (1957); Desbrun *et al.* (2008); Xiao *et al.* (2015a),
this scheme preserves the gauge symmetry of the electromagnetic field,
and the pressure field is naturally derived from the discrete internal
energy. The whole system is solved using the Hamiltonian splitting
method discovered by He et al. He *et al.* (2015a),
which was been successfully adopted in constructing symplectic particle-in-cell
schemes Xiao *et al.* (2015a). Because of its structure preserving
and explicit nature, this algorithm is especially suitable for large-scale
simulations for physics problems that are multi-scale and require
long-term fidelity and accuracy. The algorithm is verified via two
tests: studies of the dispersion relation of waves in a two-fluid
plasma system and the oscillating two-stream instability.

###### pacs:

52.65.Rr, 52.25.Dg## I Introduction

The ideal two-fluid model, a basic non-dissipative model of plasma physics, has been widely used to study fusion and astrophysical plasmas. In this model, the electrons and ions are treated as ideal fluids separately, with coupling to the electromagnetic fields through the charge and current carried by them. Although this system is easily generalized to any number of different charged species, the terminology “two-fluid” will be used here in lieu of “multi-fluid”, as is typically done. Because the ideal two-fluid system has noncanonical Hamiltonian form Morrison (1982, 1998), as was shown in Ref. Spencer and Kaufman (1982), its dynamics preserves geometric structure and there is no dissipation of invariants such as the total energy and momentum in the system. Conventional numerical algorithms for the ideal two-fluid system generally do not preserve geometric structure and thus the truncation error can accumulate coherently over simulation time-steps. This is a serious drawback when solving most electron-ion systems whose behaviors are naturally multi-scale. For example, the ion cyclotron period is thousands of times longer than that of the electron.

Symplectic methods, discovered in the 1980s Lee (1982); Ruth (1983); Feng (1985, 1986); Lee (1987); Veselov (1988); Yoshida (1990); Forest and Ruth (1990); Channell and Scovel (1990); Candy and Rozmus (1991); Tang (1993); Sanz-Serna and Calvo (1994); Shang (1999); Marsden and West (2001); Hairer *et al.* (2002); Feng and Qin (2010),
have proven to be efficient for solving finite-dimensional canonical
Hamiltonian systems. Such methods preserve the symplectic geometric
structure (2-form) associated with the original canonical Hamiltonian
system, and the numerical error of all invariants can be globally
bounded by small values throughout simulations Hairer *et al.* (2006).
In plasma physics, accelerator physics, and fluid dynamics, many of
the finite-dimensional Hamiltonian systems and most of the infinite-dimensional
Hamiltonian systems are noncanonical; for example, this is the case
for guiding center dynamics Littlejohn (1979, 1981, 1983),
the Euler fluid and magnetohydrodynamics (MHD) equations Morrison and Greene (1980),
the Vlasov-Maxwell and Vlasov-Poisson systems Morrison (1980, 1982); Weinstein and Morrison (1981); Marsden and Weinstein (1982); Morrison (1998),
and drift and gyrokinetic theories Morrison (2013); Squire *et al.* (2013); Morrison *et al.* (2013); Burby *et al.* (2015); Brizard *et al.* (2016).
The development of geometric algorithms for these systems can be challenging.
However, recently significant advances have been achieved in the development
of structure preserving geometric algorithms for charged particle
dynamics Qin and Guan (2008); Qin *et al.* (2009); Guan *et al.* (2010); Squire *et al.* (2012a); Qin *et al.* (2013); Liu *et al.* (2014); Zhang *et al.* (2014, 2015); Ellison *et al.* (2015a); He *et al.* (2015b, c); Ellison *et al.* (2015b); Liu *et al.* (2015); He *et al.* (2016),
the Vlasov-Maxwell systems Squire *et al.* (2012b, c); Xiao *et al.* (2013); Kraus (2013); Evstatiev and Shadwick (2013); Shadwick *et al.* (2014); Xiao *et al.* (2015b, a); Crouseilles *et al.* (2015); Qin *et al.* (2015); He *et al.* (2015a); Qin *et al.* (2016); Webb (2016),
compressible ideal MHD Zhou *et al.* (2014, 2016),
and incompressible fluids Pavlov *et al.* (2011); Gawlik *et al.* (2011).
All of these methods have demonstrated unparalleled long-term numerical
accuracy and fidelity compared with conventional methods. As a side
note, we point out that for infinite-dimensional Hamiltonian systems,
an alternative viewpoint is to treat them as multi-symplectic systems
Bridges (1997); Marsden *et al.* (1998), and corresponding multi-symplectic
algorithms Reich (2000); Sun and Qin. (2000); Bridges and Reich (2001); Wang and Qin (2001); Hong and Qin (2002); Chen *et al.* (2002); Guo and Wu (2003)
have also been developed.

In the present work, an explicit, high-order, noncanonical symplectic
algorithm for integrating the compressible ideal two-fluid system
is developed. We discretize the fluid as particles in the Lagrangian
description, which naturally guarantees conservation of the density.
The electromagnetic fields and internal energy are discretized over
a cubic mesh by using the theory of discrete exterior calculus (DEC)
Hirani (2003). High-order Whitney interpolating forms
Xiao *et al.* (2015a) are used to ensure the gauge symmetry of
Maxwell’s equations. The discrete Poisson bracket for the ideal two-fluid
system is obtained by the similar technique that is used in obtaining
the discrete Vlasov-Maxwell bracket Xiao *et al.* (2015a), and
the final numerical scheme is constructed by the powerful Hamiltonian
split method He *et al.* (2015a, c); Xiao *et al.* (2015a).
We note that for the existing structure preserving method for the
compressible fluid Zhou *et al.* (2014), all fields are discretized
over a moving mesh, which does not apply to cases where the mesh deforms
significantly during the evolution, such as in a rotating fluid. This
difficulty is overcome by using a fixed mesh rather than a moving
one for discretizing the electromagnetic and internal energy fields
in our method. The conservation of symplectic structure guarantees
that the numerical errors of all invariants such as the total energy
and momentum are bounded within a small value during the simulations
Hairer *et al.* (2006). Therefore, this method is most suitable
for solving long-term multi-scale problems.

The paper is organized as follows. In Sec. II the Hamiltonian theory of the ideal two-fluid system is reviewed and the geometric structure preserving method is developed. Two numerical examples, the dispersion relation of waves in an ideal two-fluid system and the oscillating two-stream instability, are given in Sec. III. Finally, in Sec. IV we conclude.

## Ii Structure preserving discretization for ideal two-fluid systems

The starting point of our development is the Lagrangian of the ideal
two-fluid system, written in terms of Lagrangian variables, which
is quite similar to the Lagrangian for the Vlasov-Maxwell system except
for the addition of internal energy terms (see e.g. Charidakos *et al.* (2014)).
This Lagrangian is given as follows:

(1) | |||||

where , , and are the mass, charge, and initial number density distribution of species , respectively, and are current position and velocity of fluid elements for species labeled by , which we take to be the initial value of in the configuration space, is the Jacobian of the coordinate transformation from the initial value to , is the internal energy per unit mass for species , and is the electromagnetic vector potential. In the arguments of the fields and we suppress the time variable. In this Lagrangian, we have ignored the entropy term in the internal energy, assuming barotropic fluids, and adopted the temporal gauge with . The permittivity and permeability are set to unity for simplicity.

Evolution equations are obtained upon variation of the action as in Hamilton’s principle

(2) |

giving rise the equations of motion via

(3) |

which yield

(4) | |||||

(5) |

where and are the electromagnetic fields. These equations are exactly the ideal two-fluid equations in the Lagrangian variable description.

Now we discretize the Lagrangian using a method very similar to that
for the discretization of the Vlasov-Maxwell system in Ref. Xiao *et al.* (2015a).
The electromagnetic fields and internal energy are sampled over a
cubic mesh, while the fluid is discretized into finite-sized smooth
particle Monaghan (1992); Xiao *et al.* (2013, 2015a)
moving between mesh grids. Modeling fluids using a set of Lagrangian
particles is also the key idea of the smoothed-particle-hydrodynamics
(SPH) method Monaghan (1992); Bonet and Lok (1999); Price and Monaghan (2004).
However, the difference is that our internal energy fields are calculated
on fixed mesh grids. Therefore, the method developed in the present
study more closely resembles the structure-preserving symplectic particle-in-cell
method of Ref. Xiao *et al.* (2015a) The resulting discrete Lagrangian
is

(6) | |||||

where

(7) |

Here, the subscript denotes the -th particle of species
, and are Whitney interpolating
maps for discrete 0-forms and 1-forms Hirani (2003); Whitney (1957); Desbrun *et al.* (2008); Xiao *et al.* (2015a),
is discrete internal energy per unit volume for species ,
is the discrete curl operator that is defined in Eq. (15),
are indices for the discrete 0-form,
1-form, 2-form, respectively. To simplify the notation, the grid size
has been set to unity. The Whitney maps are defined as
follows:

where the one-dimensional interpolation function is chosen in this paper to be

(8) |

The equations of motion arising from the action with Lagrangian of (6) are the following:

(9) | |||||

(10) |

Next we introduce two discrete fields
and , which are discrete electromagnetic
fields. We will make use of the following properties of the interpolating
forms Xiao *et al.* (2015a); Hirani (2003); Whitney (1957),

(11) | |||||

(12) | |||||

(13) | |||||

(14) | |||||

(15) | |||||

(16) | |||||

which hold for any , , and . With these identities, Eqs. (9) and (10) can be expressed as

(17) | |||||

(18) | |||||

(19) |

The continuity equations for the densities are automatically satisfied, as can be shown by directly calculating the time derivative of ,

(20) | |||||

(21) | |||||

(22) | |||||

(23) |

where can be viewed as the discrete momentum density over mesh grids, and is a discrete version of . So, Eq. (23) is essentially a kind of discrete continuity equation.

To construct the geometric structure preserving algorithm, the Hamiltonian
theory for the discretized system is considered. Note that the only
difference between the two-fluid Lagrangian and the Vlasov-Maxwell
Lagrangian is the internal energy term, which can be written as a
function of . Thus, the discrete Poisson structure of
the ideal two-fluid system can be chosen to be the same as that for
the Vlasov-Maxwell system Xiao *et al.* (2015a), which is

And the two-fluid Hamiltonian is

(24) |

It is straightforward to check that the following Hamiltonian equations are identical to Eqs. (17-19),

(25) | |||||

(26) | |||||

(27) | |||||

(28) |

Now the discrete algorithm can be developed. Using a Hamiltonian splitting
technique similar to that in Ref. He *et al.* (2015a); Xiao *et al.* (2015a),
can be split into 6 parts

(29) |

where

(30) | |||||

(31) | |||||

(32) | |||||

(33) |

It turns out that the exact solutions for all sub-systems can be found
and computed explicitly. The exact solutions for , ,
, and have been derived in Ref. Xiao *et al.* (2015a).
They are

(34) | |||||

(35) | |||||

(36) |

The solutions and are similar to . For , the exact evolution equations are

(37) | |||||

(38) | |||||

(39) | |||||

(40) |

Using the property Eq. (11) of the Whitney interpolating forms, the exact solution can be written as

(41) |

Here, can be interpreted as a discrete version of the continuous Newton’s second law,

(42) |

At first look, it seems different from Newton’s law in the Lagrangian form derived in Ref. Morrison (1998), i.e.,

(43) |

This is because the in Eq. (43) is defined as the internal energy per mass. Upon letting , the relation between and is

(44) |

Consequently,

which is identical to the pressure term in the right hand side of Eq. (43). Therefore, the pressure for the species can be defined to be Morrison (1998)

The final geometric structure-preserving scheme can be constructed from these exact solutions. For example, a first-order scheme can be chosen as

(45) |

and a second-order scheme can be constructed as

(46) | |||||

The th-order scheme can be derived from the th-order scheme by using

(47) | |||||

(48) | |||||

(49) |

## Iii Numerical examples

To verify the practicability of our explicit high-order noncanonical symplectic algorithm for ideal two-fluid systems, we apply it to two physics problems. In the first problem we examine the dispersion relation of an electron-deuterium plasma, while the second concerns the oscillation two-stream instability.

For the electron-deuterium plasma, parameters of the unperturbed uniform plasma are chosen as follows,

(50) | |||||

(51) | |||||

(52) | |||||

(53) | |||||

(54) | |||||

(55) | |||||

(56) | |||||

(57) |

where , is the speed of light in the vacuum, is the constant external magnetic field. This plasma supports both electron waves and ion waves, and their frequencies are very different since the deuterium ion is much heavier than the electron. The simulation is carried out in a mesh, and the periodical boundary condition is adopted in all directions. The grid size is choosen to be , the time step is set to . The simulation is initialized with stationary fluid particles being equally spaced with 4 particles per grid cell.

To numerically obtain the dispersion relation, the simulation is carried out with a small random perturbation. The space-time dependence of one field component is transformed into space, and a contour plot the field component in the space is used to make correspondence with the linear dispersion relation of the discrete system. In Fig. 1, such a contour plot is compared with the theoretical dispersion relation in both high frequency and low frequency ranges. We can see that the dispersion relation obtained by our geometric two-fluid algorithm agrees very well with the theory over the frequency range of the simulation. As expected, the total energy of the system is bounded to be within a interval of its intitial value of for all simulation time-steps, which is plotted in Fig. 2.

The second example is the well-known oscillating two-stream instability
Nishikawa (1968); Morales *et al.* (1974). We consider
the case of an unmagnetized cold two-fluid model, and compare with
stability condition that was previously studied in Ref. Qin and Davidson (2014).
We simulate an electron-positron plasma, with system parameters given
as follows:

(58) | |||||

(59) | |||||

(60) |

wiht the relative drift velocity between electrons and positrons chosen to be . The simulation domain is a mesh. Initial perturbations with two different wave numbers, and , are tested. According to the theoretical prediction of Ref. Qin and Davidson (2014), the mode with is stable while that with is unstable. Both of these predictions are confirmed by the simulation using our algorithm, as seen in Fig. 3. The evolution of the perturbed electrostatic field of the unstable mode is plotted in Fig. 4, which displays the space-time dependence of during the nonlinear evolution of the instability.