Package tparton

tParton: Evolution of transversity parton distribution functions

tParton is a Python package for evolving transversity PDFs using two complementary methods:

1. Direct integration (Hirai method): Numerically solves the DGLAP equations using discretized grids in x and Q².

2. Mellin moment method (Vogelsang method): Uses Mellin transforms and inverse transform via Cohen contour integration for faster, more accurate evolution.

Basic Usage

Command-line interface::

# Using Mellin moment method
python -m tparton m input.dat 3.1 10.6 --morp plus -o output.dat

# Using direct integration method  
python -m tparton t input.dat 3.1 10.6 --morp plus -o output.dat

Python API::

from tparton.m_evolution import evolve as m_evolve
from tparton.t_evolution import evolve as t_evolve

# Evolve using Mellin method (faster)
result = m_evolve(input_pdf, Q0_squared=3.1, Q_squared=10.6, 
                  morp='plus', order='NLO')

# Evolve using direct integration (more control over discretization)
result = t_evolve(input_pdf, Q0_squared=3.1, Q_squared=10.6,
                  morp='plus', order='NLO')

Navigation

• Home: https://mikesha2.github.io/tParton/
• Examples & Tutorials: https://mikesha2.github.io/tParton/examples.html
• API Documentation: https://mikesha2.github.io/tParton/api/tparton/

See Also

• GitHub repository: https://github.com/mikesha2/tParton
• Paper (arXiv): https://arxiv.org/abs/2409.00221
• Paper (JOSS): https://mikesha2.github.io/tParton/paper.pdf

Sub-modules

tparton.constants

QCD constants and parameters used in PDF evolution …

tparton.m_evolution

Mellin moment method for transversity PDF evolution (Vogelsang method) …

tparton.t_evolution

Direct integration method for transversity PDF evolution (Hirai method) …

Functions

def m_evolve(pdf: numpy.ndarray,
Q0_2: float = 0.16,
Q2: float = 5.0,
l_QCD: float = 0.25,
n_f: int = 5,
CG: float = 3,
morp: str = 'minus',
order: int = 2,
n_x: int = 200,
verbose: bool = False,
Q0_2_a: float = 8315.178393760001,
a0: float = mpf('0.0093901416424218252'),
alpha_num: bool = True,
degree: int = 5) ‑> numpy.ndarray

Expand source code

