Bivariate incomplete-Bessel kernels for the first nonlinear Vlasov-Maxwell response

This paper extends a recently developed incomplete-Bessel function formulation from linear to nonlinear plasma physics by deriving the first nonlinear Vlasov-Maxwell response in terms of a newly introduced bivariate incomplete-Bessel kernel, which naturally arises from the characteristic integral and offers a more compact alternative to traditional double cyclotron-harmonic expansions.

The Gist

Imagine a plasma as a giant, chaotic dance floor filled with charged particles (electrons and ions). These particles are constantly spinning in circles because of a magnetic field, much like children spinning on a playground merry-go-round.

When scientists want to understand how these particles react to a wave of energy (like a radio wave or a laser) passing through the plasma, they have to solve a very complex mathematical puzzle called the Vlasov–Maxwell equation.

The Old Way: The "Infinite List" Problem

For decades, the standard way to solve this puzzle was to break the spinning motion down into a list of "harmonics" (like musical notes).

• The Linear Problem (One Wave): When just one wave hits the plasma, the math produces a long list of terms involving special functions called Bessel functions. It's like trying to describe a spinning top by listing every single wobble it makes. It works, but the list gets very long and messy.
• The Nonlinear Problem (Two Waves): Things get much worse when two waves interact. The math requires a "double list." You have to sum up the interactions of the first wave and the second wave simultaneously. The paper describes this as a "double cyclotron-harmonic expansion."
• The Result: The formulas become incredibly long, with nested sums and complicated denominators. It's like trying to write a recipe that requires you to list every single grain of salt and every possible variation of a spice before you can even start cooking. It's hard to read and even harder to use on a computer.

The New Way: The "Orbit Kernel" Shortcut

In this paper, the author, Roberto Ricci, proposes a smarter way to look at the problem. Instead of breaking the spinning motion into a list of notes right away, he keeps the motion as a continuous "orbit" (the actual path the particle takes).

