Two Classic Math Methods Nail the Masses of Bottom-Quark Mesons to Within Half a Percent

Inside every bottom-flavored meson, a heavy quark and its partner are bound by the strong force in a regime that quantum chromodynamics, the theory of the strong interaction, cannot yet crack with pen and paper alone. A new theoretical study published in The European Physical Journal C shows that two of the oldest workhorses of quantum mechanics, the matrix Numerov method and the variational method, can be combined with a deceptively simple potential to reproduce the measured masses of the B, Bs, and Bc mesons with deviations of less than 0.5 percent from the experimental averages compiled by the Particle Data Group. The result, from Pritom Kumar Nath, Rishan Dev Choudhury, and Jugal Lahkar at Tezpur University in India, offers a fresh benchmark for how far inexpensive, non-relativistic calculations can be pushed before the full machinery of lattice QCD or QCD sum rules becomes indispensable.

The three mesons in the study span an instructive range of quark configurations. The B meson pairs a bottom antiquark with a light down quark, the Bs substitutes a heavier strange quark, and the Bc is a true heavyweight duel, binding a bottom antiquark to a charm quark. Because the bottom quark weighs nearly 5 GeV, roughly twenty-five times the intrinsic QCD energy scale of about 0.2 GeV, these systems allow physicists to separate the slow motion of the heavy quark from the frenetic dynamics of the light degrees of freedom. That separation is what makes effective descriptions, and potential models in particular, viable at all. The team modeled the quark-antiquark interaction with the Killingbeck potential, a compact expression containing a quadratic term, a linear term, and an inverse-distance Coulombic term, each carrying an unknown coefficient.

The choice of a harmonic quadratic term might look arbitrary, but the authors ground it in serious field theory. At finite temperature or in a complex non-perturbative vacuum, the famous Cornell potential, with its short-range Coulomb piece from one-gluon exchange and its long-range linear confinement from the QCD string, becomes screened by Debye shielding of the gluonic fields. Expanding the screened Karsch-Mehr-Satz potential in powers of the quark separation and truncating at second order produces exactly the Killingbeck form, with the quadratic coefficient emerging from the interplay between screening corrections and the thermal attenuation of confinement. In other words, the harmonic term is not a mathematical convenience but a footprint of how the gluonic vacuum responds to a bound heavy quark.

To fix the unknown Killingbeck parameters, the researchers turned to ordinary least squares regression, the same statistical tool used across the sciences for curve fitting. They evaluated the Cornell potential, whose parameters are well established in the literature, on a grid of radial points spanning 0 to 2.0 inverse GeV, a range that comfortably encloses the charge radii of all three mesons. Setting up the problem as a multivariate linear regression, they minimized the sum of squared residuals between the two potentials and solved the resulting normal equations analytically. The calibrated Killingbeck potential then served as the input interaction for every subsequent calculation, ensuring that both numerical methods operated on identical physics.

The first computational engine was the matrix Numerov method, a discretization scheme celebrated for its fourth-order global accuracy, meaning errors shrink as the fourth power of the grid spacing. The method recasts the radial Schrödinger equation into a matrix eigenvalue problem: the kinetic energy operator becomes a finite-difference matrix, the potential is averaged over neighboring grid points with the characteristic Numerov weights, and Dirichlet boundary conditions enforce the bound-state requirement that the wavefunction vanishes at the edges of the domain. Using 500 grid points over a range of 30 inverse GeV, the team extracted ground-state energies and wavefunctions for each meson, from which the probability density at the origin follows directly.

The second engine was the variational method, armed with a Gaussian trial wavefunction whose width is a single adjustable parameter. Minimizing the expectation value of the Hamiltonian with respect to that parameter yields an upper bound on the true ground-state energy. The Gaussian choice is not merely convenient here; because the Killingbeck potential contains a harmonic term, a Gaussian is in fact the exact ground-state solution of a pure harmonic oscillator, making it an exceptionally well-matched ansatz. The agreement between the two independent methods was striking: the mass spectra showed parametric stability with standard deviations of just 0.028 GeV for the Numerov approach and 0.031 GeV for the variational one, and both landed within half a percent of the measured B, Bs, and Bc masses of 5.279, 5.366, and 6.274 GeV respectively.

