Series Solutions for Atom Guides 

by William M. Golding 

ARL-TR-5335 September 2010 

Approved for public release; distribution unlimited. 
NOTICES 

Disclaimers 

The findings in this report are not to be construed as an official Department of the Army position 

unless so designated by other authorized documents. 

Citation of manufacturer’s or trade names does not constitute an official endorsement or 

approval of the use thereof. 

Destroy this report when it is no longer needed. Do not return it to the originator. 
Army Research Laboratory 

Adelphi, MD 20783-1197 

ARL-TR-5335 September 2010 

Series Solutions for Atom Guides 

William M. Golding 

Sensors and Electron Devices Directorate, ARL 

Approved for public release; distribution unlimited. 
ii

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing the 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 

September 2010 

2. REPORT TYPE 

Final 

3. DATES COVERED (From - To) 

October 2009 to August 2010 

4. TITLE AND SUBTITLE 

Series Solutions for Atom Guides 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

William M. Golding 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

U.S. Army Research Laboratory 

ATTN: RDRL-SEE-O 

2800 Powder Mill Road 

Adelphi MD 20783-1197 

8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

ARL-TR-5335 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution unlimited. 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

In an atomic waveguide based on magnetic interactions, the forces on different spin components have opposite signs. Thus, if 

one component of a spinor is bound near the guiding center, another component of the spinor is repelled from the center. This is 

simply an example of the Stern-Gerlach effect commonly used to separate particles in different spin states. One difference is 

that, in the atomic guide, the spatial degrees of freedom must be quantized as well as the spin degrees of freedom. In order to

accurately treat the complete quantum problem both inside and outside the guiding region, a power series solution valid for the

whole real axis is developed. Thus, in principle, both the bound and unbound components of the full spinor solution can be 

treated within the same framework. The power series converges and can be readily summed over a large range of the radial 

coordinate by making use of arbitrary precision arithmetic. This technique produces high precision eigenvalues as well as 

detailed information about the behavior of the spinor wavefunction at both large and small distances from the origin, making it

possible to perform detailed studies of atom propagation without resorting to the adiabatic approximation. 

15. SUBJECT TERMS

Atomic waveguides, arbitrary precision arithmetic, atomic physics, quantum mechanics, atom interferometry, atom chips 

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION 

 OF 

 ABSTRACT 

UU 

18. NUMBER 

 OF 

 PAGES 

32 

19a. NAME OF RESPONSIBLE PERSON

William M. Golding 

a. REPORT 

Unclassified 

b. ABSTRACT 

Unclassified

c. THIS PAGE 

Unclassified

19b. TELEPHONE NUMBER (Include area code) 

(301) 394-1535 

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18 
iii

Contents 

List of Figures iv

List of Tables v

1. Introduction 1

2. Development of Radial Equations 3

2.1 Hamiltonian for a Magnetic Atom ..................................................................................3

2.2 Eigenstates Common to H and Λz ...................................................................................5

3. Solution by the Frobenius Technique 7

3.1 Uncouple to Produce Fourth Order Differential Equations.............................................7

3.2 Frobenius Series Approach .............................................................................................8

4. Series Evaluations Using Arbitrary Precision Arithmetic 12

4.1 Series Calculations ........................................................................................................12

4.2 Eigenvalues: Shooting or Summing? ............................................................................16

4.3 Eigenfunctions ...............................................................................................................19

5. Conclusions 21

6. References 22

Distribution List 23
iv

List of Figures 

Figure 1. Cross-sectional field plot of the general quadrupole guiding field used in this work. 

Magnetic fields come inwards along x and go outwards along y, resulting in a quadrupole 

null at the center that can be used for atom guiding. .................................................................4

Figure 2. Term plot of the log of the magnitude of the individual terms used in a calculation 

of the first series in equation 27. The target value of radius is 40 and the maximum 

eigenvalue is 30. The largest term in the sequence is roughly 10100 and therefore 120 

digits needs to be used to do an adequate job of summing this series for any value of 

radius up to 40 for eigenvalues less than 30. ...........................................................................13

Figure 3. Calculation by series summation of the lowest two bound components using 15 

digits of precision. This kind of result is representative of hardware floating point. 

Beyond a certain point the calculation breaks down. As the number of digits used in the 

calculation is increased, the target radius can be significantly increased. In this way, 

hardware floating point limits the ability to obtain accurate eigenvalues. ..............................14

Figure 4. The same calculation as in figure 3. The only difference is that 25 digits of 

precision are used. One can see that the calculation can proceed to much larger target 

values as the number of digits is increased. .............................................................................15

Figure 5. The same calculation as displayed in figures 3 and 4, except that 100 digits were 

used in this calculation. The effective target radius can be pushed out to 40 with very 

clean results. This technique lets us calculate significantly more accurate eigenvalues as 

well as perform checks on asymptotic series expansions for future work. ..............................15

Figure 6. This is a blown up version of the tail of the 100-digit calculation displayed in 

figure 5. The plot is from radius of 0 to 40 and the scale is expanded by a factor of 

10,000. Very clean damped oscillations are observed. This is an example of the power of 

the basic arbitrary precision techniques. In the same region where we saw complete 

breakdown of the calculations at 15 and 25 digits in figures 3 and 4, extremely useful 

information about the asymptotic behavior is obtained by simply summing the series 

using 100 digits of precision. This cannot be obtained directly using hardware floating 

point techniques. ......................................................................................................................16

Figure 7. Plot of the first 52 eigenvalues. The red points are mode type 1 eigenvalues and the 

blue points are mode type 2 eigenvalues; an approximate fit to the eigenvalues as a 

function of n is given in the text. .............................................................................................18

