Quantum commutation relationship for photonic orbital angular momentum

Orbital Angular Momentum (OAM) of photons are ubiquitously used for numerous applications. However, there is a fundamental question whether photonic OAM operators satisfy standard quantum mechanical commutation relationship or not; this also poses a serious concern on the interpretation of an optical vortex as a fundamental quantum degree of freedom. Here, we examined canonical angular momentum operators defined in cylindrical coordinates, and applied them to Laguerre-Gauss (LG) modes in a graded index (GRIN) fibre. We confirmed the validity of commutation relationship for the LG modes and found that ladder operators also work properly with the increment or the decrement in units of the Dirac constant (ℏ). With those operators, we calculated the quantum-mechanical expectation value of the magnitude of angular momentum, which includes contributions from both intrinsic and extrinsic OAM. The obtained results suggest that OAM characterised by the LG modes exhibits a well-defined quantum degree of freedom.

1 Introduction

Quantum commutation relationship between operators is an indispensable characteristic in connection with measurements of physical observables [1–4]. Angular momentum operators are especially important as generators of rotation for states with angular momentum state, |m⟩, which is characterised with a quantised integer or a half-integer, m, along a certain direction (say, z) [1–4]. The states pointing directions such as x- and y-directions or any other directions in three-dimensional (3D) space are described by superposition states of orthogonal basis states [1–4], and commutation relationship is used to rotate the states. The spin, $\widehat{\mathbf{S}}$, of one-half in units of the Dirac constant, ℏ = h/(2π), where h is the Plank constant, for an electron is perfectly described with Pauli’s spin matrices [1–4]. The Orbital Angular Momentum (OAM), $\widehat{\mathbf{L}}$, for an electron trapped in spherical potential, e.g., an electron in an atom, is also characterised with an integer quantum number in units of ℏ, where an orbital, such as s, p, d⋯, described by spherical harmonics, $Y_l^m$, satisfies ${\widehat{\mathbf{L}}}^2|l,m\rangle =l(l+1){\hslash }^2|l,m\rangle$, with an integer l associated with the magnitude of angular momentum [1–4]. Furthermore, the total angular momentum, $\widehat{\mathbf{J}}=\widehat{\mathbf{S}}+\widehat{\mathbf{L}}$, of their simple sum is also known as a well-defined quantum angular momentum operator. In general, those angular momentum operators are described by elegant mathematics—Lie algebra, which is looked upon as a triumph in the mathematical formulation of quantum mechanics [1–4].

However, this elegant theoretical framework of angular momentum is less trivial in applying to a photon [5–13], since a photon is usually uni-directionally propagating (say, along z). Therefore, the apparent spherical rotational symmetry is absent for a photon, such that the situation is rather different from that of an electron confined in a spherically symmetric 3D space. The propagation direction of a photon sets a natural quantisation axis; along the propagation direction, the angular momentum operator, ${\widehat{L}}_z$, is successfully obtained by the classical analogue using the Poynting vector [5–8, 11–14]. However, angular momentum operators along the directions perpendicular to the propagation direction (x-y plane) are not defined uniquely in a gauge-invariant way [6–8]. It was argued that the spin and OAM operators would commute [6]. It is now generally believed that the total angular momentum operator $\widehat{\mathbf{J}}$ for photons is well-defined, while it is impossible to split into $\widehat{\mathbf{S}}$ and $\widehat{\mathbf{L}}$ in a local gauge-invariant way [7, 8, 15], despite a big challenge against this claim [16]. More recently, it was demonstrated that spin and orbital angular momentum were successfully separated by considering the SO(3) symmetry [17]. This is outside the scope of this paper, and we will revisit this issue in a forthcoming paper.

The nature of structured light is attracting significant attention, these days [18–23]. The increased bandwidth in fibre optic communication would be one of the most promising near-term application [20, 24]. Another attractive application will be for quantum technologies, where classical entanglement with orthogonal spin and orbital angular momentum states are correlated over variable space and time [18–23, 25–30].

Nevertheless, there are several naive questions, which should be addressed: 1) The spin of a photon could be proper angular momentum, which would satisfy commutation relationship, at least in the absence of OAM. The selection rules of absorption and excitation of a photon in materials [2–4, 9, 31] are clear evidence to expect that the spin of a photon is transferred to angular momenta of an electron. If the spin of a photon is conserved, while accepting well-defined OAM of an electron, why is it regarded as a classical degree of freedom, which is described by a commutable operator? It should be treated with an appropriate quantum commutation relationship, if spin of a photon is a proper quantum mechanical degree of freedom. Moreover, what about the relationship between spin of a photon and the polarisation [32, 33]? Detailed discussion about spin of photons will be provided in a separate paper [28]. 2) A photon with OAM carries quantised angular momentum of ℏm along the direction of propagation, which was successfully described by a Laguerre-Gauss (LG) mode in cylindrical coordinates [5]. Here, some questions arise: Why ℏ appeared in OAM, which is usually the evidence of quantisation? If the OAM is a classical degree of freedom, described by a commutable operator, we normally expect that ℏ would not appear, which is clearly not the case. Can we define standard quantum mechanical canonical angular momentum operators and apply them to LG modes? Moreover, what happens if we define ladder operators for raising or lowering angular momentum in a standard way, such as ${\widehat{L}}_{\pm }={\widehat{L}}_x\pm i{\widehat{L}}_y$, and apply them to the LG modes? Can we show the increment or the decrement of quantised angular momentum in units of ℏ? These are non-trivial questions. In this paper, we will answer to those questions on OAM by directly calculating the matrix elements.

Here is the outline of this paper: In Section 2, we derive fundamental principles and equations, which we are relying on, and explain our model. We are interested in photonics that can be applied to communication technologies and low-energy condensed-matter physics. Consequently, we will not deal with high-energy physics nor those issues related to Lorentz invariance [7, 8, 15, 16] in this paper. We are considering monochromatic coherent light from lasers, such that the incoherent unpolarised light will not be considered, either. The incoherent unpolarised light corresponds to a beam that the radius of the Stokes parameters (S₁, S₂, S₃) does not coincide with its intensity, $S_0\ne \sqrt{S_1^2+S_2^2+S_3^2}$, and thus, it is a typically mixture of beams with a misalignment, such as radiation from the Sun [9]. Here, we consider a coherent beam with properly aligned along the direction of the propagation. We have also included detailed appendices to define associate Laguerre functions and have shown various mathematical formulas to make this paper self-contained. Some of the formulas in appendices are newly derived, particularly for evaluating the ladder operations. We hope the appendices will help readers, since definitions depend on literatures for factors and signs.

In Section 3, we explain our methods to evaluate OAM operators. We clarify the challenges to apply OAM operators to plane-waves, which are not successful. Nevertheless, this would help to understand the problem, which we would like to address. There, we show that the main problem of the plane-waves is a lack of a node at the core of the waveguide, which is also called as topological charge. We show that this problem is solved by using LG modes. In Section 4, we show our main calculation results of various matrix elements for OAM and discuss their implications. Our results show that the LG modes actually satisfy the quantum canonical commutation relationship of angular momentum. We also obtained expectation values of the magnitude of OAM. Finally, in Section 5, we conclude that OAM is indeed a genuine quantum-mechanical observable at least in a graded index (GRIN) fibre [9, 34] satisfying some conditions.

2 Principles and models

2.1 Maxwell’s equations

We start from Maxwell’s equations [9, 10],

$\nabla \times \mathcal{E}=-\frac{\partial \mathcal{B}}{\partial t}(1)$

$\nabla \times \mathcal{H}=\frac{\partial \mathcal{D}}{\partial t}(2)$

$\nabla \cdot \mathcal{D}=0(3)$

$\nabla \cdot \mathcal{B}=0,(4)$

in a non-magnetic transparent material of the dielectric constant ϵ and the permittivity μ₀ without charges and currents. As usual for describing electromagnetic fields [9, 10], we use complex oscillation fields in SI units for electric field $\mathcal{E}$ (V/m), displacement $\mathcal{D}$ (C/m²), magnetic field $\mathcal{H}$ (A/m), and induction $\mathcal{B}$ (Wb/m²), with materials equations $\mathcal{D}=ϵ\mathcal{E}$ and $\mathcal{B}={\mu }_0\mathcal{H}$. All physical observables must be expected to be real without the imaginary part [3] such that experimentally observable fields should be considered by taking the real part, such as for electric field $\mathit{E}=(\mathcal{E}{+\mathcal{E}}^{*})/2$, after the calculation using a complex field of $\mathcal{E}$. While μ₀ is approximately the same as that in vacuum, ϵ is different in a material from the value of ϵ₀ in vacuum [9, 10]. In a non-uniform material, ϵ = ϵ(r) depends on a position r = (x, y, z). If we assume that the profile of ϵ is sufficiently uniform (∇ϵ ≃ 0) compared with the size of wavelength, λ, of a photon, we obtain Helmholtz equation,

${\nabla }^2\mathcal{E}={\mu }_0ϵ\frac{{\partial }^2}{\partial t^2}\mathcal{E},(5)$

which is valid in a completely uniform material and in vacuum. In this paper, we will examine the Helmholtz equation in more detail as follows, but we will not examine its validity any further. The only source of deviations would be arising from the case with significant non-uniformity (∇ϵ ≃ 0). Therefore our analysis will not be valid if ϵ is significantly changed in nano-metre-scale such as for photonic crystals [35–38] and other inhomogeneous systems [39, 40].

2.2 Mapping to Schrödinger equation

First, we see the qualitative feature of the Helmholtz equation in a uniform material by assuming a solution for a linearly polarised monochromatic plane wave

$\mathcal{E}\left(x,y,z\right)=E_0\psi \left(x,y,z\right){\mathrm{e}}^{i\left(kz-\omega t\right)}\widehat{\mathbf{n}},(6)$

where E₀ is the magnitude of the electric field, k is the wavenumber in the material, ω is the angular frequency, the unit vector $\widehat{\mathbf{n}}$ is the direction of the polarisation in the (x, y) plane due to the transversality of the electromagnetic wave, and ψ(x, y, z) is the envelope wavefunction of a photon. Inserting $\mathcal{E}$ into the Helmholtz equation, we obtain

$\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}+2ik\frac{\partial }{\partial z}-k^2\right)\psi =-{\mu }_0ϵ{\omega }^2\psi .(7)$

For various practical applications in laser optics, rays from laser sources are sufficiently collimated, such that the rays can be regarded as paraxial beams [9]. In such a case, we can use slowly varying approximation [9] to neglect the second derivative,

$\frac{{\partial }^2\psi }{\partial z^2}\ll k\frac{\partial \psi }{\partial z},k^2\psi ,(8)$

and obtain [41]

$i\overline{\lambda }\frac{\partial }{\partial z}\psi =-\frac{{\overline{\lambda }}^2}{2n}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)\psi ,(9)$

where we have used the dispersion relationship, ω = vk, with the velocity $v=1/\sqrt{{\mu }_0ϵ}=c/n$ of a photon in a material of a refractive index of n. The velocity of a photon in vacuum is $c=1/\sqrt{{\mu }_0{ϵ}_0}$ and $\overline{\lambda }$ = λ/(2π) is an angular wavelength. Eq. 9 is exactly the same form with a standard non-relativistic Schrödinger equation [2, 3, 31, 41, 42].

$i\hslash \frac{\partial }{\partial t}\psi =-\frac{{\hslash }^2}{2m}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)\psi (10)$

for a particle of mass, m, in a 2D xy-plane at time t. The correspondence is summarised in Table 1.

TABLE 1

TABLE 1. Mapping of Helmholtz equation to a non-relativistic Schrödinger equation.

This implies the propagation of a photon along a paraxial optical path follows the same equation, which describes the dynamics of a massive quantum-mechanical particle [41]. In fact, electron vortices similar to photonic ones were observed [43], and essentially the same mathematical and physical techniques were applicable to both electronic and photonic systems. It is also intuitive to recognise n for a photon corresponds to m of the particle, such that v is low in a material with large n similar to the low velocity of a heavy particle.

2.3 Coherent state for photons

On the contrary to the above similarity between electrons and photons, the important difference is coming from the nature of statistics between Fermions and Bosons. For photons, we are considering monochromatic rays from lasers, which are considered to be in a macroscopic coherent state, exhibiting Bose-Einstein condensation, enabled by the Bose statistics for spin integer particles [1, 4, 11, 44]. On the other hand, electrons are Fermions due to their spin 1/2 characteristics [1, 4], such that a macroscopic coherence is not expected, except for ordered states, such as a superconducting state, which is similar to Bose-Einstein condensation of Cooper pairs [45]. Here, we consider a coherent state for photons to understand the quantum-mechanical state with certain polarisation and orbital angular momentum. A coherent state cannot be described by a fixed number state only due to the phase coherence. Instead, a coherent state is described by a fixed phase, while allowing fluctuation in the number of photons from their average value by a superposition of states with different number of states [11, 44]. Specifically, a coherent state for σ-polarisation in a uniform material is described by

$\begin{array}{ccc}\hfill |{\alpha }_{\sigma }\rangle & =\hfill & {\mathrm{e}}^{-\frac{|{\alpha }_{\sigma }{|}^2}{2}}{\mathrm{e}}^{{\alpha }_{\sigma }{\widehat{a}}_{\sigma }^{†}}|0\rangle \hfill \\ \hfill & =\hfill & {\mathrm{e}}^{-\frac{|{\alpha }_{\sigma }{|}^2}{2}}{\displaystyle \sum _{n_{\sigma }=0}^{\infty }}\frac{{\left({\alpha }_{\sigma }{\widehat{a}}_{\sigma }^{†}\right)}^{n_{\sigma }}}{n_{\sigma }!}|n_{\sigma }\rangle ,\hfill \end{array}(11)$

where σ = H for horizontally polarised state and σ = V for vertically polarised state. α_σ is a complex number, accounting for the macroscopic wavefunction [28]. We use horizontally or vertically polarised states as basis states, for simplicity, but in general we can take other orthonormal bases such as left/right polarised states and diagonal/anti-diagonal states [9]. ${\widehat{a}}_{\sigma }^{†}$ and ${\widehat{a}}_{\sigma }$ are creation and annihilation operators, satisfying Bose commutation relationship [4, 11, 44]

$\left[{\widehat{a}}_{\sigma },{\widehat{a}}_{{\sigma }^{\prime }}^{†}\right]={\delta }_{\sigma ,{\sigma }^{\prime }},(12)$

where δ is the Kronecker delta. A coherent state is best characterised by the fact that it is an eigenstate of an annihilation operator, which can be directly confirmed by the commutation relationship as

$\begin{array}{ccc}\hfill {\widehat{a}}_{\sigma }|{\alpha }_{\sigma }\rangle & =\hfill & {\mathrm{e}}^{-\frac{|{\alpha }_{\sigma }{|}^2}{2}}{\displaystyle \sum _{n_{\sigma }=1}^{\infty }}\frac{{\alpha }^{n_{\sigma }}}{\sqrt{n_{\sigma }!}}\sqrt{n_{\sigma }}|n_{\sigma }-1\rangle \hfill \\ \hfill & =\hfill & {\alpha }_{\sigma }{\mathrm{e}}^{-\frac{|{\alpha }_{\sigma }{|}^2}{2}}{\displaystyle \sum _{n_{\sigma }=0}^{\infty }}\frac{{\alpha }^{n_{\sigma }}}{\sqrt{n_{\sigma }!}}|n_{\sigma }\rangle \hfill \\ \hfill & =\hfill & {\alpha }_{\sigma }|{\alpha }_{\sigma }\rangle .\hfill \end{array}(13)$

We can also confirm that |α_σ⟩ is normalised as ⟨α_σ|α_σ⟩ = 1. Then, we can calculate the average number of photons for both polarised states as

$\langle {\widehat{a}}_{\text{H}}^{†}{\widehat{a}}_{\text{H}}\rangle =|{\alpha }_{\text{H}}{|}^2=N_{\text{H}}=N{\cos }^2\alpha (14)$

$\langle {\widehat{a}}_{\text{V}}^{†}{\widehat{a}}_{\text{V}}\rangle =|{\alpha }_{\text{V}}{|}^2=N_{\text{V}}=N{\sin }^2\alpha ,(15)$

where N_H and N_V are the average number of photons for horizontally and vertically polarisation, respectively, N = N_H + N_V is the total number of photons, and α is the auxiliary angle $(0\le \alpha \le \frac{\pi }{2})$ to describe the polarisation. Alternatively, α_σ is determined by the polarisation state as

${\alpha }_{\text{H}}=\sqrt{N}\cos \alpha (16)$

${\alpha }_{\text{V}}=\sqrt{N}\sin \alpha {\mathrm{e}}^{i\delta },(17)$

where δ ∈ (0, 2π) is the phase of the polarisation. The overall coherent state is described by the direct product as

$|{\alpha }_{\text{H}},{\alpha }_{\text{V}}\rangle =|{\alpha }_{\text{H}}\rangle |{\alpha }_{\text{V}}\rangle .(18)$

Now, we have prepared the coherent state, and the next step is to consider the quantum many-body description of the electromagnetic field, which is achieved by considering the following complex electric field operator

$\widehat{\mathcal{E}}\left(z,t\right)=\sqrt{\frac{2\hslash \omega }{ϵV}}{\mathrm{e}}^{i\beta }\left({\widehat{a}}_{\text{H}}\widehat{\mathbf{x}}+{\widehat{a}}_{\text{V}}\widehat{\mathbf{y}}\right),(19)$

where V is the volume, $\widehat{\mathbf{x}}$ and $\widehat{\mathbf{y}}$ are unit vectors along x and y, and β = kz − ωt + β₀ describes the trivial time and space evolution with a global U (1) phase of β₀. If we apply $\widehat{\mathcal{E}}(z,t)$ to the coherent state |α_H, α_V⟩ as,

$\widehat{\mathcal{E}}\left(z,t\right)|{\alpha }_{\text{H}},{\alpha }_{\text{V}}\rangle =\mathcal{E}\left(z,t\right)|{\alpha }_{\text{H}},{\alpha }_{\text{V}}\rangle ,(20)$

we realise the state is an eigenstate of $\widehat{\mathcal{E}}(z,t)$ with the complex eigenvalue of

$\mathcal{E}\left(z,t\right)=E_0{\mathrm{e}}^{i\beta }\left(\cos \alpha \widehat{\mathbf{x}}+{\mathrm{e}}^{i\delta }\sin \alpha \widehat{\mathbf{y}}\right),(21)$

where $E_0=\sqrt{2\hslash \omega N/(ϵV)}$. This means that the coherent state is a simultaneous eigenstate for diagonalising $\widehat{\mathcal{E}}(z,t)$ and ${\widehat{a}}_{\sigma }$. Alternatively, we can consider the coherent state, $\widehat{\mathcal{E}}(z,t)|{\alpha }_{\text{H}},{\alpha }_{\text{V}}\rangle$, describes the photonic state of the system, since the multiplication of the operator does not change the state. In fact, $\mathcal{E}(z,t)$ actually corresponds to the spinor description of the polarisation state as

$\left(\begin{array}{c}\hfill {\mathcal{E}}_x\hfill \\ \hfill {\mathcal{E}}_y\hfill \end{array}\right)=E_0{\mathrm{e}}^{i\beta }\left(\begin{array}{c}\hfill \cos \alpha \hfill \\ \hfill \sin \alpha {\mathrm{e}}^{i\delta }\hfill \end{array}\right),(22)$

where the matrix part is nothing but a Jones vector [9] to describe the polarisation state of a photon. Moreover, the overall factor of the complex electric field of E₀e^iβ corresponds to the orbital part of the wavefunction for photons in a uniform material, which is the solution of Helmholtz equation.

More generally, for describing a coherent monochromatic ray from a laser source propagating in a waveguide or a fibre, the complex electric field operator must be defined as

$\widehat{\mathcal{E}}\left(\mathbf{r},t\right)=\mathcal{E}\left(\mathbf{r},t\right)\left({\widehat{a}}_{\text{H}}\widehat{\mathbf{x}}+{\widehat{a}}_{\text{V}}\widehat{\mathbf{y}}\right),(23)$

where the scalar complex electric field, $\mathcal{E}(\mathbf{r},t)=E_0\psi (\mathbf{r}){\mathrm{e}}^{i\beta }$, describes the orbital part. The state, $\widehat{\mathcal{E}}(\mathbf{r},t)|{\alpha }_{\text{H}},{\alpha }_{\text{V}}\rangle$, describes the entire photonic state, including polarisation. In this case, $\mathcal{E}(\mathbf{r},t)$ becomes

$\left(\begin{array}{c}\hfill {\mathcal{E}}_x\hfill \\ \hfill {\mathcal{E}}_y\hfill \end{array}\right)=E_0\mathrm{\Psi }\left(\mathbf{r}\right){\mathrm{e}}^{i\beta }\left(\begin{array}{c}\hfill \cos \alpha \hfill \\ \hfill \sin \alpha {\mathrm{e}}^{i\delta }\hfill \end{array}\right),(24)$

where the orbital part is described by Ψ(r) = ψ(r)e^iβ. Ψ(r) is determined by the scalar Helmholtz equation

