# Decoupled, linear, unconditionally energy stable and charge-conservative
finite element method for a inductionless magnetohydrodynamic phase-field
model

Xiaorong Wang  LSEC, NCMIS, Academy of Mathematics and Systems Science,
Chinese Academy of Sciences; School of Mathematical Science, University of
Chinese Academy of Sciences, Beijing 100190, China.  Xiaodi Zhang  Henan
Academy of Big Data, Zhengzhou University, Zhengzhou 450001, China.
<EMAIL_ADDRESS>(X. Zhang<EMAIL_ADDRESS>(X. Wang)

###### Abstract

In this paper, we consider the numerical approximation for a diffuse interface
model of the two-phase incompressible inductionless magnetohydrodynamics
problem. This model consists of Cahn-Hilliard equations, Navier-Stokes
equations and Poisson equation. We propose a linear and decoupled finite
element method to solve this highly nonlinear and multi-physics system. For
the time variable, the discretization is a combination of first order Euler
semi-implicit scheme, several first order stabilization terms and implicit-
explicit treatments for coupling terms. For the space variables, we adopt the
finite element discretization, especially, we approximate the current density
and electric potential by inf-sup stable face-volume mixed finite element
pairs. With these techniques, the scheme only involves a sequence of decoupled
linear equations to solve at each time step. We show that the scheme is
provably mass-conservative, charge-conservative and unconditionally energy
stable. Numerical experiments are performed to illustrate the features,
accuracy and efficiency of the proposed scheme.

###### keywords:

Inductionless MHD equations, Cahn-Hilliard equation, Mixed finite element
method, Decouple scheme, Energy stable, Charge-conservative

††journal: J. COMPUT. APPL. MATH.

##  1 Introduction

In recent years, the phase field method, also called the diffuse interface
method, has been widely used to simulate the motion of multiphase
incompressible immiscible fluids for both numerical and theoretical research
[1, 2, 3], such as in fields like material sciences, fracture mechanics and
fluid mechanics. Unlike the traditional sharp interface model, the phase field
method describes the interface by a balance of molecular forces in a very thin
layer rather than a free curve that evolves over time. To indicate different
phases, an auxiliary function $\phi$ called phase field or order parameter is
introduced [4, 5]. The phase function takes a distinct constant value in each
bulk phase, and varies smoothly in the interfacial region. In the diffuse
interface theory, the surface motion of the order parameter is driven by a
gradient flow. There are two popular gradient flow-type governing equations:
Allen-Cahn equations [6] and Cahn-Hilliard equations [7]. Comparing with the
classical sharp interface model, one principal advantage of phase field method
is its ability to capture the interface implicitly and automatically, and
allow for topological changes such as self-intersection, pinch-off and
splitting. Thus, phase field method has become an attractive technique, for
extensive study on the phase field approach, we refer to the recent reviews
[8, 9] and the references therein.

Magnetohydrodynamic (MHD) equations are used to describe the behavior of
electrically conductive fluids under the influence of magnetic fields and
electric currents, which couple the Navier-Stokes equations and Maxwell
equations. The theoretical analysis and numerical simulation of incompressible
MHD equations is an area of research currently undergoing intense study [10,
11, 12]. However, in most industrial and laboratory cases, the magnetic
Reynolds number is small, then the induced magnetic field can be neglected
compared with the applied magnetic field. In these cases, Maxwell equations
are replaced by Poisson equation for the electric scalar potential, which
yields the inductionless magnetohydrodynamic (IMHD) equations [13, 14, 15]. It
is well-known that numerical simulation of multiphase flow is a challenging
problem in computational fluid dynamics. This problem becomes more difficult
for MHD due to the interaction between multiple physical fields. The
multiphase MHD flow often is involved in a wide range of scientific and
industrial problems, such as astrophysics, aluminum electrolysis, liquid metal
magnetic pumps, MHD power generators, fusion reactor blankets [16, 17]. Thus
developing fast and effective algorithms for incompressible multiphase MHD
equations has great theoretical significance and applicable value.

In this paper, we focus on the two-phase incompressible IMHD problem in a
general Lipschitz domain, which is frequently performed to simulate the
movement in aluminum electrolysis cell [11, 18] and explore the theoretical
study and numerical experiments of liquid lithium-lead cladding [19].
Traditionally, most of the existing models for multiphase incompressible IMHD
problem are devoted to sharp interface models, for instance, the front-
tracking method [20, 21], the volume-of-fluid method [22, 23, 24], the level-
set method [25, 26]. There are few works have focused on the two-phase
incompressible IMHD equations by diffuse-interface method in the literature.
In 2014, Ding et al. [27] presented a two-phase incompressible inductionless
magnetohydrodynamics model and simulated the deformation of melt interface in
an aluminum electrolytic cell. In 2020, Chen and his collaborators [28]
proposed a linear, decoupled, unconditionally energy stable and second order
time-marching schemes to simulate the two-phase incompressible Cahn-Hilliard-
IMHD (CHIMHD) conducting flow. Zhang [29] derived a diffuse interface model
for IMHD fluids systematically and analyzed sharp interface limits for
different choice of the mobilities. Mao et al. [30] analyzed the fully
discrete finite element approximation of a three-dimensional CHIMHD model and
proved the well-posedness of weak solution to the phase field model by using
the classical compactness method.

Studying efficient numerical methods for solving incompressible CHIMHD flows
is the main focus of this article. The CHIMHD model is a complex system that
involves three coupled physical processes: the phase field model (Cahn-
Hilliard equation), fluid dynamics (Navier-Stokes equations), and electric
field (Possion equation). Due to the highly nonlinear with multi-physics
fields coupled in this system, it is relatively complicated and time-consuming
to solve it by nonlinear and coupled scheme [28, 30]. Our primary goal is to
propose a decoupled and efficient fully discrete mixed finite element method
to solve the system. The time discretization is an Euler semi-implicit scheme
with some extra first order stabilization terms in Cahn-Hilliard equations and
Poisson equation, and implicit-explicit treatments for coupling terms. The
subtle time splitting technique allows us to fully decouple system into three
processes at each time step. For spacial discretization, we adopt classical
inf-sup stable finite element to discretize the phase field and chemical
potential, standard velocity-pressure stable finite element to discretize the
hydrodynamic unknowns, and stable face-volume finite element pairs to
discretize the current density and electric potential. With these techniques,
we obtain an efficient fully discrete scheme, which is decoupled, linear while
still mass-conservative, charge-conservative and unconditionally energy
stable. Numerical experiments are performed to demonstrate the features,
accuracy and efficiency of the proposed scheme.

We emphasize that the proposed scheme is charge-conservative and
unconditionally energy stable. The charge conservation law, namely,
$\rm{div}\boldsymbol{J}=0$ is a physical law in electromagnetics, which plays
an important role in keeping the calculation accuracy for the simulation of
MHD fluid. It was believed that only when the numerical schemes numerically
maintains the physical conservation laws, such as momentum conservation and
charge conservation, can we get accurate results for MHD flows at high
Hartmann number [31, 15]. Due to the rapid changes near the interface, the
non-compliance of energy dissipation laws of the scheme may lead to spurious
numerical solutions [32, 33]. Thus, the design of unconditionally energy
stable schemes are of great importance for solving phase field models. To the
best of our knowledge, the scheme presented in this paper is the first
decoupled, linear and unconditionally energy stable scheme for a phase field
model of two-phase inductionless MHD flows.

The rest of this paper is organized as follows. In section 2, we present the
CHIMHD model, derive a weak formulation and the dissipative energy law. In
section 3, we propose a fully discrete charge-conservation finite element
method and prove its unconditional energy stability. In section 4, we present
some numerical simulations to validate our schemes. Some concluding remarks
are given in section 5.

##  2 The diffuse interface model

In this section, we first introduce the CHIMHD model, then present a weak
formulation, and finally show the dissipative energy law in the PDE level.

###  2.1 Model system

We consider a phase-field model for mixture of two immiscible, incompressible
and conducting fluids with the same density in a bounded domain $\Omega$ with
Lipschitz-continuous boundary $\Sigma\coloneqq\partial\Omega$ in
$\mathbb{R}^{{\rm d}},{\rm d}=2,3$. Following [30, 29], the dynamic behavior
of the flow is governed by the coupling model of the Cahn-Hilliard equations
and the IMHD equations. In this paper, we consider the incompressible CHIMHD
system as follows:

$\displaystyle\rho\left(\boldsymbol{u}_{t}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}\right)-2\nabla\cdot\left(\eta\left(\varphi\right)D\left(\boldsymbol{u}\right)\right)+\nabla
p+\varphi\nabla\mu$ $\displaystyle=\boldsymbol{J}\times\boldsymbol{B}\text{ in
}\Omega,$ (1) $\displaystyle\mathrm{div}\boldsymbol{u}$ $\displaystyle=0\text{
in }\Omega,$ (2)
$\displaystyle\sigma\left(\varphi\right)^{-1}\boldsymbol{J}+\nabla\phi-\boldsymbol{u}\times\boldsymbol{B}$
$\displaystyle=\boldsymbol{0}\text{ in }\Omega,$ (3)
$\displaystyle\mathrm{div}\boldsymbol{J}$ $\displaystyle=0\text{ in }\Omega,$
(4) $\displaystyle\varphi_{t}+{\rm div}(\varphi\boldsymbol{u})-M\Delta\mu$
$\displaystyle=0\text{ in }\Omega,$ (5)
$\displaystyle-\gamma\varepsilon\Delta\varphi+\frac{\gamma}{\varepsilon}f(\varphi)-\mu$
$\displaystyle=0\text{ in }\Omega.$ (6)

where $\rho$ is the density of fluids, $\boldsymbol{u},p$ denote the velocity
and pressure of fluids, $\boldsymbol{J},\phi$ represent the current density
and electric scalar potential, $\varphi,\mu$ denote the phase function and
chemical potential, and $\boldsymbol{B}$ is the applied magnetic field which
is assumed to be given.
$D(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}})$ is
the symmetric gradient tensor. In the Cahn-Hilliard equations (5)-(6),
$f(\varphi)=F^{\prime}(\varphi)$, with $F(\varphi)$ is the double-well
potential, $\gamma$ is the surface tension coefficient, and $\epsilon$ is the
interface thickness, $M$ represents the diffusional mobility related to the
relaxation time scale. In the IMHD equations (1)-(4), the viscosity of the
fluids $\eta(\varphi),$ and the electric conductivity $\sigma(\varphi)$ are
depending on the phase function $\varphi$. They are Lipschitz-continuous
functions of $\varphi$ satisfying

$\begin{array}[]{l}0<\min\left\\{\eta_{1},\eta_{2}\right\\}\leq\eta(\varphi)\leq\max\left\\{\eta_{1},\eta_{2}\right\\},\\\
0<\min\left\\{\sigma_{1},\sigma_{2}\right\\}\leq\sigma(\varphi)\leq\max\left\\{\sigma_{1},\sigma_{2}\right\\},\end{array}$

where $\eta_{i}$ and $\sigma_{i}$ (i=1,2) are the viscosity and electric
conductivity of the fluid on each phase. The system of equations are
supplemented with the following initial values and boundary conditions

$\begin{array}[]{r}\boldsymbol{u}(0)=\boldsymbol{u}_{0},\quad\varphi(0)=\varphi_{0}\quad\text{
in }\Omega,\\\
\boldsymbol{u}=0,\quad\boldsymbol{J}\cdot\boldsymbol{n}=0,\quad\partial_{\boldsymbol{n}}\varphi=\partial_{\boldsymbol{n}}\mu=0\quad\text{
on }\Sigma,\end{array}$

where $\boldsymbol{n}$ is the outer unit normal of $\Sigma$. In this paper, we
consider two-phase flows with matching density $\rho_{1}=\rho_{2}$. Without
lose of generality, we set $\rho\equiv 1$ in the above system (2.1).

######  Remark 2.1.

The first term on the right-hand side of (1),
$\boldsymbol{J}\times\boldsymbol{B}$ is the Lorentz force. The last term on
the left-hand side of (1) the continuum surface tension force in the potential
form, which denotes the phase introduced force. This term appears differently
in literature [34, 9]: $\nabla\cdot(\nabla\phi\otimes\nabla\phi),$
$-\mu\nabla\phi$ or $\phi\nabla\mu.$ It can be shown that these three
expressions are equivalent by redefining the pressure $p$.

######  Remark 2.2.

Though the model considered here is a two-phase model with matched density,
one can still employ this model with Boussinesq approximation to model the
effect of density difference by a gravitational force in case of small density
ratio. The large density ratio case is reserved for future work.

In this paper, we consider the truncated double well potential $F(\varphi)$,