def evolve(
    pdf: np.ndarray,
    Q0_2: float = 0.16,
    Q2: float = 5.0,
    l_QCD: float = 0.25,
    n_f: int = 5,
    CG: float = 3,
    morp: str = 'minus',
    order: int = 2,
    n_x: int = 200,
    verbose: bool = False,
    Q0_2_a: float = 91.1876**2,
    a0: float = 0.118 / 4 / pi,
    alpha_num: bool = True,
    degree: int = 5,
) -> np.ndarray:
    """Evolve transversity PDF using the Mellin moment method (Vogelsang).
    
    This is the main function for PDF evolution using Mellin transforms.
    Faster and less discretization-dependent than the direct integration method.
    
    This method:
    1. Computes Mellin moments (Eq. 20): M[f](s) = ∫₀^∞ x^{s-1} f(x) dx
    2. Evolves moments using Eq. (24) with splitting function moments from
       Eq. (26) for LO and Eq. (27) for NLO
    3. Reconstructs PDF via inverse Mellin transform (Eq. 36) using Cohen's
       method with accelerated alternating series convergence
    
    Parameters
    ----------
    pdf : ndarray
        Input PDF as x*f(x). Can be 1D array (values at x evenly
        spaced on [0, 1]) or 2D array ([[x0, x0*f(x0)], [x1, x1*f(x1)], ...]).
    Q0_2 : float, optional
        Initial energy scale squared in GeV² (default: 0.16).
    Q2 : float, optional
        Final energy scale squared in GeV² (default: 5.0).
    l_QCD : float, optional
        QCD scale parameter Λ in GeV (default: 0.25).
        Only used if alpha_num=False.
    n_f : int, optional
        Number of active quark flavors (default: 5).
    CG : float, optional
        Number of colors, NC (default: 3).
    morp : str, optional
        Distribution type (default: 'minus'). Options are 'plus'
        (ΔT q⁺ = ΔT u + ΔT d) or 'minus' (ΔT q⁻ = ΔT u - ΔT d).
    order : int, optional
        Perturbative order (default: 2). Use 1 for LO or 2 for NLO.
    n_x : int, optional
        Number of x grid points minus 1 for output (default: 200).
    verbose : bool, optional
        Print (x, x*pdf(x)) during evolution if True (default: False).
    Q0_2_a : float, optional
        Reference scale Q₀² where αs is known, in GeV² (default: 91.1876²).
        Only used if alpha_num=True.
    a0 : float, optional
        Reference coupling αs(Q0_2_a)/(4π) (default: 0.118/(4π)).
        Only used if alpha_num=True.
    alpha_num : bool, optional
        Use numerical ODE evolution for αs if True (default: True).
        If False, uses analytical approximation.
    degree : int, optional
        Convergence acceleration degree for inverse Mellin (default: 5).
        Higher values increase accuracy but slow computation.
    
    Returns
    -------
    ndarray
        Evolved PDF as 2D array [[x, x*f_evolved(x)], ...]. Shape: (n_x+2, 2)
        due to padding at boundaries.
    
    Notes
    -----
    Advantages over direct integration (t_evolution.evolve):
    - Typically 10-100x faster
    - Less sensitive to discretization
    - Better for smooth PDFs
    
    Disadvantages:
    - Less direct control over integration
    - May have issues with very peaked PDFs
    
    Examples
    --------
    >>> import numpy as np
    >>> from tparton.m_evolution import evolve
    >>> x = np.linspace(0, 1, 100)
    >>> pdf_in = x * (1-x)**3  # x*f(x) format
    >>> pdf_out = evolve(pdf_in, Q0_2=4.0, Q2=100.0, order=2)
    >>> x_out, xf_out = pdf_out[:, 0], pdf_out[:, 1]
    
    See Also
    --------
    t_evolution.evolve : Direct integration method (Hirai)
    evolveMoment : Evolves a single Mellin moment
    inv_mellin : Inverse Mellin transform (Cohen method)
    
    References
    ----------
    - Vogelsang, W. (1998). Phys. Rev. D 57, 1886-1894
    - Sha, C.M. & Ma, B. (2025). arXiv:2409.00221
    """

    if pdf.shape[-1] == 1:
        # If only the x*pdf(x) values are supplied, assume a linear spacing from 0 to 1
        xs = np.linspace(0, 1, len(pdf))
    else:
        # Otherwise split the input array
        xs, pdf = pdf[:, 0], pdf[:, 1]
    
    # Divide x*pdf(x) by x. 
    # In the Hirai method, the evolution of x*pdf(x) and pdf(x) are numerically identical and do not require this extra step.
    pdf = pdf / (xs + 1e-100)
    # We assume that pdf(0) = 0
    pdf[0] = 0

    # Interpolate the resulting (x, pdf(x)) pairs as a function
    pdf_fun = interp(xs, pdf, fill_value=0, assume_sorted=True)
    
    # Convert the pdf into one compatible with mpmath's internal floating point representation
    pdf = lambda x: mp.mpf(pdf_fun(float(x)).item())

    # The type of distribution determines eta in Eq. (31)
    eta = 1 if morp == 'plus' else -1

    # Calculate the color constants
    NC, CF, Tf, beta0, beta1 = constants(CG, n_f)

    if order == 1:
        # If the desired order of accuracy is LO, we simply set beta1 to 0, which reproduces the relevant LO equations
        beta1 = 0
        # For the sake of efficiency, we also redefine the NLO splitting function moment to be the zero function
        NLO_splitting_function_moment = lambda n, eta, CF, NC, Tf: 0
    
    if alpha_num:
        # Use the numerically evolved alpha_S
        alpha_S_Q0_2 = alpha_S_num(Q0_2, order, Q0_2_a, a0, beta0, beta1)
        alpha_S_Q2 = alpha_S_num(Q2, order, Q0_2_a, a0, beta0, beta1)
    else:
        # Use the approximate analytical expression for alpha_S in Eq. (4)
        alpha_S_Q0_2 = alpha_S(Q0_2, order, beta0, beta1, l_QCD)
        alpha_S_Q2 = alpha_S(Q2, order, beta0, beta1, l_QCD)

    # Choose the values of x at which the evolved pdf(x) will be evaluated
    if n_x > 0:
        xs = np.linspace(0, 1, n_x+2)
    # In all cases, we assume that xs[0] = 0 and xs[1] = 1, pdf(0) = pdf(1) = 0, so no evolution is necessary at these points.
    # Even if pdf(0) != 0, this slight change will not significantly affect the final numerical result.
    xs = xs[1:-1]

    # A function representing the Mellin transform of pdf(x), Eq. (20)
    pdf_m = lambda s: mellin(pdf, s)
    # A function representing the resulting evolved moments, Eq. (24)
    pdf_evolved_m = lambda s: mpc(evolveMoment(s, pdf_m(s), alpha_S_Q0_2, alpha_S_Q2, beta0, beta1, eta, CF, NC, Tf))
    # Perform Mellin inversion on the evolved moments, Eq. (36)
    pdf_evolved = np.array([inv_mellin(pdf_evolved_m, x, degree=degree, verbose=verbose).__complex__().real for x in xs])

    # Reinstate the endpoints x = 0 and x = 1
    xs = np.pad(xs, 1)
    xs[-1] = 1

    # Pad the evolved pdf so that pdf(0) = pdf(1) = 0
    pdf_evolved = np.pad(pdf_evolved, 1)
    # Organize the (x, x*pdf_evolved(x)) pairs into an array
    pdf_evolved = np.stack((xs, np.array(xs) * np.array(pdf_evolved)), axis=1)
    print('Done!')
    return pdf_evolved