${\nabla }^2\mathrm{\Psi }\left(\mathbf{r}\right)={\mu }_0ϵ\left(\mathbf{r}\right)\frac{{\partial }^2}{\partial t^2}\mathrm{\Psi }\left(\mathbf{r}\right),(25)$

and thus Ψ(r) is essentially a single-particle wavefunction, describing the orbital degree of freedom. The reason why a macroscopic number of photonic state can be described by a single wavefunction comes from the Bose-Einstein condensed character of a superposition state. All photons are occupying a single state with fixed ω, k, δ, and α, while allowing the fluctuation of the number of photons around its average value of N using a coherent state. The polarisation state is also described as a superposition state of two orthogonal polarisation basis states, coming from intrinsic internal degrees of freedom described by a Jones vector.

The coherent state for a laser beam can also be written as

$\begin{array}{ccc}\hfill |N,\alpha ,\delta \rangle & =\hfill & |{\alpha }_{\text{H}},{\alpha }_{\text{V}}\rangle \hfill \\ \hfill & =\hfill & {\mathrm{e}}^{-\frac{|{\alpha }_{\text{H}}{|}^2}{2}}{\mathrm{e}}^{{\alpha }_{\text{H}}{\widehat{a}}_{\text{H}}^{†}}{\mathrm{e}}^{-\frac{|{\alpha }_{\text{V}}{|}^2}{2}}{\mathrm{e}}^{{\alpha }_{\text{V}}{\widehat{a}}_{\text{V}}^{†}}|0\rangle \hfill \\ \hfill & =\hfill & {\mathrm{e}}^{-\frac{N}{2}}{\mathrm{e}}^{\sqrt{N}\left(\cos \alpha {\widehat{a}}_{\text{H}}^{†}+{\mathrm{e}}^{i\delta }\sin \alpha {\widehat{a}}_{\text{V}}^{†}\right)}|0\rangle .(26)\hfill \end{array}$

If we want to calculate the real electric field, instead of the complex field, we should use the electric field operator defined by

$\widehat{\mathbf{E}}=\frac{1}{2}\left(\widehat{\mathcal{E}}+{\widehat{\mathcal{E}}}^{†}\right),(27)$

which is an observable, such that we can calculate the expectation value, $⟨\widehat{\mathbf{E}}⟩$, quantum-mechanically, using |N, α, δ⟩.

2.4 Laguerre-Gauss mode in a uniform material

Above formalism is based on Maxwell equations, quantum statistics, and superposition principle. Therefore, it is virtually an exact consequence that Ψ(r) represents the wavefunction of coherent photonic states and satisfies the Helmholtz equation at least in a uniform material. The similarity of quantum-mechanical nature of Ψ(r) was intuitively suggested in many pioneering works [5, 6, 8, 14]. Now, it becomes clearer that the intuitive correlation is not merely a coincidence but firmly supported by a quantum many-body theory rather than classical Maxwell’s equations alone, since we cannot derive a wavefunction from classical mechanics. Our formalism contains a standard vacuum state of Quantum Electro-Dynamics (QED) theory [4] in the limit of n → 1, where Ψ(r) will become a simple plane wave, Ψ(r) → e^iβ.

However, the plane wave is not the only solution of the Helmholtz equation, since a solution of differential equation depends also on the symmetry and boundary condition of the system [5, 9]. This is especially true in condensed matter physics, because a material is usually patterned in a specific form with a certain symmetry. Here, we derive a LG mode solution in a uniform material in cylindrical coordinates by using the slowly varying approximation [5, 9]. For completeness, we will describe its full detail in this subsection.

Our starting point is the Helmholtz equation in cylindrical coordinates (r, ϕ)

$i\frac{\partial }{\partial z}\psi =-\frac{1}{2k}\left({\partial }_r^2+\frac{1}{r}{\partial }_r+\frac{1}{r^2}{\partial }_{\varphi }^2\right)\psi ,(28)$

where $r=\sqrt{x^2+y^2}$ is the radius and ϕ = tan⁻¹ (y/x) is the azimuthal angle. We will solve this equation by using a trial wavefunction

$\psi \left(r,\varphi ,z\right)={\left(\frac{r}{w\left(z\right)}\right)}^{m}f\left({\left(\frac{r}{w\left(z\right)}\right)}^2\right){\mathrm{e}}^{iP\left(z\right)+ik\frac{r^2}{2q\left(z\right)}+im\varphi +i\theta \left(z\right)},(29)$

where w(z) is the beam-waist size, P(z) is the phase shift for a beam expansion, q(z) is the complex spherical radius, and θ(z) is another phase shift for radial and azimuthal expansions. Here, we tentatively assume m ≥ 0 for simplicity, and yet we relax this condition for all integer values, including negative values, at the end of the calculation. While the trial wavefunction is inserted into the Helmholtz equation, it is useful to note that ∂_x (fg) = g (∂_xf ) + f (∂_x g) = ψ(∂_x f )/f + ψ∂_x g/g holds. Then, we obtain

$\begin{array}{ccc}\hfill & 2ki\left(i{P}^{\prime }+\frac{1}{q}\right)-\frac{k^2{r}^2}{q^2}\left(1-q^{\prime }\right)-2k{\theta }^{\prime }\hfill & \hfill \\ \hfill & +\frac{4}{w^2}{\left(\frac{r}{w}\right)}^2\frac{f^{″}}{f}+\frac{4ki}{q}{\left(\frac{r}{w}\right)}^2\frac{f^{\prime }}{f}+\frac{2\left(2m+1\right)}{w^2}\frac{f^{\prime }}{f}+\frac{2}{w^2}\frac{f^{\prime }}{f}\hfill & \hfill \\ \hfill & -2ki\frac{m{w}^{\prime }}{w}-4ki\frac{r^2{w}^{\prime }}{w^3}\frac{f^{\prime }}{f}+2ki\frac{m}{q}=0.(30)\hfill & \hfill \end{array}$

The Gaussian mode solution is given by assuming

$\frac{\partial P}{\partial z}=\frac{i}{q}(31)$

$\frac{\partial q}{\partial z}=1,(32)$

which will give us q(z) = z + q₀ = z − iz₀ and $P(z)=i\ln \left(1+z/q_0\right)$ with the confocal parameter $z_0=k{w}_0^2/2=\pi n{w}_0^2/\lambda$ and the minimum waist of w₀. We then obtain the Gaussian factor

${\mathrm{e}}^{ik\frac{r^2}{2q\left(z\right)}}={\mathrm{e}}^{-\frac{r^2}{w^2}}{\mathrm{e}}^{i\frac{k{r}^2}{2R}}(33)$

and the phase-shift factor

${\mathrm{e}}^{iP\left(z\right)}=\frac{w_0}{w}{\mathrm{e}}^{-i\eta \left(z\right)},(34)$

where the beam waist, w(z), the beam radius, R(z), and the phase, η(z), are given by

$w\left(z\right)=w_0\sqrt{1+{\left(\frac{z}{z_0}\right)}^2}(35)$

$R\left(z\right)=z+\frac{z_0^2}{z}(36)$

$\eta \left(z\right)={\tan }^{-1}\left(\frac{z}{z_0}\right).(37)$

The focal point of the Gaussian beam is set at the origin, where the waist becomes minimum w (0) = w₀. In a uniform material or a vacuum, there is no mechanism to confine the mode, and the beam waist can be arbitrarily controlled by the use of an optical lens up to the diffraction limit. Thus, w₀, and consequently z₀, can be controlled and determined by a boundary condition. It is also useful to note that kww′ = 2z/z₀ and $k{w}^2/q=2\left(i+z/z_0\right)$ hold, and we obtain

$\begin{array}{ccc}\hfill & \hfill & -k{w}^2{\theta }^{\prime }-2m\hfill \\ \hfill & \hfill & -4{\left(\frac{r}{w}\right)}^2\frac{f^{\prime }}{f}+2\left(m+1\right)\frac{f^{\prime }}{f}+2{\left(\frac{r}{w}\right)}^2\frac{f^{″}}{f}=0.\hfill \end{array}(38)$

We now focus on the last three terms of this equation, which can be rewritten by exchanging valuables subsequently using ρ = r/w, a = ρ², and b = 2a as

$\begin{array}{ccc}\hfill & \hfill & -4{\rho }^2\frac{f^{\prime }}{f}+2\left(m+1\right)\frac{f^{\prime }}{f}+2{\rho }^2\frac{f^{″}}{f}\hfill \\ \hfill & =\hfill & \frac{2}{f}\left[-2a\frac{d}{da}+\left(m+1\right)\frac{d}{da}+a\frac{d^2}{d{a}^2}\right]f\hfill \\ \hfill & =\hfill & \frac{4}{f}\left[b\frac{d^2}{d{b}^2}+\left(m+1-b\right)\frac{d}{db}\right]f=-4p,(39)\hfill \end{array}$

where we used the differential equation for the associate Laguerre function, $L_p^m$ (Section 4),

$\left[b\frac{d^2}{d{b}^2}+\left(m+1-b\right)\frac{d}{db}\right]L_p^m\left(b\right)=-p{L}_p^m\left(b\right),(40)$

such that we can obtain $f(b)=L_p^m(b)=L_p^m\left(2{\left(\frac{r}{w}\right)}^2\right)$. The rest of the Helmholtz equation is

$k{w}^2{\theta }^{\prime }=-2\left(2p+m\right),(41)$

which gives the phase-shift

$\theta =-\left(2p+m\right){\tan }^{-1}\left(z/z_0\right).(42)$

Finally, we obtain

$\begin{array}{ccc}\hfill \psi \left(r,\varphi ,z\right)=& \hfill & \frac{w_0}{w}{\left(\frac{\sqrt{2}r}{w}\right)}^m{L}_p^m\left(2{\left(\frac{r}{w}\right)}^2\right)\hfill \\ \hfill & \hfill & {\mathrm{e}}^{-\frac{r^2}{w^2}}{\mathrm{e}}^{ik\frac{r^2}{2R}}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{-i\left(2p+m+1\right){\tan }^{-1}\left(z/z_0\right)},\hfill \end{array}(43)$

which is not normalised, yet. The norm N_norm of the wavefunction is obtained by

$\begin{array}{ccc}\hfill N_{\mathrm{norm}}^2& =\hfill & {\int }_0^{\infty }dr2\pi r|\psi \left(r,\varphi ,z\right){|}^2\hfill \\ \hfill & =\hfill & \frac{\pi }{2}{w}_0^2{\int }_0^{\infty }db{b}^m{\left|L_p^m\left(b\right)\right|}^2{\mathrm{e}}^{-b}\hfill \\ \hfill & =\hfill & \frac{\pi }{2}{w}_0^2\frac{\left(p+m\right)!}{p!},(44)\hfill \end{array}$

where we used the orthogonality condition (Section 4)

${\int }_0^{\infty }db{\mathrm{e}}^{-b}{b}^m{L}_l^m\left(b\right)L_n^m\left(b\right)=\frac{\left(l+m\right)!}{l!}{\delta }_{l,n}.(45)$

Thus, we obtain

$N_{\mathrm{norm}}=w_0\sqrt{\frac{\pi }{2}{w}_0^2\frac{\left(p+m\right)!}{p!}}.(46)$

Now, we consider the case for a negative value of m. The only source of the azimuthal dependence in the Helmholtz equation is coming from

${\partial }_{\varphi }^2\psi =-\frac{m^2}{r^2}\psi ,(47)$

such that the solution does not depend on the sign of m.

Therefore, the final normalised wavefunction becomes

$\begin{array}{ccc}\hfill \psi \left(r,\varphi ,z\right)=& \sqrt{\frac{2}{\pi }\frac{p!}{\left(p+|m|\right)!}}\frac{1}{w}{\left(\frac{\sqrt{2}r}{w}\right)}^{|m|}{L}_p^{|m|}\left(2{\left(\frac{r}{w}\right)}^2\right)\hfill & \hfill \\ \hfill & {\mathrm{e}}^{-\frac{r^2}{w^2}}{\mathrm{e}}^{ik\frac{r^2}{2R}}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{-i\left(2p+|m|+1\right){\tan }^{-1}\left(z/z_0\right)}.\hfill & \hfill \end{array}(48)$

Here, the wavefunction was normalised in the xy-plane as

${\int }_0^{\infty }dr2\pi r|\psi \left(r,\varphi ,z\right){|}^2=1,(49)$

since our main interests in the following sections are orbital angular momentum, and this normalisation is easier to treat. On the other hand, in the consideration of the electric field in the previous subsection, the normalisation was slightly different, since we have prioritised to have the proper definition of E₀ (V/cm) as the electric field. When the number of photons that we are considering is one, this corresponds to the zero-point fluctuation of the electric field, $e_0=\sqrt{2\hslash \omega /ϵV}$, but actually a laser beam contains a macroscopic number of photons. By comparing the factors between E₀ and ψ(r, ϕ, z), we realise that we should assume V = w(z)²L, where L is the length of the system along z. This means that the magnitude of the electric field changes upon propagation due to the change of the beam waist. If the beam expands, the electric field decreases, and vice versa. This is attributed to the change in size of the mode profile for photons. If we use the normalisation of this subsection, E₀ should be simply re-defined as $E_0=\sqrt{2\hslash \omega N/ϵL}$.

It is worth making a remark on the Gouy phase [5, 42, 46–50] of

${\varphi }_{\text{G}}=\left(2p+|m|+1\right){\tan }^{-1}\left(z/z_0\right),(50)$

which is the same as a geometrical phase of Pancharatnam-Berry. This term appears due to the focusing of the beam, which will change $({\mathcal{E}}_x,{\mathcal{E}}_y)$ at z → −∞ to $(-{\mathcal{E}}_x,-{\mathcal{E}}_y)$ at z → +∞ upon crossing the focal point at z = 0. This change is taken into account within our orbital wavefunction ψ(r, ϕ, z), where the polarisation state is not changed. In the absence of OAM, corresponding to m = 0 and p = 0, this change just accompanies a phase-shift of π. With OAM, the extra phase factor of e^imϕ will contribute to it as an additional phase-shift of πm, because the focusing corresponds to rotating the phase from ϕ = 0 to ϕ = π. This global change of the orbital due to focussing together with the local rotation of the phase by OAM is responsible for the Gouy phase. In addition, the radial distribution due to the mode shape described by a Laguerre function will also contribute, in a similar way. For the radial profile, we have p-nodes along the radial direction, where the focussing corresponds to change the phase front located at r = w(z) to r = −w(z). During this change, the nodes along r will go across the origin, while adding a phase-shift of πp. In addition, there are nodes at −r, or equivalently (r, ϕ = −π), along the opposite radial direction. Therefore, the total contribution to the phase-shift from radial oscillation is 2πp. The actual change of the phase is not abrupt, and it adiabatically changes in the length scale of z₀. For the propagation in a GRIN fibre, this Gouy phase is, fortunately, not so important because of the absence of focusing, as we shall see in the next subsection.

We should be careful for the interpretation of the radial quantum number, p, which describes the number of nodes along the radial direction. This value is different from the quantum number, l, to describe the magnitude of OAM in a spherical symmetric system by $Y_l^m$, for which the value of m is limited to be m = l, l − 1, …, 0, ⋯ − (l − 1), − l. On the other hand, there is no such restriction to $L_p^{|m|}$, where the so-called LG₀₁ mode at p = 0 and m = 1 can be well-defined, for example. This also means that p cannot be the proper quantum number to be assigned as the magnitude of OAM. In fact, the LG mode is not the simultaneous eigenstate of the magnitude of OAM and the component of OAM along the quantised axis (z), albeit the expectation value of the magnitude depending on p. It is reasonable to expect that there exists the simultaneous eigenstates, according to the general theory of OAM, but LG modes do not diagonalise the operator for the magnitude of OAM. Henceforth, we also use n for the radial quantum number later, instead of the popular use of p to avoid unnecessary confusion to the momentum, p = ℏk.

2.5 Laguerre-Gauss mode in a graded index fibre

As another example, for which the Helmholtz equation can be solved exactly, we also discuss the propagation of a coherent monochromatic laser beam in the GRIN fibre [9, 34], which has the refractive index n(r) dependence given by $n{(r)}^2=ϵ(r)/{ϵ}_0=n_0^2\left(1-g^2{r}^2\right)$, where the graded index parameter, g = 2π/Λ, has the dimension of inverse length, and g must be small to justify the derivation of the Helmholtz equation such that the index profile is sufficiently gentle (Λ ≫ λ). We also define the wavenumber in vacuum as k₀ = 2π/λ for a laser beam emitted from the waveguide. The energy of the photon will not be changed upon the emission, such that ω = ck₀ is valid, while the dispersion relationship, ω = ω(k), in the waveguide is highly non-trivial, and we will obtain this from the Helmholtz equation. We also define a constant wavenumber parameter, $k_{n_0}=2\pi /{\lambda }_{n_0}=\frac{2\pi n_0}{\lambda }=k_0{n}_0$, since the waveguide is mostly determined by the core refractive index, n₀ = n (0), and this is simply a constant parameter, where $k_{n_0}$ is different from the true wavenumber, k, responsible for describing photon momentum of p = ℏk.

With those parameters, we can rewrite ${\mu }_0ϵ{\omega }^2=n{(r)}^2/c^2{(c{k}_0)}^2=k_0^{2}n{(r)}^2=k_{n_0}^2(1-g^2{r}^2)$, and the Helmholtz equation becomes

$\left({\nabla }^2+2ik\frac{\partial }{\partial z}+\left(k_{n_0}^2-k^2\right)+k_{n_0}^2{g}^2{r}^2\right)\psi =0.(51)$

in the cylindrical coordinate, we can convert this equation to

$\begin{array}{ccc}\hfill & \hfill & \left({\partial }_r^2+\frac{1}{r}{\partial }_r+2ik\partial z+\partial z^2\right.\hfill \\ \hfill & \hfill & \left.+\frac{1}{r^2}{\partial }_{\varphi }^2+\left(k_{n_0}^2-k^2\right)+k_{n_0}^2{g}^2{r}^2\right)\psi =0\hfill \end{array}(52)$

One of the conceptual advantages to consider the GRIN waveguide is that we do not have to worry about the paraxial slowly varying approximation at all, because the second derivative along z vanishes. This can be verified by confirming that the trial wavefunction

$\psi \left(r,\varphi ,z\right)={\left(\frac{r}{w_0}\right)}^{m}f\left({\left(\frac{r}{w_0}\right)}^2\right){\mathrm{e}}^{i{k}_{n_0}\frac{r^2}{2q}+im\varphi },(53)$

with the constant waist w₀ and the constant complex radius, q, becomes the solution. Again, we tentatively assume m ≥ 0. By inserting Eq. 53 into Eq. 52, we obtain

$\begin{array}{ccc}\hfill & 2k_{n_0}i\frac{1}{q}+k_{n_0}^2{r}^2\left(-\frac{1}{q^2}+g^2\right)+\left(k_{n_0}^2-k^2\right)\hfill & \hfill \\ \hfill & +\frac{4}{w_0^2}{\left(\frac{r}{w_0}\right)}^2\frac{f^{″}}{f}+\frac{4k_{n_0}i}{q}{\left(\frac{r}{w_0}\right)}^2\frac{f^{\prime }}{f}+\frac{2\left(2m+1\right)}{w^2}\frac{f^{\prime }}{f}+\frac{2}{w_0^2}\frac{f^{\prime }}{f}\hfill & \hfill \\ \hfill & +2k_{n_0}i\frac{m}{q}=0.\hfill & \hfill \end{array}(54)$

we then obtain a stable Gaussian form by noting

$q=-i\frac{i}{g}(55)$

$w_0=\sqrt{\frac{2}{g{k}_{n_0}}}.(56)$

we also use useful identities

$2k_{n_0}i\frac{1}{q}=-4\frac{4}{w_0^2}(57)$

$g{w}_0^2=\frac{2}{k_{n_0}},(58)$

and the Helmholtz equation then becomes

$\begin{array}{ccc}\hfill & \hfill & -2\left(m+1\right)+\left(k_{n_0}^2-k^2\right)\frac{w_0^2}{2}\hfill \\ \hfill & +\hfill & \left[2{\left(\frac{r}{w_0}\right)}^2\frac{f^{″}}{f}+2\left(m+1\right)\frac{f^{\prime }}{f}-4{\left(\frac{r}{w_0}\right)}^2\frac{f^{\prime }}{f}\right]=0.\hfill \\ \hfill & \hfill & \hfill \end{array}(59)$

By noticing that the last three terms in the left-hand side of Eq. 59 can be described by the associated Laguerre function, we obtain

$f=L_n^m\left(2{\left(\frac{r}{w_0}\right)}^2\right),(60)$

and the remaining equation becomes

$k^2=k_{n_0}^2-2g{k}_{n_0}\left(2n+m+1\right).(61)$

This provides the solution [9]

$k=\frac{\omega }{v_0}\sqrt{1-2\frac{\delta {\omega }_0}{\omega }\left(2n+m+1\right)},(62)$

where we have defined a phase velocity at the core as v₀ = c/n₀ and a frequency shift as δω₀ = v₀g. By solving Eq. 62 with respect to ω, we obtain the dispersion relationship

$\omega =\sqrt{v_0^2{k}^2+\delta {\omega }_0^2{\left(2n+m+1\right)}^2}+\delta {\omega }_0\left(2n+m+1\right).(63)$