$F(\varphi)=\frac{1}{4}\left\\{\begin{array}[]{ll}4(\varphi+1)^{2}&\text{ if
}\varphi<-1,\\\ \left(\varphi^{2}-1\right)^{2}&\text{ if }-1\leq\varphi\leq
1,\\\ 4(\varphi-1)^{2}&\text{ if }\varphi>1.\end{array}\right.$

The original definition of the potential is logarithmic [7], which guarantees
the phase field $\varphi$ stays within (-1,1). The Ginzburg-Landau potential
$F_{\rm{GL}}(\varphi)=\frac{1}{4}(1-\varphi^{2})^{2}$ is an popular
alternative to approximate the one suggested by Cahn and Hilliard. Following
[34, 35], we further restrict the growth of the potential to quadratic away
from the range [-1,1]. It is proved by [36] that the truncated $F(\varphi)$
also ensures the boundness of $\varphi$ in the Cahn-Hilliard equation.
Therefore, it is a common practice to consider the Cahn-Hilliard equation with
the truncated double-well potential $F(\varphi)$ [37, 34]. It is clear that
the second derivative of $F(\varphi)$ is continuous and bounded,

$L:=\max_{\varphi\in\mathbb{R}}\left|F^{\prime\prime}(\varphi)\right|=2.$

This property is of great help in handling the nonlinear double-well potential
by using stabilization method. For the treatment of the general double-well
potential, we give some comments in Remark 3.7.

###  2.2 Weak formulation and energy estimate

Firstly, we introduce some useful notations and Sobolev spaces. Let
$L^{2}(\Omega)$ be the space of square-integrable functions that equipped with
the inner product and norm:

$(f,g):=\int_{\Omega}fg\,{\rm d}x,\quad\|f\|_{0}:=(f,f)^{1/2},\quad\forall
f,g\in L^{2}(\Omega).$

Its subspace with mean zero over $\Omega$ is written as $L_{0}^{2}(\Omega)$.
Let $H^{1}(\Omega),\boldsymbol{H}({\rm div},\Omega)$ be the subspaces of
$L^{2}(\Omega)$ with square integrable gradients and square integrable
divergences respectively. The equipped norm in $\boldsymbol{H}({\rm
div},\Omega)$ is defined by

$\left\|\boldsymbol{J}\right\|_{{\rm
div},\Omega}=\left(\|\boldsymbol{J}\|_{0,\Omega}^{2}+\|{\rm
div}\boldsymbol{J}\|_{0,\Omega}^{2}\right)^{1/2}.$

Their subspaces with vanishing traces and vanishing normal traces on $\Sigma$
are denoted by $H_{0}^{1}(\Omega),\boldsymbol{H}_{0}(\mathrm{div},\Omega)$. We
refer to [38] for their definitions and inner products.

For convenience, we introduce some notations for function spaces

$\displaystyle\boldsymbol{V}$
$\displaystyle\coloneqq\boldsymbol{H}_{0}^{1}(\Omega),\quad X\coloneqq
H^{1}(\Omega),\quad Q=L_{0}^{2}(\Omega),$ $\displaystyle\boldsymbol{D}$
$\displaystyle\coloneqq\boldsymbol{H}_{0}(\mathrm{div},\Omega),\quad
S=L_{0}^{2}(\Omega).$

To derive the weak formulation of (2.1), we define a trilinear form:

$\begin{array}[]{l}O(\boldsymbol{w},\boldsymbol{u},\boldsymbol{v})=\frac{1}{2}\left(\boldsymbol{w}\nabla\boldsymbol{u},\boldsymbol{v}\right)-\frac{1}{2}\left(\boldsymbol{w}\nabla\boldsymbol{v},\boldsymbol{u}\right),\end{array}$

for any $\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}\in\boldsymbol{V}$.

Armed with the above notation, a weak formulation of the system (2.1) amounts
to find $(\boldsymbol{u},p,\boldsymbol{J},\phi,\varphi,\mu)\in$
$\boldsymbol{V}\times Q\times\boldsymbol{D}\times S\times X\times X$ such that

$\displaystyle\left\langle\boldsymbol{u}_{t},\boldsymbol{v}\right\rangle+2(\eta(\varphi)D(\boldsymbol{u}),D(\boldsymbol{v}))+\mathcal{O}(\boldsymbol{u},\boldsymbol{u},\boldsymbol{v})-(p,{\rm
div}\boldsymbol{v})+(\varphi\nabla\mu,\boldsymbol{v})-\left(\boldsymbol{J}\times\boldsymbol{B},\boldsymbol{v}\right)$
$\displaystyle=\boldsymbol{0},$ (7) $\displaystyle({\rm div}\boldsymbol{u},q)$
$\displaystyle=0,$
$\displaystyle\left(\sigma(\varphi)^{-1}\boldsymbol{J},\boldsymbol{K}\right)-\left(\boldsymbol{u}\times\boldsymbol{B},\boldsymbol{K}\right)-(\phi,{\rm
div}\boldsymbol{K})$ $\displaystyle=\boldsymbol{0},$ $\displaystyle({\rm
div}\boldsymbol{J},\chi)$ $\displaystyle=0,$
$\displaystyle\left\langle\varphi_{t},\psi\right\rangle-(\varphi\boldsymbol{u},\nabla\psi)+M(\nabla\mu,\nabla\psi)$
$\displaystyle=0,$
$\displaystyle\lambda\varepsilon(\nabla\varphi,\nabla\chi)+\frac{\lambda}{\varepsilon}(f(\varphi),\chi)-(\mu,\chi)$
$\displaystyle=0,$

for all
$(\boldsymbol{v},q,\boldsymbol{K},\theta,\psi,\chi)\in\boldsymbol{V}\times
Q\times\boldsymbol{D}\times S\times X\times X$.

Now, we are in a position to establish the energy estimate for the CHIMHD
system.

######  Theorem 2.3.

Let $(\boldsymbol{u},p,\boldsymbol{J},\phi,\varphi,\mu)$ be the solution of
(7). The following energy dissipation law holds

$\mathrm{E}(t)+\int_{0}^{t}\mathrm{P}(s)\mathrm{d}s=\mathrm{E}(0)$

where

$\displaystyle\mathrm{E}(t)$
$\displaystyle:=\int_{\Omega}\left(\frac{1}{2}|\boldsymbol{u}|^{2}+\frac{\lambda\epsilon}{2}|\nabla\varphi|^{2}+\frac{\lambda}{\epsilon}F(\varphi)\right){\rm
d}x,$ $\displaystyle\mathrm{P}(t)$
$\displaystyle:=M\|\nabla\mu\|_{0}^{2}+2\|\sqrt{\eta(\varphi)}D(\boldsymbol{u})\|_{0}^{2}+\left\|\sqrt{\sigma(\varphi)^{-1}}\boldsymbol{J}\right\|_{0}^{2}.$

###### Proof.

Taking
$(\boldsymbol{v},q,\boldsymbol{K},\theta,\psi,\chi)=(\boldsymbol{u},p,\boldsymbol{J},\phi,\mu,\varphi_{t})$
in (7), we have

$\begin{array}[]{l}\left(\boldsymbol{u}_{t},\boldsymbol{u}\right)+2\|\sqrt{\eta(\varphi)}D(\boldsymbol{u})\|_{0}^{2}-(\boldsymbol{J}\times\boldsymbol{B},\boldsymbol{u})+(\varphi\nabla\mu,\boldsymbol{u})=0,\\\
\left(\sigma(\varphi)^{-1}\boldsymbol{J},\boldsymbol{J}\right)-(\boldsymbol{u}\times\boldsymbol{B},\boldsymbol{J})=0,\\\
\left(\varphi_{t},\mu\right)-(\varphi\boldsymbol{u},\nabla\mu)+M\|\nabla\mu\|_{0}^{2}=0,\\\
\lambda\varepsilon(\nabla\varphi,\nabla\varphi_{t})+\frac{\lambda}{\varepsilon}(f(\varphi),\varphi_{t})-(\mu,\varphi_{t})=0.\end{array}$

Combining the these equalities, we obtain

$\left(\boldsymbol{u}_{t},\boldsymbol{u}\right)+\lambda\varepsilon(\nabla\varphi,\nabla\varphi_{t})+\frac{\lambda}{\varepsilon}\left(f(\varphi),\varphi_{t}\right)+2\|\sqrt{\eta(\varphi)}D(\boldsymbol{u})\|_{0}^{2}+\left\|\sqrt{\sigma(\varphi)^{-1}}\boldsymbol{J}\right\|_{0}^{2}+M\|\nabla\mu\|_{0}^{2}=0.$

Integrating both sides over $(0,t)$ yields the theorem. ∎

The energy law describes the evolution of the total energy caused by energy
conversion. Since the induced magnetic field can be neglected and the electric
field is considered to be quasi-static, the total energy $\mathrm{E}$ only
consists of the fluid kinetic energy $\frac{1}{2}\|u\|_{0}^{2}$ and the Cahn-
Hilliard free energy. The dissipation of $\mathrm{E}$ stems from the friction
losses $\|\sqrt{\eta(\varphi)}D(\boldsymbol{u})\|_{0}^{2}$, the Ohmic losses
$\left\|\sqrt{\sigma(\varphi)^{-1}}\boldsymbol{J}\right\|_{0}^{2}$ and the
diffusion transport term $M\|\nabla\mu\|_{0}^{2}$.

##  3 Decoupled energy stable finite element method

In this section, we propose a decoupled, energy stable, mixed finite element
scheme for continuous problem (7).

Let $\mathcal{T}_{h}$ be a shape-regular simplex subdivision of $\Omega$. As
usual, we introduce the local mesh size $h_{K}=\mathrm{diam}\left(K\right)$
and global mesh size $h:=\underset{K\in\mathcal{T}_{h}}{\max}h_{K}$. Here we
choose conforming finite element space pairs
$\left(\boldsymbol{V}_{h},Q_{h}\right)\subset\left(\boldsymbol{V},Q\right)$ to
discretize velocity $\boldsymbol{u}$ and pressure $p$,
$\left(\boldsymbol{D}_{h},S_{h}\right)\subset\left(\boldsymbol{D},S\right)$ to
approximate current density $\boldsymbol{J}$ and electric potential $\varphi$,
and $\left(X_{h},X_{h}\right)$ to discretize the phase field function $\phi$
and chemical potential $\mu$. In addition, we assume these spaces satisfy the
following inf-sup conditions.

######  Proposition 3.4 (inf-sup condition).

The finite element pairs $\left(\boldsymbol{V}_{h},Q_{h}\right)$,
$\left(\boldsymbol{D}_{h},S_{h}\right)$ and $\left(X_{h},X_{h}\right)$ satisfy
the following uniform inf-sup conditions:

$\displaystyle\inf_{q_{h}\in
Q_{h}}\sup_{\boldsymbol{v}_{h}\in\boldsymbol{V}_{h}}\frac{(q_{h},{\rm
div}\boldsymbol{v}_{h})}{\left\|\nabla\boldsymbol{v}_{h}\right\|_{0}\left\|q_{h}\right\|_{0}}$
$\displaystyle\geq\beta_{s},$ (8) $\displaystyle\inf_{\theta_{h}\in
S_{h}}\sup_{\boldsymbol{K}_{h}\in\boldsymbol{D}_{h}}\frac{(\psi_{h},{\rm
div}\boldsymbol{K}_{h})}{\left\|\boldsymbol{K}_{h}\right\|_{{\rm{}_{div}}}\left\|\theta_{h}\right\|_{0}}$
$\displaystyle\geq\beta_{m},$ (9) $\displaystyle\inf_{\chi_{h}\in
X_{h}}\sup_{\psi_{h}\in
X_{h}}\frac{(\nabla\chi_{h},\nabla\psi_{h})}{\left\|\chi_{h}\right\|_{1}\left\|\psi_{h}\right\|_{1}}$
$\displaystyle\geq\beta_{c},$ (10)

where $\beta_{s}$, $\beta_{m}$ and $\beta_{s}$ only depend on $\Omega$.

With these discrete spaces, the semi-discretization formulation of the system
(7) is to find
$(\boldsymbol{u}_{h},p_{h},\boldsymbol{J}_{h},\phi_{h},\varphi_{h},\mu_{h})\in\boldsymbol{V}_{h}\times
Q_{h}\times\boldsymbol{D}_{h}\times S_{h}\times X_{h}\times X_{h}$ such that

