11institutetext: Université de Reims Champagne Ardenne, CReSTIC, Reims, France
22institutetext: Univ Lyon, INSA‐Lyon, Université Claude Bernard Lyon 1, UJM-
Saint Etienne, CNRS, Inserm, CREATIS UMR 5220, U1294, F‐69100, Lyon, France
22email<EMAIL_ADDRESS>
# Cascaded multitask U-Net using topological loss for vessel segmentation and
centerline extraction††thanks: This work was supported by the French _Agence
Nationale de la Recherche (Grant ANR-20-CE45-0011)_ as part of the PreSPIN
project. This work was granted access to the HPC resources of IDRIS under the
allocation 2022-AD011013610 made by GENCI

Pierre Rougé 1122 Nicolas Passat 11 Odyssée Merveille 22

###### Abstract

Vessel segmentation and centerline extraction are two crucial preliminary
tasks for many computer-aided diagnosis tools dealing with vascular diseases.
Recently, deep-learning based methods have been widely applied to these tasks.
However, classic deep-learning approaches struggle to capture the complex
geometry and specific topology of vascular networks, which is of the utmost
importance in many applications. To overcome these limitations, the clDice
loss, a topological loss that focuses on the vessel centerlines, has been
recently proposed. This loss requires computing the skeletons of both the
ground truth and the predicted segmentation in a differentiable manner.
Currently, differentiable skeletonization algorithms are either inaccurate or
computationally demanding. In this paper, we propose to use a U-Net which
computes the vascular skeleton directly from the segmentation and the MRA
image. This method is naturally differentiable and provides a good trade-off
between accuracy and computation time. The resulting cascaded multitask U-Net
is trained with the clDice loss to embed topological constraints during the
segmentation. In addition to this topological guidance, this cascaded U-Net
also benefits from the inductive bias generated by the skeletonization during
the multitask training. This model is able to predict both the vessel
segmentation and centerlines with a more accurate topology and a low
computation time. Code is available at :
https://github.com/PierreRouge/Cascaded-U-Net-for-vessel-segmentation

###### Keywords:

vessel segmentation, U-Net, topological loss, deep-learning.

##  1 Introduction

Vascular diseases include all the alterations that occur in vessels, such as
stenosis, aneurysm, thrombosis or embolism. Consequences of these diseases,
such as stroke, are an important cause of death and disability worldwide. The
accurate vessel segmentation from angiographic images, actively investigated
over the last 30 years [5, 9], is an important step for the diagnosis and
treatment of vascular diseases.

In the last decade, deep-learning has allowed significant progress in medical
imaging [16] and especially in segmentation. In particular, it is increasingly
applied to the segmentation of vascular structures [10]. In this context, one
of the first approaches dealing with 3D angiographic images was proposed by
Kitrungrotsakul et al. [4] who used a multi pathway convolutional network
processing each of the three orthogonal planes to segment the hepatic vessels
in 3D. The same year, Livne et al. [6] and Sanches et al. [15] used for the
first time architectures building upon 2D and 3D U-Net, respectively, to
segment the brain vessels. From this time on, efforts were geared towards
designing networks dedicated to the segmentation of curvilinear structures.
Notably, Mou et al. [11] proposed a convolutional network with two attention
modules to encode both spatial and channel relationship. This method was
tested on six different imaging modalities and nine datasets (2D and 3D). In
the meantime, Tetteh et al. [18] designed an architecture aiming to perform
simultaneously vessel segmentation, centerline prediction and bifurcation
detection in angiographic images. Recently, Dang et al. [2] tackled the
problem of data annotation and proposed a weakly-supervised deep-learning
framework. The annotated patches were obtained using a classifier
discriminating between vessel and non-vessel patches and K-means algorithm.

Despite these efforts, automatic segmentation of vascular networks remains a
challenging topic, especially due to the complex topological and geometrical
properties of vessels, and their sparseness in the images. By contrast to many
anatomical structures, vessels do not constitute a compact volume at a
specific position and scale. They are organized as a multiscale network (from
large vessels to thin ones, close to the resolution of the image) in the whole
image. This represents a challenge for deep-learning methods, especially when
a topologically correct result is required for subsequent tasks such as blood
flow modeling [8].

To overcome these difficulties, Shit et al. [17] recently proposed a novel
metric specifically designed to evaluate the quality of tubular structure
segmentation. This metric, named clDice (for “centerline Dice”), mainly relies
on the skeleton of the tubular structures instead of their whole volume,
therefore focusing on topological information. To use this new metric as a
loss function, it is necessary to compute the skeleton of the predicted
segmentation in a differentiable manner. Therefore, the authors proposed a
differentiable soft-skeleton algorithm. However, the resulting skeletons do
not preserve the topology of the structures of interest. In a following work,
Menten et al. [7] proposed two new differentiable skeletonization algorithms
to overcome this limitation. They showed that using the clDice loss with this
soft-skeleton algorithm provides better and more connected segmentation of 2D
tubular structures, for instance on retinal images and on 3D tubular
structures on the Vessap dataset [19], a dataset of mice brain vascular
networks acquired at a very high resolution and with a research protocol.
However, such approach has not yet been tested on Magnetic Resonance (MR) or
X-ray Computed Tomography (CT) angiographic dataset acquired in clinical
conditions (i.e. images with more noise, artifacts and with a lower
resolution).

Figure 1: Architecture of the proposed cascaded U-Net (see Section 2.2).

Our work is inspired by Pan et al. [12] who proposed to perform the
skeletonization via a multitask architecture with shared features fusion
blocks which jointly perform centerline extraction and vessel segmentation.
Here, the skeleton is computed based on the initial image and benefits only
indirectly from the segmentation.

In this work, we propose a cascaded network with a U-Net backbone [13] that
first computes the segmentation, and then uses this segmentation and the
initial image to perform the skeletonization task. In this architecture, the
skeletonization tasks directly benefits from the segmentation, and may also
incorporate information from the initial image to produce a better skeleton.
Finally, the skeleton output, obtained in a differentiable manner, will be
used to compute the clDice loss and supervise the whole network, resulting in
a more topologically correct segmentation.

We evaluate our method against four standard U-Net, trained with the Dice or
clDice losses using the skeletonization methods introduced in [17] and [7]. We
show that our method provides segmentations with a more accurate topology and
with a lower computation time.

In this context, our main contributions are the following:

  * •

we evaluate the performance of the clDice loss for 3D vascular segmentation of
the human brain;

  * •

we propose an efficient way of performing the skeletonization operation to
compute the clDice;

  * •

we propose a cascaded multitask U-Net architecture which segments brain
vascular network with a more accurate topology and with a lower runtime.

The remainder of the article is organized as follows. In Section 2, we
describe our methodological contribution. In Section 3, we present experiments
that assess the relevance of our method. In Section 4, we discuss our results
and conclude.

##  2 Method

Table 1: Results of our skeletonization network vs. others skeletonization methods Model | Runtime (ms) $\downarrow$ |  $\chi$ error $\downarrow$ |  $\beta_{0}$ error $\downarrow$ |  $\beta_{1}$ error $\downarrow$
---|---|---|---|---
soft-skeleton algorithm |  $5$ $\pm$ $1$ |  $1313$ $\pm$ $257$ |  $1168$ $\pm$ $241$ |  $144$ $\pm$ $30$
Euler |  $558$ $\pm$ $13$ |  $1$ $\pm$ $1$ |  $2$ $\pm$ $2$ |  $1$ $\pm$ $1$
Boolean |  $1022$ $\pm$ $34$ |  $1$ $\pm$ $1$ |  $1$ $\pm$ $1$ |  $0$ $\pm$ $0$
Skeletonization network |  $9$ $\pm$ $2$ |  $412$ $\pm$ $69$ |  $294$ $\pm$ $48$ |  $118$ $\pm$ $26$
Table 2: Results of the Cascaded Multitask U-Net: mean $\pm$ standard deviation values ($p$-value). The $p$-values are computed with respect to the U-Net (Dice) results. Model | $DSC\uparrow$ | $clDice\uparrow$ |  $\chi$ error $\downarrow$ |  $\beta_{0}$ error $\downarrow$ |  $\beta_{1}$ error $\downarrow$ | Runtime (s) $\downarrow$
---|---|---|---|---|---|---
U-Net (Dice) |  0.76 $\pm$ $0.02$ |  0.85 $\pm$ $0.02$ |  $37.8$ $\pm$ $20.3$ |  $9.5$ $\pm$ $6.5$ |  $33.3$ $\pm$ $24.4$ |
U-Net (Dice + clDice Soft) |  $0.75$ $\pm$ $0.02$ |  0.85 $\pm$ $0.02$ |  $33.7$ $\pm$ $19.3$ |  $8.2$ $\pm$ $6.8$ |  $30.0$ $\pm$ $22.6$ |  $128$ $\pm$ $1$
U-Net (Dice + clDice Euler) |  $0.75$ $\pm$ $0.02$ |  0.85 $\pm$ $0.02$ |  $26.6$ $\pm$ $15.2$ |  $7.8$ $\pm$ $6.0$ |  25.0 $\pm$ $17.4$ |  $198$ $\pm$ $1$
U-Net (Dice + clDice Boolean) |  0.76 $\pm$ $0.01$ |  0.85 $\pm$ $0.02$ |  $24.5$ $\pm$ $15.9$ |  7.0 $\pm$ $6.7$ |  $25.8$ $\pm$ $16.8$ |  $301$ $\pm$ $1$
Cascaded U-Net (ours) |  $0.75$ $\pm$ $0.02$ |  $0.84$ $\pm$ $0.02$ |  23.5 $\pm$ $14.8$ |  7.0 $\pm$ $6.7$ |  $25.5$ $\pm$ $16.5$ |  $170$ $\pm$ $1$
$p$-value w.r.t U-Net (Dice) | $p<0.05$ | $p=0.3$ | $p<0.05$ | $p<0.05$ | $p<0.05$ |

###  2.1 clDice loss and differentiable skeletonizations algorithms

The clDice [17] derives from two metrics called topology precision
($T_{prec}$) and topology sensitivity ($T_{sens}$) in reference to the usual
precision and sensitivity metrics. These metrics are defined as follows:

$T_{prec}(C_{P},S_{G})=\frac{|C_{P}\cap S_{G}|}{|C_{P}|}\quad
T_{sens}(C_{G},S_{P})=\frac{|C_{G}\cap S_{P}|}{|C_{G}|}$

where $C_{P}$, $C_{G}$ and $S_{P}$, $S_{G}$ are the predicted and ground truth
centerlines and segmentations, respectively. The clDice is defined as the
harmonic mean of $T_{prec}$ and $T_{sens}$:

$clDice(S_{P},S_{G},C_{P},C_{G})=2\cdot\frac{T_{prec}(C_{P},S_{G})\cdot
T_{sens}(C_{G},S_{P})}{T_{prec}(C_{P},S_{G})+T_{sens}(C_{G},S_{P})}$

By leveraging the skeleton representation, the clDice avoids being biased by
large vessels and thus better focuses on topological information. Most of the
methods designed to extract a skeleton are not differentiable. Therefore, Shit
et al. proposed a differentiable soft-skeleton algorithm to use the clDice for
training a neural network. This algorithm uses min and max filters to perform
dilation and erosion on the predicted segmentation. Preliminary experiments
(see Section 3.2) showed that the results from this soft-skeletonization are
not sufficiently accurate for 3D vascular segmentation, in particular
regarding topology.

In a following work, Menten et al. [7] proposed two new differentiable
skeletonization algorithms. These algorithms remove simple points [14] in the
image, ensuring that the topology is not affected by the skeletonization. The
identification of simple points is done either using the Euler characteristics
or through a Boolean characterization. In the following, we will then refer to
these methods as Euler and Boolean methods, respectively. These methods have
the advantage of generating a nearly perfect skeleton in terms of topology,
but at the cost of an important computation time.

