A Three-Function Variational Principle for Stationary Nonbarotropic Magnetohydrodynamics

Abstract

The current paper is devoted to the introduction of simpler Eulerian variational principles from which all the relevant equations of nonbarotropic stationary magnetohydrodynamics can be derived for magnetic fields that lie on surfaces. A variational principle is given in terms of three independent variables for stationary nonbarotropic magnetohydrodynamic flows. This is a smaller number of variables than the eight variables that appear in the standard equations of nonbarotropic magnetohydrodynamics, which are the magnetic field, the velocity field, the specific entropy, and the density. We further investigate the case in which the flow along magnetic lines is not ideal.

1. Introduction

Variational principles for magnetohydrodynamics (MHD) were introduced previously in both Lagrangian and Eulerian approaches. Sturrock [1] studied in his book a Lagrangian variational formalism for MHD. Vladimirov and Moffatt [2] in a list of works studied a Eulerian variational approach for incompressible MHD. However, their variational approach contained three more variables in addition to the seven functions that appear in the standard equations of incompressible MHD, which are the magnetic field $\overrightarrow{B}$, velocity field $\overrightarrow{v}$, and pressure P. Kats [3] extended Moffatt’s work to compressible nonbarotropic flows, but without minimizing the number of variables and thus the computational load. Moreover, Kats demonstrated that the functions he suggested can be used to describe the motion of discontinuity surfaces [4,5]. Sakurai [6] suggested a two-function Eulerian variational principle for force-free magnetohydrodynamics and used it as a basis of a numerical scheme; his method was discussed in [1]. A method of analyzing the equations for those two variables was suggested in [7]. Yahalom and Lynden-Bell [8] combined the work of Sturrock [1] with the work of Sakurai [6] to obtain a Eulerian variational principle for barotropic MHD, which depends on only six variables. The variational derivative of the suggested action produced all the equations needed to describe barotropic MHD without additional constraints. The equations obtained were similar to the ones of Frenkel, Levich, and Stilman [9] (and also those of [10]).Yahalom [11] showed that for barotropic MHD, four functions suffice. Moreover, it was demonstrated that the discontinuities of some of those variables [12] are topological local conserved quantities.

Previous work has been concerned with barotropic MHD in which the pressure and internal energy depend only on the density and not on the specific entropy and temperature. Variational principles of nonbarotropic MHD in which the internal energy and pressure are entropy dependent and thus temperature dependent effects can be described can be found in the work of Bekenstein and Oron [13] in terms of 15 functions and A. V. Kats [3] in terms of 20 functions. We mention also the cases of incompressible flows in which density does not change in time or space and isochoric flows in which the volume of a system remains constant, which are incompressible if the specific volume in the flow is also constant; however, those cases are beyond the scope of the current work. A. V. Kats in an outstanding paper [14] (Section IV, E) demonstrated that there is a large symmetry group (which may be considered a type of gauge freedom) associated with the choice of variables, and it follows that the number of functions can be reduced.Yahalom [15,16] showed that five variables suffice to describe nonbarotropic MHD for the case in which we require a Sakurai [6] form for the magnetic field. Morrison [17] introduced a Hamiltonian formalism, but this also depends on eight canonical variables (see Table 2 in [17]).The paper of Yahalom [15] was concerned with general MHD nonstationary flows. A subsequent paper [18] was concerned with stationary flows and introduced and eight-variable stationary variational principle; here, we shall attempt to improve on this and obtain a three-variable stationary variational principle for nonbarotropic MHD. This will be done for a general case in which the magnetic field lines need not lie on entropy surfaces; for the restricted case in which the magnetic field lines lie on entropy surfaces, see [19].

Applications of this paper may arise for both linear and nonlinear stability analysis of stationary nonbarotropic MHD flows [20,21] and for designing numerical algorithms for integrating the equations of MHD [22,23,24]. Another possible application is connected to obtaining new analytic solutions in terms of the variational functions [25], as will be described below.

The plan of this work is as follows: We introduce the standard notations and equations of nonbarotropic magnetohydrodynamics for the stationary and nonstationary cases. Then, we introduce the notions of load and metage. The variational principles follow.

2. Standard Formulation of Nonbarotropic Magnetohydrodynamics

The standard set of equations solved for nonbarotropic magnetohydrodynamics are given below:

$$\frac{\partial \overrightarrow{B}}{\partial t}=\overrightarrow{\nabla }\times (\overrightarrow{v}\times \overrightarrow{B}),$$

(1)

$$\overrightarrow{\nabla }·\overrightarrow{B}=0,$$

(2)

$$\frac{\partial \rho }{\partial t}+\overrightarrow{\nabla }·(\rho \overrightarrow{v})=0,$$

(3)

$$\rho \frac{d\overrightarrow{v}}{dt}=\rho (\frac{\partial \overrightarrow{v}}{\partial t}+(\overrightarrow{v}·\overrightarrow{\nabla })\overrightarrow{v})=-\overrightarrow{\nabla }p(\rho ,s)+\frac{(\overrightarrow{\nabla }\times \overrightarrow{B})\times \overrightarrow{B}}{4\pi }.$$

(4)

$$\frac{ds}{dt}=0.$$

(5)

Here $\frac{\partial }{\partial t}$ is the partial temporal derivative, $\frac{d}{dt}$ is the material derivative, and $\overrightarrow{\nabla }$ is the standard nabla operator of vector calculus. $\overrightarrow{B}$ is the magnetic field. $\overrightarrow{v}$ is the velocity field. $\rho$ is the fluid mass density. s is the specific entropy (entropy per unit mass). Finally, $p(\rho ,s)$ is the pressure, which depends on both the density and entropy (and not just the density) through the thermodynamic equation of state (the nonbarotropic case).

The justification for these equations can be found in standard books on MHD [1]. The above is valid for a collision-dominated plasma in local thermodynamic equilibrium. Such conditions are not always fulfilled by real physical plasmas, certainly not in astrophysics or in fusion-relevant magnetic confinement studies. Yet, it is thought that the fastest macroscopic instabilities obey the above equations [12], while instabilities associated with viscous or finite conductivity terms take longer. According to a theorem by Bateman [26], every system can be described by a variational principle (including viscous plasma); the challenge is to discover a compact action functional that depends on a small amount of variational functions. This paper discusses only ideal MHD, while viscous MHD will be left for future endeavors.

Equation (1) describes magnetic field lines that are moving with the fluid (“frozen” lines). Equation (2) dictates a solenoidal field. Equation (3) dictates mass conservation. Equation (4) is the Euler equation for which pressure and the Lorentz forces apply. The term:

$$\overrightarrow{J}=\frac{\overrightarrow{\nabla }\times \overrightarrow{B}}{4\pi },$$

(6)

is the electric current density. Equation (5) is due to the lack of heat generation (null viscosity and null resistivity) in ideal nonbarotropic MHD and a lack of heat conduction; thus, the only remaining thermal process is heat convection. One needs to solve for eight variables ($\overrightarrow{v},\overrightarrow{B},\rho ,s$) eight equations (Equations (1) and (3)–(5)). Equation (2) is an initial condition on the $\overrightarrow{B}$ field and is satisfied automatically later due to Equation (1). For the stationary case, we obtain:

$$\overrightarrow{\nabla }\times (\overrightarrow{v}\times \overrightarrow{B})=0,$$

(7)

$$\overrightarrow{\nabla }·\overrightarrow{B}=0,$$

(8)

$$\overrightarrow{\nabla }·(\rho \overrightarrow{v})=0,$$

(9)

$$\rho (\overrightarrow{v}·\overrightarrow{\nabla })\overrightarrow{v}=-\overrightarrow{\nabla }p(\rho ,s)+\frac{(\overrightarrow{\nabla }\times \overrightarrow{B})\times \overrightarrow{B}}{4\pi }.$$

(10)

$$\overrightarrow{v}·\overrightarrow{\nabla }s=0.$$

(11)

3. Variational Principle of Nonbarotropic MHD

Here, we generalize the analysis of [8] for the nonbarotropic case. Let the action be:

$$\begin{array}{ccc}\hfill A& \equiv & \int \mathcal{L}{d}^{3}xdt,\hfill \\ \hfill \mathcal{L}& \equiv & {\mathcal{L}}_1+{\mathcal{L}}_2,\hfill \\ \hfill {\mathcal{L}}_1& \equiv & \rho (\frac{1}{2}{\overrightarrow{v}}^2-\epsilon (\rho ,s))+\frac{{\overrightarrow{B}}^2}{8\pi },\hfill \\ \hfill {\mathcal{L}}_2& \equiv & \nu [\frac{\partial \rho }{\partial t}+\overrightarrow{\nabla }·(\rho \overrightarrow{v})]-\rho \alpha \frac{d\chi }{dt}-\rho \beta \frac{d\eta }{dt}-\rho \sigma \frac{ds}{dt}\hfill \\ & -& \frac{\overrightarrow{B}}{4\pi }·\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta .\hfill \end{array}$$

