  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.