In our work, we chose to use a standard U-Net to perform the skeletonization.
This method is, by nature, differentiable and presents a good trade-off
between accuracy and computation time.

###  2.2 Model architecture

The backbone model used is a standard U-Net [13] with a depth of $4$, using
$2$-stride convolution for down-sampling, instance normalization and leakyReLU
activation function.

Our cascaded U-Net architecture is presented in Figure 1. It is composed of a
first U-Net taking as input an MRA image and performing the segmentation. This
task is supervised by a Dice loss and will be referred to as the segmentation
network. The output of this network is concatenated with the MRA image and fed
to a second U-Net performing the skeletonization task, also supervised by a
Dice loss. This part of the architecture will be referred to as the
skeletonization network. Similarly to Pan et al. [12], the training of these
two networks is also supervised by the clDice loss. The cascaded U-Net final
loss is then defined by:

$Loss(S_{P},S_{G},C_{P},C_{G})=Dice(S_{P},S_{G})+\\\ \lambda_{1}\cdot
Dice(C_{P},C_{G})+\lambda_{2}\cdot clDice(S_{P},S_{G},C_{P},C_{G})$ (1)

where $\lambda_{1},\lambda_{2}\in\mathbb{R}$ are two weight parameters.

This architecture presents several advantages. First, the skeletonization
network performs the skeletonization in a differentiable manner, which allows
using the clDice loss and enforcing topological constraint on the segmentation
task. In addition, unlike the other skeletonization methods, the
skeletonization network, that takes as input both the output of the
segmentation network and the original MRA image, can correct segmentation
errors and therefore limits error propagation. Finally, by jointly learning
segmentation and skeletonization, the segmentation network will benefit of an
inductive bias that will encourage learning topologically correct
segmentations. Indeed, the skeletonization task gives similar importance to
all vessels independently of their thickness; so the multitask learning can
help to enforce the importance of small vessels in the segmentation task.

###  2.3 Training configuration

All the MRA volumes were first normalized by Z-score. During training, one
batch is composed of $2$ patches of size $192\times 192\times 64$, each
randomly located in an MRA volume. One epoch consists of $250$ batches.

To ensure a fair comparison across our experiments, we used the same data
augmentation strategy for all trained networks, which is directly inspired by
nnUnet [3].

Sliding windows with $25\%$ overlap were used to reconstruct the full volume
at inference time. A Gaussian kernel was applied on overlapping parts to
reduce the weight of the voxels far from the center of the patches. We used
stochastic gradient descent with Nesterov momentum with an initial learning
rate set to $0.01$. A linear learning rate decay was applied so that learning
rate was equal to $0$ at the last epoch.

For the cascaded multitask U-Net, we chose to first pretrain the two U-Nets
separately before fine-tuning the network in the scheme described in Section
2.2. We initialized the weights of the segmentation network with the weights
of the U-Net trained with the Dice loss, and we pretrained the skeletonization
network using as input the MRA image and the segmentation ground truth.
Baseline models were trained during 750 epochs; for the cascaded multitask
U-Net we first pretrained the two networks during $750$ epochs, then fine-
tuned it during $250$ epochs.

##  3 Evaluation

In this section, we present the setup and results of the experiments conducted
to evaluate our skeletonization network and cascaded multitask U-Net. All the
results presented for deep-learning models were obtained through a $5$-fold
cross-validation.

Figure 2: Segmentation results obtained with the different approaches for
patient 8. Color boxes indicate the corresponding zoomed area.

###  3.1 Dataset and metrics

Dataset: For this study, we used a publicly available
dataset111https://public.kitware.com/Wiki/TubeTK/Data [1] containing $34$
time-of-flight (TOF) MRA volumes of the brain. All the volumes present a voxel
resolution of $0.513\times 0.513\times 0.800$ mm3 and a size of $448\times
448\times 128$. For all the volumes, the associated vessel segmentation ground
truth was available. To produce the skeleton ground truth, we used the
skeletonize function from the python package scikit-image.

Metrics: In the evaluation, we used clDice (see Section 2.1) and Dice
similarity coefficient ($DSC$). Beyond these quantitative metrics, we also
evaluated the topological quality of our segmentation by computing the
following topological scores: mean absolute error of the first two Betti
numbers $\beta_{0}$ (number of connected components) and $\beta_{1}$ (number
of tunnels), and mean absolute error of the Euler characteristics
$\chi=\sum_{i}(-1)^{i}\beta_{i}$ with $\beta_{i}$ the successive Betti
numbers. Before computing the metrics, we applied a post-processing on all the
results to remove the connected components with a volume lower than $100$
voxels.

###  3.2 Skeletonization using U-Net

To evaluate our skeletonization network, we compared it to the soft-skeleton
algorithm introduced in [17], the skeletonization based on the Boolean
characterization of simple points and the skeletonization based on the Euler
Characteristic introduced in [7]. We computed, for all methods, the mean time
required to perform the skeletonization per patch as well as the topological
metrics introduced in Section 3.1. These results, presented in Table 1, show
that our method produces better skeletons than the soft-skeleton algorithm
with regard to the topological metrics in a comparable computation time. Euler
and Boolean methods produce nearly perfect skeletons but at the cost of a more
important runtime. Our simple skeletonization network then provides a good
trade-off between accuracy and computation time. Qualitatively, the skeletons
produced by the soft-skeleton algorithm present many disconnections and a
thickness of several voxels compared to other skeletonization methods. By
contrast, our method provides more accurate skeletons.

###  3.3 Cascaded Multitask U-Net

####  3.3.1 Hyperparameters optimization

The goal of the cascaded multitask U-Net is to improve the results of the
segmentation network thanks to the clDice loss. As stated in Section 2.2, we
have to tune the hyperparameters $\lambda_{1}$ and $\lambda_{2}$ in order to
handle the trade-off between the skeletonization loss and the clDice loss. We
also tested two training configurations: a first one where the skeletonization
network weights are frozen; a second where they are updated during the
cascaded U-Net training. In both configurations, the skeletonization network
was first pre-trained. Based on these experiments, we found that the best
cascaded multitask U-Net training policy consists of freezing the weights of
the skeletonization network and set the loss weights to
$(\lambda_{1},\lambda_{2})=(0.5,0.5)$.

####  3.3.2 Validation and comparisons

We compared the segmentation of our cascaded U-Net with four methods:

  * •

a U-Net trained with a Dice loss;

  * •

a U-Net trained with the combination of Dice loss and clDice loss but with the
soft-skeleton, Euler or Boolean algorithms.

For all these methods, the trade-off between Dice and clDice loss was
optimized through a grid search.

Table 2 presents the results of the different methods. All methods present
similar results in terms of Dice and clDice. However, we can observe that
every method using the clDice loss succeeds to improve the topological
correctness of the segmentations, as indicated by the improvement of the
topological metrics. This behavior is also illustrated in Figure 2.

Besides, we can see that the model using soft-skeleton algorithm presents the
worst topological metrics, confirming that the quality of the skeletonization
is an important factor to use the clDice loss. Our method succeeds to give a
similar topological improvement with respect to the models using the Euler and
Boolean skeletonization methods but with an improved runtime (average epoch
duration). We attribute the good behavior of our method—despite an imperfect
skeletonization—to two factors. First, unlike deterministic methods, a learned
skeletonization network can correct segmentation error during training
(notably using the information of the MRA image) and therefore limits error
propagation. Second, in our method, the skeletonization acts as a
complementary task. Therefore, thanks to this joint multitask learning, the
segmentation task can benefit of an inductive bias that helps to learn a
topologically correct segmentation of its own.

##  4 Conclusion

In this article, we proposed to use a U-Net to learn the skeletonization
operation required to compute the clDice loss. This method provides a good
trade-off between the topological correctness of the skeletons and the
computation time. We then proposed a cascaded multitask U-Net to learn a
vascular segmentation with a topological constraint enforced by the clDice
loss. This cascaded multitask U-Net jointly learns vessel segmentation and
skeletonization and can then benefit of the inductive bias induced by the
skeletonization task.

In this study, we showed that the clDice loss improves the topological
correctness of vascular segmentation from MRA images. Moreover, we
demonstrated that our method gives the best topological correctness with the
least computation time. Our method may certainly be improved by investigating
dedicated skeletonization networks or improved information sharing between the
segmentation and skeletonization tasks. Such improvements will be part of
future works.

##  5 Compliance with ethical standards

This research study was conducted retrospectively using human subject data
made available in open access by Kitware at the following link:
https://public.kitware.com/Wiki/TubeTK/Data. Ethical approval was not required
as confirmed by the license attached with the open access data

##  6 Acknowledgements and conflicts of interests

This work was supported by the French Agence Nationale de la Recherche (Grant
ANR-20-CE45-0011). This work was granted access to the HPC resources of IDRIS
under the allocation 2022-AD011013610 made by GENCI.

The authors have no relevant financial or non-financial interests to disclose.

## References

  * [1] Aylward, S.R., Bullitt, E.: Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction. IEEE Transactions on Medical Imaging 21, 61–75 (2002)
  * [2] Dang, V.N., Galati, F., Cortese, R., Di Giacomo, G., Marconetto, V., Mathur, P., Lekadir, K., Lorenzi, M., Prados, F., Zuluaga, M.A.: Vessel-CAPTCHA: An efficient learning framework for vessel annotation and segmentation. Medical Image Analysis 75, 102263 (2022)
  * [3] Isensee, F., Jaeger, P.F., Kohl, S.A.A., Petersen, J., Maier-Hein, K.H.: nnU-Net: A self-configuring method for deep learning-based biomedical image segmentation. Nature methods 18, 203–211 (2021)
  * [4] Kitrungrotsakul, T., Han, X., Iwamoto, Y., Lin, L., Foruzan, A.H., Xiong, W., Chen, Y.: VesselNet: A deep convolutional neural network with multi pathways for robust hepatic vessel segmentation. Computerized Medical Imaging and Graphics 75, 74–83 (2019)
  * [5] Lesage, D., Angelini, E.D., Bloch, I., Funka-Lea, G.: A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes. Medical Image Analysis 13, 819–845 (2009)
  * [6] Livne, M., Rieger, J., Aydin, O.U., Taha, A.A., Akay, E.M., Kossen, T., Sobesky, J., Kelleher, J.D., Hildebrand, K., Frey, D., Madai, V.I.: A U-Net deep learning framework for high performance vessel segmentation in patients with cerebrovascular disease. Frontiers in Neuroscience 13, 97 (2019)
  * [7] Menten, M.J., Paetzold, J.C., Zimmer, V.A., Shit, S., Ezhov, I., Holland, R., Probst, M., Schnabel, J.A., Rueckert, D.: A skeletonization algorithm for gradient-based optimization. arXiv preprint arXiv:2309.02527 (2023)
  * [8] Miraucourt, O., Salmon, S., Szopos, M., Thiriet, M.: Blood flow in the cerebral venous system: Modeling and simulation. Computer Methods in Biomechanics and Biomedical Engineering 20, 471–482 (2017)
  * [9] Moccia, S., Momi, E.D., Hadji, S.E., Mattos, L.S.: Blood vessel segmentation algorithms - Review of methods, datasets and evaluation metrics. Computer Methods and Programs in Biomedicine 158, 71–91 (2018)
  * [10] Mookiah, M.R.K., Hogg, S., MacGillivray, T.J., Prathiba, V., Pradeepa, R., Mohan, V., Anjana, R.M., Doney, A.S., Palmer, C.N.A., Trucco, E.: A review of machine learning methods for retinal blood vessel segmentation and artery/vein classification. Medical Image Analysis 68, 101905 (2021)
  * [11] Mou, L., Zhao, Y., Fu, H., Liu, Y., Cheng, J., Zheng, Y., Su, P., Yang, J., Chen, L., Frangi, A.F., Akiba, M., Liu, J.: CS2-Net: Deep learning segmentation of curvilinear structures in medical imaging. Medical Image Analysis 67, 101874 (2021)
  * [12] Pan, L., Zhang, Z., Zheng, S., Huang, L.: MSC-Net: Multitask learning network for retinal vessel segmentation and centerline extraction. Applied Sciences 12, 403 (2021)
  * [13] Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional networks for biomedical image segmentation. In: MICCAI, Proceedings. pp. 234–241 (2015)
  * [14] Saha, P.K., Strand, R., Borgefors, G.: Digital topology and geometry in medical imaging: A survey. IEEE Transactions on Medical Imaging 34(9), 1940–1964 (2015)
  * [15] Sanches, P., Meyer, C., Vigon, V., Naegel, B.: Cerebrovascular network segmentation of MRA images with deep learning. In: ISBI, Proceedings. pp. 768–771 (2019)
  * [16] Shen, D., Wu, G., Suk, H.I.: Deep learning in medical image analysis. Annual Review of Biomedical Engineering 19, 221–248 (2017)
  * [17] Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P.W., Bauer, U., Menze, B.H.: clDice - A novel topology-preserving loss function for tubular structure segmentation. In: CVPR, Proceedings. pp. 16560–16569 (2021)
  * [18] Tetteh, G., Efremov, V., Forkert, N.D., Schneider, M., Kirschke, J., Weber, B., Zimmer, C., Piraud, M., Menze, B.H.: Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-D angiographic volumes. Frontiers in Neuroscience 14, 592352 (2020)
  * [19] Todorov, M.I., Paetzold, J.C., Schoppe, O., Tetteh, G., Shit, S., Efremov, V., Todorov-Völgyi, K., Düring, M., Dichgans, M., Piraud, M., Menze, B., Ertürk, A.: Machine learning analysis of whole mouse brain vasculature. Nature methods 17, 442–449 (2020)


