  Background: On the mean-field level, pairing correlations are incorporated
through the Bogoliubov-Valatin transformation, whereupon the particle degrees
of freedom are replaced by quasiparticles. This approach leads to a spontaneous
breaking of the particle-number symmetry and mixing of states with different
particle numbers. In order to restore the particle number, various methods have
been employed, which are based on projection approaches before or after
variation. Approximate variation-after-projection (VAP) schemes, utilizing the
Lipkin method, have mostly been used within the Lipkin-Nogami prescription.
  Purpose: Without recurring to the Lipkin-Nogami prescription, and using
instead states rotated in the gauge space, we derive the Lipkin method of
particle-number restoration up to sixth order and we test the convergence and
accuracy of the obtained expansion.
  Methods: We perform self-consistent calculations using the higher-order
Lipkin method to restore the particle-number symmetry in the framework of
superfluid nuclear energy-density functional theory. We also apply the Lipkin
method to a schematic exactly solvable two-level pairing model.
  Results: Calculations performed in open-shell tin and lead isotopes show that
the Lipkin method converges at fourth order and satisfactorily reproduces the
VAP ground-state energies and energy kernels. Near closed shells, the
higher-order Lipkin method cannot be applied because of a non-analytic kink in
the ground-state energies in function of the particle number.
  Conclusions: In open-shell nuclei, the higher-order Lipkin method provides a
good approximation to the exact VAP energies. The method is computationally
inexpensive, making it particularly suitable, for example, for future
optimizations of the nuclear energy-density functionals and simultaneous
restoration of different symmetries.
