So the key point to understanding this problem is to understand that it is the modes that contain information about the physical parameters of your photons (such as the momentum or angular momentum), and quantization is just a description of excitation of those modes.
For instance the canonical quantization of the plane-wave expansion which you've referenced starts off by expanding the vector potential $\mathbf A$ (in the Coulomb gauge) as
$$\mathbf A(\mathbf r) \propto \sum_{\lambda, \mathbf k} \epsilon_\lambda a_\lambda(\mathbf k) e^{i(\mathbf{k\cdot r}-\omega t)} + c.c.
\rightarrow
\hat{\mathbf{A}}(\mathbf r) \propto \sum_{\lambda, \mathbf k} \epsilon_\lambda \hat a_\lambda(\mathbf k) e^{i(\mathbf{k\cdot r}-\omega t)} + h.c.$$
So now each mode designated by the polarization $\lambda$ and the wavevector $\mathbf k$ represents an eigenmode of linear momentum i.e. this mode has momentum $\mathbf p = \hbar\mathbf k$ per photon, and the number of excitations are the only quantum parts i.e. the (average) number of photons in this mode is
$$N_{\lambda,\mathbf k}=\left<\hat a^\dagger_\lambda(\mathbf k) \hat a_\lambda(\mathbf k)\right>.$$
In order to find projections on orbital angular momentum (OAM) eigenstates, or any other type of eigenstate, you need a new representation which means a new modal expansion of $\hat{\mathbf{A}}$. This also means you'll have different operators $\hat a$ that in general can be written as a linear combination of $\hat a_\lambda(\mathbf k)$. For instance if you wanted to treat polarization in the H/V basis (where $\epsilon_{H/V}\propto\epsilon_{\lambda}\pm i\epsilon_{-\lambda}$) then
$$\hat a'_{H/V}\propto\hat a_\lambda\pm i\hat a_{-\lambda}.$$
For a paraxial beam around $\mathbf k = \mathbf k_0$ with transverse wavevector $\mathbf{k_\perp}=(\rho,\phi)$ an OAM mode will have the form $e^{i\ell\phi}$[1], so
$$\hat a(\mathbf{k_\perp,k_0})=\sum_{\ell,p} e^{i\ell\phi}f_{\ell,p}(\rho)\hat a_{\ell,p},$$
and where $f_{\ell,p}(\rho)$ is the radial modes (as an example one could use Laguerre Gaussian modes to have a complete set of radial modes or depending on what you want to compute just remember to sum/average/etc over the $\rho$ component in the end)). Inverting gives
$$\hat a_{\ell,p} = \int \text d\mathbf{k_\perp} e^{-i\ell\phi}f^*_{\ell,p}(\rho)\hat a(\mathbf{k_\perp,k_0}).$$
For some of the details on QFT of OAM see [2]. For nonparaxial beams $\mathbf L$ and $\mathbf S$ don't commute i.e. are coupled which is discussed in [2]. A more complete discussion of this problem is given in [3].
Sources
[1] Allen, Les, et al. "Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes." PRA 45
[2] Calvo, G. F., et al. "Quantum field theory of photons with orbital angular momentum." PRA 73 (2006)
[3] Bliokh, Konstantin Y., et al. "Angular momenta and spin-orbit interaction of nonparaxial light in free space." PRA 82 (2010) arXiv:1006.3876