# A Lagrangian perspective on the stability of ideal MHD equilibria with flow

Yao Zhou  Plasma Physics Laboratory, Princeton University, Princeton, New
Jersey 08543, USA  J. W. Burby  Plasma Physics Laboratory, Princeton
University, Princeton, New Jersey 08543, USA  Hong Qin  Plasma Physics
Laboratory, Princeton University, Princeton, New Jersey 08543, USA  Department
of Modern Physics, University of Science and Technology of China, Hefei, Anhui
230026, China

###### Abstract

We take a careful look at two approaches to deriving stability criteria for
ideal MHD equilibria. One is based on a tedious analysis of the linearized
equations of motion, while the other examines the second variation of the MHD
Hamiltonian computed with proper variational constraints. For equilibria
without flow, the two approaches are known to be fully consistent. However,
for equilibria with flow, the stability criterion obtained from the
constrained variation approach was claimed to be stronger than that derived
using the linearized equations of motion. We show this claim is incorrect by
deriving and comparing both criteria within the same framework. It turns out
that the criterion obtained from the constrained variation approach has
stricter requirements on the initial perturbations than the other. Such
requirements naturally emerge in our new treatment of the constrained
variation approach using the Euler-Poincaré structure of ideal MHD, which is
more direct and simple than the previous derivation from the Poisson
perspective.

##  I Introduction

The stability properties of an ideal MHD equilibrium determine whether or not
such a configuration can be found in nature or reproduced in the laboratory.
Therefore, theoretical tools for assessing the stability of a given
equilibrium stand as a cornerstone in the theory of magnetized plasma physics.
For _absolute equilibria_ (zero flow), where all fluid elements are
stationary, a particularly well-developed tool is based on the energy
principleBernstein  _et al._ (1958); Freidberg (1987); Kulsrud (1964); Laval
_et al._ (1965), which states that an absolute equilibrium will be linearly
stable if and only if the perturbed potential energy
$\delta^{2}W_{0}(\bm{\xi},\bm{\xi})$ is non-negative definite as a functional
of the the fluid displacement $\bm{\xi}$. For _relative equilibria_ (non-zero
flow), where fluid elements move along a time-independent Eulerian velocity
field, tools for determining stability are not as well-developed. Existing
stability criteria only provide sufficient conditions for stability. The
oldest condition was found by Frieman and RotenbergFrieman and Rotenberg
(1960) from analyzing the linearized equations of motion. More recently,
HameiriHameiri (2003) derived a sufficient condition for stability against so-
called dyanmically accessible perturbations, which we will discuss at greater
length in what follows.

Methods for deriving these stability criteria fall into two categories. In the
first categoryBernstein  _et al._ (1958); Frieman and Rotenberg (1960);
Freidberg (1987) are those methods based on careful, direct analysis of the
linearized MHD equations of motion. These methods proceed by first proving
directly that the linearized force operator is symmetric, and then deriving a
stability criterion as a result of this property. The second category consists
of methods based on energy conservation Freidberg (1987); Kulsrud (1964);
Laval  _et al._ (1965); Almaguer  _et al._ (1988); Isichenko (1998); Hameiri
(2003, 1998). Because the non-linear MHD equations possess a conserved energy
functional, the linearized equations of motion inherit a conserved quadratic
functional. This fact can be exploited to indirectly prove that the linearized
force operator is symmetricKulsrud (1964). It also immediately leads to a
sufficient condition for linear stability because the conserved quadratic
functional provides a norm that bounds linear solutions for all time whenever
it is positive definiteHolm  _et al._ (1985).

When applied to absolute equilibria, either type of method eventually leads to
the energy principle, which, again, gives a necessary and sufficient condition
for linear stability. However, neither approach gives a necessary and
sufficient condition for the linear stability of relative equilibria; the
power of each type of method is reduced to only providing sufficient stability
criteria. For example, Frieman and Rotenberg’s conditionFrieman and Rotenberg
(1960) asserts that a relative equilibrium will be linearly stable if
$\delta^{2}W(\bm{\xi},\bm{\xi})=-\int\bm{\xi}\cdot\mathbf{F}(\bm{\xi})\,\mathrm{d}^{3}x$
is positive definite. Here, $\mathbf{F}$ is the force operator appearing in
the ideal MHD equations linearized about a flowing equilibrium,

$\displaystyle\rho\ddot{\bm{\xi}}+2\rho\mathbf{v}\cdot\nabla\dot{\bm{\xi}}=\mathbf{F}(\bm{\xi}),$
(1)

where $\rho$ is the mass density and $\mathbf{v}$ is the Eulerian velocity
field of the equilibrium. The root of the difficulty in proving _necessity_ of
this condition, a difficulty which is not present when analyzing absolute
equilibria, is the presence of a first time derivative of $\bm{\xi}$ in Eq.
(1).

Another interesting subtlety ushered in by the presence of a non-zero
equilibrium flow is that the stability criteria derived from the two
approaches no longer coincide. Hameiri’s energy-based methodHameiri (2003) for
deriving a condition for linear stability of relative equilibria against
dynamically accessible perturbations by leveraging the Poisson structure of
ideal MHDMorrison and Greene (1980), which is a refined version of Isichenko’s
approachIsichenko (1998), achieves a different result from the much earlier
condition derived by Frieman and Rotenberg. The key idea is that the
degeneracy of the field-theoretic Poisson bracket implies that the time
derivative of a solution to the non-linear MHD equations cannot point in an
arbitrary direction in the (infinite-dimensional) MHD phase space; instead the
time-derivative must point along the level sets of the Casimir functionals.
The so-called dynamically accessible variationsMorrison (1998); Morrison  _et
al._ (1989); Hameiri (2003) span these special directions in the phase space.
Using this idea, Hameiri argued that a sufficient condition for linearized
stability against dynamically accessible variations is the positive
definiteness of the functional

$\displaystyle\delta^{2}H(\bm{\xi},\beta,\alpha,\bm{\zeta})=$ (2)
$\displaystyle\int[\rho(\delta\mathbf{v}(\bm{\xi},\beta,\alpha,\bm{\zeta})-\mathbf{v}\cdot\nabla\bm{\xi}+\bm{\xi}\cdot\nabla\mathbf{v})^{2}-\bm{\xi}\cdot\mathbf{F}(\bm{\xi})]\,\mathrm{d}^{3}x,$

where

$\displaystyle\delta\mathbf{v}(\bm{\xi},\beta,$
$\displaystyle\alpha,\bm{\zeta})=$ (3)
$\displaystyle\bm{\xi}\times(\nabla\times\mathbf{v})+{(\nabla\times\bm{\zeta})\times\mathbf{B}}/{\rho}+\alpha\nabla
s-\nabla\beta.$

Here, $(\bm{\xi},\beta,\alpha,\bm{\zeta})$ are known as Clebsch
variablesMarsden  _et al._ (1984); Morrison (1998) and they parameterize the
dynamically accessible variations. Hameiri noticed that the positive
definiteness of $\delta^{2}W$ is sufficient but not necessary for that of
$\delta^{2}H$. He then concluded that his stability criterion is stronger than
Frieman and Rotenberg’s.

The proof of the validity of Hameiri’s stability criterion comes from the
general theory of dynamically accessible variationsMorrison (1998); Morrison
_et al._ (1989); Hameiri (2003), which applies to general Poisson dynamical
systems. This theory gives a powerful method for finding stability criteria in
various finite and infinite dimensional Hamiltonian systems, and its
importance cannot be overstated. However, arguments along these lines, and
Hameiri’s argument in particular, may be difficult to grasp without a
familiarity with infinite dimensional Poisson dynamics. Moreover, the physical
origin of the Casimir invariants, which plays a critical role in deriving
dynamically accessible stability criteria, may seem mysterious when one starts
from a generic Poisson formulation. In particular, an obvious relationship
between these invariants and symmetry properties of the underlying physical
system is lacking.

In this paper, we will give a new energy-based derivation of Hameiri’s
stability criterion for relative ideal MHD equilibria that is more direct and
tangible than previous derivations. We will do this by way of an approach
based on an Euler-Poincaré Holm  _et al._ (1998, 2009) counterpart to the
method of dynamically accessible variations which we developed for this
purpose. Then we will derive Frieman and Rotenberg’s criterion within the same
framework and compare the two stability criteria. This comparison will yield a
subtle surprise which could also be deduced from Morrison’s description of
dynamically accessible variations on the Poisson sideMorrison (1998); Morrison
_et al._ (1989); Hameiri’s criterion requires the initial perturbation to
preserve the Noether invariants implied by the particle relabeling symmetry
Cotter and Holm (2012), while Frieman and Rotenberg’s does not. This means
Frieman and Rotenberg’s stability criterion is not strictly weaker than
Hameiri’s. Thus, the significance of this work is threefold; we provide a
simplified means for understanding Hameiri’s stability criterion; we outline a
method for deriving dynamically accessible stability criteria for general
Euler-Poincaré fluids; and we clarify that Hameiri’s stability criterion
applies only to dynamically accessible perturbations while Frieman and
Rotenberg’s criterion applies to a more general class of perturbations.

We will begin with two sections describing the Lagrangian counterpart to the
method of dynamically accessible variations that we have developed. This
general, abstract theory, while illuminating, is not strictly necessary to
follow our simplified derivation for both stability criteria, and therefore
may be skipped by those readers only interested in the difference between the
two criteria. In Sec. II, following Ref. Holm  _et al._ , 1998, we will
briefly review the notion of an Euler-Poincaré fluid. In Sec. III, we will
study perturbations to relative equilibria that preserve the Noether
invariants implied by the particle relabeling symmetry in a general Euler-
Poincaré fluid. In particular, we will demonstrate (a) how to derive Clebsch
variables that parameterize these perturbations and (b) the existence of a
quadratic functional conserved under the evolution of these perturbations.
Then in Sec. IV we will derive both stability criteria based on the Euler-
Poincaré theory. Finally, we discuss the relationship between our derivation
and previous ones, as well as the difference between Hameiri’s criterion and
Frieman and Rotenberg’s in Sec. V.

##  II Euler-Poincaré Fluids

