Skip to content

Commit

Permalink
Add CoulombTwoBody
Browse files Browse the repository at this point in the history
  • Loading branch information
ohno committed Apr 20, 2024
1 parent b4c10b5 commit 334e930
Show file tree
Hide file tree
Showing 11 changed files with 2,866 additions and 286 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,6 @@ E(H)
Helium cation has symbol $\mathrm{He}^+$ and atomic number 2 ($Z=2$). Therefore the ground state ($n=1$) energy is $-2 E_\mathrm{h}$.

```julia
using Antique
He⁺ = HydrogenAtom(Z=2)
E(He⁺)
# output> -2.0
Expand All @@ -51,6 +50,7 @@ There are more examples on each model page.
- [Morse Potential](https://ohno.github.io/Antique.jl/stable/MorsePotential/) `MorsePotential`
- [Rigid Rotor](https://ohno.github.io/Antique.jl/stable/RigidRotor/) `RigidRotor`
- [Hydrogen Atom](https://ohno.github.io/Antique.jl/stable/HydrogenAtom/) `HydrogenAtom`
- [Coulomb 2-Body System](https://ohno.github.io/Antique.jl/stable/HydrogenAtom/) `CoulombTwoBody`

## Future Works

Expand Down
262 changes: 9 additions & 253 deletions developer/dev.ipynb

Large diffs are not rendered by default.

21 changes: 11 additions & 10 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,16 +21,17 @@ makedocs(;
),
pages=[
"Home" => "index.md",
"Delta Potential" => "DeltaPotential.md" ,
"Infinite Potential Well" => "InfinitePotentialWell.md",
"Harmonic Oscillator" => "HarmonicOscillator.md" ,
"Morse Potential" => "MorsePotential.md" ,
"Pöschl-Teller Potential" => "PoschlTeller.md" ,
"Hydrogen Atom" => "HydrogenAtom.md" ,
"Rigid Rotor" => "RigidRotor.md" ,
"Infinite Potential Well 3D" => "InfinitePotentialWell3D.md",
# "Harmonic Oscillator 3D" => "HarmonicOscillator3D.md" ,
# "API reference" => "API.md" ,
"Delta Potential" => "DeltaPotential.md" ,
"Infinite Potential Well" => "InfinitePotentialWell.md" ,
"Harmonic Oscillator" => "HarmonicOscillator.md" ,
"Morse Potential" => "MorsePotential.md" ,
"Pöschl-Teller Potential" => "PoschlTeller.md" ,
"Rigid Rotor" => "RigidRotor.md" ,
"Hydrogen Atom" => "HydrogenAtom.md" ,
"Coulomb 2-Body System" => "CoulombTwoBody.md" ,
# "Infinite Potential Well 3D" => "InfinitePotentialWell3D.md",
# "Harmonic Oscillator 3D" => "HarmonicOscillator3D.md" ,
# "API reference" => "API.md" ,
],
)

Expand Down
274 changes: 274 additions & 0 deletions docs/src/CoulombTwoBody.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,274 @@
```@meta
CurrentModule = Antique
```

# Coulomb 2-Body System

This is the model of two particles interacting through Coulomb forces such as positronium, muonium, hydrogen atoms, deuterium atoms, etc.

## Definitions

This model is described with the time-independent Schrödinger equation
```math
\hat{H} \psi(\pmb{r}) = E \psi(\pmb{r}),
```
and the Hamiltonian
```math
\hat{H} = - \frac{\hbar^2}{2\mu} \nabla^2 + V(r),
```
where $\mu=\left(\frac{1}{m_1}+\frac{1}{m_2}\right)^{-1}$ is the reduced mass of particle 1 and particle 2. The potential includes only Coulomb interaction and it does not include fine or hyperfine interactions in this model. Parameters are specified with the following struct.

#### Parameters
```@docs; canonical=false
Antique.CoulombTwoBody
```

#### Potential
```@docs; canonical=false
Antique.V(::CoulombTwoBody, ::Any)
```

#### Eigen Values
```@docs; canonical=false
Antique.E(::CoulombTwoBody)
```

#### Eigen Functions
```@docs; canonical=false
Antique.ψ(::CoulombTwoBody, ::Any, ::Any, ::Any)
```

#### Radial Functions
```@docs; canonical=false
Antique.R(::CoulombTwoBody, ::Any)
```

#### Associated Laguerre Polynomials
```@docs; canonical=false
Antique.L(::CoulombTwoBody, ::Any)
```

#### Spherical Harmonics
```@docs; canonical=false
Antique.Y(::CoulombTwoBody, ::Any, ::Any)
```

#### Associated Legendre Polynomials
```@docs; canonical=false
Antique.P(::CoulombTwoBody, ::Any)
```

## Usage & Examples

[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `CoulombTwoBody` and several parameters `z₁`, `z₂`, `m₁`, `m₂`, `mₑ`, `a₀`, `Eₕ` and `` are set as optional arguments.

```@example CTB
using Antique
Ps = CoulombTwoBody(z₁=-1, z₂=1, m₁=1.0, m₂=1.0, mₑ=1.0, a₀=1.0, Eₕ=1.0, ℏ=1.0)
; #hide
```

#### Parameters

```@repl CTB
Ps.z₁
Ps.z₂
Ps.m₁
Ps.m₂
Ps.mₑ
Ps.a₀
Ps.Eₕ
Ps.ℏ
```

#### Eigen Values

Examples of [positronium](https://en.wikipedia.org/wiki/Positronium#Energy_levels):

```@repl CTB
E(Ps, n=1)
E(Ps, n=2)
```

#### Mass and Charge Dependence

The values of masses are cited from the 2018 CODATA recommended values, [E. Tiesinga, et al., Rev. Mod. Phys. 93, 025010 (2021)](https://doi.org/10.1103/RevModPhys.93.025010).

```@example CTB
me = 1.0 # me #
mµ = 206.7682830 # me # https://physics.nist.gov/cgi%2Dbin/cuu/Value?mmusme
mp = 1836.15267343 # me # https://physics.nist.gov/cgi%2Dbin/cuu/Value?mpsme
md = 3670.48296788 # me # https://physics.nist.gov/cgi%2Dbin/cuu/Value?mdsme
mt = 5496.92153573 # me # https://physics.nist.gov/cgi%2Dbin/cuu/Value?mtsme
mh = 5495.88528007 # me # https://physics.nist.gov/cgi-bin/cuu/Value?mhsme
ma = 7294.29954142 # me # https://physics.nist.gov/cgi-bin/cuu/Value?malsme
Ps = CoulombTwoBody(m₁=me, m₂=me)
Mu = CoulombTwoBody(m₁=me, m₂=mµ)
H = CoulombTwoBody(m₁=me, m₂=mp)
D = CoulombTwoBody(m₁=me, m₂=md)
T = CoulombTwoBody(m₁=me, m₂=mt)
BO = CoulombTwoBody(m₁=me, m₂=Inf)
He3⁺ = CoulombTwoBody(z₁=-1, z₂=+2, m₁=me, m₂=mh)
He4⁺ = CoulombTwoBody(z₁=-1, z₂=+2, m₁=me, m₂=ma)
He∞⁺ = CoulombTwoBody(z₁=-1, z₂=+2, m₁=me, m₂=Inf)
pμ = CoulombTwoBody(z₁=-1, z₂=+1, m₁=mµ, m₂=mp)
dμ = CoulombTwoBody(z₁=-1, z₂=+1, m₁=mµ, m₂=md)
tμ = CoulombTwoBody(z₁=-1, z₂=+1, m₁=mµ, m₂=mt)
bμ = CoulombTwoBody(z₁=-1, z₂=+1, m₁=mµ, m₂=Inf)
hμ = CoulombTwoBody(z₁=-1, z₂=+2, m₁=mµ, m₂=mh)
αμ = CoulombTwoBody(z₁=-1, z₂=+2, m₁=mµ, m₂=ma)
println(" \tE / Eₕ")
println("Ps \t", E(Ps))
println("Mu \t", E(Mu))
println("H \t", E(H))
println("D \t", E(D))
println("T \t", E(T))
println("∞H \t", E(BO))
println("³He⁺\t", E(He3⁺))
println("⁴He⁺\t", E(He4⁺))
println("∞He⁺\t", E(He∞⁺))
println("pμ \t", E(pμ))
println("dμ \t", E(dμ))
println("tμ \t", E(tμ))
println("bμ \t", E(bμ))
println("hμ \t", E(hμ))
println("αμ \t", E(αμ))
```

```@example CTB
println(" \t<δ³(r)> / a₀⁻³")
println("1/8π =\t", 1/8/π)
println("Ps \t", abs(ψ(Ps,0,0,0))^2)
println("Mu \t", abs(ψ(Mu,0,0,0))^2)
println("H \t", abs(ψ(H ,0,0,0))^2)
println("D \t", abs(ψ(D ,0,0,0))^2)
println("T \t", abs(ψ(T ,0,0,0))^2)
println("∞H \t", abs(ψ(BO,0,0,0))^2)
println("1/π = \t", 1/π)
```

#### Lifetime of Positronium

The lifetime $\tau$ of positronium (Ps, $\mathrm{e}^+\mathrm{e}^-$) is written as

```math
\tau = \frac{1}{\Gamma},
```

```math
\Gamma = 4 \pi \alpha^4 c {a_0}^2 \langle\delta^3(\pmb{r})\rangle.
```

where $\langle\delta^3(\pmb{r})\rangle = \langle\psi|\delta^3(\pmb{r})|\psi\rangle = |\psi(\pmb{0})|^2 = \frac{1}{8\pi} a_0^{-3} \simeq 2.685\times10^{29}~\mathrm{m}^{-3}$ is the value of probability density at the origin ($r=0$). Reference:
- (7.169) in D. J. Griffiths, Introduction to Elementary Particles (John Wiley & Sons, Inc. 1987) ISBN 0-471-60386-4
- [S. Berko, H. N. Pendleton, Annual Review of Nuclear and Particle Science, 30, 543 (1980)](https://doi.org/10.1146/annurev.ns.30.120180.002551)
- [A. M. Frolov, S. I. Kryuchkov, and V. H. Smith, Jr., Phys. Rev. A, 51, 4514 (1995)](https://doi.org/10.1103/PhysRevA.51.4514)

```@example CTB
α = 7.2973525693e-3 # # https://physics.nist.gov/cgi-bin/cuu/Value?alph
c = 299792458 # m s-1 # https://physics.nist.gov/cgi-bin/cuu/Value?c
a₀ = 5.29177210903e-11 # m # https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
Ps = CoulombTwoBody(z₁=1, z₂=-1, m₁=1.0, m₂=1.0, mₑ=1.0, a₀=1.0, Eₕ=1.0, ℏ=1.0)
δ = abs(ψ(Ps,0,0,0))^2 * a₀^(-3)
Γ = 4 * π * α^4 * c * a₀^2 * δ
τ = 1/Γ
println("<δ> = ", abs(ψ(Ps,0,0,0))^2, " a₀⁻³")
println(" = ", δ, " m⁻³")
println("Γ = ", Γ / 1e9, " GHz")
println("τ = ", τ / 1e-12, " ps")
```

#### Hyperfine Splitting

The hyperfine splitting of hydrogen atoms is given as

```math
\Delta E (\mathrm{H}) = -\frac{2}{3} \mu_0 \gamma_\mathrm{p} \gamma_\mathrm{e} \hbar^2 \langle\delta^3(\pmb{r})\rangle
```

in [Griffiths(1982)](https://doi.org/10.1119/1.12733). This fomula is not available for deuterium (D) and positronium (Ps). Because of the different spin between the proton and the deuteron for D, the contribution of positron-electron pair annihilation for Ps. Note the definition of gyromagnetic ratio. The mass of protons is used for all nucleons and nuclei:

```math
\begin{aligned}
&\gamma_\mathrm{e} = \frac{-e}{2 m_\mathrm{e}} g_\mathrm{e}, &
&\gamma_\mathrm{e^+} = \frac{+e}{2 m_\mathrm{e}} g_\mathrm{e}, &
&\gamma_\mathrm{\mu} = \frac{-e}{2 m_\mathrm{\mu}} g_\mathrm{\mu}, \\
&\gamma_\mathrm{p} = \frac{+e}{2 m_\mathrm{p}} g_\mathrm{p}, &
&\gamma_\mathrm{d} = \frac{+e}{2 m_\mathrm{p}} g_\mathrm{d}, &
&\gamma_\mathrm{t} = \frac{+e}{2 m_\mathrm{p}} g_\mathrm{t}, &
&\gamma_\mathrm{h} = \frac{+2e}{2 m_\mathrm{p}} g_\mathrm{h}. &
\end{aligned}
```

The value of probability density at the origin is $\langle\delta^3(\pmb{r})\rangle = \langle\psi|\delta^3(\pmb{r})|\psi\rangle = |\psi(\pmb{0})|^2 \simeq \frac{1}{\pi} a_0^{-3} \simeq 2.148\times10^{30}~\mathrm{m}^{-3}$ in Mu, H, D and T. This values are very different in Ps, $^3\mathrm{He}^+$ and muonic hydrogen ($\mathrm{p\mu}$) due to the difference of reduced masses and charges. The energy can be converted to frequency (Hz) by $v = \Delta E / h$.

```@example CTB
a₀ = 5.29177210903e-11 # m # https://physics.nist.gov/cgi-bin/cuu/Value?bohrrada0
Eₕ = 4.3597447222071e-18 # J # https://physics.nist.gov/cgi-bin/cuu/Value?hr
ℏ = 1.054571817e-34 # J s # https://physics.nist.gov/cgi-bin/cuu/Value?hbar
me = 9.1093837015e-31 # kg # https://physics.nist.gov/cgi-bin/cuu/Value?me
mµ = 1.883531627e-28 # kg # https://physics.nist.gov/cgi-bin/cuu/Value?mmu
mp = 1.67262192369e-27 # kg # https://physics.nist.gov/cgi-bin/cuu/Value?mp
md = 3.3435837724e-27 # kg # https://physics.nist.gov/cgi-bin/cuu/Value?md
mt = 5.0073567446e-27 # kg # https://physics.nist.gov/cgi-bin/cuu/Value?mt
mh = 5.0064127796e-27 # kg # https://physics.nist.gov/cgi-bin/cuu/Value?mh
e = 1.602176634e-19 # C # https://physics.nist.gov/cgi-bin/cuu/Value?e
µ₀ = 1.25663706212e-6 # N A-2 # https://physics.nist.gov/cgi-bin/cuu/Value?mu0
h = 6.62607015e-34 # J Hz-1 # https://physics.nist.gov/cgi-bin/cuu/Value?h
eV = 1.602176634e-19 # J # https://physics.nist.gov/cgi-bin/cuu/Value?evj
ge = 2.00231930436256 # https://physics.nist.gov/cgi-bin/cuu/Value?gem
gµ = 2.0023318418 # https://physics.nist.gov/cgi-bin/cuu/Value?gmum
gp = 5.5856946893 # https://physics.nist.gov/cgi-bin/cuu/Value?gp
gd = 0.8574382338 # https://physics.nist.gov/cgi-bin/cuu/Value?gdn
gt = 5.957924931 # https://physics.nist.gov/cgi-bin/cuu/Value?gtn
gh = -4.255250615 # https://physics.nist.gov/cgi-bin/cuu/Value?ghn
Ps = CoulombTwoBody(z₁=-1, z₂=+1, m₁=me, m₂=me, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
Mu = CoulombTwoBody(z₁=-1, z₂=+1, m₁=me, m₂=mµ, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
H = CoulombTwoBody(z₁=-1, z₂=+1, m₁=me, m₂=mp, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
D = CoulombTwoBody(z₁=-1, z₂=+1, m₁=me, m₂=md, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
T = CoulombTwoBody(z₁=-1, z₂=+1, m₁=me, m₂=mt, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
he = CoulombTwoBody(z₁=-1, z₂=+2, m₁=me, m₂=mh, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
pμ = CoulombTwoBody(z₁=-1, z₂=+1, m₁=mµ, m₂=mp, mₑ=me, a₀=a₀, Eₕ=Eₕ, ℏ=ℏ)
ΔE_H = 2/3 / 4 * µ₀ * ℏ^2 * e^2 * ge / me * gp / mp * abs(ψ(H,0,0,0))^2
ΔE_D = 1 / 4 * µ₀ * ℏ^2 * e^2 * ge / me * gd / mp * abs(ψ(D,0,0,0))^2
ΔE_T = 2/3 / 4 * µ₀ * ℏ^2 * e^2 * ge / me * gt / mp * abs(ψ(H,0,0,0))^2
ΔE_Ps = 7/6 / 4 * µ₀ * ℏ^2 * e^2 * ge / me * ge / me * abs(ψ(Ps,0,0,0))^2
ΔE_Mu = 2/3 / 4 * µ₀ * ℏ^2 * e^2 * ge / me * gµ / mµ * abs(ψ(Mu,0,0,0))^2
ΔE_he = 2/3 / 4 * µ₀ * ℏ^2 * e^2 * ge / me * gh / mp * abs(ψ(he,0,0,0))^2
ΔE_pµ = 2/3 / 4 * µ₀ * ℏ^2 * e^2 * gµ / mµ * gp / mp * abs(ψ(pµ,0,0,0))^2
println("H \t", ΔE_H / h / 1e6, " MHz (Antique.jl + CODATA2018)")
println(" \t", "1420.405751768(1) MHz (https://doi.org/10.48550/arXiv.hep-ph/0109128)")
println("D \t", ΔE_D / h / 1e6, " MHz (Antique.jl + CODATA2018)")
println(" \t", "327.384352522(2) MHz (https://doi.org/10.48550/arXiv.hep-ph/0109128)")
println("T \t", ΔE_T / h / 1e6, " MHz (Antique.jl + CODATA2018)")
println(" \t", "1516.701470773(8) MHz (https://doi.org/10.48550/arXiv.hep-ph/0109128)")
println("Ps\t", ΔE_Ps / h / 1e6, " MHz (Antique.jl + CODATA2018)")
println(" \t", "203391.7(6) MHz (https://doi.org/10.48550/arXiv.hep-ph/0310099)")
println("Mu\t", ΔE_Mu / h / 1e6, " MHz (Antique.jl + CODATA2018)")
println(" \t", "4463.30278(5) MHz (https://doi.org/10.48550/arXiv.hep-ph/0109128)")
println("³He⁺\t", ΔE_he / h / 1e6, " MHz (Antique.jl + CODATA2018)")
println(" \t", "-8665.649867(10) MHz (https://doi.org/10.48550/arXiv.hep-ph/0109128)")
println("µp\t", ΔE_pµ / h / 1e12, " THz (Antique.jl + CODATA2018)")
println(" \t", 0.182725*eV / h / 1e12 ," THz (https://doi.org/10.1119/1.12733, https://doi.org/10.1016/j.nimb.2012.04.001)")
```

## Testing

Unit testing and Integration testing were done using computer algebra system ([Symbolics.jl](https://symbolics.juliasymbolics.org/stable/)) and numerical integration ([QuadGK.jl](https://juliamath.github.io/QuadGK.jl/stable/)). The test script is [here](https://github.com/ohno/Antique.jl/blob/main/test/CoulombTwoBody.jl).

```@eval
using Markdown
using Antique
Markdown.parse(Antique.load("../../test/result/CoulombTwoBody.log"))
```
2 changes: 1 addition & 1 deletion docs/src/HydrogenAtom.md
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ Potential energy curve:
using CairoMakie
f = Figure()
ax = Axis(f[1,1], xlabel=L"$r~/~a_0$", ylabel=L"$V(r)~/~E_\mathrm{h}0$", limits=(0.0,15.0,-2.0,0.2))
ax = Axis(f[1,1], xlabel=L"$r~/~a_0$", ylabel=L"$V(r)~/~E_\mathrm{h}$", limits=(0.0,15.0,-2.0,0.2))
lines!(ax, 0.1:0.01:20, r -> V(H, r))
f
```
Expand Down
28 changes: 9 additions & 19 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,6 @@ E(H)
Helium cation has symbol $\mathrm{He}^+$ and atomic number 2 ($Z=2$). Therefore the ground state ($n=1$) energy is $-2 E_\mathrm{h}$.

```julia
using Antique
He⁺ = HydrogenAtom(Z=2)
E(He⁺)
# output> -2.0
Expand All @@ -46,12 +45,6 @@ There are more examples on each model page.

```@raw html
<div class="catalog">
<div class="item">
<a target="_blank" href="./DeltaPotential">
<img src="assets/fig/DeltaPotential.png" alt="DeltaPotential"/>
</a>
<code>DeltaPotential</code>
</div>
<div class="item">
<a target="_blank" href="./InfinitePotentialWell">
<img src="assets/fig/InfinitePotentialWell.png" alt="InfinitePotentialWell"/>
Expand All @@ -76,21 +69,18 @@ There are more examples on each model page.
</a>
<code>MorsePotential</code>
</div>
<div class="item">
<a target="_blank" href="./RigidRotor">
<img src="assets/fig/RigidRotor.png" alt="RigidRotor"/>
</a>
<code>RigidRotor</code>
</div>
<div class="item">
<a target="_blank" href="./HydrogenAtom">
<img src="assets/fig/HydrogenAtom.png" alt="HydrogenAtom"/>
</a>
<code>HydrogenAtom</code>
</div>
</div>
```

- [Delta Potential](https://ohno.github.io/Antique.jl/stable/DeltaPotential/) `DeltaPotential`
- [Infinite Potential Well](https://ohno.github.io/Antique.jl/stable/InfinitePotentialWell/) `InfinitePotentialWell`
- [Harmonic Oscillator](https://ohno.github.io/Antique.jl/stable/HarmonicOscillator/) `HarmonicOscillator`
- [PoschlTeller](https://ohno.github.io/Antique.jl/stable/PoschlTeller/) `PoschlTeller`
- [Morse Potential](https://ohno.github.io/Antique.jl/stable/MorsePotential/) `MorsePotential`
- [Rigid Rotor](https://ohno.github.io/Antique.jl/stable/RigidRotor/) `RigidRotor`
- [Hydrogen Atom](https://ohno.github.io/Antique.jl/stable/HydrogenAtom/) `HydrogenAtom`
- [Coulomb 2-Body System](https://ohno.github.io/Antique.jl/stable/HydrogenAtom/) `CoulombTwoBody`

## Future Works

[List of quantum-mechanical systems with analytical solutions](https://en.wikipedia.org/wiki/List_of_quantum-mechanical_systems_with_analytical_solutions)
Expand Down
3 changes: 2 additions & 1 deletion src/Antique.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@ module Antique
:PoschlTeller,
:HarmonicOscillator3D,
:RigidRotor,
:InfinitePotentialWell3D
:InfinitePotentialWell3D,
:CoulombTwoBody,
]

# for Julia 1.1
Expand Down
Loading

0 comments on commit 334e930

Please sign in to comment.