Recursion Operators and Tri-Hamiltonian Structure of the First Heavenly Equation of Pleba\'nski

We present first heavenly equation of Pleba\'nski in a two-component evolutionary form and obtain Lagrangian and Hamiltonian representations of this system. We study all point symmetries of the two-component system and, using the inverse Noether theorem in the Hamiltonian form, obtain all the integrals of motion corresponding to each variational (Noether) symmetry. We derive two linearly independent recursion operators for symmetries of this system related by a discrete symmetry of both the two-component system and its symmetry condition. Acting by these operators on the first Hamiltonian operator $J_0$ we obtain second and third Hamiltonian operators. However, we were not able to find Hamiltonian densities corresponding to the latter two operators. Therefore, we construct two recursion operators, which are either even or odd, respectively, under the above-mentioned discrete symmetry. Acting with them on $J_0$, we generate another two Hamiltonian operators $J_+$ and $J_-$ and find the corresponding Hamiltonian densities, thus obtaining second and third Hamiltonian representations for the first heavenly equation in a two-component form. Using P. Olver's theory of the functional multi-vectors, we check that the linear combination of $J_0$, $J_+$ and $J_-$ with arbitrary constant coefficients satisfies Jacobi identities. Since their skew symmetry is obvious, these three operators are compatible Hamiltonian operators and hence we obtain a tri-Hamiltonian representation of the first heavenly equation. Our well-founded conjecture applied here is that P. Olver's method works fine for nonlocal operators and our proof of the Jacobi identities and bi-Hamiltonian structures crucially depends on the validity of this conjecture.