This dispersion relationship can be intuitively understood as follows: Since ω ≠ 0 at k = 0, this implies an opening of an energy gap in a band diagram, meaning that the dispersion is “massive.” The emergence of an energy gap is reminiscent of the theory of superconductivity [45] and the Nambu-Anderson-Goldstone-Higgs theory of a broken symmetry [51–54]. We infer that a similar symmetry principle is hidden in our system. We will discuss this in a subsequent paper [55]. Here, we can recognise the increase of the gap by increasing the radial quantum number, p, and the quantised OAM number, m, because the discrete photon energy is related to the confinement degrees of freedom of photons rather than the free propagation of photons along z.

Finally, we relax the condition for m to allow negative integers without breaking the formalism. Normalising the wavefunction as before, we obtain an exact solution

$\begin{array}{ccc}\hfill \psi \left(r,\varphi ,z\right)=& \hfill & \sqrt{\frac{2}{\pi }\frac{n!}{\left(n+|m|\right)!}}\frac{1}{w_0}{\left(\frac{\sqrt{2}r}{w_0}\right)}^{|m|}\hfill \\ \hfill & \hfill & L_n^{|m|}\left(2{\left(\frac{r}{w_0}\right)}^2\right){\mathrm{e}}^{-\frac{r^2}{w_0^2}}{\mathrm{e}}^{im\varphi },\hfill \end{array}(64)$

in a GRIN fibre without the slowly varying paraxial approximation.

3 Methods

3.1 Orbital angular momentum for photons

We have confirmed the fundamental principle on how to treat a coherent laser beam on the basis of Maxwell’s equations and a quantum many-body theory. In particular, we have understood why we can describe a macroscopically coherent laser by a single-particle wavefunction, Ψ(r), due to Bose-Einstein statistics, while the entire many-body state is described by a coherent state with both orbital and spin degrees of freedom. Photons are quantum-mechanical particles with a wave nature; we can also describe them with Maxwell’s equations together with a quantum many-body theory. For coherent photons from a laser, it was less obvious how we can treat the ray quantum-mechanically; however, a laser produces indistinguishable photons with the same phase by the stimulated emission process, in which existing photons in a cavity induce recombinations of electron-hole pairs to make clones of photons as a results of a chain-reaction. Thus, we can describe a coherent monochromatic ray by the single mode of Ψ(r). If the waveguide contains several modes, it is straightforward to allow the superposition of these macroscopically coherent beams.

The fundamental equation for describing the orbital character of Ψ(r) is the Helmholtz equation, instead of the Schrödinger equation, although we have a significant similarity to a paraxial wave (Table 1). Unlike in vacuum without a material, where Ψ(r) is a simple plane-wave, the mode profile of Ψ(r) can be highly non-trivial in materials, depending on the symmetries of the waveguides and the actual profile of the refractive index, n(r). In the previous section, we have obtained LG modes in a GRIN fibre with a uniform material. Here, the LG modes (LG_nm) were clearly labelled by the radial quantum number n and the quantum optical orbital angular momentum number m along the propagation direction. The central theme of this paper is to examine the validity of this interpretation that m is indeed a proper quantum index to describe the optical OAM, and thus the angular momentum of the orbital is ℏm. Under the assumption that a standard quantum mechanical treatment is also applicable to photons, described by Ψ(r), we examine the impacts of OAM operations in the following sections.

3.2 Canonical orbital angular momentum operators

The most well-established quantum mechanical treatment is canonical commutation relationship for the position $\widehat{\mathbf{r}}=(\widehat{x},\widehat{y},\widehat{z})$ and the momentum $\widehat{\mathbf{p}}=({\widehat{p}}_x,{\widehat{p}}_y,{\widehat{p}}_z)$ operators [1–3]: $[\widehat{x},{\widehat{p}}_x]=i\hslash$, $[\widehat{y},{\widehat{p}}_y]=i\hslash$, $[\widehat{z},{\widehat{p}}_z]=i\hslash$, and commutable relationship among other combinations. The OAM operator, $\widehat{\mathbf{l}}$, is defined by use of those $\widehat{\mathbf{r}}$ and $\widehat{\mathbf{p}}$ as $\widehat{\mathbf{l}}=\widehat{\mathbf{r}}\times \widehat{\mathbf{p}}$, such that each component becomes

${\widehat{l}}_x=\frac{\hslash }{i}\left(y{\partial }_z-z{\partial }_y\right)(65)$

${\widehat{l}}_y=\frac{\hslash }{i}\left(z{\partial }_x-x{\partial }_z\right)(66)$

${\widehat{l}}_z=\frac{\hslash }{i}\left(x{\partial }_y-y{\partial }_x\right),(67)$

respectively. In a system with a spherical symmetry, the eigenstate of these operators are described by Y_lm (θ, ϕ) = ⟨θ, ϕ|l, m⟩, where θ and ϕ are polar and azimuthal angels, respectively, l and m are quantum numbers for the magnitude of OAM and the OAM component along the quantisation axis, respectively [1–3]. Our goal is to obtain a similar relationship for a system with a cylindrical symmetry, described by the LG modes. In this section, we will obtain the operator representation of $\widehat{\mathbf{l}}$ in cylindrical coordinates.

The unit vectors of cylindrical coordinates are defined for a rotation in a 3D Cartesian coordinate along the z-axis,

$\left(\begin{array}{c}\hfill \widehat{\mathbf{r}}\hfill \\ \hfill \widehat{\mathbf{\Phi }}\hfill \end{array}\right)=\left(\begin{array}{cc}\hfill \cos \varphi \hfill & \hfill \sin \varphi \hfill \\ \hfill -\sin \varphi \hfill & \hfill \cos \varphi \hfill \end{array}\right)\left(\begin{array}{c}\hfill \widehat{\mathbf{x}}\hfill \\ \hfill \widehat{\mathbf{y}}\hfill \end{array}\right),(68)$

where $\widehat{\mathbf{r}}$ is the unit vector along r and $\widehat{\mathbf{\Phi }}$ is a unit vector along the azimuthal direction, while z is unchanged. The important point in this coordinate is ϕ dependence of these unit vectors, i.e., $\widehat{\mathbf{r}}=\widehat{\mathbf{r}}(\varphi )$ and $\widehat{\mathbf{\Phi }}=\widehat{\mathbf{\Phi }}(\varphi )$. Therefore, the nabla operator is of the form

$\nabla =\left(\frac{\partial }{\partial r},\frac{1}{r}\frac{\partial }{\varphi },\frac{\partial }{\partial z}\right)=\frac{\partial }{\partial r}\widehat{\mathbf{r}}+\frac{1}{r}\frac{\partial }{\partial \varphi }\widehat{\mathbf{\Phi }}+\frac{\partial }{\partial z}\widehat{\mathbf{z}},(69)$

and the Laplacian becomes

${\nabla }^2=\frac{{\partial }^2}{\partial r^2}+\frac{{\partial }^2}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^2}{\partial {\varphi }^2}+\frac{1}{r}\frac{\partial }{\partial r},(70)$

where we have used

$\frac{\partial }{\partial \varphi }\widehat{\mathbf{r}}\left(\varphi \right)=\widehat{\mathbf{\Phi }}.(71)$

This ϕ dependence of $\widehat{\mathbf{r}}(\varphi )$ and $\widehat{\mathbf{\Phi }}(\varphi )$ makes it difficult to define and treat the OAM operators for ${\widehat{l}}_r$ and ${\widehat{l}}_{\varphi }$. Instead, we will keep using ${\widehat{l}}_x$ and ${\widehat{l}}_y$, defined above, and we will express them by (r, ϕ, z). In order to use (r, ϕ) instead of (x, y), we obtain

${\partial }_x=\cos \varphi {\partial }_r-\frac{1}{r}\sin \varphi {\partial }_{\varphi },(72)$

${\partial }_y=\sin \varphi {\partial }_r+\frac{1}{r}\cos \varphi {\partial }_{\varphi }.(73)$

By inserting these into $\widehat{\mathbf{l}}$, we obtain

${\widehat{l}}_x=\frac{\hslash }{i}\left(-z\sin \varphi {\partial }_r-\frac{z}{r}\cos \varphi {\partial }_{\varphi }+r\sin \varphi {\partial }_z\right),(74)$

${\widehat{l}}_y=\frac{\hslash }{i}\left(z\cos \varphi {\partial }_r-\frac{z}{r}\sin \varphi {\partial }_{\varphi }-r\cos \varphi {\partial }_z\right),(75)$

${\widehat{l}}_z=\frac{\hslash }{i}{\partial }_{\varphi }.(76)$

We can readily confirm the original commutation relationship

$\left[{\widehat{l}}_x,{\widehat{l}}_y\right]=i\hslash {\widehat{l}}_z(77)$

and its cyclic exchanges

$\left[{\widehat{l}}_y,{\widehat{l}}_z\right]=i\hslash {\widehat{l}}_x,(78)$

$\left[{\widehat{l}}_z,{\widehat{l}}_x\right]=i\hslash {\widehat{l}}_y(79)$

are all valid in the cylindrical coordinate (r, ϕ, z).

We also obtain raising and lowering operators, respectively, as

${\widehat{l}}_{+}={\widehat{l}}_x+i{\widehat{l}}_y(80)$

$=\hslash {\mathrm{e}}^{i\varphi }\left(z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }-r{\partial }_z\right),(81)$

${\widehat{l}}_{-}={\widehat{l}}_x-i{\widehat{l}}_y(82)$

$=\hslash {\mathrm{e}}^{-i\varphi }\left(-z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }+r{\partial }_z\right).(83)$

3.3 Application to plane waves and problems

So far, it was straightforward to develop a theory for photonic OAM. In this subsection, we will apply our canonical OAM operator to plane waves to clarify that problems arise. Specifically, we consider a plane wave with OAM in the simplest form:

$\mathrm{\Psi }\left(r,\varphi ,z\right)={\mathrm{e}}^{ikz}{\mathrm{e}}^{im\varphi },(84)$

which is not the solution of the Helmholtz equation at m ≠ 0. Nevertheless, it is useful to clarify the potential issue and to explain what we should address in the following sections.

First, multiplying it by $\widehat{\mathbf{l}}$, we obtain

${\widehat{l}}_{+}\left(r,\varphi ,z\right)\mathrm{\Psi }=\hslash {\mathrm{e}}^{i\varphi }\left(-\frac{z}{r}m-ikr\right)\mathrm{\Psi },(85)$

${\widehat{l}}_{-}\left(r,\varphi ,z\right)\mathrm{\Psi }=\hslash {\mathrm{e}}^{-i\varphi }\left(-\frac{z}{r}m+ikr\right)\mathrm{\Psi },(86)$

${\widehat{l}}_z\left(r,\varphi ,z\right)\mathrm{\Psi }=\hslash m\mathrm{\Psi },(87)$

which means that Ψ is indeed an eigenstate for l_z and that ${\widehat{l}}_{\pm }$ is effectively working to raise and lower the eigenvalue of the angular momentum component along the direction of the propagation. If we multiply them by Ψ* from the left, we obtain

${\mathrm{\Psi }}^{*}{\widehat{l}}_{+}\left(r,\varphi ,z\right)\mathrm{\Psi }=\hslash {\mathrm{e}}^{i\varphi }\left(-\frac{z}{r}m-ikr\right),(88)$

${\mathrm{\Psi }}^{*}{\widehat{l}}_{-}\left(r,\varphi ,z\right)\mathrm{\Psi }=\hslash {\mathrm{e}}^{-i\varphi }\left(-\frac{z}{r}m+ikr\right),(89)$

${\mathrm{\Psi }}^{*}{\widehat{l}}_z\left(r,\varphi ,z\right)\mathrm{\Psi }=\hslash m.(90)$

by averaging these over space, we obtain

$⟨{\widehat{l}}_x\left(r,\varphi ,z\right)⟩=0,(91)$

$⟨{\widehat{l}}_y\left(r,\varphi ,z\right)⟩=0,(92)$

$⟨{\widehat{l}}_z\left(r,\varphi ,z\right)⟩=\hslash m.(93)$

therefore, the expectation values are reasonably well-defined.

However, if we calculate the complex conjugate of ${\widehat{l}}_{-}\mathrm{\Psi }$ simply by taking its complex conjugate as

${\mathrm{\Psi }}^{*}{\left({\widehat{l}}_{-}\left(r,\varphi ,z\right)\right)}^{†}=\hslash {\mathrm{e}}^{-i\left(m-1\right)\varphi }\left(-\frac{z}{r}m-ikr\right){\mathrm{e}}^{-ikz},(94)$

and multiply this by ${\widehat{l}}_{-}\mathrm{\Psi }$ from the right to calculate the norm, we obtain

${\mathrm{\Psi }}^{*}{\left({\widehat{l}}_{-}\left(r,\varphi ,z\right)\right)}^{†}{\widehat{l}}_{-}\left(r,\varphi ,z\right)\mathrm{\Psi }={\hslash }^2\left(\frac{z^2}{r^2}{m}^2+k^2{r}^2\right),(95)$

which is a positive real value. On the other hand, the direct calculation of ${\mathrm{\Psi }}^{*}{\widehat{l}}_{+}{\widehat{l}}_{-}\mathrm{\Psi }$ becomes

${\mathrm{\Psi }}^{*}{\widehat{l}}_{+}\left(r,\varphi ,z\right){\widehat{l}}_{-}\left(r,\varphi ,z\right)\mathrm{\Psi }={\hslash }^2\left(\frac{z^2}{r^2}{m}^2+2ikz+{\left(kr\right)}^2+m\right).(96)$

This implies that

${\left({\widehat{l}}_{-}\left(r,\varphi ,z\right)\right)}^{†}\ne {\widehat{l}}_{+}\left(r,\varphi ,z\right),(97)$

which means that the l_± is not observable for the Hilbert space spanned by the plane waves with OAM. We can also confirm the conjugate relationships

${\mathrm{\Psi }}^{*}{\left({\widehat{l}}_{+}\left(r,\varphi ,z\right)\right)}^{†}{\widehat{l}}_{+}\left(r,\varphi ,z\right)\mathrm{\Psi }={\hslash }^2\left(\frac{z^2}{r^2}{m}^2+k^2{r}^2\right)(98)$

${\mathrm{\Psi }}^{*}{\widehat{l}}_{-}\left(r,\varphi ,z\right){\widehat{l}}_{+}\left(r,\varphi ,z\right)\mathrm{\Psi }={\hslash }^2\left(\frac{z^2}{r^2}{m}^2+2ikz+{\left(kr\right)}^2-m\right),(99)$

which also imply

${\left({\widehat{l}}_{+}\left(r,\varphi ,z\right)\right)}^{†}\ne {\widehat{l}}_{-}\left(r,\varphi ,z\right).(100)$

We also see

$\begin{array}{ccc}\hfill & \hfill & {\mathrm{\Psi }}^{*}{\left({\widehat{l}}_{-}\left(r,\varphi ,z\right)\right)}^{†}{\widehat{l}}_{-}\left(r,\varphi ,z\right)\mathrm{\Psi }\hfill \\ \hfill & \hfill & ={\mathrm{\Psi }}^{*}{\left({\widehat{l}}_{+}\left(r,\varphi ,z\right)\right)}^{†}{\widehat{l}}_{+}\left(r,\varphi ,z\right)\mathrm{\Psi },\hfill \end{array}(101)$

showing a classical result without providing commutation relationship. This is a remarkable difference from the standard quantum mechanics [3], which shows the canonical commutation relationship upon the calculation of the norm for ${\widehat{l}}_{\pm }{Y}_l^m(\theta ,\varphi )$.

On the other hand, the direct calculation shows

$\begin{array}{ccc}\hfill & {\mathrm{\Psi }}^{*}{\widehat{l}}_{+}\left(r,\varphi ,z\right){\widehat{l}}_{-}\left(r,\varphi ,z\right)\mathrm{\Psi }-{\mathrm{\Psi }}^{*}{\widehat{l}}_{-}\left(r,\varphi ,z\right){\widehat{l}}_{+}\left(r,\varphi ,z\right)\mathrm{\Psi }\hfill & \hfill \\ \hfill & =2\hslash {\mathrm{\Psi }}^{*}{\widehat{l}}_z\left(r,\varphi ,z\right)\mathrm{\Psi },\hfill & \hfill \end{array}(102)$

such that the commutation relationship

$\left[{\widehat{l}}_{+}\left(r,\varphi ,z\right),{\widehat{l}}_{-}\left(r,\varphi ,z\right)\right]=2\hslash {\widehat{l}}_z\left(r,\varphi ,z\right)(103)$

is indeed satisfied on the average.

These apparent contradiction and inconsistency are coming from the assumption of the ill-defined plane-wave wavefunction, Ψ(r, ϕ, z) = e^ikze^imϕ. This is confirmed by calculating the magnitude of the OAM along the radial direction as

${\mathrm{\Psi }}^{*}\left({\widehat{l}}_x^2+{\widehat{l}}_y^2\right)\mathrm{\Psi }={\hslash }^2\left(m^2\frac{z^2}{r^2}+k^2{r}^2+2ikz\right),(104)$

which gives an imaginary part, thus showing that the magnitude is not observable. If we take the average over z ∈ (0, L), we obtain

${\int }_0^L\frac{dz}{L}{\mathrm{\Psi }}^{*}\left({\widehat{l}}_x^2+{\widehat{l}}_y^2\right)\mathrm{\Psi }={\hslash }^2\left(m^2\frac{L^2}{3r^2}+k^2{r}^2+ikL\right),(105)$

which is still a complex value. Further average over r ∈ (0, R), where R is the radius of a cylindrical waveguide, gives

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^R\frac{2\pi rdr}{\pi R^2}{\int }_0^L\frac{dz}{L}{\mathrm{\Psi }}^{*}\left({\widehat{l}}_x^2+{\widehat{l}}_y^2\right)\mathrm{\Psi }\hfill \\ \hfill & \hfill & ={\hslash }^2\left(\frac{2m^2{L}^2}{3R^2}\left(\ln R-\ln 0\right)+\frac{1}{2}{k}^2{R}^2+ikL\right),\hfill \end{array}(106)$

which diverges at the origin. We can also integrate over z ∈ (L/2, L/2), and obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^R\frac{2\pi rdr}{\pi R^2}{\int }_{-L/2}^{L/2}\frac{dz}{L}{\mathrm{\Psi }}^{*}\left({\widehat{l}}_x^2+{\widehat{l}}_y^2\right)\mathrm{\Psi }\hfill \\ \hfill & \hfill & ={\hslash }^2\left(\frac{m^2{L}^2}{6R^2}\left(\ln R-\ln 0\right)+\frac{1}{2}{k}^2{R}^2\right),\hfill \end{array}(107)$

which becomes a real value, but still diverges at the origin.

The position-dependent average of the radial magnitude suggests that it contains extrinsic contributions of OAM. For both coordinates, we could not avoid the ultraviolet divergences at the origin, which are coming from the finite amplitude of the wavefunction at the origin. Without having a node at the origin, the magnitude of the OAM required to sustain the phase described by e^imϕ is impossible to exist.

However, for the LG modes, which always have nodes at the centre (r = 0) of the waveguide for ^∀m ≠ 0, there is a chance that the OAM can be well-defined quantum-mechanically. Our main purpose of this work is to confirm the validity of this concept of OAM, using canonical orbital angular momentum operators defined in this section, for the LG modes in a cylindrical GRIN fibre. In the next section, we will confirm positive results, including the observable nature of the magnitude and the commutation relationship for OAM.

4 Mathematical formulas

4.1 Laguerre function

We describe full details of Laguerre and associate Laguerre functions and related formulas in this section [56–58]. First, we consider a differential equation

$\left[a\frac{d^2}{d{a}^2}+\left(1-a\right)\frac{d}{da}+n\right]f=0,(108)$

which will be solved by assuming a Taylor series expansion

$f={\displaystyle \sum _{j=0}^{\infty }}{a}_j{a}^j,(109)$

which gives

${\displaystyle \sum _{j=0}^{\infty }}{a}^j\left[{\left(j+1\right)}^2{a}_{j+1}+\left(n-j\right)a_j\right]=0.(110)$

This provides a recurrence formula

$a_{j+1}=-\frac{{\left(n-j\right)}^2}{{\left(j+1\right)}^2}{a}_j(111)$

for j = 0, 1, 2, …, n, and a_j = 0 for j > n. Therefore, we obtain

$a_j={\left(-1\right)}^j\frac{n!}{{\left(j!\right)}^2\left(n-j\right)!}{a}_0.(112)$

The differential equation cannot be determined uniquely without providing a boundary condition. The same is true for a special function, such that there exists a room to choose the arbitrary value of a₀, while it is a standard rule to choose a₀ > 0 in mathematics. Our definition in this paper is a₀ = 1, but other people are also using a₀ = n! as an alternative definition. In this way, we obtain the solution, f = f_n(a), as

$L_n\left(a\right)={\displaystyle \sum _{j=0}^n}\frac{{\left(-1\right)}^j}{j!}\frac{n!}{\left(j!\right)\left(n-j\right)!}{a}^j(113)$

$={\displaystyle \sum _{j=0}^n}\frac{{\left(-1\right)}^j}{j!}{}_n{C}_j{a}^j,(114)$

where _nC_j is a binomial coefficient.