Figure 8. The first 10 type 1 mode radial eigenfunctions. The black curve with no zero 

crossings is the ground state wave function. The first excited state is the blue curve. ............20

Figure 9. The first 10 type 2 mode radial eigenfunctions. The black curve is the lowest 

energy type 2 mode but is actually the first excited state of the system. .................................20
v

List of Tables 

Table 1. First 52 energy eigenvalues for our system with alignment 1/2, zero bias field and a 

transverse gradient of b1=1. The sequence of eigenvalues alternates between type 1 and 

type 2 eigenvalues as n is increased. Because of this double counting the index n in the 

table is not the same as the index n used in the plot of figure 7. .............................................19
vi

INTENTIONALLY LEFT BLANK
1

1. Introduction 

A basic type of magnetic guiding structure used on atom chips is formed by placing a long, 

straight current-carrying wire in an externally produced uniform magnetic field that points in a 

direction perpendicular to the track of the wire (Fortágh, 2007). In this configuration, the roughly 

circular magnetic field of the wire is cancelled along a line that parallels the wire, creating a long 

straight line of zero magnetic field. Cold magnetic atoms can be trapped in quantum states in this 

type of magnetic field minimum near to the wire and be guided along parallel to the wire 

direction. This forms a type of magnetic waveguide for atoms that can, in principle, be used for 

constructing atom interferometers and other devices based on atom guiding. 

In an atomic waveguide based on magnetic interactions, the forces on different spin components 

have opposite signs. Thus, if one component of a spinor is bound near the guiding center, another 

component of the spinor is repelled from the center. This is simply an example of the SternGerlach effect commonly used to separate particles in different spin states. One difference 

between the simple Stern-Gerlach experiment, in which only the spin degrees of freedom are 

quantized and the atom waveguide, is that the spatial degrees of freedom must be quantized as 

well as the spin in the atom guide problem. 

In order to accurately treat the complete quantum problem both inside and outside the guiding 

region, a power series solution expanded about the guide center, yet valid for the whole real axis, 

has been developed. Using this series, both the bound and unbound components of the full spinor 

solution can be treated within the same framework. The power series converges quickly at very 

small radii but converges very slowly at relatively small distances from the guide center. 

However, this limitation is seen to be due to the fixed hardware floating point resolution 

available in modern computers. The power series solutions can be readily summed over a wide 

range of the radial coordinate by making use of software floating point techniques implemented 

in arbitrary precision arithmetic libraries. This technique produces high precision eigenvalues as 

well as detailed information about the behavior of the spinor wavefunction at both large and 

small distances from the origin, making it possible to perform detailed studies of atom 

propagation without resorting to the adiabatic approximation (Sukumar, 1997). 

Numerical solutions for the quantum modes and eigenvalues of spin 1/2 and spin 1 particles in a 

quadrupole guide have been published elsewhere (Hinds, 2000; Potvliege, 2001; Lesanovsky, 

2004; Bill, 2006); however, in all of these works, the boundary conditions used at the guide 

center have been chosen incorrectly and do not properly treat the regular singularity that exists 

there. This error in the boundary conditions clearly results from treating the two second-order 

radial equations (our equations 8) as though they uncouple at the origin before trying to perform 

a power series analysis on the full system of radial equations. The power series analysis can be 

performed without first decoupling the system but one must be careful not to assume that 
2

relations like R R        hold as   0. This type of reasoning seems to have caused 

problems with the boundary conditions used in these earlier works. In Blumel (1991), the 

Frobenius method was used on uncoupled equations to properly handle similar types of 

singularities in the problem of a neutron in the magnetic field of a rectilinear current, a different 

but related problem. 

Failure to handle the boundary conditions properly leads to an incorrect description of important 

details of waveguide properties such as the energy spectrum and the ground state of the system. 

In addition, a whole class of modes goes undetected in these calculations simply because the 

boundary conditions are not treated carefully. The quantum properties of an atom guide cannot 

be treated consistently using an inconsistent basis set. Since the quantum system is linear, using 

incorrect boundary conditions in the numerical solution leads to apparent modes that are actually 

superpositions of the true modes of the system. These mixed modes cannot generally be 

stationary states of the system; they are mixtures of nondegenerate states. Various lifetime 

estimates have been performed by several authors (Hinds, 2000; Sukamar, 1997; Potvliege, 

2001; Lesanovsky, 2004; Bill, 2006) using these nonstationary modes. All such estimates are 

questionable, since the mixed modes used for those calculations would be time dependent simply 

due to the normal time evolution of a conservative Hamiltonian system and this point is not made 

clear in those works. Thus, any time-dependent properties calculated using these nonstationary 

modes might actually be due to the natural unitary evolution of the conservative system and not 

representative of an inherent limitation of magnetic guiding. This needs to be clarified before 

attributing a lifetime to the states. 

In a previous work (Golding, 2009), series solutions to the radial differential equations were used 

simply to provide accurate initial conditions appropriate for the regular singularity at the origin. 

However, in that work, the initial conditions were simply used to initialize various hardware 

precision numerical integration routines that are routinely used for solving ordinary differential 

equations (ODE). The shooting technique, along with a root solving approach, was used to solve 

for the eigenvalues of the system as well as plot the radial eigenfunctions. The differential 

equations solved are stiff and it is difficult to judge the accuracy of the eigenvalue results 

obtained using the shooting technique. This is a well known problem with the shooting technique 

that occurs when an exponentially growing solution is present along with an exponentially 

decaying solution. When the value of the decaying solution is too small, round off errors can 

cause the growing exponential solution to become dominant, limiting the resolution of the 

calculated eigenvalues. In second order systems, this problem is often handled by integrating the 