(12)

$\epsilon$ is the internal energy per mass. We recall the thermodynamic identities:

$$\begin{array}{ccc}\hfill d\epsilon & =& Tds-Pd\frac{1}{\rho }=Tds+\frac{P}{{\rho }^2}d\rho \hfill \\ & & \frac{\partial \epsilon }{\partial s}=T,\frac{\partial \epsilon }{\partial \rho }=\frac{P}{{\rho }^2}\hfill \\ \hfill w& =& \epsilon +\frac{P}{\rho }=\epsilon +\frac{\partial \epsilon }{\partial \rho }\rho =\frac{\partial (\rho \epsilon )}{\partial \rho }\hfill \\ \hfill dw& =& d\epsilon +d(\frac{P}{\rho })=Tds+\frac{1}{\rho }dP\hfill \end{array}$$

(13)

In the above, T is the temperature and w is the specific enthalpy. $\epsilon$ is the specific internal energy. T is the temperature. w is the specific enthalpy. A special case of the equation of state is the polytropic equation of state [27]:

$$p=K{\rho }^{\gamma }$$

(14)

K and $\gamma$ may depend on the specific entropy s. Hence:

$$\frac{\partial \epsilon }{\partial \rho }=K{\rho }^{\gamma -2}\Rightarrow \epsilon =\frac{K}{\gamma -1}{\rho }^{\gamma -1}=\frac{p}{\rho (\gamma -1)}\Rightarrow \rho \epsilon =\frac{p}{\gamma -1}$$

(15)

the last identity is up to a function dependent on s. Obviously, $\nu ,\alpha ,\beta ,\sigma$ are Lagrange multipliers, which are inserted in such a way that the variational principle yields the following equations:

$$\begin{array}{ccc}& & \frac{\partial \rho }{\partial t}+\overrightarrow{\nabla }·(\rho \overrightarrow{v})=0,\hfill \\ & & \rho \frac{d\chi }{dt}=0,\hfill \\ & & \rho \frac{d\eta }{dt}=0.\hfill \\ & & \rho \frac{ds}{dt}=0.\hfill \end{array}$$

(16)

$\nu ,\alpha ,\beta ,\sigma$ may be multiple-valued. If $\rho$ does not vanish, these are just the continuity Equation (3), entropy conservation, and a declaration that Sakurai’s functions are moving with the flow. Varying with respect to $\overrightarrow{B}$, we obtain:

$$\overrightarrow{B}=\widehat{\overrightarrow{B}}\equiv \overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta .$$

(17)

Hence, $\overrightarrow{B}$ is given in Sakurai’s form and thus respects Equation (2). It can be shown that if $\overrightarrow{B}$ is in Sakurai’s form and Equation (16) is fulfilled, then it follows that also Equation (1) is respected.

We showed that the entire set of equations of nonbarotropic MHD is obtained from the above action except Euler’s equations. We now take care of that. Let us vary the above action with respect to $\overrightarrow{v}$:

$$\begin{array}{ccc}\hfill {\delta }_{\overrightarrow{v}}A& =& \int dt\{\int d^{3}xdt\rho \delta \overrightarrow{v}·[\overrightarrow{v}-\overrightarrow{\nabla }\nu -\alpha \overrightarrow{\nabla }\chi -\beta \overrightarrow{\nabla }\eta -\sigma \overrightarrow{\nabla }s]\hfill \\ & +& \oint d\overrightarrow{S}·\delta \overrightarrow{v}\rho \nu +\int d\overrightarrow{\Sigma }·\delta \overrightarrow{v}\rho [\nu ]\}.\hfill \end{array}$$

(18)

$\oint d\overrightarrow{S}·\delta \overrightarrow{v}\rho \nu$ is zero in generic cases. For astrophysical scenarios, $\rho =0$ on the boundary of the flow domain; in the case of a fluid contained in some pipe, the conditions $\delta \overrightarrow{v}·\widehat{n}=0$ are imposed ($\widehat{n}$ denotes the unit vector normal to the boundary). $\int d\overrightarrow{\Sigma }$ on the cut of $\nu$ is zero if $\nu$ is single-valued and $[\nu ]=0$, which is the case for generic topologies. For $\nu$, which is multiple-valued, only a Kutta velocity variation [28] that is parallel to the cut will remove the cut integral. An arbitrary velocity variation on the cut dictates $\rho =0$ on said surface, which contradicts the notion that a cut is arbitrary, as is the zero line of the azimuth. Later, we demonstrate that the “cut” is comoving; thus, it may become quite involved. This difficult reality may be more convenient to handle in symmetrical cases.

If the surface integrals vanish and also we have ${\delta }_{\overrightarrow{v}}A=0$ for an arbitrary velocity variation, we find that:

$$\overrightarrow{v}=\widehat{\overrightarrow{v}}\equiv \overrightarrow{\nabla }\nu +\alpha \overrightarrow{\nabla }\chi +\beta \overrightarrow{\nabla }\eta +\sigma \overrightarrow{\nabla }s.$$

(19)

This resembles the Clebsch form in nonmagnetic fluids [29,30]. Taking the variation with respect to $\rho$, we have:

$$\begin{array}{ccc}\hfill {\delta }_{\rho }A& =& \int d^{3}xdt\delta \rho [\frac{1}{2}{\overrightarrow{v}}^2-w-\frac{\partial \nu }{\partial t}-\overrightarrow{v}·\overrightarrow{\nabla }\nu ]\hfill \\ & +& \int dt\oint d\overrightarrow{S}·\overrightarrow{v}\delta \rho \nu +\int dt\int d\overrightarrow{\Sigma }·\overrightarrow{v}\delta \rho [\nu ]\hfill \\ & +& \int d^3{x\nu \delta \rho |}_{t_0}^{t_1}.\hfill \end{array}$$

(20)

$w=\frac{\partial (\epsilon \rho )}{\partial \rho }$ is the enthalpy per mass. Provided $\oint d\overrightarrow{S}·\overrightarrow{v}\delta \rho \nu$ vanishes on the boundary and $\int d\overrightarrow{\Sigma }·\overrightarrow{v}\delta \rho [\nu ]$ vanishes on the cut of $\nu$ for $\nu$, which is multiple-valued (thus, either a Kutta-type condition for the velocity in contradiction with the “cut” being an arbitrary surface or a vanishing density variation on the same) and in the beginning and the end, the following equation results:

$$\frac{d\nu }{dt}=\frac{1}{2}{\overrightarrow{v}}^2-w,$$

(21)

Since the right-hand side of the equation is single-valued as it is composed of physical variables that are not potentials, we obtain:

$$\frac{d[\nu ]}{dt}=0.$$

(22)

Thus, the cut is comoving with the flow and may become complicated. This situation may be less restrictive when the flow is symmetrical.

Finally, we vary with respect to both $\chi$ and $\eta$, leading to the following expressions:

$$\begin{array}{ccc}\hfill {\delta }_{\chi }A& =& \int d^{3}xdt\delta \chi [\frac{\partial (\rho \alpha )}{\partial t}+\overrightarrow{\nabla }·(\rho \alpha \overrightarrow{v})-\overrightarrow{\nabla }\eta ·\overrightarrow{J}]\hfill \\ & +& \int dt\oint d\overrightarrow{S}·[\frac{\overrightarrow{B}}{4\pi }\times \overrightarrow{\nabla }\eta -\overrightarrow{v}\rho \alpha ]\delta \chi \hfill \\ & +& \int dt\int d\overrightarrow{\Sigma }·[\frac{\overrightarrow{B}}{4\pi }\times \overrightarrow{\nabla }\eta -\overrightarrow{v}\rho \alpha ][\delta \chi ]\hfill \\ & -& \int d^3{x\rho \alpha \delta \chi |}_{t_0}^{t_1},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\delta }_{\eta }A& =& \int d^{3}xdt\delta \eta [\frac{\partial (\rho \beta )}{\partial t}+\overrightarrow{\nabla }·(\rho \beta \overrightarrow{v})+\overrightarrow{\nabla }\chi ·\overrightarrow{J}]\hfill \\ & +& \int dt\oint d\overrightarrow{S}·[\overrightarrow{\nabla }\chi \times \frac{\overrightarrow{B}}{4\pi }-\overrightarrow{v}\rho \beta ]\delta \eta \hfill \\ & +& \int dt\int d\overrightarrow{\Sigma }·[\overrightarrow{\nabla }\chi \times \frac{\overrightarrow{B}}{4\pi }-\overrightarrow{v}\rho \beta ][\delta \eta ]\hfill \\ & -& \int d^3{x\rho \beta \delta \eta |}_{t_0}^{t_1}.\hfill \end{array}$$