In this section, we briefly review the general notion of an Euler-Poincaré
fluid in order to prepare for a discussion of linearized stability in this
setting in the next section. The material in this section comes directly from
the definitive Refs. Holm  _et al._ , 1998 and Holm  _et al._ , 2009, which
develop the theory of Euler-Poincaré dynamical systems in great detail. For a
detailed explanation of how it applies to the Maxwell-Vlasov system see Ref.
Cendra  _et al._ , 1998. Also see Ref. Squire  _et al._ , 2013 for an
application of the theory to gyrokinetics. As a heuristic rule, conservative
fluid systems that exhibit some sort of particle-relabeling symmetry tend to
exhibit this type of structure. In order to understand our proof of previous
stability criteria later in the paper, knowledge of this formalism is not
strictly necessary. However, knowledge of the theory illuminates the structure
of our proof.

Consider a fluid confined within a three-dimensional domain
$\mathcal{D}\subset\mathbb{R}^{3}$. Label points in $\mathcal{D}$ with the
variable $X\in\mathcal{D}$ and let the time-dependent mapping
$\mathbf{x}_{t}:\mathcal{D}\rightarrow\mathcal{D}$ give the time-$t$ position
of a fluid element $\mathbf{x}_{t}(X)$ assuming it began at $X$ when $t=0$. If
we assume, as we will, that the mapping $\mathbf{x}_{t}$ is smooth with smooth
inverse $\mathbf{x}_{t}^{-1}$ for all $t$ (note that this may not always be
true due to phenomena such as cavitation and shock formation), then we may
regard the configuration of the fluid as being specified by an element of the
diffeomorphism group $Q\equiv\text{Diff}(\mathcal{D})$, which is the set of
all smooth mappings $\mathcal{D}\rightarrow\mathcal{D}$ with smooth inverses.

Let $V=\Omega_{k}(\mathcal{D})$ be the space of differential $k$-forms over
$\mathcal{D}$ with typical element $u$, i.e. the space of degree-$k$ covariant
antisymmetric tensor fields. Let $V^{*}=\Omega_{3-k}(\mathcal{D})$ be the
$(3-k)$-forms with typical element $a$. Finally, let
$\mathfrak{X}(\mathcal{D})$ be the space of vector fields on $\mathcal{D}$
tangent to $\mathcal{D}$’s boundary with typical element $\mathbf{v}$. Given a
Lagrangian function $L:Q\times\mathfrak{X}(\mathcal{D})\times
V^{*}\rightarrow\mathbb{R}$ satisfying the symmetry property Note (1)
$L(\mathbf{x}\circ\mathbf{y},\mathbf{v},\mathbf{y}^{*}a)=L(\mathbf{x},\mathbf{v},a)$
for all $\mathbf{y}\in\text{Diff}(\mathcal{D})$, the associated _Euler-
Poincaré fluid with advected parameter_ $a_{o}\in V^{*}$ is a fluid whose
paths in configuration space $t\mapsto\mathbf{x}_{t}$ are extremals of the
action functional

$\displaystyle\mathcal{A}(\mathbf{x}_{\cdot})=\int_{t_{1}}^{t_{2}}L\left(\mathbf{x}_{t},\dot{\mathbf{x}}_{t}\circ\mathbf{x}_{t}^{-1},a_{o}\right)\mathrm{d}t$
(4)

over the space of curves in $Q$ with fixed endpoints. Note that
$\mathbf{x}_{\cdot}$ denotes the curve $t\mapsto\mathbf{x}_{t}$, which is
distinct from $\mathbf{x}_{t}\in Q$.

This variational principle leads to a second-order differential equation
governing the evolution of the fluid configuration $\mathbf{x}_{t}\in Q$.
Because $\mathbf{x}_{t}$ gives the position of every fluid element in
$\mathcal{D}$ at time $t$, this is a Lagrangian, as opposed to Eulerian
description of the fluid (unfortunately, we have a Lagrangian Lagrangian
description of our fluid). The great utility of the assumed structure of an
Euler-Poincaré fluid, and in particular the symmetry property possessed by the
Lagrangian $L$, is the systematic method it offers for passing to an Eulerian
description of the fluid, which we now describe.

Define the projection map $\pi:Q\times\mathfrak{X}(\mathcal{D})\times
V^{*}\rightarrow\mathfrak{X}(\mathcal{D})\times V^{*}$ using the formula
$\pi(\mathbf{x},\mathbf{v},a)=(\mathbf{v},\mathbf{x}^{-1*}a)$. Due to the
symmetry property
$L(\mathbf{x}\circ\mathbf{y},\mathbf{v},\mathbf{y}^{*}a)=L(\mathbf{x},\mathbf{v},a)$,
there is a unique and well-defined functional
$l:\mathfrak{X}(\mathcal{D})\times V^{*}\rightarrow\mathbb{R}$ satisfying
$l\circ\pi=L$, i.e.

$\displaystyle l(\mathbf{v},a)=L(\mathbf{x},\mathbf{v},\mathbf{x}^{*}a).$ (5)

Using this _reduced Lagrangian_ , we can form a constrained variational
principle on $\mathfrak{X}(\mathcal{D})\times V^{*}$ equivalent to the
unconstrained variational principle on $Q\times\mathfrak{X}(\mathcal{D})\times
V^{*}$ defined by the action functional $\mathcal{A}$. As proved in Ref. Holm
_et al._ , 1998, the curve $\mathbf{x}_{\cdot}$ with
$\mathbf{x}_{t_{1}}=e_{1}$ and $\mathbf{x}_{t_{2}}=e_{2}$ is an extremal of
the functional $\mathcal{A}$ if and only if it solves the initial value
problem
$\dot{\mathbf{x}}_{t}=\mathbf{v}_{t}\circ\mathbf{x}_{t};~{}\mathbf{x}_{t_{1}}=e_{1}$
for a $\mathbf{v}_{t}$ that solves the constrained variational problem

$\displaystyle\delta\int_{t_{1}}^{t_{2}}\\!l(\mathbf{v}_{t},a_{t})\,\mathrm{d}t=0.$
(6)

The precise definition of this constrained variational problem is as follows.
Let $\Gamma_{t}=(u_{t},a_{t})\in\mathfrak{X}(\mathcal{D})\times V^{*}$, where
$a_{t}$ is the unique solution to the initial value problem
$\dot{a}_{t}=-\mathfrak{L}_{\mathbf{v}_{t}}a_{t};~{}a_{t_{1}}=e_{1}^{-1*}a_{o}$.
Let $\Gamma_{t,\epsilon}$ define a two parameter curve in
$\mathfrak{X}(\mathcal{D})\times V^{*}$ with $\Gamma_{t,0}=\Gamma_{t}$,
$\Gamma_{t_{1},\epsilon}=(\mathbf{v}_{t_{1}},a_{t_{1}})$,
$\Gamma_{t_{2},\epsilon}=(\mathbf{v}_{t_{2}},a_{t_{2}})$, and whose variation
has the form

$\displaystyle\delta\Gamma_{t}\equiv\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\Gamma_{t,\epsilon}=(\dot{\bm{\xi}}_{t}+\mathfrak{L}_{\mathbf{v}_{t}}\bm{\xi}_{t},-\mathfrak{L}_{\bm{\xi}_{t}}a_{t}),$
(7)

for an arbitrary time dependent vector field
$\bm{\xi}_{t}\in\mathfrak{X}(\mathcal{D})$ with
$\bm{\xi}_{t_{1}}=\bm{\xi}_{t_{2}}=0$. Then $\mathbf{v}_{t}$ is a solution to
the constrained variational problem (6) if

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\int_{t_{1}}^{t_{2}}\\!l(\Gamma_{t,\epsilon})\,\mathrm{d}t=0$
(8)

for any choice of the time-dependent vector field $\bm{\xi}_{t}$.

From this new perspective, the constrained variational principle determines
the dynamics of the Eulerian velocity field $\mathbf{v}_{t}$ and the detailed
fluid configuration $\mathbf{x}_{t}$ follows as a sort of afterthought; if we
wanted to know how each fluid element moves, all that would be necessary is
integration of the time-dependent ODE
$\dot{\mathbf{x}}_{t}=\mathbf{v}_{t}\circ\mathbf{x}_{t}$. Therefore, we have
arrived at an Eulerian description of our fluid; whereas before the fluid
configuration $\mathbf{x}_{t}$ was regarded as the independent variable, now
$\mathbf{v}_{t}$ is the independent variable. In fact, we can write down the
PDE, known as the Euler-Poincaré equations, that must be satisfied by
$\mathbf{v}_{t}$ explicitly. We have Holm  _et al._ (1998, 2009)

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\frac{\delta
l}{\delta\mathbf{v}}=-\mathfrak{L}_{\mathbf{v}_{t}}\frac{\delta
l}{\delta\mathbf{v}}+\frac{\delta l}{\delta a}\diamond a_{t},$ (9)

which follows directly from the constrained variational principle. Also,
recall that the evolution of $a_{t}$ is assumed to be specified by the
equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}{a}_{t}=-\mathfrak{L}_{\mathbf{v}_{t}}a_{t}.$
(10)

Here the functional derivative $\delta
l/\delta\mathbf{v}\in\mathfrak{X}^{*}(\mathcal{D})$ is a one-form density,
while the diamond operator $\diamond$ obeys the formula

$\displaystyle\left\langle\frac{\delta l}{\delta a}\diamond
a_{t},\bm{\xi}\right\rangle=-\left\langle\mathfrak{L}_{\bm{\xi}}a_{t},\frac{\delta
l}{\delta a}\right\rangle.$ (11)

These notations are briefly explained in the appendix for the readers’
convenience. For more detailed discussions see Refs. Holm  _et al._ , 1998 and
Holm  _et al._ , 2009.

Euler-Poincaré fluids are a special case of the more general Euler-Poincaré
dynamical systems developed by Holm, Marsden, and Ratiu in Ref. Holm  _et al._
, 1998. In the general theory, one starts with a Lie group $G$, a vector space
$V$, a (left or right) representation of $G$ on $V$, and a (left or right)
$G$-invariant Lagrangian function $L:TG\times V^{*}\rightarrow\mathbb{R}$. $G$
is regarded as the configuration space for a mechanical system while $V^{*}$
is regarded as a space of parameters. Hamilton’s principle applied to the
action functional
$\mathcal{A}=\int_{t_{1}}^{t_{2}}\\!L(g,\dot{g},a_{o})\,\mathrm{d}t$ then
leads to the usual sort of Euler-Lagrange equations on $TG$, although they are
parameterized by the variable $a_{o}\in V^{*}$. Finally, by leveraging the
$G$-invariance of $L$, an equivalent dynamical description is obtained by
reducing Hamilton’s variational principle to a constrained variational
principle on $\mathfrak{g}\times V^{*}$, where $\mathfrak{g}$ is the Lie
algebra of $G$. The paths in $\mathfrak{g}\times V^{*}$ satisfying the
constrained variational principle obey the Euler-Poincaré equations. When
$G=\text{Diff}(\mathcal{D})$, $V$ is taken to be a space of differential
forms, and the (right) representation of $G$ on $V$ is given by pullback, the
Euler-Poincaré equations are precisely Eqs. (9) and (10).

When applied to continuum theories, the passage from $TG$ to $\mathfrak{g}$ is
precisely the passage from a Lagrangian to an Eulerian description. Thus,
Holm, Marsden, and Ratiu’s theory offers an attractive means for passing to an
Eulerian description of a continuum without losing sight of the variational
properties of the Lagrangian description. In particular, the symmetry
properties of the Lagrangian imply conservation laws for the Eulerian
equations via Noether’s theorem or the Kelvin Noether theorem, as discussed in
Refs. Cotter and Holm, 2012 and Holm  _et al._ , 1998. These conservation laws
will play a crucial role in the Lagrangian formulation of the method of
dynamically accessible variations that we describe and apply in the coming
sections. The distinction between these conservation laws and the advected
parameters is also key to understanding the subtle difference between the two
stability criteria, while in the previous works from the Poisson perspective
Hameiri (2003) all the constants of motion (Casimirs) were treated as the
same.

