Package `tparton`

tParton: Evolution of transversity parton distribution functions

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

Direct integration (Hirai method): Numerically solves the DGLAP equations using discretized grids in x and Q².
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')

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

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. (25) with splitting function moments from
       Eq. (27) for LO and Eq. (28) 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 in the tilde convention x*f(x). Accepted formats are a 1D
        array of shape (N,) or a single-column array of shape (N, 1), both
        taken as the x*f(x) values at x evenly spaced on [0, 1] inclusive, or
        a two-column array of shape (N, 2) holding
        [[x0, x0*f(x0)], [x1, x1*f(x1)], ...]. Note that the second column is
        x*f(x), not f(x).
    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
        Charge-conjugation combination of a single quark flavor (default:
        'minus'). Options are 'plus' (ΔT q⁺ = ΔT q + ΔT q̄) or 'minus'
        (ΔT q⁻ = ΔT q − ΔT q̄). This selects the continuation branch η = ±1
        of the splitting function moments; it does not refer to a flavor
        combination such as ΔT u ± ΔT d. A flavor difference such as
        ΔT u − ΔT d is itself a valid non-singlet input, and is evolved with
        whichever branch its charge-conjugation structure calls for.
    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 a 2D array [x, x*f_evolved(x)], arranged by row, so
        that result[0] is the x grid and result[1] the evolved x*f(x) values.
        The shape is (2, n_x+2) when n_x > 0, the two extra points being the
        padded endpoints x = 0 and x = 1; when n_x <= 0 the input x grid is
        reused and the shape is (2, N).
    
    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]
    >>> len(x_out) == len(xf_out)
    True
    
    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
    """

    # Normalize every supported input format to a pair of 1D arrays
    xs, pdf = split_pdf_input(pdf)
    
    # 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. (32)
    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 set beta1 to 0 so that the
        # LO running coupling is used. The NLO splitting function moment is
        # dropped inside evolveMoment() via its `order` argument; it cannot be
        # suppressed by a local rebinding here, because evolveMoment() resolves
        # the module-level function, nor by rebinding the module global, which
        # would corrupt later NLO calls in the same process.
        beta1 = 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. (21)
    pdf_m = lambda s: mellin(pdf, s)
    # A function representing the resulting evolved moments, Eq. (25)
    pdf_evolved_m = lambda s: mpc(evolveMoment(s, pdf_m(s), alpha_S_Q0_2, alpha_S_Q2, beta0, beta1, eta, CF, NC, Tf, order))
    # Perform Mellin inversion on the evolved moments, Eq. (37)
    # mp.re() takes the real part for both mpf and mpc returns; mpf.__complex__
    # was removed in mpmath 1.4, so do not rely on it here.
    pdf_evolved = np.array([float(mp.re(inv_mellin(pdf_evolved_m, x, degree=degree, verbose=verbose))) 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
    # Row-major (2, N): row 0 is x, row 1 is x*f_evolved(x). The notebooks and the
    # stored .npz reference data index this as result[0] / result[1].
    pdf_evolved = np.stack((xs, np.array(xs) * np.array(pdf_evolved)))
    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. (25) with splitting function moments from Eq. (27) for LO and Eq. (28) 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 in the tilde convention xf(x). Accepted formats are a 1D array of shape (N,) or a single-column array of shape (N, 1), both taken as the xf(x) values at x evenly spaced on [0, 1] inclusive, or a two-column array of shape (N, 2) holding [[x0, x0f(x0)], [x1, x1f(x1)], …]. Note that the second column is x*f(x), not f(x).
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: Charge-conjugation combination of a single quark flavor (default: 'minus'). Options are 'plus' (ΔT q⁺ = ΔT q + ΔT q̄) or 'minus' (ΔT q⁻ = ΔT q − ΔT q̄). This selects the continuation branch η = ±1 of the splitting function moments; it does not refer to a flavor combination such as ΔT u ± ΔT d. A flavor difference such as ΔT u − ΔT d is itself a valid non-singlet input, and is evolved with whichever branch its charge-conjugation structure calls for.
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 a 2D array [x, xf_evolved(x)], arranged by row, so that result[0] is the x grid and result[1] the evolved xf(x) values. The shape is (2, n_x+2) when n_x > 0, the two extra points being the padded endpoints x = 0 and x = 1; when n_x <= 0 the input x grid is reused and the shape is (2, N).

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]
>>> len(x_out) == len(xf_out)
True

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 in the tilde convention x*f(x). Accepted formats are a 1D
        array of shape (N,) or a single-column array of shape (N, 1), both
        taken as the x*f(x) values at x evenly spaced on [0, 1] inclusive, or
        a two-column array of shape (N, 2) holding
        [[x0, x0*f(x0)], [x1, x1*f(x1)], ...]. Note that the second column is
        x*f(x), not f(x).
    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
        Charge-conjugation combination of a single quark flavor (default:
        'plus'). Options are 'plus' (ΔT q⁺ = ΔT q + ΔT q̄) or 'minus'
        (ΔT q⁻ = ΔT q − ΔT q̄). This selects which transversity splitting
        function ΔT P_q± is used; it does not refer to a flavor combination
        such as ΔT u ± ΔT d. A flavor difference such as ΔT u − ΔT d is itself
        a valid non-singlet input, and is evolved with whichever branch its
        charge-conjugation structure calls for.
    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 a 2D array [x, x*f_evolved(x)], arranged by row, so
        that result[0] is the x grid and result[1] the evolved x*f(x) values.
        The shape is (2, N), where N is the number of input points: this
        method evolves the input grid in place and does not resample it.
    
    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]
    >>> len(x_out) == len(xf_out)
    True
    
    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
    """
    # Normalize every supported input format to a pair of 1D arrays
    xs, pdf = split_pdf_input(pdf)

    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.
    # These are the n_t left endpoints of n_t equal sub-intervals, consistent with the
    # Euler step dt = (tmax - tmin) / n_t used below.
    ts = np.linspace(tmin, tmax, n_t + 1)[:-1]
    
    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

    if alpha_num:
        # Use the numerically evolved alpha_S / 2 pi, integrating outward from
        # the reference scale t_ref = ln(Q0_2_a) at which alpha_S is known.
        #
        # SciPy's odeint requires a monotonic time array, whereas the Euler grid
        # ts above is decreasing whenever the evolution runs downward in Q2. So
        # the sample points below and above t_ref are integrated separately,
        # each sorted by increasing distance from t_ref, and the results are
        # scattered back onto the original ordering of ts. Sorting ts itself and
        # evolving on the sorted grid would instead misalign the coupling with
        # the Euler steps.
        t_ref = np.log(Q0_2_a)
        alp2pi_use = np.empty_like(ts)
        below = ts < t_ref
        for mask, direction in ((below, -1.0), (~below, 1.0)):
            idx = np.flatnonzero(mask)
            if len(idx) == 0:
                continue
            # Order these sample points so that t moves monotonically away from
            # t_ref, which is the first point of the integration
            idx = idx[np.argsort(direction * ts[idx])]
            sol = odeint(ode, a0, np.concatenate(([t_ref], ts[idx])), tfirst=True).flatten()
            alp2pi_use[idx] = sol[1:] * 2
    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
    # Row-major (2, N): row 0 is x, row 1 is x*f_evolved(x). The notebooks and the
    # stored .npz reference data index this as result[0] / result[1].
    return np.stack((xs, res))

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 in the tilde convention xf(x). Accepted formats are a 1D array of shape (N,) or a single-column array of shape (N, 1), both taken as the xf(x) values at x evenly spaced on [0, 1] inclusive, or a two-column array of shape (N, 2) holding [[x0, x0f(x0)], [x1, x1f(x1)], …]. Note that the second column is x*f(x), not f(x).
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: Charge-conjugation combination of a single quark flavor (default: 'plus'). Options are 'plus' (ΔT q⁺ = ΔT q + ΔT q̄) or 'minus' (ΔT q⁻ = ΔT q − ΔT q̄). This selects which transversity splitting function ΔT P_q± is used; it does not refer to a flavor combination such as ΔT u ± ΔT d. A flavor difference such as ΔT u − ΔT d is itself a valid non-singlet input, and is evolved with whichever branch its charge-conjugation structure calls for.
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 a 2D array [x, xf_evolved(x)], arranged by row, so that result[0] is the x grid and result[1] the evolved xf(x) values. The shape is (2, N), where N is the number of input points: this method evolves the input grid in place and does not resample it.

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]
>>> len(x_out) == len(xf_out)
True

References

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

Basic Usage

Navigation

See Also

Sub-modules

Functions

Parameters

Returns

Notes

Examples

See Also

References

Parameters

Returns

Examples

See Also

References