system in both directions and matching the two solutions in a region where both are judged to be 

accurate. However, in the fourth-order system studied here, there are two oscillatory solutions in 

addition to the two exponential solutions that contribute significantly at large distances from the 

center. The oscillatory solutions cause practical difficulties in establishing satisfactory boundary 

conditions for an inward integration of the system since one cannot simply assume that the 

solution simply decays to zero far from the guide center. 
3

In this work, the solutions to the quadrupole guide problem are formulated using a power series 

approach. The behavior at the regular singular point is properly accounted for by using the 

Frobenius series method (Ince, 1956; Bender, 1978) on the fourth-order equations derived from 

the coupled radial spinor equations. This is, in principle, an exact answer but leaves us to solve 

the problem of actually summing the series. Power series generally converge slowly and they are 

often used only to obtained local behavior, as in our previous work. However, the power series 

solutions can easily be summed at large distances from the origin if arbitrary precision arithmetic 

is used throughout the calculation. In the computer algebra programs, Maple and Mathematica 

arbitrary precision arithmetic is implemented using the GNU Multiple Precision (GMP) Library 

(GNU, 2010). The calculations in this work were done using Maple 14. Using arbitrary 

precision, all of the terms in the required sums, both very big and very small, are accurately 

included in the calculations resulting in very accurate solutions for both the guide wave functions 

and eigenvalues that can be confidently used in further detailed investigations of atom guide 

behavior. 

2. Development of Radial Equations 

The equations for the radial wave functions of the quadrupole guide are derived by making use 

of the angular symmetry of the quadrupole magnetic field. This symmetry is expressed as the 

conservation of alignment,  z zz L S (Hinds, 2001; Lesanovsky, 2004; Bergeman, 1989; 

Golding, 2009). Our approach is to find eigenfunctions that are common to both the 

Hamiltonian, H and the alignment z . This approach effectively separates the angular and radial 

dependence of the problem, leaving only a coupled system of radial equations to solve. These 

equations are solved exactly by series techniques and the resulting series are summed directly 

using arbitrary precision arithmetic. 

2.1 Hamiltonian for a Magnetic Atom 

The quantum description of an atom moving in a spatially dependent magnetic field must include 

both internal and external degrees of freedom. This means that the momentum, position, and spin 

degrees of freedom must all be treated as quantum operators. Using m for the mass of the atom 

and M



 for the magnetic moment, the Hamiltonian is written as the sum of the kinetic energy, 

2

2

p

m

 and the interaction energy of the magnetic moment and the field B



. The result is 

 

2

2

p H M Br

m      . (1) 

The magnetic field B r    is independent of the z-direction and is taken to be an ideal quadrupole 

field (equation 2) extending to infinity in the x-y plane and uniform along z. An additional 
4

uniform bias field B0 is added in the z-direction to help control possible spin-dependent losses 

that may occur at the zero field point at the center of the guide. The spatial dependence of the 

ideal quadrupole field is given by 

B r B xx yy B z      1 0  ˆ ˆ ˆ   . (2) 

A cross-sectional plot of this field is shown in figure 1. The guiding center of this field is at the 

origin where the transverse field is zero. The quantity B1 is the magnitude of the transverse field 

gradient. It is taken as greater than or equal to zero in this work, although changing the sign of B1

simply changes the quadrupole field configuration from one with the field point inward along the 

x-axis to one with the field pointing inward along the y-axis. 

Figure 1. Cross-sectional field plot of the general quadrupole guiding field 

used in this work. Magnetic fields come inwards along x and go 

outwards along y, resulting in a quadrupole null at the center that 

can be used for atom guiding. 

The gradient of the transverse field is what provides the trapping force that keeps the atom 

confined in the transverse direction. The transverse field is zero at the center of the guide and the 

only nonzero field component at the very center is the longitudinal bias field B0. The magnitude 

of the transverse field is 2 2 B   Bx y B 1 1 , where  is the radial coordinate in the 

cylindrical system. Even though the field is varying in direction as one goes round the guide, the 

contours of constant field magnitude are simply circles. 
5

The potential is independent of the position along the guide. Therefore, pz is a constant of the 

motion and is ignored here. In future problems, when variations in the potential along the zdirection are considered, the variations in the longitudinal momentum will be included. 

Variations in the transverse potential should show up as variations in the guide propagation 

constant and this will affect the accuracy of sensitive measurements based on interferometry. 

In the model of the atom described by equation 1, the magnetic moment of the atom is just the 

magnetic moment of a single outer electron. The atomic model used here can be thought of as an 

alkali atom with a spin zero nucleus and a total mass m. The magnetic moment is M   S  

, 

where  is the gyromagnetic ratio of the level considered and S



 is the spin angular momentum of 

the atom. Since a spin 1⁄2 system is being considered, the spin angular momentum is proportional 

to the Pauli matrices, . 2

S     

The Schrödinger equation for the guided atom eigenstate E is 

     

2

2

2 E EE r SBr r E r

M

       

     

.

 (3) 

Equation 3 is made dimensionless by choosing a length scale such that x  x , where the new x

is dimensionless. Then    

2 22

2 2

2 2 E E

k

r r

M M

       and the recoil energy is defined as 

2 2

2 R

k E

M   , where 1 k   is the wavenumber of the optical field used to cool the atom. This is a 

rather arbitrary choice at this stage, but the recoil energy is a convenient scale for atoms that 

have been laser cooled. After dividing through by the recoil energy, pulling out the relevant 

factors, and defining the dimensionless transverse field parameter 1

1 2 R

B b

E  

  , the longitudinal 

field parameter 0

0 2 R

B b

E  

 , and the scaled energy 

R

E

E

  . 

