Quantum Commutation Relationship for Photonic Orbital Angular Momentum

Orbital Angular Momentum (OAM) of photons are already ubiquitously being 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 a cylindrical coordinate, 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 decrement in units of the Dirac constant ($\hbar$). 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.

Orbital Angular Momentum (OAM) of photons are already ubiquitously being 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 a cylindrical coordinate, 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 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.

I. INTRODUCTION
Quantum commutation relationship between operators is an indispensable characteristic in connection with measurements of physical observables [1][2][3][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][2][3][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][2][3][4], and commutation relationship is used to rotate the states.The spin, Ŝ, with one-half in units of the Dirac constant, = h/(2π), where h it the Plank constant, for an electron is perfectly described with Pauli's spin matrices [1][2][3][4].The Orbital Angular Momentum (OAM), 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 m l , satisfies L2 |l, m = l(l +1) 2 |l, m , with an integer l associated with the magnitude of angular momentum [1][2][3][4].Furthermore, the total angular momentum, Ĵ = Ŝ+ 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][2][3][4].
However, this elegant theoretical framework of angular momentum is less trivial in applying to a photon [5][6][7][8][9][10][11][12], 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 propa-gation direction of a photon sets a natural quantisation axis; in the propagation direction, the angular momentum operator, Lz , is successfully obtained by the classical analogue using the Poynting vector [5][6][7][8][11][12][13].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][7][8].It was argued that the spin and OAM operators would commute [6].It is now generally believed that the total angular momentum operator Ĵ for photons is welldefined, while it is impossible to split into Ŝ and L in a local gauge-invariant way [7,8,14], despite a big challenge against this claim [15].This is outside the scope of this paper, and we will revisit this issue in a forthcoming paper.
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, 16] 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 could the spin of a photon be regarded as a classical degree of freedom, which is described by a commutable operator?It should be treated with 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 [17,18]?Detailed discussion about this will be provided in a separate paper.(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 a cylindrical coordinate [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 momen-tum 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 L± = Lx ± i Ly , and apply them to the LG modes?Can we show the increment or 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 Sec.II, we review 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,14,15] in this paper.We are considering monochromatic coherent light sources from lasers, such that the incoherent unpolarised lights will not be considered, either.We have also included detailed appendices to define associate Laguerre functions and have shown various mathematical formulas to make this paper self-contained.In the calculations, well-known formulas and equations appear; yet the others are newly derived, particularly in evaluating the ladder operations.We hope the appendices will help readers, because definitions depend on literatures for factors and signs.
In Sec.III, 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 Sec.IV, 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 Sec.V, we conclude that OAM is indeed a genuine quantum-mechanical observable at least in a graded index (GRIN) fibre [9,19] satisfying some conditions.

A. Maxwell's equations
We start from Maxwell's equations [9,10], in a non-magnetic transparent material of the dielectric constant ǫ and the permittivity µ 0 without charges and currents.As usual for describing electromagnetic fields [9,10], we use complex oscillation fields in SI units for electric field E (V/m), displacement D (C/m 2 ), magnetic field H (A/m), and induction B (Wb/m 2 ), with materials equations D = ǫE and B = µ 0 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 E = (E + E * )/2, after the calculation using a complex field of E. While µ 0 is approximately the same as that in a vacuum, ǫ is different in a material from the value of ǫ 0 in a 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, which is valid in a perfectly uniform material and in a 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 approximations are the sufficient uniformity (∇ǫ ≃ 0).Therefore our analysis will not be valid if ǫ is significantly changed in nano-metre-scale such as for photonic crystals [20][21][22][23] and other inhomogeneous systems [24,25].

B. 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 E(x, y, z) = E 0 ψ(x, y, z)e i(kz−ωt) n, (6) where E 0 is the magnitude of the electric field, k is the wavenumber in the material, ω is the angular frequency, the unit vector 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 E into the Helmholtz equation, we obtain 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, and obtain [26] iλ where we have used the dispersion relationship, ω = vk, with the velocity v = 1/ √ µ 0 ǫ = c/n of a photon in a material of a refractive index of n.The velocity of a photon in a vacuum is c = 1/ √ µ 0 ǫ 0 and λ = λ/(2π) is an angular wavelength.Equation ( 9) is exactly the same form with a standard non-relativistic Schrödinger equation [2,3,16,26,27] for a particle of mass, m, in a 2D xy-plane at time t.The correspondence is summarised in Table I.This implies the propagation of a photon along a paraxial optical path can be described by the same mechanism with the dynamics of a massive quantummechanical particle [26].In fact, electron vortices similar to photonic ones were observed [28], and essentially the same mathematical and physical techniques are 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 of large n similar to low velocity of a heavy particle.

C. 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,29].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 [30].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,29].Specifically, a coherent state for σ-polarisation in a uniform material is described by where σ = H for horizontally polarised state and σ = V for vertically polarised state.α σ is a complex number.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].â † σ and âσ are creation and annihilation operators, satisfying Bose commutation relationship [4,11,29] 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 We can also confirm that |α σ is normalised as α σ |α σ = 1.Then, we can calculate the average number of photons for both polarised states as 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 ≤ α ≤ π 2 ) to describe the polarisation.Alternatively, α σ is determined by the polarisation state as where δ ∈ (0, 2π) is the phase of the polarisation.The overall coherent state is described by the direct product as 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 where V is the volume, x and ŷ are unit vectors along x and y, and β = kz − ωt + β 0 describes the trivial time and space evolution with a global U (1) phase of β 0 .If we apply Ê(z, t) to the coherent state |α H , α V as, we realise the state is an eigenstate of Ê(z, t) with the complex eigenvalue of E(z, t) = E 0 e iβ cos αx + e iδ sin αŷ , where E 0 = 2 ωN/(ǫV ).This means that the coherent state is a simultaneous eigenstate for diagonalising Ê(z, t) and âσ .Alternatively, we can consider the coherent state, Ê(z, t)|α H , α V , describes the photonic state of the system, since the multiplication of the operator does not change the state.In fact, E(z, t) actually corresponds to the spinor description of the polarisation state as 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 0 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 where the scalar complex electric field, E(r, t) = E 0 ψ(r)e iβ , describes the orbital part.The state, Ê(r, t)|α H , α V , describes the entire photonic state, including polarisation.In this case, E(r, t) becomes where the orbital part is described by Ψ(r) = ψ(r)e iβ .Ψ(r) is determined by the scalar Helmholtz equation 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 states 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 If we want to calculate the real electric field, instead of the complex field, we should use the electric field operator defined by which is an observable, such that we can calculate the expectation value, Ê , quantum-mechanically, using |N, α, δ .
D. 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,13].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 a cylindrical coordinate 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 a cylin- where r = x 2 + y 2 is the radius and φ = tan −1 (y/x) is the azimuthal angle.We will solve this equation by using a trial wavefunction 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 The Gaussian mode solution is given by assuming which will give us q(z) = z + q 0 = z − iz 0 and P (z) = i ln (1 + z/q 0 ) with the confocal parameter z 0 = kw 2 0 /2 = πnw 2 0 /λ and the minimum waist of w 0 .We then obtain the Gaussian factor and the phase-shift factor where the beam waist, w(z), the beam radius, R(z), and the phase, η(z), are given by The focal point of the Gaussian beam is set at the origin, where the waist becomes minimum w(0) = w 0 .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 0 , and consequently z 0 , can be controlled and determined by a boundary condition.It is also useful to note that kww ′ = 2z/z 0 and kw 2 /q = 2 (i + z/z 0 ) hold, and we obtain We now focus on the last three terms of this equation, which can be rewritten by exchanging valuables subsequently using ρ = r/w, a = ρ 2 , and b = 2a as where we used the differential equation for the associate Laguerre function, L m p , (Appendices) 40) The rest of the Helmholtz equation is which gives the phase-shift Finally, we obtain which is not normalised, yet.The norm N norm of the wavefunction is obtained by where we used the orthogonality condition (Appendices) Thus, we obtain 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 such that the solution does not depend on the sign of m.
Therefore, the final normalised wavefunction becomes Here, the wavefunction was normalised in the xy-plane as 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 0 (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 = 2 ω/ǫV , but actually a laser beam contains a macroscopic number of photons.By comparing the factors between E 0 and ψ(r, φ, z), we realise that we should assume V = w(z) 2 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 0 should be simply re-defined as E 0 = 2 ωN/ǫL.It is worth making a remark on the Gouy phase [5,27,[31][32][33][34][35] of which is the same as a geometrical phase of Pancharatnam-Berry.This term appears due to the focusing of the beam, which will change (E x , E y ) at z → −∞ to (−E x , −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 0 .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 m l , for which the value of m is limited to be m On the other hand, there is no such restriction to L |m| p , where the so-called LG 01 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.Therefore, 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.E. 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,19], which has the refractive index n(r) dependence given by n(r) 2 = ǫ(r)/ǫ 0 = n 2 0 1 − g 2 r 2 , 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 a vacuum as k 0 = 2π/λ for a laser beam emitted from the waveguide.The energy of the photon will not be changed upon the emission, such that ω = ck 0 is valid, while the dispersion relationship, ω = ω(k), in the waveguide is highly nontrivial, and we will obtain this from the Helmholtz equation.We also define a constant wavenumber parameter, k n0 = 2π/λ n0 = 2πn0 λ = k 0 n 0 , since the waveguide is mostly determined by the core refractive index, n 0 = n(0), and this is simply a constant parameter, where k n0 is different from the true wavenumber, k, responsible for describing photon momentum of p = k.
With those parameters, we can rewrite ), and the Helmholtz equation becomes In the cylindrical coordinate, we can convert this equation to 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 with the constant waist w 0 and the constant complex radius, q, becomes the solution.Again, we tentatively assume m ≥ 0. By inserting Eq. (53) into Eq.( 52), we obtain We then obtain a stable Gaussian form by noting We also use useful identities and the Helmholtz equation then becomes 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 and the remaining equation becomes This provides the solution [9] where we have defined a phase velocity at the core as v 0 = c/n 0 and a frequency shift as δω 0 = v 0 g.By solving Eq. ( 62) with respect to ω, we obtain the dispersion relationship 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 [30] and the Nambu-Anderson-Goldstone-Higgs theory of a broken symmetry [36][37][38][39].We infer that a similar symmetry principle is hidden in our system.We will discuss this in a subsequent paper.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 in a GRIN fibre without the slowly varying paraxial approximation.

A. 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 rays.
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 I).Unlike in a vacuum without a material, where Ψ(r) is a simple planewave, 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.

