  The bound states around a vortex in anisotropic superconductors is a
longstanding yet important issue. In this work, we develop a variational theory
on the basis of the Andreev approximation to obtain the energy levels and wave
functions of the low-energy quantized bound states in superconductors with
anisotropic pairing on arbitrary Fermi surface. In the case of circular Fermi
surface, the effective Schr\"odinger equation yielding the bound state energies
gets back to the theory proposed by Volovik and Kopnin many years ago. Our
generalization here enables us to prove the equidistant energy spectrum inside
a vortex in a broader class of superconductors. More importantly, we are now
able to obtain the wave functions of these bound states by projecting the
quasiclassical wave function on the eigenmodes of the effective Schr\"odinger
equation, going beyond the quasiclassical Eilenberger results, which, as we
find, are sensitive to the scattering rate. For the case of isotropic Fermi
surface, the spatial profile of the low-energy local density of states is
dominated near the vortex center and elongates along the gap antinode
directions, in addition to the ubiquitous Friedel oscillation arising from the
quantum inteference neglected in the Eilenberger theory. Moreover, as a
consequence of the pairing anisotropy, the quantized wave functions develop a
peculiar distribution of winding number, which reduces stepwise towards the
vortex center. Our work provides a flexible way to study the vortex bound
states in the future.
