A first study of strong isospin breaking effects in lattice QCD using truncated polynomials

This paper introduces a new approach utilizing automatic differentiation to compute high-order derivatives of lattice QCD observables, specifically focusing on strong isospin-breaking effects and the propagation of these derivatives through the conjugate gradient algorithm.

The Gist

Imagine you are a master chef trying to perfect a complex recipe for a giant, cosmic cake (a subatomic particle called a kaon). The recipe calls for two very specific ingredients: "Up-quark flour" and "Down-quark flour."

In the standard version of the universe, these two flours are almost identical twins. They weigh exactly the same, and they taste the same. However, in the real universe, the Down-quark is just a tiny, tiny bit heavier than the Up-quark. This tiny difference is called Strong Isospin Breaking.

Because this difference is so small, it's incredibly hard to measure its effect on the final cake. If you try to bake the cake with the "real" slightly different flours from scratch, the process is slow, expensive, and computationally heavy (like trying to bake a cake in a kitchen that is constantly shaking).

The Old Way: The "Manual Calculator" Method

Previously, scientists used a method called RM123. Think of this as a manual calculator.

1. They bake a cake using the "perfect twin" flours (where Up and Down are identical).
2. Then, they manually calculate what would happen if they swapped just a pinch of Up-flour for a pinch of Down-flour.
3. They do this math by hand (or by writing specific computer code) for the first pinch, then the second pinch, then the third.
4. The Problem: Every time they wanted to go one step further (to calculate the effect of the second or third pinch), they had to write a whole new, complicated set of instructions. It was like having to re-invent the wheel every time you wanted to add a new spice.

The New Way: The "Smart Blender" (Truncated Polynomials)

This paper introduces a new tool: Truncated Polynomials combined with Automatic Differentiation.

Imagine instead of a manual calculator, you have a Smart Blender.

• You don't just put in the flour; you put in the flour plus a tiny, invisible "tag" that says, "If I change the weight of this ingredient by 0.0001%, here is exactly how the cake changes."
• When you blend the ingredients (solve the complex math equations), the blender doesn't just give you the cake; it automatically calculates the entire chain reaction of changes for you.
• It doesn't matter if you want to know the effect of one pinch or a thousand pinches. The blender does the heavy lifting automatically.

The Big Challenge: The "Conveyor Belt"

The most difficult part of baking this cosmic cake is a step called the Conjugate Gradient Algorithm. Imagine this as a conveyor belt that has to sort millions of tiny crumbs perfectly before the cake can be finished.

• Usually, the conveyor belt stops when the crumbs are "good enough" sorted.
• The authors of this paper had to figure out how to make the conveyor belt stop perfectly even when it's carrying those invisible "change tags" (the derivatives).
• They proved that their new "Smart Blender" setup works perfectly. The conveyor belt stops at the right time, and the final result is just as accurate as the old manual method, but it happens automatically.

What Did They Find?

The team tested this new method on a specific particle (the Kaon).

1. They used their "Smart Blender" to calculate how the particle's mass changes due to the tiny weight difference between Up and Down quarks.
2. They compared their result to the old "Manual Calculator" method (RM123).
3. The Result: The numbers matched almost perfectly (the difference was less than one part in ten million!).

Why Does This Matter?

This is a huge breakthrough for a few reasons:

• Speed and Ease: Scientists no longer need to write complex, custom code for every new level of precision. They just turn on the "Automatic Differentiation" switch.
• Future Proofing: This method can be used for any parameter, not just quark masses. It can help fix mistakes in simulations or even calculate the effects of electromagnetism on particles.
• Precision: It allows physicists to understand the universe with a level of detail that was previously too difficult or expensive to achieve.

In short: The authors built a "universal translator" for physics simulations. Instead of manually translating every single sentence of a complex math problem, they built a tool that translates the whole story at once, automatically, and with perfect accuracy. This makes studying the subtle secrets of the universe much faster and easier.

Technical Summary

1. Problem Statement

In Lattice Quantum Chromodynamics (LQCD), calculating physical observables often requires accounting for strong isospin breaking (IB) effects, primarily arising from the mass difference between up ($m_u$) and down ($m_d$) quarks. These effects are crucial for precision calculations of hadronic observables (e.g., mass splittings, the Hadronic Vacuum Polarization contribution to $g-2$).

The challenge lies in efficiently computing these effects:

• Exact methods (simulating non-degenerate quarks directly) are computationally expensive, requiring complex algorithms like Rational Hybrid Monte Carlo (RHMC).
• Approximate methods (like the RM123 method) treat isospin breaking as a Taylor expansion in the mass splitting $\Delta m = (m_d - m_u)/2$. While efficient, the RM123 method requires the manual derivation, implementation, and testing of specific Feynman diagrams for each order of the expansion. As the truncation order increases, the number of diagrams grows combinatorially, making higher-order calculations tedious and error-prone.