Evolve transversity PDF using the Mellin moment method (Vogelsang).

This is the main function for PDF evolution using Mellin transforms. Faster and less discretization-dependent than the direct integration method.

This method: 1. Computes Mellin moments (Eq. 20): Mf = ∫₀^∞ x^{s-1} f(x) dx 2. Evolves moments using Eq. (24) with splitting function moments from Eq. (26) for LO and Eq. (27) for NLO 3. Reconstructs PDF via inverse Mellin transform (Eq. 36) using Cohen's method with accelerated alternating series convergence

Parameters

pdf : ndarray

Input PDF as xf(x). Can be 1D array (values at x evenly spaced on [0, 1]) or 2D array ([[x0, x0f(x0)], [x1, x1*f(x1)], …]).

Q0_2 : float, optional

Initial energy scale squared in GeV² (default: 0.16).

Q2 : float, optional

Final energy scale squared in GeV² (default: 5.0).

l_QCD : float, optional

QCD scale parameter Λ in GeV (default: 0.25). Only used if alpha_num=False.

n_f : int, optional

Number of active quark flavors (default: 5).

CG : float, optional

Number of colors, NC (default: 3).

morp : str, optional

Distribution type (default: 'minus'). Options are 'plus' (ΔT q⁺ = ΔT u + ΔT d) or 'minus' (ΔT q⁻ = ΔT u - ΔT d).

order : int, optional

Perturbative order (default: 2). Use 1 for LO or 2 for NLO.

n_x : int, optional

Number of x grid points minus 1 for output (default: 200).

verbose : bool, optional

Print (x, x*pdf(x)) during evolution if True (default: False).

Q0_2_a : float, optional

Reference scale Q₀² where αs is known, in GeV² (default: 91.1876²). Only used if alpha_num=True.

a0 : float, optional

Reference coupling αs(Q0_2_a)/(4π) (default: 0.118/(4π)). Only used if alpha_num=True.

alpha_num : bool, optional

Use numerical ODE evolution for αs if True (default: True). If False, uses analytical approximation.

degree : int, optional

Convergence acceleration degree for inverse Mellin (default: 5). Higher values increase accuracy but slow computation.

