Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want to program in the regular and irregular spherical Coulomb wavefunctions $F_\ell(\gamma,kr)$ and $G_\ell(\gamma,kr)$, respectively, which are defined in terms of the regular and irregular solutions to Kummer's differential equation

$$\Big[z\frac{d^2}{dz^2}+(\beta-z)\frac{d}{dz}-\alpha\Big]\,f(z)=0\,.$$

The regular solution is given by the standard Confluent hypergeometric function ${}_1F_1(\alpha|\,\beta|\,z)$. I can access this function in Mathematica using Hypergeometric1F1, out of which I can construct the regular spherical Coulomb wavefunction $F_\ell(\gamma,kr)$.

But I am unable to find the two irregular solutions, typically denoted $W_1$ and $W_2$. I need these to construct the irregular spherical Coulomb wavefunction. The $W$'s and the regular ${}_1F_1$ are connected by the relations:

$${}_1F_1(\alpha|\,\beta|\,z)=W_1(\alpha|\,\beta|\,z)+W_2(\alpha|\,\beta|\,z)$$

Does anyone know how I can access both irregular solutions?

share|improve this question
2  
Mathworld says the second solution is implemented as HypergeometricU. You can also uncover this by solving the equation directly, DSolve[z f''[z] + (\[Beta] - z) f'[z] - \[Alpha] f[z] == 0, f[z], z]. –  whuber Mar 29 '13 at 20:04
add comment

2 Answers 2

up vote 3 down vote accepted

The formulae here also allow for an expression in terms of the irregular Whittaker function WhittakerW[]:

CoulombPhase[ℓ_, η_] := (LogGamma[1 + ℓ + I η] - LogGamma[1 + ℓ - I η])/(2 I)

CoulombHPlus[ℓ_, η_, ρ_] :=
                 Exp[π (η - I ℓ)/2 + I CoulombPhase[ℓ, η]] WhittakerW[-I η, ℓ + 1/2, -2 I ρ]

CoulombHMinus[ℓ_, η_, ρ_] :=
                   Exp[π (η + I ℓ)/2 - I CoulombPhase[ℓ, η]] WhittakerW[I η, ℓ + 1/2, 2 I ρ]

Compare with OP's implementation (which has a different argument convention):

N[{CoulombHPlus[2/3, 3/4, 5], CoulombHplus[2/3, 4/3, 5 3/4]}]
   {-0.93049 + 0.597603 I, -0.93049 + 0.597603 I}

N[{CoulombHMinus[2/3, 3/4, 5], CoulombHminus[2/3, 4/3, 5 3/4]}]
   {-0.93049 - 0.597603 I, -0.93049 - 0.597603 I}

These are of course the complex-valued solutions (they are related to the real solutions in a way similar to the relationship of the Hankel functions with the Bessel functions); if one wants the real-valued irregular solution, this is easily done. For completeness, I will also include an implementation of the regular Coulomb function:

CoulombFactor[ℓ_, η_] :=
      Exp[-π (η + I (ℓ + 1))/2] Sqrt[Gamma[1 + ℓ - I η] Gamma[1 + ℓ + I η]]/(2 (2 ℓ + 1)!)

CoulombF[ℓ_, η_, ρ_] := CoulombFactor[ℓ, η] WhittakerM[I η, ℓ + 1/2, 2 I ρ]

CoulombG[ℓ_, η_, ρ_] := (CoulombHPlus[ℓ, η, ρ] + CoulombHMinus[ℓ, η, ρ])/2

Compare:

N[{(CoulombHPlus[2/3, 3/4, 5] - CoulombHMinus[2/3, 3/4, 5])/(2 I),
   CoulombF[2/3, 3/4, 5]}] // Chop
   {0.597603, 0.597603}

N[{(CoulombHPlus[2/3, 3/4, 5] + CoulombHMinus[2/3, 3/4, 5])/2,
   CoulombG[2/3, 3/4, 5]}] // Chop
   {-0.93049, -0.93049}

Replicate the plots here:

With[{ℓ = 0, η = Sqrt[15/2]}, 
     Plot[{CoulombF[ℓ, η, ρ], CoulombG[ℓ, η, ρ], Abs[CoulombHPlus[ℓ, η, ρ]]}, {ρ, 0, 30}, 
          PlotRange -> {-2, 5}, PlotStyle -> {Blue, Red, Green}]]

Coulomb wave function plots

With[{ℓ = 5, η = 0}, 
     Plot[{CoulombF[ℓ, η, ρ], CoulombG[ℓ, η, ρ], Abs[CoulombHPlus[ℓ, η, ρ]]}, {ρ, 0, 30}, 
          PlotRange -> {-2, 5}, PlotStyle -> {Blue, Red, Green}]]

Coulomb wave function plots

and the plots here:

(* from http://dlmf.nist.gov/help/vrml/aboutcolor *)
DLMFHeightColor[u_?NumericQ] := Blend[{Blue, Cyan, Green, Yellow, Red}, u] /; 0 <= u <= 1

Plot3D[CoulombF[0, η, ρ], {ρ, 0, 5}, {η, -2, 2}, BoundaryStyle -> None,
       ColorFunction -> DLMFHeightColor, Mesh -> False, PlotPoints -> 75,
       WorkingPrecision -> 20]

surface plot of regular Coulomb wave function

Plot3D[CoulombG[0, η, ρ], {ρ, 0, 5}, {η, -2, 2}, BoundaryStyle -> None,
       ClippingStyle -> None, ColorFunction -> DLMFHeightColor,
       Mesh -> False, PlotPoints -> 55, WorkingPrecision -> 20]

surface plot of irregular Coulomb wave function

I needed to use arbitrary precision evaluation for these plots since the evaluation of the Coulomb functions with these definitions at machine precision seemed unstable for some combinations of arguments (try them yourself). There should be a way to reformulate these formulae for better stability (there are already a number of papers on the numerical evaluation of Coulomb functions), but I don't have time to look through those at the moment.

share|improve this answer
add comment

I stumbled upon the NIST Digital Library of Mathematical Functions online, and the relations given in section 33.2 is sufficient to program the useful Coulomb wavefunctions into Mathematica. My implementation is as follows:

The Sommerfeld parameter:

γ[k_] = 1/k

The irregular Coulomb functions $H^{\pm}_\ell(\gamma,\,kr)$

coulθ = (k r - γ[k] Log[2 k r] - ℓ π/2 + Arg[Gamma[ℓ + 1 + I γ[k]]]);

CoulombHplus[ℓ_, k_, r_] = E^(I coulθ) (-2 I k r)^(ℓ + 1 + I γ[k])
   HypergeometricU[ℓ + 1 + I γ[k],  2 ℓ + 2, -2 I k r];
CoulombHminus[ℓ_, k_, r_] = E^(-I coulθ) (2 I k r)^(ℓ + 1 - I γ[k])
   HypergeometricU[ℓ + 1 - I γ[k], 2 ℓ + 2, 2 I k r];

Then the standard regular and irregular Coulomb functions follow: $F_\ell(\gamma,\,kr)$ and $G_\ell(\gamma,\,kr)$

CoulombF[ℓ_, k_, r_] = 1/(2 I) (CoulombHplus[ℓ, k, r] - CoulombHminus[ℓ, k, r]);
CoulombG[ℓ_, k_, r_] = 1/2 (CoulombHplus[ℓ, k, r] + CoulombHminus[ℓ, k, r]);

A good consistency check is to compare against their known asymptotic behavior as $r\rightarrow\infty$:

asympCoulombF[ℓ_, k_, r_] = Sin[k r - π ℓ/2 - γ[k] Log[2 k r] + Arg[Gamma[ℓ + 1 + I γ[k]]]];
asympCoulombG[ℓ_, k_, r_] = Cos[k r - π ℓ/2 - γ[k] Log[2 k r] + Arg[Gamma[ℓ + 1 + I γ[k]]]];
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.