B. Canonical orbital angular momentum operators
The most well-established quantum mechanical treatment is canonical commutation relationship for the position r = (x, ŷ, ẑ) and the momentum p = (p x , py , pz ) operators [1][2][3]: [x, px ] = i , [ŷ, py ] = i , [ẑ, pz ] = i , and commutable relationship among other combinations.The OAM operator, l, is defined by use of those r and p as l = r × p, such that each component becomes 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][2][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 l in a cylindrical coordinate.
The unit vectors of a cylindrical coordinate are defined for a rotation in a 3D Cartesian coordinate along the z axis, where r is the unit vector along r and Φ 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., r = r(φ) and Φ = Φ(φ).Therefore, the nabla operator is of the form and the Laplacian becomes where we have used This φ dependence of r(φ) and Φ(φ) makes it difficult to define and treat the OAM operators for lr and lφ .Instead, we will keep using lx and ly , defined above, and we will express them by (r, φ, z).In order to use (r, φ) instead of (x, y), we obtain By inserting these into l, we obtain lx = i −z sin φ∂ r − z r cos φ∂ φ + r sin φ∂ z , (74) We can readily confirm the original commutation relationship lx , ly = i lz (77) and its cyclic exchanges lz , lx = i ly (79) are all valid in the cylindrical coordinate (r, φ, z).
We also obtain raising and lowering operators, respectively, as l+ = lx + i ly (80)

C. Application to plane waves & 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: 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 l, we obtain which means that Ψ is indeed an eigenstate for l z and that l± 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 By averaging these over space, we obtain lx (r, φ, z) = 0, (91) ly (r, φ, z) = 0, (92) Therefore, the expectation values are reasonably welldefined.
However, if we calculate the complex conjugate of l− Ψ simply by taking its complex conjugate as and multiply this by l− Ψ from the right to calculate the norm, we obtain which is a positive real value.On the other hand, the direct calculation of Ψ * l + l− Ψ becomes This implies that l− (r, φ, z) † = l+ (r, φ, z), 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 which also imply We also see 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 l± Y m l (θ, φ).On the other hand, the direct calculation shows such that the commutation relationship l+ (r, φ, z), l− (r, φ, z) = 2 lz (r, φ, z) 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 ikz e imφ .This is confirmed by calculating the magnitude of the OAM along the radial direction as which gives an imaginary part, thus showing that the magnitude is not observable.If we take the average over z ∈ (0, L), we obtain which is still a complex value.Further average over r ∈ (0, R), where R is the radius of a cylindrical waveguide, gives which diverges at the origin.We can also integrate over z ∈ (L/2, L/2), and obtain 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 quantummechanically.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.

IV. 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][9][10]] 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 [40][41][42], 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 2 0 to simplify calculations.We use various formulas for associate Laguerre functions, which are summarised in Appendixes.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, Lz , along z for a coherent state is obtained by multiplying N as Lz = mN .