$\displaystyle\left((\boldsymbol{u}_{h})_{t},\boldsymbol{v}_{h}\right)+2(\eta(\varphi_{h})D(\boldsymbol{u}_{h}),D(\boldsymbol{v}_{h}))+\mathcal{O}(\boldsymbol{u}_{h},\boldsymbol{u}_{h},\boldsymbol{v}_{h})-(p_{h},{\rm
div}\boldsymbol{v}_{h})$ (11)
$\displaystyle+(\varphi_{h}\nabla\mu_{h},\boldsymbol{v}_{h})-\left(\boldsymbol{J}_{h}\times\boldsymbol{B}_{h},\boldsymbol{v}_{h}\right)$
$\displaystyle=\boldsymbol{0},$ $\displaystyle({\rm
div}\boldsymbol{u}_{h},q_{h})$ $\displaystyle=0,$
$\displaystyle\left(\sigma(\varphi_{h})^{-1}\boldsymbol{J}_{h},\boldsymbol{K}_{h}\right)-\left(\boldsymbol{u}_{h}\times\boldsymbol{B}_{h},\boldsymbol{K}_{h}\right)-(\phi_{h},{\rm
div}\boldsymbol{K}_{h})$ $\displaystyle=\boldsymbol{0},$ $\displaystyle({\rm
div}\boldsymbol{J}_{h},\chi_{h})$ $\displaystyle=0,$
$\displaystyle\left((\varphi_{h})_{t},\psi_{h}\right)-(\varphi_{h}\boldsymbol{u}_{h},\nabla\psi_{h})+M(\nabla\mu_{h},\nabla\psi_{h})$
$\displaystyle=0,$
$\displaystyle\lambda\varepsilon(\nabla\varphi_{h},\nabla\chi_{h})+\frac{\lambda}{\varepsilon}(f(\varphi_{h}),\chi_{h})-(\mu_{h},\chi_{h})$
$\displaystyle=0,$

for all
$(\boldsymbol{v}_{h},q_{h},\boldsymbol{K}_{h},\theta_{h},\psi_{h},\chi_{h})\in\boldsymbol{V}_{h}\times
Q_{h}\times\boldsymbol{D}_{h}\times S_{h}\times X_{h}\times X_{h}$.

With similar arguments in Theorem 2.3, one can easily get the energy stability
of the semi-discrete scheme, thus the details are omitted here.

######  Theorem 3.5.

Let $(\boldsymbol{u},p,\boldsymbol{J},\phi,\varphi,\mu)$ be the solution of
(7). The energy dissipation law holds

$\mathrm{E}(t)+\int_{0}^{t}\mathrm{P}(s)\mathrm{d}s=\mathrm{E}(0)$

where

$\displaystyle\mathrm{E}(t)$
$\displaystyle:=\int_{\Omega}\left(\frac{1}{2}|\boldsymbol{u}_{h}|^{2}+\frac{\lambda\epsilon}{2}|\nabla\varphi_{h}|^{2}+\frac{\lambda}{\epsilon}F(\varphi_{h})\right){\rm
d}x,$ $\displaystyle\mathrm{P}(t)$
$\displaystyle:=M\|\nabla\mu_{h}\|_{0}^{2}+2\|\sqrt{\eta(\varphi_{h})}D(\boldsymbol{u}_{h})\|_{0}^{2}+\left\|\sqrt{\sigma(\varphi_{h})^{-1}}\boldsymbol{J}_{h}\right\|_{0}^{2}.$

Let $\left\\{t_{n}=n\tau:\,\,n=0,1,\cdots,N\right\\},\tau=T/N,$ be an
equidistant partition of the time interval $[0,T].$ For any time dependent
function $\omega\left(x,t\right)$, the full discrete approximation to
$\omega\left(t_{n}\right)$ will be denoted by $\omega_{h}^{n}$. A fully
discrete mixed finite element scheme for problem (7) reads as follows:

Given the initial guess $\boldsymbol{u}^{0}$ and $\phi^{0}$, we compute

$\left(\boldsymbol{u}_{h}^{n+1},p_{h}^{n+1},\boldsymbol{J}_{h}^{n+1},\varphi_{h}^{n+1},\phi_{h}^{n+1},\mu_{h}^{n+1}\right),\,n=0,1,\cdots,N-1,$

by the following three steps.
Step 1: Find $\left(\varphi_{h}^{n+1},\mu_{h}^{n+1}\right)\in X_{h}\times
X_{h}$ such that

$\begin{cases}\left(\delta_{t}\varphi_{h}^{n+1},\psi_{h}\right)+M\left(\nabla\mu_{h}^{n+1},\nabla\psi_{h}\right)+\tau\left(\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\varphi_{h}^{n}\nabla\psi_{h}\right)&=\left(\varphi_{h}^{n}\boldsymbol{u}_{h}^{n},\nabla\psi_{h}\right),\\\
\lambda\varepsilon\left(\nabla\varphi_{h}^{n+1},\nabla\chi_{h}\right)+\left(\frac{\lambda}{\epsilon}(\varphi_{h}^{n+1}-\varphi_{h}^{n}),\chi_{h}\right)\\\
+\left(\frac{\lambda}{\varepsilon}f(\varphi^{n}),\chi_{h}\right)-\left(\mu_{h}^{n+1},\chi_{h}\right)&=0,\end{cases}$
(12)

for all $(\psi_{h},\chi_{h})\in X_{h}\times X_{h}$.
Step 2: Find
$\left(\boldsymbol{J}_{h}^{n+1},\phi_{h}^{n+1}\right)\in\boldsymbol{D}_{h}\times
S_{h}$ such that

$\begin{cases}\left(\sigma\left(\varphi_{h}^{n+1}\right)^{-1}\boldsymbol{J}_{h}^{n+1},\boldsymbol{K}_{h}\right)+\tau\left(\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B},\boldsymbol{K}_{h}^{n+1}\times\boldsymbol{B}\right)\\\
-\left(\phi_{h}^{n+1},\mathrm{div}\boldsymbol{K}_{h}\right)-\tau\left(\varphi_{h}^{n+1}\nabla\mu_{h}^{n+1}\times\boldsymbol{B},\boldsymbol{K}_{h}^{n+1}\right)&=\left(\boldsymbol{u}_{h}^{n}\times\boldsymbol{B},\boldsymbol{K}_{h}^{n+1}\right),\\\
\left(\mathrm{div}\boldsymbol{J}_{h}^{n+1},\theta_{h}\right)&=0,\end{cases}$
(13)

for all $(\boldsymbol{K}_{h},\theta_{h})\in\boldsymbol{D}_{h}\times S_{h}$.
Step 3: Find
$\left(\boldsymbol{u}^{n+1},p^{n+1}\right)\in\boldsymbol{V}_{h}\times Q_{h}$
such that

$\begin{cases}(\delta_{t}\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h})+O(\boldsymbol{u}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h})-2\left(\eta\left(\varphi_{h}^{n+1}\right)D\left(\boldsymbol{u}_{h}^{n+1}\right),D\left(\boldsymbol{v}_{h}^{n+1}\right)\right)\\\
\qquad-\left(p_{h}^{n+1},{\rm
div}\boldsymbol{v}_{h}\right)+\left(\varphi_{h}^{n+1}\nabla\mu_{h}^{n+1},\boldsymbol{v}_{h}\right)+\left(\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B},\boldsymbol{v}_{h}\right)&=\boldsymbol{0},\\\
\left(\mathrm{div}\boldsymbol{u}^{n+1},q_{h}\right)&=0,\end{cases}$ (14)

for all $(\boldsymbol{v}_{h},q_{h})\in\boldsymbol{V}_{h}\times Q_{h}$.

Before we get into the discussion on the properties of this scheme, several
remarks about the scheme are given in order.

######  Remark 3.6.

To decouple the nonlinear coupled multiphysics system, we introduce two first
-order stabilization terms. The first stabilization term
$\tau(\varphi^{n}_{h}\nabla\mu^{n}_{h},\varphi^{n}_{h}\nabla\psi_{h})$ in step
1 is to decouple Cahn–Hilliard equations and Navier–Stokes equations [35]. The
second stabilization term
$\tau\left(\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B},\boldsymbol{K}_{h}\times\boldsymbol{B}\right)$
in step 2 is to decouple Poisson equation and Navier–Stokes equations [39].
These extra two stabilization terms are vital to keep the couping term
explicitly while preserving the energy stability, see the proof of Theorem
3.5.

######  Remark 3.7.

In step 1, we employ the stabilized method [35] to treat the double-well
potential $F(\varphi)$ explicitly without suffering from any time step
constraint. Note that this stabilizing term
$\frac{\lambda}{\epsilon}(\varphi^{n+1}_{h}-\varphi^{n}_{h})$ introduces an
extra consistent error of order $\tau$, which is the same order as the overall
truncation error of the scheme. Note that there are many other efficient
methods on constructing energy stable schemes for the Cahn-Hilliard equation,
such as, convex-splitting method [40], invariant energy quadratization method
[41], and scalar auxiliary variable method [42]. Here we adopt the stabilized
explicit method only for simplicity, and it is easy to extend the scheme to
the methods mentioned above.

######  Remark 3.8.

In this paper, we only focus on the decoupling of multiphysics problem rather
the decoupling of all variable in each physical problem. The velocity and
pressure in step 3 can be further decoupled by using the first order pressure
correction scheme [34], and we leave it to the interested readers.

It is clear that the scheme given by (12)-(14) is a decoupled, linear scheme.
Next we want to show that the scheme is mass-conservative and charge-
conservative.

######  Proposition 3.9.

Let
$\left(\boldsymbol{u}_{h}^{m},p_{h}^{m},\boldsymbol{J}_{h}^{m},\phi_{h}^{m},\varphi_{h}^{m},\mu_{h}^{m}\right)$
solve (12)-(14) for any $1\leq m\leq N$, then the scheme is mass-conservative
and charge-conservative, namely,

$\int_{\Omega}\varphi_{h}^{m}\mathrm{\leavevmode\nobreak\
d}x=\int_{\Omega}\varphi_{h}^{0}\mathrm{\leavevmode\nobreak\ d}x,\quad{\rm
div}\boldsymbol{J}_{h}^{m}=0.$

###### Proof.

Letting $\psi_{h}=1$ in the first equation of (12), we have mass conservation
$\int_{\Omega}\varphi_{h}^{m}\mathrm{\leavevmode\nobreak\
d}x=\int_{\Omega}\varphi_{h}^{0}\mathrm{\leavevmode\nobreak\ d}x$. Then, we
note that for all $\theta_{h}\in S_{h}$, there holds

$\left(\theta_{h},{\rm div}\boldsymbol{J}_{h}^{m}\right)=0,$

and ${\rm div}\boldsymbol{J}_{h}^{m}\in S_{h}$. Taking $\theta_{h}={\rm
div}\boldsymbol{J}_{h}^{m}$, we obtain ${\rm div}\boldsymbol{J}_{h}^{m}=0$. ∎

Now, we are in a position to prove the unconditional energy stability of the
decoupled scheme as follows, which is analogous to that of the original
problem in Theorem 2.3.

######  Theorem 3.10.

The decoupled scheme is unconditionally energy stable in the sense that the
following energy estimate

$\delta_{t}\mathrm{E}^{n+1}+\mathrm{P}^{n}\leq 0\quad\forall n\geq 0,$ (15)

holds, where

$\displaystyle\mathrm{E}^{n+1}$
$\displaystyle:=\frac{1}{2}\|\boldsymbol{u}_{h}^{n+1}\|_{0}^{2}+\frac{\lambda\epsilon}{2}\|\nabla\varphi_{h}^{n+1}\|_{0}^{2}+\frac{\lambda}{\epsilon}\left(F\left(\varphi_{h}^{n+1}\right),1\right),$
$\displaystyle\mathrm{P}^{n}$
$\displaystyle:=M\|\nabla\mu_{h}^{n+1}\|_{0}^{2}+2\|\sqrt{\eta(\varphi_{h}^{n+1})}D(\boldsymbol{u}_{h}^{n+1})\|_{0}^{2}+\left\|\sqrt{\sigma(\varphi_{h}^{n+1})^{-1}}\boldsymbol{J}_{h}^{n+1}\right\|_{0}^{2}.$

###### Proof.

Letting
$\left(\psi_{h},\chi_{h}\right)=\left(\mu_{h}^{n+1},\delta_{t}\varphi_{h}^{n+1}\right)$
in (12), we have

$\displaystyle\left(\delta_{t}\varphi_{h}^{n+1},\mu_{h}^{n+1}\right)+\tau\|\varphi_{h}^{n}\nabla\mu_{h}^{n+1}\|^{2}+M\|\nabla\mu_{h}^{n+1}\|^{2}$
$\displaystyle=(\varphi_{h}^{n}\boldsymbol{u}_{h}^{n},\nabla\mu_{h}^{n+1}),$
(16)
$\displaystyle\lambda\varepsilon(\nabla\varphi_{h}^{n+1},\nabla\delta_{t}\varphi_{h}^{n+1})+\frac{\lambda}{\varepsilon}(f(\varphi_{h}^{n}),\delta_{t}\varphi_{h}^{n+1})+\frac{1}{\epsilon\tau}\|\varphi_{h}^{n+1}-\varphi_{h}^{n}\|^{2}$
$\displaystyle=(\mu_{h}^{n+1},\delta_{t}\varphi_{h}^{n+1}).$ (17)

Next, setting
$\left(\boldsymbol{K}_{h},\theta_{h}\right)=\left(\boldsymbol{J}_{h}^{n+1},\phi_{h}^{n+1}\right)$,
we get

$\left\|\sqrt{\sigma(\varphi_{h}^{n+1})^{-1}}\boldsymbol{J}_{h}^{n+1}\right\|_{0}^{2}+\tau\left\|\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}\right\|_{0,\Omega}^{2}=-(\boldsymbol{u}_{h}^{n},\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B})+\tau(\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}).$
(18)