Beyond masses, the calculations produced decay constants, the quantities that quantify the overlap of the quark and antiquark wavefunctions at the origin and govern how readily a charged meson annihilates into a lepton-neutrino pair through a virtual W boson. Here the picture grew more nuanced. For the B meson, both methods systematically overestimated the leptonic decay constant, with the Numerov value of 0.365 GeV exceeding experiment by roughly 84 percent and the variational value of 0.314 GeV by about 58 percent. For the Bs meson, agreement improved markedly, with deviations shrinking to around 10 to 15 percent. For the Bc, where both constituents are heavy, the calculated upper bounds aligned well with existing theoretical predictions. The pattern is physically telling: the lighter the companion quark, the more relativistic its motion, and the less appropriate a purely non-relativistic framework becomes.

The study also computed the oscillation frequencies of neutral B and Bs mesons, which spontaneously transform into their antiparticles through weak-interaction mixing, a phenomenon that has become one of the most sensitive probes of new physics at LHCb. Because the mass difference between the mixed states depends quadratically on the decay constant, the overestimated B-meson decay constant propagates into the oscillation prediction, while the Bs results bracket the experimental value more tightly. Additionally, the team extracted the slope and curvature of the Isgur-Wise function, the universal form factor describing how the light-quark cloud survives a sudden velocity change of the heavy quark during semileptonic decay. The B-meson slope of roughly 0.338 satisfies the model-independent Bjorken lower bound of one quarter but falls short of the stronger Uraltsev bound of three quarters, a shortfall the authors attribute squarely to missing relativistic corrections. The Bc slope of about 3.70 and curvature near 4.0 are dramatically larger, reflecting the ultra-compact, tightly bound geometry of the double-heavy system, where the uncertainty principle dictates a broad momentum-space wavefunction exquisitely sensitive to recoil.

The authors are candid about the limitations. The overestimation of the wavefunction at the origin in light-quark systems stems from the rigid, localized behavior of the potential at zero separation, and they argue that incorporating the Darwin term, which smears the singular interaction over the quantum fluctuations of the light quark, together with a mock-meson mass formulation, would be the natural remedy. Relativistic corrections, they conclude, are the most probable source of the remaining discrepancies in both decay constants and Isgur-Wise parameters, rather than any flaw in the confining potential’s functional form.

Perhaps the most forward-looking implication lies in hot matter. Because the screened Coulomb potential in a quark-gluon plasma reduces precisely to the Killingbeck form with a temperature-dependent quadratic coefficient, the same framework could track how meson bound states dissolve as the medium heats past the deconfinement transition, yielding critical dissociation temperatures for comparison with heavy-ion collision data. What began as a comparative exercise in numerical methods thus doubles as a template for probing quark confinement under the most extreme conditions the universe has to offer, from the first microseconds after the Big Bang to the fireballs recreated at colliders today.

Subject of Research: Non-relativistic potential model calculations of the spectroscopic and structural properties of bottom-flavored B, Bs, and Bc mesons

Article Title: Spectroscopic and structural properties of B, (B_s), and (B_c) mesons within a non-relativistic potential model: a comparative analysis via matrix Numerov and variational methods

Article References: Kumar Nath, P., Dev Choudhury, R., & Lahkar, J. (2026). Spectroscopic and structural properties of B, $$B_s$$, and $$B_c$$ mesons within a non-relativistic potential model: a comparative analysis via matrix Numerov and variational methods. The European Physical Journal C, 86(10), Article 1159. https://doi.org/10.1140/epjc/s10052-026-16433-1

Image Credits: AI Generated

DOI: 10.1140/epjc/s10052-026-16433-1

Keywords: B mesons, Bs meson, Bc meson, Killingbeck potential, Cornell potential, matrix Numerov method, variational method, decay constants, Isgur-Wise function, quarkonium spectroscopy, particle-antiparticle mixing, quark-gluon plasma