B. 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.

Rising operation for m > 0
We assume m > 0 and calculate l+ Ψ m n = e i(m+1)φ e ikz 1 w 0 where the last factor becomes where we have used Therefore, we obtain The most significant part of this expression is that we confirm l+ Ψ m n ∝ Ψ m+1 n , 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 l+ .Nevertheless, the main role of l+ to increase m was successfully confirmed for m > 0.

Rising operation for m = 0
We then continue to calculate for m = 0 as where the last factor becomes and thus we obtain 120) which is exactly the same expression with that of m > 0.

Rising operation for m < 0
We obtain a similar result by directly calculating l+ Ψ m n = e i(m+1)φ e ikz 1 w 0 where the last factor becomes where we have used Thus, we obtain Therefore, the raising operator is successfully working to increment m, independently of the value and the sign of m.

Lowering operation for m < 0
Next, we apply the lowering operator to the case for m < 0, and obtain l− Ψ m n = e i(m−1)φ e ikz 1 w 0 where the last factor becomes . (127) This also shows that the lowering operator can actually lower m to m − 1.

Lowering operation for m = 0
Similarly, we obtain for m = 0: where the last factor becomes which yields This is the same formula as that we obtained for m < 0.

Lowering operation for m > 0
Finally, we apply the lowering operator to the case for m > 0, and obtain where the last factor becomes 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 .