2.2 Eigenstates Common to H and Λz 

The common eigenstates of H and z are found using the eigenstates of Lz, Sz in the following 

way 

, ,

, ,. z

H    

  



  (4) 

When the longitudinal bias field, 0 b is zero the Hamiltonian is invariant under various symmetry 

operations. This symmetry is the cause of a degeneracy of the eigenstates at 0 b  0 . The two 
6

degenerate states for each value of energy have opposite alignment. The operation Py x  takes an 

eigenstate , to its degenerate partner , , 

, , Py x     , (5) 

where Py takes y to –y and  x exchanges the spinor components when working in the position 

representation. In this representation, the spinor eigenfunction is given by   r r ,        . 

The component form is 

   

 

 

 

/2

/2

,

,

i

i

i

R e

r e

R e



  

 

 

  

  

 



 

   

         

 (6) 

and the degenerate partner is found as follows 

   

 

 

 

 

/2

/2

/2

/2

( )

i

i

yx yx i

i

i

i

R e

P r Pe

R e

R e

e

R e

r



 

 





 





 



 













 







 

   

 

    

 

 





. (7) 

Once the eigenstate is found for positive alignment, the degenerate partner of negative alignment 

is easily constructed. For this reason, positive alignment is assumed throughout the calculations 

when 0 b is zero and degeneracy is not an issue. 

When the assumed form (equation 6) is used in the Schrödinger equation (equation 3) the 

angular dependence drops out leaving only the following radial equations for the specific 

alignment  , 

 

 

2

2

2 0 1

2

2

2 0 1

1 1/2

() ,

1 1/2

() .

L R R R R bR b R

L R R R R bR b R

 

 



 

 



 

 

     

     



      

       

 (8) 

These two equations define the differential operators L and L that will be used in section 3. In 

the case of a spin half particle with integral orbital angular momentum, the pair must be solved 

for the allowed half-integral values of the alignment 135 , ,

222      , although the solutions 

for negative alignment are simply derived using symmetry when the bias field, 0 b is zero. 
7

3. Solution by the Frobenius Technique 

The series solution of differential equations near singular points requires the use of the Frobenius 

method (Ince, 1956; Bender, 1978). This is essentially a generalized power series approach that 

allows an extra degree of freedom for proper handling of the system behavior near regular 

singular points. Even though the differential equation contains singularities, the Frobenius 

method produces some solutions that are well behaved at the singular points and some solutions 

that have logarithmic or other types of singular behavior that do not usually satisfy the physical 

requirements of a desired solution. 

 The differential operator form of the coupled radial equations is useful for understanding the 

series solution technique that follows. The operators L/ defined in equation 8 are explicitly 

written as 

 

 

2

2

2 0

2

2

2 0

1 1/2

( ),

1 1/2

( ).

L b

L b

 

 





 





 







     

      

 (9) 

   

Using operator notation the coupled system is just written as the pair of differential equations 

1

1

,

.

LR b R

LR b R





  

  





 (10) 

3.1 Uncouple to Produce Fourth Order Differential Equations 

In order to apply the Frobenius method to this system, the equations are uncoupled producing 

two fourth-order equations that are central to this technique. By solving for both R R and  

directly on the right-hand side in the above system (equation 10) and substituting back into the 

other equation, we find the pair of fourth-order differential equations 

1 ,

1

1 .

1

1

1

L L R bR

b

L L R bR

b









  

  

    

 

    

 

 (11) 
8

These equations are uncoupled in the sense that only one unknown function is involved in each. 

Both equations are fourth order, homogeneous linear differential equations and therefore four 

independent solutions are required for the general solution of each. Both equations have a regular 

singularity at the origin and an irregular singularity at infinity. Equations of this type are called 

Hamburger equations (Ince, 1956). 

The general solution of each equation in equation 11 is actually only the general solution for one 

component of the two-component spinors defined in equation 7. The other component is found 

by substitution back into the original system (equation 10). Thus there are two distinct ways of 

finding the general solution of the system (equation 10). One could either solve the first equation 

in equation 11 for R and use the second equation of equation 10 to calculate the corresponding 

R components or use the reverse process to get R from R . Because the Frobenius method 

always produces a log-free solution for the largest indicial index found, these two reversed 

solution processes determine the two distinct types of physically acceptable modes in this 

system. 

3.2 Frobenius Series Approach 

The Frobenius series technique is applied by assuming a modified power series solution of the 

form 

1

1

0

n

p n

n

R cp  

 



  , (12) 

where the coefficients n cp and the index  1 are unknown except that 0 cp  1. In both the 

subscript p1and in the coefficient cp , the p refers to the fact that this solution process starts 

with the spin-up (or plus) radial component. The number 1 in the subscripts labels these solutions 

as modes of type 1. The reverse solution process will produce modes of type 2 and the letter m

will indicate a spin-down or (minus) component. Setting R R   p1 in the lower equation of 

equation 11 results in 

1 1

1

11 1

1

2 1

1

0 0

2

1 1

0

1

1

1 0

p p

n n

n n

n n

n

n n

n

L L R bR

b

cp L L b cp

cp L L b cp

 







 







 

    

 

 

 

 



    

 

      

       

 



. (13) 

The result of letting the operator 1 L L

   act on 1 n    is 
9

    1 1 1 3 15 1

3 5

1 , , n n P Pn Pn

L L      

  

 

 

             

 (14) 

and the functions Pi are defined here 

  

   

   

  

22 2

3 0 0 00

43 2

2 2

2 2

5

32 2 2

43 2 2 2 4 2

1 0

( , ) 2 4 2 2 2 2 2 /2 2 2

( , ) 4 6 17 / 2 18 2 6

4 18 17 4 3 / 2 6 6