4.1.1 Rising operator

By directly calculating the derivative, we obtain

$\left[a\frac{d}{da}-a\right]L_n\left(a\right)={\displaystyle \sum _{j=0}^{n+1}}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n+1\right)!}{j!\left(n+1-j\right)!}j{a}^j,(115)$

which can be combined with this identity

$\left(n+1\right)L_n\left(a\right)={\displaystyle \sum _{j=0}^{n+1}}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n+1\right)!}{j!\left(n+1-j\right)!}\left(n+1-j\right)a^j,(116)$

to obtain

$\left[a\frac{d}{da}-a+\left(n+1\right)\right]L_n\left(a\right)=\left(n+1\right)L_{n+1}\left(a\right).(117)$

This formula works as a raising operator to increase the radial index, n.

4.1.2 Lowering operator

Quite similarly, we can also obtain the lowering operator. Calculating a derivative,

$a\frac{d}{da}{L}_n\left(a\right)={\displaystyle \sum _{j=0}^n}\frac{{\left(-1\right)}^j}{j!}\frac{n!}{j!\left(n-j\right)!}j{a}^j(118)$

$\begin{array}{ccc}\hfill & =\hfill & {\displaystyle \sum _{j=0}^{n-1}}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n-1\right)!}{j!\left(n-1-j\right)!}\frac{n}{\left(n-j\right)}j{a}^j\hfill \\ \hfill & \hfill & +n\frac{{\left(-1\right)}^n}{n!}{a}^n\hfill \end{array}(119)$

together with the identity

$-n{L}_n\left(a\right)={\displaystyle \sum _{j=0}^{n-1}}\frac{{\left(-1\right)}^j}{j!}\frac{n!}{j!\left(n-1-j\right)!}\left(-n\right)\frac{n}{n-j}{a}^j-\frac{{\left(-1\right)}^n}{n!}{a}^n,(120)$

we obtain the lowering operation formula,

$\left[a\frac{d}{da}-n\right]L_n\left(a\right)=-n{L}_{n-1}\left(a\right).(121)$

By summing up raising and lowering operators, we also obtain the recurrence relationship

$\left(2n+1-a\right)L_n\left(a\right)=\left(n+1\right)L_{n+1}\left(a\right)+n{L}_{n-1}\left(a\right),(122)$

which correlate 3 successive Laguerre functions.

4.1.3 Generating function

The generating function is defined as a function, whose coefficients of series expansion are Laguerre functions. Therefore, it is defined as

$G\left(t,\tau \right)={\displaystyle \sum _{i=0}^{\infty }}{L}_i\left(t\right){\tau }^i.(123)$

By inserting the series expansion form of L_i(t), we obtain

$G\left(t,\tau \right)={\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=0}^i}\frac{{\left(-1\right)}^j}{j!}\frac{i!}{j!\left(i-j\right)!}{t}^j{\tau }^i(124)$

$={\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{k=0}^{\infty }}\frac{{\left(-1\right)}^j}{j!}\frac{\left(k+j\right)!}{j!k!}{t}^j{\tau }^{k+j},(125)$

where we used k = i − j in the 2nd line. Together with the binomial theorem

${\left(\frac{1}{1-\tau }\right)}^{j+1}={\displaystyle \sum _{k=0}^{\infty }}\frac{\left(k+j\right)!}{k!j!}{\tau }^k,(126)$

we finally obtain the analytic formula for the generating function as

$G\left(t,\tau \right)=\frac{1}{1-\tau }{\displaystyle \sum _{j=0}^{\infty }}\frac{1}{j!}{\left(-\frac{t\tau }{1-\tau }\right)}^j(127)$

$=\frac{1}{1-\tau }\exp \left(-\frac{t\tau }{1-\tau }\right).(128)$

Using this generating function, we can obtain the orthogonality relationship, which is used to confirm the orthogonality against modes with different radial numbers and calculate the normalisation factors. In order to derive it, we evaluate the following sum of the integrals,

$\begin{array}{ccc}\hfill & \hfill & {\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{i^{\prime }=0}^{\infty }}{\tau }^i{\tau }^{\prime i^{\prime }}{\int }_0^{\infty }dt{\mathrm{e}}^{-t}{L}_i\left(t\right)L_{i^{\prime }}\left(t\right)\hfill \\ \hfill & =\hfill & \frac{1}{\left(1-\tau \right)\left(1-{\tau }^{\prime }\right)}{\int }_0^{\infty }dt\exp \left(-t\frac{1-\tau {\tau }^{\prime }}{\left(1-\tau \right)\left(1-{\tau }^{\prime }\right)}\right)\hfill \\ \hfill & =\hfill & \frac{1}{1-\tau {\tau }^{\prime }}{\left[\exp \left(-t\frac{1-\tau {\tau }^{\prime }}{\left(1-\tau \right)\left(1-{\tau }^{\prime }\right)}\right)\right]}_0^{\infty }\hfill \\ \hfill & =\hfill & \frac{1}{1-\tau {\tau }^{\prime }}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _i^{\infty }}{\left(\tau {\tau }^{\prime }\right)}^i\hfill \\ \hfill & =\hfill & {\displaystyle \sum _i^{\infty }}{\displaystyle \sum _{i^{\prime }}^{\infty }}{\delta }_{i,i^{\prime }}{\tau }^i{\tau }^{\prime i^{\prime }}.(129)\hfill \end{array}$

Comparing the first term and the last one, we obtain

${\int }_0^{\infty }dt{\mathrm{e}}^{-t}{L}_i\left(t\right)L_{i^{\prime }}\left(t\right)={\delta }_{i,i^{\prime }}.(130)$

4.1.4 Rodrigues formula

Rodrigues formula is an operator form of the representation of Laguerre function, which will be suitable for quantum mechanics. In order to obtain it, we just need to evaluate the following function

$\begin{array}{ccc}\hfill {\mathrm{e}}^t\frac{d^n}{d{t}^n}\left(t^n{\mathrm{e}}^{-t}\right)& =\hfill & {\displaystyle \sum _{j=0}^n}\frac{n!}{j!\left(n-j!\right)}\frac{n!}{\left(n-j!\right)}{\left(-1\right)}^{n-j}{t}^{n-j}\hfill \\ \hfill & =\hfill & n!{\displaystyle \sum _{j=0}^n}\frac{{\left(-1\right)}^j}{j!}\frac{n!}{j!\left(n-j!\right)}{t}^j\hfill \\ \hfill & =\hfill & n!L_n\left(t\right),(131)\hfill \end{array}$

and thus, we obtain

$L_n\left(t\right)=\frac{1}{n!}{\mathrm{e}}^t\frac{d^n}{d{t}^n}\left(t^n{\mathrm{e}}^{-t}\right).(132)$

4.2 Associated Laguerre function

The associated Laguerre function is defined as

$L_n^m\left(t\right)={\left(-1\right)}^m\frac{d^m}{d{t}^m}{L}_{n+m}\left(t\right).(133)$

The factor of (−1)^m guarantees the first Taylor series expansion coefficient of a₀ to be positive a₀ > 0 as a mathematical convention.

The differential equation for the associated Laguerre function is derived from that of the Laguerre function

$\left[t\frac{d^2}{d{t}^2}+\left(1-t\right)\frac{d}{dt}+n\right]L_n\left(t\right)=0,(134)$

by the mth derivative of this equation,

$\left[\frac{d^m}{d{t}^m}\left(t\frac{d^2}{d{t}^2}\right)+\frac{d^m}{d{t}^m}\left(\left(1-t\right)\frac{d}{dt}\right)+n\frac{d^m}{d{t}^m}\right]L_n\left(t\right)=0,(135)$

which becomes,

$\left[t\frac{d^2}{d{t}^2}+\left(m+1-t\right)\frac{d}{dt}+\left(n-m\right)\right]{\left(-1\right)}^m\frac{d^m}{d{t}^m}{L}_n\left(t\right)=0.(136)$

By exchanging n → n + m, we obtain

$\left[t\frac{d^2}{d{t}^2}+\left(m+1-t\right)\frac{d}{dt}+n\right]L_n^m\left(t\right)=0.(137)$

We also obtain the Taylor series expansion of $L_n^m(t)$ by direct calculation. Inserting the series expression for L_n(t) into the definition, we obtain

$L_n^m\left(t\right)={\left(-1\right)}^m\frac{d^m}{d{t}^m}{\displaystyle \sum _{j=0}^{n+m}}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n+m\right)!}{j!\left(n+m-j\right)!}{t}^j(138)$

$\begin{array}{ccc}\hfill & =\hfill & {\left(-1\right)}^m{\displaystyle \sum _{j=m}^{n+m}}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n+m\right)!}{j!\left(n+m-j\right)!}\frac{j!}{\left(j-m\right)!}{t}^{j-m}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{j=0}^n}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n+m\right)!}{\left(j+m\right)!\left(n-j\right)!}{t}^j,(139)\hfill \end{array}$

which shows that the term at j = 0 is indeed positive.

4.2.1 Generating function

The generating function for the associated Laguerre function is defined by

$G\left(t,\tau \right)={\displaystyle \sum _{i=0}^{\infty }}{L}_i^m\left(t\right){\tau }^i(140)$

$\begin{array}{ccc}\hfill & =\hfill & {\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=0}^n}\frac{{\left(-1\right)}^j}{j!}\frac{\left(n+m\right)!}{\left(j+m\right)!\left(n-j\right)!}{t}^j{\tau }^i\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{j=0}^{\infty }}{\displaystyle \sum _{k=0}^{\infty }}\frac{{\left(-1\right)}^j}{j!}\frac{\left(j+k+m\right)!}{\left(j+m\right)!k!}{t}^j{\tau }^{j+k}.\hfill \end{array}(141)$

Using the binomial theorem,

${\left(\frac{1}{1-\tau }\right)}^{j+m+1}={\displaystyle \sum _{k=0}^{\infty }}\frac{\left(k+j+m\right)!}{k!\left(j+m\right)!}{x}^k,(142)$

we obtain

$G\left(t,\tau \right)={\displaystyle \sum _{j=0}^{\infty }}\frac{{\left(-1\right)}^j}{j!}{\left(\frac{1}{1-\tau }\right)}^{j+m+1}{\left(t\tau \right)}^j(143)$

$={\left(\frac{1}{1-\tau }\right)}^{m+1}\sum _{j=0}^{\infty }\frac{{\left(-1\right)}^j}{j!}{\left(-\frac{t\tau }{1-\tau }\right)}^j(144)$

$={\left(\frac{1}{1-\tau }\right)}^{m+1}\exp \left(-\frac{t\tau }{1-\tau }\right).(145)$

4.2.2 Recurrence relationship

We obtain the recurrence relationship for the associated Laguerre function. By calculating the derivative of the generating function by τ, we obtain

$\begin{array}{ccc}\hfill & \hfill & {\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right)n{\tau }^{n-1}\hfill \\ \hfill & =\hfill & \left(m+1\right)\frac{{\mathrm{e}}^{-\frac{t\tau }{1-\tau }}}{{\left(1-\tau \right)}^{m+2}}-\frac{{\mathrm{e}}^{-\frac{t\tau }{1-\tau }}}{{\left(1-\tau \right)}^{m+1}}\frac{t\left(1-\tau \right)+t\tau }{{\left(1-\tau \right)}^2}\hfill \\ \hfill & =\hfill & \frac{m+1}{1-\tau }{\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^n-\frac{t}{{\left(1-\tau \right)}^2}{\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^n.\hfill \end{array}(146)$

Then, we obtain

$\begin{array}{ccc}\hfill & \hfill & {\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right)n{\tau }^{n-1}{\left(1-\tau \right)}^2\hfill \\ \hfill & =\hfill & \frac{m+1}{1-\tau }{\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^n\left(1-\tau \right)-t{\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^n,\hfill \end{array}(147)$

from which we obtain the recurrence relationship

$\begin{array}{ccc}\hfill & \hfill & \left(n+1\right)L_{n+1}^m\left(t\right)-\left(2n+m+1-t\right)L_n^m\left(t\right)\hfill \\ \hfill & +\left(n+m\right)L_{n-1}^m\left(t\right)=0\hfill \end{array}(148)$

for n ≥ 1.

4.2.3 Ladder operators for radial quantum number

For obtaining ladder operators, we calculate the derivative of the generating function by t as

${\displaystyle \sum _{n=0}^{\infty }}\frac{d}{dt}{L}_n^m\left(t\right){\tau }^n=-\frac{\tau }{1-\tau }{\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^n,(149)$

which becomes

${\displaystyle \sum _{n=0}^{\infty }}\frac{d}{dt}{L}_n^m\left(t\right){\tau }^n-{\displaystyle \sum _{n=0}^{\infty }}\frac{d}{dt}{L}_n^m\left(t\right){\tau }^{n+1}=-{\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^{n+1}.(150)$

Then, we obtain the identity for lowering n

$\frac{d}{dt}{L}_n^m\left(t\right)=\left[\frac{d}{dt}-1\right]L_{n-1}^m\left(t\right).(151)$

However, this expression is not perfect, since the derivative operator remained in the right-hand side, which will be removed later.

Next, we construct the raising operator by calculating the derivative of the recurrence equation by t as

$\begin{array}{ccc}\hfill & \hfill & \left(n+1\right)\left[\frac{d}{dt}-1\right]L_n^m\left(t\right)+L_n^m\left(t\right)\hfill \\ \hfill & \hfill & -\left(2n+m+1-t\right)\frac{d}{dt}{L}_n^m\left(t\right)\hfill \\ \hfill & \hfill & +\left(n+m\right)\frac{d}{dt}{L}_{n-1}^m\left(t\right)=0,\hfill \end{array}(152)$

which becomes

$\left[\left(n+m-t\right)\frac{d}{dt}+n\right]L_n^m\left(t\right)=\left[\left(n+m\right)\frac{d}{dt}\right]L_{n-1}^m\left(t\right).(153)$

by using the lowering identity, this becomes

$\begin{array}{ccc}\hfill & \hfill & \left(n+m-t\right)\left[\frac{d}{dt}-1\right]L_{n-1}^m\left(t\right)+n{L}_n^m\left(t\right)\hfill \\ \hfill & \hfill & =\left(n+m\right)\frac{d}{dt}{L}_{n-1}^m\left(t\right),\hfill \end{array}(154)$

which gives the raising operator

$\left[t\frac{d}{dt}-t+n+m+1\right]L_n^m\left(t\right)=\left(n+1\right)L_{n+1}^m\left(t\right)(155)$

for n ≥ 1. This expression is preferable, since the derivative operation appeared only in the left-side. Together with this raising operator and the recurrence formula, we can eliminate $L_{n+1}^m(t)$ to obtain the lowering operator

$\left[t\frac{d}{dt}-n\right]L_n^m\left(t\right)=-\left(n+m\right)L_{n-1}^m\left(t\right),(156)$

while keeping m unchanged. These ladder operations for the associated Laguerre functions are consistent with those for Laguerre function in the limit of m = 0.

4.2.4 Ladder operators for orbital angular momentum

The above formulas for raising and lowering the radial quantum numbers are known in literatures [56–58], while we could not find appropriate formulas for raising and lowering quantum number m for orbital angular momentum without affecting the radial quantum number of n. Here, we derived these by direct calculations. First, we obtain the raising operator by calculating

$L_n^{m+1}\left(t\right)={\left(-1\right)}^{m+1}\frac{d^{m+1}}{d{t}^{m+1}}{L}_{n+m+1}\left(t\right)(157)$

$=-\frac{d}{dt}{\left(-1\right)}^m\frac{d^m}{d{t}^m}{L}_{\left(n+1\right)+m}\left(t\right)(158)$

$=-\frac{d}{dt}{L}_{n+1}^m\left(t\right)(159)$

$=-\left[\frac{d}{dt}-1\right]L_n^m\left(t\right).(160)$

Thus, the raising operator is described as

$\left[\frac{d}{dt}-1\right]L_n^m\left(t\right)=-L_n^{m+1}\left(t\right).(161)$

It was less straightforward to obtain the lowering operator. As for preparations, we recognised several useful recurrence formulas

$L_n^{m+1}\left(t\right)=-\left[\frac{d}{dt}-1\right]L_n^m\left(t\right),(162)$

$L_n^{m+1}\left(t\right)=-\frac{d}{dt}{L}_{n+1}^m\left(t\right),(163)$

$L_n^{m+1}\left(t\right)=L_{n-1}^{m+1}\left(t\right)+L_n^m\left(t\right),(164)$

$L_n^m\left(t\right)=L_{n-1}^m\left(t\right)+L_n^{m-1}\left(t\right).(165)$

By using the recurrence formula, we obtain

$L_{n-1}^m\left(t\right)=\frac{2n+m+1-t}{n+m}{L}_n^m\left(t\right)-\frac{n+1}{n+m}{L}_{n+1}^m\left(t\right).(166)$

Next, we use the raising operator for n to obtain

$L_{n+1}^m\left(t\right)=\frac{1}{n+1}\left[t\frac{d}{dt}-t+n+m+1\right]L_n^m\left(t\right).(167)$

By combining these equations, we obtain

$L_{n-1}^m\left(t\right)=\frac{1}{n+m}\left[n-t\frac{d}{dt}\right]L_n^m\left(t\right).(168)$

Inserting this into Eq. (B33), we obtain the lowering operator

$\left[t\frac{d}{dt}+m\right]L_n^m\left(t\right)=\left(n+m\right)L_n^{m-1}\left(t\right).(169)$

4.2.5 Orthogonality relationship

The orthogonality relationship is obtained in a similar way by using the generating function

${\displaystyle \sum _{n=0}^{\infty }}{L}_n^m\left(t\right){\tau }^n=\frac{{\mathrm{e}}^{-\frac{t\tau }{1-\tau }}}{{\left(1-\tau \right)}^{m+1}},(170)$

and calculating the sum

$\begin{array}{ccc}\hfill & \hfill & {\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \sum _{n^{\prime }=0}^{\infty }}{\tau }^n{\tau }^{\prime n^{\prime }}{\int }_0^{\infty }dt{\mathrm{e}}^{-t}{t}^m{L}_n^m\left(t\right)L_{n^{\prime }}^m\left(t\right)\hfill \\ \hfill & \hfill & =\frac{1}{{\left(1-\tau \right)}^{m+1}{\left(1-{\tau }^{\prime }\right)}^{m+1}}{\int }_0^{\infty }dt{t}^m{\mathrm{e}}^{-t\frac{1-\tau {\tau }^{\prime }}{\left(1-\tau \right)\left(1-{\tau }^{\prime }\right)}}\hfill \\ \hfill & \hfill & =\frac{m!}{{\left(1-\tau {\tau }^{\prime }\right)}^{m+1}}\hfill \\ \hfill & \hfill & ={\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \sum _{n^{\prime }=0}^{\infty }}\frac{\left(n+m\right)!}{n!}{\delta }_{n,n^{\prime }}{\tau }^n{\tau }^{\prime n^{\prime }},\hfill \end{array}(171)$

where we used the binomial theorem

${\left(\frac{1}{1-\tau {\tau }^{\prime }}\right)}^{m+1}={\displaystyle \sum _{n=0}^{\infty }}\frac{\left(n+m\right)!}{n!m!}{\tau }^n{\tau }^{\prime n}.(172)$

Thus, we obtain the orthogonality relationship

${\int }_0^{\infty }dt{\mathrm{e}}^{-t}{t}^m{L}_n^m\left(t\right)L_{n^{\prime }}^m\left(t\right)=\frac{\left(n+m\right)!}{n!}{\delta }_{n,n^{\prime }}.(173)$

4.2.6 Rodrigues formula

The Rodrigues formula for the associated Laguerre function is obtained by the direct calculations. First, we calculate

$\begin{array}{ccc}\hfill L_n^m\left(t\right)& =\hfill & {\left(-1\right)}^m\frac{d^m}{d{t}^m}{L}_{n+m}\left(t\right)\hfill \\ \hfill & =\hfill & {\left(-1\right)}^m\frac{d^m}{d{t}^m}\left[\frac{1}{\left(n+m\right)!}{\mathrm{e}}^t\frac{d^{n+m}}{d{t}^{n+m}}\left(t^{n+m}{\mathrm{e}}^{-t}\right)\right]\hfill \\ \hfill & =\hfill & \frac{d^m}{d{t}^m}\left[{\displaystyle \sum _{k=0}^{n+m}}\frac{1}{\left(n+m-k\right)!}{t}^{n+m-k}\right.\left.{\left(-1\right)}^{n-k}\frac{\left(n+m\right)!}{\left(n+m-k\right)!k!}\right]\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{k=0}^n}\frac{1}{\left(n+m-k\right)!}\frac{\left(n+m-k\right)!}{\left(n-k\right)!}{t}^{n-k}\hfill \\ \hfill & \hfill & {\left(-1\right)}^{n-k}\frac{\left(n+m\right)!}{\left(n+m-k\right)!k!}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{k=0}^n}\frac{\left(n+m\right)!}{\left(n+m-k\right)!k!\left(n-k\right)!}{\left(-t\right)}^{n-k}\hfill \end{array}(174)$

On the other hand, we calculate

$\begin{array}{ccc}\hfill & \hfill & \frac{1}{n!}{t}^{-m}{\mathrm{e}}^t\frac{d^n}{d{t}^n}\left(t^{n+m}{\mathrm{e}}^{-t}\right)\hfill \\ \hfill & ={\displaystyle \sum _{k=0}^n}\frac{\left(n+m\right)!}{k!\left(n-k\right)!\left(n+m-k\right)!}{\left(-t\right)}^{n-k}\hfill \end{array}(175)$

By comparison, we obtain the Rodrigues formula

$L_n^m\left(t\right)=\frac{1}{n!}{t}^{-m}{\mathrm{e}}^t\frac{d^n}{d{t}^n}\left(t^{n+m}{\mathrm{e}}^{-t}\right).(176)$

4.2.7 Integration formulas

We also obtained integration formulas for the associated Laguerre functions:

${\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^m{L}_n^m\left(a\right)L_{n^{\prime }}^m\left(a\right)=\frac{\left(n+m\right)!}{n!}{\delta }_{n,n^{\prime }},(177)$

${\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^m{L}_n^m\left(a\right)L_n^{m-1}\left(a\right)=\frac{\left(n+m\right)!}{n!},(178)$

${\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{m+1}{L}_n^{m+1}\left(a\right)L_n^m\left(a\right)=\frac{\left(n+m+1\right)!}{n!},(179)$

${\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{m+1}{L}_n^m\left(a\right)L_n^m\left(a\right)=\frac{\left(n+m\right)!}{n!}\left(2n+m+1\right).(180)$

These are useful to calculate the matrix elements. We also obtained an identity,

$\begin{array}{ccc}\hfill & \hfill & m{\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{m-1}{L}_n^m\left(a\right)L_n^m\left(a\right)+\frac{\left(n+m\right)!}{n!}\hfill \\ \hfill & \hfill & =2{\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^m{L}_n^m\left(a\right)L_n^{m+1}\left(a\right).\hfill \end{array}(181)$

5 Results and discussions

We consider applications of the canonical OAM operators to the LG modes in a GRIN fibre. We use the normalised LG mode [5, 6, 8–10]

${\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)=⟨r,\varphi ,z|{\mathrm{\Psi }}_n^m⟩(182)$

$\begin{array}{ccc}\hfill & \hfill & =\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+|m|\right)!}}{\left(\frac{\sqrt{2}r}{w_0}\right)}^{|m|}{L}_n^{|m|}\left(2{\left(\frac{r}{w_0}\right)}^2\right){\mathrm{e}}^{-\frac{r^2}{w_0^2}}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz},\hfill \end{array}(183)$

where n is the radial quantum number and m is the quantum number of OAM along the direction of propagation. In principle, we should also consider a superposition state made of the LG modes with different quantum numbers [59–61], but we will not consider this in this paper for simplicity; yet our formalism works well. We define a normalised cross-sectional area as $a=2r^2/w_0^2$ to simplify calculations. We use various formulas for associate Laguerre functions, which are summarised in the previous section.

5.1 Expectation value

First, we have checked the expectation values of $\widehat{\mathbf{l}}$ by use of the LG modes. This was straightforward by noting ${\widehat{l}}_{\pm }{\psi }_n^m\propto {\psi }_n^{m\pm 1}\propto {\mathrm{e}}^{i(m\pm 1)\varphi }$ and

${\int }_0^{2\pi }\frac{d\varphi }{2\pi }{\mathrm{e}}^{\pm i\varphi }=0,(184)$

and thus we obtain

$⟨{\widehat{l}}_x⟩=0,(185)$

$⟨{\widehat{l}}_y⟩=0,(186)$

$⟨{\widehat{l}}_z⟩=\hslash m.(187)$

Therefore, the quantum-mechanical expectation value of OAM is well-defined for all directions. This is a single particle expectation value, and the total angular momentum, ${\widehat{L}}_z$, along z for a coherent state is obtained by multiplying N as $⟨{\widehat{L}}_z⟩=\hslash mN$.

5.2 Ladder operations

We will evaluate ladder operations to the LG modes. The calculations are straightforward but tedious, so that we will split them into several sections.

5.2.1 Rising operation for m > 0

We assume m > 0 and calculate

$\begin{array}{ccc}\hfill {\widehat{l}}_{+}{\mathrm{\Psi }}_n^m=& \hfill & \hslash {\mathrm{e}}^{i\left(m+1\right)\varphi }{\mathrm{e}}^{ikz}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{\left(\frac{\sqrt{2}r}{w_0}\right)}^m{\mathrm{e}}^{-\frac{r^2}{w_0^2}}\hfill \\ \hfill & \hfill & \cdot \left(z\left(\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}-ikr\right)L_n^m\left(a\right),\hfill \end{array}(188)$

Where the last factor becomes

$\begin{array}{ccc}\hfill & \hfill & \left(z\left(\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}-ikr\right)L_n^m\left(a\right)\hfill \\ \hfill & \hfill & =2\sqrt{2}\frac{z}{w_0}\left(\frac{\sqrt{2}r}{w_0}\right)\left(\frac{d}{da}{L}_n^m\left(a\right)-\frac{1}{2}{L}_n^m\left(a\right)-\frac{ik{w}_0^2}{4z}{L}_n^m\left(a\right)\right)\hfill \\ \hfill & \hfill & =-2\sqrt{2}\frac{z}{w_0}\left(\frac{\sqrt{2}r}{w_0}\right)\left(L_n^{m+1}\left(a\right)-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)L_n^m\left(a\right)\right),\hfill \end{array}(189)$

Where we have used

$\left[\frac{d}{da}-1\right]L_n^m\left(a\right)=-L_n^{m+1}\left(a\right).(190)$

Therefore, we obtain

$\begin{array}{ccc}\hfill {\widehat{l}}_{+}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)=& -2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n+m+1}{\mathrm{\Psi }}_n^{m+1}\left(r,\varphi ,z\right)\hfill & \hfill \\ \hfill & \cdot \left(1-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)\frac{L_n^m\left(a\right)}{L_n^{m+1}\left(a\right)}\right).\hfill & \hfill \end{array}(191)$

The most significant part of this expression is that we confirm ${\widehat{l}}_{+}{\mathrm{\Psi }}_n^m\propto {\mathrm{\Psi }}_n^{m+1}$, which means that the raising operator properly works to increase the quantum number m of OAM. Unfortunately, the coefficient is not a constant, which depends on both z and r through a. Therefore, the shape of the orbital would be significantly distorted upon the application of ${\widehat{l}}_{+}$. Nevertheless, the main role of ${\widehat{l}}_{+}$ to increase m was successfully confirmed for m > 0.

5.2.2 Rising operation for m = 0

We then continue to calculate for m = 0 as

$\begin{array}{ccc}\hfill {\widehat{l}}_{+}{\mathrm{\Psi }}_n^0=& \hslash {\mathrm{e}}^{i\varphi }{\mathrm{e}}^{ikz}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+1\right)!}}\sqrt{n+1}{\mathrm{e}}^{-\frac{r^2}{w_0^2}}\hfill & \hfill \\ \hfill & \cdot \left(z\left(4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r}{w_0^2}\right)-ikr\right)L_n^0\left(a\right),\hfill & \hfill \end{array}(192)$