(24)

If the temporal and boundary conditions are satisfied with respect to the variations $\delta \chi$ and $\delta \eta$ on the domain boundary and on the cuts for the case that some (or all) of the relevant variables are multiple-valued, we obtain the following equations:

$$\frac{d\alpha }{dt}=\frac{\overrightarrow{\nabla }\eta ·\overrightarrow{J}}{\rho },\frac{d\beta }{dt}=-\frac{\overrightarrow{\nabla }\chi ·\overrightarrow{J}}{\rho },$$

(25)

in which Equation (3) is used. By suitable temporal conditions, we require that $\delta \eta$ and $\delta \chi$ vanish at the initial and final times. Cases making the boundary term null include a boundary located at infinity in which both $\overrightarrow{B}$ and $\rho$ are null or a boundary that is impermeable and perfectly conducting. For the integral over the “cuts” to become null, one can use $\delta \eta$ and $\delta \chi$, which are single-valued. It is shown that $\chi$ can always be single-valued; hence, taking $\delta \chi$ to be single-valued is not a restriction. In some topologies, $\eta$ is multiple-valued, and for those cases, a single-valued $\delta \eta$ suffices to make the cut term vanish.

Finally, we take a variation of the action with respect to s:

$$\begin{array}{ccc}\hfill {\delta }_{s}A& =& \int d^{3}xdt\delta s[\frac{\partial (\rho \sigma )}{\partial t}+\overrightarrow{\nabla }·(\rho \sigma \overrightarrow{v})-\rho T]\hfill \\ & +& \int dt\oint d\overrightarrow{S}·\rho \sigma \overrightarrow{v}\delta s-\int d^3{x\rho \sigma \delta s|}_{t_0}^{t_1},\hfill \end{array}$$

(26)

in which the temperature is $T=\frac{\partial \epsilon }{\partial s}$. We notice that according to Equation (19), $\sigma$ is single-valued, and hence, no cuts are needed. Taking into account the continuity Equation (3) we obtain for locations in which the density $\rho$ is not null, the result:

$$\frac{d\sigma }{dt}=T,$$

(27)

provided that ${\delta }_{s}A$ vanishes for an arbitrary $\delta s$.

4. Euler’s Equations

We shall now show that a velocity field given by Equation (19), such that the equations for $\alpha ,\beta ,\chi ,\eta ,\nu ,\sigma ,s$ satisfy the corresponding Equations (16), (21), (25), and (27) must satisfy Euler’s equations. Let us calculate the material derivative of $\overrightarrow{v}$:

$$\begin{array}{ccc}\hfill \frac{d\overrightarrow{v}}{dt}& =& \frac{d\overrightarrow{\nabla }\nu }{dt}+\frac{d\alpha }{dt}\overrightarrow{\nabla }\chi +\alpha \frac{d\overrightarrow{\nabla }\chi }{dt}+\frac{d\beta }{dt}\overrightarrow{\nabla }\eta +\beta \frac{d\overrightarrow{\nabla }\eta }{dt}\hfill \\ & +& \frac{d\sigma }{dt}\overrightarrow{\nabla }s+\sigma \frac{d\overrightarrow{\nabla }s}{dt}.\hfill \end{array}$$

(28)

It can be easily shown that:

$$\begin{array}{ccc}\hfill \frac{d\overrightarrow{\nabla }\nu }{dt}& =& \overrightarrow{\nabla }\frac{d\nu }{dt}-\overrightarrow{\nabla }{v}_k\frac{\partial \nu }{\partial x_k}=\overrightarrow{\nabla }(\frac{1}{2}{\overrightarrow{v}}^2-w)-\overrightarrow{\nabla }{v}_k\frac{\partial \nu }{\partial x_k},\hfill \\ \hfill \frac{d\overrightarrow{\nabla }\eta }{dt}& =& \overrightarrow{\nabla }\frac{d\eta }{dt}-\overrightarrow{\nabla }{v}_k\frac{\partial \eta }{\partial x_k}=-\overrightarrow{\nabla }{v}_k\frac{\partial \eta }{\partial x_k},\hfill \\ \hfill \frac{d\overrightarrow{\nabla }\chi }{dt}& =& \overrightarrow{\nabla }\frac{d\chi }{dt}-\overrightarrow{\nabla }{v}_k\frac{\partial \chi }{\partial x_k}=-\overrightarrow{\nabla }{v}_k\frac{\partial \chi }{\partial x_k},\hfill \\ \hfill \frac{d\overrightarrow{\nabla }s}{dt}& =& \overrightarrow{\nabla }\frac{ds}{dt}-\overrightarrow{\nabla }{v}_k\frac{\partial s}{\partial x_k}=-\overrightarrow{\nabla }{v}_k\frac{\partial s}{\partial x_k}.\hfill \end{array}$$

(29)

$x_k$ is a Cartesian coordinate, and a summation convention is assumed. Inserting the result from Equations (16) and (29) into Equation (28) yields:

$$\begin{array}{ccc}\hfill \frac{d\overrightarrow{v}}{dt}& =& -\overrightarrow{\nabla }{v}_k(\frac{\partial \nu }{\partial x_k}+\alpha \frac{\partial \chi }{\partial x_k}+\beta \frac{\partial \eta }{\partial x_k}+\sigma \frac{\partial s}{\partial x_k})\hfill \\ & +& \overrightarrow{\nabla }(\frac{1}{2}{\overrightarrow{v}}^2-w)+T\overrightarrow{\nabla }s\hfill \\ & +& \frac{1}{\rho }((\overrightarrow{\nabla }\eta ·\overrightarrow{J})\overrightarrow{\nabla }\chi -(\overrightarrow{\nabla }\chi ·\overrightarrow{J})\overrightarrow{\nabla }\eta )\hfill \\ & =& -\overrightarrow{\nabla }{v}_k{v}_k+\overrightarrow{\nabla }(\frac{1}{2}{\overrightarrow{v}}^2-w)+T\overrightarrow{\nabla }s\hfill \\ & +& \frac{1}{\rho }\overrightarrow{J}\times (\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )\hfill \\ & =& -\frac{\overrightarrow{\nabla }p}{\rho }+\frac{1}{\rho }\overrightarrow{J}\times \overrightarrow{B}.\hfill \end{array}$$

(30)

Equations (17) and (19) are used. This shows that the nonbarotropic MHD Euler equations can be obtained from the action of Equation (12), and thus, all the equations of nonbarotropic MHD can be derived from the action for arbitrary volume variations restricted only on the relevant boundaries and cuts. Boundary variations serve to obtain boundary and initial conditions for the equations.

5. Simplified Action

It might be argued that the above claims are misleading. A simplified action for nonbarotropic MHD is presented, and instead, six more functions $\alpha ,\beta ,\chi ,\eta ,\nu ,\sigma$ are added to $\overrightarrow{B},\overrightarrow{v},\rho ,s$. The action described in Equation (12) in a pedagogical form can be simplified. It is clear that the Lagrangian density of Equation (12) can be written as:

$$\begin{array}{ccc}\hfill \mathcal{L}& =& -\rho [\frac{\partial \nu }{\partial t}+\alpha \frac{\partial \chi }{\partial t}+\beta \frac{\partial \eta }{\partial t}+\sigma \frac{\partial s}{\partial t}+\epsilon (\rho ,s)]\hfill \\ & +& \frac{1}{2}\rho [{(\overrightarrow{v}-\widehat{\overrightarrow{v}})}^2-{(\widehat{\overrightarrow{v}})}^2]\hfill \\ & +& \frac{1}{8\pi }[{(\overrightarrow{B}-\widehat{\overrightarrow{B}})}^2-{(\widehat{\overrightarrow{B}})}^2]+\frac{\partial (\nu \rho )}{\partial t}+\overrightarrow{\nabla }·(\nu \rho \overrightarrow{v}).\hfill \end{array}$$

(31)

$\widehat{\overrightarrow{v}}$ is a shorthand notation for $\overrightarrow{\nabla }\nu +\alpha \overrightarrow{\nabla }\chi +\beta \overrightarrow{\nabla }\eta +\sigma \overrightarrow{\nabla }s$ (see Equation (19)), and $\widehat{\overrightarrow{B}}$ is a shorthand notation for $\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta$ (see Equation (17)). Thus, $\mathcal{L}$ has four contributions:

$$\begin{array}{ccc}\hfill \mathcal{L}& =& \widehat{\mathcal{L}}+{\mathcal{L}}_{\overrightarrow{v}}+{\mathcal{L}}_{\overrightarrow{B}}+{\mathcal{L}}_{boundary},\hfill \\ \hfill \widehat{\mathcal{L}}& \equiv & -\rho \left[\frac{\partial \nu }{\partial t}+\alpha \frac{\partial \chi }{\partial t}+\beta \frac{\partial \eta }{\partial t}+\sigma \frac{\partial s}{\partial t}\right.\hfill \\ & +& \left.\epsilon (\rho ,s)+\frac{1}{2}{(\overrightarrow{\nabla }\nu +\alpha \overrightarrow{\nabla }\chi +\beta \overrightarrow{\nabla }\eta +\sigma \overrightarrow{\nabla }s)}^2\right]\hfill \\ & -& \frac{1}{8\pi }{(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )}^2\hfill \\ \hfill {\mathcal{L}}_{\overrightarrow{v}}& \equiv & \frac{1}{2}\rho {(\overrightarrow{v}-\widehat{\overrightarrow{v}})}^2,\hfill \\ \hfill {\mathcal{L}}_{\overrightarrow{B}}& \equiv & \frac{1}{8\pi }{(\overrightarrow{B}-\widehat{\overrightarrow{B}})}^2,\hfill \\ \hfill {\mathcal{L}}_{boundary}& \equiv & \frac{\partial (\nu \rho )}{\partial t}+\overrightarrow{\nabla }·(\nu \rho \overrightarrow{v}).\hfill \end{array}$$

(32)

The only term containing $\overrightarrow{v}$ is (${\mathcal{L}}_{boundary}$ depends on $\overrightarrow{v}$, but being a boundary term in space and time, it does not affect the derived volume equations) ${\mathcal{L}}_{\overrightarrow{v}}$. After we equate to zero the variational derivative with respect to $\overrightarrow{v}$, we obtain Equation (19), but this will otherwise have no effect on the other variational derivatives. The only term depending on $\overrightarrow{B}$ is ${\mathcal{L}}_{\overrightarrow{B}}$, and this term will lead, after we set the variational derivative to zero, to Equation (17), but will not change other variations. ${\mathcal{L}}_{boundary}$ contains only complete partial derivatives and thus does not contribute to the equations, although it can influence the boundary conditions. Thus, Equations (16), (21), (25) and (27) are derived from:

$$\begin{array}{ll}\widehat{\mathcal{L}}[\alpha ,\beta ,\chi ,\eta ,\nu ,\rho ,\sigma ,s]& =-\rho [\frac{\partial \nu }{\partial t}+\alpha \frac{\partial \chi }{\partial t}+\beta \frac{\partial \eta }{\partial t}+\sigma \frac{\partial s}{\partial t}+\epsilon (\rho ,s)\\ & +\frac{1}{2}{(\overrightarrow{\nabla }\nu +\alpha \overrightarrow{\nabla }\chi +\beta \overrightarrow{\nabla }\eta +\sigma \overrightarrow{\nabla }s)}^2]\\ & -\frac{1}{8\pi }{(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )}^2\end{array}$$

(33)

in which $\widehat{\overrightarrow{v}}$ replaces $\overrightarrow{v}$ and $\widehat{\overrightarrow{B}}$ replaces $\overrightarrow{B}$. After integrating Equations (16), (21), (25), and (27), we are allowed to insert the potentials $\alpha ,\beta ,\chi ,\eta ,\nu ,\sigma ,s$ into Equations (17) and (19), thus deriving the vector fields $\overrightarrow{v}$ and $\overrightarrow{B}$. To summarize, we showed that the general ideal nonbarotropic MHD problem is reduced from eight equations (Equations (1) and (3)–(5)) and the additional constraint (Equation (2)) to a problem of eight first-order (in the temporal derivative) unconstrained equations. This set of equations can be derived from the Lagrangian density $\widehat{\mathcal{L}}$.

6. Stationary Nonbarotropic MHD

Stationary configurations are unique to Eulerian fluid dynamics with no counterpart in the Lagrangian description of fluid dynamics. Stationary configurations are defined by the fact that the fields $\overrightarrow{v},\overrightarrow{B},\rho ,s$ are independent of time. This does not imply that the potentials $\alpha ,\beta ,\chi ,\eta ,\nu ,\sigma$ are functions of only the spatial coordinates. Indeed, choosing the potentials in such a way will lead to erroneous results such that the stationary equations of motion can not be derived from $\widehat{\mathcal{L}}$ of Equation (32). This problem can be solved as follows. Let us choose $\alpha ,\beta ,\chi ,\nu ,\sigma$ to depend on the spatial coordinates alone. $\eta$ is chosen such that:

$$\eta =\overline{\eta }-t,$$

(34)

in which $\overline{\eta }$ depends on the spatial coordinates alone. $\widehat{\mathcal{L}}$ of Equation (32) will become:

$$\widehat{\mathcal{L}}=\rho (\beta -\epsilon (\rho ,s))-\frac{1}{2}\rho {(\overrightarrow{\nabla }\nu +\alpha \overrightarrow{\nabla }\chi +\beta \overrightarrow{\nabla }\overline{\eta }+\sigma \overrightarrow{\nabla }s)}^2-\frac{1}{8\pi }{(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\overline{\eta })}^2.$$

(35)

This can be compared with [2] (Equation 6.12) for incompressible flows in which their I is analogous to our $\beta$. Notice that $\beta$ is not a conserved quantity, but I is conserved.

Varying the Lagrangian $\widehat{L}=\int \widehat{\mathcal{L}}{d}^{3}x$ with respect to $\nu ,\alpha ,\beta ,\chi ,\eta ,\rho ,\sigma ,s$ leads to the following equations:

$$\begin{array}{ccc}& & \overrightarrow{\nabla }·(\rho \widehat{\overrightarrow{v}})=0,\hfill \\ & & \rho \widehat{\overrightarrow{v}}·\overrightarrow{\nabla }\chi =0,\hfill \\ & & \rho (\widehat{\overrightarrow{v}}·\overrightarrow{\nabla }\overline{\eta }-1)=0,\hfill \\ & & \widehat{\overrightarrow{v}}·\overrightarrow{\nabla }\alpha =\frac{\overrightarrow{\nabla }\overline{\eta }·\widehat{\overrightarrow{J}}}{\rho },\hfill \\ & & \widehat{\overrightarrow{v}}·\overrightarrow{\nabla }\beta =-\frac{\overrightarrow{\nabla }\chi ·\widehat{\overrightarrow{J}}}{\rho },\hfill \\ & & \beta =\frac{1}{2}{\widehat{\overrightarrow{v}}}^2+w,\hfill \\ & & \rho \widehat{\overrightarrow{v}}·\overrightarrow{\nabla }s=0,\hfill \\ & & \rho \widehat{\overrightarrow{v}}·\overrightarrow{\nabla }\sigma =\rho T.\hfill \end{array}$$

(36)

These lead to the stationary nonbarotropic MHD equations:

$$\overrightarrow{\nabla }\times (\widehat{\overrightarrow{v}}\times \widehat{\overrightarrow{B}})=0,$$

(37)

$$\rho (\widehat{\overrightarrow{v}}·\overrightarrow{\nabla })\widehat{\overrightarrow{v}}=-\overrightarrow{\nabla }p(\rho ,s)+\frac{(\overrightarrow{\nabla }\times \widehat{\overrightarrow{B}})\times \widehat{\overrightarrow{B}}}{4\pi }.$$

(38)

In what follows, we attempt to reduce the number of variational variables from eight to three.

7. Load and Metage

The following section follows closely a similar section in [8]. Consider a thin tube surrounding a magnetic field line as described in Figure 1. The magnetic flux contained within the tube is:

$$\Delta \Phi =\int \overrightarrow{B}·d\overrightarrow{S}$$

(39)

and the mass contained within the tube is:

$$\Delta M=\int \rho d\overrightarrow{l}·d\overrightarrow{S},$$

(40)

in which $dl$ is a short distance along the tube. Since the magnetic field lines are comoving as dictated by Equations (1) and (3), both the quantities $\Delta \Phi$ and $\Delta M$ will not change during the motion of the tube. Since the tube is thin, we may define the comoving magnetic load:

$$\lambda =\frac{\Delta M}{\Delta \Phi }=\oint \frac{\rho }{B}dl,$$

(41)