35 6 2 17 / 2 3 / 2 6 6 9 19 / 2

16

Pn n b n b n b b

Pn n n n

n

P b

     

   

    

  



   

    

          

      

     

        







 (15). 

Inserting equation 14 into equation 13 produces 

    1 1 11 135 2

1 3 1 5 1 11

0

,, 0 n n nn

nn n n

n

cp P cp P n cp P n b cp       

    





    . (16) 

By adjusting the indices in the terms of equation 16, an equation in which each term of the sum 

contains only one power of  can be obtained. In this form, the coefficients n cp are defined to be 

zero for negative n and it is easy to pick out both the indicial equation and the recurrence 

relations for the coefficients, 

     1 2 5

5 1 23 1 41 1 6

0

0

, 2, 0

0 for 0

1

n

n n nn

n

n

cp P n cp P n cp P b cp

cp n

cp

  

  

 



   

 





. (17) 

This equation is solved by setting the factor multiplying 1 n  5    to zero for each n . This 

produces the following system of equations that defines both the indicial equations and the 

recurrence relations for the coefficients, 

    2

5 1 23 1 41 1 6 , 2, 0, 0 n n nn cp P n cp P n cp P b cp n            . (18) 

The first several of these equations have been written out explicitly here, 
10

 

 

 

  

 

 

 

 

05 1

15 1

25 1 03 1

35 1 13 1

45 1 23 1 01

5 5 1 3 3 1 11

2

65 1 43 1 21 1 0

75 1

0, 0, 0,

1, 1, 0,

2, 2, 0, 0,

3, 3, 1, 0,

4, 4, 2, 0,

5, 5, 3, 0,

6, 6, 4, 0,

7, 7,

n cp P

n cp P

n cp P cp P

n cp P cp P

n cp P cp P cp P

n cp P cp P cp P

n cp P cp P cp P b cp

n cp P





 

 

 

 

 



 

 

 

 

  

  

   

  

 

2

53 1 31 1 1

2

85 1 63 1 41 1 2

5, 0,

8, 8, 6, 0,

cp P cp P b cp

n cp P cp P cp P b cp



 

  

    

 (19) 

Since 0 cp  1, the n  0 equation determines the allowed values of  1 in terms of the alignment. 

This is the indicial equation and is given by 5 1 P (0, ) 0   or

   43 2 2 2 4 2

11 1 1

35 6 2 17 / 2 3 / 2 6 6 9 19 / 2 0

16                     (20) 

By inspecting equation 19, it can be determined that all of the odd n coefficients n cp are zero 

since P5 1   1, 0   and all odd coefficients are proportional to 1 cp . 

The roots of the indicial equation in equation 20 determine the form of the solutions for Rp1

generated by the Frobenius technique. The four roots are 

1       [ 1/ 2, 7 / 2, 5 / 2, 1/ 2]    . (21) 

The allowed values for  are positive half-integers so all of the values of 1 are integers. In the 

Frobenius method, repeated roots and roots separated by integers generally introduce logarithmic 

singularities into the solutions and negative values of  1 produce solutions with poles at the 

origin. All of the singular solutions must be excluded in the physical problem. In the Frobenius 

method, the solution corresponding to the largest indicial root is always free of logarithmic 

terms. For   1/2 , the largest allowed value of  1 is always   5/2 . 

This calculation is essentially repeated for the type 2 modes by solving the first of equation in 

equation 11 using 

2

2

0

n

m n

n

R R cm  

 





   (22) 

to develop the slightly different indicial equation for  2 , 

   43 2 2 2 4 2

22 2 2

35 6 2 17 / 2 3 / 2 6 6 9 19 / 2 0

16                     . (23) 
11

This indicial equation can be obtained directly from equation 20 by simply changing the sign of 

 and the solutions are easily obtained in the same way from equation 21, 

2        [ 5 / 2, 1/ 2 , 7 / 2,1/ 2 ]    . (24) 

However, the largest index for positive  is now   7/2 and since the largest index is always 

log-free, the solution (equation 22) using this value of  2 is the second physically acceptable 

solution of the system and is referred to as mode 2. 

It is straightforward to see that for   7/2 there are at least two solutions to equation 11 that are 

unacceptable because of the   singularity at the origin. For smaller values of  , it is not as 

obvious because there are some cases in which three values of  are unique integers. In these 

cases, the three solutions need to be checked more carefully to be certain they are not all logfree. We have checked this using the computer algebra program, Maple, and it turns out that two, 

and only two, log-free modes can be found for any value of  in this system. 

Once the allowed indices have been determined, the coefficients of the series expansions 

equations 12 and 22 must be found to complete the solution. The n cp are determined directly 

from equation 18 as 

 

 

2

23 1 41 1 6

5 1

2, 2 and even ,

n nn

n

cp P n cp P b cp

cp n P n





       (25) 

and a similar expression is obtained by the same technique for the n cm . The recurrence relations 

defining the two physically acceptable solutions, determined by 1     5/2 for the n cp and 

2     7/2 for the n cm are 

   

   

    

   

2 2

1 60 4 0 2

2 22

1 60 4 0 2

0

2

0

(2 ) , 22 32 1

2 12 , 4 2 32 1

1, 1,

0, 0 for all even and for all 0

2 nn n

n

nn n

n

n n

b cp b cp n b cp

cp

nn n n

b cm b cm n n b cm

cm

nn n n

cp cm

cp cm

n

n n

 

 

  







 

 

     



 

       



 

 

. (26) 

Using these recurrence relations, the following power series solutions for the two types of 

physical modes can be used directly for calculating eigenvalues, eigenfunctions, and the quantum 