Introduction
Plebański had demonstrated how the Einstein field equations in the complex domain under (anti-)self-duality condition reduce to a single scalar PDE, which he called the first and the second heavenly equations [26]. Earlier we had shown that the second heavenly equation is a bi-Hamiltonian system [19]. We had also obtained bi-Hamiltonian representations for all other equations of the heavenly type [20,32,33,34] except the first heavenly equation.
Here we show that the first heavenly equation of Plebański presented in the two-component evolutionary form where t =t +ỹ, x =t −ỹ, possesses two linearly independent recursion operators for symmetries. (From now on, letter subscripts denote partial derivatives with respect to corresponding variables.) It should be noted here that recursion operators and Hamiltonian structures can also be found in the initial, non-evolutionary form of a given PDE by the methods of [10].
In Section 2, we present the Lagrangian for the system (1.2). Lagrangian of this system turns out to be degenerate and therefore we use the Dirac's theory of constraints [5]. The Poisson bracket of Dirac constraints provides us a symplectic operator and symplectic 2-form. The inverse of the symplectic operator is the first Hamiltonian operator J 0 for the two-component system. We find the corresponding Hamiltonian density and present our system in a Hamiltonian form.
In Section 3, we obtain all point Lie symmetries of the first heavenly system (FHS) together with the table of commutators of the basis Lie algebra generators and show how the corresponding Hamiltonians of the variational symmetry flows can be obtained from the inverse Noether theorem in a Hamiltonian form. These Hamiltonians are integrals of the motion along the first heavenly flow for all the variational symmetry flows that commute with FHS. Each Hamiltonian is also conserved along all the symmetry flows that commute with the flow generated by this Hamiltonian.
In Section 4, we show that the symmetry condition for the first heavenly equation (1.1) in a one component form can be compactly expressed in terms of Lax operators and have a two dimensional total divergence form and therefore defines a potential variable. We prove that the potential of a symmetry is again a symmetry and thereby a recursion relation for symmetries is derived. These are partner symmetries which we introduced in [15] for the complex Monge-Ampère equation and later classified heavenly equations possessing partner symmetries in [30], with particular applications in [19] for the Plebański second heavenly equation and in [33] for the mixed heavenly and Husain systems. In [16,31] we have demonstrated how to use partner symmetries for constructing noninvariant solutions of the complex Monge-Ampère equation and corresponding (anti-)self-dual Ricci-flat metrics without Killing vectors which are important for discovering an explicit metric of the gravitational instanton K3.
We have obtained the recursion operator R 1 for symmetries in a two-component form. We also discover a second recursion operator R 2 by applying to R 1 a discrete symmetry of the system (1.2) and its symmetry condition. We also construct the recursion operators R + = R 1 + R 2 and R − = R 1 − R 2 that are even and odd, respectively, under the discrete symmetry. The idea of constructing recursion operators from the Lax-type representation in the context of partner symmetries appeared in our paper [15] and later, in a more general context, in the paper [28]. Here we derive recursion operators using our approach of partner symmetries while A. Sergyeyev in [28] uses a more general approach for deriving recursion operators which works even when partner symmetries do not exist (see, e.g., our paper [32]).
In Section 5, acting by these recursion operators on J 0 we obtain several alternatives for the second Hamiltonian operator for the two-component system, We have failed to find the Hamiltonian densities for J 1 and J 2 such that provide bi-Hamiltonian representations of the two-component system. On the contrary, for the even and odd Hamiltonian operators J ε we have found the appropriate Hamiltonian densities and obtained two bi-Hamiltonian representations of the first heavenly equation of Plebański in a two-component form (1.2). Using jointly all three Hamiltonian operators J 0 , J + and J − we obtain tri-Hamiltonian representation of the system (1.2), provided that the Poisson compatibility of these operators is proved.
In Section 6, we prove that the Jacobi identities are satisfied for an arbitrary linear combination of the three operators and hereby prove that they are indeed Hamiltonian and compatible. For this purpose we use P. Olver's theory of functional multi-vectors [21]. P. Olver's theory was formulated for matrix-differential Hamiltonian operators while our operators J ε are nonlocal. We show how this theory works nicely even for these nonlocal operators, so that the applicability of P. Olver's method to nonlocal Hamiltonian operators seems to be a well-founded conjecture though a rigorous generalization of his theory to nonlocal Hamiltonian operators is still awaited. Under this conjecture, we have shown that equations (1.2) yield a tri-Hamiltonian completely integrable system in the sense of Magri [13,14].
It should be noted that there is a number of tri-Hamiltonian systems (and even multi-Hamiltonian systems) in the literature [1,2,3] and [11]. Our results show that a tri-Hamiltonian system can be found also in the multi-dimensional case.
We also note that the Olver's method for for checking Jacobi identity is not the only one that works for nonlocal operators. Many authors use the formula by Gel'fand, Dorfman and co-workers. Recent examples of (1+1)-dimensional systems treated in this way can be found in [12].

Lagrangian, symplectic and Hamiltonian structures
We obtain a Lagrangian for this two-component system which is degenerate in the following sense. The canonical momenta which satisfy canonical Poisson brackets, cannot be inverted for velocities. Following Dirac's theory of constraints [5], we impose (2.2) as constraints and calculate the Poisson bracket of the left-hand sides of the constraints (2.3) The result is expressed in the matrix-operator form as which is obviously skew symmetric. From now on D ξ denotes total derivative operator with respect to ξ. The corresponding 2-form Ω is obtained by integrating the density which is closed up to a total divergence so that Ω = V ωdV is closed under appropriate boundary conditions and hereby Ω is a symplectic form, while K in (2.4) is a symplectic operator. Taking the inverse of K we obtain Hamiltonian operator which is obviously skew symmetric, while the Jacobi identities for J 0 are satisfied as a consequence of the closure of the symplectic 2-form Ω. The corresponding Hamiltonian density H 1 is found by using (2.1) and (2.2) in the relation Thus, we present the two-component system (1.2) in the Hamiltonian form where δ u , δ v denote Euler-Lagrange operators related to variational derivatives of Hamiltonian functional [21].

