The First-Order Euler-Lagrange equations and some of their uses

In many nonlinear field theories, relevant solutions may be found by reducing the order of the original Euler-Lagrange equations, e.g., to first order equations (Bogomolnyi equations, self-duality equations, etc.). Here we generalise, further develop and apply one particular method for the order reduction of nonlinear field equations which, despite its systematic and versatile character, is not widely known.


I. INTRODUCTION
Nonlinear field theories are ubiquitous in the description of physical systems from particle physics [1] - [4] to condensed matter systems [5] - [7] and cosmology [8], where any genuine interaction is generally related to the nonlinearity of the underlying field theory. In these theories, one powerful strategy to obtain solutions of physical importance is to reduce the order of the original field equations (the Euler-Lagrange (EL) equations) of the system. The resulting equations of lower order -Bogomolnyi equations, self-duality equations, Bäcklund transformations, etc. -are easier to solve and allow to obtain a large number of relevant solutions with particular characteristics, like solitons, nonlinear waves, vortices, monopoles, instantons, etc. There exist several known methods to achieve this reduction of order, where the best-known one is probably the Bogomolnyi trick [9], [10], [11] of completing a square. To consider an example, let us assume that we have the energy functional of a field theory which for static fields may be expressed as a sum of two terms, E = d d x(A 2 + B 2 ) (typically, A depends on first derivatives, whereas B only depends on the fields). This may trivially be rewritten as (I.1) If, in addition, Q is a homotopy invariant (i.e., AB is locally a total derivative), then it does not contribute to the EL equations, and its value only depends on the boundary conditions imposed on the fields. As a consequence, E andĒ lead to the same EL equations. Further,Ē is non-negative, so E obeys the inequality E ≥ |Q| (Bogomolnyi bound) which is saturated by solutions to the reduced-order (usually, first-order) equation A = ±B (Bogomolnyi equation or BPS equation).
Recently it has been observed [12] that it can be useful to partly invert the logic of this construction. That is to say, let us assume that we have two functionals (functions of the fields, their derivatives, and possibly also of the coordinates x µ ) A, B which are in some sense "duals" of each other, and which are such that the product AB is locally a total derivative (the integral Q = 2 d d xAB is a homotopy invariant). This automatically implies that the "energy functional" E = d d x(A 2 + B 2 ) is a BPS action, and the "self-duality equations" (BPS equations) A = ±B provide global minima of this action. This construction is useful, because it immediately allows for some simple generalisations (to give just one example, A g = gA and B g = g −1 B have the same homotopy invariant Q and, therefore, lead to the new BPS action E g = d d x(g 2 A 2 + g −2 B 2 ) and BPS equations gA = ±g −1 B; here, g can be a rather arbitrary function of fields and coordinates).
We remark that in this paper we are mainly interested in the (local) order-reduced field equations and not so much in global considerations. We shall, therefore, use the notions of "homotopy invariant" and of "total derivative" interchangably.
The Bogomolnyi trick is very simple in simple cases (e.g. one field in one dimension), but it is not completely obvious how to generalise it to more fields and higher dimensions. More fields require, in general, to complete more squares, where frequently it is not obvious which terms should be paired into squares, so applying the method requires some guesswork. Further, the "mixed" (AB type) terms still have to add up to a homotopy invariant, which is not obvious, either. In other words, the Bogomolnyi trick does not provide a criterion as to whether it can be applied, or whether the theory under consideration has a nontrivial BPS sector (nontrivial first-order solutions), at all.
A second method is known under the name of "first-order formalism" [13]- [18]. It essentially consists in identifying a first integral of the field theory under consideration and is, therefore, especially well adapted for one-dimensional systems, where it can easily handle the case of several fields. It can also be used in theories which are effectively one-dimensional, e.g., because the considered field configurations are co-dimension one defects, or (in some cases) because of a symmetry reduction (assuming, e.g., spherical symmetry). But in the most general higher-dimensional case, the method, again, does not provide a criterion as to whether it can be applied, i.e., whether the required first integrals can be found.
A third, rather recent method was called "on-shell method" by its inventors [19]- [21]. As it was developed up to now, the method can only be applied to effectively one-dimensional systems, where it, however, can handle the multiple-field case. To explain the method, let us consider as a specific example a theory of several fields in one dimension with energy functional E = dxE(φ a , φ a′ ) where a = 1 . . . m and φ a′ ≡ ∂ x φ a . The method then consists of the following two steps. Firstly, one tries to re-express the m EL equations where D x is the total x derivative, acting both on explicit and on implicit functions of x (e.g. Further, the functions f a and g a may, in principle, depend both on x and on the fields, but not on derivatives of the fields. For simplicity, we assume from now on that f = f (φ a ), g = g(φ a ) do not depend on x. The second step then consists in adding and subtracting m functions X a (φ a ) in the following way, The following pair of first-order equations are then sufficient conditions for the original EL equations, The applicability of the method is restricted i) by the fact that, right now, it only works in one dimension (or in effectively one-dimensional systems), and ii) by the condition that the g a must not depend on the φ a′ , which cannot always be fulfilled. Very recently, some generalisations of the method have been developed, where this last condition can be weakened [20], [21]. Before presenting a fourth method, which will be the main theme of this paper, for illustrative purposes we want to apply the methods presented so far to the simplest possible system, namely a real scalar field in 1 + 1 dimensional space-time with the standard lagrangian density where the potential U is non-negative, as always. Further, we assume for the moment that U has two zeros at the vacuum values φ = φ 1 , φ 2 (φ 2 > φ 1 without loss of generality). For static field configurations, this leads to the energy functional The Bogomolnyi trick just requires to complete the square, and the mixed term Q is indeed a homotopy invariant, as it must be. The value of Q depends on the imposed boundary conditions. Finite energy requires that both φ + and φ − take one of the two vacuum values φ 1 or φ 2 , which leads to the values Q = 0 (trivial or vacuum solution), . The corresponding BPS equation just reads φ ′ = ± √ 2U = W ,φ , and W is usually called the superpotential. The first-order formalism for the simple system at hand just boils down to the observation that the EL equation φ ′′ = U ,φ may be integrated (after multiplication by φ ′ ) to the equation and we recover the BPS equation. Finally, the on-shell method introduces the function X(φ) by adding and subtracting D x X in the EL equation Inserting φ ′ from the first equation into the second and integrating the last equation leads to U = (1/2)X 2 + const, but finite energy requires const = 0, so X = ± √ 2U and φ ′ = ± √ 2U , and we recover the BPS equation, again.
The fourth method we want to consider was proposed under the names of "strong necessary conditions" or "Bogomolnyi decomposition" by its inventors [22]- [30]. It is the main purpose of the present paper to generalise and further develop this method for the order reduction of Euler-Lagrange equations, to review some known applications, and to apply it to new nonlinear systems. For reasons which will become rather obvious in a moment, we prefer to call this method the "First-Order Euler-Lagrange formalism" (FOEL formalism) and the resulting order-reduced field equations the "First-Order Euler-Lagrange equations" (FOEL equations). We want to emphasize already at this point that the FOEL method i) is completely general, i.e, it may be applied to all systems which allow for a reduction of order and, ii) is systematic, i.e., requires (almost) no guesswork. In particular, it provides an alternative -and much more systematic -derivation of all known Bogomolnyi equations of nonlinear soliton-supporting field theories, as well as Bäcklund transformations of certain 1+1 dimensional field theories, among other results, thereby demonstrating both its usefulness and its versatile character. This holds true despite the fact that the method is based on a combination of two very simple (in fact, almost trivial) observations, as we shall explain in the next section.
The paper is organised as follows. In Section 2, we introduce the FOEL formalism in its most general form. In Section 3 we consider various examples for its application, for field theories in 1+1 dimensions, 2+1 dimensions, as well as for field theories coupled to gravity. Section 4 contains our conclusions. We always assume the speed of light equal to one, c = 1. In Minkowski space, we use the metric sign convention ds 2 = dt 2 − d x 2 . Further, in all examples we assume that some units of length and energy (or action) have been fixed, such that both our coordinates x µ and our fields φ a are dimensionless. All coupling constants which may appear in some examples are, therefore, dimensionless, as well.