calculated along the field line. The parts of the magnetic line that go out of the flow region to parts of space in which $\rho =0$ have zero contribution. $\lambda$ is single-valued, which can be experimentally measured in principle. Since $\lambda$ is comoving, it satisfies:

$$\frac{d\lambda }{dt}=0.$$

(42)

Furthermore, surfaces of constant magnetic load are comoving and are made from magnetic field lines. Hence, the gradient to such surfaces must be normal to the surface, hence orthogonal to the field lines:

$$\overrightarrow{\nabla }\lambda ·\overrightarrow{B}=0.$$

(43)

Consider an arbitrary comoving point on the field line denoted by i, and consider another comoving point on the field line and mark it as r. The integral:

$$\mu (r)={\int }_i^r\frac{\rho }{B}dl+\mu (i),$$

(44)

is also a comoving quantity, which we denote after Lynden-Bell and Katz [31] as the magnetic metage. $\mu (i)$ is a number that can be chosen differently for each line. Thus:

$$\frac{d\mu }{dt}=0.$$

(45)

By differentiating along the field line, we derive:

$$\overrightarrow{\nabla }\mu ·\overrightarrow{B}=\rho .$$

(46)

$\mu$ will be generally multiple-valued.

We now have two comoving coordinates of flow, $\lambda ,\mu$; obviously, in a three-dimensional flow, we should also have a third. However, before defining the third function, we first change $\lambda$ to a specific function of $\lambda$. Consider the magnetic flux within a load surface $\Phi (\lambda )$. The magnetic flux is a comoving quantity and depends only on $\lambda$ of its boundary. Define:

$$\chi =\frac{\Phi (\lambda )}{2\pi }.$$

(47)

$\chi$ satisfies the equations:

$$\frac{d\chi }{dt}=0,\overrightarrow{B}·\overrightarrow{\nabla }\chi =0.$$

(48)

Let us define ${\eta }^{*}$. Since $\overrightarrow{\nabla }\mu$ is not orthogonal to $\overrightarrow{B}$, we can choose $\overrightarrow{\nabla }{\eta }^{*}$ to be orthogonal to $\overrightarrow{B}$ and not be in the direction of $\overrightarrow{\nabla }\chi$. We chose ${\eta }^{*}$ not to depend uniquely on $\chi$. Both $\overrightarrow{\nabla }{\eta }^{*}$ and $\overrightarrow{\nabla }\chi$ are orthogonal to $\overrightarrow{B}$, so it must be that:

$$\overrightarrow{B}=A\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }{\eta }^{*}.$$

(49)

However, using Equation (2), we have:

$$\overrightarrow{\nabla }·\overrightarrow{B}=\overrightarrow{\nabla }A·(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }{\eta }^{*})=0.$$

(50)

This implies that A is a function of $\chi ,{\eta }^{*}$. Now, we can define a new comoving function $\eta$ such that:

$$\eta ={\int }_0^{{\eta }^{*}}A(\chi ,{\eta }^{{}^{\prime }*})d{\eta }^{{}^{\prime }*},\frac{d\eta }{dt}=0.$$

(51)

In terms of this function, we obtain the Sakurai (Euler potentials) presentation:

$$\overrightarrow{B}=\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta .$$

(52)

We thus showed how $\chi ,\eta$ can be practically defined for a given $\overrightarrow{B},\rho$. $\eta$ is defined in a nonunique way since one can redefine $\eta$ by performing: $\eta \to \eta +f(\chi )$ in which $f(\chi )$ is any arbitrary function. $\chi ,\eta$ serve as labels of the field lines. Moreover, we obtain the following expression for the magnetic flux:

$$\Phi =\int \overrightarrow{B}·d\overrightarrow{S}=\int d\chi d\eta .$$

(53)

In the case that the surface integral is performed inside a load contour, we obtain:

$$\Phi (\lambda )={\int }_{\lambda }d\chi d\eta =\chi {\int }_{\lambda }d\eta =\left\{\begin{array}{c}\chi [\eta ]\\ \chi ({\eta }_{max}-{\eta }_{min})\end{array}\right.$$

(54)

In one case, the load surfaces are topological cylinders. $\eta$ is thus multiple-valued, and we obtain the upper value. In the second case, the surfaces of constant load are topological spheres, and $\eta$ is single-valued and has both minimal ${\eta }_{min}$ and maximal ${\eta }_{max}$ values. Hence the lower value of is obtained. In some cases, $\eta$ is identical to twice the latitude angle $\theta$. In those cases, ${\eta }_{min}=0$ (value at the “north pole”) and ${\eta }_{max}=2\pi$ (value at the “south pole”).

Comparing with Equation (47), it follows that $\eta$ can be either single-valued or multiple-valued and that its discontinuity across its cut in the multiple-valued scenario is $[\eta ]=2\pi$.

So far, the discussion has not differentiated the cases of stationary and nonstationary flows. It should be noted that even for stationary flows, one can have a nonstationary $\eta$ coordinate as the magnetic field depends only on the gradient of $\eta$ (see Equation (52)), in particular if $\eta$ is stationary, then $\eta +g(t)$, which is clearly not stationary, will produce according to Equation (52) a stationary magnetic field. In what follows, it is advantageous to use the form of $\eta$ given in Equation (34) in which $\overline{\eta }$ is stationary.

$\chi ,\eta ,\mu$ are sufficient to label any fluid element in three dimensions. However, for a nonbarotropic ideal flow, there is also another possible label, s, which is comoving due to Equation (5). The question then arises about the relation of s to the other labels. Since we need to decide regarding the preferred set of labels, it may seem that the physical labels are $\chi ,\eta ,s$, which are the surfaces on which the magnetic fields lie and the entropy. Thus, $\mu$ is a function of $\chi ,\eta ,s$. If the magnetic field lines lie on the entropy surface, $\mu$ is an independent label. The density is:

$$\rho =\frac{\partial \mu }{\partial s}\frac{\partial (\chi ,\eta ,s)}{\partial (x,y,z)}.$$

(55)

$\mu$ can be defined for each line separately according to Equation (44), and it is thus obvious that such a choice exists in which $\mu$ is uniquely a function of s. One may also think of the entropy s as functions $\chi ,\eta ,\mu$. In what follows, we shall ignore the status of s as a label and consider it as a variational function, which only attains the status of a label at the variational extremum.

8. A Simpler Variational Principle of Stationary Nonbarotropic Magnetohydrodynamics

In a previous paper [18], we showed that stationary nonbarotropic magnetohydrodynamics can be described in terms of eight first-order differential equations and by an action principle from which those equations can be derived. Below, we show that one can do better for the case in which the magnetic field lines lie on an entropy surface; in this case three functions suffice to describe stationary nonbarotropic magnetohydrodynamics.

Consider Equation (48); for a stationary flow, it takes the form:

$$\overrightarrow{v}·\overrightarrow{\nabla }\chi =0.$$

(56)

Hence, $\overrightarrow{v}$ can take the form:

$$\overrightarrow{v}=\frac{\overrightarrow{\nabla }\chi \times \overrightarrow{K}}{\rho }.$$

(57)

However, the velocity field must satisfy the stationary mass conservation Equation (3):

$$\overrightarrow{\nabla }·(\rho \overrightarrow{v})=0.$$

(58)

A sufficient condition (although not necessary) for $\overrightarrow{v}$ to be a solution of Equation (58) is that $\overrightarrow{K}$ is of the form $\overrightarrow{K}=\overrightarrow{\nabla }N$; here, N is arbitrary. Thus, $\overrightarrow{v}$ can be written as:

$$\overrightarrow{v}=\frac{\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }N}{\rho }.$$

(59)

Let us now calculate $\overrightarrow{v}\times \overrightarrow{B}$ in which $\overrightarrow{B}$ is given by Sakurai’s presentation Equation (52):

$$\begin{array}{ccc}\hfill \overrightarrow{v}\times \overrightarrow{B}& =& (\frac{\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }N}{\rho })\times (\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )\hfill \\ & =& \frac{1}{\rho }\overrightarrow{\nabla }\chi (\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }N)·\overrightarrow{\nabla }\eta .\hfill \end{array}$$

(60)

N can be at most a function of the three coordinates $\chi ,\mu ,\overline{\eta }$; hence:

$$\overrightarrow{\nabla }N=\frac{\partial N}{\partial \chi }\overrightarrow{\nabla }\chi +\frac{\partial N}{\partial \mu }\overrightarrow{\nabla }\mu +\frac{\partial N}{\partial \overline{\eta }}\overrightarrow{\nabla }\overline{\eta }.$$

(61)

Inserting Equation (61) into Equation (60) yields:

$$\overrightarrow{v}\times \overrightarrow{B}=\frac{1}{\rho }\overrightarrow{\nabla }\chi \frac{\partial N}{\partial \mu }(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\mu )·\overrightarrow{\nabla }\overline{\eta }.$$

(62)

Rearranging terms and using Equation (52), we can write:

$$\overrightarrow{v}\times \overrightarrow{B}=-\frac{1}{\rho }\overrightarrow{\nabla }\chi \frac{\partial N}{\partial \mu }(\overrightarrow{\nabla }\mu ·\overrightarrow{B}).$$

(63)

However, using Equation (46), this simplifies to the form:

$$\overrightarrow{v}\times \overrightarrow{B}=-\overrightarrow{\nabla }\chi \frac{\partial N}{\partial \mu }.$$

(64)

Inserting Equation (64) into Equation (7) leads to the equation:

$$\overrightarrow{\nabla }(\frac{\partial N}{\partial \mu })\times \overrightarrow{\nabla }\chi =0.$$

(65)

However, since N is at most a function of $\chi ,\mu ,\overline{\eta }$, it follows that $\frac{\partial N}{\partial \mu }$ is some function of $\chi$:

$$\frac{\partial N}{\partial \mu }=-F(\chi ).$$

(66)

This can be easily integrated to yield:

$$N=-\mu F(\chi )+G(\chi ,\overline{\eta }).$$

(67)

Inserting N into Equation (59) results in:

$$\overrightarrow{v}=\frac{\overrightarrow{\nabla }\chi \times (-F(\chi )\overrightarrow{\nabla }\mu +\frac{\partial G}{\partial \overline{\eta }}\overrightarrow{\nabla }\overline{\eta })}{\rho }.$$

(68)

We shall replace $\chi ,\overline{\eta }$ with ${\chi }^{\prime },{\overline{\eta }}^{\prime }$ defined as:

$${\chi }^{\prime }=\int F(\chi )d\chi ,{\overline{\eta }}^{\prime }=\frac{\overline{\eta }}{F(\chi )}.$$

(69)

This does not affect Equation (52) as:

$$\overrightarrow{B}=\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta =\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\overline{\eta }=\overrightarrow{\nabla }{\chi }^{\prime }\times \overrightarrow{\nabla }{\overline{\eta }}^{\prime }.$$

(70)

However, the velocity takes a simpler form:

$$\overrightarrow{v}=\frac{\overrightarrow{\nabla }{\chi }^{\prime }\times \overrightarrow{\nabla }(-\mu +G^{\prime }({\chi }^{\prime },{\overline{\eta }}^{\prime }))}{\rho },$$

(71)

$G^{\prime }=\frac{G}{F}$. We recall that $\mu$ is defined in Equation (44) up to a constant, which may vary between field lines. As the lines are labeled by their ${\chi }^{\prime },{\overline{\eta }}^{\prime }$ values, we can add a function of ${\chi }^{\prime },{\overline{\eta }}^{\prime }$ to $\mu$ without affecting its desired properties. We thus define ${\mu }^{\prime }$ in the form:

$${\mu }^{\prime }=\mu -G^{\prime }({\chi }^{\prime },{\overline{\eta }}^{\prime }).$$

(72)

${\mu }^{\prime }$ can be nonsingle-valued. Inserting Equation (72) into Equation (71) leads to:

$$\overrightarrow{v}=\frac{\overrightarrow{\nabla }{\mu }^{\prime }\times \overrightarrow{\nabla }{\chi }^{\prime }}{\rho }.$$

(73)

The primes on $\chi ,\mu ,\overline{\eta }$ will be ignored from now on. This is analogous to [2] (Equation 7.11) for incompressible flows; our $\mu$ and $\chi$ resemble their A and $\Psi$. $\overrightarrow{v}$ satisfies the following:

$$\overrightarrow{v}·\overrightarrow{\nabla }\mu =0,\overrightarrow{v}·\overrightarrow{\nabla }\chi =0,\overrightarrow{v}·\overrightarrow{\nabla }\overline{\eta }=1,$$

(74)

where we used Equation (45) and (52). $\mu ,\chi$ are both comoving and stationary. $\overline{\eta }$ satisfies Equation (34). If:

$$basis=(\overrightarrow{\nabla }\chi ,\overrightarrow{\nabla }\overline{\eta },\overrightarrow{\nabla }\mu ),$$

(75)

is a local vector basis at any point in space, then their exists a dual basis:

$$dualbasis=\frac{1}{\rho }(\overrightarrow{\nabla }\overline{\eta }\times \overrightarrow{\nabla }\mu ,\overrightarrow{\nabla }\mu \times \overrightarrow{\nabla }\chi ,\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\overline{\eta })=(\frac{\overrightarrow{\nabla }\overline{\eta }\times \overrightarrow{\nabla }\mu }{\rho },\overrightarrow{v},\frac{\overrightarrow{B}}{\rho }),$$

(76)

such that:

$$basi{s}_i·dualbasi{s}_j={\delta }_{ij},i,j\in [1,2,3],$$

(77)

${\delta }_{ij}$ is Kronecker’s delta. The surfaces $\chi ,\mu ,\overline{\eta }$ generate a local vector basis for space, while the physical fields of interest $\overrightarrow{v},\overrightarrow{B}$ are contained in the dual basis. We can now construct a vector product of $\overrightarrow{v}$ and $\overrightarrow{B}$, and taking into account Equations (52) and (73), we arrive at the equation:

$$\overrightarrow{v}\times \overrightarrow{B}=\overrightarrow{\nabla }\chi ;$$

(78)

thus, $\overrightarrow{v}$ and $\overrightarrow{B}$ lie on $\chi$ surfaces and are a vector basis for this two-dimensional surface. This can be compared with [2] (Equation 5.6) for incompressible flows; the J appearing in that paper is analogous to $\chi$.

9. A Three-Function Variational Principle for Stationary MHD

Previously, we demonstrated that, provided that the $\overrightarrow{v}$ is given by Equation (73) and $\overrightarrow{B}$ is given by Equation (52), then Equations (7)–(9) are satisfied. To complete the set of equations we show how the Euler Equation (4) can be derived from the action:

$$\begin{array}{ccc}\hfill A& \equiv & \int \mathcal{L}{d}^{3}xdt,\hfill \\ \hfill \mathcal{L}& \equiv & \rho (\frac{1}{2}{\overrightarrow{v}}^2-\epsilon (\rho ,s))-\frac{{\overrightarrow{B}}^2}{8\pi },\hfill \end{array}$$

(79)

in which both $\overrightarrow{v}$ and $\overrightarrow{B}$ are given by Equations (52) and (73), respectively, and the density $\rho$ is given by Equation (45):

$$\rho =\overrightarrow{\nabla }\mu ·\overrightarrow{B}=\overrightarrow{\nabla }\mu ·(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )=\frac{\partial (\chi ,\eta ,\mu )}{\partial (x,y,z)}.$$

(80)

The Lagrangian density of Equation (79) takes the more explicit form:

$$\mathcal{L}[\chi ,\eta ,\mu ]=\rho \left(\frac{1}{2}{(\frac{\overrightarrow{\nabla }\mu \times \overrightarrow{\nabla }\chi }{\rho })}^2-\epsilon (\rho ,s(\chi ,\eta ,\mu ))\right)-\frac{{(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )}^2}{8\pi }$$

(81)

and can be seen explicitly to depend on only three functions. We underline that the magnetic field lines lie on entropy surfaces. s must be a function of $\chi ,\eta$ only and does not depend on $\mu$. Let us make arbitrary small variations $\delta {\alpha }_i=(\delta \chi ,\delta \eta ,\delta \mu )$ of the functions ${\alpha }_i=(\chi ,\eta ,\mu )$. Let us define a $\Delta$ variation that does not modify the ${\alpha }_i$’s, such that:

$$\Delta {\alpha }_i=\delta {\alpha }_i+(\overrightarrow{\xi }·\overrightarrow{\nabla }){\alpha }_i=0,$$

(82)

in which $\overrightarrow{\xi }$ is the Lagrangian displacement; thus:

$$\delta {\alpha }_i=-\overrightarrow{\nabla }{\alpha }_i·\overrightarrow{\xi }.$$

(83)

This leads to the equation:

$$\overrightarrow{\xi }\equiv -\frac{\partial \overrightarrow{r}}{\partial {\alpha }_i}\delta {\alpha }_i.$$

(84)

Making a variation of $\rho$ given in Equation (80) with respect to ${\alpha }_i$ yields:

$$\delta \rho =-\overrightarrow{\nabla }·(\rho \overrightarrow{\xi }).$$