Symmetries and integrals of motion
Point Lie symmetries of the system (1.2) are determined by the following basis generators where a, b and f , g are arbitrary smooth functions of one and two variables, respectively and subscripts denote partial derivatives. We note that the obvious translation invariance symmetries are generated by X 1 , X 2 and also by Y a=1 ± Z b=1 . In particular, X = (Z b=1 − Y a=1 )/2 = −∂ t is the generator of the first heavenly flow (1.2) itself.
The Lie algebra of point symmetries is determined by the Table 1 of commutators of the basis generators where the commutator [X i , X j ] stands at the intersection of ith row and jth column. Here we use the following shorthand notatioñ and the primes denote derivatives of arbitrary functions a and b of a single variable.
We need symmetry characteristics determining symmetries in evolutionary form [21] with the independent variables not being transformed under symmetry transformations. For the point symmetry generator of the form X = ξ i ∂ x i + η α ∂ u α , where the summation over repeated indices is used, the symmetry characteristics are defined as ϕ α = η α − u α i ξ i with the subscripts i denoting derivatives with respect to x i . In our problem, α = 1, 2, u 1 = u, u 2 = v, η 1 = η u , determine the transformation of u and v, respectively. We also use u t = v and v t = Q where Q is the right-hand side of the second of our equations (1.2). Symmetry characteristics have the form Applying the formula (3.2) to the generators in (3.1), we obtain characteristics of these symmetries In accordance with our previous remark that X = (Z b=1 − Y a=1 )/2 = −∂ t generates the first heavenly flow, we note that (ϕ b=1 − ϕ a=1 )/2 = v and (ψ b=1 − ψ a=1 )/2 = Q are characteristics of this symmetry with the group parameter τ = t.
First Hamiltonian structure provides a link between symmetries in evolutionary form and integrals of motion conserved by the Hamiltonian flow (2.7). Replacing time t by the group parameter τ in (2.7) and using u τ = ϕ, v τ = ψ for symmetries in the evolutionary form, we obtain the Hamiltonian form of the Noether theorem for any conserved density H of an integral of motion To determine the integral H that corresponds to a known symmetry with the characteristic (ϕ, ψ) we use the inverse Noether theorem Here (3.5) is obtained by applying the operator K to both sides of (3.4).
We now apply the formula (3.5) to determine integrals H i corresponding to all variational symmetries with characteristics (ϕ i , ψ i ) from the list (3.3). Using the expression (2.4) for K, we rewrite the formula (3.5) in an explicit form which supplies the formulas for determining integrals H i for the known symmetries (ϕ i , ψ i ) Here K 11 is determined from (2.4) to be K 11 = (vz + u xz )Dx + (vx − u xx )Dz + vxz. We always start with solving the second equation in (3.6) in which, due to the fact that ϕ i never contain derivatives of v, δ v H i is reduced to the partial derivative H i v , so that the equation H i v = uxzϕ i is easily integrated with respect to v with the "constant of integration" h i [u] depending only on u and its derivatives. Then the operator δ u is applied to the resulting H i , which involves the unknown δ u h i [u], and then it is equated to δ u H i following from the first equation in (3.6) to determine δ u h i [u]. Finally, we reconstruct h i [u] and hence H i . If we encounter a contradiction, then this particular symmetry is not a variational one and does not lead to an integral.
This solution algorithm for the symmetries listed in (3.1) with characteristics (3.3) yields the following results. X 5 generates a non-variational symmetry. For all other symmetries the corresponding integrals read with arbitrary smooth functions a(x), b(z), f (t + x,x), g(t − x,z). According to our previous remarks, the Hamiltonian H 1 of the first heavenly flow, obtained in (2.6), is also contained in the set (3.7) as a linear combination H 1 = (H b=1 − H a=1 )/2. Since the symmetry generators X 1 , X 2 , X 3 , Y a , Z b commute with the generator −∂ t of the first heavenly flow, the corresponding Hamiltonian densities H 1 , H 2 , H 3 , H a , Z b are densities of integrals of motion for the flow (1.2), subject to suitable boundary conditions. Moreover, each H i is the density of the integral of motion along all the symmetry flows with characteristics (3.3) that commute with the corresponding X i .
Converting the 2-component first heavenly system back to the first heavenly equation in the one-component form, one could obtain its integrals as reductions of the integrals for the system.