behavior of the quadrupole waveguide. The mode 1 solutions for both spinor components are 

    0

5/2

1

0

1/2

3/2

1

2

1 0 1 0

4 μ 2 μ 6 5 

,

 

n

p n

n

n n

m n n

n n

R cp

b

R cp cp n nn

b b





 



 

  

 



   

 



 



  



 

 (27) 
12

and the equivalent solutions for the mode 2 spinor components are 

  2   0

7/2

2

0

1/2

5/2

2

1 1 0 0

8 μ 2 μ

,

 12 57 

n

m n

n

n n

p n n

n n

b

n nn

b

R cm

R cm cm

b







 

 



 

 



   

 

    







 

. (28) 

4. Series Evaluations Using Arbitrary Precision Arithmetic 

The waveguide modes are completely defined by the coefficients in equation 26 and the series 

solutions presented in equations 27 and 28. In order to find the eigenvalues, the series solutions 

are evaluated at some large value of  and  is adjusted until the bound component of the spinor 

at that large radius is near zero. The bound component is just the combination R R p  m (Hinds, 

2000; Golding, 2009) for the specific mode type, either 1 or 2. Asymptotic analysis shows that 

the bound components die off exponentially, but there is a small quickly decaying oscillatory 

component that can have an effect on the eigenvalues calculated using this technique and this 

possibility must be studied further. 

The series solutions converge but they cannot be easily summed using the standard floating point 

hardware available in most computers. The solution adopted here is to use an existing 

implementation of arbitrary precision arithmetic for these calculations. The computer algebra 

program Maple makes use of the GNU Multiple Precision (GMP) Library (GNU, 2010) to 

perform the basic arithmetic calculations addition, subtraction, division, multiplication, and 

exponentiation. These operations are all that are needed to sum the power series solutions in 

equations 27 and 28 to perform the calculations reported in section 4.1. 

4.1 Series Calculations 

The basic difficulty in direct summation of the series solutions in equations 27 and 28 is not that 

the series diverge or that they do not converge quickly, it is really just that the intermediate terms 

in the sum are very large and cannot be represented accurately in hardware floating point. 

Arbitrary precision arithmetic techniques eliminate this constraint at the expense of some speed 

and efficiency. The tradeoff is extremely useful when extreme accuracy is needed and other 

techniques are suspect. 

In order to do these calculations effectively, the number of digits needed to represent the largest 

number in the sum to at least the precision needed in the final result should be estimated. In 

addition, the total number of terms needed for convergence at the desired accuracy needs to be 

known. For example, in the solution of an eigenvalue problem, the value of an eigenfunction 

might be required to 20 significant digits at   40 . In order to guarantee that the sums will be 

accurately calculated, individual terms of the series should be calculated for a reasonable test 
13

case and the maximum term value obtained. If the maximum value is 300 10 and the desired 

precision of the zero is 20 10 then at least 320 digits are needed to avoid losing precision in the 

final sum due to a loss of precision of any single term in the series. In practice, some guard digits 

should be added and convergence should be checked by changing the number of digits and the 

number of terms used in the calculation. 

Both the number of digits and the number of terms required to evaluate a particular sum 

accurately can be visualized by creating a term plot. This is just a plot of the log base 10 of the 

magnitude of the terms in the sequence to be summed. The term plot for the series defining Rp1

is shown in figure 2. The summation calculated is the first sum in equation 27. The terms plotted 

are just the individual 10 log n

n cp  for n from zero to the largest value needed. The number of 

digits needed for the sum is just the difference between the maximum value and the minimum 

value assuming the plot is terminated when the term with the desired precision is reached. 

. 

Figure 2. Term plot of the log of the magnitude of the individual terms 

used in a calculation of the first series in equation 27. The target 

value of radius is 40 and the maximum eigenvalue is 30. The 

largest term in the sequence is roughly 10100 and therefore 

120 digits needs to be used to do an adequate job of summing 

this series for any value of radius up to 40 for eigenvalues less than 30. 

If the number of required digits is underestimated, the calculation will breakdown numerically at 

values of the radial coordinate smaller than the desired target value. The results of using different 

numbers of digits in the summations are shown in figures 3, 4, and 5. In figure 3, the number of 

digits used is 15, which is approximately the number of digits used in floating point hardware 

calculations on many machines. The calculation breaks down at around  12 . The range of the 
14

calculation can be extended significantly by using 25 digits, as shown in figure 4. In figure 5, 

100 digits were used in the calculation and there are no obvious problems with the calculations 

over the whole range from zero to 40. In figure 6, the tail of the 100-digit calculation is blown up 

by a factor of 10,000 and a very clean but heavily damped oscillation is observed. This particular 

oscillation can be predicted by asymptotic techniques but it has been difficult to verify in 

previous calculations performed at hardware precision. Understanding this behavior in detail can 

help avoid certain types of systematic errors in the eigenvalue calculations that involve 

accurately locating zeros in the tails of the eigenfunctions. 

Figure 3. Calculation by series summation of the lowest two bound components using 15 digits of 

precision. This kind of result is representative of hardware floating point. Beyond a certain point 

the calculation breaks down. As the number of digits used in the calculation is increased, the 

target radius can be significantly increased. In this way, hardware floating point limits the 

ability to obtain accurate eigenvalues. 
15

Figure 4. The same calculation as in figure 3. The only difference is that 25 digits of precision are used. 

One can see that the calculation can proceed to much larger target values as the number of digits 

is increased. 

Figure 5. The same calculation as displayed in figures 3 and 4, except that 100 digits were used in this 

calculation. The effective target radius can be pushed out to 40 with very clean results. This 

technique lets us calculate significantly more accurate eigenvalues as well as perform checks 