Then, taking
$\left(\boldsymbol{v}_{h},q_{h}\right)=\left(\boldsymbol{u}_{h}^{n+1},p_{h}^{n+1}\right)$,
and using the identity

$\left(\delta_{t}\boldsymbol{u}_{h}^{n+1},\boldsymbol{u}_{h}^{n+1}\right)=\frac{1}{2\tau}\left(\left\|\boldsymbol{u}_{h}^{n+1}\right\|_{0,\Omega}^{2}-\left\|\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}+\left\|\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}\right),$

it yields

$\displaystyle\frac{1}{2\tau}\left(\left\|\boldsymbol{u}_{h}^{n+1}\right\|_{0,\Omega}^{2}-\left\|\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}+\left\|\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}\right)+2\|\sqrt{\eta(\varphi_{h}^{n+1})}D(\boldsymbol{u}_{h}^{n+1})\|_{0}^{2}$
$\displaystyle\qquad=\left(\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h},\boldsymbol{u}_{h}^{n+1}\right)-\left(\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\boldsymbol{u}_{h}^{n+1}\right).$
(19)

By making the summations of (16)-(19), we obtain

$\displaystyle\frac{1}{2\tau}\left(\left\|\boldsymbol{u}_{h}^{n+1}\right\|_{0,\Omega}^{2}-\left\|\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}+\left\|\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}\right)+2\|\sqrt{\eta(\varphi_{h}^{n+1})}D(\boldsymbol{u}_{h}^{n+1})\|_{0}^{2}$
$\displaystyle+\left\|\sqrt{\sigma(\varphi_{h}^{n+1})^{-1}}\boldsymbol{J}_{h}^{n+1}\right\|_{0}^{2}+\tau\left\|\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}\right\|_{0,\Omega}^{2}+\tau\|\varphi_{h}^{n}\nabla\mu_{h}^{n+1}\|^{2}+M\|\nabla\mu_{h}^{n+1}\|^{2}$
$\displaystyle=\left(\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h}-\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right)$
$\displaystyle\quad+\tau(\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h})$
$\displaystyle\quad-\frac{\lambda}{\varepsilon}(f(\varphi_{h}^{n}),\delta_{t}\varphi_{h}^{n+1})-\frac{1}{\epsilon\tau}\|\varphi_{h}^{n+1}-\varphi_{h}^{n}\|^{2}.$
(20)

Using Young inequality, we derive the right hand side of the first term of
(20) has the following estimate,

$\displaystyle\left(\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h}-\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right)$
$\displaystyle\leq\left\|\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h}-\varphi_{h}^{n}\nabla\mu_{h}^{n+1}\right\|_{0,\Omega}\left\|\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}$
$\displaystyle\leq\frac{\tau}{2}\left\|\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h}-\varphi_{h}^{n}\nabla\mu_{h}^{n+1}\right\|_{0,\Omega}^{2}+\frac{1}{2\tau}\left\|\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}$
$\displaystyle=\frac{\tau}{2}\left\|\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h}\right\|_{0,\Omega}^{2}-\tau(\varphi_{h}^{n}\nabla\mu_{h}^{n+1},\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}_{h})$
$\displaystyle\quad+\frac{\tau}{2}\left\|\varphi_{h}^{n}\nabla\mu_{h}^{n+1}\right\|_{0,\Omega}^{2}+\frac{1}{2\tau}\left\|\boldsymbol{u}_{h}^{n+1}-\boldsymbol{u}_{h}^{n}\right\|_{0,\Omega}^{2}.$
(21)

For the last two term in (20), using the Taylor expansion,

$F\left(\varphi_{h}^{n+1}\right)-F\left(\varphi_{h}^{n}\right)=\left(f\left(\varphi_{h}^{n}\right),\varphi_{h}^{n+1}-\varphi_{h}^{n}\right)+\frac{f^{\prime}\left(\xi_{h}^{n}\right)}{2}\|\varphi_{h}^{n+1}-\varphi_{h}^{n}\|^{2},$

then we have

$\displaystyle-\frac{\lambda}{\varepsilon}(f(\varphi_{h}^{n}),\delta_{t}\varphi_{h}^{n+1})-\frac{\lambda}{\epsilon\tau}\|\varphi_{h}^{n+1}-\varphi_{h}^{n}\|^{2}$
$\displaystyle=-\frac{\lambda}{\epsilon\tau}\left(F\left(\varphi_{h}^{n+1}\right)-F\left(\varphi_{h}^{n}\right)\right)+\frac{\lambda}{\epsilon\tau}\left(\frac{f^{\prime}\left(\xi_{h}^{n}\right)}{2}-1\right)\|\varphi_{h}^{n+1}-\varphi_{h}^{n}\|^{2}.$
(22)

Plugging (21) and (22) into (20), and dropping the positive term
$\frac{\tau}{2}\left\|\boldsymbol{J}_{h}^{n+1}\times\boldsymbol{B}\right\|_{0,\Omega}^{2}$
and $\frac{\tau}{2}\|\varphi_{h}^{n}\nabla\mu_{h}^{n+1}\|^{2}$, we get the
required estimate (15). The proof is thus complete. ∎

##  4 Numerical Experiments

In this section, we present a series of 2D numerical experiments to illustrate
the features of the proposed algorithms. The finite element method is
implemented on the finite element software FreeFEM developed by [43]. For any
integer $k\geq 0,$ let $P_{k}(K)$ be the space of polynomials of degree $k$ on
element $K$, and denote $\boldsymbol{P}_{k}(K)=P_{k}(K)^{2}$. We employ the
Mini-element [44] to approximate the velocity and pressure

$\boldsymbol{V}_{h}=\boldsymbol{P}_{1,h}^{b}\cap\boldsymbol{V},\quad
Q_{h}=\left\\{q_{h}\in Q:\left.q_{h}\right|_{K}\in P_{1}(K),\,\forall
K\in\mathcal{T}_{h}\right\\},$

where $P_{1,h}^{b}=\left\\{v_{h}\in C^{0}(\Omega):\left.v_{h}\right|_{K}\in
P_{1}(K)\oplus{\rm span}\\{\hat{b}\\},\,\forall K\subset T_{h}\right\\}$,
$\hat{b}$ is a bubble function. We choose the lowest-order Raviart-Thomas [45]
element space

$\boldsymbol{D}_{h}=\left\\{\boldsymbol{K}_{h}\in\boldsymbol{D}:\left.\boldsymbol{K}_{h}\right|_{K}\in\boldsymbol{P}_{0}(K)+\boldsymbol{x}P_{0}(K),\,\forall
K\in\mathcal{T}_{h}\right\\},$

to approximate the current density, the discontinuous and piecewise constant
finite element space to approximate the electric potential

$S_{h}=\left\\{\psi_{h}\in S:\left.\psi_{h}\right|_{K}\in P_{0}(K),\,\forall
K\in\mathcal{T}_{h}\right\\}.$

The phase field $\varphi$ and chemical potential $\mu$ are discretized by
first order Lagrange finite element space $\left(X_{h},X_{h}\right)$, where

$X_{h}=\left\\{\chi_{h}\in L_{0}^{2}:\left.\chi_{h}\right|_{K}\in
P_{1}(K),\,\forall K\in\mathcal{T}_{h}\right\\}.$

######  Example 4.11 (Convergence and accuracy).

The first example is used to verify the convergence rates in both time and
space. The computational domain is set as $\Omega=\left(0,1\right){}^{2}$, and
the external magnetic field is $\boldsymbol{B}=\left(0,0,1\right)^{{\rm T}}$.
The physical parameters are given by $R_{e}=\kappa=1$ with terminal time
$T=1$. The right-hand side functions, initial conditions and Dirichlet
boundary conditions are chosen such that the given solutions satisfy the
system.

Let the approximation errors at the final time $t=T$ be denoted by

$e_{\omega}=\omega(T)-\omega_{h}^{N}\quad\omega\in\left\\{\boldsymbol{u},p,\boldsymbol{J},\phi,\varphi,\mu\right\\}.$

First, we test the temporal convergence orders. The analytic solutions are
chosen as

$\displaystyle\boldsymbol{u}$
$\displaystyle=\left(y\exp\left(-t\right),x\cos\left(t\right)\right),\quad
p=\sin\left(t\right),$ $\displaystyle\boldsymbol{J}$
$\displaystyle=\left(\sin\left(t\right),\cos\left(t\right)\right),\quad\phi=1,$
$\displaystyle\varphi$ $\displaystyle=(x+y)\exp\left(-t\right),\quad\mu=x\cos
t.$

Note that the exact solutions are linear or constant in space, the main error
comes from the discretization of the time variable. We fix a mesh size with
$h=1/10$ and test the convergence rate with respect to the time step. Then the
errors and orders are displayed in Tables 3-3. From these tables, we observe
that the errors of all variable decrease as the mesh is refined, with
convergence order of $O(\tau)$, which accords with our theoretical analysis
completely.

Table 1: Time convergence rates of the scheme for
$\left(\boldsymbol{u},p\right)$

$\tau$ | $\|e_{\boldsymbol{u}}\|_{0}$ | $\|\nabla e_{\boldsymbol{u}}\|_{0}$ | $\|e_{p}\|_{0}$
---|---|---|---
0.2 | 1.7227e-4(—) | 1.3384e-3(—) | 3.4903e-2(—)
0.1 | 7.5918e-05(1.18) | 5.9037e-4(1.18) | 1.6788e-2(1.06)
0.05 | 3.5535e-05(1.10) | 2.7649e-4(1.09) | 8.1861e-3(1.04)
0.025 | 1.7206e-05(1.05) | 1.3391e-4(1.05) | 4.0362e-3(1.02)
0.0125 | 8.4693e-06(1.02) | 6.5926e-05(1.02) | 2.0033e-3(1.01)
0.00625 | 4.2022e-06(1.01) | 3.2712e-05(1.01) | 9.9790e-4(1.01)

Table 2: Time convergence rates of the scheme for
$\left(\boldsymbol{J},\phi\right)$

$\tau$ | $\|e_{\boldsymbol{J}}\|_{\text{div }}$ | $\|e_{\phi}\|_{0}$ | $\|\nabla\cdot\boldsymbol{J}_{h}^{N}\|_{0}$
---|---|---|---
0.2 | 6.9951e-3(—) | 3.6822e-2(—) | 3.68218e-12
0.1 | 3.6099e-4(0.95) | 1.9026e-2(0.95) | 1.90247e-12
0.05 | 1.8227e-3(0.99) | 9.6627e-3(0.98) | 9.66331e-13
0.025 | 9.1439e-4(1.00) | 4.8671e-3(0.99) | 4.86857e-13
0.0125 | 4.5781e-4(1.00) | 2.4422e-3(0.99) | 2.44362e-13
0.00625 | 2.2904e-4(1.00) | 1.2233e-3(1.00) | 1.22111e-13

Table 3: Time convergence rates of the scheme for $\left(\varphi,\mu\right)$

$\tau$ | $\|e_{\varphi}\|_{0}$ | $\|\nabla e_{\varphi}\|_{0}$ | $\|e_{\mu}\|_{0}$ | $\|\nabla e_{\mu}\|_{0}$
---|---|---|---|---
0.2 | 1.1298e-1(—) | 1.8314e-2(—) | 1.1941e-1 (—) | 6.1937e-3(—)
0.1 | 5.1171e-2(1.14) | 5.9434e-3(1.62) | 7.0687e-2(0.77) | 3.2132e-3(0.95)
0.05 | 2.4156e-2(1.08) | 2.2442e-3(1.40) | 3.6871e-2(0.94) | 1.5619e-3(1.04)
0.025 | 1.1710e-2(1.04) | 9.6111e-4(1.22) | 1.8674e-2(0.98) | 7.6263e-4(1.03)
0.0125 | 5.7621e-3(1.02) | 4.4380e-4(1.11) | 9.3800e-3(0.99) | 3.7607e-4(1.02)
0.00625 | 2.8576e-3(1.01) | 2.1321e-4(1.06) | 4.6987e-3(1.00) | 1.8666e-5(1.01)

Next, we aim to check the spatial approximation orders. To this end, we choose
the exact solution

$\displaystyle\boldsymbol{u}$
$\displaystyle=\left(\sin\left(y\right)\exp\left(-t\right),x^{2}\cos\left(t\right)\right),\quad
p=y\sin\left(t\right),$ $\displaystyle\boldsymbol{J}$
$\displaystyle=\left(y^{2}\sin\left(t\right),\sin\left(x\right)\cos\left(t\right)\right),\quad\phi=x\exp\left(-t\right),$
$\displaystyle\varphi$
$\displaystyle=\sin(x)\exp(-t),\quad\mu=\cos(y)\cos(t).$

With initial time step and mesh width $h_{0}=2\tau_{0}=1/2$, we simultaneously
refine time and space size such that the relation $h=2\tau$ holds. The
corresponding convergent results are demonstrated in Tables 5-5 and a first
order convergence of the proposed numerical scheme can be observed
asymptotically, which agrees with our expected results.

