The Helmholtz wave equation is linearized using the Feshbach–Villars procedure used for linearizing the Klein–Gordon equation, based on the close algebraic analogy between the Helmholtz equation and the Klein–Gordon equation for a spin-0 particle. The Foldy–Wouthuysen iterative diagonalization technique is then applied to the linearized Helmholtz equation to obtain a Hamiltonian description for a system with varying refractive index. The Hamiltonian has a wavelength-dependent part absent in the traditional descriptions. Besides reproducing all the traditional quasi-paraxial terms, our method leads to additional contributions dependent on the wavelength. Applied to the axially symmetric graded-index fiber, this method results in wavelength-dependent modifications of the paraxial behavior and the aberration coefficients to all orders. Explicit expression for the modified aberration coefficients to the …