on asymptotic series expansions for future work. 
16

Figure 6. This is a blown up version of the tail of the 100-digit calculation displayed in figure 5. The 

plot is from radius of 0 to 40 and the scale is expanded by a factor of 10,000. Very clean 

damped oscillations are observed. This is an example of the power of the basic arbitrary 

precision techniques. In the same region where we saw complete breakdown of the 

calculations at 15 and 25 digits in figures 3 and 4, extremely useful information about the 

asymptotic behavior is obtained by simply summing the series using 100 digits of 

precision. This cannot be obtained directly using hardware floating point techniques. 

4.2 Eigenvalues: Shooting or Summing? 

In order to solve eigenvalue problems of the type studied here, a variation of the shooting 

method is used. The shooting method is the technique in which a differential equation is 

integrated from a boundary point where the solution is known to a target point where the solution 

should attain a target value. In practice, the eigenvalue or some other parameter can be adjusted 

so that the unknown function has the desired behavior at the target point. When the solution is 

acceptable, the eigenvalue or desired value of the parameter has been found. One limitation of 

this technique is that the differential equation solver must be able to integrate the system 

accurately from the known point to the target point. This can be difficult if the system is stiff or 

unstable or otherwise very sensitive to the exact value of the eigenvalue. 

To determine the eigenvalues using the series approach, a variation of the shooting technique is 

used. In this approach, the solution does not depend on numerical solutions of the differential 

equations from an initial point to a target point. The only requirement is to accurately sum the 

series at the target point for various trial eigenvalues. The initial condition is built in to the series 

solution. The eigenvalue must be varied and the summation repeated until the target value for the 

unknown function is determined. This type of effective shooting technique can be accomplished 

over a large range of distances in systems that support arbitrary precision arithmetic and have 

series solutions. 
17

Because the sum of the two components R R p  m is an apparent bound state for both mode types 

1 and 2 (Hinds, 2000; Golding, 2009), the leading behavior of this bound component dies off 

exponentially with distance from the center of the guide. Using this fact, the target behavior for 

the eigenvalue solution is just that 0 as R R p m     . The eigenvalues are then readily 

found using a root solving routine that uses arbitrary precision arithmetic to look for zeros of the 

sums defined in equations 27 and 28 at large radii. We have found that in this way the target 

point can be much farther from the origin than it can using the similar ODE based shooting 

techniques. This results in much more precise calculations of the eigenvalues. This enhanced 

resolution is expected to be useful as these functions are used for further calculations on the 

detailed properties of atom guides. 

The results of eigenvalue calculations for 0   1/ 2, 0 b  and b1 1  are shown in figure 7. The 

first 52 eigenvalues are shown. The eigenvalue data can be approximated by the simple equation

0.7355 1.893790658 1.270128257  n . In table 1, the first 26 values of each type of mode are 

presented. Careful inspection shows that eigenvalues form an interleaved series in the sense that 

the first eigenvalue is a type 1 mode and second is a type 2 mode; this pattern holds up as far as 

we have checked. In figure 7, mode type 1 points are red and mode type 2 points are blue. 
18

Figure 7. Plot of the first 52 eigenvalues. The red points are mode type 1 eigenvalues and the blue points are 

mode type 2 eigenvalues; an approximate fit to the eigenvalues as a function of n is given in the 

text. 
19

Table 1. First 52 energy eigenvalues for our system with alignment 1/2, 

zero bias field and a transverse gradient of b1=1. The sequence of 

eigenvalues alternates between type 1 and type 2 eigenvalues as n 

is increased. Because of this double counting the index n in the 

table is not the same as the index n used in the plot of figure 7. 

n Type 1 Eigenvalue Type 2 Eigenvalue 

1 1.893790658 2.770959115 

2 3.735223647 4.436734695 

3 5.211244681 5.827843267 

4 6.504703191 7.067475471 

5 7.681587325 8.20585412 

6 8.774989674 9.269728735 

7 9.80443751 10.27548299 

8 10.78265521 11.23407047 

9 11.71851157 12.15327128 

10 12.61849734 13.03886265 

11 13.48754237 13.89528345 

12 14.32950179 14.72603767 

13 15.14746211 15.53395267 

14 15.9439429 16.32135181 

15 16.72103469 17.09017376 

16 17.48049623 17.84205754 

17 18.22382493 18.57840462 

18 18.952309 19.30042528 

19 19.66706684 20.00917379 

20 20.36907731 20.70557569 

21 21.05920333 21.39044912 

22 21.73821055 22.06452178 

23 22.40678229 22.72844459 

24 23.06553163 23.38280277 

25 23.7150113 24.02812491 

26 24.35572179 24.66489048 

4.3 Eigenfunctions 

Once the eigenvalues are determined, the corresponding eigenfunctions are found by a 

straightforward summation of the series solutions equations 27 and 28 at any desired value of  . 

The eigenfunctions determined in this way display detailed behavior in regions beyond those 

accessible using the hardware precision ODE shooting techniques. This type of behavior is 

shown in figures 8 and 9. It is clear that the bound component does not simply decay away 

exponentially but has a small and quickly decaying oscillatory component as well. This behavior 

is explained by a detailed asymptotic solution of the problem. The ability to compare the direct 

calculations of the series technique with the asymptotic solutions over wide ranges will be very 

helpful in checking various aspects of the asymptotic behavior of these solutions. 
20

Figure 8. The first 10 type 1 mode radial eigenfunctions. The black curve with no zero crossings is the ground state 

wave function. The first excited state is the blue curve. 

Figure 9. The first 10 type 2 mode radial eigenfunctions. The black curve is the lowest energy type 2 mode but is 

actually the first excited state of the system. 
21