Table 4: Full discretization convergence rates of the scheme for
$\left(\boldsymbol{u},p,\boldsymbol{J}\right)$

$\left(\tau,h\right)$ | $\|\nabla e_{\boldsymbol{u}}\|_{0}$ | $\|e_{p}\|_{0}$ | $\|e_{\boldsymbol{J}}\|_{\text{div }}$ | $\|\nabla\cdot\boldsymbol{J}_{h}^{N}\|_{0}$
---|---|---|---|---
$\left(\tau_{0},h_{0}\right)$ | 1.1685e-1(—) | 5.4282e-2(—) | 1.5322e-1(—) | 1.02906e-11
$\left(\tau_{0},h_{0}\right)/2$ | 5.7591e-2(1.02) | 1.8775e-2(1.53) | 7.7659e-2(0.98) | 1.03095e-11
$\left(\tau_{0},h_{0}\right)/4$ | 2.8635e-2(1.01) | 7.4194e-3(1.34) | 3.9010e-2(0.99) | 1.04248e-11
$\left(\tau_{0},h_{0}\right)/8$ | 1.4287e-2(1.00) | 3.3507e-3(1.15) | 1.9536e-2(1.00) | 1.05131e-11
$\left(\tau_{0},h_{0}\right)/16$ | 7.1370e-3(1.00) | 1.6213e-3(1.05) | 9.7733e-3(1.00) | 1.05643e-11
$\left(\tau_{0},h_{0}\right)/32$ | 3.5671e-3(1.00) | 8.0384e-4(1.01) | 4.8877e-3(1.00) | 1.05915e-11

Table 5: Full discretization convergence rates of the scheme for
$\left(\phi,\varphi,\text{$\mu$}\right)$

$\left(\tau,h\right)$ | $\|e_{\phi}\|_{0}$ | $\|\nabla e_{\varphi}\|_{0}$ | $\|\nabla e_{\mu}\|_{0}$
---|---|---|---
$\left(\tau_{0},h_{0}\right)$ | 4.5526e-2(—) | 2.8424e-2(—) | 6.3754e-2(—)
$\left(\tau_{0},h_{0}\right)/2$ | 2.2281e-2(1.00) | 1.4898e-2(0.93) | 3.2766e-2(0.96)
$\left(\tau_{0},h_{0}\right)/4$ | 1.1450e-2(0.99) | 7.5740e-3(0.98) | 1.6567e-2(0.98)
$\left(\tau_{0},h_{0}\right)/8$ | 5.7815e-3(0.99) | 3.8063e-3(0.99) | 8.3157e-3(0.99)
$\left(\tau_{0},h_{0}\right)/16$ | 2.900e-3(1.00) | 1.9060e-3(1.00) | 4.1632e-3(1.00)
$\left(\tau_{0},h_{0}\right)/32$ | 1.4520e-3(1.00) | 9.5336e-4(1.00) | 2.0825e-3(1.00)

Finally, we verify the exactly divergence-free property of the discrete
current density. From the last column of Tables 3 and 5, the approximate
solutions yields $\left\|\nabla\cdot\boldsymbol{J}_{h}\right\|_{0}$ in the
order of $10^{-11}\sim 10^{-13}$, which is almost divergence-free. These tiny
errors mainly result from the numerical integral errors and rounding errors.

######  Example 4.12 (Shape relaxation).

In this example, we simulate the evolution of a square shaped bubble and two
kissing circular bubbles in the domain $\Omega=(0,1)^{2}$. We set the external
magnetic field as $\boldsymbol{B}=\left(0,0,1\right)^{{\rm T}}$ and the
physical parameters

$\eta=\sigma=1,\quad\epsilon=0.01,\quad\gamma=M=0.1.$

The initial velocity is taken as zero vector. For a square shaped bubble, the
initial profile of phase function $\varphi$ is chosen to be

$\varphi_{0}=\tanh\left(\frac{|x+y-1|+|x-y|-0.4}{\sqrt{2}\varepsilon}\right).$

For two kissing circular bubbles, $\varphi_{0}$ is taken as

$\varphi_{0}=1-\tanh\left(\frac{\left\|x-x_{o_{1}}\right\|-r_{1}}{\sqrt{2}\varepsilon}\right)-\tanh\left(\frac{\left\|x-x_{o_{2}}\right\|-r_{2}}{\sqrt{2}\varepsilon}\right)$

where$\left\|x-x_{o_{i}}\right\|$ is the Eulerian distance between the points
$x$ and $x_{o_{i}}$. The two points $x_{o_{1}}=\left(0.3,0.5\right)$ and
$x_{o_{2}}=\left(0.7,0.5\right)$ , are the center of two bubbles,
$r_{1}=r_{2}=0.2$ are their radius.

With the prescribed data, we conduct the numerical experiments with the mesh
size $h=1/64$ and time step $\tau=0.01$, and calculate the total energy ${\rm
E}$ and the mass $\int_{\Omega}\phi_{h}^{n}{\rm d}x$ at each time step.

For a square shaped bubble, we choose the terminal time $T=1$. Fig. 1 shows
the evolution curve of the energy ${\rm E}$ and the mass
$\int_{\Omega}\phi_{h}^{n}{\rm d}x$. We notice that energy decay monotonically
and the mass $\int_{\Omega}\phi_{h}^{n}{\rm d}x$ remains a constant. This
confirms that our decoupled scheme is energy stable and mass conservative.
Snapshots of the phase evolution at different time are presented in Fig. 2. We
observe that the isolated square shape relaxes to a circular shape, due to the
effect of surface tension.

(a) The energy

(b) The mass

Figure 1: Time evolution of the energy and mass for a square shaped bubble.

(a) $t=0$

(b) $t=0.01$

(c) $t=0.21$

(d) $t=0$.36

(e) $t=0.51$

(f) $t=1.00$

Figure 2: Snapshots of the relaxation of a square shape. The black circle
indicates the zero-level set of $\varphi_{h}^{n}$.

For two kissing circular bubbles, we choose the terminal time $T=15$. The time
evolution of the energy ${\rm E}$ and the mass $\int_{\Omega}\phi_{h}^{n}{\rm
d}x$ are displayed in Fig. 3. We still observe that energy curves decay
monotonically and the mass $\int_{\Omega}\phi_{h}^{n}{\rm d}x$ remains
constant in time. Again, this confirms that the decoupled scheme is energy
stable and mass-conservative. Fig. 4 displays some snapshots of the phase
evolution. From this figure, one can see that as time evolves, the two bubbles
quickly connect together and eventually coalesces into one big bubble under
the influence of surface tension.

(a) The energy

(b) The mass

Figure 3: Time evolution of the energy and mass for two kissing circular
bubbles.

(a) $t=0$

(b) $t=0.01$

(c) $t=0.51$

(d) $t=2.51$

(e) $t=5.01$

(f) $t=7.51$

(g) $t=10.01$

(h) 12.51

(i) $t=15.00$

Figure 4: Snapshots of the relaxation of two kissing circular bubbles. The
black circle indicates the zero-level set of $\varphi_{h}^{n}$.

######  Example 4.13 (Kelvin–Helmholtz instability).

In this example, we simulate the Kelvin–Helmholtz instability [46], which is
one of the most fundamental instabilities in incompressible fluids.
Considering the model problem on a rectangular domain
$\Omega=(0,0.5)\times(0,1)$, we choose
$\boldsymbol{B}=\left(0,0,1\right)^{{\rm T}}$,
$\eta=0.0002,\,\sigma=1,\,\epsilon=M=0.01,$ and $\gamma=0.001.$ The initial
conditions for the phase function $\varphi$ and velocity $\boldsymbol{u}$ are
given by

$\phi_{0}=\tanh\left(\frac{6\left(y-y_{c}\right)}{\sqrt{2}\epsilon}\right),\quad\boldsymbol{u}_{0}=\left(\tanh\left(50\left(y-y_{c}\right)\right),0\right)^{{}^{{\rm
T}}},$

where $y_{c}=0.5+0.005\sin\left(4\pi x\right)$. The boundary condition for
velocity is

$\boldsymbol{u}=\left(\pm 1,0\right)^{{\rm T}}\quad\text{ at }y=\pm
1,\quad\boldsymbol{u}=\left(0,0\right)^{{\rm T}}\quad\text{ at }x=\pm 1.$

A uniform mesh with the step size $h=1/128$ and a uniform time partition with
the time step size $\tau=0.01$ are used in this simulation. Fig. 5 shows
several snapshots of the phase field to illustrate the evolution of the
shearing interface. We observe that the flow sweeps the initial interfacial
vorticity into the center in early stage. As vorticity accumulates at the
center, the interface begins to steepen and the height of the instability gets
larger. At late stage, roll-up follows and the interface evolves into a
spiral. The dynamics of the interface is similar to those obtained in [47,
48].

(a) $t=0$

(b) $t=0.21$

(c) $t=0.41$

(d) $t=0.61$

(e) $t=0.81$

(f) $t=1.01$

(g) $t=1.21$

(h) $t=1.61$

(i) $t=2.00$

Figure 5: Snapshots of the relaxation of the shearing interface. The black
circle indicates the zero-level set of $\varphi_{h}^{n}$.

######  Example 4.14 (Gravity-driven flow).

This example is to study the effect of the gravity on the two-phase fluid.
Similar to [28], we supplement the gravity force as a force term
$\boldsymbol{f}_{g}$ on the right-hand side of the momentum equation (1),

$\boldsymbol{f}_{g}=\boldsymbol{g}\frac{H\left(\varphi\right)+1}{2},$

where
$H\left(\varphi\right)=\frac{1}{1+\exp\left(-\frac{\varphi}{\epsilon}\right)}$
is a regularized approximation of the Heaviside step function, and
$\boldsymbol{g}=(0,10)$ is the gravity. The physical parameters are given as
$\eta=\epsilon=0.01,\sigma=100,\,M=0.001,\ T=2.5.$ The initial data of
$\boldsymbol{u}$ is $\boldsymbol{0}$ and of $\varphi$ is set as

$\varphi_{0}=-\tanh\left(\frac{\left\|x-x_{o}\right\|-r}{\sqrt{2}\varepsilon}\right),$

where $x_{o}=\left(0.5,0.8\right)$, and $r=0.1$.

We simulate the CHIMHD system using the proposed scheme with $h=0.01$ and
$\tau=0.005$. In Fig. 6, we plot some snapshots of the phase-field profile
$\varphi$ with $\gamma=0.01,0.005,0.001$. We observe that the effect of the
gravity the circle bubble changes shape as it falls, and it transforms to be
flat when it approaches to the boundary. The smaller $\gamma$ is, the greater
deformation. Since the Ginzburg–Landau energy models adhesion forces, it can
be expected that a reduction of the parameter $\gamma$ reduces adhesion forces
and leads to instabilities. Thus, the results verify this point.

|  |  |
---|---|---|---
|  |  |
|  |  |
$t=1.01$ | $t=1.51$ | $t=2.01$ | $t=2.5$
Figure 6: Profile of the phase fields for $\gamma=0.01,0.005,0.001$.(From top
to bottom.)

##  5 Summary

In this paper, we propose a linear and decoupled finite element method for the
CHIMHD model. This full discrete scheme is based on first order Euler semi-
implicit scheme with some first order stabilization terms and implicit-
explicit treatments for time discretization, and stable mixed finite element
approximation for space discretization. In particular, we solve the current
density and electric potential simultaneously by using stable face-volume
mixed finite element pairs to ensure the discrete current density are
divergence-free exactly. The scheme is proved to be mass-conservative, charge-
conservative and unconditionally energy stable. We performed some numerical
tests to verify the features, accuracy and efficiency of the proposed scheme.
In the further work, we will study the highly efficient scheme for the two-
phase IMHD flows with large density ratio.

## References