C. 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 Appendixes.
Then, we obtain 2. Norm of l+Ψ m n for m < 0 Similarly, we obtain Next, we calculate the norm of l− Ψ m n for m ≤ 0 as

Summary of the norm of l±Ψ m n
The above direct calculations show that we obtain the same norm for l± Ψ m n : 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 ( kw 0 ) 2 /2 = (pw 0 ) 2 /2 with the momentum p = k is intuitively understandable, since pw 0 has the dimension of the angular momentum.As for the case of the plane wave, however, the fact that we obtain 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 l+ l− and l− l+ , directly, using the LG modes, which are shown in the next subsection.

D. Validity of commutation relationship for LG modes
First, we must calculate the wavefunction of l+ l− Ψ m n and l− l+ Ψ m n and calculate the expectation value.This is straightforward but tedious.Again, we will split calculations for positive and negative values of m.

Ladder operators
Before calculating the matrix elements, we will describe operators, l+ l− and l− l+ , by using a cylindrical coordinate (r, φ, z), as and l− l+ Thus, we also confirmed the commutation relationship for l± , l+ , l− = 2 lz . (142) Therefore, if we apply these operators to a LG mode, we must confirm that this identity is always valid.This is useful to check the validity of the calculation.
2. l+ l− operation for m ≥ 0 First, we evaluate various terms as follows: and finally By summing up some of these terms, we obtain By using the integration formulas (Appendixes), we finally obtain 3. l− l+ operation for m ≥ 0 The only source of the difference between l+ l− and l− l+ operations is coming from the sign of i∂ φ .Therefore, it is straightforward to obtain This result, together with the previous result for l− l+ , confirms that the commutation relationship over average indeed satisfies l+ , l− = 2 lz , where lz = m for m ≥ 0. We can proceed for m ≤ 0 in a similar way.First, we evaluate the factors: and the other factors are the same ones for m ≥ 0.Then, we obtain Therefore, we finally obtain 5. l− l+ operation for m ≤ 0 Again, the only source of the change between l+ l− and l− l+ operations is coming from the sign of i∂ φ , such that we obtain Therefore, the commutation relationship over average is also valid for m ≤ 0.