Returns

ndarray

Evolved PDF as 2D array [[x, x*f_evolved(x)], …]. Shape: (n_x+2, 2) due to padding at boundaries.

Notes

Advantages over direct integration (t_evolution.evolve): - Typically 10-100x faster - Less sensitive to discretization - Better for smooth PDFs

Disadvantages: - Less direct control over integration - May have issues with very peaked PDFs

Examples

>>> import numpy as np
>>> from tparton.m_evolution import evolve
>>> x = np.linspace(0, 1, 100)
>>> pdf_in = x * (1-x)**3  # x*f(x) format
>>> pdf_out = evolve(pdf_in, Q0_2=4.0, Q2=100.0, order=2)
>>> x_out, xf_out = pdf_out[:, 0], pdf_out[:, 1]

See Also

evolve()

Direct integration method (Hirai)

evolveMoment

Evolves a single Mellin moment

inv_mellin

Inverse Mellin transform (Cohen method)

References

• Vogelsang, W. (1998). Phys. Rev. D 57, 1886-1894
• Sha, C.M. & Ma, B. (2025). arXiv:2409.00221

def t_evolve(pdf: numpy.ndarray,
Q0_2: float = 0.16,
Q2: float = 5.0,
l_QCD: float = 0.25,
n_f: int = 5,
CG: float = 3,
n_t: int = 100,
n_z: int = 500,
morp: str = 'plus',
order: int = 2,
logScale: bool = False,
verbose: bool = False,
Q0_2_a: float = 8315.178393760001,
a0: float = 0.009390141642421825,
alpha_num: bool = True) ‑> numpy.ndarray

Expand source code