##  III Clebsch Variables and Energy Conservation

Now we will use the general Euler-Poincaré equations to study linear
perturbations about relative equilibria that preserve the Noether invariants
implied by the particle relabeling symmetry. By definition, a relative
equilibrium $(\mathbf{v}_{o},a_{o})$ is a stationary solution to the combined
system of PDE’s defined by equations (9) and (10), which we will refer to as
the Euler-Poincaré system. Moreover, as discussed in Ref. Cotter and Holm,
2012, the Noether invariants are given by

$\displaystyle\left<\frac{\delta
l}{\delta\mathbf{v}},\bm{\eta}_{t}\right>=const.$ (12)

where $\bm{\eta}_{t}$ is the unique solution to the initial value problem
$\dot{\bm{\eta}}_{t}=-\mathfrak{L}_{\mathbf{v}_{t}}\bm{\eta}_{t};~{}\bm{\eta}_{t_{1}}=\bm{\eta}_{o}$,
$\bm{\eta}_{o}$ is any vector field that generates an infinitesimal symmetry
of $a_{t_{1}}$, i.e. $\mathfrak{L}_{\bm{\eta}_{o}}a_{t_{1}}=0$, and $t_{1}$ is
the initial time, which we will assume to be zero henceforth. Therefore, we
would like to consider dynamical solutions to the Euler-Poincaré system that
start out infinitesimally close to $(\mathbf{v}_{o},a_{o})$ and that all have
the same values of the Noether quantities given in Eq. (12). As we will see
explicitly when applying the general theory we are about to develop to ideal
MHD, fixing the values of the Noether quantities is the Lagrangian analogue of
fixing the values of the Casimirs when working with Poisson systems. In fact,
these Noether invariants, together with the advected parameters which are
intrinsic to the Euler-Poincaré system, appear to be fully equivalent to the
Casimir invariants modulo a Legendre transformation, as will be discussed in a
forthcoming article. Thus, we are considering the behavior of what Morrison
has dubbed dynamically accessible variations in the Poisson context from the
Lagrangian perspective.

The goal of our analysis is twofold. First of all, we would like to find a
parameterization of these dynamically accessible linear perturbations that
reflects the fact that any solution to the Euler-Poincaré system on
$\mathfrak{X}(\mathcal{D})\times V^{*}$ is the projection of a solution to the
variational principle defined by $\mathcal{A}$. The parameterizing variables
are known as Clebsch variables. Second, we would like to find the equations of
motion for these Clebsch variables and an associated conservation law.

Let $\Gamma_{t,\epsilon}=(\mathbf{v}_{t,\epsilon},a_{t,\epsilon})$ be an
$\epsilon$-dependent family of solutions to the Euler-Poincaré system with
$\Gamma_{t,0}=\Gamma_{o}=(\mathbf{v}_{o},a_{o})$ the stationary solution. If
we define the operator
$\delta=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\big{|}_{0}$, then
$\delta\Gamma_{t}=(\delta\mathbf{v}_{t},\delta
a_{t})\in\mathfrak{X}(\mathcal{D})\times V^{*}$ is a linearized solution to
the Euler-Poincaré system. It must satisfy the linearized equations of motion:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\left(\delta\frac{\delta
l}{\delta\mathbf{v}}\right)$ $\displaystyle=$
$\displaystyle-\mathfrak{L}_{\delta\mathbf{v}_{t}}\frac{\delta
l}{\delta\mathbf{v}}-\mathfrak{L}_{\mathbf{v}_{o}}\left(\delta\frac{\delta
l}{\delta\mathbf{v}}\right)$ $\displaystyle+\left(\delta\frac{\delta l}{\delta
a}\right)\diamond a_{o}+\frac{\delta l}{\delta a}\diamond\delta a_{t}$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\delta a_{t}$ $\displaystyle=$
$\displaystyle-\mathfrak{L}_{\delta\mathbf{v}_{t}}a_{o}-\mathfrak{L}_{\mathbf{v}_{o}}\delta
a_{t}.$ (14)

In addition, $\delta\Gamma_{t}=(\delta\mathbf{v}_{t},\delta a_{t})$ must have
the form

$\displaystyle\delta\mathbf{v}_{t}$
$\displaystyle=\dot{\bm{\xi}}_{t}+\mathfrak{L}_{\mathbf{v}_{o}}\bm{\xi}_{t}$
(15) $\displaystyle\delta a_{t}$
$\displaystyle=-\mathfrak{L}_{\bm{\xi}_{t}}a_{o},$ (16)

for some time-dependent vector field $\bm{\xi}_{t}$. This follows from the
fact that for each $\epsilon$, $\Gamma_{t,\epsilon}$ must be the projection of
a solution to the original variational principle defined by $\mathcal{A}$.
That is,
$\Gamma_{t,\epsilon}=\pi(\mathbf{x}_{\epsilon,t},\dot{\mathbf{x}}_{\epsilon,t}\circ\mathbf{x}_{\epsilon,t}^{-1},a_{o})$
(incidentally, this also explains the form of the constrained variations in
Eq. (7)). Finally, because we assume that the Noether invariants (12)
evaluated on $\Gamma_{t,\epsilon}$ are independent of $t$ and $\epsilon$, we
have

$\displaystyle\left<\delta\frac{\delta
l}{\delta\mathbf{v}}+\mathfrak{L}_{\bm{\xi}_{t}}\frac{\delta
l}{\delta\mathbf{v}},\bm{\eta}\right>=0,$ (17)

for all vector fields $\bm{\eta}$ generating infinitesimal symmetries of
$a_{o}$, $\mathfrak{L}_{\bm{\eta}}a_{o}=0$.

These three conditions on $\delta\Gamma_{t}$ contain a great deal of
information, which we will now unravel. Combining Eq. (III) with Eq. (15) and
using both the identity $\mathfrak{L}_{\mathbf{v}_{o}}(u\diamond
a)=(\mathfrak{L}_{\mathbf{v}_{o}}u)\diamond
a+u\diamond(\mathfrak{L}_{\mathbf{v}_{o}}a)$ and the equilibrium condition
$\mathfrak{L}_{\mathbf{v}_{o}}\frac{\delta l}{\delta\mathbf{v}}=\frac{\delta
l}{\delta a}\diamond a_{o}$, we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mu_{t}+\mathfrak{L}_{\mathbf{v}_{o}}\mu_{t}=\omega_{t}\diamond
a_{o},$ (18)

where $\mu_{t}=\delta\frac{\delta
l}{\delta\mathbf{v}}+\mathfrak{L}_{\bm{\xi}_{t}}\frac{\delta
l}{\delta\mathbf{v}}$ and $\omega_{t}=\delta\frac{\delta l}{\delta
a}+\mathfrak{L}_{\bm{\xi}_{t}}\frac{\delta l}{\delta a}$. This equation may be
regarded as a linear, inhomogeneous partial differential equation for the one-
form density $\mu_{t}$. It admits a simple formal solution in terms of the the
Lagrangian trajectory $\mathbf{y}_{t}\in Q$ of the equilibrium flow
$\mathbf{v}_{o}$. Indeed, using the fact that
$\dot{\mathbf{y}}_{t}=\mathbf{v}_{o}\circ\mathbf{y}_{t}$, the left-hand side
of Eq. (18) can be written
$\mathbf{y}_{t}^{-1*}\frac{\mathrm{d}}{\mathrm{d}t}(\mathbf{y}_{t}^{*}\mu_{t})$,
which implies

$\displaystyle\mu_{t}=\mu_{0}+\left(\int_{0}^{t}\\!(\mathbf{y}_{\tau}\circ\mathbf{y}_{t}^{-1})^{*}\omega_{\tau}\,\mathrm{d}\tau\right)\diamond
a_{o}.$ (19)

Moreover, in light of equation (17), $\mu_{0}=u_{0}\diamond a_{o}$ for some
$u_{0}\in V$. Therefore, by Eq. (19), for all $t$, there is some $u_{t}\in V$
and some $\bm{\xi}_{t}\in\mathfrak{X}(\mathcal{D})$ such that Note (2)

$\displaystyle\delta\frac{\delta l}{\delta\mathbf{v}}$
$\displaystyle=-\mathfrak{L}_{\bm{\xi}_{t}}\frac{\delta
l}{\delta\mathbf{v}}+u_{t}\diamond a_{o}$ (20) $\displaystyle\delta a_{t}$
$\displaystyle=-\mathfrak{L}_{\bm{\xi}_{t}}a_{o},$ (21)

and

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}{\bm{\xi}_{t}}$ $\displaystyle=$
$\displaystyle-\mathfrak{L}_{\mathbf{v}_{o}}\bm{\xi}_{t}+\delta\mathbf{v}_{t},$
(22) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}u_{t}$ $\displaystyle=$
$\displaystyle-\mathfrak{L}_{\mathbf{v}_{o}}u_{t}+\delta\frac{\delta l}{\delta
a}+\mathfrak{L}_{\bm{\xi}_{t}}\frac{\delta l}{\delta a}.$ (23)

Provided that the Lagrangian is non-degenerate, Eqs. (20) and (21) tell us how
to complete the first goal of our analysis; these equations identify the
appropriate Clebsch variables Marsden and Ratiu (1999); Morrison (1998). The
idea behind Clebsch variables in general is to represent quantities subject to
constraints in a way that automatically satisfies the constraints. For
instance, a magnetic field $\mathbf{B}$ is subject to the constraint
$\nabla\cdot\mathbf{B}=0$. If $\mathbf{B}$ is expressed as
$\nabla\alpha\times\nabla\beta$, then it clearly automatically satisfies the
divergence-free constraint. Conversely, locally every magnetic field can be
represented in this way. In this case $\alpha,\beta$ are known as (local)
Clebsch variables. In the setting of this paper, Clebsch variables should give
a representation of linearized solutions to the Euler-Poincaré system that
guarantee they satisfy the constraint of preserving the Noether invariants (to
first order), i.e. Eq. (17). To see that Eqs. (20) and (21) do indeed identify
Clebsch variables, first it is necessary to define what it means for the
Lagrangian to be non-degenerate. $l$ is non-degenerate when the linear
operator
$\frac{\delta^{2}l}{\delta\mathbf{v}^{2}}(\mathbf{v}_{o},a_{o}):\mathfrak{X}(\mathcal{D})\rightarrow\mathfrak{X}(\mathcal{D})^{*}$
given by the formula

$\displaystyle\frac{\delta^{2}l}{\delta\mathbf{v}^{2}}(\mathbf{v}_{o},a_{o})[\mathbf{w}]=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v}_{o}+\epsilon\mathbf{w},a_{o})$ (24)

is invertible. When the Lagrangian is non-degenerate, Eqs. (20) and (21) can
be used to explicitly express $\delta\mathbf{v}_{t}$ and $\delta a_{t}$ in
terms of $\bm{\xi}_{t}$ and $u_{t}$ using the mapping
$\mathfrak{C}:\mathfrak{X}(\mathcal{D})\times
V\rightarrow\mathfrak{X}(\mathcal{D})\times V^{*}$ given by

$\displaystyle\mathfrak{C}(\bm{\xi},u)=\bigg{(}\frac{\delta^{2}l}{\delta\mathbf{v}^{2}}^{-1}\left[u\diamond
a_{o}-\mathfrak{L}_{\bm{\xi}}\frac{\delta
l}{\delta\mathbf{v}}+\frac{\delta^{2}l}{\delta
a\delta\mathbf{v}}[\mathfrak{L}_{\bm{\xi}}a_{o}]\right]$ (25)
$\displaystyle,-\mathfrak{L}_{\bm{\xi}}a_{o}\bigg{)},$

where $\frac{\delta^{2}l}{\delta
a\delta\mathbf{v}}[a^{\prime}]=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\big{|}_{0}\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v}_{o},a_{o}+\epsilon a^{\prime})$. Indeed, Eqs.
(20) and (21) can be written in terms of $\mathfrak{C}$ as