E. Summary of expectation values of l+ l− and l− l+ operations
The above results are summarised as follows, is also valid for ∀ m.

F. Magnitude of OAM
The above calculations have confirmed the quantum commutation relationship, [ l+ , l− ] = 2 lz , 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 l+ l− and l− l+ are useful to evaluate the magnitude of the OAM, because the inner product of the vectorial OAM operator is described as We can confirm that the expectation value is independent of whether we are using l+ l− or l− l+ , and we obtain 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 while if we average over z ∈ (−L/2, L/2), it becomes 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 δω 0 m ≫ v 0 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, 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 lx , ly , 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 where w 1 = w 0 / √ 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 .

G. Transfer matrix element
Finally, we calculate the matrix elements of l± 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 l± .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 for m > 0, and 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 for m ≥ 0, and for m < 0. For both cases, the relationship, is always valid for ∀ m.This implies l † + = l− (181) l † − = l+ (182) 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 [43][44][45][46][47][48][49].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.

V. CONCLUSIONS
A photon, an elementary particle with the internal spin degree of freedom, can have an orbital degree of freedom [5,6,[8][9][10][11][12].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,29].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 ray 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.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][9][10][11][12].
We have defined canonical OAM operators in a cylindrical coordinate and have applied them to the LG modes.We have confirmed that the OAM is quantised along the direction of propagation and that the quantummechanical 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 origindependent 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 ray of photons.for n ≥ 1.

Ladder operators for radial quantum number
For obtaining ladder operators, we calculate the derivative of the generating function by t as However, this expression is not perfect, since the derivative operator remained in the right-hand side, which will be removed later.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.

Ladder operators for orbital angular momentum
The above formulas for raising and lowering the radial quantum numbers are known in literatures [50][51][52], 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 (B37)

Orthogonality relationship
The orthogonality relationship is obtained in a similar way by using the generating function and calculating the sum

First, we have
checked the expectation values of l by use of the LG modes.This was straightforward by noting l± ψ m n ∝ ψ m±1 n ∝ e i(m±1)φ and

1 )
m d m dt m L (n+1)+m (t) straightforward to obtain the lowering operator.As for preparations, we recognised several useful recurrence formulas use the raising operator for n to obtainL m n+1 (t) = 1 n + 1 t d dt − t + n + m + 1 L m n (t).

TABLE I .
Mapping of Helmholtz equation to a nonrelativistic Schrödinger equation.
Next, we construct the raising operator by calculating the derivative of the recurrence equation by t as ≥ 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 m n+1 (t) to obtain the lowering operator