II. THE FIRST-ORDER EULER-LAGRANGE FORMALISM
To explain the two simple observations which provide the starting point of the method, let us, for the moment, consider a theory of real scalar fields φ a with an action functional where m is the dimension of field space, d is the dimension of physical space (or space-time), and ∂ µ ≡ (∂/∂x µ ). The lagrangian density (energy density in the static case) is restricted to depend only on the fields and their first derivatives. The necessary generalizations for the inclusion of gauge and/or gravitational fields will be presented when required.
then this new action leads to the same EL equations as the old action S. Here, the functions J µ are, in general, functions of the coordinates x µ , the fields φ a and their first derivatives φ a ,µ , and the D µ are total derivatives, Here, repeated indices are summed over (Einstein summation convention). The important point is that, in contrast to the second-order EL equations, the first-order EL (FOEL) equations (II.3) are not invariant, so by appropriately choosing J µ (i.e.,L in (II.3)), we may obtain nontrivial FOEL equations (e.g., Bogomolnyi equations) with nontrivial solutions (e.g., BPS solitons). From Eq. (II.5) it seems that the new lagrangianL will contain second derivatives, which we do not want to permit. If min(d, m) = 1, this is indeed the case, so J must be restricted to depend only on x and φ a (or J µ only on x µ and φ), but not on first derivatives of the fields. For min(d, m) > 1, on the other hand, there exist certain antisymmetric combinations of first derivatives such that the unwanted second derivatives φ a ,µν cancel. Let us consider the simplest nontrivial case m = d = 2 more explicitly. The most general expression for the functions J µ is (using x 1 ≡ x, x 2 ≡ y and the summation convention w.r.t. b) leading to (here, ∇ ≡ (∂ x , ∂ y )) and, indeed, terms containing either φ a ,xy or φ 1 ,x φ 1 ,y , etc., have cancelled. Here, several comments are in order. 1) If H b = H b (φ a ) does not depend explicitly on x µ , then ∇H b = 0 and the above expression simplifies. Further, H 1 and H 2 only enter in the combination G(φ a ) = (H 2 ,φ 1 − H 1 ,φ 2 ), so the above total divergence simplifies to and, with the restriction H b = H b (φ a ), this is the most general total derivative term which may be added to a lagrangian density for m = d = 2.
2) The expression G(φ a ) φ 1 ,x φ 2 ,y − φ 1 ,y φ 2 ,x is precisely (proportional to) the topological charge density of two-dimensional nonlinear field theories supporting topological solitons. So it is not surprising that this term will be important in the derivation of the Bogomolnyi equations of said theories. 3) Before generalizing to higher dimensions, it is useful to introduce a more compact notation. Defining K 1 = H 2 , K 2 = −H 1 , J µ may be expressed in the compact notation (it turns out that the K a are more suitable for generalizations than the H a ). The total derivative then is (assuming, for the moment, general functions K a (x µ , φ a )) (here and below we use the notation φ a ,µ ≡ ∂ x µ φ a and K c ,b ≡ ∂ φ b K c ). Now, the generalization to higher dimensions d and m is like follows. The most general expression for J µ reads where both F µ and the K's in general depend on x µ and φ a . Further, the K's are antisymmetric tensors both in physical space and in field space. If we assume, in addition, that the K's only depend on the fields φ a and not explicitly on the coordinates x µ (as will be the case in all our applications), then the total divergence of J µ may be expressed like where the F µ 1 ...µ j a 1 ...a j (φ a ) are tensors which are completely antisymmetric both in the coordinate space and in the field space indices. In general, the expression for the F µ 1 ...µ j a 1 ...a j tensors in terms of the K's is rather complicated and given by where the sign σ(d, m, j) (which is irrelevant for applications -we just show it for completeness) is (the derivation is relegated to appendix A). Here, the subindex ,a 1 means the ∂ φ a 1 derivative of the K's, and the bracket means antisymmetrisation w.r.t. the enclosed indices (but remember that the ǫ tensor is already antisymmetric, so the antisymmetrisation is only w.r.t. a 1 ). Fortunately, in the simplest case d = m = j (which is the case which is relevant, e.g., for topological solitons), the expression for F is very simple, where G is an arbitrary function of the fields (formally, in terms of , as easily follows from the general formula (II.14) and the Schouten identity). The above expression is, in fact, the most general completely antisymmetric tensor of maximal rank in both spaces (antisymmetric tensors of maximal rank are essentially given by one function, multiplied by the corresponding ǫ tensors).
The possibility to express the total derivative D µ J µ (locally) by an arbitrary antisymmetric tensor (without having to bother about its relation to the K's) continues, in fact, to hold for j = m, even for d ≥ m, i.e., where F µ 1 ...µm a 1 ...am (φ a ) is an arbitrary tensor-valued function of φ a which is completely antisymmetric both in the coordinate and in the field space indices. This is proven in appendix B. In all our explicit applications, the total derivatives we need to consider are of the above type (II.17), so we never have to worry about the cumbersome formula (II.14).
We further remark that, in principle, already the slightly more general equations are sufficient conditions for the EL equations, where the C µ a are some constants. These equations may, however, be generated from the standard FOEL equations (II.3) by the addition of the further ,µ , so this case is, in fact, covered by the standard FOEL equations.
Finally, let us remark that there is one significant difference between d = m = 1 and max(d, m) > 1. For d = m = 1, the number of FOEL equations (two) equals the number of unknowns φ and F , therefore we always expect to find at least local solutions (which may or may not be extendable to the desired global solutions). For max(d, m) > 1, on the other hand, the number of FOEL equations is, in general, bigger than the number of unknowns φ a , F µ and F µ 1 ...µm a 1 ...am . To find solutions one, therefore, has to assume that not all FOEL equations are independent, which introduces certain additional constraints. The FOEL method produces nontrivial solutions precisely for those field theories where these additional constraints can be imposed consistently.

III. APPLICATIONS OF THE FOEL METHOD
A. 1 + 1 dimensional field theories In a first example, for illustrative purposes, we apply the FOEL formalism to the simple case of one static standard scalar field. Then we consider the generalisations to generalised dynamics and to several scalar fields, providing an explicit example for each case. Finally, we briefly review the simple derivation of Bäcklund transformations using the FOEL formalism.

Real scalar field
First of all, we want to apply the method to the simplest case, that is, the standard field theory of one real scalar field, (I.6), which, obviously, has been done before [26]. If we calculated the FOEL equations directly for the energy density of the static energy functional (I.7), we would just find that is, the trivial solution of a field sitting in one of the extrema of U (one of the vacua if the condition of finite energy is imposed) for all x. Instead, we add a total derivative term −D x F to the static energy functional (I.7), where, for simplicity, we assume that F only depends on φ and not on x (the minus sign in front of the total derivative is for convenience). The two resulting FOEL equations are Inserting the first equation into the second leads to Finite energy requires the constant to be zero, C = 0, leading to which is just the Bogomolnyi equation. Further, F may be identified with the superpotential, F = W . Finally, for the on-shell value of the energy (i.e., for the energy evaluated for a FOEL solution) we find (the vertical bar indicates evaluation at the FOEL solution) and, therefore, for the original energy, As a simple, explicit example, we choose the well-known φ 4 kink with potential U = (1/2)(1 − φ 2 ) 2 with two vacua at φ ± = ±1. Eq. (III.4) then leads to which provides the kink/antikink solutions φ = ± tanh(x − x 0 ) (here, the integration constant x 0 provides the kink position). Further, F = φ − (1/3)φ 3 , leading to the well-known energy result

Generalised Dynamics
We continue with the case of one real scalar field in 1+1 dimensions where now we allow, however, for lagrangian densities L(X, φ) which are rather general functions of the scalar field φ and the Poincare-invariant combination Theories of this type are known under the names of "generalised dynamics" or "k field theories" (k stands for kinetic). For simplicity, we shall again only consider the static case, such that the energy density is E(Y, φ) = −L(−X, φ), where we use the new kinetic variable Y ≡ −X = (1/2)φ ′ 2 for convenience. As always, we add a total derivative to the energy density, leading to the FOEL equations where we used φ ′ = √ 2Y . This equation may be integrated once to give , and the Y partial derivative of this expression coincides with the first FOEL equation only for C = 0). Eliminating F ,φ from Eqs. (III.10), (III.12) leads to which is just the first integral of the first-order formalism for generalised dynamics [14], [17]. Physically, this relation is known as the "zero pressure condition" or the "zero strain condition" [14], [17], [31], because the l.h.s expression in Eq. (III.13) is the pressure component of the energymomentum tensor (equally, the only strain component) in 1+1 dimensions [32]. Finally, the energy density for FOEL solutions is E = √ 2Y F ,φ = φ ′ F ,φ , leading to the simple energy expression in terms of the function F , as in the case of standard dynamics. We remark that the simplicity and the systematic character of the FOEL method is borne out in this case by the simple derivation of the first-order equations and the energy expression. The explicit solution of the first-order equations for a particular model of generalised dynamics, on the other hand, is as difficult in the FOEL formalism as it is in any other first-order method. The first-order equations are, after all, equivalent in the different approaches. In the FOEL formalism, the solution strategy is like follows. Firstly, interpret Eq. (III.10) as an algebraic equation for φ ′ (remember that for generalised dynamics E Y depends on Y , i.e., on φ ′ = √ 2Y ). This will, in general, produce 2R roots where the Y r (F ,φ ) are R given functions (roots) of F ,φ . Secondly, for a given root r insert the corresponding Y r (F ,φ ) instead of Y in Eq. (III.12) and solve for F ,φ = F r,φ . Thirdly, insert this F r,φ back into Eq. (III. 15) and now consider this equation as a first-order ODE. The whole method is, obviously, first order, but can still be quite complicated, due to the algebraic equations (III.10) and (III.12). As a simple example, we consider the case of the simplest k field theory leading to compactons (kinks with a compact domain) [33]. The static energy density is so the potential is just the φ 4 theory potential with its two vacua at φ = ±1, but the kinetic term is the square of the standard one. The first FOEL equation is (remember (III. 17) and the once-integrated second equation (III.12) is Inserting this back into the first equation leads to with the compacton solution (we assume that the integration constant (kink position) x 0 = 0, for simplicity) Here, x c = 4 3 4 π 2 is the compacton boundary. Finally, for the function F we get and to the compacton energy (III.23)

Several Fields
Now we consider the case of several real scalar fields, where for simplicity we only consider theories which have a standard (quadratic) kinetic term but may have a non-cartesian target space metric (i.e., field theories of the nonlinear sigma model type). Adding a total derivative −D x F to the static energy density, we get where G −1 is the inverse metric. The second set of FOEL equations is which simplifies to and may be integrated to (the integration constant must be zero, as always). If we identify F with the superpotential W from other first-order approaches, then the above is the superpotential equation relating the potential U and the superpotential W . In other approaches, this equation must essentially be guessed, whereas here it is a completely straight-forward result of the FOEL method. Finally, the energy for FOEL solutions is As one particular example, we consider the kinks in a massive nonlinear sigma model originally found in [34]. The energy functional for static configurations (we are still in 1+1 dimensions!) reads where φ = (φ 1 , φ 2 , φ 3 ) is a unit vector field, φ 2 = 1, taking values in the two-sphere. The kinetic (non-linear sigma model) term is invariant under general rotations of the field vector. For ǫ = 0, the potential breaks this symmetry down to rotations about the third axis in field space, whereas for ǫ = 0, only a discrete subgroup of the target space rotations remains. It is useful to parametrise the unit vector field by two fields (longitude and latitude) like φ = (sin θ cos φ, sin θ sin φ, cos θ).
The energy density, shifted by the usual total derivative then reads and (after inserting for θ ′ , φ ′ from above) Finally, the superpotential equation (the first integral of the last two FOEL equations) is We display both the (unintegrated) FOEL equations and the superpotential equation, because the former are slightly more general than the latter (i.e., (III.36) implies (III.34) and (III.35) but not the other way round), which will be important for the case ǫ = 0. In a first step, we consider the case ǫ = 0. Then it is sufficient to consider the three equations (III.32), (III.33) and (III.36). As φ does not show up in Eq. (III.36), it is consistent to assume F = F (θ) ⇒ F ,φ = 0, which immediately leads to φ = φ 0 = const. (III.36) then gives F ,θ = ±m sin θ ⇒ F = ∓m cos θ (III. 37) which, for the plus sign (kink) immediately leads to interpolating between the vacua θ = 0 (north pole) at x = −∞ and θ = π (south pole) at x = ∞. Finally, the kink energy is Next, we assume ǫ = 0. We shall find that the only topological soliton (kink) solutions will again have a constant φ; i.e., φ ′ = 0. It is, in fact, easy to deduce this fact directly from the potential. The form of the potential implies that any topologically nontrivial field configuration with finite energy must interpolate between the north pole and the south pole (e.g. θ(−∞) = 0, θ(∞) = π for a kink-like configuration). But the suppression factor sin 2 θ in the potential then implies that the field φ may take any values at the boundaries x = ±∞. Any nontrivial φ configuration may, therefore, be deformed continuously into the configuration φ ′ = 0, which obviously lowers the energy.
We shall find, however, that for ǫ = 0 not all values φ 0 are allowed, and the allowed solutions are isolated solutions from the point of view of the FOEL equations. Indeed, the assumption F = F (θ) is incompatible with the superpotential equation (III.36), because the r.h.s. explicitly depends on φ. So to find these isolated solutions, we have to use, instead, the un-integrated FOEL equations (before replacing φ ′ by F ,φ / sin 2 θ). We find that Eq. (III.34) is compatible with φ ′ = 0 for any value of φ = φ 0 . Eq. (III.35), on the other hand, is compatible with φ ′ = 0 only for sin φ 0 cos φ 0 = 0, i.e., for φ 0 = 0, π/2, π, 3π/2. Integrating Eq. (III.34) then leads to F 2 ,θ = m 2 (1 + ǫ 2 cos 2 φ 0 ) sin 2 θ. The resulting equation for F is exactly like in the ǫ = 0 case (see Eq. (III.37)) for φ 0 = 0, π, leading to the same kink solution and energy. For φ 0 = π/2, 3π/2, instead, the equation for F reads F ,θ = ±m ′ sin θ , m ′ = m 1 + ǫ 2 (III. 40) so the corresponding solution and energy may be found by the replacement m → m ′ . As m ′ > m, it follows that the solutions for φ 0 = 0, π are true global minima, whereas the solutions for φ 0 = π/2, 3π/2 are sphaleron-type solutions, i.e., saddle points which are local maxima in the φ 0 direction, whereas they are minima w.r.t. all other directions in the (infinite-dimensional) configuration space. It is interesting to note that, in this case, the FOEL method is able to find both the minima and the sphalerons. We end this example by remarking that in this model there also exist non-topological kinks which take the same value (e.g. the north pole) for x → ±∞ [34]. Obviously, the FOEL method (or any other first-order method) is not able to find these non-topological kinks, because the corresponding energy expression is zero for non-topological kink configurations, only allowing for the trivial solution.

Bäcklund transformations
The FOEL formalism also allows for a simple derivation of Bäcklund transformations [25]. As this is rather surprising, we want to briefly review this result where, for simplicity, we consider the Sine-Gordon (SG) example with Lagragian density (for a more general discussion beyond the SG example we refer to [28]). Taking light-cone coordinates x ± = 1 2 (x ± t) we have the following Lagrangian density and EL equation Bäcklund transformations are relevant for obtaining time-dependent solutions, so our system is no longer effectively one-dimensional, which will add some further constraints (the number of FOEL equations grows rapidly with the number of dimensions). The basic idea for the derivation of Bäcklund transformations in the FOEL formalism is to duplicate the system by adding a second Sine-Gordon Lagrangian depending on a second real scalar field ψ, L = L SG (φ) + λL SG (ψ) (here λ is a real parameter). As Bäcklund transformations relate different solutions of the same SG equation, this is a rather natural step.
If we now add a total derivative of the form then, alltogether, we havē The FOEL equations resulting from the variations w.r.t. φ and ψ are whereas the variations w.r.t. the field derivatives give In particular, we find that λ must be negative. Choosing λ = −1 for simplicity, we get with the general solution Expressing everything in terms of the η ± , we are left with the following four FOEL equations, and Adding and subtracting them, and using the addition theorems for trigonometric functions, we get the two equations with the common first integral (the analog of the superpotential equation) and the obvious solution The separation constant β is usually called the Bäcklund parameter. If we insert these solutions into Eqs. (III.57), (III.58) and re-express everything in terms of φ and ψ, then we just obtain the well-known Bäcklund transformations Once again, we want to emphasize the systematic character of the FOEL calculation. Indeed, after the reduction of the number of independent equations, the remaining steps are exactly as before, i.e., replace the field derivatives η +,x + etc., by the F + ,η + , etc., and then find the first integral (the "superpotential equation") of the resulting equations.

B. 2+1 dimensional field theories
In this section, we shall consider two examples, namely the baby Skyrme model and its submodels, on the one hand, and the generalised Maxwell-Higgs model, on the other hand. The FOEL formalism (under a different name) has already been applied to the baby Skyrme model [35] (as well as its gauged version [36], which under certain conditions permits an order reduction, too [37]), whereas for the generalised Maxwell-Higgs model this calculation is new.

The baby Skyrme model
Here we review the calculation of Bogomolnyi topological solitons (baby Skyrmions) for the baby Skyrme model and its submodels, using the FOEL formalism, for details we refer to [35], [36]. The field of the baby Skyrme model takes values in the two-sphere, so may be parametrised by a unit three-vector φ. Here we prefer to use a complex scalar field w = u + iv which is related to the unit vector via stereographic projection, In terms of the real and imaginary parts u and v, the energy functional of the baby Skyrme model (here σ and τ are non-negative real constants). It turns out that, in order to find the BPS solitons, it is enough to add the topological density term as a total derivative, The resulting FOEL equations are Starting from these equations, we now want to consider different submodels and special cases. In all cases, these equations cannot be all independent, because we have 6 equations for 3 unknowns.

The CP(1) model:
The CP(1) model or nonlinear sigma model consists of the quadratic kinetic term only. In our notation, it is defined by σ = 1, τ = 0, and U = 0. In this case, adding (III.72) and (III.75), we get whereas subtracting (III.74) from (III.73) gives The BPS baby Skyrme model: The BPS baby Skyrme model is the baby Skyrme model without the quadratic term, σ = 0., In addition, we set τ = 1. Eqs. (III.72) -(III.75) are now non-linear in the field derivatives, and to make them linear we impose the following non-linear first-order equation where K is a (at the moment unknown) function of u and v. But now the four equations (III.72) -(III.75) boil down to just one equation Inserting this back inte Eqs. (III.70), (III.71) we get with the common first integral (III.83) Eliminating G, we, therefore, end up with the single nonlinear first-order equation As we have just one equation for the two unknowns u and v, there exists an infinite-dimensional solution space for each winding number, which is related to the infinitely many symmetries (the area-preserving diffeomorphisms) of the energy functional (III.68) for σ = 0. For details we refer to [38], [39], [40]. , which makes them overdetermined and, in general, no solution exists. But we still may find particular solutions for specific potentials by the following procedure. We start with a specific solution of the Cauchy-Riemann equations (a specific holomorphic function w(z)) and interpret equation (III.80) as a defining equation for K for this given holomorphic w. Then we insert the resulting K into Eqs. (III.70) and (III.71) and determine the corresponding potential U . For solutions to the Cauchy-Riemann equations it holds that u ,x v ,y − v ,x u ,y = u 2 ,x + u 2 ,y = v 2 ,x + v 2 ,y , which allows to express all kinetic terms in (III.70) and (III.71) in terms of K. Replacing also G by K, Eqs. (III.70) and (III.71) simplify to with the common first integral which now should be understood as a defining equation for U , given K.
Let us give a simple example. Choosing w = z, i.e., u = x, v = y, we get (1 + ww) 4 (III.88) that is, the so-called "holomorphic potential" [41]- [43] (holomorphic because it has the holomorphic solution w = z). Choosing w = z 2 , i.e., u = x 2 − y 2 , v = 2xy instead, we get (1 + ww) 4 (III.89) and the resulting potential has two vacua, at w = 0 (north pole) and at w = ∞ (south pole). Higher powers w = z n , n > 2 result in potentials which are no longer rational functions. Instead, they contain roots and so might not belong to the class of potentials which one wants to permit. We remark that similar BPS-type solutions on compact domains (on tori) -again leading to particular potentials -were studied in [44].

The generalised Maxwell-Higgs model
The abelian Higgs model (or Maxwell-Higgs model) is known to possess BPS vortex solutions, although an analytical expression for these solutions is not known. Recently, some generalisations have been studied within the first-order formalism [45] and using the on-shell method [20]. These generalisations are defined by the lagrangian density Further, ψ is a complex scalar field, and A µ is the gauge potential of Maxwell electrodynamics. We assume that the potential U takes its only vacuum value at |ψ| = 1, giving rise to the usual "Mexican hat" type spontaneous symmetry breaking. The function w is similar to the (here, diagonal) target space metric for non-linear sigma models, but now for a gauge theory. Finally, the function h is frequently called "dielectric function", because it generalises the dielectric constant to a field-dependent function. For static configurations we choose the temporal gauge A 0 = 0. We could now introduce the FOEL method directly for the two-dimensional static energy functional but, instead, we follow [45], [20] and perform a symmetry reduction to axially symmetric configurations first. Concretely, we introduce polar coordinates x = r cos θ, y = r sin θ and make the ansatz ψ = e inθ g(r) , n ∈ Z (III.91) and where the condition of finite energy requires the real functions a and g to obey the following boundary conditions, The static energy functional (divided by 2π for convenience; further, from now on we assume e = 1) then reads E 2π = drr h 2 a ,r r 2 + w g 2 ,r + g 2 a 2 r 2 + U = dy ha 2 ,y + w 2yg 2 ,y + 1 2 g 2 a 2 y + U (III.94) where we introduced the new variable y = r 2 . Subtracting a total derivative −D y F (g, a), the resulting energy density then reads We notice the explicit presence of different powers of the independent variable y in this expression, which has the consequence that in the purely algebraic part of the FOEL equations each power of y has to vanish independently. This is the trace left in the effectively one-dimensional functional of the more restrictive character of the FOEL equations in higher dimensions. Explicitly, varying w.r.t. the field derivatives we get the two first FOEL equations Varying w.r.t. g we find h ,g a 2 ,y + w ,g 2yg 2 ,y + 1 2 g 2 a 2 y + wg a 2 y + U ,g − F ,gg g ,y − F ,ag a ,y = 0 h ,g F ,a 2h 2 + w ,g y 1 8  Our results coincide with the ones of [20], but we believe that the method used here is simpler and more systematic.
As always, we want to end with some explicit examples. First of all, choosing h = 1 and w = 1, we recover the standard abelian Higgs model. Indeed, w = 1 implies K ,g = −2g ⇒ K = 1 − g 2 , leading to the standard abelian Higgs potential and h should be a function of g 2 in order to avoid a singularity at g = 0 for w. As a more explicit example, we may choose h = (1 + g 2 ) −m , leading to where m is a positive integer. In particular, the so constructed w is positive definite in the fundamental domain of the standard abelian Higgs vortex (i.e., in the interval 0 ≤ g ≤ 1 where the vortex takes its values), as it must be.

C. Self-gravitating field theories
Self-gravitating field theories, that is, field theories coupled to gravity in the standard way and with the Einstein-Hilbert term included are, in general, not reducible to lower order. But after some simplifying assumptions (e.g., symmetry reductions), such a reduction of order may be possible (i.e, a first integral of the field equations may exist). Two known examples where this happens are scalar field inflation and "thick brane world models", where the 3+1 dimensional universe is assumed to be a brane of finite thickness in a 4+1 dimensional bulk universe, and the finite thickness is the result of a finite extension of a soliton (a kink) in the fifth dimension. As scalar field inflation and thick brane world models are formally very similar, we shall consider only the first case. Finally, we will consider the case of the BPS Skyrme model in a curved space-time and rederive the conditions which must hold such that this system remains a BPS theory.

Scalar field inflation
Scalar field inflation is known to possess a first integral, where the methods to derive this first integral are known under the names of "Hamilton-Jacobi approach" [46], [47], "fake supersymmetry" (or "fake supergravity") [48], [49], the "superpotential method" [50], or the already considered first-order formalism [13], [15], [16]. Here we want to rederive this result using the FOEL formalism. Our starting point is the action where S EH is the Einstein-Hilbert action, L m is the matter (scalar field) lagrangian L m = (1/2)g µν ∂ µ φ∂ ν φ − U , g µν is the metric tensor, g = det g µν , and R is the Ricci scalar. Further, κ is a constant related to Newton's constant by κ = 4πG. The resulting EL equations (the Einstein equations) are compatible with the cosmological ansatz for a spatially flat universe, where H is the Hubble "constant" (the Hubble function). So the function G is essentially the Hubble function. Finally, insertingφ andȧ from Eqs. (III.122) and (III.123) into Eq. (III.121), we find the "superpotential equation" or "Hamilton-Jacobi equation" where G should be identified with the "superpotential" W from other approaches. Inserting, instead, Eqs. (III.122) and (III.123) into Eq. (III.120), we get the φ derivative of the superpotential equation, i.e, an identity. Our results coincide, of course, with the results from other methods. We want to emphasize, once more, the simple and systematic character of the FOEL method.

BPS Skyrmions on curved space-times
The Skyrme model is a nonlinear field theory in 3+1 dimensions which is considered to provide a mesonic low-energy effective action for Quantum Chromodynamics (QCD). Its field U takes values in the group manifold SU(2), U ∈ SU(2), and, physically, may be identified with the pions. The lagrangian density of the Skyrme model consists of a term quadratic in first derivatives (the "non-linear sigma model term") and a term which is quartic in first derivatives (the "Skyrme term"). Further, the original model may be generalised naturally to include both a potential term (supposed to give masses to the pions) and a term sextic in first derivatives (which we shall simply call the "sextic term"). Quite recently, it was found that within this class of generalised Skyrme models there exists a submodel which has the BPS property [51], i.e., both a BPS equation for static configurations and infinitely many solutions which satisfy the BPS equation and saturate the corresponding Bogomolnyi bound. As always, this BPS equation can be derived using the FOEL formalism [52]. This so-called "BPS Skyrme model" consists of the potential and the sextic term only (for details see [51], [53]), where c is a constant, L µ = U † ∂ µ U is the left-invariant chiral current and B µ is the baryon current (topological current). Further, we already introduced the generalisations necessary on curved space-times. In flat (Minkowski) space-time, and for potentials U = U (TrU), the BPS equations are compatible with the axially symmetric ansatz in spherical polar coordinates, where B is the baryon number (topological degree). Further, this ansatz leads to the spherically symmetric action S BPS = dt dr r 2 sin θ dθ dφ L BPS = −4π dt dr r 2 cB 2 2r 4 sin 4 f f ′2 + U (f ) (III.127) (we assume from now on that the potential U (f ) has its unique vacuum at f = 0). It turns out that the same axially symmetric ansatz (III.126) is compatible with the field equations of the full self-gravitating system for the Schwarzschild-type metric ansatz ds 2 = σ 2 (r)N (r)dt 2 − dr 2 N (r) − r 2 (dθ 2 + sin 2 θdφ 2 ) , N (r) ≡ 1 − 2m(r) r (III.128) (where we defined the "mass function" m(r) for later convenience). Generalising S BPS for this metric and adding the Einstein-Hilbert action for the same metric results in the total action (for self-gravitating Skyrmions in general, and for the EH action for this metric, we refer to [54], and for self-gravitating BPS Skyrmions to [55]- [58]). We now might try to add a total derivative D r F (f, m, σ) to the corresponding lagrangian density L tot and to apply the FOEL method. It turns out, however, that any assumption of a nontrivial F leads to a contradiction, so the only solution which the FOEL method is able to reproduce for the full self-gravitating system requires F = Cm (where C is a constant) and leads to the vacuum solution f = 0 for the Skyrme field and to the Schwarzschild solution for the metric, m = m ADM = const. and σ = κC = const. We still may pursue a less ambitious goal and consider the "BPS Skyrme model" in a fixed background metric (i.e., for fixed functions N (r) and σ(r)) and ask the question for which background metrics this system still admits a BPS equation and BPS solutions (i.e., is a genuine BPS Skyrme model). That is to say, we skip the EH term and add the total derivative D r F = F ,f f ′ , leading to the "energy density" where now f is the only dynamical field. The resulting FOEL equations are As the first term in this equation does not depend on r, the second term cannot be r-dependent, either, leading to the conclusion σ 2 N = const., i.e., the time-time component g tt of the metric must be constant. This precisely agrees with the result recently derived in [57]. Assuming this, the above equation may be integrated to the "superpotential equation" As in the flat space case, the functional form of f is completely determined by the potential U .

IV. CONCLUSIONS
It was the main purpose of the present paper to generalise and further develop a systematic method (which we called the First-Order Euler-Lagrange (FOEL) formalism) for the reductionof-order of EL equations of nonlinear field theories originally introduced in [22]- [30]. Further, we reviewed some known applications of the method and presented some new ones. Concretely, the FOEL equations for generalised dynamics and for the case of several fields in 1+1 dimensions, for the generalised Maxwell-Higgs system in 2+1 dimensions, as well as for all field theories coupled to gravity are new results. As said, the formalism applies in all cases where an order reduction may be performed, not just in the cases reviewed here. The self-duality equations for instantons, e.g., were already derived in [22]. It would, of course, be interesting to discover new field theories possessing a BPS sector using the FOEL formalism. Here, the most nontrivial part is the identification of a candidate field theory, because once such a candidate is found, the formalism provides a systematic way to find (or disprove) the BPS sector. Another question of interest concerns the relation of the FOEL formalism with supersymmetry (SUSY). It is well-known that theories with a BPS sector typically allow for SUSY extensions. Further, SUSY transformations produce a total derivative term when acting on the lagrangian density. So one wonders whether the total derivative term D µ J µ of the FOEL method is related to the total derivative term of SUSY transformations, and whether the current J µ of the FOEL method is related to the (bosonic part of the) supercurrent of the SUSY-extended theory. These and related questions shall be investigated in future publications.
Next, we perform the following changes of index names, j = i + 1, a k → a k+1 , c → a 1 and a further reshuffling to obtain min(m,d) Here, the expression in the first line multiplying the antisymmetric product of the φ a µ is already antisymmetric in µ 1 . . . µ j and in a 2 . . . a j , so all that is missing for the explicit expression for the F µ 1 ...µ j a 1 ...a j tensors is an antisymmetrisation w.r.t. a 1 , leading to Eq. (II.14). For the antisymmetrisation it is sufficient to sum over all cyclic permutations, i.e., T [a 1 ǫ a 2 ...a j ] = 1 j T a 1 ǫ a 2 ...a j + T a 2 ǫ a 3 ...a j a 1 + · · · + T a j ǫ a 1 ...a j−1 (A.4) because the expression is already antisymmetric in a 2 . . . a j .