Transitions and decays

GIModel.QuarkModelTransitions computes decay amplitudes and widths from the states and wavefunctions of a GIModel spectrum. It is a submodule of GIModel, so there is nothing extra to install:

using GIModel
using GIModel.QuarkModelTransitions
params, mq = load_parameters_and_quark_masses(default_parameters_path())
spec = compute_spectrum(params, Meson(mq, :c, :c);
    levels = spectrum_levels(2; L_labels = ("S", "P")), solver = OscillatorSolver())

The pattern

Every calculation follows the same four steps:

1. Resolve the states: physical_state(spectrum, label).

2. Choose an operator: photon emission, pseudoscalar emission, a leptonic current, or an annihilation.

3. Evaluate the amplitude: matrix_element(final, operator, initial).

4. Convert to a width in MeV: decay_width.

psi = physical_state(spec, "1^3S_1")
eta_c = physical_state(spec, "1^1S_0")
amplitude = matrix_element(eta_c, PhotonEmission(mq), psi)

RadiativeAmplitude(1^3S_1 -> 1^1S_0 + gamma, DirectM1, recoil_order=0, 0.6881204742164155 + 0.0im)

decay_width(amplitude)       # MeV

0.0023635549674362938

The argument order is always (final, operator, initial), like   .

Physical states

A PhysicalState is a self-contained copy of one meson: its label, mass, , and the signed components with their radial waves.

psi

PhysicalState("1^3S_1", 3.0913730064733147 GeV, 1 component)

Because it holds all components, a mixed parent or daughter contributes coherently, with interference, and nothing needs to be re-assembled by hand. Masses come from the spectrum, so decay momenta follow from the model masses. Comparing at other masses shows how to use different masses.

What can be computed

operator	process	final state for matrix_element	final state for decay_width
PhotonEmission	(E1, M1, M2)	daughter PhysicalState	the same
MultipolePhotonEmission	, every multipole, any	daughter PhysicalState	the same
PseudoscalarEmission		TwoMesonChannel(B, P)	the same
LeptonicCurrent	,	Vacuum()	MasslessLeptonPair(), LeptonNeutrinoChannel
TwoPhotonAnnihilation		TwoPhotonChannel()	the same
GluonicAnnihilation	,	not available	TwoGluonChannel(), ThreeGluonChannel()

Widths are always in MeV. The amplitude units depend on the process:

result	.value is
M1 photon amplitude	magnetic moment in nuclear magnetons
E1 and M2 photon amplitude	width amplitude in
multipole photon amplitude	.multipoles: one width amplitude in per multipole
pseudoscalar emission	dimensionless partial-wave amplitudes (see below)
leptonic current	dimensionless reduced current
two-photon	width amplitude in

Photon emission

PhotonEmission follows Appendix D of the paper. You do not name the multipole: it is inferred from the two states. With verbose = true the selection is logged:

chi_c1 = physical_state(spec, "1^3P_1")
e1 = matrix_element(psi, PhotonEmission(mq), chi_c1; verbose = true)
decay_width(e1)

┌ Info: photon transition selected
│   initial = "1^3P_1"
│   final = "1^3S_1"
│   transition_class = :AllowedE1
│   multipole = :E1
│   recoil_order = 0
└   q_GeV = 0.39208783378512374

0.21042114059605513

The classes are DirectM1 (same radial level), HinderedM1 (different radial levels), AllowedE1 (spin-conserving S–P), SpinFlipE1 and SpinFlipM2. PhotonEmission(mq; recoil_order = 2) adds the correction where it is implemented (currently M1).

The photon couples to the quark charges, so flavor must be explicit. States built from :u and :d work; the averaged :q is rejected for photon emission because it cannot distinguish   from  . States obtained from an isoscalar spectrum carry this information and are resolved into coherent u and d contributions automatically.

The spin-flip kernels ( ) use the denominators   (J=2, M2) and   (J=1, E1) and weight each emitter by its centre-of-mass position, as E1 does. The paper prints   and , which overstate these rates by a factor of 2; see the Table VI report.

Nonzero D-wave admixtures are rejected by these Appendix-D kernels with: “Not implemented yet. Please submit issue if needed, and/or PR with implementation.” No components are silently discarded. Use MultipolePhotonEmission for D waves and mixed states.

All multipoles

MultipolePhotonEmission evaluates the quark convection and spin currents of the same interaction, exact in the photon momentum, for any pair of states, and resolves the result into multipoles. Electric multipoles are taken from the charge density (Siegert’s theorem with  ), as for the leading E1. There is no smearing. In    the E1 is accompanied by an M2 from the charm magnetic moment, and multipole_fractions gives the normalized amplitudes that CLEO and BESIII measure:

chi_c2 = physical_state(spec, "1^3P_2")
multi = matrix_element(psi, MultipolePhotonEmission(mq), chi_c2)
multipole_fractions(multi)

Dict{Symbol, Float64} with 3 entries:
  :E1 => 0.995588
  :M2 => -0.0938277
  :E3 => 0.0

The M2 fraction is close to the first-order value  . An E3 needs S–D or P–F mixing, absent from these pure states. The width decay_width(multi) also includes the spin-current contribution to E1, a correction of relative order that the leading kernel omits.

Pseudoscalar emission

PseudoscalarEmission(g, h, masses) is the elementary-emission operator of Eq. (19),

where is the momentum of the emitted pseudoscalar and the momentum of the quark that emits it. The couplings g and h (in ) are not predicted by the spectrum; you supply them, typically fitted to one or two known widths.