Lax pairs and recursion operators
We start with the truncated Lax pair, with the spectral parameter λ = 0, for the first heavenly equation (1.1) (see [6]) which commute on solutions of (1.1). In the new variables t =t +ỹ, x =t −ỹ, used in the two-component system (1.2), they become Let ϕ be a symmetry characteristic [21] of the first heavenly equation (1.1) in a one component form, so that the Lie equation reads u τ = ϕ, where τ is the group parameter and independent variables do not transform under the group. The equation (1.1), transformed to the new variables, has the form with the truncated Lax pair for this equation given in (4.1). Symmetry condition for equation (4.2) reads It can be expressed in terms of Lax operators (4.1) as and hereby have a two-dimensional divergence form. Hence it defines a potential variableφ by the equations It is easy to check that the symmetry condition can also be put in the form This equation is satisfied for the potentialφ. Indeed, using its definition (4.4) we have since L 1 and L 2 commute on solutions of the equation (4.2). Thus, the potential of a symmetry is again a symmetry and we obtain a recursion relation for symmetries.
In an explicit form the definition (4.4) of the symmetry potential becomes where v = u t , and using the definition of Q in (1.2), we definẽ For the two-component form (1.2) of our equation, we introduce the two-component symmetry characteristic (ϕ, ψ), where ψ = ϕ t and the same for the potential of the symmetry (φ,ψ) with ψ =φ t , and use this notation in (4.5) with the result We solve the first equation here forφ in the form and use this in the second equation in (4.6) to obtaiñ Then we obtain our first recursion operator R 1 in 2×2 matrix form which acts on two-component symmetry characteristics with the following explicit form for R 1 We also discover a second recursion operator by noticing the discrete symmetry possessed both by the first heavenly equation in the two-component form and its symmetry condition. For equation (1.1) we choose among its discrete symmetries the symmetry of exchangingx ↔ −z andt ↔ỹ, which is inherited by the two-component system (1.2) in the form We note that this discrete symmetry also leaves invariant the symmetry condition for the first heavenly equation (4.3). The alternative Lax pair for this equation becomes and the symmetry condition is now expressed as Repeating the reasoning that has led us to the first recursion operator, we discover that the symmetry potential which is determined by the conditionŝ also satisfies the symmetry condition and hereby is also a symmetry and hence we have obtained another recursion for symmetries. The second recursion operator in a matrix 2×2 form, obtained from R 1 by the same discrete symmetry transformation, reads where we have used Dỹ = D t − D x for a two-component representation and We have added extra overall minus in R 2 to have the action of the discrete symmetry to be For future use, we will need the recursion operators either even or odd relative to the discrete symmetry, i.e., R + = R 1 +R 2 and R − = R 1 −R 2 . Joining both cases we introduce R ε = R 1 +εR 2 with the explicit form where ε = ±,

