Exactness of Bohr-Sommerfeld quantisation for two non-central potentials
David T. S. Perkins Robert A. Smith

	
	
	
	
	School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, United Kingdom
	
	Exactness of Bohr-Sommerfeld quantisation for two non-central potentials
    David T S Perkins and Robert A Smith
    
========================================================================


        In this paper we demonstrate the integrability of the Hamilton-Jacobi equation for two non-central potentials in spherical polar coordinates, and 
        present complete solutions for the classically bound orbits. We then show that the semiclassical method of Bohr-Sommerfeld quantisation exactly
        reproduces the bound state spectra of the corresponding quantum mechanical Schrödinger equations. One of these potentials has previously been
        analysed in parabolic coordinates; the results for the other are, to the authors' best knowledge, original.
	
	PACS numbers: 44.20.Jj, 03.65.Sq, 03.65.Ge
	
	Keywords: Bohr-Sommerfeld quantisation, Hamilton-Jacobi theory
	
	Submitted to: J. Phys. A: Math. Theor.
	
	

	

§ INTRODUCTION


        The Kepler-Coulomb system is of fundamental importance in physics and chemistry, and is exactly soluble both classically and quantum 
        mechanically <cit.>. Its classical and quantum solutions are related in the sense that applying Bohr-Sommerfeld quantisation to classical 
        bound orbits reproduces exactly the quantum mechanical bound state spectrum derived from the Schrödinger equation. This naturally leads us to ask 
        which other systems possess the same properties, namely classical and quantum solubility with exactness of Bohr-Sommerfeld quantisation. An 
        obvious approach is to consider three-dimensional potentials which generalise the Kepler-Coulomb system.

        In 1997 Dutt et. al. <cit.> showed that the two non-central potentials,
        
    V_A(r)=-κ r-ρcosθ r^2sinθ+γ r^2sin^2θ,    
                    V_B(r)=-κ r-ρcosθ r^2sin^2θ+γ r^2sin^2θ,

        can be solved quantum mechanically using methods from supersymmetric quantum mechanics (SUSYQM) <cit.>. After separation in 
        spherical polar coordinates, only the polar differential equation is changed from the Kepler-Coulomb form, with the substitution z=lntan(θ/2)
        yielding the exactly soluble hyperbolic Scarf or hyperbolic Rosen-Morse potentials, respectively. The exact quantum solubility of the potentials 
        V_A(r) and V_B(r) led us to consider their classical and semiclassical solubility.

        It turns out that V_B(r) has been extensively studied in the literature. Makarov et. al. <cit.> showed that V_B(r) can be
        separated classically and quantum mechanically in both spherical and parabolic coordinates. Kibler and Campigotto <cit.> later solved the 
        system both classically and quantum mechanically in parabolic coordinates, and demonstrated the exactness of Bohr-Sommerfeld quantisation. In a
        following paper,  Kibler et. al. <cit.> solved the Schrödinger equation in spherical, parabolic and prolate spheroidal coordinates, and derived
        the coefficients of the interbasis expansions. The special case where ρ=0 is known as the Hartmann potential <cit.>, and was 
        originally introduced to model ring-shaped molecules such as benzene <cit.>.

        In the following we consider potentials of the form
        
    V(r) = V_1(r) + V_2(θ) r^2,

        where we will eventually set V_1(r) equal to the Kepler-Coulomb potential, V_1(r)=-κ/r, and either 
        V_2(θ)=-ρθ+γ cosec^2 θ or 
        V_2(θ)=-ρθ cosec θ+γ cosec^2 θ, to obtain the non-central potentials V_A(r) and 
        V_B(r), respectively. In section 2 we separate the Hamilton-Jacobi equation for the general potential, V(r), and perform the radial 
        and azimuthal integrals. In section 3 we perform the polar integral for the potential V_A(r), and hence construct its classical solution. We 
        then perform Bohr-Sommerfeld quantisation of V_A(r) in section 4, and compare this with the quantum mechanical result in section 5. The
        classical solution, Bohr-Sommerfeld quantisation, and quantum solution of the Makarov-Kibler potential, V_B(r), are then presented in 
        sections 6, 7, and 8, respectively.


	

