In this paper we prove symmetry of nonnegative solutions of the integral equation \(\begin{aligned} u (\zeta ) = \int \limits _{{\mathbb {H}}^n} |\zeta ^{-1} \xi |^{-(Q-\alpha )} u(\xi )^{p} d\xi \quad 1< p \le \frac{Q+\alpha }{Q-\alpha },\ 0< \alpha <Q \end{aligned}\) on the Heisenberg group \({\mathbb {H}}^n = {\mathbb {C}}^n \times {\mathbb {R}},\) \(Q= 2n +2\) using the moving plane method and the Hardy–Littlewood–Sobolev inequality proved by Frank and Lieb for the Heisenberg group. For p subcritical, i.e., \(1< p < \frac{Q+\alpha }{Q-\alpha }\) we show nonexistence of positive solution of this integral equation, while for the critical case, \(p = \frac{Q+\alpha }{Q-\alpha }\) we prove that the solutions are cylindrical and are unique up to a Heisenberg translation and suitable scaling of the function \(\begin{aligned} u_0 (z,t) = \left( (1+ |z|^2)^2 + t^2 \right) ^{- \frac{Q-\alpha }{4}}, \quad (z,t ) \in {\mathbb {H}}^n. \end{aligned}\) As a consequence, we also obtain the symmetry and classification of nonnegative \(C^2\) solutions of the equation \(\begin{aligned} \Delta _{\mathbb {H}}u + u^{p} = 0 \quad \text {for } 1< p \le \frac{Q+\alpha }{Q-\alpha } \text { in } {\mathbb {H}}^n \end{aligned}\) without any partial symmetry assumption on the function u.
