diff --git a/docs/src/GEM.md b/docs/src/GEM.md index 6b9580d..2a5a0c6 100644 --- a/docs/src/GEM.md +++ b/docs/src/GEM.md @@ -143,6 +143,154 @@ round(eta_c.E[1] * 1000; digits=6) The 1.18 MeV difference from the paper table is consistent with using the rounded parameters and rounded variational exponent shown above. +## Examples of Hadron Spectroscopy + +The following three examples use the parameters and basis ranges specified in +[Issue #39](https://github.com/JuliaFewBody/TwoBody.jl/issues/39#issuecomment-5349154437). +Natural units (``\hbar=c=1``) are used: energies and masses are in GeV, +and Gaussian ranges are in GeV``^{-1}``. The displayed ground-state masses +are converted to MeV. All three examples use ``l=0`` Gaussian basis functions. + +### ``\Lambda_c(1/2^+)`` + +Parameters follow [Kim, Hiyama, Oka, and Suzuki (2020)](https://doi.org/10.1103/PhysRevD.102.014004). + +The charm quark and scalar diquark are treated as a two-body system, with +``\mu=M_{qq}M_c/(M_{qq}+M_c)`` and + +```math +\hat H = \frac{\boldsymbol p^2}{2\mu} + M_{qq} + M_c +- \frac{\alpha}{r} + \lambda r + C. +``` + +```@example lambda_c +using TwoBody + +Mqq = 0.725 +Mc = 1.750 +μ = inv(inv(Mqq) + inv(Mc)) +α = 0.06 / μ + +H = Hamiltonian( + RestEnergy(m=Mqq), + RestEnergy(m=Mc), + Kinetic(hbar=1, m=μ), + Coulomb(coefficient=-α), + Linear(coefficient=0.165), + Constant(constant=-0.83116597), +) +BS = GeometricBasisSet(GaussianBasis, 0.01, 9.0, 40) + +round(1000 * solve(H, BS; info=0).E[1]; digits=3) +``` + +### ``\eta_c(1S)``: Meng, Wang, and Oka + +Parameters follow [Meng, Wang, and Oka (2024)](https://doi.org/10.48550/arXiv.2404.01238). + +The AL1 Hamiltonian is + +```math +\hat H = \frac{\boldsymbol p^2}{2\mu} + m_1 + m_2 +- \frac{\kappa}{r} + \lambda r - \Lambda ++ \frac{2\pi\kappa'}{3m_1m_2} + \frac{e^{-r^2/r_0^2}}{\pi^{3/2}r_0^3} + \boldsymbol\sigma_1\cdot\boldsymbol\sigma_2. +``` + +Here ``r_0=A(2\mu)^{-B}`` and +``\langle\boldsymbol\sigma_1\cdot\boldsymbol\sigma_2\rangle=-3`` +for the spin-singlet state. + +```@example meng_charmonium +using TwoBody + +m₁ = 1.836 +m₂ = 1.836 +κ = 0.5069 +κ′ = 1.8609 +spin = -3 +μ = inv(inv(m₁) + inv(m₂)) +r₀ = 1.6553 * (2m₁ * m₂ / (m₁ + m₂))^(-0.2204) + +H = Hamiltonian( + RestEnergy(m=m₁), + RestEnergy(m=m₂), + Kinetic(hbar=1, m=μ), + Coulomb(coefficient=-κ), + Linear(coefficient=0.1653), + Constant(constant=-0.8321), + Gaussian( + coefficient=2π * κ′ * spin / (3m₁ * m₂ * (sqrt(π) * r₀)^3), + exponent=inv(r₀^2), + ), +) +BS = GeometricBasisSet(GaussianBasis, 0.1, 80.0, 20) + +round(1000 * solve(H, BS; info=0).E[1]; digits=3) +``` + +### ``\eta_c(1S)``: Arifi, Happ, Ohno, and Oka + +Parameters follow [Arifi, Happ, Ohno, and Oka (2024)](https://arxiv.org/abs/2401.07933). + +The semirelativistic Hamiltonian is + +```math +\hat H = \sqrt{m_1^2+\boldsymbol p^2}+\sqrt{m_2^2+\boldsymbol p^2} ++ a + br - \frac{4\alpha_s}{3r} ++ \frac{32\pi\alpha_s}{9m_1m_2} + \left(\frac{\lambda}{\sqrt{\pi}}\right)^3 e^{-\lambda^2r^2} + \boldsymbol S_1\cdot\boldsymbol S_2. +``` + +Here ``\lambda=\Lambda\sqrt{\mu}`` and +``\langle\boldsymbol S_1\cdot\boldsymbol S_2\rangle=-3/4``. +`RelativisticKinetic` represents ``\sqrt{m^2+\boldsymbol p^2}-m``, +so each constituent also needs a `RestEnergy` term. + +```@example arifi_charmonium +using TwoBody + +m₁ = 1.515 +m₂ = 1.515 +αs = 0.285 +spin = -3/4 +μ = inv(inv(m₁) + inv(m₂)) +λ = 1.437 * sqrt(μ) + +H = Hamiltonian( + RestEnergy(m=m₁), + RestEnergy(m=m₂), + RelativisticKinetic(m=m₁), + RelativisticKinetic(m=m₂), + Constant(constant=-0.189), + Linear(coefficient=0.092), + Coulomb(coefficient=-4αs/3), + Gaussian( + coefficient=32π * αs * (λ / sqrt(π))^3 * spin / (9m₁ * m₂), + exponent=λ^2, + ), +) +BS = GeometricBasisSet(GaussianBasis, 0.358, 2.720, 10) + +round(1000 * solve(H, BS; info=0).E[1]; digits=3) +``` + +### Comparison with the reference masses + +| Example | TwoBody.jl (MeV) | Reference quoted in Issue #39 (MeV) | +|:--|--:|--:| +| Kim, Hiyama, Oka, and Suzuki | 2286.000 | 2286 | +| Meng, Wang, and Oka | 3005.252 | 3005 | +| Arifi, Happ, Ohno, and Oka | 3020.148 | 3019 | + +The last example differs from the quoted reference by about 1.15 MeV with +the supplied parameters and basis ranges; it does not reproduce 3019 MeV +to the displayed precision. The regression tests allow 1 MeV for the Meng +example and 2 MeV for the Arifi example. These tolerances compare with the +quoted reference masses and are not estimates of numerical convergence. + ## API reference ```@docs; canonical=false @@ -152,3 +300,4 @@ TwoBody.ComplexGaussianBasisSet TwoBody.φp TwoBody.ψp ``` + diff --git a/test/Rayleigh-Ritz.jl b/test/Rayleigh-Ritz.jl index 6839826..c59f08e 100644 --- a/test/Rayleigh-Ritz.jl +++ b/test/Rayleigh-Ritz.jl @@ -57,6 +57,66 @@ @printf("%3d\t%.9f\t%.9f\t%s\n", 1, numerical, analytical, acceptance ? "✔" : "✗") @test acceptance + @testset "hadron spectroscopy" begin + Mqq = 0.725 + Mc = 1.750 + μ = inv(inv(Mqq) + inv(Mc)) + lambda_c = Hamiltonian( + RestEnergy(m=Mqq), + RestEnergy(m=Mc), + Kinetic(hbar=1, m=μ), + Coulomb(coefficient=-0.06 / μ), + Linear(coefficient=0.165), + Constant(constant=-0.83116597), + ) + lambda_c_basis = GeometricBasisSet(GaussianBasis, 0.01, 9.0, 40) + @test solve(lambda_c, lambda_c_basis; info=0).E[1] ≈ 2.286 atol=1e-6 + + m₁ = 1.836 + m₂ = 1.836 + κ = 0.5069 + κ′ = 1.8609 + spin = -3 + μ = inv(inv(m₁) + inv(m₂)) + r₀ = 1.6553 * (2m₁ * m₂ / (m₁ + m₂))^(-0.2204) + meng = Hamiltonian( + RestEnergy(m=m₁), + RestEnergy(m=m₂), + Kinetic(hbar=1, m=μ), + Coulomb(coefficient=-κ), + Linear(coefficient=0.1653), + Constant(constant=-0.8321), + Gaussian( + coefficient=2π * κ′ * spin / (3m₁ * m₂ * (sqrt(π) * r₀)^3), + exponent=inv(r₀^2), + ), + ) + meng_basis = GeometricBasisSet(GaussianBasis, 0.1, 80.0, 20) + @test solve(meng, meng_basis; info=0).E[1] ≈ 3.005 atol=1e-3 + + m₁ = 1.515 + m₂ = 1.515 + αs = 0.285 + spin = -3/4 + μ = inv(inv(m₁) + inv(m₂)) + λ = 1.437 * sqrt(μ) + arifi = Hamiltonian( + RestEnergy(m=m₁), + RestEnergy(m=m₂), + RelativisticKinetic(m=m₁), + RelativisticKinetic(m=m₂), + Constant(constant=-0.189), + Linear(coefficient=0.092), + Coulomb(coefficient=-4αs/3), + Gaussian( + coefficient=32π * αs * (λ / sqrt(π))^3 * spin / (9m₁ * m₂), + exponent=λ^2, + ), + ) + arifi_basis = GeometricBasisSet(GaussianBasis, 0.358, 2.720, 10) + @test solve(arifi, arifi_basis; info=0).E[1] ≈ 3.019 atol=2e-3 + end + @testset "Pérez-Torres STO-3G contraction" begin exponents = (0.109818, 0.405771, 2.227660) contraction = (0.444635, 0.535328, 0.154329)