§ HAMILTON-JACOBI EQUATION FOR NON-CENTRAL POTENTIALS
 


	In spherical polar coordinates, the Hamilton-Jacobi (HJ) equation for potentials of the form given in eq. <ref> can be written as,
	
    ∂ S∂ t + 1 2μ[(∂ S∂ r)^2
                    + 1 r^2(∂ S∂θ)^2 
                    + 1 r^2sin^2θ(∂ S∂ϕ)^2] 
                    + V_1(r) + V_2(θ) r^2 = 0,

	where μ is the particle's mass. Substituting a solution of the form,
	
    S = -ε t + α_ϕϕ + W_1(r) + W_2(θ),

	allows for the separation of the HJ equation into first order non-linear differential equations for W_1(r) and W_2(θ), upon introduction 
        of a separation constant α_θ^2 <cit.>. Rearranging these differential equations yields
	
        
    W_1(r) = ∫ dr √(2μ(ε-V_1(r))-α_θ^2 r^2), 
                    
       W_2(θ) = ∫ dθ √(α_θ^2-α_ϕ^2sin^2θ-2μ V_2(θ)).

        
	
        We may identify ε as the total energy, and α_ϕ as the z-component of the angular momentum. The HJ equations of motion
        (EOMs) are then given by,
	
	
    β_r = ∂ S∂ε = 
                       -t + ∫ dr μ√(2μ(ε-V_1(r))-α_θ^2 r^2),
                    
       β_θ = ∂ S∂α_θ =
                       -∫ dr α_θ r^2√(2μ(ε-V_1(r))-α_θ^2 r^2)
                       + ∫ dθ α_θ√(α_θ^2-α_ϕ^2 cosec^2 θ-2μ V_2(θ)),
                    
       β_ϕ = ∂ S∂α_ϕ = 
                      ϕ - ∫ dθ α_ϕsin^2θ√(α_θ^2-α_ϕ^2 cosec^2 θ-2μ V_2(θ)),

	
	
	where the β_i are constants. We may interpret β_r and β_ϕ as the initial values of time, -t_0, 
        and azimuthal angle, ϕ_0, respectively; we may set them equal to zero without loss of generality. Since we are concerned with bound orbits, 
        we require ε < 0.
	
	The constant α_θ can be related to the system's angular momentum, using α_ϕ=p_ϕ, as
	
    α_θ^2 = p_θ^2 + p_ϕ^2sin^2θ + 2μ V_2(θ).

	This shows that α_θ has a central piece (the first two terms), and a non-central piece (the final term). 
        Clearly α_θ no longer has the physical interpretation as the total angular momentum.
	
	For our systems V_1(r) = -κ/r, leaving the radial integrals unchanged from the Kepler-Coulomb problem. Performing these integrals, 
        we find a parametric relation between r and t in terms of an intermediary variable w <cit.>. In particular we obtain
	
	
    r =1 2(r_1+r_2)-1 2(r_2-r_1)cos w, 
                    
       t = √(μ(r_1+r_2)^3 8κ)(w-r_2-r_1 r_2+r_1sin w),

	        
	
	where r_1,2 are the minimum and maximum values of the orbital radius, respectively, given by
	
    r_1,2 = κ 2|ε|∓√(κ^2 4|ε|^2 - α_θ^2 2μ|ε|).

	Evaluating the radial integral in eq. <ref> then yields
	
    r = 2r_1r_2[r_1+r_2 + (r_2-r_1)cos(ψ-β_θ)]^-1,

	with ψ simply being the remaining polar integral in eq. <ref>, 
        
    ψ=∫α_θ dθ√(α_θ^2-α_ϕ^2 cosec^2 θ-2μ V_2(θ)).

	Finally we can set β_θ = 0 without consequence, as this simply dictates the initial radial position.


	