def evolve(
    pdf: np.ndarray,
    Q0_2: float = 0.16,
    Q2: float = 5.0,
    l_QCD: float = 0.25,
    n_f: int = 5,
    CG: float = 3,
    n_t: int = 100,
    n_z: int = 500,
    morp: str = 'plus',
    order: int = 2,
    logScale: bool = False,
    verbose: bool = False,
    Q0_2_a: float = 91.1876 ** 2,
    a0: float = 0.118 / 4 / np.pi,
    alpha_num: bool = True
) -> np.ndarray:
    """Evolve transversity PDF using the direct integration method (Hirai).
    
    This is the main function for PDF evolution using direct numerical integration
    of the DGLAP equation (Eq. 1). More robust for peaked PDFs but slower than
    Mellin method.
    
    This method:
    1. Discretizes t = ln(Q²) into n_t Euler steps
    2. At each step, computes the convolution integral (Eq. 19):
       
           f̃(x) ⊗ g(x) = ∫ dx̃ f̃(x/z) g(z)
       
       using Simpson's rule with n_z integration points
    3. Updates PDF using forward Euler: f̃(t + dt) ≈ f̃(t) + dt·f̃'(t)
    
    Parameters
    ----------
    pdf : ndarray
        Input PDF as x*f(x). Can be 1D array (values at x evenly
        spaced on [0, 1]) or 2D array ([[x0, x0*f(x0)], [x1, x1*f(x1)], ...]).
    Q0_2 : float, optional
        Initial energy scale squared in GeV² (default: 0.16).
    Q2 : float, optional
        Final energy scale squared in GeV² (default: 5.0).
    l_QCD : float, optional
        QCD scale parameter Λ in GeV (default: 0.25).
        Only used if alpha_num=False.
    n_f : int, optional
        Number of active quark flavors (default: 5).
    CG : float, optional
        Number of colors, NC (default: 3).
    n_t : int, optional
        Number of Euler time steps (default: 100).
        More steps = better accuracy but slower.
    n_z : int, optional
        Number of z points for convolution integrals (default: 500).
        More points = better accuracy but slower.
    morp : str, optional
        Distribution type (default: 'plus'). Options are 'plus'
        (ΔT q⁺ = ΔT u + ΔT d) or 'minus' (ΔT q⁻ = ΔT u - ΔT d).
    order : int, optional
        Perturbative order (default: 2). Use 1 for LO or 2 for NLO.
    logScale : bool, optional
        Use logarithmic spacing for z points (default: False).
        Recommended for peaked PDFs.
    verbose : bool, optional
        Print progress (time step count) if True (default: False).
    Q0_2_a : float, optional
        Reference scale Q₀² where αs is known, in GeV² (default: 91.1876²).
        Only used if alpha_num=True.
    a0 : float, optional
        Reference coupling αs(Q0_2_a)/(4π) (default: 0.118/(4π)).
        Only used if alpha_num=True.
    alpha_num : bool, optional
        Use numerical ODE evolution for αs if True (default: True).
        If False, uses analytical approximation.
    
    Returns
    -------
    ndarray
        Evolved PDF as 2D array [[x, x*f_evolved(x)], ...]. Shape: (n+1, 2)
        where n is the number of input points.
    
    Examples
    --------
    >>> import numpy as np
    >>> from tparton.t_evolution import evolve
    >>> x = np.linspace(0, 1, 100)
    >>> pdf_in = x * (1-x)**3  # x*f(x) format
    >>> pdf_out = evolve(pdf_in, Q0_2=4.0, Q2=100.0, n_t=200, n_z=1000)
    >>> x_out, xf_out = pdf_out[:, 0], pdf_out[:, 1]
    
    See Also
    --------
    m_evolution.evolve : Mellin moment method (Vogelsang)
    integrate : Performs convolution at a single x value
    splitting : Evaluates splitting functions
    
    References
    ----------
    - Hirai, M., Kumano, S., & Saito, N. (1998). Comput. Phys. Commun. 111, 150-160
    - Sha, C.M. & Ma, B. (2025). arXiv:2409.00221
    """
    if pdf.shape[-1] == 1:
        # If only the x*pdf(x) values are supplied, assume a linear spacing from 0 to 1
        xs = np.linspace(0, 1, len(pdf))
    else:
        # Otherwise split the input array
        xs, pdf = pdf[:, 0], pdf[:, 1]

    sign = 1 if morp == 'plus' else -1
    lnlam = 2 * np.log(l_QCD)

    # Calculate the color constants
    _, CF, Tf, beta0, beta1 = constants(CG, n_f)

    # Define the (log) starting and ending energy scales squared
    tmin = np.log(Q0_2)
    tmax = np.log(Q2)

    # Define the timepoints between those energy scales at which Eq. (1) will be integrated
    ts = np.linspace(tmin, tmax, n_t)
    
    def _beta_ode(x, a):
        """QCD beta-function ODE for running coupling.

        Parameters
        ----------
        x : float
            Log-energy variable (t = ln Q²). Not used explicitly.
        a : float
            Coupling `a = α_s/(4π)` at scale `x`.

        Returns
        -------
        float
            Time derivative da/dt according to LO/NLO beta function.
        """
        return -beta0 * a * a - (beta1 * a * a * a if order == 2 else 0)
    ode = _beta_ode

    # SciPy requires that the times be monotonically increasing or decreasing
    less = ts < np.log(Q0_2_a)
    ts_less = ts[less]
    ts_greater =ts[~less]

    # For energies strictly below the reference energy, evolve toward lower t
    alp2pi_num_less = odeint(ode, a0, [np.log(Q0_2_a)] + list(ts_less[::-1]), tfirst=True).flatten() * 2
    alp2pi_num_less = alp2pi_num_less[-1:0:-1]
    # For energies strictly above the reference energy, evolve toward higher t
    alp2pi_num_greater = odeint(ode, a0, [np.log(Q0_2_a)] + list(ts_greater), tfirst=True).flatten() * 2
    # Combine the alpha / 2 pi in increasing order of energy scale
    alp2pi_num_greater = alp2pi_num_greater[1:]
    alp2pi_num = list(alp2pi_num_less) + list(alp2pi_num_greater)
    if alpha_num:
        # Use the numerically evolved alpha_S / 2 pi
        alp2pi_use = alp2pi_num
    else:
        # Use the approximate analytical expression for alpha_S in Eq. (4)
        alp2pi_use = alp2pi(ts, lnlam, order, beta0, beta1)

    # Euler integration of Eq. (1) by a small timestep dt
    dt = (tmax - tmin) / n_t
    res = np.copy(pdf)

    for i, alp in enumerate(alp2pi_use):
        if verbose:
            print(i+1, ' of ', len(ts), 'time steps')
        # Perform the convolution at each x using (possibly log-scaled) z integration points
        inc = np.array([
            integrate(
                res,
                index,
                z_grid(xs[index], n_z, logScale),
                alp,
                order,
                CF,
                sign,
                CG,
                Tf,
                xs,
            )
            for index in range(1, len(xs) - 1)
        ])
        # Ensure that x*pdf(x) = 0 at x = 0 and x = 1
        inc = np.pad(inc, 1)
        res += dt * inc * alp
    return np.stack((xs, res), axis=1)