## References

  * [1] D. M. Anderson, G. B. McFadden, A. A. Wheeler, Diffuse-interface methods in fluid mechanics, in: Annual review of fluid mechanics, Vol. 30, Vol. 30 of Annu. Rev. Fluid Mech., Annual Reviews, Palo Alto, CA, 1998, pp. 139–165.
  * [2] H. Ding, P. D. Spelt, C. Shu, Diffuse interface model for incompressible two-phase flows with large density ratios, Journal of Computational Physics 226 (2) (2007) 2078 – 2095.
  * [3] C. Liu, J. Shen, A phase field model for the mixture of two incompressible fluids and its approximation by a Fourier-spectral method, Phys. D 179 (3-4) (2003) 211–228.
  * [4] L. Rayleigh, On the theory of surface forces. ii. compressible fluids, Philosophical Magazine Series 5 33 (1892) 209–220.
  * [5] J. D. van der Waals, The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density, Journal of Statistical Physics 20 (1893) 197–200.
  * [6] S. M. A. J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metallurgica 27 (1979) 1085–1095.
  * [7] J. W. Cahn, J. E. Hilliard, Free energy of a nonuniform system. I. Interfacial Free Energy, The Journal of Chemical Physics 28 (2) (1958) 258.
  * [8] J. Kim, Phase-field models for multi-component fluid flows, Commun. Comput. Phys. 12 (3) (2012) 613–661.
  * [9] J. Shen, Modeling and numerical approximation of two-phase incompressible flows by a phase-field approach, in: Multiscale modeling and analysis for materials simulation, Vol. 22 of Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap., World Sci. Publ., Hackensack, NJ, 2012, pp. 147–195.
  * [10] J.-F. Gerbeau, A stabilized finite element method for the incompressible magnetohydrodynamic equations, Numer. Math. 87 (1) (2000) 83–111.
  * [11] J.-F. Gerbeau, C. Le Bris, T. Lelièvre, Mathematical methods for the magnetohydrodynamics of liquid metals, Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2006.
  * [12] R. Moreau, Magnetohydrodynamics, Vol. 3 of Fluid Mechanics and its Applications, Kluwer Academic Publishers Group, Dordrecht, 1990.
  * [13] S. Badia, A. F. Martín, R. Planas, Block recursive LU preconditioners for the thermally coupled incompressible inductionless MHD problem, J. Comput. Phys. 274 (2014) 562–591.
  * [14] L. Li, M. Ni, W. Zheng, A charge-conservative finite element method for inductionless MHD equations. Part I: Convergence, SIAM J. Sci. Comput. 41 (4) (2019) B796–B815.
  * [15] M.-J. Ni, R. Munipalli, N. B. Morley, P. Huang, M. A. Abdou, A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. I. On a rectangular collocated grid system, J. Comput. Phys. 227 (1) (2007) 174–204.
  * [16] M. Abdou, T. A. TEAM, A. Ying, N. Morley, K. Gulec, S. Smolentsev, M. Kotschenreuther, S. Malang, S. Zinkle, T. Rognlien, P. Fogarty, B. Nelson, R. Nygren, K. McCarthy, M. Youssef, N. Ghoniem, D. Sze, C. Wong, M. Sawan, H. Khater, R. Woolley, R. Mattas, R. Moir, S. Sharafat, J. Brooks, A. Hassanein, D. Petti, M. Tillack, M. Ulrickson, T. Uchimoto, On the exploration of innovative concepts for fusion chamber technology, Fusion Engineering and Design 54 (2) (2001) 181–247.
  * [17] P. A. Davidson, An introduction to magnetohydrodynamics, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2001.
  * [18] D. Munger, A. Vincent, A level set approach to simulate magnetohydrodynamic-instabilities in aluminum reduction cells, J. Comput. Phys. 217 (2) (2006) 295–311.
  * [19] M. Abdou, D. Sze, C. Wong, M. Sawan, A. Ying, N. B. Morley, S. Malang, U.S. plans and strategy for ITER blanket testing, Fusion Science and Technology 47 (3) (2005) 475–487.
  * [20] R. Samulyak, J. Du, J. Glimm, Z. Xu, A numerical algorithm for MHD of free surface flows at low magnetic Reynolds numbers, J. Comput. Phys. 226 (2) (2007) 1532–1549.
  * [21] J. Zhang, M. Ni, Direct numerical simulations of incompressible multiphase magnetohydrodynamics with phase change, J. Comput. Phys. 375 (2018) 717–746.
  * [22] H. Huang, A. Ying, M. Abdou, 3d mhd free surface fluid flow simulation based on magnetic-field induction equations, Fusion Engineering and Design 63-64 (2002) 361–368.
  * [23] Y. Pan, J. Zhang, Z.-H. Wang, M.-J. Ni, Development of mhd solver based on an adaptive mesh refinement technique, Fusion Engineering and Design 87 (5-6) (2012) 630–633.
  * [24] K. Takatani, Mathematical modeling of incompressible mhd flows with free surface, ISIJ Int. 47 (4) (2007) 545–551.
  * [25] H. Ki, Level set method for two-phase incompressible flows under magnetic fields, Comput. Phys. Comm. 181 (6) (2010) 999–1007.
  * [26] W.-X. Xie, L. Cai, J.-H. Feng, Tracking entropy wave in ideal mhd equations by weighted ghost fluid method, Applied Mathematical Modelling 31 (11) (2007) 2503–2514.
  * [27] H. L. Liwei Ding, Yafei Cao, Z. Liu, Mhd numerical simulation of aluminum electrolytic cell (in chinese), Metal Materials and Metallurgy Engineering 42 (4) (2014) 8–13.
  * [28] R. Chen, H. Zhang, Second-order energy stable schemes for the new model of the Cahn-Hilliard-MHD equations, Adv. Comput. Math. 46 (6) (2020) 79.
  * [29] X. Zhang, Sharp-interface limits of the diffuse interface model for two-phase inductionless magnetohydrodynamic fluids, arXiv preprint arXiv:2106.10433 (2021).
  * [30] S. Mao, X. Wang, Fully discrete finite element approximation of a three-dimensional diffuse interface model for two-phase incompressible inductionless magnetohydrodynamic fluids, Submitted (2021).
  * [31] M.-J. Ni, J.-F. Li, A consistent and conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part III: On a staggered mesh, J. Comput. Phys. 231 (2) (2012) 281–298.
  * [32] J. Hua, P. Lin, C. Liu, Q. Wang, Energy law preserving c0 finite element schemes for phase field models in two-phase flow computations, Journal of Computational Physics 230 (19) (2011) 7115–7131.
  * [33] H. Johnston, J.-G. Liu, Accurate, stable and efficient Navier-Stokes solvers based on explicit treatment of the pressure term, J. Comput. Phys. 199 (1) (2004) 221–259.
  * [34] J. Shen, X. Yang, Decoupled, energy stable schemes for phase-field models of two-phase incompressible flows, SIAM J. Numer. Anal. 53 (1) (2015) 279–296.
  * [35] J. Shen, X. Yang, Numerical approximations of Allen-Cahn and Cahn-Hilliard equations, Discrete Contin. Dyn. Syst. 28 (4) (2010) 1669–1691.
  * [36] L. A. Caffarelli, N. E. Muler, An $L^{\infty}$ bound for solutions of the Cahn-Hilliard equation, Arch. Rational Mech. Anal. 133 (2) (1995) 129–144.
  * [37] J. Shen, X. Yang, A phase-field model and its numerical approximation for two-phase incompressible flows with different densities and viscosities, SIAM J. Sci. Comput. 32 (3) (2010) 1159–1179.
  * [38] P. Monk, Finite element methods for Maxwell’s equations, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 2003\.
  * [39] Q. D. Xiaodi Zhang, A decoupled, unconditionally energy stable and charge-conservative finite element method for inductionless magnetohydrodynamic equations, Submitted (2020).
  * [40] D. J. Eyre, An unconditionally stable one-step scheme for gradient systems, unpublished (1997) 15.
  * [41] X. Yang, J. Zhao, Q. Wang, J. Shen, Numerical approximations for a three-component Cahn-Hilliard phase-field model based on the invariant energy quadratization method, Math. Models Methods Appl. Sci. 27 (11) (2017) 1993–2030.
  * [42] J. Shen, J. Xu, J. Yang, A new class of efficient and robust energy stable schemes for gradient flows, SIAM Rev. 61 (3) (2019) 474–506.
  * [43] F. Hecht, New development in freefem++, J. Numer. Math. 20 (3-4) (2012) 251–265.
  * [44] D. N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (4) (1984) 337–344 (1985).
  * [45] P. A. Raviart, J. M. Thomas, A mixed finite element method for 2-nd order elliptic problems, in: I. Galligani, E. Magenes (Eds.), Mathematical Aspects of Finite Element Methods, Springer Berlin Heidelberg, Berlin, Heidelberg, 1977, pp. 292–315.
  * [46] P. G. Drazin, W. H. Reid, Hydrodynamic stability, 2nd Edition, Cambridge Mathematical Library, Cambridge University Press, Cambridge, 2004.
  * [47] H. D. Ceniceros, R. L. Nós, A. M. Roma, Three-dimensional, fully adaptive simulations of phase-field fluid models, J. Comput. Phys. 229 (17) (2010) 6135–6155.
  * [48] H. G. Lee, J. Kim, Two-dimensional Kelvin-Helmholtz instabilities of multi-component fluids, Eur. J. Mech. B Fluids 49 (part A) (2015) 77–88.


# Controlling Quantum Chaos: Optimal Coherent Targeting

Steven Tomsovic<EMAIL_ADDRESS>Institut für Theoretische Physik,
Universität Regensburg, D-93040 Regensburg, Germany  Department of Physics and
Astronomy, Washington State University, Pullman, WA USA  Juan Diego Urbina
Institut für Theoretische Physik, Universität Regensburg, D-93040 Regensburg,
Germany  Klaus Richter  Institut für Theoretische Physik, Universität
Regensburg, D-93040 Regensburg, Germany

###### Abstract

One of the principal goals of controlling classical chaotic dynamical systems
is known as targeting, which is the very weakly perturbative process of using
the system’s extreme sensitivity to initial conditions in order to arrive at a
predetermined target state. It is shown that a generalization to chaotic
quantum systems is possible in the semiclassical regime, but requires tailored
perturbations whose effects must undo the dynamical spreading of the evolving
quantum state. The procedure described here is applied to initially minimum
uncertainty wave packets in the quantum kicked rotor, a preeminent quantum
chaotic paradigm, to illustrate the method, and investigate its accuracy. The
method’s error can be made to vanish as $\hbar\rightarrow 0$.

quantum chaos, control, semiclassical mechanics

One of the several meanings ascribed to the term controlling chaos [1, 2] is
the concept of targeting [3, 4, 5, 6] in which the exponential instability of
a chaotic dynamical system is taken advantage of in an optimal way. In Refs.
[3, 4] small perturbations of a system parameter leading to some predetermined
final state were discussed, but a very weak perturbation applied to the
system’s trajectory (initial conditions) can be applied as well. A well known
early example is the redirecting in late 1983 of the International Sun-Earth
Explorer-3 spacecraft towards a rendezvous with the Giacobini-Zimmer comet in
1985 [7, 2].

A natural question arises as to whether it is possible to bring this concept
from the classical to the quantum realm and to achieve some form of targeting
with highly excited or far-from-equilibrium quantum chaotic systems. In this
regard, three desiderata can serve as a blueprint. First, let us assume that
the quantum system is isolated and has a well defined classical analog.
Ideally, the initial quantum and target states of interest would be as
classical as possible given uncertainty relations. Second, an optimal quantum
transport path can be identified, conceivably guided by the classical analog’s
(chaotic) dynamics. Finally, it is preferable to avoid monitoring or
measurements, and to seek applying purely coherent unitary weak perturbations
to the quantum system’s dynamics in order to achieve the optimal targeting.
Clearly, such an optimal coherent chaotic quantum targeting procedure would
fall into the general classification of a quantum control problem.

The subject of quantum control theory has a long history; a survey of the
theory and applications can be found in [8] and a more recent pedagogical
overview of optimal control theory in [9]. Within this enormous body of work,
are some concepts a bit closer to the subject of this Letter. In [10],
Sugawara describes ‘wave packet shaping’, which is applied to the integrable
Morse oscillator potential. It gives a complicated time-dependent laser field
that shifts a wave packet from its ground state location to another desired
location. Another method makes use of phase space structural implications of
Kolmogorov-Arnol’d-Moser theory [11, 12, 13]; see Ref. [14], which describes
the creation of a non-dispersive electronic wave packet following a trajectory
about which the local phase space dynamics has an approximately harmonic
nature. However, the highly excited, strongly chaotic regime remains
challenging.

Figure 1: Schematic of optimal coherent chaotic quantum targeting. The
circular zone represents the Wigner transform density of a nearly minimal
uncertainty state, which can be slightly shifted (via a unitary operator,
$U_{s}$) to be centered on an optimal trajectory starting at $(q_{0},p_{0})$.
As it propagates the density follows an optimal chaotic trajectory, but is
locally spreading, which must be counteracted by contractions, $U_{M^{-1}}$.
At the end, it can be shifted from $(q_{\tau},p_{\tau})$ to the centroid of
the target state $(q_{\beta},p_{\beta})$.

Returning to the first desideratum, for a quantum chaotic system the initial
and target states ideally would be minimum uncertainty states with well
defined mean momenta and positions, $(q_{\alpha},p_{\alpha})$ and
$(q_{\beta},p_{\beta})$, respectively. The system would be in a semiclassical
regime, i.e. $\hbar$ is sufficiently small that an asymptotic theory, such as
advanced methods applicable to chaotic systems [15, 16] and based on time-
dependent Wentzel-Kramers-Brillouin (WKB) theory [17, 18, 19, 20, 21] describe
the propagation well. The Wigner transforms [22] of these states are analogous
to localized Liouvillian phase space densities with volumes determined by
$\hbar$ as opposed to single trajectories and this adds an inescapable new
ingredient, arguably the most important one.