§ CLASSICAL MOTION IN THE COTANGENT POTENTIAL V_A(R,Θ)



        We first consider V_A(r,θ), for which V_2(θ)=-ρθ+γ cosec^2θ, and set γ=0; we may recover
        the results for γ≠ 0 by replacing α_ϕ^2→α_ϕ^2=α_ϕ^2+2μγ inside the square 
        roots of eq. <ref> and eq. <ref>. To complete the classical solution, we need to evaluate the integrals
	
	
    ψ = ∫α_θ dθ√(α_θ^2-α_ϕ^2 cosec^2 θ-2μρθ), 
                    
       ϕ = ∫α_ϕ dθsin^2θ√(α_θ^2-α_ϕ^2 cosec^2 θ-2μρθ).

	        
	
	The latter integral is performed using the substitution u = θ, which leads to the EOM between 
        θ and ϕ,
	
    [α_θ^2α_ϕ^2+μ^2ρ^2α_ϕ^4 - 1]^1 2cosϕ 
                     +μρα_ϕ^2 = θ.

	This shows that the periods of motion in θ and ϕ are identical for the case γ = 0. From eq. <ref> it 
        follows that the minimum and maximum values of θ in the motion are,
	
    θ_1,2 = ^-1(μρα_ϕ^2±[α_θ^2α_ϕ^2+μ^2ρ^2α_ϕ^4-1]^1 2).

	
	To calculate the integral for ψ, we first change variable from θ to ϕ using eq. <ref> to obtain
	
    ψ = A∫dϕ 1+(B+Ccosϕ)^2,

	with the constants A, B, and C given by
	
    A = α_θα_ϕ,      
                    B = μρα_ϕ^2,      
                    C = √(A^2+B^2-1).

	We next employ the tangent half-angle substitution s=tan(1 2ϕ) to yield
	
    ψ = 2A (B-C)^2+1∫(s^2+1) ds s^4+2Ds^2+E,

	where D and E are given by
	
    D = 2-A^2 (B-C)^2+1,      
                    E = (B+C)^2+1 (B-C)^2+1.

	Further progress is made by factorising the quartic in eq. <ref> as
	
    s^4+2Ds^2+E = (s^2+2Fs+G)(s^2-2Fs+G),

	where F and G are given by,
	
    G = √(E),      
                    F=√(1 2(G-D)).

	The integrand of eq. <ref> may then be split up into partial fractions, yielding the final result
	
    ψ(ϕ) = A (B-C)^2+1{(G-1) 4FGln(s^2-2Fs+G s^2+2Fs+G) 
        + (G+1) 2GH[tan^-1(s+F H) + tan^-1(s-F H)]},

        where H is given by
        
    H=√(G-F^2)=√(1 2(G+D)).
 
	Unfortunately, this appears to be the most concise manner in which to present this solution. We have therefore obtained a complete classical
        solution allowing us to use ϕ as the driving variable when plotting the orbits. Starting from ϕ, we can determine θ(ϕ) from
        eq. <ref>, and ψ(ϕ) from eq. <ref>; we can then determine r(ψ) and t(ψ) from
        eq. <ref> and eq. <ref>, respectively, giving a complete description of the motion.
	
	In order to ensure bound orbits, we have to impose the following restrictions on the separation constants,
	
    0 ≤α_θ^2≤κ^2μ 2|ε|,        α_ϕ^2(α_θ^2-α_ϕ^2)+μ^2ρ^2≥ 0,         α_ϕ^2≥ 0.

	The first of these inequalities is found by requiring that r_1,2 be real and positive, whilst the second is found by requiring that θ_1,2
        be real. The third inequality is self-evident as α_ϕ is still the z angular momentum of the system. The second inequality shows that 
        α_θ^2<α_ϕ^2 is now a possibility, emphasizing the fact that α_θ is no longer the total angular momentum.
        The conditions on ε, α_θ^2 and α_ϕ^2 may be rewritten as
        
    ε<0,        
                    0 ≤α_θ^2≤κ^2μ 2|ε|,        
                    0≤α_ϕ^2≤1 2[√(α_θ^4+4μ^2ρ^2)+α_θ^2].
 

        In fig. <ref> we show two representative examples of orbits in the potential V_A( r) for γ=0. These orbits are generally
        not closed since the periods of the motion in r and θ,ϕ are not commensurate except for special parameter values.
	
	

        The motion takes place on the surface defined by eq. <ref>, which upon multiplication by rsinθ and rewriting in 
        cartesian coordinates becomes
        
    Cx+B√(x^2+y^2)=z  →  
                    (C^2-B^2)x^2-B^2y^2+z^2-2Czx=0.

        This is clearly a quadric surface, and upon diagonalisation of the corresponding symmetric matrix we find that this is the elliptic cone given by
        
    [√(α_θ^4+4μ^2ρ^2)-α_θ^2α_ϕ^2]x'^2
                    +2μ^2ρ^2α_ϕ^4y'^2
                    -[α_θ^2+√(α_θ^4+4μ^2ρ^2)α_ϕ^2]z'^2=0,

	where
	
    [ x'; y'; z' ] = 
                    [  cosθ_c       0 -sinθ_c;       0       1       0;  sinθ_c       0  cosθ_c ][ x; y; z ],

	corresponding to an anticlockwise rotation about the y axis by an angle θ_c defined by
	
    tanθ_c=√(α_θ^4+4μ^2ρ^2)+α_θ^2-2α_ϕ^2 
                                            2√(α_ϕ^2(α_ϕ^2-α_θ^2)+μ^2ρ^2).

        The axis of the elliptic cone therefore has polar angles (θ_c,ϕ_c) where ϕ_c=ϕ_0 (note that we have previously
        set ϕ_0=0, but setting ϕ_0 0 would just cause rotation of the xz-plane by ϕ_0). The half-angles of the elliptic
        cone in the x'z'- and y'z'-planes, θ_x'z' and θ_y'z', are then
        
	
    tanθ_x'z' = √(√(α_θ^4+4μ^2ρ^2)+α_θ^2√(α_θ^4+4μ^2ρ^2)-α_θ^2)
                                               = √(α_θ^4+4μ^2ρ^2)+α_θ^2 2μρ,
       tanθ_y'z' = α_ϕμρ√(1 2[√(α_θ^4+4μ^2ρ^2)+α_θ^2]).

                
	
        We note that straightforward geometry implies that θ_c=1 2(θ_2-θ_1) and θ_x'z'=1 2(θ_2+θ_1);
        these can be shown to be equivalent to the previous results using trigonometrical identities.
        
        From the above results, we can see how the shape of the elliptic cone varies as a function of the parameters of the problem.
        If we set ρ=0, we find that θ_x'z'=θ_y'z'=π/2 and cosθ_c=α_ϕ/α_θ, corresponding to motion in
        the plane perpendicular to the conserved angular momentum, as expected for the Kepler-Coulomb problem. For ρ 0, θ_x'z'
        increases from π/4 to π/2 as α_θ increases from 0 to ∞; θ_x'z' does not depend upon α_ϕ.
        Alternatively, if we fix α_θ and increase ρ, the cone folds since θ_x'z' decreases from π/2 to π/4 as ρ increases
        from 0 to ∞. For fixed values of ρ and α_θ, θ_y'z' increases from 0 to θ_x'z' as α_ϕ
        increases from 0 to its maximum value α_ϕ^2=1 2[√(α_θ^4+4μ^2ρ^2)+α_θ^2].
	
        Turning our attention now to the γ≠ 0 case, we replace α_ϕ by α_ϕ=√(α_ϕ^2+2μγ) 
        as appropriate in our previous calculations. The equation for θ(ϕ) becomes
	
    [α_θ^2α_ϕ^2 + μ^2ρ^2α_ϕ^4 - 1]^1 2cos(α_ϕα_ϕϕ) + μρα_ϕ^2 = θ,

	whilst the form of ψ(ϕ) given in eq. <ref> remains the same with the changes of parameter
	
    s = tan(α_ϕα_ϕϕ 2),          
                    A = α_θα_ϕ,          
                    B = μρα_ϕ^2.

	The final change is in the second inequality of eq. <ref>, which becomes 
        α_ϕ^2(α_θ^2-α_ϕ^2)+μ^2ρ^2≥ 0.
        The conditions on ε, α_θ and α_ϕ are now
        
    -γκ^2 (4γ^2-ρ^2)<ε<0,
       max[0,μ(4γ^2-ρ^2) 2γ]≤α_θ^2≤κ^2μ/2|ε|,
       0≤α_ϕ^2≤1 2[√(α_θ^4+4μ^2ρ^2)+α_θ^2]-2μγ.

        The additional restrictions on ε and α_θ^2 only occur when γ>1 2ρ, and arise from the fact
        that α_ϕ^2 must be non-negative.
	
        In fig. <ref> we show two representative examples of orbits in the potential V_A( r) for γ≠ 0.
        These orbits are generally not closed since the periods of the motion in r, θ and ϕ are not commensurate except for special  
        parameter values. Since the motion in θ and ϕ is generally incommensurate, the orbit is not confined to a fixed surface unless
        α_ϕ/α_ϕ is rational. The typical motion has the periods of the r, θ and ϕ motions all irrationally
        related, leading to the type of orbit seen in fig. <ref>a, whilst an orbit with the θ and ϕ motions
        rationally related is shown in fig. <ref>b.
	
	
	

	

§ BOHR-SOMMERFELD QUANTISATION OF THE COTANGENT 
        POTENTIAL V_A(R,Θ)

        

        Bohr-Sommerfeld quantisation (BSQ) is an extension of Bohr's 1913 quantum theory of the hydrogen atom. It is part of the “old quantum theory”
        in which quantum conditions are imposed on the classical solution of a problem. This was superseded after 1925 by the “new quantum theory”
        of Born, Heisenberg and Schrödinger, which is the physically correct theory. In general BSQ gives an incorrect result, but there are special systems 
        such as the harmonic oscillator and hydrogen atom, for which the correct quantum mechanical result is obtained. We will demonstrate that
        V_A(r,θ) is one such special system for which BSQ is exact.
         
        To apply Bohr-Sommerfeld quantisation, we must first rewrite our classical equations in terms of action-angle variables <cit.>. 
        Starting from a Hamiltonian description of our system, with the coordinates, q_i, and momenta, p_i, each showing periodic motion, the action
        variables, J_i, are defined by
	
    J_i = 1 2π∮ p_i dq_i
                            = 1 2π∮∂ W_i∂ q_i dq_i

	where we assume the Hamilton-Jacobi equation has a separable solution
        
    S=-ε t+∑_i W_i(q_i;{α_j}).

        The J_i({α_j}) form a new set of constant momenta, which are functions of the HJ separation constants, α_i. Since the energy
        is a separation constant, we can write the Hamiltonian as a function of the J_i. Their conjugate coordinates are the angle variables, ξ_i, 
        defined by
        
    ξ_i=∂ W∂ J_i=∑_j ∂ W_j(q_j;{α_k})∂ J_i.

        The time evolution of the action variables is given by
	
    ξ_i(t) = ξ_i(0)+ω_i t,
                      where  ω_i({J_k}) = ∂ H({J_k})∂ J_i;

	ω_i is the constant frequency associated with ξ_i.

        The final step in Bohr-Sommerfeld quantisation is to set J_i=(n_i+ν_i)ħ, where n_i is a non-negative integer, and the Maslov
        index ν_i equals 0 if q_i has no turning points, and 1 2 if q_i oscillates between two turning points. <cit.>
        For our system this means that
        
    J_r = (n_r+1/2)ħ,      
                    J_θ = (n_θ+1/2)ħ,      
                    J_ϕ = n_ϕħ.


        For the general separable system, the action variables are given by
        
        
    J_r=1π∫_r_1^r_2 dr √(2μ(ε-V_1(r))-α_θ^2 r^2),
                    
       J_θ= 1π∫_θ_1^θ_2 dθ √(α_θ^2-α_ϕ^2sin^2θ-2μ V_2(θ)),
                    
       J_ϕ=1 2π∫_0^2π p_ϕ dϕ=α_ϕ.
 
        
        The integral for J_r is unchanged from the Kepler-Coulomb problem and may be written as
        
    J_r=√(2μ|ε|)π∫_r_1^r_2√((r_2-r)(r-r_1)) r dr.

        If we rewrite this as a contour integral around the branch cut between r_1 and r_2, it may be evaluated by deforming the contour and
        considering the residues of the poles at r=0 and r=∞ to obtain
        
    J_r=√(2μ|ε|)[1 2(r_1+r_2)-√(r_1r_2)]
                        =κ√(μ 2|ε|)-α_θ.

        Setting V_2(θ)=-ρθ, and making the substitution u=θ, the integral for J_θ becomes
        
    J_θ=α_ϕπ∫_u_2^u_1√((u_2-u)(u-u_1)) u^2+1 du.

        Once again we rewrite this as a contour integral around the branch cut between u_2 and u_1, and evaluate it by deforming the contour and
        considering the residues of the poles at u=± i and u=∞ to obtain
        
    J_θ=α_ϕ[Re√((u_1-i)(i-u_2))-1]
                                =√(1 2(√(α_θ^4+4μ^2ρ^2)+α_θ^2))-α_ϕ.


        We now rearrange the equations for J_r, J_θ and J_ϕ to write H≡ε as
        
    H=-|ε|=-μκ^2 2{J_r+(J_θ+J_ϕ)
                         [ 1-μ^2ρ^2(J_θ+J_ϕ)^4]^1/2}^2,
	
        and Bohr-Sommerfeld quantisation, using eq. <ref>, finally gives
        
    E(n_r, n_θ, n_ϕ)=-μκ^2 2ħ^2{n_r+1 2+(n_θ+n_ϕ+1 2)
                                                          [ 1-μ^2ρ^2ħ^4(n_θ+n_ϕ+1 2)^4]^1/2}^2.

        To obtain the results for γ 0, we replace J_ϕ≡α_ϕ by 
        J_ϕ≡α_ϕ=√(J_ϕ^2+2μγ), giving
         
    H=-|ε|=-μκ^2 2{J_r+(J_θ+√(J_ϕ^2+2μγ))
                                                [ 1-μ^2ρ^2(J_θ+√(J_ϕ^2+2μγ))^4]^1/2}^2,
	
        and BSQ proceeds as before using eq. <ref>.

        From eqs. <ref> and <ref>, we see that in the γ=0 case, the frequencies associated with 
        θ and ϕ are identical, since J_θ and J_ϕ only occur in the combination J_θ+J_ϕ. In other words,
        ω_r≠ω_θ=ω_ϕ. When γ≠ 0, J_θ+J_ϕ is replaced by J_θ+√(J_ϕ^2+2μγ),
        and it follows that ω_r≠ω_θ≠ω_ϕ. The system then has three independent frequencies, and the orbits for
        γ≠ 0 are very different from those for γ=0. If we start from the Kepler-Coulomb potential, all three frequencies are the same;
        when the θ term is then added, ω_r becomes different to ω_θ=ω_ϕ; when the cosec^2 θ 
        term is finally added, all three frequencies are different.
	

	

§ QUANTUM SOLUTION OF THE COTANGENT POTENTIAL V_A(R,Θ)
 
        

	
        The Schrödinger equation for the general potential given in eq. <ref> is
	
    -ħ^2 2μ∇^2Ψ+(V_1(r)+V_2(θ) r^2)Ψ=EΨ(𝐫).

	Separating variables in the standard manner using Ψ(r,θ,ϕ)=R(r)Θ(θ)Φ(ϕ), we obtain
	
	
    Φ(ϕ) = e^in_ϕϕ,
                    
       d^2R dr^2+2 rdR dr + (2μħ^2(E-V_1(r))-l(l+1) r^2)R = 0, 
                    
       d^2Θ dθ^2+θdΘ dθ+(l(l+1)-2μ V_2(θ)ħ^2
                        -n_ϕ^2sin^2θ)Θ = 0,

	        
	
	where n_ϕ is an integer, and l(l+1) is the common separation constant associated with R(r) and Θ(θ); l will no longer be
        a non-negative integer when V_2(θ)≠ 0.
	
	As in the classical case, the radial equation is unaffected by the non-central potential, and so has the standard hydrogen
        atom radial wavefunction,
	
    R_n_r l(r) = (Qr)^l e^-Qr L_n_r^2l+1(2Qr),
                        where    
                    Q=√(2μ|E|ħ^2)

        although l is no longer a non-negative integer. The L_n_r^2l+1(w) are associated Laguerre polynomials,  and the system has energy
	
    E=-μκ^2(n_r+l+1)^2


        To solve the polar equation for the cotangent potential, we substitute V_2(θ)=-ρθ in eq. <ref>, and change
        variable to u=θ, which yields
        
    d^2Θ du^2+u (1+u^2)dΘ du+[l(l+1) (1+u^2)^2+2μρħ^2u (1+u^2)^2
                    -n_ϕ^2 (1+u^2)]Θ=0.
	
        We next remove the double poles at u=± i by setting 
        
    Θ(u)=exp[-α 2^-1u](1+u^2)^2β-1 4χ(u),

        so that χ(u) satisfies the Romanovski equation
        
    (1+u^2)d^2χ du^2+(2β u+α)dχ du-n_θ(n_θ+2β-1)χ=0,

        where α, β and l obey the conditions
        
        
    (β-1)^2-1 4α^2=(l+1 2)^2
       α(β-1)=-2μρħ^2
       n_ϕ^2-(β-1 2)^2=n_θ(n_θ+2β-1).

        
        These can be solved to give the results for α, β and l,
        
        
    β=1 2-n_θ-n_ϕ
       α=2μρħ^2(n_θ+n_ϕ+1 2)
       l+1 2=(n_θ+n_ϕ+1 2)
                                                        [1-μ^2ρ^2ħ^4(n_θ+n_ϕ+1 2)^4]^1 2.

                
        
        The normalisable solutions of the Romanovski equation are the Romanovski polynomials, which have weight function,
        W^(α,β)(u), and corresponding Rodrigues formula,
        
        
    W^(α,β)(u)=(1+u^2)^β-1e^-α^-1u,
        R^(α,β)_n(u)=1 2^n n!1 W^(α,β)(u)d^n du^n[(1+u^2)^n W^(α,β)(u)].

        
        They are related to Jacobi polynomials of complex parameters and imaginary argument by
        
    R^(α,β)_n(u)=(-i)^n P^(β-1+iα 2,β-1-iα 2)_n(iu),

        but it is more useful to treat them as real polynomials. They were first discovered by Routh in 1884 <cit.>, and the later rediscovered by
        Romanovski in 1929 <cit.>. Their applications in physics have recently been discussed by Raposo et al <cit.> and 
        Alvarez-Castillo <cit.> and we are following their definitions. We note that the orthogonality of the polar wavefunctions for the cotangent
        potential is not the standard orthogonality with respect to the weight function occuring in the Rodrigues formula. The wavefunctions for different
        n_θ have different values of the parameters α and β. The fact that these wavefunctions are orthogonal is, however, guaranteed 
        by the Sturm-Liouville nature of the original problem.

        The unnormalised polar wavefunctions for the cotangent potential are therefore
        
    Θ(θ)=exp[-μρħ^2(n_θ+n_ϕ+1 2) θ]
                                   (sinθ)^n_θ+n_ϕ  R_n_θ^(α,1 2-n_θ-n_ϕ)(θ),

        and the energy for the complete wavefunction labelled by quantum numbers (n_r,n_θ,n_ϕ) is
        
        
    E(n_r,n_θ,n_ϕ)   =-μκ^2 2ħ^2(n_r+l+1)^2
       =-μκ^2 2ħ^2{n_r+1 2+(n_θ+n_ϕ+1 2)
                                                          [ 1-μ^2ρ^2ħ^4(n_θ+n_ϕ+1 2)^4]^1/2}^2,

        
        which agrees with the Bohr-Sommerfeld result of eq. <ref>. The case where γ≠ 0 is then obtained by replacing
        n_ϕ^2 by n_ϕ^2+2μγħ^2. It follows that Bohr-Sommerfeld quantisation exactly reproduces the quantum
        mechanical spectrum for the cotangent potential.
       


	

§ CLASSICAL MOTION IN THE MAKAROV-KIBLER POTENTIAL V_B(R,Θ)



	We now consider the Makarov-Kibler potential, where V_2(θ)=-ρ cosec θθ+γ cosec^2 θ, 
        and first set γ=0. Following the analysis of section <ref>, we need to evaluate the integrals
	
	
    ψ = ∫α_θ dθ√(α_θ^2-α_ϕ^2 cosec^2 θ
                               +2μρ cosec θθ),
                    
       ϕ = ∫α_ϕ dθsin^2θ√(α_θ^2-α_ϕ^2 cosec^2 θ 
                               +2μρ cosec θθ).

	
        The solution of the first equation is found by changing variable to v=cosθ, yielding
	
    √(1-α_ϕ^2α_θ^2+μ^2ρ^2α_θ^4)cosψ 
                    +μρα_θ^2=cosθ.

        This shows that the periods of motion in ψ and θ are the same when γ=0. From eq. <ref> this means
        that the periods of motion in r and θ are the same when γ=0. The minimum and maximum values of θ in the 
        motion are
	
    θ_1,2=cos^-1[μρα_θ^2∓√(1-α_ϕ^2α_θ^2+μ^2ρ^2α_θ^4)].


        To evaluate the integral for ϕ, we change variable from θ to ψ using eq. <ref> to obtain	
        
    ϕ=∫dψ 1-(M+Ncosψ)^2,

        with the constants M and N given by
        
    M=μρα_θ^2,        
                   N=√(1-α_ϕ^2α_θ^2+μ^2ρ^2α_θ^4).

        The tangent half-angle substitution finally gives
        
    ϕ(ψ)=1√((1+M)^2-N^2)tan^-1[√(1+M-N 1+M+N)tan(1 2ψ)]
             +1√((1-M)^2-N^2)tan^-1[√(1-M+N 1-M-N)tan(1 2ψ)].

        We have therefore obtained a complete classical solution allowing us to use ψ as the driving variable when plotting orbits. Starting from
        ψ, we can determine θ(ψ) from eq. <ref>, ϕ(ψ) from eq. <ref>, r(ψ) from eq. <ref>
        and t(ψ) from eq. <ref>, giving a complete description of the motion.
        The period of the motion in ϕ is clearly different from that in ψ, and hence different from that in r and θ.

	In order to ensure bound orbits, we have to impose the following restrictions on the separation constants,
        
    0≤α_θ^2≤κ^2μ 2|ε|,        α_θ^2(α_θ^2-α_ϕ^2)+μ^2ρ^2≥ 0,        α_ϕ^2≥ 2μρ.

        The first of these inequalities is found by requiring that r_1,2 be real and positive, whilst the second is found by requiring that θ_1,2
        be real. The third inequality is found by requiring that 1-M≥ N, which is necessary if eq. <ref> describes a periodic solution.
        The conditions on ε, α_θ^2 and α_ϕ^2 may be rewritten as
        
    ε<0,        
                    0≤α_θ^2≤κ^2μ 2|ε|,        
                    2μρ≤α_ϕ^2≤α_θ^2+μ^2ρ^2α_θ^2.

	
        The orbits in the Makarov-Kibler potential are similar to orbits in the cotangent potential, in that they lie on a quadric surface. To show that this
        is the case, we use eq. <ref> and eq. <ref> to obtain
	
    pr+qrcosθ=2r_1r_2,

        where the constants p and q are given by
        
    p=r_1+r_2-(r_2-r_1)μρ√(α_θ^4-α_θ^2α_ϕ^2+μ^2ρ^2),        
                    q=(r_2-r_1)α_θ^2√(α_θ^4-α_θ^2α_ϕ^2+μ^2ρ^2).

	Upon rewriting in cartesian coordinates and rearranging, this gives
	
    p^2x^2+p^2y^2+(p^2-q^2)z^2+4qr_1r_2z=4r_1^2r_2^2.

	Shifting the z-axis to eliminate the linear term in eq. <ref> using z' = z + 2r_1r_2q/(p^2-q^2), the equation
        for the quadric surface becomes
	
    (p^2-q^2) 4r_1^2r_2^2 (x^2+y^2)+(p^2-q^2)^2 4p^2r_1^2r_2^2 z'^2=1.

        The nature of the surface depends upon whether p^2-q^2 is positive, negative or zero, in which cases it is an ellipsoid, hyperboloid of
        two sheets or paraboloid, respectively. All three situations are possible depending upon the parameter values. 
        In fig. <ref> we show orbits in the Makarov-Kibler potential for γ=0 which lie on an ellipsoidal surface, whilst in 
        fig. <ref> we show orbits which lie on a hyperboloidal surface.
	
	
	
	
	

        Consider now the γ≠ 0 case, which we again treat by replacing α_ϕ by 
        α_ϕ=√(α_ϕ^2+2μγ) as appropriate in our previous calculations. The equation for θ(ψ)
        then becomes
        
    √(1-α_ϕ^2α_θ^2+μ^2ρ^2α_θ^4)cosψ 
                    +μρα_θ^2=cosθ,

        whilst that for ϕ(ψ) becomes
        
    ϕ(ψ)=α_ϕα_ϕ{1√((1+M)^2-N^2)tan^-1[√(1+M-N 1+M+N)tan(1 2ψ)].
             +.1√((1-M)^2-N^2)tan^-1[√(1-M+N 1-M-N)tan(1 2ψ)]},

        where M and N are now defined by
         
    M=μρα_θ^2,        
                   N=√(1-α_ϕ^2α_θ^2+μ^2ρ^2α_θ^4).

        We see that the motion in ψ and θ, and hence in r and θ, maintains the same period. The period of the motion in
        ϕ is changed when γ≠ 0, but this motion already generally has a different period from the r and θ motion.
        It follows that setting γ≠ 0 has no qualitative effect on the orbits in the Makarov-Kibler potential, with the orbits remaining
        confined to the same types of quadric surfaces.


	

§ BOHR-SOMMERFELD QUANTISATION OF THE MAKAROV-KIBLER POTENTIAL V_B(R,Θ)
 
        

        We now perform Bohr-Sommerfeld quantisation for the Makarov-Kibler potential. The results for J_r and J_ϕ are the same
        as before, and are given by eq. <ref> and eq. <ref>, respectively. The integral for J_θ in the case
        where γ=0 is
        
        
    J_θ=1π∫_θ_1^θ_2 dθ √(α_θ^2-α_ϕ^2 cosec^2 θ
                          +2μρ cosec θθ)
       =α_θπ∫_v_2^v_1√((v_1-v)(v-v_2)) dv 1-v^2,

        
        where we have made the substitution v=cosθ. This may be rewritten as a contour integral around the branch cut between v_2
        and v_1, and evaluated by considering the residues of the poles at v=± 1 and v=∞ to obtain
        
    J_θ=α_θ-1 2(√(α_ϕ^2+2μρ)+√(α_ϕ^2-2μρ))

        We now rearrange the equations for J_r, J_θ and J_ϕ to write H≡ε as
        
    H=-μκ^2 2[J_r+J_θ+1 2(√(J_ϕ^2+2μρ)+√(J_ϕ^2-2μρ) )]^2,
 
        and, from eq. <ref>, Bohr-Sommerfeld quantisation gives
        
    E(n_r, n_θ, n_ϕ)=-μκ^2 2ħ^2[ n_r+n_θ+1+1 2(√(n_ϕ^2+2μρħ^2)+
                                                         √(n_ϕ^2-2μρħ^2) )]^2.
 
        The result for γ≠ 0 is then simply found by replacing J_ϕ^2 by J_ϕ^2+2μγ to give
        
    H=-μκ^2 2[J_r+J_θ+1 2(√(J_ϕ^2+2μ(γ+ρ))
                           +√(J_ϕ^2+2μ(γ-ρ)) )]^2,

        with Bohr-Sommerfeld quantisation giving
         
    E(n_r, n_θ, n_ϕ)=-μκ^2 2ħ^2[ n_r+n_θ+1+1 2(
                                                          √(n_ϕ^2+2μ(γ+ρ)ħ^2)
                                                          +√(n_ϕ^2+2μ(γ-ρ)ħ^2) )]^2.

        These results agree with those of Kibler and Campigotto <cit.>, which were obtained by separating the classical motion in
        parabolic polar coordinates.
        
	From eq. <ref>, we see that ω_r=ω_θ≠ω_ϕ even when γ≠ 0, as already
        seen from the classical solution.
	

	

§ QUANTUM MECHANICS OF THE MAKAROV-KIBLER POTENTIAL


        To solve the polar Schrödinger equation for the Makarov-Kibler potential, we substitute V_2(θ)=-ρ cosec θθ in
        eq. <ref>, and change variable to v=cosθ, which yields
        
    d^2Θ dv^2-2v (1-v^2)dΘ dv+[l(l+1) (1-v^2)+2μρħ^2v(1-v^2)^2
                    -n_ϕ^2 (1-v^2)^2]Θ=0.

        We next remove the double poles at v=± 1 by setting Θ(v)=(1-v)^α/2(1+v)^β/2χ(v), where
        α=√(n_ϕ^2-2μρħ^2) and β=√(n_ϕ^2+2μρħ^2), and χ(v) satisfies 
        
    (1-v^2)d^2χ dv^2+[(β-α)-(α+β+2)v]dχ dv+n_θ(n_θ+α+β+1)χ=0

        with n_θ=l-1 2(α+β). The normalisable solutions of this equation are the Jacobi polynomials, 
        P_n^(α,β)(v), defined by the Rodrigues formula
        
    P_n^(α,β)(v)=(-1)^n 2^n n!(1-v)^-α(1+v)^-βd^n dv^n[(1-v)^n+α(1+v)^n+β].

        The unnormalised polar wavefunctions for the Makarov-Kibler potential are therefore
        
    Θ(θ)=(1-cosθ)^α/2(1+cosθ)^β/2P_n_θ^(α,β)(cosθ),

        and the energy for the complete wavefunction labelled by quantum numbers (n_r,n_θ,n_ϕ) is
        
    E(n_r,n_θ,n_ϕ)   =-μκ^2 2ħ^2(n_r+l+1)^2
       =-μκ^2 2ħ^2[n_r+n_θ+1+1 2(√(n_ϕ^2+2μρħ^2)
                                                        +√(n_ϕ^2-2μρħ^2))]^2,

        which agrees with the Bohr-Sommerfeld result of eq. <ref>. The case where γ≠ 0 is then
        obtained by replacing n_ϕ^2 by n_ϕ^2+2μγħ^2, which would yield eq. <ref>.
        It follows that Bohr-Sommerfeld quantization exactly reproduces the correct quantum mechanical spectrum for the Makarov-Kibler potential.

		
	

§ CONCLUSIONS

        We have shown that the cotangent and Makarov-Kibler potentials, V_A( r) and V_B( r), defined in eq. <ref>, are
        classically and quantum mechanically exactly soluble in spherical polar coordinates. Moreover, the quantum mechanical spectrum can be
        obtained from the classical solution in both cases via Bohr-Sommerfeld quantisation. However, the lifting of degeneracies of the frequencies of 
        the angle variables in the
        classical solution differs for the two potentials. Starting with the Kepler-Coulomb potential, all three frequencies are identical,
        ω_r=ω_θ=ω_ϕ. In the Makarov-Kibler potential, adding the -ρ cosec θθ/r^2 term then
        lifts the degeneracy of the ϕ-motion, so that ω_r=ω_θ≠ω_ϕ. Adding the  γ cosec^2θ/r^2
        term does not further lift the degeneracy. By contrast, for the cotangent potential, adding the -ρθ/r^2 term lifts the degeneracy of the 
        r-motion,
        so that ω_r≠ω_θ=ω_ϕ. Adding the γ cosec^2θ/r^2 term then lifts the degeneracy of the
        ϕ-motion, so that ω_r≠ω_θ≠ω_ϕ for the general motion. Another difference between the two potentials
        is that the Makarov-Kibler potential is superintegrable, being soluble in spherical polar, parabolic, and prolate spheroidal coordinates, whilst the
        cotangent potential is only soluble in spherical polar coordinates.

        A further interesting feature is that the classical orbits in the Makarov-Kibler and cotangent potentials both lie on quadric surfaces, the latter only in 
        the case γ=0. In the Makarov-Kibler potential, these surfaces are ellipsoids, parabaloids, or hyperboloids of two sheets, according to the
        values of the constants of motion. In the cotangent potential, these surfaces are elliptic cones.

        Finally the identification of a new system (the cotangent potential) for which Bohr-Sommerfeld quantisation is exact, poses the general question
        of why certain special systems have this property and others do not. The fact that the cotangent potential is not superintegrable indicates that
        this is not generally a requirement.

	
	
        Cordani
        Cordani B 2003 The Kepler Problem (Basel: Springer)

        Dutt
	Dutt R, Gangopadhyaya A and Sukhatme U P 1997 Am. J. Phys. 65 400

        SUSYQM
	Cooper F, Khare A and Sukhatme U P 2001 Supersymmetry and Quantum Mechanics  (Singapore: World Scientific)

        Makarov
        Makarov A A, Smorodinsky J A, Valiev K and Winternitz P 1967 Nuovo Cimento A 52, 1061

        Kibler1
	Kibler M and Campigotto C 1993 Int. J. Quantum Chem. 45 209

        Kibler2
        Kibler M, Mardoyan L G and Pogosyan G S 1994 Int. J. Quantum Chem. 52 1301

        Hartmann
        Hartmann H 1976 Theor. Chim. Acta 24 201
	
	Hartmann_benzene
        Hartmann H and Schuch D 1980 Int. J. Quantum Chem. 18 125
		
        Goldstein
	Goldstein H 1980 Classical Mechanics (London: Addison-Wesley)
		
	beige_book
	Pars L A 1965 A Treatise on Analytical Dynamics  (London: Heinemann)

	Routh
	Routh E J 1884 Proc. London Math. Soc. 16 245
	
	Romanovski
        Romanovski V 1929 C. R. Acad. Sci. (Paris) 188 1023
	
	Raposo
        Raposo A P, Weber H J, Alvarez-Castillo A E and Kirchbach M 2007 Centr. Eur. J. Phys. 5 253
	
	Alvarez
	Alvarez-Castillo D E and Kirchbach M 2007 Rev. Mex. Fis. E 53 (2) 143
	

	
	
	
	
	