The authors aim to develop a framework that automates the calculation of these derivatives to arbitrary orders without manual diagrammatic derivation.

2. Methodology

The paper proposes a novel approach using Automatic Differentiation (AD) via Truncated Polynomials to generalize the RM123 method.

A. Truncated Polynomials

Instead of standard floating-point numbers, the authors use truncated polynomials (a form of forward-mode AD). A variable $x$ is represented as a polynomial:
$$\tilde{x} = x^{(0)} + x^{(1)}\Delta m + x^{(2)}\Delta m^2 + \dots + x^{(K)}\Delta m^K$$
By overloading mathematical operations to act on these polynomials, any differentiable function $f(\tilde{x})$ automatically computes the Taylor expansion coefficients (derivatives) of the output up to order $K$.

B. Reweighting Framework

The expectation value of an observable $\langle O \rangle_{\Delta m}$ is computed by reweighting an isospin-symmetric ensemble ($\Delta m = 0$):
$$\langle O \rangle_{\Delta m} = \frac{\langle O(\Delta m) W_{IB}(\Delta m) \rangle_0}{\langle W_{IB}(\Delta m) \rangle_0}$$

• Valence Contribution ($O(\Delta m)$): The correlation function (e.g., Kaon) where quark masses depend on $\Delta m$.
• Sea Contribution ($W_{IB}$): The reweighting factor involving the determinant of the Dirac operator.

By replacing $\Delta m$ with a truncated polynomial $\tilde{\Delta m}$ (where the coefficient of $\Delta m$ is 1 and others are 0), the entire reweighting equation yields the full Taylor expansion of the observable automatically.

C. Integration with Conjugate Gradient (CG)

A critical technical hurdle is solving the Dirac equation $D\psi = \eta$ iteratively using the Conjugate Gradient algorithm. The authors investigate whether AD propagates correctly through this iterative solver.

• Stopping Criterion Adaptation: Standard CG stops when the residual norm is below a tolerance. With truncated polynomials, the residual becomes a polynomial. The authors propose enforcing the tolerance condition order-by-order on the polynomial coefficients to ensure convergence of all derivatives.

3. Key Contributions

1. Automation of RM123: The paper demonstrates that truncated polynomials can fully automate the generation of isospin-breaking corrections, eliminating the need for manual derivation of higher-order diagrams.
2. Validation of AD in Iterative Solvers: It provides a proof-of-concept that automatic differentiation can successfully propagate derivatives through the Conjugate Gradient algorithm, a non-trivial iterative process.
3. Implementation: The authors implemented this framework using Julia (specifically FormalSeries.jl for polynomials and LatticeGPU.jl for GPU-accelerated lattice simulations).

4. Results

The study was performed on the A654 ensemble (CLS initiative, $N_f=2+1$ Wilson fermions, $48 \times 24^3$ volume).

• Observable: The derivative of the charged Kaon mass with respect to the mass splitting, $\partial M_{K^+} / \partial \Delta m$.
• Comparison: The authors compared the first-order derivative obtained via:
  1. Truncated Polynomials (Automatic Differentiation): Solving the Dirac equation with polynomial inputs.
  2. Standard RM123: Explicitly computing the specific first-order diagram (Eq. 16 in the paper).
• Outcome:
  • The ratio of the first-order to leading-order correlation function was fitted to extract the mass derivative, yielding $M^{(1)}_{K^+} = -3.86(52)$.
  • Agreement: The relative difference between the AD method and the explicit RM123 method was found to be of the order of $10^{-7}$.
  • This confirms that the derivatives are correctly propagated through the CG algorithm and that the stopping criterion adaptation works effectively.

5. Significance and Future Outlook

• Scalability: This method allows for the calculation of isospin-breaking effects to arbitrary orders in $\Delta m$ with minimal additional coding effort, as the algebra handles the combinatorial complexity of higher-order diagrams automatically.
• General Applicability: While focused on strong isospin breaking, the framework is general and can be applied to:
  • Electromagnetic isospin breaking.
  • Correction of mistunings in bare parameters.
  • Derivatives with respect to any action parameter.
• Next Steps: The authors plan to extend the analysis to higher orders in the mass expansion, incorporate sea-quark contributions (reweighting factors), and perform a detailed cost-benefit analysis comparing the computational overhead of AD against standard RM123 implementations.

In conclusion, this work establishes a robust, automated pathway for high-precision LQCD calculations involving isospin breaking, leveraging modern automatic differentiation techniques to overcome the limitations of manual diagrammatic expansions.