Evolve transversity PDF using the direct integration method (Hirai).

This is the main function for PDF evolution using direct numerical integration of the DGLAP equation (Eq. 1). More robust for peaked PDFs but slower than Mellin method.

This method: 1. Discretizes t = ln(Q²) into n_t Euler steps 2. At each step, computes the convolution integral (Eq. 19):

   f̃(x) ⊗ g(x) = ∫ dx̃ f̃(x/z) g(z)

using Simpson's rule with n_z integration points 3. Updates PDF using forward Euler: f̃(t + dt) ≈ f̃(t) + dt·f̃'(t)

Parameters

pdf : ndarray

Input PDF as xf(x). Can be 1D array (values at x evenly spaced on [0, 1]) or 2D array ([[x0, x0f(x0)], [x1, x1*f(x1)], …]).

Q0_2 : float, optional

Initial energy scale squared in GeV² (default: 0.16).

Q2 : float, optional

Final energy scale squared in GeV² (default: 5.0).

l_QCD : float, optional

QCD scale parameter Λ in GeV (default: 0.25). Only used if alpha_num=False.

n_f : int, optional

Number of active quark flavors (default: 5).

CG : float, optional

Number of colors, NC (default: 3).

n_t : int, optional

Number of Euler time steps (default: 100). More steps = better accuracy but slower.

n_z : int, optional

Number of z points for convolution integrals (default: 500). More points = better accuracy but slower.

morp : str, optional

Distribution type (default: 'plus'). Options are 'plus' (ΔT q⁺ = ΔT u + ΔT d) or 'minus' (ΔT q⁻ = ΔT u - ΔT d).

order : int, optional

Perturbative order (default: 2). Use 1 for LO or 2 for NLO.

logScale : bool, optional

Use logarithmic spacing for z points (default: False). Recommended for peaked PDFs.

verbose : bool, optional

Print progress (time step count) if True (default: False).

Q0_2_a : float, optional

Reference scale Q₀² where αs is known, in GeV² (default: 91.1876²). Only used if alpha_num=True.

a0 : float, optional

Reference coupling αs(Q0_2_a)/(4π) (default: 0.118/(4π)). Only used if alpha_num=True.

alpha_num : bool, optional

Use numerical ODE evolution for αs if True (default: True). If False, uses analytical approximation.

Returns

ndarray

Evolved PDF as 2D array [[x, x*f_evolved(x)], …]. Shape: (n+1, 2) where n is the number of input points.

Examples

>>> import numpy as np
>>> from tparton.t_evolution import evolve
>>> x = np.linspace(0, 1, 100)
>>> pdf_in = x * (1-x)**3  # x*f(x) format
>>> pdf_out = evolve(pdf_in, Q0_2=4.0, Q2=100.0, n_t=200, n_z=1000)
>>> x_out, xf_out = pdf_out[:, 0], pdf_out[:, 1]

See Also

evolve()

Mellin moment method (Vogelsang)

integrate

Performs convolution at a single x value

splitting

Evaluates splitting functions

References

• Hirai, M., Kumano, S., & Saito, N. (1998). Comput. Phys. Commun. 111, 150-160
• Sha, C.M. & Ma, B. (2025). arXiv:2409.00221