1. The Linear Shortcut: In a previous study, the author showed that for a single wave, you can skip the long list and use a single, compact mathematical object called an incomplete-Bessel function (let's call it a "G-function"). This function captures the entire spinning path in one neat package.
2. The Nonlinear Leap: This paper takes that idea one step further. When two waves interact, the math usually requires a messy double sum. Ricci shows that if you keep the orbits unbroken, the interaction naturally creates a new, slightly more complex object: a bivariate incomplete-Bessel function.
   • Think of the linear version as a single "orbit kernel."
   • The nonlinear version is a "bivariate orbit kernel." It's a single mathematical tool that holds the information of both the outer wave and the inner wave's response, wrapped together.

The Analogy: The Nested Russian Doll

Imagine the problem is like a set of Russian nesting dolls.

• The Old Method: To understand the big doll, you have to open it, take out the next one, write down its size, open that one, write down its size, and so on. If you have two layers of interaction, you have to do this twice, creating a massive pile of notes.
• The New Method: Ricci suggests looking at the whole stack of dolls as a single, unified shape. He introduces a new "bivariate" tool that describes the entire stack at once. You don't need to open the dolls to see the pattern; the tool captures the relationship between the inner and outer layers instantly.

Why This Matters

The paper claims three main benefits of this new approach:

1. It's the "Natural" Way: The author argues that this bivariate function isn't just a trick; it's the natural mathematical shape that emerges when you follow the actual path (orbit) of the particles. The old "double sum" method is just what happens when you force this natural shape into a long list of numbers.
2. It Recovers the Old Results: The author proves that if you do expand this new bivariate function into a list, you get exactly the same complicated formulas that everyone has been using for years (the Liu–Tripathi formula). This confirms the new method is correct.
3. It's Better for Computers: While the paper doesn't promise a specific new technology (like a new fusion reactor), it suggests that for computer simulations, this new method is much cleaner. Instead of summing thousands of terms that might cancel each other out (causing numerical errors), a computer can calculate this single "bivariate kernel" and its derivatives. It's like replacing a messy spreadsheet with a single, elegant formula.

Summary

In simple terms, this paper says: "Stop trying to describe the complex dance of plasma particles by listing every single step. Instead, use a new mathematical 'lens' (the bivariate incomplete-Bessel function) that captures the whole dance in one go. It's cleaner, it's mathematically natural, and it makes the heavy lifting of calculating plasma behavior much more organized."

Technical Summary

Technical Summary: Bivariate Incomplete-Bessel Kernels for the First Nonlinear Vlasov–Maxwell Response

Problem Statement
The kinetic theory of nonlinear wave coupling in magnetised plasmas traditionally relies on a weakly nonlinear expansion of the Vlasov–Maxwell system about a homogeneous equilibrium. The standard approach expresses the plasma response as a double cyclotron-harmonic expansion. While explicit, this representation becomes algebraically cumbersome, particularly when electromagnetic waves propagate at arbitrary angles to the background magnetic field. In the first nonlinear correction, the response involves two coupled harmonic indices: one arising from the outer characteristic integration of the output mode and another from the inner linear response acting as a source. This structure leads to complex double sums of Bessel functions and nested resonance denominators, as seen in classical formulations by Liu and Tripathi (1986) and subsequent works.

Methodology
This paper proposes an alternative formulation that avoids the immediate expansion of the Larmor phase into a harmonic series. Instead, the authors retain the phase along the unperturbed orbit and evaluate the characteristic integrals directly.

1. Linear Baseline: The authors recall a recent linear formulation (Ricci 2026) where the characteristic integral is evaluated in terms of the incomplete-Bessel function $G_\mu(z, \psi)$. This function replaces the infinite sum of cyclotron harmonics before the susceptibility tensor is projected.
2. Nonlinear Extension: The first nonlinear Vlasov equation is solved using the same characteristic method. The nonlinear source term consists of the Lorentz force of one mode acting on the linear response of another. By expressing the inner linear response in its incomplete-Bessel form, the outer characteristic integral acts on a shifted $G_{\mu_q}$ response.
3. Kernel Derivation: This process naturally generates a bivariate incomplete-Bessel function, denoted as $G^{(r)}_{\mu,\nu}(z, \psi; w, \chi)$. This function serves as the orbit-resolvent kernel for the first nonlinear problem, where the integer $r$ accounts for finite harmonic shifts generated by transverse velocity factors and gyrophase derivatives.
4. Recovery of Classical Results: The authors derive the nonlinear distribution function explicitly in terms of these bivariate functions. They then demonstrate that expanding these functions using the Jacobi–Anger formula and Turkin's function recovers the traditional double-harmonic Liu–Tripathi formula, thereby validating the new formulation.
5. Susceptibility Projection: Finally, the paper derives the nonlinear susceptibility tensor by projecting the bivariate distribution onto velocity moments. This is achieved through "projected bivariate kernels," which are angular Fourier projections of the bivariate functions, avoiding the need to return to the full double-harmonic expansion.

Key Contributions and Results

• Definition of Bivariate Kernels: The paper introduces $G^{(r)}_{\mu,\nu}$ as the natural orbit integral generated by the first nonlinear characteristic problem. It establishes that this function is not merely a notational convenience but the fundamental object arising from the iteration of the characteristic method.
• Structural Reorganization: The nonlinear distribution function is expressed compactly using differential operators acting on these bivariate kernels. The formulation isolates the angular contractions into a small set of reusable special-function blocks.
• Mathematical Identities: A comprehensive set of identities for the bivariate functions is derived, including superposition properties, harmonic expansions, recurrence relations, and differential equations (e.g., the outer Bessel equation and inner Turkin recurrence).
• Consistency Verification: The paper provides a rigorous derivation showing that the bivariate formulation is mathematically equivalent to the classical double-harmonic approach. The traditional formula is recovered by expanding the bivariate function and applying Graf's addition theorem to contract the intermediate harmonic indices.
• Explicit Susceptibility Formula: The nonlinear susceptibility tensor is presented as a finite sum over nine elementary bivariate blocks (indexed by $a, b \in \{-1, 0, 1\}$), followed by a two-dimensional velocity integral. The angular dependence is encapsulated in projected kernels $K^{(r,a)}_N$ and their derivatives.

Significance and Claims
The paper claims that the primary advantage of this formulation is both structural and potentially computational:

• Structural: It identifies the specific special-function object (the bivariate incomplete-Bessel kernel) that underlies the first nonlinear cyclotron response, organizing the theory around the geometry of the characteristic flow rather than the algebraic complexity of harmonic sums.
• Computational: The authors suggest that for numerical implementation, the bivariate formulation offers a more structured target. Instead of handling long, oscillatory double sums of Bessel functions with shifted indices, one can tabulate and reuse the projected kernels and their derivatives. This flexibility may be advantageous in regimes where direct harmonic summation is slowly convergent or numerically delicate due to cancellations.
• Scalability: The authors note that this approach extends recursively to higher orders. The $n$-th order distribution is expected to be expressible in terms of $n$-variate orbit kernels, generated by adding one characteristic time integral for each nonlinear interaction.

The paper maintains a modest stance, acknowledging that the numerical advantage is not a universal reduction in computational cost but depends on specific plasma parameters (Larmor radii, detunings, damping). It explicitly states that the construction applies to the fixed-characteristic perturbative hierarchy and does not account for nonlinear orbit corrections (such as trapping or ponderomotive drifts) that would modify the characteristic operator itself.