As sketched in Fig. 1, these densities provide a geometric picture for the
task of identifying an optimal quantum transport path, the second desideratum.
For a fully chaotic Hamiltonian dynamical system, there exist directions of
maximal exponential deviation under perturbations as well as maximal
exponential compression [23]. Intuitively, separating along the direction of
maximal exponential rate from the centroid of the initial density while
simultaneously converging at the maximal exponential rate towards the centroid
of the target density, gives an optimal dynamical path from the initial to
target state. For the desired arrival time, $\tau$, there exists such a
heteroclinic trajectory segment that maximizes it’s initial conditions’,
$(q_{0},p_{0})$, weight in the Liouvillian density of the initial state, and
similarly for it’s endpoint, $(q_{\tau},p_{\tau})$, within the density for the
target state. Crudely speaking, $(q_{0},p_{0})$ and $(q_{\tau},p_{\tau})$ are
closest to $(q_{\alpha},p_{\alpha})$ and $(q_{\beta},p_{\beta})$,
respectively.

The methods for finding heteroclinic trajectory segments developed in [24, 25,
26] are sufficient for the purposes here, and it has previously been shown
that for chaotic systems the full set of heteroclinic trajectory segments can
be used to construct the full time evolution of transport coefficients
involving wave packets [27, 28, 25] and coherent states in chaotic many-body
bosonic systems [29, 26]. However, to provide an optimal transport pathway for
the quantum dynamics, here the idea is to follow just one segment, thus
evoking the third desideratum. Centering the density on the trajectory segment
and maintaining its local nature is sufficient to accomplish this. Both
operations, translations and dilations, fall into the class of linear
canonical transformations, which have precise unitary transformation
counterparts as was known very early on to Dirac [30], and was developed fully
later [31, 32]. There it was shown that these unitary transformations have at
most quadratic generators.

A unitary shift transformation, say $U_{s}$, slightly shifting the initial and
final state’s centroids towards the initial and final conditions of the
heteroclinic trajectory segment, is given by $U_{s}(\Delta p,\Delta
q)=\exp\left[i\left(\Delta p\hat{q}-\hat{p}\Delta q\right)/\hbar\right]$. In
configuration space, it is expressed as

$\braket{q}{U_{s}(\Delta p,\Delta
q)}{q^{\prime}}=\\!\exp\\!\left[\frac{i}{\hbar}\Delta
p\\!\left(q\\!-\\!\frac{\Delta
q}{2}\right)\right]\\!\delta\\!\left(q\\!-\\!q^{\prime}\\!-\\!\Delta q\right)$
(1)

where the global phase part, $-\Delta p\Delta q/(2\hbar)$, is most likely of
no concern and can be dropped. Curiously, small perturbations to guide the
system toward following a heteroclinic trajectory segment is a critical part
of classical targeting, but is not nearly as essential in the quantum context,
especially in the very strongly chaotic limit where no transport barriers are
present.

Countering the spreading of the quantum state under evolution is the most
critical operation and must be done periodically before the state is spread to
such an extent that the nonlinear dynamics becomes manifest in the shape of
the spreading Liouvillian density. The heteroclinic trajectory segment’s
stability matrix for each such time interval describes the successive wave
packet spreadings, and is describing a local linear canonical transformation.
The local unitary transformation associated with the inverse stability matrix
for that portion of the trajectory, say $U_{M^{-1}}$, unwinds the spreading
and returns the wave packet to its minimum uncertainty form.

Within the semiclassical propagation method called linearized wave packet
dynamics [33, 34] are an accounting for this spreading, and the method
effectively constructs the corresponding unitary transformation in a
configuration representation. The multidimensional expressions can be found in
[16], and for an alternative approach for suppression of spreading through
periodic nonlinear kicking see [35]. The main ingredients are the second
derivatives of the generating function, $S$, which can be expressed in terms
of the stability matrix elements (block ordering is consistent with the kicked
rotor equations ahead) [34]. To unwind this spreading using the stability
matrix follows using block $2$ x $2$ matrix inversion formulas [36] and some
algebra. One finds

$\displaystyle\left(\\!\begin{array}[]{cc}\frac{\partial^{2}S}{\partial\vec{q}\partial\vec{q}}&\frac{\partial^{2}S}{\partial\vec{q}^{\prime}\partial\vec{q}}\\\
\frac{\partial^{2}S}{\partial\vec{q}\partial\vec{q}^{\prime}}&\frac{\partial^{2}S}{\partial\vec{q}^{\prime}\partial\vec{q}^{\prime}}\end{array}\\!\right)=$
(4)
$\displaystyle\left(\\!\begin{array}[]{cc}-\bf{M^{-1}_{21}}\cdot\bf{M_{22}}&\bf{M^{-1}_{21}}\\\
\bf{M_{11}}\cdot\bf{M^{-1}_{21}}\cdot\bf{M_{22}}\\!-\\!\bf{M_{12}}&\\!\\!-\bf{M_{11}}\cdot\bf{M^{-1}_{21}}\end{array}\\!\right)\,.$
(7)

At the moment in time, $t$, that the spreading is countered, the
transformation must be centered locally about the heteroclinic trajectory’s
location $(q_{t},p_{t})$ using the coordinates $q-q_{t}$ and
$q^{\prime}-q_{t}$. In addition, centering the momentum coordinate in
configuration space is accomplished by multiplying the transformation by
$\exp\left[ip_{t}\left(q-q^{\prime}\right)/\hbar\right]$.

If the unitary dynamics of the uncontrolled quantum chaotic system of interest
is denoted by $U(t)$, then by creating a controlled quantum dynamics given by

$U_{\rm cqd}(\tau)=U_{s}\prod_{j=1}^{n}\\{U_{M^{-1}}U(t)\\}_{j}U_{s}\,,$ (8)

with $U_{M^{-1}}$ unwinding spreading, $\tau=nt$, and the shift operators
$U_{s}$, Eq. (1), the initial quantum state can be made to follow any
trajectory segment that exists in the classical analog dynamical system. This
technique of optimal coherent chaotic quantum targeting is illustrated for the
quantum kicked rotor next.

The kicked rotor has a long history of providing a simple, powerful paradigm
for both classical and quantum dynamical systems [37, 38, 39, 40]. Its quantum
version has been experimentally realized with cold atoms in a kicked optical
lattice for a variety of purposes [41, 42, 43, 44]. The classical Hamiltonian
is given by

$H(q,p)=\frac{p^{2}}{2}-\frac{K}{4\pi^{2}}\cos(2\pi
q)\sum_{n=-\infty}^{\infty}\delta(t-n)\,.$ (9)

By choosing the kicking strength greater than roughly $2\pi$, it generates a
strongly chaotic classical dynamics. The resulting mapping equations for the
version on the unit-periodicity phase space torus are

$\begin{split}&p_{n+1}=p_{n}-\frac{K}{2\pi}\sin(2\pi q_{n})\ \pmod{1}\,,\\\
&q_{n+1}=q_{n}+p_{n+1}\qquad\qquad\pmod{1}\,.\end{split}$ (10)

The single step stability matrix,

$\left(\begin{array}[]{c}\delta p_{n+1}\\\ \delta
q_{n+1}\end{array}\right)={\bf M}_{n}\left(\begin{array}[]{c}\delta p_{n}\\\
\delta q_{n}\end{array}\right)\;,$ (11)

is critical for the targeting $U_{M^{-1}}$ and given by

${\bf M}_{n}=\left(\begin{array}[]{cc}m_{11}&m_{21}\\\
m_{12}&m_{22}\end{array}\right)=\left(\begin{array}[]{cc}1&-K\cos\left(2\pi
q_{n}\right)\\\ 1&1-K\cos\left(2\pi q_{n}\right)\end{array}\right)$ (12)

and with inverse given by

${\bf M}_{n}^{-1}=\left(\begin{array}[]{cc}m_{22}&-m_{21}\\\
-m_{12}&m_{11}\end{array}\right)\ .$ (13)

The quantum dynamics are generated by a unitary or Floquet operator,

$\widehat{U}=\exp\left(\displaystyle\frac{-i\widehat{p}^{2}}{2\hbar}\right)\;\exp\left[\displaystyle\frac{iK}{4\pi\hbar^{2}}\cos
2\pi\widehat{q}\right]\;.$ (14)

In a configuration representation, $U_{jk}=\langle q_{j}\mid\widehat{U}\mid
q_{k}\rangle$, with null Bloch phases it becomes

U_jk = 1iN exp[ iπ(j-k)2N ] exp[ iNK2π cos2πkN ] ,

(15)

where $N$ is the Hilbert space dimension, $j,k=1,...N$, and Planck’s constant
is $2\pi\hbar=1/N$. The semiclassical limit of $\hbar\rightarrow 0$ is
equivalent to $N\\!\rightarrow\\!\infty$. The kicking strength
$K\\!=\\!4\pi^{2}/5$ and Hilbert space dimensionality $N=500$ are chosen for
Fig. 2-4.

The initial state $\ket{\alpha}$ reads

$\braket{q_{j}}{\alpha}=A(\hbar)\exp\\!\left[\\!-\frac{\left(q_{j}-q_{\alpha}\right)^{2}}{2\hbar}\\!+\\!\frac{i}{\hbar}p_{\alpha}\left(q_{j}\\!-\\!q_{\alpha}\right)\right]$
(16)

where $A(\hbar)$ is a normalization constant, and likewise for the target
state $\ket{\beta}$. This form shares momentum and position minimum
uncertainty equally on the unit torus. Thus, the wave packet’s phase space
analogy has a circular symmetry for all $\hbar$. The area inside the wave
packet’s two standard deviation contour is equal to $h\ (=1/N)$. The integer
value $m$ of the position $q_{j}=m+j/N$ is chosen such that
$\left|q_{j}-q_{\alpha}\right|\leq 1/2$ for all $j$. In this way,
$p_{\alpha}=\braket{\alpha}{\widehat{p}}{\alpha}$ and likewise,
$q_{\alpha}=\braket{\alpha}{\widehat{q}}{\alpha}$. The selected wave packet
centroids are $(q_{\alpha},p_{\alpha})=(0.5,0)$ and
$(q_{\beta},p_{\beta})=(0,0)$, respectively.

Figure 2: Quantum chaotic spreading of the kicked rotor initial state $|\alpha\rangle$. The real (purple solid) and absolute (dashed green) values of the one and three time propagations of the initial state, $\braket{q}{\alpha(1)}$ (top panel) and $\braket{q}{\alpha(3)}$ (bottom panel). If uncontrolled the state spreads and rapidly acquires an ergodic looking appearance. Table 1: Targeting heteroclinic trajectory. Values $n=1,...4$ correspond to position and momentum centroids of the propagated states in the four panels of Fig. 3. n | $q_{n}$ | $p_{n}$
---|---|---
0 |  0.50063152866 |  0.00056701866
1 | 0.50618488472 | 0.00555335605
2 | 0.56055984292 | 0.05437495820
3 | 1.08164081735 | 0.52108097443
4 | 0.98601176250 | -0.09562905485

In such a strongly chaotic system the initial state spreads rapidly on a
logarithmically short time scale, the Ehrenfest time,
$\tau_{E}\\!=\\!\ln(S/\hbar)/\mu$; see Fig. 2. The classical action $S$ can be
taken as unity for the unit torus kicked rotor, and the Lyapunov exponent,
$\mu\\!=\\!\ln(K/2)$. Here, $\tau_{E}\approx 4$, which implies ergodic
statistical behavior by this time. Further propagation appears random.

Figure 3: Quantum guiding along heteroclinic classical paths. The real (solid
purple) and absolute (dashed green) values of $\braket{q}{\alpha(t)}$ are
obtained using Eq. (8). At each step the propagated state has a near perfect
overlap (not discernible directly from the figure) with a wave packet of the
form of Eq. (16) with the time appropriate point of the heteroclinic
trajectory in Table 1.

To illustrate quantum targeting for the kicked rotor, a heteroclinic
trajectory segment of four iterations is chosen and given in Table 1. The
initial point is quite close to $(q_{\alpha},p_{\alpha})$ and the final point
not quite as close to $(q_{\beta},p_{\beta})$. Extending this trajectory
segment to six iterations, one prior and one latter, would make using the
before and after shifts, $U_{s}$, largely unnecessary in Eq. (8); they are
kept in this example. With the wave packet spreading countered at each time
step, as in Eq. (8), it can be seen that the initial state follows the chosen
heteroclinic trajectory exceedingly well; see Fig. 3. In fact, the absolute
squared overlap of the evolving state following Eq. (8), with Gaussian wave
packets created using Eq. (16) and the successive heteroclinic phase points,
is equivalent to the accuracy shown in Fig. 5 ahead along the entire
trajectory segment.

The important role of $U_{M^{-1}}$ is illustrated for the second iteration of
the Floquet operator in Fig. 4. With $\ket{\alpha^{\prime}(2)}$ the propagated
ket before applying $U_{M^{-1}}$ and $\ket{\alpha(2)}$ afterwards, the effect
of the inverse stability matrix of the single step of the heteroclinic
trajectory is markedly relocalizing, in fact bringing the quantum state back
to minimum uncertainty. This technique is quite flexible in the sense that one
could undo the stretching continuously in time or even double the time shown
here by preemptively undoing the stretching of the coming next iteration or
other choices.