where the last factor becomes

$\begin{array}{ccc}\hfill & \hfill & \left(z\left(+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r}{w_0^2}\right)-ikr\right)L_n^0\left(a\right)\hfill \\ \hfill & \hfill & =2\sqrt{2}\frac{z}{w_0}\left(\frac{\sqrt{2}r}{w_0}\right)\left(\frac{d}{da}{L}_n^0\left(a\right)-\frac{1}{2}{L}_n^0\left(a\right)-\frac{ik{w}_0^2}{4z}{L}_n^0\left(a\right)\right)\hfill \\ \hfill & \hfill & =-2\sqrt{2}\frac{z}{w_0}\left(\frac{\sqrt{2}r}{w_0}\right)\left(L_n^1\left(a\right)-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)L_n^0\left(a\right)\right),\hfill \end{array}(193)$

and thus we obtain

$\begin{array}{ccc}\hfill {\widehat{l}}_{+}{\mathrm{\Psi }}_n^0\left(r,\varphi ,z\right)=& -2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n+1}{\mathrm{\Psi }}_n^1\left(r,\varphi ,z\right)\hfill & \hfill \\ \hfill & \cdot \left(1-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)\frac{L_n^0\left(a\right)}{L_n^1\left(a\right)}\right),\hfill & \hfill \end{array}(194)$

which is exactly the same expression with that of m > 0.

5.2.3 Rising operation for m < 0

We obtain a similar result by directly calculating

$\begin{array}{ccc}\hfill & \hfill & {\widehat{l}}_{+}{\mathrm{\Psi }}_n^m=\hslash {\mathrm{e}}^{i\left(m+1\right)\varphi }{\mathrm{e}}^{ikz}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m\right)!}}{\left(\frac{\sqrt{2}r}{w_0}\right)}^{-m}{\mathrm{e}}^{-\frac{r^2}{w_0^2}}\hfill \\ \hfill & \hfill & \cdot \left(z\left(-\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}-ikr\right)L_n^{-m}\left(a\right),\hfill \end{array}(195)$

where the last factor becomes

$\begin{array}{ccc}\hfill & \hfill & \left(z\left(-\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}-ikr\right)L_n^{-m}\left(a\right)\hfill \\ \hfill & \hfill & =\sqrt{2}\frac{z}{w_0}\left(\frac{w_0}{\sqrt{2}r}\right)\hfill \\ \hfill & \hfill & \left(2\left(\frac{d}{da}-m\right)L_n^{-m}\left(a\right)-a{L}_n^{-m}\left(a\right)-\frac{ik{w}_0^2}{z}{L}_n^{-m}\left(a\right)\right)\hfill \\ \hfill & \hfill & =2\sqrt{2}\frac{z}{w_0}\left(n-m\right)\left(\frac{w_0}{\sqrt{2}r}\right)L_n^{-m-1}\left(a\right)\hfill \\ \hfill & \hfill & \left(1-\frac{a}{2}\frac{1}{n-m}\left(1+\frac{ik{w}_0^2}{2z}\right)\frac{L_n^{|m|}\left(a\right)}{L_n^{|m+1|}\left(a\right)}\right),\hfill \end{array}(196)$

where we have used

$\left[a\frac{d}{da}-m\right]L_n^{-m}\left(a\right)=\left(n-m\right)L_n^{-m-1}\left(a\right).(197)$

Thus, we obtain

$\begin{array}{ccc}\hfill {\widehat{l}}_{+}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)=& \hfill & 2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n-m}{\mathrm{\Psi }}_n^{m+1}\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & \cdot \left(1-\frac{a}{2}\frac{1}{n-m}\left(1+\frac{ik{w}_0^2}{2z}\right)\frac{L_n^{|m|}\left(a\right)}{L_n^{|m+1|}\left(a\right)}\right).\hfill \end{array}(198)$

therefore, the raising operator is successfully working to increment m, independently of the value and the sign of m.

5.2.4 Lowering operation for m < 0

Next, we apply the lowering operator to the case for m < 0, and obtain

$\begin{array}{ccc}\hfill {\widehat{l}}_{-}{\mathrm{\Psi }}_n^m& \hfill & =\hslash {\mathrm{e}}^{i\left(m-1\right)\varphi }{\mathrm{e}}^{ikz}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m+1\right)!}}\hfill \\ \hfill & \hfill & \sqrt{n-m+1}{\left(\frac{\sqrt{2}r}{w_0}\right)}^{-m}\hfill \\ \hfill & \hfill & \cdot {\mathrm{e}}^{-\frac{r^2}{w_0^2}}\left(-z\left(-\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)\right.\hfill \\ \hfill & \hfill & \left.-z\frac{m}{r}+ikr\right)L_n^{-m}\left(a\right),\hfill \end{array}(199)$

Where the last factor becomes

$\begin{array}{ccc}\hfill & \hfill & \left(-z\left(-\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}+ikr\right)L_n^{-m}\left(a\right)\hfill \\ \hfill & =\hfill & -4z\frac{r}{w_0^2}\left(\left(\frac{d}{da}-\frac{1}{2}\right)L_n^{-m}\left(a\right)-a{L}_n^m\left(a\right)-\frac{ik{w}_0^2}{4z}{L}_n^{-m}\left(a\right)\right)\hfill \\ \hfill & =\hfill & 2\sqrt{2}z\frac{\sqrt{2}r}{w_0^2}{L}_n^{-m+1}\left(a\right)\left(1-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)\frac{L_n^{|m|}\left(a\right)}{L_n^{|m-1|}\left(a\right)}\right),\hfill \end{array}(200)$

Which yields

$\begin{array}{ccc}\hfill {\widehat{l}}_{-}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)=& \hfill & 2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n-m+1}{\mathrm{\Psi }}_n^{m-1}\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & \cdot \left(1-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)\frac{L_n^{|m|}\left(a\right)}{L_n^{|m-1|}\left(a\right)}\right).\hfill \end{array}(201)$

this also shows that the lowering operator can actually lower m to m − 1.

5.2.5 Lowering operation for m = 0

Similarly, we obtain for m = 0:

$\begin{array}{ccc}\hfill {\widehat{l}}_{-}{\mathrm{\Psi }}_n^0=& \hslash {\mathrm{e}}^{-i\varphi }{\mathrm{e}}^{ikz}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+1\right)!}}\sqrt{n+1}{\mathrm{e}}^{-\frac{r^2}{w_0^2}}\hfill & \hfill \\ \hfill & \cdot \left(-z\left(4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r}{w_0^2}\right)+ikr\right)L_n^0\left(a\right),\hfill & \hfill \end{array}(202)$

Where the last factor becomes

$\begin{array}{ccc}\hfill & \hfill & \left(-z\left(4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r}{w_0^2}\right)-ikr\right)L_n^0\left(a\right)\hfill \\ \hfill & \hfill & =-2\sqrt{2}\frac{\sqrt{2}r}{w_0}\left(L_n^1\left(a\right)-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)L_n^0\left(a\right)\right),\hfill \end{array}(203)$

Which yields

$\begin{array}{ccc}\hfill {\widehat{l}}_{-}{\mathrm{\Psi }}_n^0\left(r,\varphi ,z\right)=& \hfill & -2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n+1}{\mathrm{\Psi }}_n^{-1}\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & \left(1-\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)\frac{L_n^0\left(a\right)}{L_n^1\left(a\right)}\right).\hfill \end{array}(204)$

This is the same formula as that we obtained for m < 0.

5.2.6 Lowering operation for m > 0

Finally, we apply the lowering operator to the case for m > 0, and obtain

$\begin{array}{ccc}\hfill & \hfill & {\widehat{l}}_{-}{\mathrm{\Psi }}_n^m\hfill \\ \hfill & \hfill & =\hslash {\mathrm{e}}^{i\left(m-1\right)\varphi }{\mathrm{e}}^{ikz}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m-1\right)!}}\frac{1}{\sqrt{n+m}}{\left(\frac{\sqrt{2}r}{w_0}\right)}^m\hfill \\ \hfill & \hfill & \cdot {\mathrm{e}}^{-\frac{r^2}{w_0^2}}\left(-z\left(\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}+ikr\right)L_n^m\left(a\right),\hfill \end{array}(205)$

where the last factor becomes

$\begin{array}{ccc}\hfill & \hfill & \left(-z\left(\frac{m}{r}+4\frac{r}{w_0^2}\frac{d}{da}-2\frac{r^2}{w_0^2}\right)-z\frac{m}{r}+ikr\right)L_n^m\left(a\right)\hfill \\ \hfill & \hfill & =-\frac{z}{r}\left(2\left(a\frac{d}{da}+m\right)-a-\frac{ik{r}^2}{z}\right)L_n^m\left(a\right)\hfill \\ \hfill & \hfill & =-2\frac{z}{r}\left(\left(n+m\right)L_n^{m-1}\left(a\right)-\frac{r^2}{w_0^2}\left(1+\frac{ik{w}_0^2}{2z}\right)L_n^m\left(a\right)\right),\hfill \end{array}(206)$

which yields

$\begin{array}{ccc}\hfill {\widehat{l}}_{-}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)=& \hfill & -2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n+m}{\mathrm{\Psi }}_n^{m-1}\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & \cdot \left(1-\frac{1}{n+m}\frac{a}{2}\left(1+\frac{ik{w}_0^2}{2z}\right)\frac{L_n^m\left(a\right)}{L_n^{m-1}\left(a\right)}\right).\hfill \end{array}(207)$

therefore, the lowering operator is also successfully working to decrement m, independently of the value and the sign of m.

Thus, in this subsection, by obtaining the wavefunction after the ladder operations, we confirmed that the ladder operations work to change the quantised OAM along the propagation direction in units of ℏ.

5.3 Norm after ladder operations

In this subsection, we obtain the norm of the wavefunctions after ladder operations. We calculate it for separately depending on the sign of m, as in the previous subsection. We have extensively used the integration formulas, which are summarised in Section 4.

5.3.1 Norm of ${\widehat{l}}_{+}{\mathrm{\Psi }}_n^m$ for m ≥ 0

First, we rewire the wavefunction by using the formula, $L_n^m(a)=L_n^{m+1}(a)-L_{n-1}^{m+1}(a)$, as

$\begin{array}{ccc}\hfill {\widehat{l}}_{+}{\mathrm{\Psi }}_n^m=& \hfill & -2\sqrt{2}\frac{z}{w_0}\hslash \sqrt{n+m+1}{\psi }_n^{m+1}\hfill \\ \hfill & \hfill & \left(\frac{1}{2}\left(1+\frac{ik{w}_0^2}{2z}\right)+\frac{1}{2}\left(1-\frac{ik{w}_0^2}{2z}\right)\frac{L_n^{m+1}\left(a\right)}{L_n^{m+1}\left(a\right)}\right).\hfill \end{array}(208)$

Then, we obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr|{\widehat{l}}_{+}{\psi }_n^m\left(r,\varphi ,z\right){|}^2\hfill \\ \hfill & \hfill & =2{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n+m+1\right)\left(1-\frac{ik{w}_0^2}{2z}\right)\left(1+\frac{ik{w}_0^2}{2z}\right)\hfill \\ \hfill & \hfill & \cdot \frac{1}{w_0^2}\left({\int }_0^{\infty }2\pi rdr\frac{2}{\pi }\frac{n!}{\left(n+m+1\right)!}{a}^{m+1}{\mathrm{e}}^{-a}{L}_n^{m+1}\left(a\right)L_n^{m+1}\left(a\right)\right.\hfill \\ \hfill & \hfill & \left.+{\int }_0^{\infty }2\pi rdr\frac{2}{\pi }\frac{n!}{\left(n+m+1\right)!}{a}^{m+1}{\mathrm{e}}^{-a}{L}_{n-1}^{m+1}\left(a\right)L_{n-1}^{m+1}\left(a\right)\right)\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+m+1\right).\hfill \end{array}(209)$

5.3.2 Norm of ${\widehat{l}}_{+}{\mathrm{\Psi }}_n^m$ for m < 0