Bi-Hamiltonian representations of the f irst heavenly equation in two-component form
Composing the recursion operator R 1 with the first Hamiltonian operator J 0 in (2.5) we obtain a candidate for a second Hamiltonian operator J 1 for the two-component system The operator J 1 is obviously skew symmetric, so the remaining task is to check Jacobi identities which will be considered in the next section. Similarly, composing the recursion operator R 2 with the Hamiltonian operator J 0 we obtain a candidate for a third Hamiltonian operator J 2 for the two-component system, which can be also obtained directly by applying the discrete symmetry transformation (4.7) to J 1 The operator J 2 is also obviously skew symmetric, so it remains only to check the Jacobi identities. However, using simple natural ansatzes, we failed to find Hamiltonian functions corresponding to the operators J 1 and J 2 such that acting by these operators on the variational derivatives of the Hamiltonian functions we could generate the two-component flow (1.2). Therefore, we considered operators J + and J − generated from J 0 by the recursion operators R + and R − , respectively, which have the additional property of being either even or odd under the discrete symmetry transformation (4.7). Using simultaneously both even and odd operators in the form J ε = R ε J 0 = J 1 + εJ 2 , where ε = ± and using the formula (4.9) for R ε , we obtain the following result The operator J ε is also obviously skew symmetric. Jacobi identities for J ε and compatibility of it with J 0 and also between J + and J − are proved in the next section, so that all of them are mutually compatible Hamiltonian operators. For the operator J ε we have found the corresponding Hamiltonian densities H 0ε using a simple and natural ansatz of the density to be linear in v and moreover where h 0ε [u] depends only on u and its partial derivatives. Then the coefficient functions h 1ε and h 0ε are determined from the requirement that the two-component system (1.2) should admit two bi-Hamiltonian representations for ε = ± with the Hamiltonian operators J ε provided that we have proved the Hamiltonian property for the operators J + and J − and compatibility of the Hamiltonian operators J 0 , J + and J − . These operators were automatically Hamiltonian if we would know that the recursion operators R + and R − are hereditary (Nijenhuis) [9,29], but this property of our recursion operators is not known. For this reason, in section 6 we check that the three operators J 0 , J + and J − are mutually compatible Hamiltonian operators (form a Poisson pencil ). For this purpose we use the technique of the functional multi-vectors from P. Olver's book [21].