Figure 4: Countering wave packet spreading. The real (solid purple) and
absolute (dashed green) values of $\braket{q}{\alpha^{\prime}(2)}$ are shown.
Applying the unitary transformation associated with the local linear canonical
transformation of the inverse one-step stability matrix undoes wave packet
spreading without shifting the position or momentum centroids and generates
the second panel from the top of Fig. 3.

As a function of a shrinking $2\pi\hbar=1/N$, Fig. 5 demonstrates the
exponentially rapid decrease in the difference from unity of the overlap
between the target state and the state propagated via quantum control, $U_{\rm
cqd}(\tau)$, for the quantum kicked rotor. This serves as a proof of principle
and demonstrates the efficiency of our method.

Figure 5: Fidelity of quantum targeting for the kicked rotor. The squared
overlap of the target state with the controlled propagated state, i.e.
$|\braket{\beta}{U_{cqd}(\tau=4)}{\alpha(0)}|^{2}$, see Eq. (8), is plotted
versus Planck’s constant. As $\hbar\rightarrow 0$, the error tends to vanish
exponentially. Dimensionalities ranging from $N=100$ to $1450$ are included.

There are two main sources of inaccuracies arising in this illustration, one
due to the nonlinearities in the local dynamics and the other due to working
with the periodic nature of the torus in this specific example. With regards
to the nonlinearities, it gives rise to curvature of the stable and unstable
manifolds. With regards to the torus, the Gaussian construction above, Eq.
(16), is not explicitly periodic. Both sources lead to vanishing errors as
$\hbar\rightarrow 0$.

In conclusion, a generalization of targeting control in classically chaotic
systems is constructed for quantum systems that is fully coherent (unitary).
Using the system’s extreme sensitivity to initial conditions, counteracting
the quantum spreading (scrambling) is demonstrated, and by means of a weak
perturbation, the quantum system is put on an optimum track in order to arrive
in an exponentially fast way and with high precision at a predetermined target
state and given time. Hence, in a sense there are no free optimization
parameters as everything is determined by the heteroclinic trajectory
segment(s) chosen and its short time stability matrices. The method is
flexible in the sense that any trajectory segment of the unperturbed chaotic
classical system can be selected and followed, plus the unwinding of the wave
packet uncertainty spreading can be accomplished in multiple ways. In contrast
to classical targeting control, the most critical element in the quantum case
is the local unitary transformation $U_{M^{-1}}$.

The underlying protocol, comprising time sequences of tailored $U_{M^{-1}}$,
can be regarded as a sophisticated Floquet engineering method. Furthermore, if
the initial state is specifically chosen as the target state, stabilized
periodic quantum dynamics (with any desired period) can be achieved.

Our proposed method is, by construction, rather general and may serve as a
conceptual platform for a broader range of applications in various branches of
quantum control. In particular, the ideas presented here naturally extend to
coherent states in bosonic many-body systems as linear canonical
transformations are related to squeezing [45, 46], and that is left for future
work.

###### Acknowledgements.

We thank Mathias Steinhuber for scientific discussions and help with devising
Fig. 1. We acknowledge support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation), project Ri681/15-1, within the Reinhart-Koselleck
Programme and Vielberth Foundation for financial support.

## References

  * Ott  _et al._ [1990] E. Ott, C. Grebogi, and J. A. Yorke, Controlling chaos, Phys. Rev. Lett. 64, 1196 (1990).
  * Ott [2006] E. Ott, Controlling chaos, Scholarpedia 1, 1699 (2006), revision #91167.
  * Shinbrot  _et al._ [1990] T. Shinbrot, E. Ott, C. Grebogi, and J. A. Yorke, Using chaos to direct trajectories to targets, Phys. Rev. Lett. 65, 3215 (1990).
  * Kostelich  _et al._ [1993] E. J. Kostelich, C. Grebogi, E. Ott, and J. A. Yorke, Higher-dimensional targeting, Phys. Rev. E 47, 305 (1993).
  * Bollt [1995] E. M. Bollt,  _Controlling chaos, targeting, and transport_ , Ph.D. thesis, University of Colorado (1995).
  * Schroer and Ott [1997] C. G. Schroer and E. Ott, Targeting in Hamiltonian systems that have mixed regular/chaotic phase spaces, Chaos 7, 512 (1997).
  * Farquhar  _et al._ [1985] R. Farquhar, D. Muhonen, and S. A. Davis, Trajectories and orbit maneuvers for the ISEE-3/ICE comet mission (1985), J. Astronaut. Sci. 33, 235 (1985).
  * Dong and Petersen [2010] D. Dong and I. R. Petersen, Quantum control theory and applications: A survey, IET Control Theory Appl. 4, 2651 (2010).
  * James [2021] M. R. James, Optimal quantum control theory, Annu. rev. control robot. auton. syst. 4, 343 (2021).
  * Sugawara [2003] M. Sugawara, General formulation of locally designed coherent control theory for quantum system, J. Chem. Phys. 118, 6784 (2003).
  * Kolmogorov [1954] A. N. Kolmogorov, On the conservation of conditionally periodic motions under small perturbation of the hamiltonian, Dokl. Akad. Nauk SSSR 98, 527 (1954).
  * Arnol’d [1963] V. I. Arnol’d, Proof of a theorem of a. n. kolmogorov on the preservation of conditionally periodic motions under a small perturbation of the hamiltonian, Russ. Math. Surv. 18, 9 (1963).
  * Moser [1962] J. Moser, On invariant curves of area-preserving mappings of an annulus, Nachr. Akad. Wiss. Göttingen, II, 673 (1962).
  * Madronero  _et al._ [2006] J. Madronero, A. Ponomarev, A. R. R. Carvalho, S. Wimberger, C. Viviescas, A. Kolovsky, K. Hornberger, P. Schlagheck, A. Krug, and A. Buchleitner, Quantum chaos, transport, and control - in quantum optics, Adv. in Atom. Mol. Opt. Phys. 53, 33 (2006).
  * Huber  _et al._ [1988] D. Huber, E. J. Heller, and R. G. Littlejohn, Generalized gaussian wave packet dynamics, schrödinger equation, and stationary phase approximation, J. Chem. Phys. 89, 2003 (1988).
  * Pal  _et al._ [2016] H. Pal, M. Vyas, and S. Tomsovic, Generalized gaussian wave packet dynamics: Integrable and chaotic systems, Phys. Rev. E 93, 012213 (2016), arXiv:1510.08051 [quant-ph].
  * Maslov and Fedoriuk [1981] V. P. Maslov and M. V. Fedoriuk,  _Semiclassical approximation in quantum mechanics_ (Reidel Publishing Company, Dordrecht, 1981).
  * Keller [1958] J. B. Keller, Corrected Bohr-Sommerfeld quantum conditions for nonseparable systems, Ann. Phys. (N.Y.) 4, 180 (1958).
  * Wentzel [1926] G. Wentzel, Eine verallgemeinerung der quantenbedingungen für die zwecke der wellenmechanik, Zeit. der Phys. 38, 518 (1926).
  * Kramers [1926] H. A. Kramers, Wellenmechanik und halbzhlige quantisierung, Zeit. der Phys. 39, 828 (1926).
  * Brillouin [1926] L. Brillouin, La mécanique ondulatoire de schrödinger: une méthode générale de resolution par approximations successives, Comptes Rendus de l’Acad. des Sci. 183, 24 (1926).
  * Wigner [1932] E. P. Wigner, On the quantum correction for thermodynamic equilibrium, Phys. Rev. 40, 749 (1932).
  * Poincaré [1899] H. Poincaré,  _Les méthodes nouvelles de la mécanique céleste_ , Vol. 3 (Gauthier-Villars et fils, Paris, 1899).
  * Li and Tomsovic [2017] J. Li and S. Tomsovic, Accurate determination of heteroclinic orbits in chaotic dynamical systems, J. Phys. A: Math. Theor. 50, 135101 (2017), arXiv:1507.06455 [nlin.CD].
  * Tomsovic and Heller [1993] S. Tomsovic and E. J. Heller, The long-time semiclassical dynamics of chaos: the stadium billiard, Phys. Rev. E 47, 282 (1993).
  * Tomsovic [2018] S. Tomsovic, Complex saddle trajectories for multidimensional quantum wave packet/coherent state propagation: application to a many-body system, Phys. Rev. E 98, 023301 (2018), arXiv:1804.10511 [cond-mat.stat-mech].
  * Tomsovic and Heller [1991] S. Tomsovic and E. J. Heller, Semiclassical dynamics of chaotic motion: Unexpected long time accuracy, Phys. Rev. Lett. 67, 664 (1991).
  * O’Connor  _et al._ [1992] P. W. O’Connor, S. Tomsovic, and E. J. Heller, Semiclassical dynamics in the strongly chaotic regime - breaking the log-time barrier, Physica D 55, 340 (1992).
  * Tomsovic  _et al._ [2018] S. Tomsovic, P. Schlagheck, D. Ullmo, J.-D. Urbina, and K. Richter, Post-Ehrenfest many-body quantum interferences in ultracold atoms far-out-of-equilibrium, Phys. Rev. A 97, 061606(R) (2018), arXiv:1711.04693v2 [quant-ph].
  * Dirac [1930] P. A. M. Dirac,  _The Principles of Quantum Mechanics_ (Oxford University Press, Oxford, 1930).
  * Collins Jr. [1970] S. A. Collins Jr., Lens-system diffraction integral written in terms of matrix optics, J. Opt. Soc. Am. 60, 1168 (1970).
  * Moshinsky and Quesne [1971] M. Moshinsky and C. Quesne, Linear canonical transformations and their unitary representations, J. Math. Phys. 12, 1772 (1971).
  * Heller [1975] E. J. Heller, Time-dependent approach to semiclassical dynamics, J. Chem. Phys. 62, 1544 (1975).
  * Heller [1991] E. J. Heller, Wavepacket dynamics and quantum chaology, in  _Chaos and Quantum Physics_ , edited by M. J. Giannoni, A. Voros, and J. Zinn-Justin (North-Holland, Amsterdam, 1991) pp. 547–663.
  * Goussev  _et al._ [2018] A. Goussev, P. Reck, F. Moser, A. Moro, C. Gorini, and K. Richter, Overcoming dispersive spreading of quantum wave packets via periodic nonlinear kicking, Phys. Rev. A 98, 013620 (2018).
  * Lu and Shiou [2002] T.-T. Lu and S.-H. Shiou, Inverses of $2\times 2$ block matrices, Comput. Math. Appl. 43, 119 (2002).
  * Chirikov [1979] B. V. Chirikov, A universal instability of many-dimensional oscillator systems, Phys. Rep. 52, 263 (1979).
  * Izrailev [1990] F. M. Izrailev, Simple models of quantum chaos: Spectrum and eigenfunctions, Phys. Rep. 196-392, 299 (1990).
  * Lakshminarayan [1997] A. Lakshminarayan, Accuracy of trace formulas, Pramana 48, 517 (1997).
  * Tomsovic and Lakshminarayan [2007] S. Tomsovic and A. Lakshminarayan, Fluctuations of finite-time stability exponents in the standard map and the detection of small islands, Phys. Rev. E 76, 036207 (2007), arXiv:0706.1494 [nlin.CD].
  * Moore  _et al._ [1995] F. L. Moore, J. C. Robinson, C. F. Bharucha, B. Sundaram, and M. G. Raizen, Atom optics realization of the quantum $\delta$-kicked rotor, Phys. Rev. Lett. 75, 4598 (1995).
  * Ammann  _et al._ [1998] H. Ammann, R. Gray, I. Shvarchuck, and N. Christensen, Quantum delta-kicked rotor: Experimental observation of decoherence, Phys. Rev. Lett. 80, 4111 (1998).
  * Ryu  _et al._ [2006] C. Ryu, M. F. Anderen, A. Vaziri, M. B. d’Arcy, J. M. Grossman, K. Helmerson, and W. D. Phillips, High-order quantum resonances observed in a periodically kicked bose-einstein condensate, Phys. Rev. Lett. 96, 160403 (2006).
  * Sajjad  _et al._ [2022] R. Sajjad, J. L. Tanlimco, H. Mas, A. Cao, E. Nolasco-Martinez, E. Q. S. ans Flávio L. N. Santos ans Patrizia Vignolo, T. Macrì, and D. M. Weld, Observation of the quantum boomerang effect, Phys. Rev. X 12, 011035 (2022).
  * Han  _et al._ [1988] D. Han, Y. S. Kim, and M. E. Noz, Linear canonical transformations of coherent and squeezed states in the wigner phase space, Phys. Rev. A 37, 807 (1988).
  * Walls [1983] D. F. Walls, Squeezed states of light, Nature 306, 141 (1983).