(85)

Making a variation of s results in:

$$\delta s=\frac{\partial s}{\partial {\alpha }_i}\delta {\alpha }_i=-\frac{\partial s}{\partial {\alpha }_i}\overrightarrow{\nabla }{\alpha }_i·\overrightarrow{\xi }=-\overrightarrow{\nabla }s·\overrightarrow{\xi }.$$

(86)

Furthermore, taking the variation of $\overrightarrow{B}$ given by Sakurai’s Equation (52) with respect to ${\alpha }_i$ yields:

$$\delta \overrightarrow{B}=\overrightarrow{\nabla }\times (\overrightarrow{\xi }\times \overrightarrow{B}).$$

(87)

It remains to calculate $\delta \overrightarrow{v}$ by varying Equation (73); this yields:

$$\delta \overrightarrow{v}=-\frac{\delta \rho }{\rho }\overrightarrow{v}+\frac{1}{\rho }\overrightarrow{\nabla }\times (\rho \overrightarrow{\xi }\times \overrightarrow{v}).$$

(88)

Varying the action results in:

$$\begin{array}{ccc}\hfill \delta A& =& \int \delta \mathcal{L}{d}^{3}xdt,\hfill \\ \hfill \delta \mathcal{L}& =& \delta \rho (\frac{1}{2}{\overrightarrow{v}}^2-w(\rho ,s))-\rho T\delta s+\rho \overrightarrow{v}·\delta \overrightarrow{v}-\frac{\overrightarrow{B}·\delta \overrightarrow{B}}{4\pi },\hfill \end{array}$$

(89)

Inserting Equations (85), (87), and (88) into Equation (89) yields:

$$\begin{array}{ccc}\hfill \delta \mathcal{L}& =& \overrightarrow{v}·\overrightarrow{\nabla }\times (\rho \overrightarrow{\xi }\times \overrightarrow{v})-\frac{\overrightarrow{B}·\overrightarrow{\nabla }\times (\overrightarrow{\xi }\times \overrightarrow{B})}{4\pi }-\delta \rho (\frac{1}{2}{\overrightarrow{v}}^2+w)+\rho T\overrightarrow{\nabla }s·\overrightarrow{\xi }\hfill \\ & =& \overrightarrow{v}·\overrightarrow{\nabla }\times (\rho \overrightarrow{\xi }\times \overrightarrow{v})-\frac{\overrightarrow{B}·\overrightarrow{\nabla }\times (\overrightarrow{\xi }\times \overrightarrow{B})}{4\pi }+\overrightarrow{\nabla }·(\rho \overrightarrow{\xi })(\frac{1}{2}{\overrightarrow{v}}^2+w)\hfill \\ & +& \rho T\overrightarrow{\nabla }s·\overrightarrow{\xi }.\hfill \end{array}$$

(90)

Using the well-known vector identity:

$$\overrightarrow{A}·\overrightarrow{\nabla }\times (\overrightarrow{C}\times \overrightarrow{A})=\overrightarrow{\nabla }·((\overrightarrow{C}\times \overrightarrow{A})\times \overrightarrow{A})+(\overrightarrow{C}\times \overrightarrow{A})·\overrightarrow{\nabla }\times \overrightarrow{A}$$

(91)

and the theorem of Gauss, we can write now Equation (89) in the form:

$$\begin{array}{ccc}\hfill \delta A& =& \int dt\{\oint d\overrightarrow{S}·[\rho (\overrightarrow{\xi }\times \overrightarrow{v})\times \overrightarrow{v}-\frac{(\overrightarrow{\xi }\times \overrightarrow{B})\times \overrightarrow{B}}{4\pi }+(\frac{1}{2}{\overrightarrow{v}}^2+w)\rho \overrightarrow{\xi }]\hfill \\ & +& \int d^{3}x\overrightarrow{\xi }·[\rho \overrightarrow{v}\times \overrightarrow{\omega }+\overrightarrow{J}\times \overrightarrow{B}-\rho \overrightarrow{\nabla }(\frac{1}{2}{\overrightarrow{v}}^2+w)+\rho T\overrightarrow{\nabla }s]\}.\hfill \end{array}$$

(92)

Time integration is of no consequence in the above. Notice that we used Equation (6) and the vorticity $\overrightarrow{\omega }=\overrightarrow{\nabla }\times \overrightarrow{v}$. If $\delta A=0$ for $\overrightarrow{\xi }$, nullifying the boundary term, but otherwise arbitrary, then:

$$\rho \overrightarrow{v}\times \overrightarrow{\omega }+\overrightarrow{J}\times \overrightarrow{B}-\rho \overrightarrow{\nabla }(\frac{1}{2}{\overrightarrow{v}}^2+w)+\rho T\overrightarrow{\nabla }s=0.$$

(93)

Using the well-known vector identity:

$$\frac{1}{2}\overrightarrow{\nabla }({\overrightarrow{v}}^2)=(\overrightarrow{v}·\overrightarrow{\nabla })\overrightarrow{v}+\overrightarrow{v}\times (\overrightarrow{\nabla }\times \overrightarrow{v})$$

(94)

and rearranging terms, we recover the stationary Euler equation:

$$\rho (\overrightarrow{v}·\overrightarrow{\nabla })\overrightarrow{v}=-\overrightarrow{\nabla }p+\overrightarrow{J}\times \overrightarrow{B}.$$

(95)

10. An Application: Helical Stratified Magnetic Field

10.1. The Magnetic Field and Related Labels

Consider [32] a magnetohydrodynamic flow of uniform density $\rho$. Furthermore, assume that the flow contains a helical stratified magnetic field:

$$\overrightarrow{B}=\left\{\begin{array}{cc}2B_{\perp }(1-\frac{R}{a})\widehat{\varphi }+B_{z0}\widehat{z}\hfill & R<a\hfill \\ 0\hfill & R>a\hfill \end{array}\right.$$

(96)

$R,\varphi ,z$ are cylindrical coordinates. $\widehat{R},\widehat{\varphi },\widehat{z}$ are unit vectors. $B_{z0},B_{\perp }$ are constants. The field is confined to a cylinder of radius a and does not depend on z. We assume that $z=0$ and $z=L$ are identified such that a topological torus is realized. The only field lines that will be closed are such that:

$$\frac{n}{m}=\frac{B_{\perp }}{\pi R{B}_{z0}}\left(1-\frac{R}{a}\right)L,n,m\mathrm{integers}$$

(97)

Lines not satisfying this condition are surface-filling. The field lines lie on cylindrical surfaces, so one can calculate the flux through a circular surface lying on the plane z and bounded by R. The flux is:

$$\Phi =\int \overrightarrow{B}·d\overrightarrow{S}=\pi R^2{B}_{z0}$$

(98)

The $\chi$ function can now be calculated according to Equation (47) to yield the value:

$$\chi =\frac{1}{2}{B}_{z0}{R}^2.$$

(99)

Solving Equation (17) for $\eta$, we obtain the following nonunique solution:

$$\eta =\varphi -\frac{2B_{\perp }}{B_{z0}}(1-\frac{R}{a})\frac{z}{R}.$$

(100)

Finally, we solve Equation (46) for $\mu$; here, we suggest the following simple and nonunique solution:

$$\mu =\frac{\rho }{B_{z0}}z.$$

(101)

Thus, $\mu$ surfaces are just z planes. Notice that since we have identified the planes $z=0$ and $z=L$, $\mu$ is nonsingle-valued. The same can be said to be on $\eta$, which is doubly nonsingle-valued in both the z and $\varphi$ directions.

10.2. The Velocity Field

A stationary velocity field $\overrightarrow{v}$ must satisfy Equations (7) and (9). Such a velocity field can be constructed using the labels $\mu$ and $\chi$ (see Equation (73)):

$$\overrightarrow{v}=k\frac{\overrightarrow{\nabla }\mu \times \overrightarrow{\nabla }\chi }{\rho }.$$

(102)

k is a dimensional constant that we chose such that $k=\frac{v_0}{a}$. Plugging in $\mu$ from Equation (101) and $\chi$ from Equation (99), we arrive at the simple expression:

$$\overrightarrow{v}=v_0\frac{R}{a}\widehat{\varphi },$$

(103)

This expression can be shown to solve Equations (1) and (3) by direct substitution. Equation (10) can be solved by the pressure function:

$$p(R)=\rho \left[\frac{B_{\perp }^2}{\pi }\left(3\frac{R}{a}-\frac{R^2}{a^2}-ln(\frac{R}{a})-2\right)+\frac{1}{2}{v}_0^2\left(\frac{R^2}{a^2}-1\right)\right],p(a)=0.$$

(104)

11. The Three-Function Action Principle for a Static Configuration

The static configuration is a stationary flow such that $\overrightarrow{v}=0$. In this case, the mass conservation Equation (7) and magnetic field Equation (9) are satisfied trivially. To complete the set of equations, we show how the static Euler Equation (4) can be derived from the action:

$$\begin{array}{ccc}\hfill A& \equiv & -\int \mathcal{L}{d}^{3}xdt,\hfill \\ \hfill \mathcal{L}& \equiv & \rho \epsilon (\rho ,s)+\frac{{\overrightarrow{B}}^2}{8\pi },\hfill \end{array}$$

(105)

in which $\overrightarrow{B}$ is given by Equation (52) and the density $\rho$ is given by Equation (80). The Lagrangian density of Equation (105) can be put in the more explicit form:

$$\mathcal{L}[\chi ,\eta ,\mu ]=\rho \epsilon (\rho ,s(\chi ,\eta ,\mu ))+\frac{{(\overrightarrow{\nabla }\chi \times \overrightarrow{\nabla }\eta )}^2}{8\pi }$$

(106)

Varying the action results in:

$$\begin{array}{ccc}\hfill \delta A& =& -\int \delta \mathcal{L}{d}^{3}xdt,\hfill \\ \hfill \delta \mathcal{L}& =& \delta \rho (w(\rho ,s))+\rho T\delta s+\frac{\overrightarrow{B}·\delta \overrightarrow{B}}{4\pi },\hfill \end{array}$$

(107)

Inserting Equations (85) and (87) into Equation (107) yields:

$$\begin{array}{ccc}\hfill \delta \mathcal{L}& =& \frac{\overrightarrow{B}·\overrightarrow{\nabla }\times (\overrightarrow{\xi }\times \overrightarrow{B})}{4\pi }+\delta \rho w+\rho T\overrightarrow{\nabla }s·\overrightarrow{\xi }\hfill \\ & =& \frac{\overrightarrow{B}·\overrightarrow{\nabla }\times (\overrightarrow{\xi }\times \overrightarrow{B})}{4\pi }-\overrightarrow{\nabla }·(\rho \overrightarrow{\xi })w+\rho T\overrightarrow{\nabla }s·\overrightarrow{\xi }.\hfill \end{array}$$

(108)

Using the well-known vector identity (Equations (91)) and the theorem of Gauss, we can write now Equation (107) in the form:

$$\begin{array}{ccc}\hfill \delta A& =& \int dt\{\oint d\overrightarrow{S}·[-\frac{(\overrightarrow{\xi }\times \overrightarrow{B})\times \overrightarrow{B}}{4\pi }+w\rho \overrightarrow{\xi }]\hfill \\ & +& \int d^{3}x\overrightarrow{\xi }·[\overrightarrow{J}\times \overrightarrow{B}-\rho \overrightarrow{\nabla }w+\rho T\overrightarrow{\nabla }s]\}.\hfill \end{array}$$

(109)

The time integration is of no consequence in the above expression. We used Equation (6). If $\delta A=0$ for an arbitrary $\overrightarrow{\xi }$, but such that the boundary term is null, then we obtain:

$$\overrightarrow{J}\times \overrightarrow{B}-\rho \overrightarrow{\nabla }w+\rho T\overrightarrow{\nabla }s=0.$$

(110)

and rearranging terms, we recover the stationary Euler equation:

$$\overrightarrow{J}\times \overrightarrow{B}-\overrightarrow{\nabla }p=0.$$

(111)

12. Transport Phenomena

In many plasmas including static configurations, heat is transferred preferably along magnetic field lines:

$${\overrightarrow{J}}_H=-\widehat{k}\overrightarrow{\nabla }T.$$

(112)

in which $\widehat{k}$ is a tensor of heat conductivity. This tensor is usually larger in the magnetic field direction and thus can be written as:

$$\widehat{k}=k_{\perp }(I-\widehat{b}\otimes \widehat{b})+k_{\parallel }\widehat{b}\otimes \widehat{b}.$$

(113)

In the above, $\widehat{b}\equiv \frac{\overrightarrow{B}}{B}$ is a unit vector in the magnetic field direction, ⊗ is the tensor product, I is the unit matrix, $k_{\perp }$ is the heat conductivity in the directions perpendicular to the magnetic field lines, and $k_{\parallel }$ is the larger heat conductivity in the direction parallel to the magnetic field lines. The equation for a stationary heat flux configuration is:

$$\overrightarrow{\nabla }·{\overrightarrow{J}}_H=0\Rightarrow \overrightarrow{\nabla }·(\widehat{k}\overrightarrow{\nabla }T)=0.$$

(114)

This equation can be derived from the heat Lagrangian and Lagrangian density:

$$L_H\equiv \int {\mathcal{L}}_H{d}^{3}x,{\mathcal{L}}_H\equiv \frac{1}{2}{(\overrightarrow{\nabla }T)}^t\widehat{k}\overrightarrow{\nabla }T=\frac{1}{2}{\partial }_{i}T{\widehat{k}}_{ij}{\partial }_{j}T.$$

(115)

In the above, ${(\overrightarrow{\nabla }T)}^t$ is the transpose of $\overrightarrow{\nabla }T$, ${\partial }_i\equiv \frac{\partial }{\partial x_i}$, and the Einstein summation convention is assumed. Taking the variation with respect to the temperature T yields:

$$\delta {\mathcal{L}}_H={(\overrightarrow{\nabla }T)}^t\widehat{k}\overrightarrow{\nabla }\delta T;$$

(116)

hence:

$$\delta L_H=\int d^{3}x\left[\overrightarrow{\nabla }·(\widehat{k}\overrightarrow{\nabla }T\delta T)-\delta T\overrightarrow{\nabla }·(\widehat{k}\overrightarrow{\nabla }T)\right],$$

(117)

and using the theorem of Gauss:

$$\delta L_H=\int d\overrightarrow{S}·(\widehat{k}\overrightarrow{\nabla }T)\delta T-\int d^{3}x\delta T\overrightarrow{\nabla }·(\widehat{k}\overrightarrow{\nabla }T);$$

(118)

thus for appropriate boundary conditions, we derive Equation (114). We notice that heat conduction is not taken into account in ideal MHD, which only assumes the convection of heat. However, provided that conduction is seen as a secondary process with respect to convection, we may obtain using the ideal variational principle a stationary or static magnetic field configuration using the appropriate variational expression given in previous sections. Then using the known magnetic field configuration, we derive the appropriate heat flux transport using $L_H$.

13. Conclusions

It is shown that stationary nonbarotropic magnetohydrodynamics can be derived from a variational principle of three functions. We showed this for both the stationary and static cases.

Possible applications include stability studies of stationary MHD and the development of numerical methods for integrating MHD equations. It may be possible to incorporate the current formalism in existing codes instead of developing a new code from scratch. Possible codes were described in [33,34,35,36,37]. We anticipate applications of this work to linear and nonlinear stability studies of known MHD configurations [20,21,38]. To achieve this, we may need to add additional constants of motion constraints to the action [39,40,41,42,43]. As for designing numerical methods for integratingthe equations of MHD, one may follow the approach described in [22,23,24,28].

Another application of the variational variables is obtaining analytic solutions for the MHD equations. Although the equations are very difficult to solve, both being partial differential equations and nonlinear, possible solutions can be found using variational functions. An example for this approach is the self-gravitating torus described in [25] and also in Section 10.

One can use continuous symmetries that appear in the action to derive through the Noether theorem new conservation laws. An example of such a derivation can be found in [32,44].

Topological invariants have always been useful, and there are such invariants in MHD flows. For example, the two helicities (magnetic helicity and cross-helicity) have long been useful in research into the problem of hydrogen fusion and in various astrophysical scenarios. In previous studied [8,12,45], relations between the helicities and symmetries of MHD were made. Furthermore, the functions of the current action are helpful for identifying and characterizing new topological invariants in MHD [32,46,47].

Although ideal MHD does not describe real plasmas fully, we showed here how processes such as heat conduction can be also described using variational analysis provided that the magnetic field configuration is given approximately by ideal variational analysis.

To conclude, we underline the limitations of the current work. First, MHD is a fluid theory and as such misses some of the processes that can only be addressed by a detailed kinetic model. Second, ideal MHD considered here neglects important nonideal processes such as resistive heating, heat conduction, and viscous effects. Moreover, we assume magnetic fields that lie on surfaces (that is Equation (17)), while for some MHD configurations, magnetic field lines are volume filling; for such cases, the current approach is not applicable.