This paper develops a Riemannian optimization approach for computing p dominant eigenpairs of an \(n\times n\) Hermitian quaternion matrix. We introduce a trace maximization model with weights on the quaternion Stiefel manifold \({{\,\textrm{St}\,}}_{\mathbb {Q}}(n,p)\) and prove that its local solutions are dominant eigenvectors or dominant invariant subspaces. A Riemannian Barzilai-Borwein method on \({{\,\textrm{St}\,}}_{\mathbb {Q}}(n,p)\) with a modified nonmonotone line search is developed, with global convergence established. Applying the algorithm to the trace optimization problem, we obtain a scalable quaternion eigensolver. Key to its scalability is the efficient implementation of retractions on the quaternion Stiefel manifold via reduction to linear systems with complex coefficients, which can be solved by state-of-the-art scientific computing libraries. Numerical experiments on synthetic and real data demonstrate the effectiveness of the proposed RBB method for dominant quaternion eigenpair computations.