multiloop_sunset_equal_masses.sage

# This files provides the routines for extracting the coefficients of the sunset graphs 
var('u')

def initialize_sunset_variables_equal_mass(looporder, Nmax):
    """
    Initialize sunset diagram variables for perturbative calculations.
    
    Parameters:
    -----------
    looporder : int
        The loop order for the calculation
    Nmax : int
        Maximum number of terms
    max_log_power : int
        Maximum logarithmic power for symbolic variables
    PFEqualMass : list or array
        Input for OA function (equal mass propagator factors)
    
    Returns:
    --------
    dict : Dictionary containing all initialized variables
    """
    # Initialize main sunset variables
    max_log_power=looporder

    Lsunset = OA(PFEqualMass[looporder])
    
    Persunset = Period([1 for i in range(looporder + 1)], Nmax)
    
    Qsunset = Q([1 for i in range(looporder + 1)], Nmax)*(-1)
    
    Rsunset = R([1 for i in range(looporder + 1)], Nmax)+I*pi
    
    # Create symbolic variables for each log power
    var_groups = {}
    for log_pow in range(max_log_power):
        var_name = chr(65 + log_pow)  # A, B, C, D, ...
        var_groups[log_pow] = [SR.var(f'{var_name}_{n}') for n in range(1, Nmax)]
    
    var_groups[looporder] = [SR.var(f'Frob_{n}') for n in range(0, looporder)]

    Frobenius_basis_tmp=(Lsunset.annihilator_of_composition(1/t).local_basis_expansions(0,Nmax))
    Frobenius_basis=sum(var_groups[looporder][n]*SR(sage_eval(str(Frobenius_basis_tmp[looporder-1-n]), locals={'t': t, 'log': log})) for n in range(looporder)).subs(t=1/t).simplify().subs(log(t)==log(t)-I*pi);

    # Return all variables in a dictionary
    return {
        'Lsunset': Lsunset,
        'Persunset': Persunset,
        'Qsunset': Qsunset,
        'Rsunset': Rsunset,
        'var_groups': var_groups,
        'Frobenius': Frobenius_basis
    }

# Example usage:
# sunset_vars = initialize_sunset_variables_equal_mass(looporder=2, Nmax=10)
# Lsunset = sunset_vars['Lsunset']
# var_groups = sunset_vars['var_groups']


def build_equation_system_from_PF(log_series, var_groups, leading_log_ser, Nmax, max_log_power, return_expressions=False):
    """
    Build system of equations for each logarithmic power.
    
    Parameters:
    -----------
    log_series : dict
        Dictionary of log series from generate_log_series()
    var_groups : dict
        Dictionary of symbolic variable groups indexed by log power
    leading_log_ser : list
        Coefficient series for the constant term
    Nmax : int
        Maximum number of terms
    max_log_power : int
        Maximum logarithmic power
    return_expressions : bool, optional
        If True, return both systems and expressions dict
    
    Returns:
    --------
    list or tuple : System of equations, optionally with expressions dictionary
    """
    systems = []
    expressions = {} if return_expressions else None
    
    for target_log_pow in range(max_log_power + 1):
        # Start with the constant term contribution
        if target_log_pow < len(leading_log_ser):
            expr = leading_log_ser[target_log_pow][0]
        else:
            expr = 0
        
        # Add contributions from each variable group
        for var_log_pow in range(max_log_power):
            vars_list = var_groups[var_log_pow]
            
            # Determine if this variable group contributes to this log power
            # A variables (log_pow=0) only contribute to target_log_pow=0
            # B variables (log_pow=1) contribute to target_log_pow=0,1
            # etc.
            if var_log_pow >= target_log_pow:
                for n in range(Nmax - 1):
                    if var_log_pow == 0:
                        # No-log terms only contribute to no-log equation
                        if target_log_pow == 0:
                            expr += vars_list[n] * log_series[var_log_pow][n]
                    else:
                        # Extract the appropriate log coefficient
                        series_data = log_series[var_log_pow][n]
                        # Find the entry with the right log power
                        for coeff_data in series_data:
                            if coeff_data[1] == target_log_pow:
                                expr += vars_list[n] * coeff_data[0]
                                break
        
        if return_expressions:
            expressions[target_log_pow] = expr
        
        # Extract coefficients of u to form equations
        equations = [x[0] for x in expr.coefficients(u)]
        systems.extend(equations)
    
    if return_expressions:
        return systems, expressions
    else:
        return systems


# Example usage:
# systems = build_equation_system_from_PF(log_series, var_groups, Lsumsunset0_coefflog_ser, 
#                                  Nmax=10, max_log_power=3)
# 
# # Or with expressions for debugging:
# systems, exprs = build_equation_system_from_PF(log_series, var_groups, Lsumsunset0_coefflog_ser, 
#                                         Nmax=10, max_log_power=3, return_expressions=True)