Proof of Jacobi identities and compatibility of Hamiltonian structures
Compatibility of three Hamiltonian structures means that the linear combination of the three Hamiltonian operators with arbitrary constant coefficients is also a Hamiltonian operator.
Since J 0 , J + and J − are obviously skew symmetric, the remaining problem is to prove the Jacobi identities for the linear combination J = aJ + + bJ − + cJ 0 of the three Hamiltonian operators with arbitrary constant coefficients a, b and c. Thus, we check simultaneously that J + and J − are indeed Hamiltonian operators, at b = c = 0 and a = c = 0 respectively, that the three Hamiltonian operators J + , J − and J 0 are compatible and the system (1.2) in the form (5.3) is tri-Hamiltonian. Since J + = J 1 + J 2 and J − = J 1 − J 2 , we obviously can equivalently change the definition of J to J = aJ 1 + bJ 2 + cJ 0 in our analysis of the Jacobi identities for J to simplify the calculations.
To prove the Jacobi identities for J, we use the technique of P. Olver's book [21], so below we give a short summary of the notation and results from this book. Let A l be the vector space of l-component differential functions that depend on independent and dependent variables of the problem and also on partial derivatives of the dependent variables up to some fixed order. for all functionals P, Q and R, where δ is the variational derivative. However, the direct verification of the Jacobi identity (6.1) is a hopelessly complicated computational task. For this reason we will use P. Olver's theory of the functional multi-vectors, in particular, his criterion: Here u i 0 = u i and in our case i, j = 1, 2, while u 1 = u and u 2 = v, and ω = (ω 1 , ω 2 ) = (η, θ) is a functional one-form corresponding to a "uni-vector" with the following property for the action of total derivatives D J (ω i ) = ω i J . By the definition of the space of functional multi-vectors, the integrals of total divergences in both Θ and pr v Dω (Θ) always vanish, so that we can integrate by parts terms in the integrands, if needed. We note also that by definition of the prolonged evolutionary vector field pr v Dω , it commutes with total derivatives and annihilates uni-vectors pr v Dω ω i = 0.
In the following it is convenient to introduce the short-hand notation µ = vz+u xz , ν = vx−u xx and A = uxz, so that 2A x = µx − νz and 2vxz = µx + νz. To check the Jacobi identities for the operator J = aJ 1 + bJ 2 + cJ 0 we set D = J, where where The bi-vector Θ in the theorem above has the form According to P. Olver's criterion (6.2) with D = J applied to Θ in (6.6), the condition for Jacobi identities to be satisfied reads where pr v Jω is the prolonged evolutionary vector field with the characteristic Jω which acts on each term in the integrand of (6.6). We will further skip the integral sign in the condition (6.7), keeping in mind the possibility of integration by parts while always omitting total divergences, and leaving out the characteristic Jω in the notation pr v for the prolonged vector field in (6.7) for brevity. The criterion (6.7) becomes where "mod DIV" (skipped in the following) means that the left-hand side of (6.8) equal to a total divergence is equivalent to zero. The action of the vector field pr v Jω involved in (6.8) is defined in terms of the matrix elements of the operator J in (6.4), (6.5) according to the following rules (with u = u 1 , v = u 2 and the letter subscripts denoting partial derivatives) pr v Jω (µ) = pr v Jω (vz + u xz ) = Dz J 21 η + J 22 θ + D x J 11 η + J 12 θ With these prolongation formulas plugged in the criterion (6.8), three groups of terms, bilinear in η and its derivatives, linear in η and without η, should separately either vanish or be reducible to a total divergence form using integration by parts when necessary. The operator J = aJ 1 + bJ 2 + cJ 0 should satisfy the Jacobi identities for arbitrary parameters a, b and c and hence vanishing of terms, up to a total divergence, should occur separately in each subgroup with distinct bilinear dependence on the coefficients a, b and c. Moreover, each subgroup is further divided into subgroups of terms with the trilinear, bilinear, linear dependence on µ, ν and its derivatives or terms without µ and ν and each such subgroup also should separately vanish or be reducible to a total divergence. Sometimes such cancelations of terms will occur only on account of the identities 2A x = µx − νz and 2vxz = µx + νz following from definitions of µ, ν and A given before the formula (6.4).
It is easy to check that all terms bilinear in η and its derivatives are canceled. For terms linear in η we consider separately groups containing a 2 , b 2 , c 2 , ac, bc and ab. For example, terms with a 2 without µ (there are no terms with ν in this group) are combined into the single term (we set a = 1 here) which is a total divergence. Terms linear in µ combine to Terms quadratic in µ are where the cancelations are obvious after expanding all derivatives of the products.
As another example, we consider terms linear in η containing ab. Terms without µ and ν do not cancel completely, the remainder being (we skip ab in the following) which should be joined with the terms linear in µ and ν. All such terms finally combine to the following expressions where we have used integration by parts skipping total divergencies and, at the large step, the relation 2A x = µx − νz. The remaining term is bilinear in µ and ν and it should be joined to other bilinear terms, linear in η and containing ab, which jointly yield which either vanish or cancel after integrating by parts.
In a similar way, we check the cancelation of linear in η terms proportional to b 2 , c 2 , ac and bc and also distinct cancelation of groups of terms without η that contain each of bilinear combinations of the coefficients a, b and c. Each such group is subdivided into subgroups with different powers of µ, ν and their derivatives. The calculations are straightforward but much too lengthy to be presented here, especially for the groups of terms without η.
Thus, the Jacobi identities are satisfied for the linear combination J = aJ 1 + bJ 2 + cJ 0 and equivalently for J = aJ + + bJ − + cJ 0 with arbitrary constant coefficients which proves that all the three operators J + , J − and J 0 are Hamiltonian and mutually compatible as soon as P. Olver's criterion is applicable. However, P. Olver's criterion is formulated for matrixdifferential operators, while our operators J + and J − are nonlocal. An extensive literature exists on the theory of nonlocal Hamiltonian operators in 1 + 1 dimensions, e.g., [4,7,8,18,27]. However, to the authors' knowledge, no effective methods exist for checking the Jacobi identities in the multi-dimensional case except O. Mokhov's method in [17] whose efficiency depends on the level of complexity of the nonlocality. As a hint to the applicability of P. Olver's method for nonlocal Hamiltonian operators, we can regard the papers [22,25], where the Jacobi identities for nonlocal symmetries were shown to be correct if nonlocal variables are included in symmetries' characteristics, which imply "ghost" terms in the commutators. Since according to P. Olver's method nonlocal terms are automatically included in the characteristic of the evolutionary vector field pr v Jω , with nonlocal J, and all such terms are canceled in the process of application of the criterion, we believe that the criterion works correctly in this more general case, which we consider as a well-founded conjecture. Thus, under this conjecture we have shown that the first heavenly equation in the two component evolutionary form (1.2), being converted to the form (5.3), is indeed a tri-Hamiltonian system.
Of course, a rigorous formulation of this method for checking the Jacobi identities for nonlocal Hamiltonian operators is still awaited and could be a very worthwhile project [24].

