We consider the model of a quantum harmonic oscillator governed by a Lindblad master equation where the typical drive and loss channels are multi-photon processes instead of single-photon ones; this implies a dissipation operator of order 2k with integer k>1 for a k-photon process. We prove that the corresponding PDE makes the state converge, for large time, to an invariant subspace spanned by a set of k selected basis vectors; the latter physically correspond to so-called coherent states with the same amplitude and uniformly distributed phases. We also show that this convergence features a finite set of bounded invariant functionals of the state (physical observables), such that the final state in the invariant subspace can be directly predicted from the initial state. The proof includes the full arguments towards the well-posedness of the corresponding dynamics in proper Banach spaces of Hermitian trace-class operators equipped with adapted nuclear norms. It relies on the Hille-Yosida theorem and Lyapunov convergence analysis.