5. Conclusions 

The use of series solution techniques along with arbitrary precision arithmetic is a very powerful 

combination for investigating the solutions of differential equations in regions where ordinary 

hardware precision numerical integration techniques break down. The use of arbitrary precision 

in differential equation solvers is implemented in Maple and these can be made to produce 

accurate solutions at large radii. However, the solution of a differential equation by the series 

technique simply requires the appropriate power series formulation of the solution and a software 

system that implements an arbitrary precision library, such as GMP (GNU, 2010). By using a 

judicious choice of the number of digits and terms needed for a particular set of calculations, the 

solutions can be calculated very quickly and with excellent resolution. We have found this 

technique to be extremely useful for calculating radial wavefunctions and eigenfunctions for the 

atomic waveguide problem as it produces accurate solutions in situations where spinor equations 

can be uncoupled to produce higher order differential equations that are difficult to solve 

analytically. 
22

6. References 

Bender, C. M.; Orszag, S. A. Advanced Mathematical Methods for Scientists and Engineers;

McGraw-Hill: New York, NY, 1978. 

Bergman, T. H.; et al. JOSA B 1989, 6 (11), 2249. 

Bill, J. Phys. Rev. A 73 2006, 053609. 

Blümel, R.; Dietrich, K. Phys. Rev. A 1991, 43, 22–28. 

Fortágh, J.; Zimmerman, C. RMP 2007, 79, 235. 

GNU Multiprecision Library documentation located at http://gmplib.org/ (accessed 2010). 

Golding, W. M. Atomic Waveguides for Atom Chips; ARL-TR-5014; U.S. Army Research 

Laboratory: Adelphi, MD, Oct. 2009. 

Hinds, E. A.; Eberlein, C. Phys. Rev. A 2000, 61, 033614. 

Hinds, E. A.; Eberlein, C. Phys. Rev. A 2001, 64, 039902(E). 

Ince, E. L. Ordinary Differential Equations; Dover: New York, NY, 1956, p. 396. 

Lesanovsky, I.; Schmelcher, P. Phys. Rev. A 2004, 70, 063604. 

Potvliege, R. M.; Zehnlé, V. Phys. Rev. A 2001, 63, 025601. 

Sukumar, C. V.; Brink, D. M. Phys. Rev. A 1997, 56, 2451. 
23

NO. OF 

COPIES ORGANIZATION 

1 PDF DEFENSE TECH INFO CTR 

 ATTN DTIC OCA 

 8725 JOHN J KINGMAN RD STE 0944 

 FT BELVOIR VA 22060-6218 

1 DARPA 

 ATTN J ABO-SHAEER 

 3701 N FAIRFAX DR 

 ARLINGTON VA 22203-1714 

1 CD OFC OF THE SECY OF DEFNS 

 ATTN ODDRE (R&AT) 

 THE PENTAGON 

 WASHINGTON DC 20301-3080 

1 US ARMY RSRCH DEV AND ENGRG 

 CMND 

 ARMAMENT RSRCH DEV & ENGRG 

 CTR 

 ARMAMENT ENGRG & TECHNLGY 

 CTR 

 ATTN AMSRD AAR AEF T 

 J MATTS 

 BLDG 305 

 ABERDEEN PROVING GROUND MD 

 21005-5001 

1 PM TIMS, PROFILER (MMS-P) 

 AN/TMQ-52 

 ATTN B GRIFFIES 

 BUILDING 563 

 FT MONMOUTH NJ 07703 

1 US ARMY INFO SYS ENGRG CMND 

 ATTN AMSEL IE TD A RIVERA 

 FT HUACHUCA AZ 85613-5300 

1 COMMANDER 

 US ARMY RDECOM 

 ATTN AMSRD AMR 

 W C MCCORKLE 

 5400 FOWLER RD 

 REDSTONE ARSENAL AL 35898-5000 

3 DEPT OF THE ARMY 

 WEAPONS SCIENCES 

 DIRECTORATE 

 ATTN AMSRD AMS WS 

 ATTN AMSRD AMS WS H EVERITT 

 ATTN AMSRD AMS WS K MYNENI 

 ATTN AMSRD AMS WS T BAHDER 

 REDSTONE ARSENAL AL 35898-5000 

NO. OF 

COPIES ORGANIZATION 

1 US GOVERNMENT PRINT OFF 

 DEPOSITORY RECEIVING SECTION 

 ATTN MAIL STOP IDAD J TATE 

 732 NORTH CAPITOL ST NW 

 WASHINGTON DC 20402 

1 US ARMY RSRCH LAB 

 ATTN RDRL CIM G T LANDFRIED 

 BLDG 4600 

 ABERDEEN PROVING GROUND MD 

 21005-5066 

1 HD US ARMY RSRCH OFC 

1 CD ATTN P REYNOLDS 

 PO BOX 12211 

 RESEARCH TRIANGLE PARK NC 

 27709-2211 

29 US ARMY RSRCH LAB 

 ATTN IMNE ALC HRR MAIL & 

 RECORDS MGMT 

 ATTN RDRL CIM L TECHL LIB 

 ATTN RDRL CIM P TECHL PUB 

 ATTN RDRL SEE O N FELL 

 ATTN RDRL SEE O P LEE 

 ATTN RDRL SEE O P PELLEGRINO 

 ATTN RDRL SEE O W M GOLDING 

 (20 COPIES) 

 ATTN RDRL SER E F CROWNE 

 ATTN RDRL SES E H BRANDT 

 ATTN RDRL SES P A EDELSTEIN 

 ADELPHI MD 20783-1197 

TOTAL: 43 (40, 2 CDS, 1 PDF) 
24

INTENTIONALLY LEFT BLANK