Similarly, we obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr|{\widehat{l}}_{+}{\psi }_n^m\left(r,\varphi ,z\right){|}^2\hfill \\ \hfill & =\hfill & 8{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n-m\right)\frac{n!}{\left(n-m-1\right)!}\frac{1}{w_0^2}{w}_0^2\hfill \\ \hfill & \hfill & \cdot \left({\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{-m-1}{L}_n^{-m-1}\left(a\right)L_n^{-m-1}\left(a\right)\right.\hfill \\ \hfill & \hfill & -\frac{1}{n-m}{\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{-m}{L}_n^{-m-1}\left(a\right)L_n^{-m}\left(a\right)\hfill \\ \hfill & \hfill & +\frac{1}{4{\left(n-m\right)}^2}\left(1+\frac{k^2{w}_0^4}{4z^2}\right)\hfill \\ \hfill & \hfill & \left.{\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{-m+1}{L}_n^{-m}\left(a\right)L_n^{-m}\left(a\right)\right)\hfill \\ \hfill & =\hfill & 8{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n-m\right)\frac{n!}{\left(n-m-1\right)!}\hfill \\ \hfill & \hfill & \cdot \left(\frac{\left(n-m-1\right)!}{n!}-\frac{\left(n-m\right)!}{n!}\frac{1}{n-m}\right.\hfill \\ \hfill & \hfill & \left.+\frac{1}{4{\left(n-m\right)}^2}\left(1+\frac{k^2{w}_0^4}{4z^2}\right)\frac{\left(n-m\right)!}{n!}\left(2n-m+1\right)\right)\hfill \\ \hfill & =\hfill & {\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n-m+1\right).\hfill \end{array}(210)$

5.3.3 Norm of ${\widehat{l}}_{-}{\mathrm{\Psi }}_n^m$ for m ≤ 0

Next, we calculate the norm of ${\widehat{l}}_{-}{\mathrm{\Psi }}_n^m$ for m ≤ 0 as

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr|{\widehat{l}}_{-}{\psi }_n^m\left(r,\varphi ,z\right){|}^2\hfill \\ \hfill & =\hfill & 8{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n-m+1\right)\frac{n!}{\left(n-m+1\right)!}\frac{1}{w_0^2}{w}_0^2\hfill \\ \hfill & \hfill & \cdot {\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{-m-1}\left(L_n^{-m-1}\left(a\right)L_n^{-m-1}\left(a\right)\right.\hfill \\ \hfill & \hfill & -\left.L_n^{-m}\left(a\right)L_n^{-m-1}\left(a\right)+\frac{1}{4}\left(1+\frac{k^2{w}_0^4}{4z^2}\right)L_n^{-m}\left(a\right)L_n^{-m}\left(a\right)\right)\hfill \\ \hfill & =\hfill & 8{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n-m+1\right)\frac{n!}{\left(n-m+1\right)!}\hfill \\ \hfill & \hfill & \cdot \left(\frac{\left(n-m-1\right)!}{n!}-\frac{\left(n-m-1\right)!}{n!}\right.\hfill \\ \hfill & \hfill & \left.+\frac{1}{4}\left(1+\frac{k^2{w}_0^4}{4z^2}\right)\frac{\left(n-m\right)!}{n!}\left(2n-m+1\right)\right)\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n-m+1\right).\hfill \end{array}(211)$

5.3.4 Norm of ${\widehat{l}}_{-}{\mathrm{\Psi }}_n^m$ for m > 0

Finally, we calculate ${\widehat{l}}_{-}{\mathrm{\Psi }}_n^m$ for m > 0 as

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr|{\widehat{l}}_{-}{\psi }_n^m\left(r,\varphi ,z\right){|}^2\hfill \\ \hfill & =\hfill & 8{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n+m\right)\frac{n!}{\left(n+m-1\right)!}\frac{1}{w_0^2}{w}_0^2\hfill \\ \hfill & \hfill & \cdot \left({\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{m-1}{L}_n^{m-1}\left(a\right)L_n^{m-1}\left(a\right)\right.\hfill \\ \hfill & \hfill & -\frac{1}{n+m}{\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^m{L}_n^m\left(a\right)L_n^{m-1}\left(a\right)\hfill \\ \hfill & \hfill & \left.+\frac{1}{4{\left(n+m\right)}^2}\left(1+\frac{k^2{w}_0^4}{4z^2}\right){\int }_0^{\infty }da{\mathrm{e}}^{-a}{a}^{m+1}{L}_n^m\left(a\right)L_n^m\left(a\right)\right)\hfill \\ \hfill & =\hfill & 8{\left(\frac{z}{w_0}\right)}^2{\hslash }^2\left(n+m\right)\frac{n!}{\left(n+m-1\right)!}\hfill \\ \hfill & \hfill & \cdot \left(\frac{\left(n+m-1\right)!}{n!}-\frac{\left(n+m\right)!}{n!}\frac{1}{n+m}\right.\hfill \\ \hfill & \hfill & \left.+\frac{1}{4{\left(n+m\right)}^2}\left(1+\frac{k^2{w}_0^4}{4z^2}\right)\frac{\left(n+m\right)!}{n!}\left(2n+m+1\right)\right)\hfill \\ \hfill & =\hfill & {\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+m+1\right).\hfill \end{array}(212)$

5.3.5 Summary of the norm of ${\widehat{l}}_{\pm }{\mathrm{\Psi }}_n^m$

The above direct calculations show that we obtain the same norm for ${\widehat{l}}_{\pm }{\mathrm{\Psi }}_n^m$:

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr|{\widehat{l}}_{\pm }{\psi }_n^m\left(r,\varphi ,z\right){|}^2\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+|m|+1\right),\hfill \end{array}(213)$

which is independent of the sign of m. The first term is coming from the origin dependent extrinsic OAM, while the second term corresponds to the contribution from the intrinsic OAM, which is always positive and finite. The obtained form of ${(\hslash k{w}_0)}^2/2={(p{w}_0)}^2/2$ with the momentum p = ℏk is intuitively understandable, since pw₀ has the dimension of the angular momentum. As for the case of the plane wave, however, the fact that we obtain $⟨|{\widehat{l}}_{+}{\psi }_n^m{|}^2⟩=⟨|{\widehat{l}}_{-}{\psi }_n^m{|}^2⟩$ means that we cannot confirm the validity of the commutation relationship, and thus, we cannot obtain the expectation value of the magnitude of OAM by simple norm calculations. This might be linked to the fact that the applications of ladder operations contain position dependent factors, such that the orbitals are significantly distorted. In fact, the LG modes are not eigenstates for the magnitude of the OAM, and they are superposition states with different magnitude of the OAM. Nevertheless, we can calculate the expectation value of the magnitude of the OAM, as we shall see below. For that purpose, it is inevitable to calculate $⟨{\widehat{l}}_{+}{\widehat{l}}_{-}⟩$ and $⟨{\widehat{l}}_{-}{\widehat{l}}_{+}⟩$, directly, using the LG modes, which are shown in the next subsection.

5.4 Validity of commutation relationship for LG modes

First, we must calculate the wavefunction of ${\widehat{l}}_{+}{\widehat{l}}_{-}{\mathrm{\Psi }}_n^m$ and ${\widehat{l}}_{-}{\widehat{l}}_{+}{\mathrm{\Psi }}_n^m$ and calculate the expectation value. This is straightforward but tedious. Again, we will split calculations for positive and negative values of m.

5.4.1 Ladder operators

Before calculating the matrix elements, we will describe operators, ${\widehat{l}}_{+}{\widehat{l}}_{-}$ and ${\widehat{l}}_{-}{\widehat{l}}_{+}$, by using cylindrical coordinates (r, ϕ, z), as

$\begin{array}{ccc}\hfill \frac{{\widehat{l}}_{+}{\widehat{l}}_{-}}{{\hslash }^2}& =\hfill & \hslash {\mathrm{e}}^{i\varphi }\left(z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }-r{\partial }_z\right)\hfill \\ \hfill & \hfill & \hslash {\mathrm{e}}^{-i\varphi }\left(-z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }+r{\partial }_z\right)\hfill \\ \hfill & =\hfill & z\left(-z{\partial }_r^2-\frac{z^2}{r^2}i{\partial }_{\varphi }+\frac{z}{r}i{\partial }_r{\partial }_{\varphi }+{\partial }_z+r{\partial }_r{\partial }_z\right)\hfill \\ \hfill & \hfill & +\frac{z}{r}i\left(-z{\partial }_{\varphi }{\partial }_r+\frac{z}{r}i{\partial }_{{\varphi }^2}+r{\partial }_{\varphi }{\partial }_z\right)\hfill \\ \hfill & \hfill & +\frac{z}{r}i\left(-i\right)\left(-z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }+r{\partial }_z\right)\hfill \\ \hfill & \hfill & -r\left(-{\partial }_r-z{\partial }_z{\partial }_r+\frac{i}{r}{\partial }_{\varphi }+\frac{z}{r}i{\partial }_z{\partial }_{\varphi }+r{\partial }_z^2\right)\hfill \\ \hfill & =\hfill & -z^2{\partial }_r^2-\frac{z^2}{r}{\partial }_r-\frac{z^2}{r^2}{\partial }_{\varphi }^2\hfill \\ \hfill & \hfill & +\left(1+2z{\partial }_z\right)r{\partial }_r-r^2{\partial }_z^2+2z{\partial }_z-i{\partial }_{\varphi },\hfill \end{array}(214)$

and

$\begin{array}{ccc}\hfill \frac{{\widehat{l}}_{-}{\widehat{l}}_{+}}{{\hslash }^2}& =\hfill & \hslash {\mathrm{e}}^{-i\varphi }\left(-z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }+r{\partial }_z\right)\hslash \hfill \\ \hfill & \hfill & {\mathrm{e}}^{i\varphi }\left(z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }-r{\partial }_z\right)\hfill \\ \hfill & =\hfill & -z\left(z{\partial }_r^2-\frac{z^2}{r^2}i{\partial }_{\varphi }+\frac{z}{r}i{\partial }_r{\partial }_{\varphi }-{\partial }_z-r{\partial }_r{\partial }_z\right)\hfill \\ \hfill & \hfill & +\frac{z}{r}i\left(z{\partial }_{\varphi }{\partial }_r+\frac{z}{r}i{\partial }_{{\varphi }^2}-r{\partial }_{\varphi }{\partial }_z\right)\hfill \\ \hfill & \hfill & -\frac{z}{r}\left(z{\partial }_r+\frac{z}{r}i{\partial }_{\varphi }-r{\partial }_z\right)\hfill \\ \hfill & \hfill & +r\left({\partial }_r+z{\partial }_z{\partial }_r+\frac{i}{r}{\partial }_{\varphi }+\frac{z}{r}i{\partial }_z{\partial }_{\varphi }-r{\partial }_z^2\right)\hfill \\ \hfill & =\hfill & -z^2{\partial }_r^2-\frac{z^2}{r}{\partial }_r-\frac{z^2}{r^2}{\partial }_{\varphi }^2\hfill \\ \hfill & \hfill & +\left(1+2z{\partial }_z\right)r{\partial }_r-r^2{\partial }_z^2+2z{\partial }_z+i{\partial }_{\varphi }.\hfill \end{array}(215)$

Thus, we also confirmed the commutation relationship for ${\widehat{l}}_{\pm }$,

$\left[{\widehat{l}}_{+},{\widehat{l}}_{-}\right]=2\hslash {\widehat{l}}_z.(216)$

therefore, if we apply these operators to an LG mode, we must confirm that this identity is always valid. This is useful to check the validity of the calculation.

5.4.2 ${\widehat{l}}_{+}{\widehat{l}}_{-}$ operation for m ≥ 0

First, we evaluate various terms as follows:

$\begin{array}{ccc}\hfill {\partial }_r{\mathrm{\Psi }}_n^m& =\hfill & \frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{\left(\frac{\sqrt{2}r}{w_0}\right)}^m{\mathrm{e}}^{-\frac{r^2}{w_0^2}}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \frac{1}{r}\left[m+2a\frac{d}{da}-a\right]L_n^m\left(a\right),\hfill \end{array}(217)$