$\displaystyle(\delta\mathbf{v}_{t},\delta
a_{t})=\mathfrak{C}(\bm{\xi}_{t},u_{t}).$ (26)

Here $\bm{\xi}$ and $u$ are the Clebsch variables; any
$(\delta\mathbf{v},\delta a)$ expressed in the form $\mathfrak{C}(\bm{\xi},u)$
will automatically satisfy Eq. (17). Conversely, by the derivation of Eqs.
(20) and (21), any $(\delta\mathbf{v}_{t},\delta a_{t})$ that is a linearized
solution preserving the Noether invariants implied by the particle relabeling
symmetry must be of the form $\mathfrak{C}(\bm{\xi}_{t},u_{t})$.

Now we turn to the first half of the second goal of this section; specifying
appropriate dynamical equations for the Clebsch variables $\bm{\xi},u$.
Whatever these Clebsch variable evolution equations might be, they must have
two properties: (a) if $(\bm{\xi}_{t},u_{t})$ is a solution to the Clebsch
variable evolution equations, then $\mathfrak{C}(\bm{\xi}_{t},u_{t})$ must be
a solution to the linearized Euler-Poincaré system (b) every linearized
solution of the Euler-Poincaré system that preserves the Noether invariants
must be given by applying $\mathfrak{C}$ to a solution to the Clebsch variable
evolution equations. To find such equations, we only have to examine Eqs. (22)
and (23). If every instance of $\delta\mathbf{v}_{t}$ and $\delta a_{t}$
appearing in these equations is expressed in terms of $\bm{\xi}_{t}$ and
$u_{t}$ using $\mathfrak{C}$, then it becomes clear that these equations
define a closed system of time evolution equations for $\bm{\xi}_{t}$ and
$u_{t}$. Moreover, it is easy to check that solutions to this system of
equations map to solutions of the linearized Euler-Poincaré system, Eqs. (III)
and (14), under $\mathfrak{C}$ and satisfy Eqs. (15) and (16). Thus, these
evolution equations satisfy property (a). Conversely, the derivation of Eqs.
(22) and (23) given earlier in this section proves that these evolution
equations satisfy property (b). Therefore, Eqs. (22) and (23) are satisfactory
evolution equations for the Clebsch variables.

Finally, we turn to finding a quadratic functional preserved by the dynamics
of the dynamically accessible perturbations. Because $l$ has no explicit time-
dependence, energy is conserved, which is given by

$\displaystyle h(\mathbf{v},a)=\left<\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v},a),\mathbf{v}\right>-l(\mathbf{v},a).$ (27)

Therefore, when this energy functional is evaluated along the
$\Gamma_{t,\epsilon}$ introduced earlier in this section, we have
$\frac{\mathrm{d}}{\mathrm{d}t}h\circ\Gamma_{t,\epsilon}=0$. Then, assuming
partial derivatives in $\epsilon$ and $t$ commute, we have a primitive form of
a linearized conservation law:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\frac{\mathrm{d}^{2}}{\mathrm{d}\epsilon^{2}}\bigg{|}_{\epsilon=0}h\circ\Gamma_{t,\epsilon}=0.$
(28)

Using the fact that $\Gamma_{t,\epsilon}$ is a solution to the Euler-Poincaré
system with fixed values of the Noether invariants associated with the
particle relabeling symmetry,
$\frac{\mathrm{d}^{2}}{\mathrm{d}\epsilon^{2}}\big{|}_{0}h\circ\Gamma_{t,\epsilon}$
can be expressed in the form

$\displaystyle\frac{\mathrm{d}^{2}}{\mathrm{d}\epsilon^{2}}\bigg{|}_{0}h\circ\Gamma_{t,\epsilon}=$
(29)
$\displaystyle\left\langle\delta\left(\mathfrak{L}_{\mathbf{v}_{t,\epsilon}}\frac{\delta
l}{\delta\mathbf{v}}-\frac{\delta l}{\delta a}\diamond
a_{t,\epsilon}\right),\bm{\xi}_{t}\right\rangle-\left\langle\delta(\mathfrak{L}_{\mathbf{v}_{t,\epsilon}}{a_{t,\epsilon}}),u_{t}\right\rangle,$

where, as before $\delta=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\big{|}_{0}$.
Using $\mathfrak{C}$, this quantity can be expressed entirely in terms of the
Clebsch variables $\bm{\xi}_{t}$ and $u_{t}$. Thus, the quadratic functional

$\displaystyle\delta^{2}h(\bm{\xi},u)=$ (30)
$\displaystyle-\bigg{<}\mathfrak{L}_{\delta\mathbf{v}}\frac{\delta
l}{\delta\mathbf{v}}+\mathfrak{L}_{\mathbf{v}_{o}}\left(\frac{\delta^{2}l}{\delta\mathbf{v}^{2}}[\delta\mathbf{v}]\right)+\mathfrak{L}_{\mathbf{v}_{o}}\left(\frac{\delta^{2}l}{\delta
a\delta\mathbf{v}}[\delta a]\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-\left(\frac{\delta^{2}l}{\delta\mathbf{v}\delta
a}[\delta\mathbf{v}]+\frac{\delta^{2}l}{\delta a^{2}}[\delta a]\right)\diamond
a_{o}-\frac{\delta l}{\delta a}\diamond\delta a,\bm{\xi}\bigg{>}$
$\displaystyle-\left<\mathfrak{L}_{\delta\mathbf{v}}a_{o}+\mathfrak{L}_{\mathbf{v}_{o}}\delta
a,u\right>,$

where $(\delta\mathbf{v},\delta a)=\mathfrak{C}(\bm{\xi},u)$, is conserved by
the Clebsch variable evolution equations. Therefore, whenever this functional
is positive definite, it defines a norm (in fact, an inner product!) on the
space of Clebsch variables, $\mathfrak{X}(\mathcal{D})\times V$, that bounds
solutions to the Clebsch variable evolution equations for all time. This means
positive definiteness of $\delta^{2}h$ is sufficient to guarantee stability of
the Clebsch variable evolution equations Holm  _et al._ (1985). Likewise, if
$||\cdot||:\mathfrak{X}(\mathcal{D})\times V^{*}$ is any norm such that the
linear mapping $\mathfrak{C}:\mathfrak{X}(\mathcal{D})\times
V\rightarrow\mathfrak{X}(\mathcal{D})\times V$ is bounded, the norm of the
evolution of a dynamically accessible $(\delta\mathbf{v}_{t},\delta a_{t})$,
$||(\delta\mathbf{v}_{t},\delta
a_{t})||\leq||\mathfrak{C}||\cdot||(\bm{\xi}_{t},u_{t})||$, will be bounded
for all time when $\delta^{2}h$ is positive definite, thus implying stability
of (linear) dynamically accessible perturbations. Note that
$\delta^{2}h(\bm{\xi},u)$ _does not_ directly define a norm on the space of
dynamically accessible variations when it is positive definite because Eq.
(29) is not a quadratic functional of $(\delta\mathbf{v},\delta a)$.

##  IV An elementary derivation of the dynamically accessible stability
criterion

Now we will apply the general theory developed in the previous two sections to
relative equilibria of the ideal MHD equations. First, we will cast the
nonlinear MHD equations in the form of an Euler-Poincaré system using the
Lagrangian

$\displaystyle
L[\mathbf{x},\mathbf{v},\rho_{o}\,\mathrm{d}^{3}x,s_{o},\mathbf{B}_{o}\cdot\mathrm{d}\mathbf{S}]=$
(31)
$\displaystyle\int_{\mathcal{D}}\left(\frac{1}{2}\mathbf{v}\cdot\mathbf{v}-\epsilon\left(\mathbf{x}^{-1*}s_{o},\mathcal{J}(\mathbf{x}^{-1})\mathbf{x}^{-1*}\rho_{o}\right)\right)\mathbf{x}^{-1*}(\rho_{o}\mathrm{d}^{3}x)$
$\displaystyle-\frac{1}{2}\int_{\mathcal{D}}\mathcal{J}(\mathbf{x}^{-1})^{2}(\mathbf{x}^{-1*}\mathbf{B}_{o})\cdot(\mathbf{x}^{-1*}\mathbf{B}_{o})\,\mathrm{d}^{3}x$

where
$\mathcal{J}(\mathbf{x})\,\mathrm{d}^{3}x\equiv\mathbf{x}^{*}\,\mathrm{d}^{3}x$
defines the Jacobian of the mapping $\mathbf{x}$,
$a\equiv(\rho_{o}\,\mathrm{d}^{3}x,s_{o},\mathbf{B}_{o}\cdot\mathrm{d}\mathbf{S})$
are the advected parameters, and the specific internal energy $\epsilon$ is
defined by the first law of thermodynamics
$\mathrm{d}\epsilon=T(s,\rho)\mathrm{d}s-p(s,\rho)\mathrm{d}(1/\rho)$, with
$T$ and $p$ being the temperature and the pressure respectively. The advected
quantities correspond to the local conservation (advection) of mass, entropy,
and magnetic flux respectively.

It is straightforward to verify that this Lagrangian satisfies the symmetry
property characteristic of Euler-Poincaré fluids,
$L(\mathbf{x}\circ\mathbf{y},\mathbf{v},\mathbf{y}^{*}a)=L(\mathbf{x},\mathbf{v},a)$.
Therefore, the reduced Lagrangian can be calculated, giving

$\displaystyle
l(\mathbf{v},a)=\int_{\mathcal{D}}\left(\frac{1}{2}\mathbf{v}\cdot\mathbf{v}-\epsilon(s,\rho)\right)\rho\,\mathrm{d}^{3}x-\frac{1}{2}\int_{\mathcal{D}}\mathbf{B}\cdot\mathbf{B}\,\mathrm{d}^{3}x,$
(32)

with $a=(\rho\,\mathrm{d}^{3}x,s,\mathbf{B}\cdot\mathrm{d}\mathbf{S})$. Using
the constrained variations

$\displaystyle\delta\mathbf{v}$ $\displaystyle=$
$\displaystyle\dot{\bm{\xi}}+\mathbf{v}\cdot\nabla\bm{\xi}-\bm{\xi}\cdot\nabla\mathbf{v},$
(33) $\displaystyle\delta\rho$ $\displaystyle=$
$\displaystyle-\nabla\cdot(\rho\bm{\xi}),$ (34) $\displaystyle\delta s$
$\displaystyle=$ $\displaystyle-\bm{\xi}\cdot\nabla s,$ (35)
$\displaystyle\delta\mathbf{B}$ $\displaystyle=$
$\displaystyle\nabla\times(\bm{\xi}\times\mathbf{B}),$ (36)

which can be derived from Eq. (7), the reduced variational principle
$\delta\int l\,\mathrm{d}t=0$ then leads to the momentum equation,

$\frac{\mathrm{d}}{\mathrm{d}t}(\rho\mathbf{v})=-\nabla\cdot(\rho\mathbf{v}\mathbf{v})+(\nabla\times\mathbf{B})\times\mathbf{B}-\nabla
p,$ (37)

which corresponds to Eq. (9) with $\delta
L/\delta\mathbf{v}=\rho\mathbf{v}\otimes\mathrm{d}^{3}x$. This equation,
together with the analogue of Eq. (10),

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\rho$ $\displaystyle=$
$\displaystyle-\nabla\cdot(\rho\mathbf{v}),$ (38)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}s$ $\displaystyle=$
$\displaystyle-\mathbf{v}\cdot\nabla s,$ (39)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{B}$ $\displaystyle=$
$\displaystyle\nabla\times(\mathbf{v}\times\mathbf{B}),$ (40)

then gives the complete ideal MHD equations as an Euler-Poincaré system. Such
a formulation was first presented by Newcomb Newcomb (1962).