Conclusion
We have converted the first heavenly equation into a two-component evolutionary form. Using two different Lax pairs derived for the one-component form of this equation, we have constructed two independent recursion operators R 1 and R 2 for symmetries of the first heavenly system (FHS) from the Lax operators. They are related by a discrete symmetry of both the FHS and its symmetry condition. We have obtained a Lagrangian for the two-component FHS. Applying to this degenerate Lagrangian the Dirac's theory of constraints, we have obtained a symplectic operator and its inverse, the latter being a Hamiltonian operator J 0 . We have found the corresponding Hamiltonian density, thus presenting the first heavenly system in a Hamiltonian form.
We have determined all local Lie point symmetries of the FHS system and, using the inverse Noether theorem in Hamiltonian form, we obtained Hamiltonians generating all the variational (Noether) point symmetries. These Hamiltonians are integrals of the motion along the first heavenly flow if the symmetry flows generated by them commute with the FHS flow. Each Hamiltonian generating a variational point symmetry flow is conserved along each point symmetry flow that commutes with the flow generated by the Hamiltonian under study. Converting the FH system back to the first heavenly equation in the one-component form, one could obtain its integrals as reductions of the integrals for the system.
Composing each of the recursion operators with J 0 , we have obtained candidates for two more Hamiltonian operators J 1 = R 1 J 0 and J 2 = R 2 J 0 which are related by the discrete symmetry transformation between the recursion operators. However, using some simple and natural ansatzes for densities we failed to find the Hamiltonian densities corresponding to J 1 and J 2 , so that these operators cannot be used to convert the first heavenly flow in new Hamiltonian forms. Therefore, we introduced operators J + = J 1 + J 2 and J − = J 1 − J 2 , even and odd with respect to the discrete symmetry, respectively. For these operators we have succeeded to find corresponding Hamiltonian densities H 0+ and H 0− by using a simple natural ansatz and hence we have obtained two more Hamiltonian representations for the first heavenly flow, provided that we prove that J + and J − are indeed Hamiltonian operators. These operators were automatically Hamiltonian if we would know that the recursion operators are hereditary (Nijenhuis), but this property of our recursion operators is not known. Therefore, since operators J 0 , J + and J − are obviously skew-symmetric, the remaining problem is to check directly the Jacobi identities for these operators and prove the compatibility of the three Hamiltonian operators J 0 , J + and J − . For this purpose we checked the Jacobi identities for the linear combination of these operators with arbitrary constant coefficients. To prove the Jacobi identities directly appears to be "a hopelessly complicated computational task" [21]. Therefore, we have applied P. Olver's theory of functional multi-vectors from his book [21]. In this way, under the well-founded conjecture of applicability of the Olver's method to nonlocal Hamiltonian operators, we have proved that all three operators are Hamiltonian (which has been already proved independently for J 0 ) and that all these operators are compatible. Hence, our main result is that the first heavenly equation in an evolutionary two-component form is a tri-Hamiltonian system. Earlier examples of tri-Hamiltonian systems, bidirectional Boussinesq system and three mutually compatible Hamiltonian structures associated with the Korteweg-de Vries and Camassa-Holm bi-Hamiltonian flows, are mentioned in P. Olver's talk [23]. There are also many other wellknown and famous multi-Hamiltonian systems (see, e.g., [1,2,3,11]). The work on symmetries and conservation laws following from these Hamiltonian properties is now in progress.