$\begin{array}{ccc}\hfill r{\partial }_r{\mathrm{\Psi }}_n^m& \hfill & =\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{a}^{m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[2\left(n+m\right)L_n^{m-1}-m{L}_n^m-a{L}_n^m\right],\hfill \end{array}(218)$

$\begin{array}{ccc}\hfill \frac{1}{r}{\partial }_r{\mathrm{\Psi }}_n^m& =\hfill & \frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{a}^{m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \frac{2}{w_0^2}\left[-2L_n^{m+1}+L_n^m+\frac{m}{a}{L}_n^m\right],\hfill \end{array}(219)$

$\begin{array}{ccc}\hfill -z^2{\partial }_r^2{\mathrm{\Psi }}_n^m=& \hfill & -2\frac{z^2}{w_0^2}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{a}^{m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[m\left(m-1\right)\frac{1}{a}{L}_n^m\left(a\right)+2L_n^{m+1}\left(a\right)\right.\hfill \\ \hfill & \hfill & -\left.\left(4n+2m+3\right)L_n^m\left(a\right)+a{L}_n^m\left(a\right)\right],\hfill \end{array}(220)$

$-\frac{z^2}{r^2}{\partial }_{\varphi }^2{\mathrm{\Psi }}_n^m=2\frac{z^2}{w_0^2}{m}^2\frac{1}{a}{\mathrm{\Psi }}_n^m,(221)$

$-r^2{\partial }_z^2{\mathrm{\Psi }}_n^m=\frac{1}{2}{k}^2{w}_0^{2}a{\mathrm{\Psi }}_n^m,(222)$

$\left(r{\partial }_r+2z{\partial }_{z}r{\partial }_r\right){\mathrm{\Psi }}_n^m=\left(1+2ikz\right)r{\partial }_r{\mathrm{\Psi }}_n^m,(223)$

$2z{\partial }_z{\mathrm{\Psi }}_n^m=2ikz{\mathrm{\Psi }}_n^m,(224)$

and finally

$-i{\partial }_{\varphi }{\mathrm{\Psi }}_n^m=m{\mathrm{\Psi }}_n^m.(225)$

By summing up some of these terms, we obtain

$\begin{array}{ccc}\hfill & \hfill & \left[-z^2{\partial }_r^2-\frac{z^2}{r}{\partial }_r-\frac{z^2}{r^2}{\partial }_{\varphi }^2\right]{\mathrm{\Psi }}_n^m\hfill \\ \hfill & \hfill & =2\frac{z^2}{w_0^2}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{a}^{m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[m\left(m-1\right)\frac{1}{a}{L}_n^m\left(a\right)+2L_n^{m+1}\left(a\right)\right.\hfill \\ \hfill & \hfill & -\left(4n+2m+3\right)L_n^m\left(a\right)+a{L}_n^m\left(a\right)\hfill \\ \hfill & \hfill & \left.+2L_n^{m+1}\left(a\right)-L_n^m\left(a\right)-\frac{m}{a}{L}_n^m\left(a\right)+\frac{m^2}{a}{L}_n^m\left(a\right)\right]\hfill \\ \hfill & =\hfill & 2\frac{z^2}{w_0^2}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n+m\right)!}}{a}^{m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[2\left(2n+m+1\right)L_n^m\left(a\right)-a{L}_n^m\left(a\right)\right].\hfill \end{array}(226)$

By using the integration formulas (Section 4), we finally obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr{\left({\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{+}\left(r,\varphi ,z\right){\widehat{l}}_{-}\left(r,\varphi ,z\right){\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & ={\hslash }^{2}2{\left(\frac{z}{w_0}\right)}^2\left(2\left(2n+m+1\right)-\left(2n+m+1\right)\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\left(1+2ikz\right)\left(2\left(n+m\right)-m-\left(2n+m+1\right)\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\frac{1}{2}{\left(k{w}_0\right)}^2\left(2n+m+1\right)+{\hslash }^{2}2ikz+{\hslash }^{2}m\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+m+1\right)+{\hslash }^2\left(m-1\right).\hfill \end{array}(227)$

5.4.3 ${\widehat{l}}_{-}{\widehat{l}}_{+}$ operation for m ≥ 0

The only source of the difference between ${\widehat{l}}_{+}{\widehat{l}}_{-}$ and ${\widehat{l}}_{-}{\widehat{l}}_{+}$ operations is coming from the sign of i∂_ϕ. Therefore, it is straightforward to obtain

$\begin{array}{ccc}\hfill & {\int }_0^{\infty }2\pi rdr{\left({\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{-}\left(r,\varphi ,z\right){\widehat{l}}_{+}\left(r,\varphi ,z\right){\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill & \hfill \\ \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+m+1\right)-{\hslash }^2\left(m+1\right)\hfill & \hfill \end{array}(228)$

This result, together with the previous result for $⟨{\widehat{l}}_{-}{\widehat{l}}_{+}⟩$, confirms that the commutation relationship over average indeed satisfies

$⟨\left[{\widehat{l}}_{+},{\widehat{l}}_{-}\right]⟩=2\hslash ⟨{\widehat{l}}_z⟩,(229)$

where $⟨{\widehat{l}}_z⟩=\hslash m$ for m ≥ 0.

5.4.4 ${\widehat{l}}_{+}{\widehat{l}}_{-}$ operation for m ≤ 0

We can proceed for m ≤ 0 in a similar way. First, we evaluate the factors:

$\begin{array}{ccc}\hfill r{\partial }_r{\mathrm{\Psi }}_n^m=& \hfill & \frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m\right)!}}{a}^{-m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[2\left(n-m\right)L_n^{-m-1}\left(a\right)+m{L}_n^{-m}\left(a\right)-a{L}_n^{-m}\left(a\right)\right],\hfill \end{array}(230)$

$\begin{array}{ccc}\hfill \frac{1}{r}{\partial }_r{\mathrm{\Psi }}_n^m=& \hfill & \frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m\right)!}}{a}^{-m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \frac{2}{w_0^2}\left[-2L_n^{-m+1}\left(a\right)+L_n^{-m}\left(a\right)-\frac{m}{a}{L}_n^{-m}\left(a\right)\right],\hfill \end{array}(231)$

$\begin{array}{ccc}\hfill -z^2{\partial }_r^2{\mathrm{\Psi }}_n^m=& \hfill & -2\frac{z^2}{w_0^2}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m\right)!}}{a}^{-m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[m\left(m+1\right)\frac{1}{a}{L}_n^{-m}\left(a\right)+2L_n^{-m+1}\left(a\right)\right.\hfill \\ \hfill & \hfill & -\left.\left(4n-2m+3\right)L_n^{-m}\left(a\right)+a{L}_n^{-m}\left(a\right)\right],\hfill \end{array}(232)$

and the other factors are the same ones for m ≥ 0. Then, we obtain

$\begin{array}{ccc}\hfill & \hfill & \left[-z^2{\partial }_r^2-\frac{z^2}{r}{\partial }_r-\frac{z^2}{r^2}{\partial }_{\varphi }^2\right]{\mathrm{\Psi }}_n^m\hfill \\ \hfill & =\hfill & 2\frac{z^2}{w_0^2}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m\right)!}}{a}^{-m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[-2L_n^{-m+1}\left(a\right)+\left(4n-2m+3\right)L_n^{-m}\left(a\right)\right.\hfill \\ \hfill & \hfill & -m\left(m+1\right)\frac{1}{a}{L}_n^{-m}\left(a\right)-a{L}_n^{-m}\hfill \\ \hfill & \hfill & 2L_n^{-m+1}\left(a\right)-L_n^{-m}\left(a\right)\hfill \\ \hfill & \hfill & \left.+\frac{m}{a}{L}_n^{-m}\left(a\right)+\frac{m^2}{a}{L}_n^{-m}\left(a\right)\right]\hfill \\ \hfill & =\hfill & 2\frac{z^2}{w_0^2}\frac{1}{w_0}\sqrt{\frac{2}{\pi }\frac{n!}{\left(n-m\right)!}}{a}^{-m/2}{\mathrm{e}}^{-a/2}{\mathrm{e}}^{im\varphi }{\mathrm{e}}^{ikz}\hfill \\ \hfill & \hfill & \cdot \left[2\left(2n+m+1\right)L_n^{-m}\left(a\right)-a{L}_n^{-m}\left(a\right)\right].\hfill \end{array}(233)$

Therefore, we finally obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr{\left({\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{+}\left(r,\varphi ,z\right){\widehat{l}}_{-}\left(r,\varphi ,z\right){\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & ={\hslash }^{2}2{\left(\frac{z}{w_0}\right)}^2\left(2\left(2n-m+1\right)-\left(2n-m+1\right)\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\left(1+2ikz\right)\left(2\left(n-m\right)+m-\left(2n-m+1\right)\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\frac{1}{2}{\left(k{w}_0\right)}^2\left(2n-m+1\right)+{\hslash }^{2}2ikz+{\hslash }^{2}m\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n-m+1\right)+{\hslash }^2\left(m-1\right).\hfill \end{array}(234)$

5.4.5 ${\widehat{l}}_{-}{\widehat{l}}_{+}$ operation for m ≤ 0

Again, the only source of the change between ${\widehat{l}}_{+}{\widehat{l}}_{-}$ and ${\widehat{l}}_{-}{\widehat{l}}_{+}$ operations is coming from the sign of i∂_ϕ, such that we obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr{\left({\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{-}\left(r,\varphi ,z\right){\widehat{l}}_{+}\left(r,\varphi ,z\right){\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n-m+1\right)-{\hslash }^2\left(m+1\right).\hfill \end{array}(235)$

Therefore, the commutation relationship over average

$⟨\left[{\widehat{l}}_{+},{\widehat{l}}_{-}\right]⟩=2\hslash ⟨{\widehat{l}}_z⟩,(236)$

is also valid for m ≤ 0.

5.5 Summary of expectation values of ${\widehat{l}}_{+}{\widehat{l}}_{-}$ and ${\widehat{l}}_{-}{\widehat{l}}_{+}$ operations

The above results are summarised as follows,

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr{\left({\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{+}\left(r,\varphi ,z\right){\widehat{l}}_{-}\left(r,\varphi ,z\right){\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+|m|+1\right)+{\hslash }^2\left(m-1\right),\hfill \end{array}(237)$

and

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr{\left({\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{-}\left(r,\varphi ,z\right){\widehat{l}}_{+}\left(r,\varphi ,z\right){\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & ={\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+|m|+1\right)-{\hslash }^2\left(m+1\right),\hfill \end{array}(238)$

which are independent of the sign of m. The commutation relationship over average

$⟨\left[{\widehat{l}}_{+},{\widehat{l}}_{-}\right]⟩=2\hslash ⟨{\widehat{l}}_z⟩,(239)$

is also valid for ^∀m.

5.6 Magnitude of OAM

The above calculations have confirmed the quantum commutation relationship, $[{\widehat{l}}_{+},{\widehat{l}}_{-}]=2\hslash {\widehat{l}}_z$, as an expectation value after the application to the LG mode. This is trivial, because the commutation relationship is valid at the operator level, such that it should be valid even after the application to the LG mode. Nevertheless, obtained matrix elements of expectation values of ${\widehat{l}}_{+}{\widehat{l}}_{-}$ and ${\widehat{l}}_{-}{\widehat{l}}_{+}$ are useful to evaluate the magnitude of the OAM, because the inner product of the vectorial OAM operator is described as

$\widehat{\mathbf{l}}\cdot \widehat{\mathbf{l}}={\widehat{l}}_x^2+{\widehat{l}}_y^2+{\widehat{l}}_z^2(240)$

$={\widehat{l}}_{+}{\widehat{l}}_{-}+{\widehat{l}}_z^2-\hslash {\widehat{l}}_z(241)$

$={\widehat{l}}_{-}{\widehat{l}}_{+}+{\widehat{l}}_z^2+\hslash {\widehat{l}}_z.(242)$

We can confirm that the expectation value is independent of whether we are using $⟨{\widehat{l}}_{+}{\widehat{l}}_{-}⟩$ or $⟨{\widehat{l}}_{-}{\widehat{l}}_{+}⟩$, and we obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }2\pi rdr{\mathrm{\Psi }}_n^m{\left(r,\varphi ,z\right)}^{*}\widehat{\mathbf{l}}\cdot \widehat{\mathbf{l}}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & =\hfill & {\int }_0^{\infty }2\pi rdr{\mathrm{\Psi }}_n^m{\left(r,\varphi ,z\right)}^{*}{\widehat{l}}_{+}{\widehat{l}}_{-}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)+{\hslash }^{2}m\left(m-1\right)\hfill \\ \hfill & =\hfill & {\int }_0^{\infty }2\pi rdr{\mathrm{\Psi }}_n^m{\left(r,\varphi ,z\right)}^{*}{\widehat{l}}_{-}{\widehat{l}}_{+}{\mathrm{\Psi }}_n^m\left(r,\varphi ,z\right)+{\hslash }^{2}m\left(m+1\right)\hfill \\ \hfill & =\hfill & {\hslash }^2\left(2{\left(\frac{z}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+|m|+1\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\left(m+1\right)\left(m-1\right).\hfill \end{array}(243)$

The expectation value does not depend on the sign of m, which is reasonable in a system with a helical symmetry. The first term contains the origin (z = 0) dependent contribution of the extrinsic OAM. If we take the average over z ∈ (0, L), we obtain

$\begin{array}{ccc}\hfill ⟨\widehat{\mathbf{l}}\cdot \widehat{\mathbf{l}}⟩& =\hfill & {\hslash }^2\left(\frac{2}{3}{\left(\frac{L}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+|m|+1\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\left(m+1\right)\left(m-1\right),\hfill \end{array}(244)$

while if we average over z ∈ (−L/2, L/2), it becomes

$\begin{array}{ccc}\hfill ⟨\widehat{\mathbf{l}}\cdot \widehat{\mathbf{l}}⟩& =\hfill & {\hslash }^2\left(\frac{1}{6}{\left(\frac{L}{w_0}\right)}^2+\frac{1}{2}{\left(k{w}_0\right)}^2\right)\left(2n+|m|+1\right)\hfill \\ \hfill & \hfill & +{\hslash }^2\left(m+1\right)\left(m-1\right).\hfill \end{array}(245)$

The other contributions are from intrinsic OAM. If m ≫ 1, the most of the energy of photons is used to sustain the rotating motion as OAM, such that ℏδω₀m ≫ ℏv₀k in the dispersion relationship. In this limit, the intrinsic OAM is dominated by the contribution from ℏm, which is consistent with the above formula.

In the opposite limit of the absence of the definite quantised OAM (n = m = 0), the beam becomes a simple Gaussian wave. Even in this case,

$⟨{\widehat{l}}_{\text{intrinsic}}^2⟩\to \frac{1}{2}{\left(\hslash k{w}_0\right)}^2-{\hslash }^2\sim \frac{1}{2}{\left(p{w}_0\right)}^2-{\hslash }^2(246)$

holds, which is an intuitive formula, because the same amount of the angular momentum magnitude with spin of ℏ is subtracted. The finite intrinsic OAM is coming from quantum mechanical fluctuations due to ${\widehat{l}}_x$, ${\widehat{l}}_y$, and quantum commutation relationship. Even if the quantum number of the mode is zero (m = 0), the quantum fluctuation is inevitable, such that the finite value of the magnitude remains as the zero-point fluctuation. If the spin component of ℏ is negligible, the OAM fluctuation is of the order of

$\sqrt{⟨{\widehat{l}}^2⟩}\sim p\frac{w_0}{\sqrt{2}}=p{w}_1,(247)$

where $w_1=w_0/\sqrt{2}$ is the effective waist for the OAM. The total amount of fluctuation as a coherent laser beam is obtained by multiplying the number of photons, N.

5.7 Transfer matrix element

Finally, we calculate the matrix elements of ${\widehat{l}}_{\pm }$ among the LG modes with different m. In order to allow these to couple, the phase matching condition must be satisfied, such that the energy, ℏω, and the momentum, p = ℏk, would be conserved through the operation of ${\widehat{l}}_{\pm }$. Otherwise, the transfer matrix elements would vanish due to the destructive interference upon propagation. Strictly speaking, such a condition would not be satisfied during the propagation in the GRIN fibre, since the value of k would be different among modes with different m due to the dispersion relationship. Therefore, the coupling is expected only in the limit of g → 0 as a Gaussian beam, where the dispersion is almost negligible and the material is considered to be almost uniform. We also assume that the beam is sufficiently collimated, such that the impact of the Gouy phase is negligible. We then obtain

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }dr2\pi r{\left({\psi }_n^{m+1}\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{+}{\psi }_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & =-\hslash \frac{z}{w_0}\sqrt{2\left(n+m+1\right)}\left(1+i\frac{k{w}_0^2}{2z}\right)\hfill \end{array}(248)$

for m ≥ 0,

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }dr2\pi r{\left({\psi }_n^{m+1}\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{+}{\psi }_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & =\hslash \frac{z}{w_0}\sqrt{2\left(n+|m|\right)}\left(1-i\frac{k{w}_0^2}{2z}\right)\hfill \end{array}(249)$

for m < 0,

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }dr2\pi r{\left({\psi }_n^{m-1}\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{-}{\psi }_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & =-\hslash \frac{z}{w_0}\sqrt{2\left(n+m\right)}\left(1-i\frac{k{w}_0^2}{2z}\right)\hfill \end{array}(250)$

for m > 0, and

$\begin{array}{ccc}\hfill & \hfill & {\int }_0^{\infty }dr2\pi r{\left({\psi }_n^{m-1}\left(r,\varphi ,z\right)\right)}^{*}{\widehat{l}}_{-}{\psi }_n^m\left(r,\varphi ,z\right)\hfill \\ \hfill & \hfill & =\hslash \frac{z}{w_0}\sqrt{2\left(n+|m|+1\right)}\left(1+i\frac{k{w}_0^2}{2z}\right)\hfill \end{array}(251)$

for m ≤ 0, respectively, If we take the average over z ∈ (0, T), assuming a thickness of T for the optical plate to increment or decrement the value of m, these results show

$\begin{array}{ccc}\hfill {\left(\left.⟨m+1\right|{\widehat{l}}_{+}\left|m\right.⟩\right)}^{*}& =\hfill & -\hslash \sqrt{2\left(n+m+1\right)}\frac{1}{2}\left(\frac{T}{w_0}-ik{w}_0\right)\hfill \\ \hfill & =\hfill & \left.⟨m\right|{\widehat{l}}_{-}\left|m+1\right.⟩\hfill \end{array}(252)$

for m ≥ 0, and

$\begin{array}{ccc}\hfill {\left(\left.⟨m+1\right|{\widehat{l}}_{+}\left|m\right.⟩\right)}^{*}& =\hfill & -\hslash \sqrt{2\left(n+|m|\right)}\frac{1}{2}\left(\frac{T}{w_0}+ik{w}_0\right)\hfill \\ \hfill & =\hfill & \left.⟨m\right|{\widehat{l}}_{-}\left|m+1\right.⟩\hfill \end{array}(253)$

for m < 0. For both cases, the relationship,

${\left(\left.⟨m+1\right|{\widehat{l}}_{+}\left|m\right.⟩\right)}^{*}=\left.⟨m\right|{\widehat{l}}_{-}\left|m+1\right.⟩,(254)$

is always valid for ^∀m. This implies

${\widehat{l}}_{+}^{†}={\widehat{l}}_{-}(255)$

${\widehat{l}}_{-}^{†}={\widehat{l}}_{+}(256)$

for a system described at least in a Hilbert space spanned by the LG modes. These results also suggest that OAM can be a proper quantum mechanical observable, satisfying the commutation relationship, at least for a system described by the LG modes. Experimentally, there are many successful demonstrations to control m [62–68]. In this paper, we have provided a theoretical justification for enabling the increment or decrement of the quantum number, m, thus confirming quantum-mechanical description of the OAM.

5.8 Application to numerical calculations

As an application of this theory, we have numerically calculated how the mode profile is changed up by the ladder operations (Figures 1, 2). Here, we assume a Gaussian profile in the input beam with λ = 1.55 μm in a GRIN fibre with the maximum refractive index at the core of n₀ = 1.50. The refractive index of the core is assumed to be changed dn = 0.05 over the core radius of 25 μm, which corresponds to the GRIN parameter of g = 0.05/25 = 0.002 [(μm)⁻¹]. We assume a horizontally polarised mode, such that the complex electric field has only x-component, as $({\mathcal{E}}_x,{\mathcal{E}}_y)=({\mathrm{\Psi }}_0^0(r,\varphi ,z),0)$. The wavefunction ${\mathrm{\Psi }}_0^0(r,\varphi ,z)$ is normalised to be unity, upon the integration of $|{\mathrm{\Psi }}_0^0(r,\varphi ,z)|$ over space. Figures 1A, B show the amplitude, $|{\mathcal{E}}_x|$, and the phase, ϕ_x, of the calculated input mode, respectively. There is no phase singularity in ϕ_x, and the amplitude is peaked at the centre of the core in the waveguide, as we expected. Then, we have applied the ladder operation to the input, as ${\widehat{l}}_{+}{\mathrm{\Psi }}_0^0(r,\varphi ,z)$, and we have calculated the amplitude and the phase of the output, as shown in Figures 1C, D. We confirmed the typical doughnut shape in the intensity [5, 6, 8–13], while the left vortex is confirmed for increasing the phase along the counter-clock-wise direction. The phase singularity is found at the centre of the mode, where the intensity vanishes. This confirms that our ladder operator worked properly to increment the quantum number of the orbital angular momentum.

FIGURE 1

FIGURE 1. Ladder operation to increment the quantum number for orbital angular momentum. (A) The amplitude and (B) the phase of the input Gaussian beam without a vortex. (C) The amplitude and (D) the phase of the output beam after the ladder operation. The left vortex with topological charge of 1 is generated upon the ladder operation.

FIGURE 2

FIGURE 2. Generation of the left vortex upon passing through a vortex lens. The vortex lens was assumed to be located at z = 0, where the Gaussian input beam is converted to be the left vortex. (A) Real part and (B) imaginary part of the complex electric field, ${\mathcal{E}}_x$, whose intensity was normalised to be unity upon integrating over space. The isosurfaces of $\mathfrak{R}({\mathcal{E}}_x)=0.025$ and −0.025 are shown in red and blue, while the isosurfaces of $\mathfrak{I}({\mathcal{E}}_x)=0.025$ and −0.025 are shown in magenta and cyan. The phase is rotating along the counter-clock-wise direction, seen from the direction of the propagation (z > 0) in the detector side.

We have also confirmed the generation of the left vortex along the direction of the propagation (Figure 2). We have assumed the same parameters to calculate the mode profiles in Figure 1, and considered that a vortex lens [5, 23, 59–61, 67, 69–73], is located at the middle of the direction of the propagation (z) at z = 0. In the calculation, we have used the ladder operation of ${\widehat{l}}_{+}(r,\varphi ,z=0){\mathrm{\Psi }}_0^0(r,\varphi ,z)$, which allows to consider the intrinsic orbital angular momentum with the vanishing effective thickness of the vortex T → 0. It is clear that the vortex lens successfully converted the Gaussian beam to the left vortex. It is straightforward to apply our analytic formulas for more complex mode with n ≠ 0 and m ≠ 0.

5.9 Extension and limitation of the theory

So far, we have discussed applications of canonical orbital angular momentum operators to Laguerre-Gauss modes in a GRIN fibre, and found that the ladder operations worked properly and analytically calculated the magnitude of orbital angular momentum. This is consistent with a view that orbital angular momentum is a well-defined quantum observable at least in a GRIN fibre. However, if orbital angular momentum is properly defined only in a GRIN fibre, applications of the present theory is limited. Finally, we discuss about the possibility to extend our theory for a more generic case and its limitation.

The reason why we have considered the GRIN fibre is that we can treat Laguerre-Gauss modes virtually exactly within the paraxial approximation. This allows us to consider the collimated beam, propagating in the GRIN fibre, for the fixed waist, w₀. The GRIN fibre contains a limit of a uniform material at g = 0, where the most important application would be the vacuum at n₀ = 1. In the vacuum, however, we need to consider the propagation dependent waist of w = w(z) and the spherical beam radius of R = R(z), as we derived in Section 2.4. It is also important to consider the Gouy phase [5, 42, 46–50]. These factors will add extra contributions during the applications of ladder operations. Fortunately, these extra contributions are negligible for the sufficiently large confocal parameter, known also as Rayleigh length, $z_0=\pi n_0{w}_0^2/\lambda \gg \lambda$. For example, if w₀ = 1.0 mm, z₀ = 2.0 m for λ = 1.55 μm, such that the present work is mostly varied for experiments on an optical bench. In other words, the present theory is applicable for a sufficiently collimated beam. For the diffraction limited narrow beam, z₀ ∼ w₀ ∼ λ, a direction application of our theory must be carefully considered. For the intermediate case, z₀ > w₀ > λ, we can apply our theory near the beam waist, z ≪ z₀. Usually, the thickness (T) of the region, where the refractive index is changed locally in a vortex lens, is of the order of λ, and such that we can consider T ≪ z₀, which justifies the increment and the decrement of orbital angular momentum upon applications of a vortex lens.

Finally, we would like to discuss briefly whether our theory can applicable to other families of structured light, such as Hermite-Gauss [5, 9], Bessel-Gauss [74–77], and Ince-Gauss [78] beams. The Bessel-Gauss beams are quite attractive, since it prevents the diffraction and the mode shape is preserved upon propagation [74–77]. It is beyond the scope of the present work to develop ladder operations to these special functions in general rather than Laguerre-Gauss beams, but Laguerre-Gauss beams have correlation to Hermite-Gauss beams [5, 9]. The Hemite-Gauss beams in the GRIN fibre is given by

${\mathrm{\Phi }}_{l,m}\left(r,\varphi ,z\right)=⟨r,\varphi ,z|{\mathrm{\Phi }}_{l,m}⟩(257)$

$\begin{array}{ccc}\hfill & \hfill & =\frac{1}{w_0}\sqrt{\frac{1}{\pi 2^{l+m-1}l!m!}}{H}_l\left(\frac{\sqrt{2}x}{w_0}\right)H_m\left(\frac{\sqrt{2}y}{w_0}\right)\hfill \\ \hfill & \hfill & {\mathrm{e}}^{-\frac{r^2}{w_0^2}}{\mathrm{e}}^{ikz},\hfill \end{array}(258)$

where H_l(x) and H_m(y) are Hermite functions with l and m as quantum numbers for horizontal and vertical directions [5, 9]. We are interested in the relationship to the Laguerre-Gauss beams ${\mathrm{\Psi }}_n^m(r,\varphi ,z)=⟨r,\varphi ,z|{\mathrm{\Psi }}_n^m⟩$ for the fundamental mode, and we use the formulas of $L_0^m(x)=1$, H₀(x) = 1, and H₁(x) = 2x. As we have seen above, the Laguerre-Gauss beam of $|{\mathrm{\Psi }}_0^m\rangle$ gives the expectation value of $⟨{\widehat{L}}_z⟩=\hslash mN$. This means that the direction of the propagation becomes the quantisation axis for orbital angular momentum and the left vortexed state, $|\mathrm{L},m\rangle =|{\mathrm{\Psi }}_0^m\rangle$, carries orbital angular momentum of ℏm per a photon. We also realise that the mode is degenerate with another mode with the opposite chirality of the right vortexed state, $|\mathrm{R},m\rangle =|{\mathrm{\Psi }}_0^{-m}\rangle$. Then, according to the quantum mechanical superposition principle, we consider the superposition state between orthogonal states of |L, m⟩ and ‖R, m⟩, which forms the SU(2) states [2, 3, 9, 9, 10, 13, 36–38, 79–88]. The generators of rotation in $\mathfrak{s}u(2)$ Lie algebra [89–92] become

${\widehat{\ell }}_1=\hslash m{\sigma }_1(259)$

${\widehat{\ell }}_2=\hslash m{\sigma }_2(260)$

${\widehat{\ell }}_3=\hslash m{\sigma }_3,(261)$

where σ_i (i = 1, 2, 3) is the Pauli matrix, and the commutation relationship becomes

$\left[{\widehat{\ell }}_i,{\widehat{\ell }}_j\right]=2\hslash m{\displaystyle \sum _{k=1}^3}{ϵ}_{\mathit{ijk}}{\widehat{\ell }}_k,(262)$

where ϵ_ijk is a complete anti-symmetric tensor. Then, we obtain the horizontal state in orbital angular momentum state as

$\begin{array}{ccc}\hfill & \hfill & \langle r,\varphi ,z|\mathrm{H},m\rangle \hfill \\ \hfill & \hfill & =\frac{1}{\sqrt{2}}\left(\langle r,\varphi ,z|{\mathrm{\Psi }}_0^m\rangle +\langle r,\varphi ,z|{\mathrm{\Psi }}_0^{-m}\rangle \right)\hfill \end{array}(263)$

$=\langle r,\varphi ,z|{\mathrm{\Phi }}_{\mathrm{1,0}}\rangle ,(264)$

where we have used the formula of r (e^imϕ + e^−imϕ) = 2r cos(ϕ) = 2x. This means that the Laguerre-Gauss beam is converted into the Hermite-Gauss beam [5]. Similarly, the vertical state in orbital angular momentum becomes

$\begin{array}{ccc}\hfill & \hfill & \langle r,\varphi ,z|\mathrm{V},m\rangle \hfill \\ \hfill & \hfill & =\frac{1}{\sqrt{2}}\left(-\langle r,\varphi ,z|{\mathrm{\Psi }}_0^m\rangle +\langle r,\varphi ,z|{\mathrm{\Psi }}_0^{-m}\rangle \right)\hfill \end{array}(265)$

$=-i\langle r,\varphi ,z|{\mathrm{\Phi }}_{\mathrm{0,1}}\rangle ,(266)$

where we have used the formula of r (−e^imϕ + e^−imϕ) = −2ir sin(ϕ) = −i2y. Therefore, both horizontal and vertical states in orbital angular momentum are described by Hermite-Gauss beams [5]. We can calculate the average orbital angular momentum per a photon as $({\ell }_1,{\ell }_2,{\ell }_3)=(⟨{\widehat{\ell }}_1⟩,⟨{\widehat{\ell }}_2⟩,⟨{\widehat{\ell }}_3⟩)$, which can be shown on the higher-order Poincaré sphere [59–61, 69, 93], as shown in Figure 3. The left and right Laguerre-Gauss beams with topological charge of m correspond to the north pole at ℏ(0, 0, m) and the south pole at ℏ(0, 0, − m), respectively. The superposition state of these states form the higher-order Poincaré sphere. The horizontal and vertical states in orbital angular momentum correspond to the points at ℏ(m, 0, 0) and ℏ(−m, 0, 0), respectively. The mode profiles of these beams are horizontal and vertical dipoles, respectively, and the superposition state among these states by a rotator along the ℓ₃ axis correspond to rotate the direction of the dipole. Therefore, the diagonal and anti-diagonal states can also be described by the Hermite-Gauss beams, such that the states along the equator in the higher-order Poincaré sphere can be described by the Hermite-Gauss beams. The ladder operators ${\widehat{l}}_{\pm }$ corresponds to increment and decrement the radius of ± ℏ to move among spheres with different radius. For coherent photons with N photons, the radius is multiplied with N, which increases upon increasing the power.

FIGURE 3

FIGURE 3. The higher-order Poincaré sphere for orbital angular momentum. By considering the superposition states between the left and the right vortexed states with topological charge of m and −m, respectively, the expectation values of the orbital angular momentum, (ℓ₁, ℓ₂, ℓ₃), are calculated, which form the sphere with the radius of ℏm per a photon. (A) The amplitude and (B) the phase of superposition states are calculated at typical points. The ladder operators of ${\widehat{l}}_{\pm }$ change the radius of the sphere by incrementing and decrementing m.

6 Conclusion

A photon, an elementary particle with the internal spin degree of freedom, can have an orbital degree of freedom [5, 6, 8–13]. In a vacuum, a photon travels at the speed of light, c, and is described by a plane wave [10]. The many-body state for photons can be described by a QED theory [1, 4, 11, 12, 44]. On the other hand, for a photon confined to a region with a larger refractive index, i.e., a waveguide, the mode is described as a confined mode, which is a bound state with a discrete energy level. The nature of this mode is completely different from that described by a plane wave allowing a continuous energy spectrum. We have shown that the fundamental equation to describe the orbital wavefunction of photons in a waveguide is Helmholtz equation [9, 10] for a monochromatic coherent beam emitted from a laser, where the spin degree of freedom is described by a Jones vector. We must solve the Helmholtz equation in a material including the refractive index profile and the symmetry of the system [88]. The reason why the many-body photonic state can be described by a single wavefunction is based on the Bose-Einstein condensation nature of a coherent state that allows a macroscopic number of photons to occupy the lowest energy state, because of the Bose statistics due to the integer spin. As a specific example, we have considered a GRIN fibre with a cylindrical symmetry, for which we could solve the Helmholtz equation exactly by using LG modes [5, 6, 8–12].

We have defined canonical OAM operators in cylindrical coordinates and have applied them to the LG modes. We have confirmed that the OAM is quantised along the direction of propagation and that the quantum-mechanical expectation value is indeed obtained as ℏm, while the average values along the directions perpendicular to the propagation vanish. We have found that the ladder operators to increase or decrease m work successfully to increment or decrement in units of ℏ. We could also calculate the quantum-mechanical average of the magnitude of OAM as a function of the radial quantum numbers of n and m. We have also confirmed the contributions from the intrinsic OAM and the origin-dependent extrinsic OAM. Finally, we have calculated the matrix elements of the ladder operators and have confirmed that the angular momentum operators are observable at least in the Hilbert space spanned by the LG modes. From those results, we conclude that the OAM is a proper quantum-mechanical degree of freedom and that a standard quantum-mechanical treatment is applicable to a monochromatic coherent beam of photons.

Data availability statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Author contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Funding

This work is supported by JSPS KAKENHI Grant Number JP 18K19958.

Acknowledgments

The author would like to express sincere thanks to Prof I. Tomita for continuous discussions and encouragements.

Conflict of interest

Author SS was employed by Hitachi, Ltd.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Dirac PAM. The principle of quantum mechanics. Oxford: Oxford University Press (1930).

Google Scholar

2. Baym G. Lectures on quantum mechanics. New York: Westview Press (1969).

Google Scholar

3. Sakurai JJ, Napolitano JJ. Modern quantum mechanics. Edinburgh: Pearson (2014).

Google Scholar

4. Sakurai JJ. Advanced quantum mechanics. New York: Addison-Wesley Publishing Company (1967).

Google Scholar

5. Allen L, Beijersbergen MW, Spreeuw RJC, Woerdman JP. Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes. Phys Rev A (1992) 45:8185–9. doi:10.1103/PhysRevA.45.8185

PubMed Abstract | CrossRef Full Text | Google Scholar

6. v Enk SJ, Nienhuis G. Commutation rules and eigenvalues of spin and orbital angular momentum of radiation fields. J Mod Opt (1994) 41:963–77. doi:10.1080/09500349414550911

CrossRef Full Text | Google Scholar

7. Leader E, Lorcé C. The angular momentum controversy: What’s it all about and does it matter? Phys Rep (2014) 541:163–248. doi:10.1016/j.physrep.2014.02.010

CrossRef Full Text | Google Scholar

8. Barnett SM, Allen L, Cameron RP, Gilson CR, Padgett MJ, Speirits FC, et al. On the natures of the spin and orbital parts of optical angular momentum. J Opt (2016) 18:064004. doi:10.1088/2040-8978/18/6/064004

CrossRef Full Text | Google Scholar

9. Yariv Y, Yeh P. Photonics: Optical electronics in modern communications. Oxford: Oxford University Press (1997).

Google Scholar

10. Jackson JD. Classical electrodynamics. New York: John Wiley and Sons (1999).

Google Scholar

11. Grynberg G, Aspect A, Fabre C. Introduction to quantum optics: From the semi-classical approach to quantized light. Cambridge: Cambridge University Press (2010).

Google Scholar

12. Bliokh KY, Rodríguez-Fortuño FJ, Nori F, Zayats AV. Spin-orbit interactions of light. Nat Photon (2015) 9:796–808. doi:10.1038/NPHOTON.2015.201

CrossRef Full Text | Google Scholar

13. Al-Attili AZ, Burt D, Li Z, Higashitarumizu N, Gardes F, Ishikawa Y, et al. Chiral germanium micro-gears for tuning orbital angular momentum. Sci Rep (2022) 12:7465. doi:10.1038/s41598-022-11245-1

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Allen L, Padgett MJ. The poynting vector in Laguerre-Gaussian beams and the interpretation of their angular momentum density. Opt Comm (2000) 184:67–71. doi:10.1016/S0030-4018(00)00960-3

CrossRef Full Text | Google Scholar

15. Ji X. Comment on “Spin and orbital angular momentum in gauge theories: Nucleon spin structure and multipole radiation revisited”. Phys Rev Lett (2010) 104:039101. doi:10.1103/PhysRevLett.104.039101

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Chen XS, Lü XF, Sun WM, Wang F, Goldman T. Spin and orbital angular momentum in gauge theories: Nucleon spin structure and multipole radiation revisited. Phys Rev Lett (2008) 100:232002. doi:10.1103/PhysRevLett.100.232002

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Yang LP, Khosravi F, Jacob Z. Quantum field theory for spin operator of the photon. Phys Rev Res (2022) 4:023165. doi:10.1103/PhysRevResearch.4.023165

CrossRef Full Text | Google Scholar

18. Forbes A, d Oliveira M, Dennis MR. Structured light. Nat Photon (2021) 15:253–62. doi:10.1038/s41566-021-00780-4

CrossRef Full Text | Google Scholar

19. Nape I, Sephton B, Ornelas P, Moodley C, Forbes A. Quantum structured light in high dimensions. APL Photon (2023) 8:051101. doi:10.1063/5.0138224

CrossRef Full Text | Google Scholar

20. Ma M, Lian Y, Wang Y, Lu Z. Generation, transmission and application of orbital angular momentum in optical fiber: A review. Front Phys (2021) 9:773505. doi:10.3389/fphy.2021.773505

CrossRef Full Text | Google Scholar

21. Rosen GFQ, Tamborenea PI, Kuhn T. Interplay between optical vortices and condensed matter. Rev Mod Phys (2022) 94:035003. doi:10.1103/RevModPhys.94.035003

CrossRef Full Text | Google Scholar

22. Shen Y, Rosales-Guzmán C. Nonseparable states of light: From quantum to classical. Laser Photon Rev (2022) 16:2100533. doi:10.1002/lpor.202100533

CrossRef Full Text | Google Scholar

23. Cisowski C, Götte JB, Franke-Arnold S. Colloquium: Geometric phases of light: Insights from fiber bundle theory. Rev Mod Phys (2022) 94:031001. doi:10.1103/revmodphys.94.031001

CrossRef Full Text | Google Scholar

24. Liu R, Zhang J, Liu J, Lin Z, Li Z, Lin Z, et al. 1-pbps orbital angular momentum fibre-optic transmission. Light Sci Appl (2022) 11:202. doi:10.1038/s41377-022-00889-3

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Spreeuw BJC. A classical analogy of entanglement. Found Phys (1998) 28:361–74. doi:10.1023/A:1018703709245

CrossRef Full Text | Google Scholar

26. Shen Y. Rays, waves, su(2) symmetry and geometry: Toolkits for structured light. J Opt (2021) 23:124004. doi:10.1088/2040-8986/ac3676

CrossRef Full Text | Google Scholar

27. Saito S. Poincaré rotator for vortexed photons. Front Phys (2021) 9:646228. doi:10.3389/fphy.2021.646228

CrossRef Full Text | Google Scholar

28. Saito S. Spin of photons: Nature of polarisation (2023). arXiv (2023) 2303.17112. doi:10.48550/arXiv.2303.17112

CrossRef Full Text | Google Scholar

29. Saito S. SU(2) symmetry of coherent photons and application to poincaré rotator (2023). arXiv (2023) 2303.18199. doi:10.48550/arXiv.2303.18199

CrossRef Full Text | Google Scholar

30. Saito S. Macroscopic singlet, triplet, and colour-charged states of coherent photons (2023). arXiv (2023) 2304.01216. doi:10.48550/arXiv.2304.01216

CrossRef Full Text | Google Scholar

31. Chuang SL. Physics of photonic devices. New York: John Wiley and SonsWiley (2009).

Google Scholar

32. Goldstein DH. Polarized light. London: CRC Press (2011). doi:10.1201/b10436

CrossRef Full Text | Google Scholar

33. Gil JJ, Ossikovski R. Polarized light and the mueller matrix approach. London: CRC Press (2016). doi:10.1201/b19711

CrossRef Full Text | Google Scholar

34. Kawakami S, Nishizawa J. An optical waveguide with the optimum distribution of the refractive index with reference to waveform distortion. IEEE Trans Microw Theor Techn. (1968) 16:814–8. doi:10.1109/TMTT.1968.1126797

CrossRef Full Text | Google Scholar

35. Joannopoulos JD, Johnson SG, Winn JN, Meade RD. Photonic crystals: Molding the flow og light. New York: Princeton Univ. Press (2008).

Google Scholar

36. Sotto M, Tomita I, Debnath K, Saito S. Polarization rotation and mode splitting in photonic crystal line-defect waveguides. Front Phys (2018) 6:85. doi:10.3389/fphy.2018.00085

CrossRef Full Text | Google Scholar

37. Sotto M, Debnath K, Khokhar AZ, Tomita I, Thomson D, Saito S. Anomalous zero-group-velocity photonic bonding states with local chirality. J Opt Soc Am B (2018) 35:2356–63. doi:10.1364/JOSAB.35.002356

CrossRef Full Text | Google Scholar

38. Sotto M, Debnath K, Tomita I, Saito S. Spin-orbit coupling of light in photonic crystal waveguides. Phys Rev A (2019) 99:053845. doi:10.1103/PhysRevA.99.053845

CrossRef Full Text | Google Scholar

39. Bliokh KY, Bekshaev AY, Nori F. Optical momentum, spin, and angular momentum in dispersive media. Phys Rev Lett (2017) 119:073901. doi:10.1103/PhysRevLett.119.073901

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bliokh KY, Bekshaev AY, Nori F. Optical momentum and angular momentum in complex media: From the abraham-minkowski debate to unusual properties of surface plasmon-polaritons. New J Phys (2017) 19:123014. doi:10.1088/1367-2630/aa8913

CrossRef Full Text | Google Scholar

41. Barnett SM, Babiker M, Padgett MJ. Optical orbital angular momentum. Phil Trans R Soc A (2016) 375:20150444. doi:10.1098/rsta.2015.0444

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Simon R, Mukunda N. Bargmann invariant and the geometry of the Güoy effect. Phys Rev Lett (1993) 70:880–3. doi:10.1103/PhysRevLett.70.880

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Lloyd SM, Babiker M, Thirunavukkarasu G, Yuan Y. Electron vortices: Beams with orbital angular momentum. Rev Mod Phys (2017) 89:035004. doi:10.1103/RevModPhys.89.035004

CrossRef Full Text | Google Scholar

44. Fox M. Quantum optics: An introduction. Oxford: Oxford University Press (2006).

Google Scholar

45. Schrieffer JR. Theory of superconductivity. Boca Raton: CRC Press (1971). doi:10.1201/9780429495700

CrossRef Full Text | Google Scholar

46. Pancharatnam S. Generalized theory of interference, and its applications. Proc Indian Acad Sci Sect A (1956) XLIV:247–62. doi:10.1007/BF03046050

CrossRef Full Text | Google Scholar

47. Berry MV. Quantual phase factors accompanying adiabatic changes. Proc R Sco Lond A (1984) 392:45–57. doi:10.1098/rspa.1984.0023

CrossRef Full Text | Google Scholar

48. Tomita A, Cao RY. Observation of Berry’s topological phase by use of an optical fiber. Phys Rev Lett (1986) 57:937–40. doi:10.1103/PhysRevLett.57.937

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Hamazaki J, Mineta Y, Oka K, Morita R. Direct observation of gouy phase shift in a propagating optical vortex. Opt Exp (2006) 14:8382–92. doi:10.1364/OE.14.008382

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Bliokh K. Geometrodynamics of polarized light: Berry phase and spin Hall effect in a gradient-index medium. J Opt A: Pure Appl Opt (2009) 11:094009. doi:10.1088/1464-4258/11/9/094009

CrossRef Full Text | Google Scholar

51. Nambu Y. Quasi-particles and gauge invariance in the theory of superconductivity. Phys Rev (1960) 117:648–63. doi:10.1103/PhysRev.117.648

CrossRef Full Text | Google Scholar

52. Anderson PW. Random-phase approximation in the theory of superconductivity. Phys Rev (1958) 112:1900–16. doi:10.1103/PhysRev.112.1900

CrossRef Full Text | Google Scholar

53. Goldstone J, Salam A, Weinberg S. Broken symmetries. Phy Rev (1962) 127:965–70. doi:10.1103/PhysRev.127.965

CrossRef Full Text | Google Scholar

54. Higgs PW. Broken symmetries and the masses of gauge bosons. Phys Lett (1962) 12:508–9. doi:10.1103/PhysRevLett.13.508

CrossRef Full Text | Google Scholar

55. Saito S. Dirac equation for photons: Origin of polarisation (2023). arXiv (2023) 2303.18196. doi:10.48550/arXiv.2303.18196

CrossRef Full Text | Google Scholar

56. Arfken GB, Weber HJ. Mathematical methods for physicists. London: Academic Press (2005).

Google Scholar

57. Whittaker ET, Watson GN. A course of modern analysis. Cambridge: Cambridge University Press (1962). doi:10.1017/9781009004091

CrossRef Full Text | Google Scholar

58. Bateman H. Higher transcendental functions [volumes I-III]. New York: MacGraw-Hill Book Company (1953).

Google Scholar

59. Padgett MJ, Courtial J. Poincaré-sphere equivalent for light beams containing orbital angular momentum. Opt Lett (1999) 24:430–2. doi:10.1364/OL.24.000430

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Milione G, Sztul HI, Nolan DA, Alfano RR. Higher-order poincaré sphere, Stokes parameters, and the angular momentum of light. Phys Rev Lett (2011) 107:053601. doi:10.1103/PhysRevLett.107.053601

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Liu Z, Liu Y, Ke Y, Liu Y, Shu W, Luo H, et al. Generation of arbitrary vector vortex beams on hybrid-order Poincaré sphere. Photon Res (2017) 5:15–21. doi:10.1364/PRJ.5.000015

CrossRef Full Text | Google Scholar

62. Marrucci L, Manzo C, Paparo D. Optical spin-to-orbital angular momentum conversion in inhomogeneous anisotropic media. Phys Rev Lett (2006) 96:163905. doi:10.1103/PhysRevLett.96.163905

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Machavariani G, Lumer Y, Moshe I, Meir A, Jackel S. Efficient extracavity generation of radially and azimuthally polarized beams. Opt Lett (2007) 32:1468–70. doi:10.1364/OL.32.001468

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Lai WJ, Lim BC, Phua PB, Tiaw KS, Teo HH, Hong MH. Generation of radially polarized beam with a segmented spiral varying retarder. Opt Exp (2008) 16:15694–9. doi:10.1364/OE.16.015694

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Guan B, Scott RP, Qin C, Fontaine NK, Su T, Ferrari C, et al. Free-space coherent optical communication with orbital angular, momentum multiplexing/demultiplexing using a hybrid 3d photonic integrated circuit. Opt Exp (2013) 22:145–56. doi:10.1364/OE.22.000145

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Sun J, Moresco M, Leake G, Coolbaugh D, Watts MR. Generating and identifying optical orbital angular momentum with silicon photonic circuits. Opt Lett (2014) 39:5977–80. doi:10.1364/OL.39.005977

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Naidoo D, Roux FS, Dudley A, Litvin I, Piccirillo B, Marrucci L, et al. Controlled generation of higher-order Poincaré sphere beams from a laser. Nat Photon (2016) 10:327–32. doi:10.1038/NPHOTON.2016.37

CrossRef Full Text | Google Scholar

68. Dorney KM, Rego L, Brooks NJ, San Román J, Liao CT, Ellis JL, et al. Controlling the polarization and vortex charge of attosecond high-harmonic beams via simultaneous spin-orbit momentum conservation. Nat Photon (2019) 13:123–30. doi:10.1038/s41566-018-0304-3

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Erhard M, Fickler R, Krenn M, Zeilinger A. Twisted photons: New quantum perspectives in high dimensions. Light: Sci Appl (2018) 7:17146. doi:10.1038/lsa.2017.146

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Andrews DL. Symmetry and quantum features in optical vortices. Symmetry (2021) 13:1368. doi:10.3390/sym.13081368

CrossRef Full Text | Google Scholar

71. Angelsky OV, Bekshaev AY, Dragan GS, Maksimyak PP, Zenkova CY, Zheng J. Structured light control and diagnostics using optical crystals. Front Phys (2021) 9:715045. doi:10.3389/fphy.2021.715045

CrossRef Full Text | Google Scholar

72. Agarwal GS. SU(2) structure of the poincaré sphere for light beams with orbital angular momentum. J Opt Soc A A (1999) 16:2914–6. doi:10.1364/JOSAA.16.002914

CrossRef Full Text | Google Scholar

73. Golub MA, Shimshi L, Davidson N, Friesem AA. Mode-matched phase diffractive optical element for detecting laser modes with spiral phases. Appl Opt (2007) 46:7823–8. doi:10.1364/AO.46.007823

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Gori F, Guattari G, Padovani C. Bessel gauss beams. Opt Commun (1987) 64:491–5. doi:10.1016/0030-4018(87)90276-8

CrossRef Full Text | Google Scholar

75. Durnin J, Miceli JJJ, Eberly JH. Diffraction-free beams. Phys Rev Lett (1987) 58:1499–501. doi:10.1103/physrevlett.58.1499

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Wang W, Ye T, Wu Z. Probability property of orbital angular momentum distortion in turbulence. Opt Exp (2021) 29:44157–73. doi:10.1364/OE.445175

CrossRef Full Text | Google Scholar

77. Wang W, Zhang G, Ye T, Wu Z, Bai L. Scintillation of the orbital angular momentum of a bessel Gaussian beam and its application on multi-parameter multiplexing. Opt Exp (2023) 31:4507–20. doi:10.1364/OE.478127

CrossRef Full Text | Google Scholar

78. Bandres MA, Gutiérrez-Vega JC. Ince-Gaussian beams. Opt Lett (2004) 29:144–6. doi:10.1364/OL.29.000144

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Jones RC. A new calculus for the treatment of optical systems i. description and discussion of the calculus. J Opt Soc Am (1941) 31:488–93. doi:10.1364/JOSA.31.000488

CrossRef Full Text | Google Scholar

80. Payne WT. Elementary spinor theory. Am J Phys (1952) 20:253–62. doi:10.1119/1.1933190

CrossRef Full Text | Google Scholar

81. Born M, Wolf E. Principles of optics. Cambridge: Cambridge University Press (1999). doi:10.1017/9781108769914

CrossRef Full Text | Google Scholar

82. Collett E. Stokes parameters for quantum systems. Am J Phys (1970) 38:563–74. doi:10.1119/1.1976407

CrossRef Full Text | Google Scholar

83. Luis A. Degree of polarization in quantum optics. Phys Rev A (2002) 66:013806. doi:10.1103/PhysRevA.66.013806

CrossRef Full Text | Google Scholar

84. Luis A. Polarization distributions and degree of polarization for quantum Gaussian light fields. Opt Comm (2007) 273:173–81. doi:10.1016/j.optcom.2007.01.016

CrossRef Full Text | Google Scholar

85. Björk G, Söderholm J, Sánchez-Soto LL, Klimov AB, Ghiu I, Marian P, et al. Quantum degrees of polarization. Opt Comm (2010) 283:4440–7. doi:10.1016/j.optcom.2010.04.088

CrossRef Full Text | Google Scholar

86. d Castillo Gft , García IR. The Jones vector as a spinor and its representation on the Poincaré sphere. Rev Mex Fis (2011) 57:406–13. doi:10.48550/arXiv.1303.4496

CrossRef Full Text | Google Scholar

87. Al-Attili AZ, Kako S, Husain MK, Gardes FY, Higashitarumizu N, Iwamoto S, et al. Whispering gallery mode resonances from ge micro-disks on suspended beams. Front Mat (2015) 2:43. doi:10.3389/fmats.2015.00043

CrossRef Full Text | Google Scholar

88. Saito S, Tomita I, Sotto M, Debnath K, Byers J, Al-Attili AZ, et al. Si photonic waveguides with broken symmetries: Applications from modulators to quantum simulations. Jpn J Appl Phys (2020) 59:SO0801. doi:10.35848/1347-4065/ab85ad

CrossRef Full Text | Google Scholar

89. Fulton W, Harris J. Representation theory: A first course. New York: Springer (2004).

Google Scholar

90. Hall BC. Lie groups, Lie algebras, and representations; an elementary introduction. Switzerland: Springer (2003).

Google Scholar

91. Pfeifer W. The Lie Algebras su(N) An Introduction. Berlin: Springer Basel AG (2003).

Google Scholar

92. Georgi H. Lie algebras in particle physics: From isospin to unified theories (Frontiers in physics). Massachusetts: Westview Press (1999).

Google Scholar

93. Holleczek A, Aiello A, Gabriel C, Marquardt C, Leuchs G. Classical and quantum properties of cylindrically polarized states of light. Opt Exp (2011) 19:9714–36. doi:10.1364/OE.19.009714

PubMed Abstract | CrossRef Full Text | Google Scholar