Using this formulation, we will now present Frieman and Rotenberg’s stability
criterion Frieman and Rotenberg (1960) for a relative equilibrium
$(\mathbf{v}_{o},a_{o})$. First, we linearize the ideal MHD equations
(37)-(40) and obtain the linearized MHD equations for perturbations
$(\delta\mathbf{v},\delta a)$. Then, by expressing $(\delta\mathbf{v},\delta
a)$ in terms of $(\bm{\xi},\dot{\bm{\xi}})$ using Eqs. (33)-(36), we have the
linearized equation of motion Frieman and Rotenberg (1960)

$\displaystyle\rho_{o}\ddot{\bm{\xi}}+2\rho_{o}\mathbf{v}_{o}\cdot\nabla\dot{\bm{\xi}}=\mathbf{F}(\bm{\xi}),$
(41)

where

$\displaystyle\mathbf{F}(\bm{\xi})=\nabla\cdot(\rho_{o}\bm{\xi}\mathbf{v}_{o}\cdot\nabla\mathbf{v}_{o}-\rho_{o}\mathbf{v}_{o}\mathbf{v}_{o}\cdot\nabla\bm{\xi})+\mathbf{F}_{0}(\bm{\xi}),$
(42)

and $\mathbf{F}_{0}$ is the linearized force operator for equilibria without
flow Freidberg (1987). Using Eq. (41) and the self-adjointness of
$\mathbf{F}$, it is simple to verify that the quantity

$\displaystyle
H_{2}(\bm{\xi},\dot{\bm{\xi}})=\int_{\mathcal{D}}\dot{\bm{\xi}}\cdot\dot{\bm{\xi}}\rho_{o}\,\mathrm{d}^{3}x-\int_{\mathcal{D}}\bm{\xi}\cdot\mathbf{F}(\bm{\xi})\,\mathrm{d}^{3}x$
(43)

is conserved Frieman and Rotenberg (1960); Hameiri (2003). A simple argument
Freidberg (1987) then shows that this implies the relative equilibrium is
stable when
$\delta^{2}W(\bm{\xi})=-\int\bm{\xi}\cdot\mathbf{F}(\bm{\xi})\,\mathrm{d}^{3}x$
is positive definite. This is exactly Frieman and Rotenberg’s sufficient
condition for linear stability.

Next, we will derive Hameiri’s sufficient condition for linear stability of a
relative equilibrium against dynamically accessible perturbations. Our method
is essentially that discussed earlier in the general context of Euler-Poincaré
fluids, but all of the manipulations we will perform can be verified easily
without recourse to the abstract theory. Consider a linear perturbation
$(\delta\mathbf{v},\delta a)$ to equilibrium $(\mathbf{v}_{o},a_{o})$ of the
dynamically accessible form Hameiri (2003); Morrison (1998)

$\displaystyle\delta\mathbf{v}$
$\displaystyle=\bm{\xi}\times(\nabla\times\mathbf{v}_{o})+{(\nabla\times\bm{\zeta})\times\mathbf{B}_{o}}/{\rho_{o}}+\alpha\nabla
s_{o}-\nabla\beta$ (44) $\displaystyle\delta\rho$
$\displaystyle=-\nabla\cdot(\rho_{o}\bm{\xi})$ (45) $\displaystyle\delta s$
$\displaystyle=-\bm{\xi}\cdot\nabla s_{o}$ (46)
$\displaystyle\delta\mathbf{B}$
$\displaystyle=\nabla\times(\bm{\xi}\times\mathbf{B}_{o})$ (47)

and regard it as an initial condition to the MHD equations linearized about
the equilibrium.

First we will show that, as this initial perturbation evolves under the
linearized MHD equations, it will remain in this form. To see this, suppose
that $(\bm{\xi},\bm{\zeta},\alpha,\beta)$ satisfied the evolution equations

$\displaystyle\dot{\bm{\xi}}+\mathbf{v}_{o}\cdot\nabla\bm{\xi}-\bm{\xi}\cdot\nabla\mathbf{v}_{o}$
$\displaystyle=\delta\mathbf{v}$ (48)
$\displaystyle\dot{\beta}+\mathbf{v}_{o}\cdot\nabla\beta$
$\displaystyle={\delta
p}/{\rho_{o}}+\bm{\xi}\cdot{(\nabla\times\mathbf{B}_{o})\times\mathbf{B}_{o}}/{\rho_{o}}$
(49) $\displaystyle\dot{\alpha}+\mathbf{v}_{o}\cdot\nabla\alpha$
$\displaystyle=\delta T+\bm{\xi}\cdot\nabla T_{o}$ (50)
$\displaystyle\dot{\bm{\zeta}}+\mathbf{v}_{o}\cdot\nabla\bm{\zeta}+\nabla{\mathbf{v}_{o}}\cdot\bm{\zeta}$
$\displaystyle=\delta\mathbf{B}+\bm{\xi}\cdot\nabla\mathbf{B}_{o}+\nabla\bm{\xi}\cdot\mathbf{B}_{o},$
(51)

where every instance of $\delta\mathbf{v}$ or $\delta a$ is to be expressed in
terms of $(\bm{\xi},\bm{\zeta},\alpha,\beta)$ using Eqs. (44)-(47). Then it is
a straightforward exercise in algebra to verify that the resulting
$(\delta\mathbf{v}_{t},\delta a_{t})$ satisfies the linearized MHD equations.
But this $(\delta\mathbf{v}_{t},\delta a_{t})$ is obviously in the form
specified by Eqs. (44)-(47), which verifies the claim because of uniqueness of
solutions to the linearized equations.

Next, we will derive a condition for stability of the solutions to the
evolution equation specified in Eqs. (48)-(51). As a first step, notice that
the variable $\bm{\xi}$ satisfies Eq. (41). Thus,
$H_{2}(\mathbf{\xi},\dot{\mathbf{\xi}})$ is constant in time. But this means
that the quadratic functional

$\displaystyle\delta^{2}h(\bm{\xi},\bm{\zeta},\alpha,\beta)=$ (52)
$\displaystyle\int_{\mathcal{D}}(\delta\mathbf{v}-\mathbf{v}_{o}\cdot\nabla\bm{\xi}+\bm{\xi}\cdot\nabla\mathbf{v}_{o})^{2}\rho_{o}\,\mathrm{d}^{3}x-\int_{\mathcal{D}}\bm{\xi}\cdot\mathbf{F}(\bm{\xi})\,\mathrm{d}^{3}x,$

where $\delta\mathbf{v}$ is expressed in terms of
$(\bm{\xi},\bm{\zeta},\alpha,\beta)$ using Eq. (44), is conserved along
solutions to the evolution equations (48)-(51). Therefore, whenever the
quadratic functional $\delta^{2}h$ is positive definite, it defines a norm
that bounds solutions to the evolution equations for
$(\bm{\xi},\bm{\zeta},\alpha,\beta)$ for all time, which implies stability.

Finally, we can give a sufficient condition for stability of the relative
equilibrium $(\mathbf{v}_{o},a_{o})$ to linear perturbations of the form of
Eqs. (44)-(47). When the functional $\delta^{2}h$ is positive definite, we
have just found that solutions to the evolution equations for
$(\bm{\xi},\bm{\zeta},\alpha,\beta)$ are bounded for all time by the norm
defined by $\delta^{2}h$. But we also know that every solution to the
linearized MHD equations with this type of initial condition can be expressed
in terms of a solution to the $(\bm{\xi},\bm{\zeta},\alpha,\beta)$ evolution
equations using Eqs. (44)-(47). Therefore, given any norm on the space of
$(\delta\mathbf{v},\delta a)$’s such that the linear mapping specified by Eqs.
(44)-(47) is bounded, $(\delta\mathbf{v}_{t},a_{t})$ will be bounded for all
time in that norm. That is to say, whenever $\delta^{2}h$ is positive
definite, linear perturbations to the relative equilibrium
$(\mathbf{v}_{o},a_{o})$ of the form specified by Eqs. (44)-(47) will be
stable. This is precisely the condition derived by HameiriHameiri (2003).
Clearly, it only directly applies to a restricted class of the initial
perturbations appearing in Frieman and Rotenberg’s stability criterion.

Note that only simple algebra, coupled with an apparently lucky guess for the
appropriate evolution equations for $(\bm{\xi},\bm{\zeta},\alpha,\beta)$, was
necessary to derive this condition. In order to understand how the evolution
equations for $(\bm{\xi},\bm{\zeta},\alpha,\beta)$ were deduced, one should
refer to the general theory for Euler-Poincaré fluids developed in the
previous section; the evolution equations follow from Eqs. (22) and (23).
Also, it is not a coincidence the quadratic functional $\delta^{2}h$ is
conserved, as it can be obtained directly Hameiri (2003) from Eq. (30), with
its conservation guaranteed by Eq. (28). The self-adjointness of
$\mathbf{F}(\bm{\xi})$ can, in turn, be proved as a result of the conservation
of $\delta^{2}h$.

Moreover, in light of the general theory, the special form of perturbations
considered in Hameiri’s stability criterion are the most general linear
perturbations that preserve the Noether invariants implied by the particle
relabeling symmetry. Indeed, Eqs. (44)-(47) follow from the expression for the
mapping $\mathfrak{C}$ defined in Eq. (25). On the other hand, from the
Poisson perspective used by Hameiri Hameiri (2003) and Morrison Morrison
(1998), this type of perturbation was found to be the most general form of
linear perturbation that preserves the Casimir invariants. Thus, we conclude
that in order to study dynamically accessible perturbations to systems that
admit both a Poisson and an Euler-Poincaré formulation, one can either use the
Poisson formulation to find perturbations that conserve the Casimirs, or the
Euler-Poincaré formulation to find perturbations that preserve the Noether
invariants implied by particle relabeling symmetry.

##  V Discussion

Eqs. (44)-(47) were previously obtained using the non-canonical Poisson
bracket for the ideal MHD and referred to as dynamically accessible variations
by MorrisonMorrison (1998). They were constructed to be the most general
linear perturbations that preserve all Casimirs of the system. From the point
of view adopted in this paper, this form of perturbation follows from Eq.
(25), and therefore gives the most general form of linear perturbation that
preserves the Noether invariants implied by the particle relabeling symmetry
Cotter and Holm (2012) and the advected parameters. We expect that this
correspondence can be proven to hold more generally in systems that admit both
a Poisson and Euler-Poincaré description using the fact the Poisson
formulation can be derived from the Euler-Poincaré formulation using the
Legendre transformHolm  _et al._ (2009).

Using Eqs. (44)-(47) as constraints, Hameiri Hameiri (2003) obtained
$\delta^{2}H$ as in Eq. (2) with $\delta\mathbf{v}$ expressed in terms of
$(\bm{\xi},{\bm{\zeta}},{\alpha},\beta)$ using Eq. (I). Hameiri then claimed
that positivity of $\delta^{2}H$ as a functional of
$(\bm{\xi},{\bm{\zeta}},{\alpha},\beta)$ was sufficient for linear stability
of the equilibrium to dynamically accessible variations. The justification of
this claim was not clearly explained. In fact, the justification can be
deduced from relevant discussions in Refs. Hameiri, 2003; Morrison, 1998.
Alternatively, the argument presented in this paper provides another proof.
More importantly, our derivation clearly shows that this criterion only
applies to perturbations that preserve the Noether invariants.

As mentioned in the introduction, Hameiri claimed that his stability criterion
is stronger than what Frieman and Rotenberg obtained Frieman and Rotenberg
(1960) via direct analysis of the linearized equations of motion (1). But in
this paper, we have shown that such a claim is incorrect because the two
stability criteria apply to different types of perturbations. In the
derivation of Frieman and Rotenberg’s criterion, only Eqs. (45)-(48) were
used, so the perturbations preserve the advected parameters but not
necessarily the Noether invariants, which is in line with the traditional
energy principle for absolute equilibria Freidberg (1987); Bernstein  _et al._
(1958). On the other hand, Hameiri’s criterion further requires the initial
$\delta\mathbf{v}$ to be of the form of Eq. (44) to preserve the Noether
invariants. However, if one only considers the dynamically accessible
perturbations, then Frieman and Rotenberg’s condition is indeed weaker than
Hameiri’s condition, i.e. the Hameiri condition could indicate stability
against dynamically accessible perturbations while the Frieman and Rotenberg
condition is inconclusive. It should also be noted that for initial
perturbations that do not conserve the advected parameters, the applicability
of both criteria, and even the traditional energy principle, becomes
questionable.