The final channel is ordered: TwoMesonChannel(surviving, emitted), and the emitted meson must be a state. The result holds one amplitude per allowed partial wave :

# K*⁺ (u s̄) → K⁰ (d s̄) + π⁺ (u d̄), with unit direct coupling g and h = 0.
S1, S0 = BasisState(1, "S", 3, 1), BasisState(1, "S", 1, 0)
solver = OscillatorSolver()
us = compute_spectrum(params, Meson(mq, :u, :s); levels = [S1], solver)
ds = compute_spectrum(params, Meson(mq, :d, :s); levels = [S0], solver)
ud = compute_spectrum(params, Meson(mq, :u, :d); levels = [S0], solver)
K_star, K_zero, pion = physical_state(us, S1), physical_state(ds, S0), physical_state(ud, S0)
strong = matrix_element(TwoMesonChannel(K_zero, pion), PseudoscalarEmission(1.0, 0.0, mq), K_star)

TransitionAmplitude
  transition     1^3S_1 -> 1^1S_0 + 1^1S_0
  operator       PseudoscalarEmission
  momentum       0.31207606274646865 GeV
  normalization  RelativisticTwoBodyNormalization
  helicities     1
  partial waves  1
    (L,S)=(1,0)  -0.40378226230698816 + 0.0im
  terms          1

[w => strong[w] for w in partial_waves(strong)]

1-element Vector{Pair{PartialWave, ComplexF64}}:
 PartialWave(1, 0) => -0.40378226230698816 + 0.0im

A vector decaying to two pseudoscalars has one P-wave amplitude. The flavor algebra is exact: a channel that violates charge or flavor, such as  , gives exactly zero. See Strong decays beyond the paper for neutral isospin states, axial mesons with two partial waves, and a calibration of g and h.

This operator uses the full solved wavefunctions of parent and daughter, with all radial nodes and each state’s own size. The paper’s Table V instead used a single oscillator scale   GeV for every meson. The 1985 paper section explains how that table is reproduced separately.

Leptonic decays

A vector meson annihilates into a lepton pair through the electromagnetic current. LeptonicCurrent(:electromagnetic, masses) takes the charges from the flavor components:

current = LeptonicCurrent(:electromagnetic, mq)
decay_width(MasslessLeptonPair(), current, psi)     # Γ(J/ψ → e⁺e⁻) in MeV

0.010511053477744904

matrix_element(Vacuum(), current, psi) returns the dimensionless reduced current instead of a width. Weak decays of pseudoscalars use a LeptonNeutrinoChannel with an explicit CKM element.

Annihilation into photons and gluons

TwoPhotonAnnihilation takes explicit charge factors as AnnihilationTerms: for the effective squared charge is  .

two_photon = TwoPhotonAnnihilation(mq, AnnihilationTerm((:c, :c), 4 / 9))
decay_width(TwoPhotonChannel(), two_photon, eta_c)  # Γ(η_c → γγ) in MeV

0.0075250209715686505

GluonicAnnihilation takes the strong coupling, either as a number or as a function of the meson mass:

gluonic = GluonicAnnihilation(mq, M -> alpha_s_q(M))
(eta_c = decay_width(TwoGluonChannel(), gluonic, eta_c),
 psi = decay_width(ThreeGluonChannel(), gluonic, psi))

(eta_c = 24.583472509534875, psi = 0.1897757984435488)

These annihilation widths depend on the smeared wavefunction at the origin (Eq. (17)), which is sensitive to the short-distance part of the wave. Check their convergence when you quote them.

Comparing at other masses

matrix_element and decay_width use the model masses stored in the states. To see how an amplitude would change at an experimental mass or momentum, keeping the wavefunctions fixed, use mass_correction_factor:

factor = mass_correction_factor(TwoPhotonChannel(), two_photon, eta_c; target_mass = 2.984)

1.0085644209633375 + 0.0im

The factor multiplies the amplitude (.value), not the width, except for gluonic annihilation, where no amplitude exists. The required keyword depends on the operator:

operator	keyword
LeptonicCurrent, TwoPhotonAnnihilation, GluonicAnnihilation	target_mass (GeV)
PhotonEmission	target_momentum (GeV)
PseudoscalarEmission	target_momentum (GeV) and partial_wave

States you build yourself

A PhysicalState can also be constructed by hand from any radial waves, for example a simple oscillator function for an analytic check, or a linear combination of solved states for an isospin eigenstate:

wave = OscillatorWave(0, 0.5, [1.0])     # one oscillator function, β = 0.5 GeV
toy = PhysicalState("toy ψ", 3.10, [(
    basis = BasisState(1, "S", 3, 1; flavors = (:c, :c)), coefficient = 1.0, wave = wave)])
toy.mass_GeV, toy.J

(3.1, 1)

The package does not check such states: their masses and waves are your assumptions.

Charge radii

charge_radius_squared computes the mean-square charge radius of a meson from its radial wave (Table VII(d) of the paper), in :

pion_wave = radial_wave(ud, "1^1S_0")
r2 = charge_radius_squared(pion_wave, mq["u"], 2 / 3, mq["d"], 1 / 3)
r2 * 0.1973^2            # fm²

0.42333886261382325

What is not covered

The operators above are the Godfrey–Isgur ones. The module does not implement quark-pair-creation () decay models, emission of mesons other than pseudoscalars, weak hadronic or semileptonic transitions, or coupled-channel effects. See Scope and limitations.