Finally, we will conclude with a summary of our results. By presenting and
comparing two existing stability criteria for ideal MHD equilibria with flow
within the same framework of the Euler-Poincaré theory, we have shown that
they actually apply to different types of perturbations. We also have
formulated a Lagrangian analogue to the method of dynamically accessible
perturbations for stability analysis that applies to general Euler-Poincaré
fluids.

###### Acknowledgements.

The authors would like to thank C. Liu and J. Squire for helpful discussions.
This work was supported by the U.S. Department of Energy under Contract No.
DE-AC02-09CH11466.

##  Appendix A Some Notation

_One-form densities_ Holm  _et al._ (1998, 2009): The dual space associated to
$\mathfrak{X}(\mathcal{D})$, $\mathfrak{X}(\mathcal{D})^{*}$, is naturally
identified with the set of one-form densities,
$\mathfrak{X}(\mathcal{D})^{*}=\Gamma(\Omega^{1}(\mathcal{D})\otimes\Omega^{3}(\mathcal{D}))$.
Thus, a typical element of $\mathfrak{X}(\mathcal{D})^{*}$ assigns to each
point $X\in\mathcal{D}$ a tensor of the form $\theta(X)\otimes\sigma_{o}$,
where $\theta$ is some one-form and $\sigma_{o}=\mathrm{d}^{3}x$ is the
(constant) standard volume form on $\mathcal{D}$. Given a vector field
$\mathbf{w}$ on $\mathcal{D}$, the Lie derivative of such a one-form density
is given by
$\mathfrak{L}_{\mathbf{w}}(\theta\otimes\sigma_{o})=(\mathfrak{L}_{\mathbf{w}}\theta)\otimes\sigma_{o}+\theta\otimes(\mathfrak{L}_{\mathbf{w}}\sigma_{o})=(\mathfrak{L}_{\mathbf{w}}\theta+\theta\nabla\cdot\mathbf{w})\otimes\sigma_{o}$.
If we denote the pairing between a vector field and a one-form density as
$\left<\theta\otimes\sigma_{o},\mathbf{w}\right>\equiv\int\theta(\mathbf{w})\sigma_{o}$,
we then see that
$\left<\theta\otimes\sigma_{o},\mathfrak{L}_{\mathbf{v}}\mathbf{w}\right>=-\left<\mathfrak{L}_{\mathbf{v}}(\theta\otimes\sigma_{o}),\mathbf{w}\right>$.
Note that, in order to derive this identity, one must utilize the fact that
vector fields are tangent to the boundary of $\mathcal{D}$.

_Functional derivatives_ : Given a functional
$l:\mathfrak{X}(\mathcal{D})\times V^{*}\rightarrow\mathbb{R}$, its partial
functional derivatives are defined by the equation

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}$
$\displaystyle l(\mathbf{v}+\epsilon\delta\mathbf{v},a+\epsilon\delta a)=$
(53) $\displaystyle\left<\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v},a),\delta\mathbf{v}\right>+\left<\delta
a,\frac{\delta l}{\delta a}(\mathbf{v},a)\right>.$

Here, $\left<\cdot,\cdot\right>$ denotes the pairing between a vector space
and its dual and $(\delta\mathbf{v},\delta
a)\in\mathfrak{X}(\mathcal{D})\times V^{*}$. We therefore see that for each
$(\mathbf{v},a)$, $\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v},a)\in\mathfrak{X}(\mathcal{D})^{*}$ is a one-
form density and $\frac{\delta l}{\delta a}(\mathbf{v},a)\in V^{**}=V$ is a
$k$-form.

We will also have occasion to use higher functional derivatives in this
article. Specifically, in Eq. (30), we employ the notations
$\frac{\delta^{2}l}{\delta
a\delta\mathbf{v}},\frac{\delta^{2}l}{\delta\mathbf{v}\delta
a},\frac{\delta^{2}}{\delta a^{2}},$ and
$\frac{\delta^{2}l}{\delta\mathbf{v}^{2}}$. The relevant definitions are

$\displaystyle\frac{\delta^{2}l}{\delta
a\delta\mathbf{v}}(\mathbf{v},a)[\delta a]$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v},a+\epsilon\delta a)$ (54)
$\displaystyle\frac{\delta^{2}l}{\delta\mathbf{v}\delta
a}(\mathbf{v},a)[\delta\mathbf{v}]$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\frac{\delta
l}{\delta a}(\mathbf{v}+\epsilon\delta\mathbf{v},a)$ (55)
$\displaystyle\frac{\delta^{2}l}{\delta a^{2}}(\mathbf{v},a)[\delta a]$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\frac{\delta
l}{\delta a}(\mathbf{v},a+\epsilon\delta a)$ (56)
$\displaystyle\frac{\delta^{2}l}{\delta\mathbf{v}^{2}}(\mathbf{v},a)[\delta\mathbf{v}]$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\bigg{|}_{0}\frac{\delta
l}{\delta\mathbf{v}}(\mathbf{v}+\epsilon\delta\mathbf{v},a).$ (57)

An unfortunate consequence of this notation is that $\frac{\delta^{2}l}{\delta
a\delta\mathbf{v}}(\mathbf{v},a)\neq\frac{\delta^{2}l}{\delta\mathbf{v}\delta
a}(\mathbf{v},a)$. That this is true is obvious from the above definitions;
$\frac{\delta^{2}l}{\delta a\delta\mathbf{v}}(\mathbf{v},a)$ is linear map
from $V^{*}$ into $\mathfrak{X}(\mathcal{D})^{*}$ while
$\frac{\delta^{2}l}{\delta\mathbf{v}\delta a}(\mathbf{v},a)$ is a linear map
from $\mathfrak{X}(\mathcal{D})$ into $V$.

_Diamond product_ Holm  _et al._ (1998, 2009): Given $(u,a)\in V\times V^{*}$,
consider the linear functional on $\mathfrak{X}(\mathcal{D})$,
$\mathbf{w}\mapsto\left<\mathfrak{L}_{\mathbf{w}}a,u\right>$. By definition,
this linear functional specifies an element, $\tilde{\theta}(u,a)$, of the
dual space $\mathfrak{X}(\mathcal{D})^{*}$ that satisfies

$\displaystyle\left<\mathfrak{L}_{\mathbf{w}}a,u\right>=\left<\tilde{\theta}(u,a),\mathbf{w}\right>.$
(58)

Note that on the left-hand side of this expression, the pairing is between $V$
and $V^{*}$, while on the right-hand side it is between
$\mathfrak{X}(\mathcal{D})$ and $\mathfrak{X}(\mathcal{D})^{*}$. Because
$\tilde{\theta}(u,a)$ must be a bilinear function of its arguments, it defines
a $\mathfrak{X}(\mathcal{D})^{*}$-valued product $(u,a)\mapsto u\diamond a$
that satisfies

$\displaystyle\left<u\diamond
a,\mathbf{w}\right>=-\left<\mathfrak{L}_{\mathbf{w}}a,u\right>.$ (59)

## References

  * Bernstein  _et al._ (1958) I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M. Kulsrud, “An energy principle for hydromagnetic stability problems,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 244, 17 (1958).
  * Freidberg (1987) J. P. Freidberg,  _Ideal magnetohydrodynamics_ (Plenum Press, New York, NY, 1987).
  * Kulsrud (1964) R. Kulsrud, “General stability theory in plasma physics,” in  _Advanced Plasma Theory_ , Vol. 1 (1964) p. 54.
  * Laval  _et al._ (1965) G. Laval, C. Mercier, and R. Pellat, “Necessity of the energy principles for magnetostatic stability,” Nuclear Fusion, 5, 156 (1965).
  * Frieman and Rotenberg (1960) E. Frieman and M. Rotenberg, “On hydromagnetic stability of stationary equilibria,” Rev. Mod. Phys., 32, 898 (1960).
  * Hameiri (2003) E. Hameiri, “Dynamically accessible perturbations and magnetohydrodynamic stability,” Physics of Plasmas, 10, 2643 (2003).
  * Almaguer  _et al._ (1988) J. A. Almaguer, E. Hameiri, J. Herrera, and D. D. Holm, “Lyapunov stability analysis of magnetohydrodynamic plasma equilibria with axisymmetric toroidal flow,” Physics of Fluids, 31, 1930 (1988).
  * Isichenko (1998) M. B. Isichenko, “Nonlinear hydrodynamic stability,” Phys. Rev. Lett., 80, 972 (1998).
  * Hameiri (1998) E. Hameiri, “Variational principles for equilibrium states with plasma flow,” Physics of Plasmas, 5, 3270 (1998).
  * Holm  _et al._ (1985) D. D. Holm, J. E. Marsden, T. Ratiu, and A. Weinstein, “Nonlinear stability of fluid and plasma equilibria,” Physics Reports, 123, 1 (1985).
  * Morrison and Greene (1980) P. J. Morrison and J. M. Greene, “Noncanonical hamiltonian density formulation of hydrodynamics and ideal magnetohydrodynamics,” Physical Review Letters, 45, 790 (1980).
  * Morrison (1998) P. J. Morrison, “Hamiltonian description of the ideal fluid,” Rev. Mod. Phys., 70, 467 (1998).
  * Morrison  _et al._ (1989) P. J. Morrison, D. Pfirsch, _et al._ , “Free-energy expressions for vlasov equilibria,” Physical Review A, 40, 3898 (1989).
  * Marsden  _et al._ (1984) J. E. Marsden, T. Raţiu, and A. Weinstein, “Semidirect products and reduction in mechanics,” Transactions of the American Mathematical Society, 281, 147 (1984).
  * Holm  _et al._ (1998) D. D. Holm, J. E. Marsden, and T. S. Ratiu, “The euler-Poincaré equations and semidirect products with applications to continuum theories,” Advances in Mathematics, 137, 1 (1998), ISSN 0001-8708.
  * Holm  _et al._ (2009) D. Holm, T. Schmah, and C. Stoica,  _Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions_, Oxford Texts in Applied and Engineering Mathematics (OUP Oxford, 2009) ISBN 9780199212903.
  * Cotter and Holm (2012) C. J. Cotter and D. D. Holm, “On noether’s theorem for the euler–poincaré equation on the diffeomorphism group with advected quantities,” Foundations of Computational Mathematics, 1 (2012).
  * Cendra  _et al._ (1998) H. Cendra, D. D. Holm, M. J. Hoyle, and J. E. Marsden, “The maxwell–vlasov equations in euler–poincaré form,” Journal of Mathematical Physics, 39, 3138 (1998).
  * Squire  _et al._ (2013) J. Squire, H. Qin, W. Tang, and C. Chandre, “The hamiltonian structure and euler-poincaré formulation of the vlasov-maxwell and gyrokinetic systems,” Physics of Plasmas, 20, 022501 (2013).
  * Note (1) $\mathbf{y}^{*}a$ is the standard notation for the pullback of the differential form $a$ by the mapping $\mathbf{y}$. See Ref. Marsden for details.
  * Note (2) In this paper we focus on stationary solution $(\mathbf{v}_{o},a_{o})$, while in fact Eqs. (20)-(23) can be valid for non-stationary solution $(\mathbf{v}_{t},a_{t})$ as well.
  * Marsden and Ratiu (1999) J. E. Marsden and T. S. Ratiu,  _Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems_ , Vol. 17 (Springer, 1999).
  * Newcomb (1962) W. A. Newcomb, “Lagrangian and hamiltonian methods in magnetohydrodynamics,” Nuclear Fusion Supplement, 2, 451 (1962).