# Analytical Aspects of the Compressible Magnetohydrodynamic Equations under Leray Scaling

Edoardo Niccolai

► To cite this version:

Edoardo Niccolai. Analytical Aspects of the Compressible Magnetohydrodynamic Equations under Leray Scaling. 2025. hal-05328121

HAL Id: hal-05328121

https://hal.science/hal-05328121v1

Preprint submitted on 23 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical Aspects of the
Compressible Magnetohydrodynamic Equations
under Leray Scaling

Edoardo Niccolai

**ABSTRACT.** This paper examines the analytical framework of the compressible magnetohydrodynamic (MHD) equations, emphasizing their connection with Leray's self-similar formulation of the Navier-Stokes system. Within this perspective, the work surveys key developments primarily in the weak and variational setting, global regularity, and spectral decay for viscous and resistive flows. The analysis aims to consolidate classical and modern contributions into a unified mathematical narrative, clarifying how Leray-type self-similarity and scaling arguments persist and influence the emergence of singular or self-organized regimes.

**KEYWORDS:** Compressible magnetohydrodynamics, Leray self-similar solutions, weak solutions, global regularity, spectral decay, singularity formation.

## CONTENTS

<table><tr><td>1. Framework and Notational Setup</td><td>4</td></tr><tr><td>2. From Leray to Modern Analysis: a Quick Look</td><td>4</td></tr><tr><td>3. Full Compressible MHD Equations</td><td>4</td></tr><tr><td>4. Fully Non-local Compressible MHD Model</td><td>5</td></tr><tr><td>5. A Variational Methodology</td><td>7</td></tr><tr><td>5.1. Recasting in Weighted Lebesgue Spaces</td><td>8</td></tr><tr><td>6. Criteria for Finite-Time Singularities</td><td>9</td></tr><tr><td>7. Dissipative Regime and Asymptotic Stabilization</td><td>12</td></tr><tr><td>8. Liouville Constraints and Entropy-Based Estimates</td><td>14</td></tr><tr><td>8.1. Some Details and Nuances</td><td>16</td></tr><tr><td>8.1.1. Spectral Linearization Around the Momentum Equation</td><td>16</td></tr><tr><td>8.1.1.1. Green Closure of the Momentum Balance</td><td>19</td></tr><tr><td>8.1.1.2. Iterated Functional Hierarchy and K<sub>P,λ</sub>-Cascade</td><td>20</td></tr><tr><td>8.1.2. Regularized Determinants of Integral Operators</td><td>21</td></tr><tr><td>8.1.2.1. Resolvent and Function-operator Nested Integral Cascade</td><td>22</td></tr><tr><td>8.1.2.2. Kernel-Based Analysis and Spectral Projectors</td><td>23</td></tr><tr><td>I. Iterated Chain Representation and Fredholm Determinant Series</td><td>23</td></tr><tr><td>II. Double Composition</td><td>25</td></tr><tr><td>III. Spectral Projectors and Contour Integral Sketch</td><td>26</td></tr><tr><td>IV. Faà di Bruno hierarchies</td><td>27</td></tr><tr><td>V. Higher-Order Parametric Differentiation</td><td>28</td></tr><tr><td>VI. Vue d'Ensemble</td><td>29</td></tr><tr><td>VII. Pseudo-differential Symbolics and Microlocal Propagation</td><td>30</td></tr><tr><td>8.1.2.3. Radial Reductions, Spherical Decoupling and Mellin-Bessel via Evaluation</td><td>32</td></tr><tr><td>8.1.2.4. Decomposition of the Fredholm Trace-log</td><td>34</td></tr><tr><td>8.1.3. Two Supporting Propositions</td><td>35</td></tr><tr><td>8.1.3.1. I. On the Magnetic Divergence-Free Condition</td><td>35</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>8.1.3.2. II. On the Maxwell-Transport Cancellation</td>
      <td>35</td>
    </tr>
    <tr>
      <td>9. Leray-Type Self-Similar Ansatz</td>
      <td>35</td>
    </tr>
    <tr>
      <td>9.1. Geometric Interpretation via Lie Derivatives</td>
      <td>36</td>
    </tr>
    <tr>
      <td>9.2. Energetic and Entropic Scalings</td>
      <td>37</td>
    </tr>
    <tr>
      <td>9.3. Entropy and Energy Scaling</td>
      <td>38</td>
    </tr>
    <tr>
      <td>9.4. Self-Similar Reduction of the Full System</td>
      <td>38</td>
    </tr>
    <tr>
      <td>9.5. Further Remarks</td>
      <td>39</td>
    </tr>
    <tr>
      <td>10. Derivation of the Self-Similar ODE System</td>
      <td>40</td>
    </tr>
    <tr>
      <td>10.1. Radial Symmetry and Model Simplification</td>
      <td>40</td>
    </tr>
    <tr>
      <td>10.2. Final Form of the Self-Similar System</td>
      <td>40</td>
    </tr>
    <tr>
      <td>10.2.1. Various Enlightenment</td>
      <td>40</td>
    </tr>
    <tr>
      <td>10.2.1.1. Non-local Integral Chain</td>
      <td>41</td>
    </tr>
    <tr>
      <td>10.2.1.2. Schur Complement Expansion</td>
      <td>44</td>
    </tr>
    <tr>
      <td>10.3. Algebraic Power-Law Profiles</td>
      <td>50</td>
    </tr>
    <tr>
      <td>10.4. Barotropic Closure Conditions</td>
      <td>51</td>
    </tr>
    <tr>
      <td>10.5. Qualitative Behavior and Invariant Regions</td>
      <td>53</td>
    </tr>
    <tr>
      <td>11. Center Manifold Dynamics</td>
      <td>53</td>
    </tr>
    <tr>
      <td>11.1. Symmetry Groups and Invariant Structures</td>
      <td>53</td>
    </tr>
    <tr>
      <td>11.2. Dynamical Reformulation</td>
      <td>54</td>
    </tr>
    <tr>
      <td>11.3. Continuation Criterion in F<sub>g</sub><sup>c</sup></td>
      <td>55</td>
    </tr>
    <tr>
      <td>11.4. Taylor Expansion at r = 0</td>
      <td>56</td>
    </tr>
    <tr>
      <td>11.5. Jacobian and Local Spectral Features</td>
      <td>57</td>
    </tr>
    <tr>
      <td>11.6. Existence of a Local Center Manifold</td>
      <td>57</td>
    </tr>
    <tr>
      <td>11.7. Autonomous Dynamics and Invariant Geometry</td>
      <td>57</td>
    </tr>
    <tr>
      <td>11.8. Stable Modes and End-State Analysis</td>
      <td>58</td>
    </tr>
    <tr>
      <td>11.9. Global Matching via Invariant Manifolds</td>
      <td>58</td>
    </tr>
    <tr>
      <td>12. Stability via Lyapunov Functionals</td>
      <td>58</td>
    </tr>
    <tr>
      <td>12.1. Anatomy of the Linear Operator</td>
      <td>59</td>
    </tr>
    <tr>
      <td>12.2. Spectral Properties and Gaps</td>
      <td>59</td>
    </tr>
    <tr>
      <td>12.3. Lyapunov Control of Non-linear Behavior</td>
      <td>60</td>
    </tr>
    <tr>
      <td>12.4. Essential and Pseudo-Spectral Analysis</td>
      <td>60</td>
    </tr>
    <tr>
      <td>12.5. Sobolev Space with Weights</td>
      <td>60</td>
    </tr>
    <tr>
      <td>12.6. Variational Principles and Lyapunov Selection</td>
      <td>61</td>
    </tr>
    <tr>
      <td>13. Far-Field Behavior and Matching Theory</td>
      <td>61</td>
    </tr>
    <tr>
      <td>13.1. Calderón-Zygmund-Type Elliptic Estimate</td>
      <td>62</td>
    </tr>
    <tr>
      <td>13.2. Sharp Algebraic Decay under Finite Similarity Energy</td>
      <td>62</td>
    </tr>
    <tr>
      <td>13.3. Asymptotics for r → ∞</td>
      <td>64</td>
    </tr>
    <tr>
      <td>13.4. Spectral Admissibility</td>
      <td>64</td>
    </tr>
    <tr>
      <td>13.5. Gap Set-up for the Linearized Operator</td>
      <td>64</td>
    </tr>
    <tr>
      <td>13.5.1. Two Analytical Perspectives</td>
      <td>66</td>
    </tr>
    <tr>
      <td>13.5.1.1. I. Hardy Inequality in L<sub>ω</sub><sup>2</sup></td>
      <td>66</td>
    </tr>
    <tr>
      <td>13.5.1.2. II. Exact Drift Identity</td>
      <td>66</td>
    </tr>
    <tr>
      <td>13.6. Curvature-Driven Coercivity</td>
      <td>66</td>
    </tr>
    <tr>
      <td>13.7. Anisotropic Fractional Smoothing via Decaying-Adapted Techniques</td>
      <td>67</td>
    </tr>
    <tr>
      <td>13.8. Microlocal Smoothing under Degenerate Weights</td>
      <td>68</td>
    </tr>
    <tr>
      <td>14. Numerical Strategies for Solving the ODE System</td>
      <td>69</td>
    </tr>
    <tr>
      <td>14.1. Compactification and Domain Mapping</td>
      <td>70</td>
    </tr>
    <tr>
      <td>14.2. Shooting Method</td>
      <td>70</td>
    </tr>
    <tr>
      <td>14.3. Convergence Metrics and Residuals</td>
      <td>70</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>14.4. Simulated Profiles</td>
      <td>70</td>
    </tr>
    <tr>
      <td>15. On the Existence of Non-Trivial Self-Similar Solutions</td>
      <td>71</td>
    </tr>
    <tr>
      <td>15.1. Non-Trivial Profiles</td>
      <td>71</td>
    </tr>
    <tr>
      <td>15.2. Level 1. Explicit Numerical Construction</td>
      <td>71</td>
    </tr>
    <tr>
      <td>15.3. Level 2. Existence via Fixed-Point Scheme</td>
      <td>71</td>
    </tr>
    <tr>
      <td>15.3.1. Functional Analytic Setting</td>
      <td>72</td>
    </tr>
    <tr>
      <td>15.4. Level 3. Inner-Outer Expansion and Global Solvability via Matched Asymptotics</td>
      <td>73</td>
    </tr>
    <tr>
      <td>15.4.1. Addendum on Radial Resolvent and Residue Identities</td>
      <td>74</td>
    </tr>
    <tr>
      <td>15.4.2. Involvements and Limitations</td>
      <td>75</td>
    </tr>
    <tr>
      <td>16. Geometric Blow-Up Theory and Log-Lattice Fluid Singularities</td>
      <td>75</td>
    </tr>
    <tr>
      <td>16.1. Comparative Insights</td>
      <td>75</td>
    </tr>
    <tr>
      <td>16.2. A Mechanism for Unification</td>
      <td>76</td>
    </tr>
    <tr>
      <td>16.2.1. Gelfand-Type Blow-Up and Pole Coalescence Process</td>
      <td>76</td>
    </tr>
    <tr>
      <td>16.2.1.1. Gaussian-weighted Modal Amplification of the Radial Equation</td>
      <td>78</td>
    </tr>
    <tr>
      <td>16.2.1.2. Exact algebraic reduction and Hamiltonian structure</td>
      <td>79</td>
    </tr>
    <tr>
      <td>16.2.2. Abstract Arrangement of Blow-Up Manifolds</td>
      <td>83</td>
    </tr>
    <tr>
      <td>16.2.3. Geometric Patterns and Flow Dynamics</td>
      <td>85</td>
    </tr>
    <tr>
      <td>Bibliography · List of Works Mentioned in the Article</td>
      <td>87</td>
    </tr>
  </tbody>
</table>
# 1. Framework and Notational Setup

The compressible magnetohydrodynamic (MHD) equations describe the dynamics of a conducting viscous fluid under magnetic effects, coupling hydrodynamic motion and electromagnetic induction in a non-linear PDE setting. This Section presents the analytical and notational conventions adopted throughout, forming the basis for the subsequent discussion of self-similar formulations and regularity mechanisms.

Let $(\vec{u}, \vec{B}, \rho, p)$ denote the velocity field, magnetic field, mass density, and pressure, respectively. The relevant physical constants are:

<table><tr><td>$\vec{u}(\vec{x}, t)$</td><td>velocity field (m/s)</td></tr><tr><td>$\vec{B}(\vec{x}, t)$</td><td>magnetic field (Tesla, or Alfvén units)</td></tr><tr><td>$\rho(\vec{x}, t)$</td><td>fluid density (kg/m³)</td></tr><tr><td>$p$</td><td>fluid pressure (Pa)</td></tr><tr><td>$P_B = \frac{|\vec{B}|^2}{2\mu_0}$</td><td>magnetic pressure (Pa)</td></tr><tr><td>$\mu$</td><td>kinematic viscosity (m²/s)</td></tr><tr><td>$\nu$</td><td>magnetic diffusivity (m²/s)</td></tr><tr><td>$\mu_0$</td><td>magnetic permeability of vacuum (H/m)</td></tr></table>

# 2. From Leray to Modern Analysis: a Quick Look

Leray's pioneering architecture [48] introduced backward self-similar solutions of Navier-Stokes [57] [80] via

$$ \vec{u}(\vec{x}, t) = (-t)^{-\frac{1}{2}} \vec{U} \left( \frac{\vec{x}}{\sqrt{-t}} \right), \quad (1) $$

$$ p(\vec{x}, t) = (-t)^{-1} P \left( \frac{\vec{x}}{\sqrt{-t}} \right), \quad (2) $$

establishing the foundation for modern singularity analysis in non-linear fluid models and MHD.

Subsequent analyses in compressible MHD settled weak and strong existence results [38], [9], and local smoothness in the inviscid limit [82]. Regularized and filtered models [41] demonstrated global regularity, while decay properties of strong solutions [26] were characterized through spectral estimates. In [81] local strong solutions to the full 3D compressible MHD system are shown to extend globally in time, provided the density stays bounded away from vacuum and concentration. Structure-preserving discretizations [15], global decay inequalities [8], and non-resistive formulations with general boundary data [49] integrate seamlessly within the same analytical idealization, to which the variational mathematization in [21] adds a complementary perspective on weak-variational and boundary-driven dynamics.

In the incompressible limit, the non-existence of finite-energy self-similar solutions [58] and its refinements [83] delineate sharp constraints on admissible similarity profiles. Boundary-induced non-trivial states [35] and structural analyses [44] clarify how over-restrictive decay and compactness assumptions can suppress physically meaningful configurations.

Analogous patterns observed in simplified MHD models [67], [47] reinforce the plausibility of self-similar dynamics consistent with empirical scaling exponents and asymptotic regimes.

# 3. Full Compressible MHD Equations

The full system of equations in 3D reads

$$ \partial_t \rho + \nabla \cdot (\rho \vec{u}) = 0, \quad (3) $$

$$ \rho \left( \partial_t \vec{u} + (\vec{u} \cdot \nabla) \vec{u} \right) = -\nabla \left( p + \frac{|\vec{B}|^2}{2\mu_0} \right) + \frac{1}{\mu_0} (\vec{B} \cdot \nabla) \vec{B} + \mu \Delta \vec{u}, \quad (4) $$

$$ \partial_t \vec{B} = \nabla \times (\vec{u} \times \vec{B}) + \nu \Delta \vec{B}, \quad (5) $$

$$ \nabla \cdot \vec{B} = 0. \quad (6) $$
We interpret this system in the weak formulation on Sobolev spaces [79] $\vec{u}, \vec{B} \in W^s(\mathbb{R}^3)$, $\rho \in L^\infty \cap W^s$, $p \in W^{s-1}$, with $s \gg \frac{3}{2}$ (Theorem 6.1), where

$$W^s(\mathbb{R}^3) = \{\varphi \in S_c^*(\mathbb{R}^3) | (1 + |\xi|^2)^{\frac{3}{2}} \hat{\varphi}(\xi) \in L^2(\mathbb{R}^3)\}. \quad (7)$$

Here $S_c$ denotes the Schwartz space [32] [16] [77] of rapidly decreasing smooth functions, and $S_c^*$ its continuous dual, viz. the space of tempered distributions.

The system conserves the total mechanical and magnetic energy in the ideal limit $\mu = \nu = 0$ (inviscid dynamics), via the time-dependent functional,

$$\mathcal{E}[t] = \int_{\mathbb{R}^3} \left( \frac{1}{2} \rho |\vec{u}|^2 + \frac{1}{2\mu_0} |\vec{B}|^2 + \rho e(\rho) \right) dx, \quad (8)$$

which reflects the absence of viscous and resistive dissipation. The internal energy,

$$e'(\rho) = \frac{p(\rho)}{\rho^2}, \quad e(\rho) = e(1) + \int_1^\rho \frac{p(s)}{s^2} ds, \quad (9)$$

encodes the thermodynamic equation of state, and $\mathcal{E}[\cdot]$ plays a central role in both variational formulations and weak compactness arguments for the compressible MHD system.

Weak solutions $(\rho, \vec{u}, \vec{B})$ to the compressible MHD system are expected to satisfy the global in time energy inequality

$$\mathcal{E}[t_2] + \int_{t_1}^{t_2} \mathcal{D}[\rho, \vec{u}, \vec{B}](t) dt \le \mathcal{E}[t_1], \quad \text{for all } 0 \le t_1 \le t_2. \quad (10)$$

This inequality expresses the dissipative nature of the MHD dynamics in the non-ideal regime $\mu, \nu > 0$, and serves as a Lyapunov-type bound controlling the energy evolution over time. Eq. (10) is typically recovered by testing the weak formulation of the governing PDEs with the solution variables themselves, under suitable integrability and regularity assumptions (e.g. $\rho \vec{u} \in L_t^\infty L_x^2$, $\vec{B} \in L_t^\infty L_x^2$, $\nabla \vec{u}, \nabla \vec{B} \in L_{t,x}^2$).

A more advanced functional analytic framework involves working with the system in the form of a quasi-linear symmetric hyperbolic-parabolic system, under the Friedrichs-Lax-Kato theory [23] [46] [40], which guarantees local in time existence in Banach spaces with time weights.

## 4. Fully Non-local Compressible MHD Model

The local system (3) (4) (5) (6) admits an equivalent formulation in terms of singular integral operators, by convolving the equations against the Newton potential and by rewriting differential operators in integral form. This perspective unveils the long-range structure inherent in the classical partial differential formulation and provides a pathway towards fractional and non-local generalizations.

A classical step consists in expressing the weak formulation of the system against smooth, compactly supported test fields. For any smooth vector field $\vec{\psi}$ with compact support, the momentum balance can be written as

$$\int_0^T \int_{\mathbb{R}^3} \left\{ \rho \vec{u} \cdot \partial_t \vec{\psi} + \rho (\vec{u} \otimes \vec{u}) : \nabla \vec{\psi} + \left( p + \frac{|\vec{B}|^2}{2\mu_0} \right) \nabla \cdot \vec{\psi} - \frac{1}{\mu_0} (\vec{B} \otimes \vec{B}) : \nabla \vec{\psi} \right\} d\vec{x} dt = \mu \int_0^T \int_{\mathbb{R}^3} \nabla \vec{u} : \nabla \vec{\psi} d\vec{x} dt, \quad (11)$$

and analogous distributional identities hold for the induction equation when tested against divergence-free vector fields.

The Newton potential $G(\vec{x}) = \frac{1}{4\pi|\vec{x}|}$ satisfies $-\Delta G = \delta$, so that for any sufficiently regular field $\varphi$ we have

$$\varphi(\vec{x}) = \int_{\mathbb{R}^3} \nabla G(\vec{x} - \vec{y}) \cdot \nabla \varphi(\vec{y}) d\vec{y}. \quad (12)$$

The differential operators appearing in the classical MHD equations can be written as convolution-type singular integral operators acting on the solution fields.
This integral suggests a systematic generalization: one can replace the Newton kernel with a fractional Riesz-type kernel of order $3 + \alpha$ and replace classical gradients with fractional difference quotients. In this sense, the compressible MHD equations in Section 3 arise as the limiting case $\alpha \to 0$ of a broader family of integro-differential systems with non-local transport, pressure, and magnetic tension.

Let $0 < \alpha < 2$ be fixed. For any locally integrable scalar or vector field $\varphi$, we define the principal-value fractional integral by p.v. $\int_{\mathbb{R}^3} \frac{\varphi(\vec{y})}{|\vec{x}-\vec{y}|^{3+\alpha}} d\vec{y} = \lim_{\varepsilon \to 0^+} \int_{|\vec{x}-\vec{y}|>\varepsilon} \frac{\varphi(\vec{y})}{|\vec{x}-\vec{y}|^{3+\alpha}} d\vec{y}$.

Let

$$K_H(\vec{x}, \vec{y}) = \frac{\vec{x}-\vec{y}}{|\vec{x}-\vec{y}|^{3+\alpha}} \text{ be the vector kernel,}$$

$$K(\vec{x}, \vec{y}, \vec{z}) = \frac{(\vec{x}-\vec{y})\otimes(\vec{x}-\vec{z})}{|\vec{x}-\vec{y}|^{3+\alpha}|\vec{x}-\vec{z}|^{3+\alpha}} \text{ the symmetric double kernel, and}$$

$$\Phi(\vec{x}, \vec{y}) = \frac{1}{2} \frac{(\vec{x}-\vec{y}) \cdot (\vec{u}(\vec{x},t) + \vec{u}(\vec{y},t))}{|\vec{x}-\vec{y}|^{3+\alpha}} \text{ the scalar exchange kernel (symmetric in } x, y).$$

We shall denote the Newton potential by $G(\vec{x}) = -\frac{1}{4\pi|\vec{x}|}$, $\nabla G(\vec{x}) = \frac{\vec{x}}{4\pi|\vec{x}|^3}$, and the Helmholtz-Leray projector by $\mathbb{P}_{HL}$.

### A. Continuity equation with non-local mass exchange.

$$\begin{aligned} \partial_t \rho(\vec{x}, t) + \nabla \cdot (\rho \vec{u})(\vec{x}, t) = & -\text{p.v.} \int_{\mathbb{R}^3} \Phi(\vec{x}, \vec{y}) (\rho(\vec{y}, t) - \rho(\vec{x}, t)) d\vec{y} \quad (13) \\ & - \text{p.v.} \int_{\mathbb{R}^3} \frac{\vec{x} - \vec{y}}{|\vec{x} - \vec{y}|^{3+\alpha}} \cdot (\rho \vec{u})(\vec{y}, t) d\vec{y}. \end{aligned}$$

The first non-local integral models fractional, long-range mass exchange controlled by the local velocity field. The second term is for a fractional divergence-type transport coupling distant cells.

### B. Integro-differential momentum equation. The integro-differential equality of (4) will be

$$\begin{aligned} \rho(\vec{x}, t) (\partial_t \vec{u} + (\vec{u} \cdot \nabla) \vec{u})(\vec{x}, t) = & -\nabla \left( p + \frac{|\vec{B}|^2}{2\mu_0} \right) (\vec{x}, t) + \frac{1}{\mu_0} (\vec{B} \cdot \nabla) \vec{B}(\vec{x}, t) + \mu \Delta \vec{u}(\vec{x}, t) \quad (14) \\ & + \text{p.v.} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} K(\vec{x}, \vec{y}, \vec{z}) [\vec{u}(\vec{y}, t) - \vec{u}(\vec{z}, t)] d\vec{y} d\vec{z} \\ & + \frac{1}{\mu_0} \text{p.v.} \int_{\mathbb{R}^3} \frac{(\vec{x} - \vec{y}) \times (\nabla_{\vec{y}} \times \vec{B}(\vec{y}, t))}{|\vec{x} - \vec{y}|^3} d\vec{y}. \end{aligned}$$

The double principal-value integral exemplifies a strongly non-local momentum transfer (fractional, anisotropic) and couples distant fluid elements in a singular manner. The last integral expresses a non-local magnetic-tension-like contribution written via the curl of $\vec{B}$ and the Newton-type kernel.

The magnetic tension can be simultaneously written as a non-local convolution with $G$,

$$\frac{1}{\mu_0} (\vec{B} \cdot \nabla) \vec{B}(\vec{x}, t) = \frac{1}{\mu_0} \text{p.v.} \int_{\mathbb{R}^3} \nabla_{\vec{x}} G(\vec{x} - \vec{y}) \cdot (\vec{B}(\vec{x}, t) \cdot \vec{B}(\vec{y}, t)) d\vec{y}, \quad (15)$$

whilst the global energy balance with fractional dissipation and non-local transfer will have this form,

$$\begin{aligned} \frac{d}{dt} \int_{\mathbb{R}^3} \left( \frac{1}{2} \rho |\vec{u}|^2 + \frac{1}{2\mu_0} |\vec{B}|^2 + \rho e(\rho) \right) d\vec{x} = & -\mu \int_{\mathbb{R}^3} |\nabla \vec{u}|^2 d\vec{x} - \nu \int_{\mathbb{R}^3} |\nabla \vec{B}|^2 d\vec{x} \\ & - \frac{1}{2} \iint_{\mathbb{R}^3 \times \mathbb{R}^3} \frac{|\vec{u}(\vec{x}, t) - \vec{u}(\vec{y}, t)|^2}{|\vec{x} - \vec{y}|^{3+\alpha}} d\vec{x} d\vec{y} \\ & - \frac{1}{2\mu_0} \iint_{\mathbb{R}^3 \times \mathbb{R}^3} \frac{|\vec{B}(\vec{x}, t) - \vec{B}(\vec{y}, t)|^2}{|\vec{x} - \vec{y}|^{3+\alpha}} d\vec{x} d\vec{y}, \end{aligned} \quad (16)$$

where the fractional dissipation constants are positive, so that $\frac{d}{dt} \mathcal{E}[t] \le 0$.
### C. Induction with non-local and convolutional reconnection terms.

$$ 
\begin{aligned}
\partial_t \vec{B}(\vec{x}, t) = & \nabla \times (\vec{u} \times \vec{B})(\vec{x}, t) + \nu \triangle \vec{B}(\vec{x}, t) \\
& - \nabla \times \text{p.v.} \int_{\mathbb{R}^3} K_H(\vec{x}, \vec{y}) \times (\nabla_{\vec{y}} \times \vec{B}(\vec{y}, t)) d\vec{y} \\
& + \nabla \times \left( \text{p.v.} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} K(\vec{x}, \vec{y}, \vec{z}) [\vec{B}(\vec{y}, t) \times \vec{u}(\vec{z}, t)] d\vec{y} d\vec{z} \right).
\end{aligned} 
\quad (17) $$

The second line adds a non-local Hall-like induction term modeled by a principal-value convolution of the curl of $\vec{B}$ with a singular vector kernel. The double convolution illustrates non-local transfer of magnetic topology mediated by distant velocity and magnetic structures.

### D. Constraint and its non-local enforcement via the Helmholtz-Leray projector. Let

$$ \mathbb{P}_{\text{HL}} \vec{F}(\vec{x}) = \vec{F}(\vec{x}) + \nabla \int_{\mathbb{R}^3} G(\vec{x} - \vec{y}) \nabla_{\vec{y}} \cdot \vec{F}(\vec{y}) d\vec{y}, \quad (18) $$

which guarantees $\nabla \cdot (\mathbb{P}_{\text{HL}} \vec{F}) = 0$ for sufficiently regular $\vec{F}$. So

$$ \nabla \cdot \vec{B}(\vec{x}, t) = 0, \quad (19) $$

$$ \vec{B}(\vec{x}, t) = \mathbb{P}_{\text{HL}} \vec{B}(\vec{x}, t). \quad (20) $$

Eq. (20) enforces the constraint by projecting any candidate magnetic field onto the divergence-free subspace through a non-local Newton-potential.

## 5. A Variational Methodology

The time-dependent energy functional $\mathcal{E}[t]$ naturally induces a weak variational structure for the compressible MHD system. To make this precise, we introduce this functional space,

$$ V_{\mathcal{E}} = \left\{ (\vec{U}, \vec{B}, R) \begin{array}{l} \left| \begin{array}{l} \vec{U} \in W_{\text{rad}}^{2}(\mathbb{R}^{3}; \mathbb{R}^{3}) \\ \vec{B} \in W_{\text{rad}}^{2}(\mathbb{R}^{3}; \mathbb{R}^{3}) \\ R \in W_{\text{rad}}^{1}(\mathbb{R}^{3}; \mathbb{R}) \end{array} \right. \end{array} \right\}, \quad (21) $$

with radial symmetry enforced in the Sobolev norm,

$$ \| (R, \vec{U}, \vec{B}) \|_{V_{\mathcal{E}}}^{2} = \int_{\mathbb{R}^3} \left[ R^2 + |\vec{U}|^2 + |\vec{B}|^2 + |\nabla \vec{U}|^2 + |\nabla \vec{B}|^2 + |\nabla R|^2 \right] w(x) dx, \quad (22) $$

where $w(x)$ is a radial weight function such as $w(x) = (1 + |x|^2)^{-\alpha}$, $\alpha > 2$, to ensure integrability at infinity. Here $W_{\text{rad}}^k(\mathbb{R}^3)$ denotes the subspace of the Sobolev space $W^k(\mathbb{R}^3)$ consisting of radially symmetric functions, that is, functions depending only on $r = |\vec{x}|$. The radial symmetry reduces the complexity of the variational formulation to scalar equations in one dimension.

Each component is equipped with a weighted Sobolev norm

$$ \|\varphi\|_{W_{\text{rad}}^k}^2 = \sum_{j=0}^{k} \int_{0}^{\infty} |r^j \varphi^{(j)}(r)|^2 r^{2-2j} dr, \quad (23) $$

which balances smoothness (derivatives up to order $k$) and decay (via radial weights).

This choice of space underpins

(1) sufficient regularity for performing integration by parts in the weak formulation,

(2) appropriate decay of profiles as $r \to \infty$,

(3) compatibility with the scaling induced by the dilation group $G_\lambda$ (Section 11.1),

(4) well-posedness of the weak formulation with test functions in $V_{\mathcal{E}}$.

Moreover, $\mathcal{E}[t]$ provides a natural functional framework to carry out energy-based arguments, compactness and variational minimization techniques, and connects directly with the energy functional $\mathcal{E}[t]$, which acts as a Lyapunov functional (Section 12) for the inviscid MHD system.

We derive the weak formulation by testing each equation of the self-similar compressible MHD system against smooth, compactly supported, radially symmetric test functions $(\vec{\phi}, \vec{\theta}, \chi) \in V_{\mathcal{E}}$.
Integrating over $\mathbb{R}^3$ and performing integration by parts (justified by the decay provided by the weighted norm), one obtains coupled bilinear and trilinear forms,

$$a[(\vec{U}, \vec{B}, R), (\vec{\phi}, \vec{\theta}, \chi)] + b[(\vec{U}, \vec{B}, R), (\vec{\phi}, \vec{\theta}, \chi)] = 0, \quad (24)$$

where $a[\cdot, \cdot]$ is bilinear and controls the dissipative and linear scaling terms, while $b[\cdot, \cdot]$ contains the non-linear advection and magnetic coupling terms. This structure is amenable to fixed-point arguments (Section 15.3.1).

We now write out the weak formulation of the reduced, self-similar compressible MHD equations. Testing the mass, momentum, and magnetic equations with $(\chi, \vec{\phi}, \vec{\theta}) \in V_g$ and integrating by parts yields the following system:

(1) mass conservation (tested with $\chi$),

$$\int_{\mathbb{R}^3} [-\gamma R\chi - \beta(\vec{X} \cdot \nabla R)\chi + (\vec{U} \cdot \nabla R)\chi + R(\nabla \cdot \vec{U})\chi] dx = 0, \quad (25)$$

(2) momentum equation (tested with $\vec{\phi}$),

$$\int_{\mathbb{R}^3} R \left( -\alpha \vec{U} - \beta(\vec{X} \cdot \nabla)\vec{U} + (\vec{U} \cdot \nabla)\vec{U} \right) \cdot \vec{\phi} dx + \mu \int_{\mathbb{R}^3} \nabla \vec{U} : \nabla \vec{\phi} dx \\ = - \int_{\mathbb{R}^3} P(R) \nabla \cdot \vec{\phi} dx - \frac{1}{\mu_0} \int_{\mathbb{R}^3} (\vec{B} \cdot \nabla) \vec{B} \cdot \vec{\phi} dx, \quad (26)$$

(3) magnetic field equation (tested with $\vec{\theta}$),

$$\int_{\mathbb{R}^3} \left\{ -\alpha \vec{B} - \beta(\vec{X} \cdot \nabla)\vec{B} + (\vec{U} \cdot \nabla)\vec{B} \right\} \cdot \vec{\theta} dx \\ = \nu \left\{ \int_{\mathbb{R}^3} \nabla \vec{B} : \nabla \vec{\theta} dx + \int_{\mathbb{R}^3} \left[ (\vec{B} \cdot \nabla)\vec{U} - (\nabla \cdot \vec{U})\vec{B} \right] \cdot \vec{\theta} \right\} dx. \quad (27)$$

This variational structure allows for a precise definition of weak solutions in the energy space $V_g$, where non-linearity and degeneracy are handled by the trilinear structure and radial symmetry.

The coercivity of the bilinear form $a[\cdot, \cdot]$ under the weighted norm, together with appropriate bounds on $b[\cdot, \cdot]$, enables the use of compactness methods, via Lions-Aubin lemma [4] [51] in weighted Sobolev spaces, and monotonicity techniques. The weak formulation is also the natural setting for constructing approximate solutions via Galerkin approximation [25] and for proving convergence to energy-dissipating self-similar profiles.

## 5.1. Recasting in Weighted Lebesgue Spaces

Alternatively, one can recast $V_g$ in terms of weighted Lebesgue spaces with radial measure

$$L_{\delta}^{2}(\mathbb{R}^{3}) = \left\{ \varphi : \mathbb{R}^{3} \rightarrow \mathbb{R} \left| \int_{\mathbb{R}^{3}} |\varphi(x)|^{2}(1 + |x|^{2})^{\delta} dx < \infty \right. \right\}, \quad (28)$$

with typical choices of $\delta < -\frac{3}{2}$ to ensure decay and integrability at infinity.

In this framework, one specifies

$$V_g^w = L_{\delta_1}^2(\mathbb{R}^3; \mathbb{R}^3) \times L_{\delta_2}^2(\mathbb{R}^3; \mathbb{R}^3) \times L_{\delta_3}^2(\mathbb{R}^3), \quad (29)$$

with appropriate exponents $(\delta_1, \delta_2, \delta_3)$ controlling decay of velocity, magnetic field and density, respectively.

Such weighted Lebesgue spaces allow energy inequalities of the form

$$\frac{d}{dt} \mathcal{E}[t] \leq -\mathcal{D}[t], \quad (30)$$

where the dissipation functional

$$\mathcal{D}[t] = \mathcal{D}_{\mu}[t] + \mathcal{D}_{\nu}[t] \quad (31)$$
accounts for viscous dissipation $\mathcal{D}_{\mu}[t]$ and magnetic diffusion dissipation $\mathcal{D}_{\nu}[t]$, specified respectively as

$$ \mathcal{D}_{\mu}[t] = \mu \int_{\mathbb{R}^3} |\nabla \vec{U}(x, t)|^2 (1 + |x|^2)^{\delta_1} dx, \quad (32) $$

$$ \mathcal{D}_{\nu}[t] = \nu \int_{\mathbb{R}^3} |\nabla \vec{B}(x, t)|^2 (1 + |x|^2)^{\delta_2} dx. \quad (33) $$

These functionals capture the energy decay mechanisms and are essential in weak convergence and compactness theory, especially when dealing with vanishing dissipation limits.

## 6. Criteria for Finite-Time Singularities

We begin by formulating a localized criterion for finite-time singularity formation in the compressible magnetohydrodynamic flow, based on weighted energy concentration and Sobolev norm amplification.

**Theorem 6.1** (Local blow-up characterization for compressible MHD). Let

$$ \vec{U}_{\text{phys}}(\vec{x}, t) = (\rho, \vec{u}, \vec{B}, p)^{\text{T}} \quad (34) $$

be a classical solution to the full compressible magnetohydrodynamics system in dimension $d = 3$, posed on $\mathbb{R}^3 \times [0,T)$ with $T < \infty$, where the superscript T is the matrix transposition. The vector $\vec{B}$ denotes the magnetic flux density, related to the magnetic field intensity $\vec{H}$ via the constitutive relation

$$ \vec{B} = \mu \vec{H}, \quad (35) $$

where $\mu > 0$ is the magnetic permeability of the medium. In this formulation, $\vec{H}$ is implicitly defined and does not explicitly appear as an independent unknown in the system.

Assume the initial data obeys

$$ \rho_0 \in L^\infty \cap \mathfrak{F}^s(\mathbb{R}^3), \quad \vec{u}_0, \vec{B}_0 \in \mathfrak{F}^s(\mathbb{R}^3), \quad p_0 \in \mathfrak{F}^{s-1}(\mathbb{R}^3), \quad \text{for some } s > \frac{d}{2} + 1. \quad (36) $$

If the solution develops a singularity at $t = T$, then

$$ \int_0^T \left( \|\nabla \vec{u}(\cdot, t)\|_{L^\infty} + \|\nabla \vec{B}(\cdot, t)\|_{L^\infty} \right) dt = \infty. \quad (37) $$

That is, finite-time blow-up implies divergence in the Lipschitz seminorm of either the velocity field $\vec{u}$ or the magnetic field $\vec{B}$.

**Proof.** Let $s > \frac{d}{2} + 1$ so that $\mathfrak{F}^s(\mathbb{R}^3)$ embeds continuously into the Lipschitz class $W^{1,\infty}(\mathbb{R}^3)$. The compressible MHD system is locally well-posed in $\mathfrak{F}^s$ for such $s$, and the existence time depends on the norm

$$ \|(\rho, \vec{u}, \vec{B}, p)(t)\|_{\mathfrak{F}^s \times \mathfrak{F}^s \times \mathfrak{F}^s \times \mathfrak{F}^{s-1}}. $$

Suppose the integral

$$ \int_0^T \left( \|\nabla \vec{u}(t)\|_{L^\infty} + \|\nabla \vec{B}(t)\|_{L^\infty} \right) dt < \infty. \quad (38) $$

Then, by standard commutator and energy estimates (see Kato-Ponce-type inequalities), one has that the $\mathfrak{F}^s$ norms remain finite on $[0,T)$. The energy functional

$$ \mathcal{E}_s[t] = \|\rho(t)\|_{\mathfrak{F}^s}^2 + \|\vec{u}(t)\|_{\mathfrak{F}^s}^2 + \|\vec{B}(t)\|_{\mathfrak{F}^s}^2 + \|p(t)\|_{\mathfrak{F}^{s-1}}^2 \quad (39) $$

satisfies a differential inequality

$$ \frac{d}{dt} \mathcal{E}_s[t] \le \left( \|\nabla \vec{u}(t)\|_{L^\infty} + \|\nabla \vec{B}(t)\|_{L^\infty} \right) \mathcal{E}_s[t]. \quad (40) $$

Applying Grönwall's lemma [31], we deduce

$$ \mathcal{E}_s[t] \le \mathcal{E}_s[0] \exp \left\{ \int_0^t \left( \|\nabla \vec{u}(\tau)\|_{L^\infty} + \|\nabla \vec{B}(\tau)\|_{L^\infty} \right) d\tau \right\}. \quad (41) $$
If

$$ \int_{0}^{T} \left( \|\nabla \vec{u}(\tau)\|_{L^{\infty}} + \|\nabla \vec{B}(\tau)\|_{L^{\infty}} \right) d\tau < \infty, \qquad (42) $$

then $\mathcal{E}_s[t]$ remains bounded on $[0, T)$, which permits continuation of the solution beyond $T$. This contradicts the blow-up assumption, hence the integral must diverge as $t \nearrow T$, i.e.,

$$ \int_{0}^{T} \left( \|\nabla \vec{u}(t)\|_{L^{\infty}} + \|\nabla \vec{B}(t)\|_{L^{\infty}} \right) dt = \infty. \qquad (43) $$

This is an MHD analogue of the Beale-Kato-Majda blow-up condition [7] for the incompressible Euler system, where the vorticity $\omega = \nabla \times \vec{u}$ meets

$$ \int_{0}^{T} \|\omega(t)\|_{L^{\infty}} dt = \infty \Rightarrow \text{loss of regularity at } t = T. \qquad (44) $$

In the compressible case, $\omega$ no longer controls the full dynamics due to compressibility and coupling with $\vec{B}$. Nonetheless, the critical mechanism—blow-up triggered by growth in the Lipschitz norm—remains valid. The extension of the BKM framework to MHD systems is justified by similar structural properties of the non-linearities, including the transport-type structure and boundedness of commutators in high regularity spaces. $\square$

Scholium 6.1 (Scaling-critical energy control in $\mathcal{F}^s$). To make the argument quantitative, observe that the energy functional $\mathcal{E}_s[t]$ satisfies the differential inequality

$$ \frac{d}{dt} \mathcal{E}_s[t] \le C_s \left( \|\nabla \vec{u}(t)\|_{L^\infty} + \|\nabla \vec{B}(t)\|_{L^\infty} \right) \mathcal{E}_s[t], \qquad (45) $$

where $C_s > 0$ depends on $s$ through the commutator constants in the estimates

$$ \|[\Lambda^s, \varphi]g\|_{L^2} \lesssim \|\nabla\varphi\|_{L^\infty} \|g\|_{\mathcal{F}^{s-1}}, \quad s > 1, \qquad (46) $$

letting $\Lambda^s = (1 - \Delta)^{\frac{s}{2}}$ be the standard pseudo-differential operator via the Fourier transform $F$ as

$$ \widehat{\Lambda^s \varphi}(\xi) = (1 + |\xi|^2)^{\frac{s}{2}} \hat{\varphi}(\xi), \qquad (47) $$

where

$$ \hat{\varphi}(\xi) = F[\varphi](\xi) = \int_{\mathbb{R}^3} \varphi(x) e^{-ix \cdot \xi} dx. \qquad (48) $$

This allows $\Lambda^s$ to be viewed as a Fourier multiplier operator with symbol $m_s(\xi) = (1 + |\xi|^2)^{\frac{s}{2}}$, a smooth and positive function that grows polynomially at infinity, thus defining an isomorphism between the Sobolev space $\mathcal{F}^s(\mathbb{R}^3)$ and $L^2(\mathbb{R}^3)$,

$$ \|\varphi\|_{\mathcal{F}^s} = \|\Lambda^s \varphi\|_{L^2} = \left( \int_{\mathbb{R}^3} (1 + |\xi|^2)^s |\hat{\varphi}(\xi)|^2 d\xi \right)^{\frac{1}{2}}. \qquad (49) $$

The commutator estimate relies on the symbolic calculus of pseudo-differential operators and the boundedness of multiplier operators on Sobolev spaces, which guarantees the control of non-linearities arising in the MHD system at this level of regularity.

When $\mathcal{E}_s[0]$ is finite and the integral of $\|\nabla \vec{u}\|_{L^\infty} + \|\nabla \vec{B}\|_{L^\infty}$ is bounded over $[0,T)$, we can infer from Grönwall's lemma that

$$ \mathcal{E}_s[t] \le \mathcal{E}_s[0] \exp \left( C_s \int_0^t \left( \|\nabla \vec{u}(\tau)\|_{L^\infty} + \|\nabla \vec{B}(\tau)\|_{L^\infty} \right) d\tau \right), \qquad (50) $$

for all $t \in [0, T)$.

Moreover, this estimate is sharp in the scaling sense: the threshold $s > \frac{d}{2} + 1$ is critical to ensure that the Sobolev algebra structure $\mathcal{F}^s \cdot \mathcal{F}^s \subset \mathcal{F}^s$ holds, and that non-linear terms in the MHD system preserve regularity.

The blow-up mechanism is controlled by the growth of the Lipschitz norms of $\vec{u}$ and $\vec{B}$, consistent with the quasi-linear nature of the system and the paradifferential structure of the equations. $\diamond$
The preceding blow-up criterion admits a complementary continuation result: smooth solutions persist as long as the integral of the Lipschitz norms remains finite. We formulate this as an analytic well-posedness theorem for the full viscous-resistive MHD system.

**Theorem 6.2** (Analytic well-posedness and continuation for the full MHD system). Let $d = 3$ and consider the compressible magnetohydrodynamic Eqq. (3) (4) (5) (6) with viscosity $\mu > 0$ and resistivity $\nu > 0$. Fix initial data,

$$ \rho_0 \in L^\infty(\mathbb{R}^3) \cap \mathfrak{F}^s(\mathbb{R}^3), \quad (51) $$

$$ \vec{u}_0, \vec{B}_0 \in \mathfrak{F}^s(\mathbb{R}^3), \quad p_0 \in \mathfrak{F}^{s-1}(\mathbb{R}^3), \quad (52) $$

with $s > \frac{d}{2} + 1$, $\rho_0 > 0$, and $\nabla \cdot \vec{B}_0 = 0$. Then there exists $T > 0$ such that

$$ (\rho, \vec{u}, \vec{B}, p) \in C([0,T]; (L^\infty \cap \mathfrak{F}^s) \times \mathfrak{F}^s \times \mathfrak{F}^{s-1}) \quad (53) $$

has a unique classical solution on $[0,T]$. If $T^* < \infty$ is the maximal time of existence, then necessarily

$$ \int_0^{T^*} \left( \|\nabla \vec{u}(t)\|_{L^\infty} + \|\nabla \vec{B}(t)\|_{L^\infty} \right) dt = \infty. \quad (54) $$

*Proof.* First, $\triangle$ acts componentwise on vector fields,

$$ \triangle \vec{v} = (\partial_{x_1 x_1} v_1 + \partial_{x_2 x_2} v_1 + \partial_{x_3 x_3} v_1, \partial_{x_1 x_1} v_2 + \partial_{x_2 x_2} v_2 + \partial_{x_3 x_3} v_2, \partial_{x_1 x_1} v_3 + \partial_{x_2 x_2} v_3 + \partial_{x_3 x_3} v_3)^T. \quad (55) $$

Accordingly the viscous and resistive terms are $\mu \triangle \vec{u}$ and $\nu \triangle \vec{B}$ (applied componentwise). All spatial integrals below are over $\mathbb{R}^3$; decay at infinity justifies integration by parts.

Let $\partial^\alpha$ be any spatial derivative with $|\alpha| \le s$. Differentiate the momentum and induction equations by $\partial^\alpha$, multiply the differentiated momentum equation by $\partial^\alpha \vec{u}$ and the differentiated induction equation by $\partial^\alpha \vec{B}$, integrate and sum over $|\alpha| \le s$. Using integration by parts and $\nabla \cdot \vec{B} = 0$ we have the exact identity,

$$ \frac{1}{2} \frac{d}{dt} (\|\vec{u}\|_{\mathfrak{F}^s}^2 + \|\vec{B}\|_{\mathfrak{F}^s}^2) + \mu \sum_{|\alpha| \le s} \langle \triangle \partial^\alpha \vec{u}, \partial^\alpha \vec{u} \rangle_{L^2} + \nu \sum_{|\alpha| \le s} \langle \triangle \partial^\alpha \vec{B}, \partial^\alpha \vec{B} \rangle_{L^2} = \Sigma, \quad (56) $$

where $\Sigma$ collects the transport and commutator contributions. Integration by parts recovers the coercive form of the dissipative terms (boundary terms vanish by decay),

$$ \langle \triangle \partial^\alpha \vec{u}, \partial^\alpha \vec{u} \rangle_{L^2} = - \|\nabla \partial^\alpha \vec{u}\|_{L^2}^2, \quad \langle \triangle \partial^\alpha \vec{B}, \partial^\alpha \vec{B} \rangle_{L^2} = - \|\nabla \partial^\alpha \vec{B}\|_{L^2}^2. \quad (57) $$

From here the dissipative contribution is

$$ \mu \|\nabla \vec{u}\|_{\mathfrak{F}^s}^2 + \nu \|\nabla \vec{B}\|_{\mathfrak{F}^s}^2. \quad (58) $$

We now estimate $\Sigma$. Decompose each transport term by commutator,

$$ \partial^\alpha (\vec{u} \cdot \nabla \vec{u}) = \vec{u} \cdot \nabla \partial^\alpha \vec{u} + [\partial^\alpha, \vec{u} \cdot \nabla] \vec{u}, \quad (59) $$

and analogous decompositions for the magnetic terms. Integrating the principal transport piece $\vec{u} \cdot \nabla \partial^\alpha \vec{u}$ against $\partial^\alpha \vec{u}$ produces a skew-symmetric term that vanishes under integration, up to terms controlled by $\|\nabla \vec{u}\|_{L^\infty}$-norm. The remaining commutators are estimated by the sharp Moser-Kato bounds [55] [56] [39]. For $s > \frac{d}{2}$,

$$ \| [\partial^\alpha, \varphi] g \|_{L^2} \le C (\|\nabla \varphi\|_{L^\infty} \|g\|_{\mathfrak{F}^{s-1}} + \|\varphi\|_{\mathfrak{F}^s} \|\nabla g\|_{L^\infty}). \quad (60) $$

Applying this with $\varphi = \vec{u}$ or $\varphi = \vec{B}$ and summing over $|\alpha| \le s$ implies the concrete bound,

$$ |\Sigma| \le C (\|\nabla \vec{u}\|_{L^\infty} + \|\nabla \vec{B}\|_{L^\infty}) (\|\vec{u}\|_{\mathfrak{F}^s}^2 + \|\vec{B}\|_{\mathfrak{F}^s}^2). \quad (61) $$

For the density, commute $\partial^\alpha$ with $\partial_t \rho + \vec{u} \cdot \nabla \rho = -\rho \nabla \cdot \vec{u}$ and take $L^2$-inner products to derive

$$ \frac{d}{dt} \|\rho\|_{\mathfrak{F}^s}^2 \le C \|\nabla \vec{u}\|_{L^\infty} \|\rho\|_{\mathfrak{F}^s}^2. \quad (62) $$
Collecting the two differential inequalities we arrive at

$$ \frac{d}{dt} (\|\rho\|_{\mathfrak{F}^s}^2 + \|\vec{u}\|_{\mathfrak{F}^s}^2 + \|\vec{B}\|_{\mathfrak{F}^s}^2) \le C (\|\nabla \vec{u}\|_{L^\infty} + \|\nabla \vec{B}\|_{L^\infty}) (\|\rho\|_{\mathfrak{F}^s}^2 + \|\vec{u}\|_{\mathfrak{F}^s}^2 + \|\vec{B}\|_{\mathfrak{F}^s}^2). \quad (63) $$

To justify the linear compensation, note that the linear operators,

$$ \vec{u} \mapsto \mu \Delta \vec{u}, \quad \vec{B} \mapsto \nu \Delta \vec{B}, \quad (64) $$

are sectorial on $\mathfrak{F}^s(\mathbb{R}^3)$ and generate analytic heat semigroups $e^{t\mu \Delta}$ and $e^{t\nu \Delta}$. On $\mathbb{R}^3$, these semigroups yield parabolic smoothing rather than exponential decay of the full energy: for each multi-index $|\alpha| \le s$ and $t > 0$,

$$ \|\nabla^k e^{t\mu \Delta} \partial^\alpha \varphi\|_{L^2} \le C_k t^{-\frac{k}{2}} \|\partial^\alpha \varphi\|_{L^2}, \quad (65) $$

$$ \|\nabla^k e^{t\nu \Delta} \partial^\alpha g\|_{L^2} \le C_k t^{-\frac{k}{2}} \|\partial^\alpha g\|_{L^2}, \quad (66) $$

for constants $C_k > 0$ depending only on $k, \mu, \nu$. Thus $e^{t\mu \Delta}$ and $e^{t\nu \Delta}$ are smoothing but not exponentially decaying in $t$ on the unbounded domain $\mathbb{R}^3$.

The exponential factor in Grönwall-type inequalities arises solely from non-linear estimates, not from the linear semigroup itself,

$$ \begin{aligned} \|\rho(t)\|_{\mathfrak{F}^s}^2 &+ \|\vec{u}(t)\|_{\mathfrak{F}^s}^2 + \|\vec{B}(t)\|_{\mathfrak{F}^s}^2 \\ &\le (\|\rho_0\|_{\mathfrak{F}^s}^2 + \|\vec{u}_0\|_{\mathfrak{F}^s}^2 + \|\vec{B}_0\|_{\mathfrak{F}^s}^2) \\ &\quad \times \exp \left( C \int_0^t (\|\nabla \vec{u}\|_{L^\infty} + \|\nabla \vec{B}\|_{L^\infty}) d\tau \right), \end{aligned} \quad (67) $$

for $t \in [0, T]$.

The $\mathfrak{F}^s$ norms remain finite while the integral in the exponential is finite, and the local existence scheme of the source (Picard iteration in the class determined by these bounds) extends as long as $\int_0^T (\|\nabla \vec{u}\|_{L^\infty} + \|\nabla \vec{B}\|_{L^\infty}) dt$ remains finite. The continuation criterion follows immediately. $\square$

## 7. Dissipative Regime and Asymptotic Stabilization

We inspect how smallness of the initial perturbation and the combined viscous-resistive mechanism enforce quantitative alignment between the velocity and magnetic fields, leading to asymptotic stabilization, algebraic on $\mathbb{R}^3$, exponential only under additional spectral-gap or compactness hypotheses.

**Theorem 7.1 (Magneto-kinematic alignment for dissipative MHD).** Let the hypotheses of Theorem 6.2 hold and assume in addition $\mu > 0, \nu > 0$. There exists $\varepsilon_1 = \varepsilon_1(s, \mu, \nu) > 0$ such that if the initial perturbation satisfies

$$ \|\rho_0 - 1\|_{\mathfrak{F}^s} + \|\vec{u}_0\|_{\mathfrak{F}^s} + \|\vec{B}_0\|_{\mathfrak{F}^s} \le \varepsilon_1, \quad (68) $$

then the solution is global and the velocity and magnetic field align asymptotically in the following quantified sense. For all times such that $\|\vec{B}(t)\|_{L^2} > 0$, define

$$ \vartheta(t) = \frac{\langle \vec{u}(t), \vec{B}(t) \rangle_{L^2}}{\|\vec{B}(t)\|_{L^2}^2}. \quad (69) $$

If $\vec{B}(t) = 0$ for some $t$, the alignment is trivially satisfied with $\vec{u}(t) = 0$, and by temporal continuity $\vartheta(t)$ is extended as the limit of the above expression when $\|\vec{B}(t)\|_{L^2} \to 0$.ª Then there exist constants $C > 0$ and $\lambda > 0$, depending only on $(s, \mu, \nu)$, such that for all $t \ge 0$

$$ \|\vec{u}(t) - \vartheta(t)\vec{B}(t)\|_{L^2} \le C\mathcal{E}_s[0]^{\frac{1}{2}}(1+t)^{-\lambda}. \quad (70) $$

The magneto-kinematic misalignment energy decays algebraically in time on $\mathbb{R}^3$, and exponentially fast on compact spatial subsets or under spectral localization of low frequencies.

<footnote>ª The definition is meaningful whenever $\|\vec{B}(t)\|_{L^2} > 0$. If $\vec{B}(t) = 0$ for some $t$, the alignment is trivially satisfied with $\vec{u}(t) = 0$, and $\vartheta(t)$ is extended by temporal continuity at such points.</footnote>
*Proof.* Under the smallness assumption on the initial data, global existence is ensured by the analytic well-posedness framework. Set $\vec{w}(t) = \vec{u}(t) - \vartheta(t)\vec{B}(t)$. The time derivative of $\vartheta(t)$ is taken only where $\|\vec{B}(t)\|_{L^2} > 0$, and extended continuously elsewhere. Differentiating $\vartheta(t)$ in time and substituting the MHD momentum and induction equations, integration by parts using $\nabla \cdot \vec{B} = 0$ gives

$$ \frac{1}{2} \frac{d}{dt} \|\vec{w}\|_{L^2}^2 + \mu \|\nabla \vec{w}\|_{L^2}^2 \le -\vartheta(t) \nu \langle \Delta \vec{B}, \vec{w} \rangle_{L^2} + \mathcal{Q}(t), \quad (71) $$

where $\mathcal{Q}(t)$ collects the quadratic and commutator contributions, estimated by the Moser-Kato bounds as

$$ |\mathcal{Q}(t)| \le C \left( \|\rho - 1\|_{\mathcal{F}^s} + \|\vec{u}\|_{\mathcal{F}^s} + \|\vec{B}\|_{\mathcal{F}^s} \right) \left( \|\nabla \vec{u}\|_{L^2}^2 + \|\nabla \vec{B}\|_{L^2}^2 \right). \quad (72) $$

For $\varepsilon_1$ sufficiently small, the non-linear contribution $\mathcal{Q}(t)$ can be absorbed by the dissipative terms at each fixed time. Using Cauchy-Schwarz and Young inequalities to treat the mixed term $-\vartheta\nu\langle\Delta\vec{B}, \vec{w}\rangle_{L^2}$ entails a differential inequality of the form

$$ \frac{d}{dt} \|\vec{w}\|_{L^2}^2 + 2\alpha(t) \|\vec{w}\|_{L^2}^2 \le 0, \quad (73) $$

where the time-dependent coefficient $\alpha(t) > 0$ arises from the dissipative terms and, by parabolic smoothing of higher norms, satisfies a lower bound of algebraic type for large $t$, for instance $\alpha(t) \gtrsim (1+t)^{-1}$ under the present small-data regime. Integration of the above inequality therefore implies an algebraic decay,

$$ \|\vec{w}(t)\|_{L^2} \le C \mathcal{E}_s[0]^{\frac{1}{2}} (1+t)^{-\gamma}, \quad (74) $$

for some exponent $\gamma > 0$ depending only on $(s, \mu, \nu)$. The magneto-kinematic misalignment decays algebraically in time on $\mathbb{R}^3$, while exponential alignment can be recovered under additional spectral-gap hypotheses or on compact spatial domains. $\square$

The exponential alignment estimate controls the relative geometry of $\vec{u}$ and $\vec{B}$, but does not yet quantify the regularizing influence of dissipation at higher derivatives. The next theorem identifies the precise parabolic smoothing effect inherent in the viscous-resistive MHD flow.

**Theorem 7.2 (Dissipative smoothing and parabolic gain).** Under the same hypotheses and smallness regime as in Theorem 7.1, the global solution satisfies for each integer $m \ge 0$ with $s + 2m$ finite,

$$ \int_{\frac{t}{2}}^{t} \left( \|\vec{u}(\tau)\|_{\mathcal{F}^{s+2m}}^2 + \|\vec{B}(\tau)\|_{\mathcal{F}^{s+2m}}^2 \right) d\tau \le C_m t^{-2m} \mathcal{E}_s[0], \quad (75) $$

and consequently,

$$ \|\vec{u}(t)\|_{\mathcal{F}^{s+2m}} + \|\vec{B}(t)\|_{\mathcal{F}^{s+2m}} \le C'_m t^{-m} \mathcal{E}_s[0]^{\frac{1}{2}}. \quad (76) $$

*Proof.* Apply $\partial^\alpha$ with $|\alpha| \le s + 2m$ to the MHD equations, multiply respectively by $\partial^\alpha \vec{u}$ and $\partial^\alpha \vec{B}$, integrate over $\mathbb{R}^3$, and sum. After integrating by parts the viscous and resistive terms, each involving $\Delta$, one obtains

$$ \frac{1}{2} \frac{d}{dt} \mathcal{E}_{s+2m}[t] + \mu \|\nabla \vec{u}\|_{\mathcal{F}^{s+2m}}^2 + \nu \|\nabla \vec{B}\|_{\mathcal{F}^{s+2m}}^2 = \mathcal{J}_{s+2m}(t), \quad (77) $$

where $\mathcal{J}_{s+2m}(t)$ denotes the commutator and transport coupling terms. Using Moser-Kato estimates valid for $s + 2m > \frac{d}{2}$,

$$ |\mathcal{J}_{s+2m}(t)| \le C_{\text{com}} \mathcal{E}_s[0] \mathcal{E}_{s+2m}[t], \quad (78) $$

where $C_{\text{com}} > 0$ is the commutator constant depending only on $s, m$. Introduce a small parameter $\eta > 0$ and write

$$ |\mathcal{J}_{s+2m}(t)| \le \eta \mathcal{E}_{s+2m}[t] + \frac{C_{\text{com}}^2}{4\eta} \mathcal{E}_s[0]^2. \quad (79) $$

Choose $\eta$ so that $\eta \le \frac{1}{2}c$ (where $c = \min\{\mu, \nu\}$ is the dissipation constant),

$$ \eta = \min \left\{ \frac{c}{2}, \frac{c}{2C_{\text{com}}} \right\}. \quad (80) $$
Then

$$|\mathcal{F}_{s+2m}(t)| \le \frac{c}{2}\mathcal{E}_{s+2m}[t] + C^*\mathcal{E}_s[0]^2, \qquad (81)$$

and substituting into the differential inequality, it is determined that

$$\frac{d}{dt}\mathcal{E}_{s+2m}[t] + \frac{c}{2}\mathcal{E}_{s+2m}[t] \le C^*\mathcal{E}_s[0]^2, \qquad (82)$$

which closes under the smallness of $\mathcal{E}_s[0]$. Integrating over $(\frac{t}{2}, t)$ the time-averaged bound is

$$\int_{\frac{t}{2}}^{t} \mathcal{E}_{s+2m}[\tau] d\tau \le C_m t^{-2m} \mathcal{E}_s[0], \qquad (83)$$

and interpolation between $\mathcal{F}^s$ and $\mathcal{F}^{s+2m}$ provides the instantaneous estimate,

$$\|\vec{u}(t)\|_{\mathcal{F}^{s+2m}} + \|\vec{B}(t)\|_{\mathcal{F}^{s+2m}} \le C'_m t^{-m} \mathcal{E}_s[0]^{\frac{1}{2}}. \qquad (84)$$

□

The smoothing effect induces a hierarchy of time-decaying Sobolev norms, ensuring enhanced spatial regularity as $t$ increases. Combining this with the magneto-kinematic alignment mechanism we can write a corollary on the complete asymptotic isotropization of the dissipative MHD dynamics.

**Corollary 7.1** (Asymptotic isotropization and magnetic relaxation). Under the hypotheses of Theorems 7.1 and 7.2, the global dissipative MHD flow satisfies complete asymptotic isotropization. Namely, there exist constants $C, \lambda > 0$, depending only on $(s, \mu, \nu)$, such that

$$\|\vec{u}(t)\|_{\mathcal{F}^s} + \|\vec{B}(t)\|_{\mathcal{F}^s} + \|\rho(t) - 1\|_{\mathcal{F}^s} \le C\mathcal{E}_s[0]^{\frac{1}{2}}(1+t)^{-\lambda}, \quad t \ge 0. \qquad (85)$$

In particular, the velocity and magnetic fields decay algebraically to zero in $\mathcal{F}^s(\mathbb{R}^3)$, while exponentially fast on compact subsets or under frequency localization,

$$\frac{\vec{u}(t)}{\|\vec{u}(t)\|_{L^2}} - \frac{\vec{B}(t)}{\|\vec{B}(t)\|_{L^2}} \to 0 \text{ in } L^2(\mathbb{R}^3) \text{ as } t \to \infty. \qquad (86)$$

Hence, dissipative magnetohydrodynamic flows relax algebraically (at the rate stated above) toward a spatially homogeneous, magnetically aligned equilibrium on $\mathbb{R}^3$; exponential relaxation is achieved only under additional spectral hypotheses or on compact domains.

*Proof.* From Theorem 7.1 we have exponential decay of the misalignment energy $\|\vec{u} - \vartheta\vec{B}\|_{L^2}^2$, and Theorem 7.2 confers uniform parabolic gain and algebraic damping of higher Sobolev energies. Summing these inequalities and invoking Grönwall's lemma (together with the algebraic decay of $\|\vec{u} - \vartheta\vec{B}\|_{L^2}^2$ and the parabolic gain for higher Sobolev norms) yields algebraic decay of $\mathcal{E}_s[t]$ as stated. Exponential decay would require additional spectral-gap assumptions. □

## 8. Liouville Constraints and Entropy-Based Estimates

Let us move on to a rigidity result showing that any finite-energy self-similar profile to the compressible MHD equations must be trivial.

**Theorem 8.1** (Liouville-type rigidity under finite energy). Let $(\vec{U}, \vec{H}, R)$ be a smooth solution to the stationary similarity system derived from the backward self-similar ansatz,

$$\vec{u}(\vec{x}, t) = \frac{1}{\sqrt{-t}} \vec{U} \left( \frac{\vec{x}}{\sqrt{-t}} \right), \quad \vec{B}(\vec{x}, t) = \frac{1}{\sqrt{-t}} \vec{H} \left( \frac{\vec{x}}{\sqrt{-t}} \right), \quad \rho(\vec{x}, t) = R \left( \frac{\vec{x}}{\sqrt{-t}} \right), \qquad (87)$$

with $R(\vec{X}) > 0$ and $(\vec{U}, \vec{H}, R) \in \mathcal{C}^2(\mathbb{R}^3)$ solving the reduced compressible MHD system,

$$-\frac{1}{2}\vec{U} - \frac{1}{2}(\vec{X} \cdot \nabla)\vec{U} + (\vec{U} \cdot \nabla)\vec{U} + \frac{1}{R}\nabla P(R) = \mu \Delta \vec{U} + \frac{1}{\mu_0 R}(\nabla \times \vec{H}) \times \vec{H}, \qquad (88)$$

$$-\frac{1}{2}\vec{H} - \frac{1}{2}(\vec{X} \cdot \nabla)\vec{H} + (\vec{U} \cdot \nabla)\vec{H} - (\vec{H} \cdot \nabla)\vec{U} = \nu \Delta \vec{H}, \qquad (89)$$

$$\nabla \cdot (R\vec{U}) = 0, \qquad (90)$$

where $P(R)$ is a smooth barotropic pressure law, $\frac{dP}{dR}(R) > 0$.
Let the total similarity energy be finite,

$$ \int_{\mathbb{R}^3} \left( \frac{1}{2} R |\vec{U}|^2 + \frac{|\vec{H}|^2}{2\mu_0} + \mathcal{P}(R \| R_\infty) \right) d\vec{X} < \infty, \qquad (91) $$

where

$$ \mathcal{P}(R \| R_\infty) = \int_{R_\infty}^{R} \frac{P(s) - P(R_\infty)}{s^2} ds \qquad (92) $$

is the relative internal energy with respect to a constant reference density $R_\infty > 0$.

Then

$$ \vec{U} = \vec{0}, \quad \vec{H} = \vec{0}, \quad R = R_\infty. \qquad (93) $$

**Scholium 8.1 (On the relative entropy).** The function $\mathcal{P}(R \| R_\infty)$ is the relative entropy (energy) functional that measures the deviation of a state $R$ from a constant reference density $R_\infty$. It arises naturally when comparing internal energies in compressible barotropic flows.

Formally, if $P(R)$ is a thermodynamic pressure law, then the internal energy per unit mass $e(R)$ fulfills

$$ \frac{d}{dR}(Re(R)) = \frac{P(R)}{R} \Rightarrow Re(R) = \int_{0}^{R} \frac{P(s)}{s} ds. \qquad (94) $$

The relative energy then admits a second-order Taylor expansion form,

$$ \mathcal{P}(R \| R_\infty) = e(R) - e(R_\infty) - e'(R_\infty)(R - R_\infty), \qquad (95) $$

and under standard regularity assumptions on $P$, one can show that this expression reduces to (92).

For a polytropic law $P(R) = aR^\gamma$ with $\gamma > 1$, we can compute

$$ \mathcal{P}(R \| R_\infty) = a \int_{R_\infty}^{R} \frac{s^\gamma - R_\infty^\gamma}{s^2} ds, \qquad (96) $$

which is positive for $R \neq R_\infty$, vanishes at $R = R_\infty$, and grows quadratically near $R_\infty$. ◇

**Proof.** Multiply the momentum equation (88) by $\vec{U}$ and integrate over $\mathbb{R}^3$. The linear scaling and transport terms satisfy

$$ \int_{\mathbb{R}^3} \left( -\frac{1}{2} \vec{U} \cdot \vec{U} - \frac{1}{2} (\vec{X} \cdot \nabla) \vec{U} \cdot \vec{U} \right) d\vec{X} = -\frac{1}{2} \int |\vec{U}|^2 d\vec{X} + \frac{3}{4} \int |\vec{U}|^2 d\vec{X} = \frac{1}{4} \int |\vec{U}|^2 d\vec{X}. \qquad (97) $$

The non-linear convective term vanishes,

$$ \int_{\mathbb{R}^3} (\vec{U} \cdot \nabla) \vec{U} \cdot \vec{U} = \frac{1}{2} \int_{\mathbb{R}^3} \vec{U} \cdot \nabla |\vec{U}|^2 = -\frac{1}{2} \int_{\mathbb{R}^3} (\nabla \cdot \vec{U}) |\vec{U}|^2, \qquad (98) $$

which is zero whenever $\nabla \cdot \vec{U} = 0$ or when $\vec{U}$ decays sufficiently fast at infinity.

Viscous and magnetic diffusion terms are non-negative,

$$ \mu \int_{\mathbb{R}^3} |\nabla \vec{U}|^2 + \nu \int_{\mathbb{R}^3} |\nabla \vec{H}|^2 \geq 0. \qquad (99) $$

For the pressure term,

$$ \int_{\mathbb{R}^3} \frac{1}{R} \nabla P(R) \cdot \vec{U} = - \int_{\mathbb{R}^3} \frac{P(R)}{R} (\nabla \cdot \vec{U}) = \int_{\mathbb{R}^3} \nabla \cdot \left( \frac{P(R)}{R} \vec{U} \right) - \int_{\mathbb{R}^3} \vec{U} \cdot \nabla \left( \frac{P(R)}{R} \right). \qquad (100) $$

The divergence term vanishes at infinity by integrability of the energy. By means of a smooth radial cutoff $\zeta_{\rho_{\text{rad}}}$ equal to 1 on $\bar{\mathbb{B}}_{\rho_{\text{rad}}}$ and 0 outside $\bar{\mathbb{B}}_{2\rho_{\text{rad}}}$, one arrives at

$$ \int_{\mathbb{R}^3} \frac{\nabla P(R)}{R} \cdot \vec{U} \zeta_{\rho_{\text{rad}}} = - \int_{\mathbb{R}^3} \left( P(R) - P(R_\infty) \right) \nabla \cdot \left( \frac{\vec{U}}{R} \zeta_{\rho_{\text{rad}}} \right). \qquad (101) $$

Passing to the limit $\rho_{\text{rad}} \to \infty$ and using decay implied by finite energy,

$$ \vec{U} \in L^2(\mathbb{R}^3), \quad P(R) - P(R_\infty) \in L^1(\mathbb{R}^3), $$
one recovers

$$ \int_{\mathbb{R}^3} \nabla P(R) \cdot \frac{\vec{U}}{R} = - \int_{\mathbb{R}^3} \left( P(R) - P(R_\infty) \right) \nabla \cdot \left( \frac{\vec{U}}{R} \right). \qquad (102) $$

The right-hand side vanishes because $\nabla \cdot (R\vec{U}) = 0$ in the weak sense.

For the induction equation, multiplying (89) by $\vec{H}$ and integrating by parts gives

$$ \frac{1}{4} \|\vec{H}\|_{L^2}^2 + \nu \|\nabla \vec{H}\|_{L^2}^2 = 0, \qquad (103) $$

hence $\vec{H} = 0$.

Returning to (88) and testing again with $\vec{U}$, ellipticity of the viscous term and the finite-energy condition imply

$$ \mu \|\nabla \vec{U}\|_{L^2}^2 = 0 \Rightarrow \vec{U} \text{ is constant.} \qquad (104) $$

Since $\vec{U} \in L^2(\mathbb{R}^3)$, it follows that $\vec{U} = 0$.

The continuity equation reduces to

$$ \nabla \cdot (R\vec{U}) = 0 \Rightarrow \vec{U} = 0 \Rightarrow \nabla R = 0, \qquad (105) $$

so $R$ is constant. The finite internal-energy constraint forces $R = R_\infty$.

By the weighted Poincaré inequality [65],

$$ \|R - R_\infty\|_{L^2_\omega} \lesssim \|\nabla R\|_{L^2_\omega}, \qquad (106) $$

with $\omega$ a radial weight adapted to the energy decay, it can be concluded that $R = R_\infty$.

Ergo $(\vec{U}, \vec{H}, R)$ coincides with the trivial constant profile, completing the proof. □

## 8.1. Some Details and Nuances

It is appropriate to take a closer look at some technical specificities that integrate further theoretical material to the topic seen in the previous Theorem 8.1.

### 8.1.1. Spectral Linearization Around the Momentum Equation

If we start from Eq. (88), and adopt the modal ansatz in similarity time as above, we recover a precise integral identity for the density-resolvent and the resulting non-local momentum operator.

Let

$$ R(\vec{X}, \tau) = R(\vec{X}) + e^{\lambda\tau} r(\vec{X}), \qquad (107) $$

$$ \vec{U}(\vec{X}, \tau) = \vec{U}(\vec{X}) + e^{\lambda\tau} \vec{u}(\vec{X}), \qquad (108) $$

$$ \vec{H}(\vec{X}, \tau) = \vec{H}(\vec{X}) + e^{\lambda\tau} \vec{h}(\vec{X}). \qquad (109) $$

The scalar transport-dilation resolvent is

$$ S_{f,\lambda}[u] = \lambda u + \frac{1}{2} \nabla \cdot (\vec{X} u) + \nabla \cdot (u \vec{U}), \qquad (110) $$

with a generic $u$-placeholder.

Linearized continuity becomes

$$ S_{f,\lambda}[r] = -\nabla \cdot (R\vec{u}). \qquad (111) $$

So, under spectral invertibility of $S_{f,\lambda}$,

$$ r = -S_{f,\lambda}^{-1} \left( \nabla \cdot (R\vec{u}) \right). \qquad (112) $$

From which the induction linearization,

$$ \begin{aligned} \lambda \vec{h} = & -\frac{1}{2} \vec{h} - \frac{1}{2} (\vec{X} \cdot \nabla) \vec{h} + (\vec{U} \cdot \nabla) \vec{h} + (\vec{u} \cdot \nabla) \vec{H} \\ & - (\vec{H} \cdot \nabla) \vec{u} - (\vec{h} \cdot \nabla) \vec{U} + \nu \triangle \vec{h}. \end{aligned} \qquad (113) $$
This allows us to outline the momentum linearization with explicit substitution $r = -S_{\mathcal{F}\lambda}^{-1}$,

$$ 
\begin{aligned}
\lambda \vec{u} = & -\frac{1}{2} \vec{u} - \frac{1}{2} (\vec{X} \cdot \nabla) \vec{u} + (\vec{U} \cdot \nabla) \vec{u} + (\vec{u} \cdot \nabla) \vec{U} \\
& + \underbrace{\frac{1}{R} \nabla (P'(R)r) - \frac{r}{R^2} \nabla P(R)}_{\text{linearized pressure with R-prefactor}} \\
& + \underbrace{\frac{1}{\mu_0 R} ((\nabla \times \vec{H}) \times \vec{h} + (\nabla \times \vec{h}) \times \vec{H}) - \frac{r}{\mu_0 R^2} ((\nabla \times \vec{H}) \times \vec{H})}_{\text{linearized Lorentz force}} \\
& + \mu \left( \frac{1}{R} \Delta \vec{u} - \frac{r}{R^2} \Delta \vec{U} \right).
\end{aligned} 
\quad (114) $$

Substitute $r$ from (112) into the momentum equation to obtain a non-local formula,

$$ 
\begin{aligned}
\lambda \vec{u} = & \left(-\frac{1}{2} - \frac{1}{2}(\vec{X} \cdot \nabla) + (\vec{U} \cdot \nabla)\right) \vec{u} + (\vec{u} \cdot \nabla) \vec{U} + \mu \frac{1}{R} \Delta \vec{u} \\
& - \frac{1}{R} \nabla \left(P'(R) S_{\mathcal{F}\lambda}^{-1} (\nabla \cdot (R \vec{u}))\right) + \frac{1}{\mu_0 R} \left((\nabla \times \vec{H}) \times \vec{h} + (\nabla \times \vec{h}) \times \vec{H}\right) \\
& + \frac{1}{\mu_0 R^2} \left(S_{\mathcal{F}\lambda}^{-1} (\nabla \cdot (R \vec{u}))\right) \left((\nabla \times \vec{H}) \times \vec{H}\right) - \mu \frac{S_{\mathcal{F}\lambda}^{-1} (\nabla \cdot (R \vec{u}))}{R^2} \Delta \vec{U}.
\end{aligned} 
\quad (115) $$

And here we go to the coupled spectral system,

$$ 
\left\{
\begin{array}{ll}
(i) & S_{\mathcal{F}\lambda}[r] + \nabla \cdot (R\vec{u}) = 0, \\
(ii) & \lambda\vec{h} + \frac{1}{2}\vec{h} + \frac{1}{2}(\vec{X} \cdot \nabla)\vec{h} - (\vec{U} \cdot \nabla)\vec{h} + (\vec{H} \cdot \nabla)\vec{u} + (\vec{h} \cdot \nabla)\vec{U} - \nu \Delta \vec{h} = (\vec{u} \cdot \nabla)\vec{H}, \\
(iii) & \text{Eq. (115) namely, momentum with } r \text{ eliminated via } S_{\mathcal{F}\lambda}^{-1}.
\end{array}
\right. 
\quad (116) $$

We arrive at the same result working with the characteristic $Y$-flow and the Jacobian determinant,

$$ \frac{d}{ds} Y(s; \vec{X}) = v(Y(s; \vec{X})), \quad v(\vec{X}) = -\frac{1}{2}\vec{X} - \vec{U}(\vec{X}), \quad Y(0; \vec{X}) = \vec{X}, \quad (117) $$

$$ J(s; \vec{X}) = \det \left( \partial_{\vec{X}} Y(s; \vec{X}) \right), \quad \frac{d}{ds} J(s; \vec{X}) = (\nabla \cdot v) (Y(s; \vec{X})) J(s; \vec{X}), \quad J(0; \vec{X}) = 1. \quad (118) $$

The resolvent is

$$ S_{\mathcal{F}\lambda}^{-1}[g](\vec{X}) = \int_{0}^{\infty} e^{-\lambda s} J(s; \vec{X}) g(Y(s; \vec{X})) ds, \quad (119) $$

which is the classical variation of constants representation for the scalar transport-dilation operator $S_{\mathcal{F}\lambda}[\cdot]$ along the flow $Y$. The integral is formal and converges under the usual spectral and decay hypotheses on $\Re\lambda$ and on the profile $(R, \vec{U})$.

The explicit elimination of $r$ via the resolvent is

$$ r(\vec{X}) = -S_{\mathcal{F}\lambda}^{-1}(\nabla \cdot (R\vec{u}))(\vec{X}) = -\int_{0}^{\infty} e^{-\lambda s} J(s; \vec{X})^{-1}(\nabla \cdot (R\vec{u}))(Y(s; \vec{X})) ds. \quad (120) $$

The non-local pressure-term in integral form becomes

$$ 
\begin{aligned}
\frac{1}{R(\vec{X})} \nabla (P'(R(\vec{X}))r(\vec{X})) &= -\frac{1}{R(\vec{X})} \nabla_{\vec{X}} \left\{ P'(R(\vec{X})) \int_{0}^{\infty} e^{-\lambda s} J(s; \vec{X})^{-1} (\nabla \cdot (R\vec{u})) (Y(s; \vec{X})) ds \right\} \\
&= -\int_{\mathbb{R}^3} K_{P,\lambda}(\vec{X}, \vec{Z}) \cdot \vec{u}(\vec{Z}) d\vec{Z},
\end{aligned} 
\quad (121) $$

with the pressure kernel given by the distributional identity,

$$ K_{P,\lambda}(\vec{X}, \vec{Z}) = \int_{0}^{\infty} e^{-\lambda s} J(s; \vec{X}) \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \nabla_{\vec{Z}} \cdot (R(\vec{Z})\delta(Y(s; \vec{X}) - \vec{Z})) \right\} ds. \quad (122) $$

Here $\delta$ is the Dirac delta [17]; the kernel is formal but explicit: integration in $s$ pushes the non-locality onto a bilinear kernel acting on $\vec{u}$.
The Lorentz magnetic non-local contribution is

$$ \frac{1}{\mu_0 R(\vec{X})} \left( (\nabla \times \vec{H}) \times \vec{h} - (\nabla \times \vec{h}) \times \vec{H} \right) - \frac{r(\vec{X})}{\mu_0 R(\vec{X})^2} \left( (\nabla \times \vec{H}) \times \vec{H} \right) \\ = \int_{\mathbb{R}^3} K_{M,\lambda}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z} + \int_{\mathbb{R}^3} \tilde{K}_M(\vec{X}, \vec{Z}) \vec{h}(\vec{Z}) d\vec{Z}, \quad (123) $$

where the magnetic kernel $K_{M,\lambda}$ contains the same resolvent weight $\int_0^\infty e^{-\lambda s} J^{-1}$, explicitly inferred through the substitution in (120), whilst $\tilde{K}_M$ is local in $\lambda$, and it comes from linear action of induction on $\vec{h}$.

The final non-local eigenproblem will be

$$ \lambda \vec{u}(\vec{X}) = \left( -\frac{1}{2} - \frac{1}{2}(\vec{X} \cdot \nabla) + (\vec{U} \cdot \nabla) \right) \vec{u}(\vec{X}) + (\vec{u} \cdot \nabla) \vec{U}(\vec{X}) + \mu \frac{1}{R(\vec{X})} \Delta \vec{u}(\vec{X}) \\ - \int_{\mathbb{R}^3} K_{P,\lambda}(\vec{X}, \vec{Z}) \cdot \vec{u}(\vec{Z}) d\vec{Z} + \int_{\mathbb{R}^3} K_{M,\lambda}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z} + \int_{\mathbb{R}^3} \tilde{K}_M(\vec{X}, \vec{Z}) \vec{h}(\vec{Z}) d\vec{Z}. \quad (124) $$

Ergo (116).

To connect the displayed Lorentz magnetic non-local decomposition with the iterated-kernel cascade, substitute the resolvent Eq. (119) into each pressure and magnetic contribution.

Let $\vec{Z}_0 = \vec{X}$. This anchor identifies the initial spatial argument for the characteristic chain and is used consistently below to indicate the evaluation point of the first resolvent iteration.

With this substitution every occurrence of $K_{P,\lambda}$ and $K_{M,\lambda}$ is expressed as a similarity-time weighted local density transported along the characteristic sequence $Y(s_j; \vec{Z}_{j-1})$, $j = 1, \dots, n$. The $n$-fold kernels that follow are therefore the natural compositions of these transported densities.

Let

$$ K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}) = \int_{0}^{\infty} ds_{1} \cdots \int_{0}^{\infty} ds_{n} \int_{\mathbb{R}^3} d\vec{Z}_1 \cdots \int_{\mathbb{R}^3} d\vec{Z}_{n-1} \prod_{j=1}^{n} e^{-\lambda s_j} J(s_j; \vec{Z}_{j-1})^{-1} \\ \times \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \nabla_{\vec{Z}_1} \cdot (R(\vec{Z}_1)\delta(Y(s_1; \vec{X}) - \vec{Z}_1)) \right. \\ \times \prod_{j=2}^{n-1} \left[ \nabla_{\vec{Z}_j} \cdot (R(\vec{Z}_j)\delta(Y(s_j; \vec{Z}_{j-1}) - \vec{Z}_j)) \right] \\ \left. \times \left[ \nabla_{\vec{Z}_n} \cdot (R(\vec{Z}_n)\delta(Y(s_n; \vec{Z}_{n-1}) - \vec{Z}_n)) \right] \right\}, \quad (125) $$

and

$$ - \int_{\mathbb{R}^3} K_{P,\lambda}(\vec{X}, \vec{Z}) \cdot \vec{u}(\vec{Z}) d\vec{Z} = - \sum_{n=1}^{\infty} \int_{\mathbb{R}^3} K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}) \cdot \vec{u}(\vec{Z}) d\vec{Z}. \quad (126) $$

Analogously, construct $K_{M,\lambda}^{(n)}(\vec{X}, \vec{Z})$ by replacing the pressure-type divergence factors in (125) with the magnetic integrands appearing in $K_{M,\lambda}$,

$$ \lambda \vec{u}(\vec{X}) = \left( -\frac{1}{2} - \frac{1}{2}(\vec{X} \cdot \nabla) + (\vec{U} \cdot \nabla) \right) \vec{u}(\vec{X}) + (\vec{u} \cdot \nabla) \vec{U}(\vec{X}) + \mu \frac{1}{R(\vec{X})} \Delta \vec{u}(\vec{X}) \\ + \sum_{n=1}^{\infty} \int_{\mathbb{R}^3} K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}) \cdot \vec{u}(\vec{Z}) d\vec{Z} + \sum_{n=1}^{\infty} \int_{\mathbb{R}^3} K_{M,\lambda}^{(n)}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z} \\ + \int_{\mathbb{R}^3} \tilde{K}_M(\vec{X}, \vec{Z}) \vec{h}(\vec{Z}) d\vec{Z}. \quad (127) $$

In addition,

$$ K_{P,\lambda}^{(1)}(\vec{X}, \vec{Z}) = \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \int_{0}^{\infty} e^{-\lambda s} J(s; \vec{X})^{-1} \nabla_{\vec{Z}} \cdot (R(\vec{Z})\delta(Y(s; \vec{X}) - \vec{Z})) ds \right\}, \quad (128) $$
and

$$ K_{P,\lambda}^{(2)}(\vec{X}, \vec{Z}) = \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \int_{0}^{\infty} ds_{1} \int_{0}^{\infty} ds_{2} \int_{\mathbb{R}^{3}} d\vec{Z}_{1} \right. \\ \times e^{-\lambda(s_{1}+s_{2})} J(s_{1}; \vec{X})^{-1} J(s_{2}; \vec{Z}_{1})^{-1} \\ \times \nabla_{\vec{Z}_{1}} \cdot (R(\vec{Z}_{1}) \delta(Y(s_{1}; \vec{X}) - \vec{Z}_{1})) \\ \left. \times \nabla_{\vec{Z}} \cdot (R(\vec{Z}) \delta(Y(s_{2}; \vec{Z}_{1}) - \vec{Z})) \right\}. \quad (129) $$

The resulting s-integrals and spatial convolutions in (127) converge in the same operator topology that validates the Neumann expansion of $(\mathbb{I} - T(\lambda))^{-1}$, where $\mathbb{I}$ is the identity matrix.

### 8.1.1.1. Green Closure of the Momentum Balance

In order to exhibit an explicit non-local closure for the similarity momentum operator, we invert the viscous Laplacian appearing in the momentum balance and express the profile $\vec{U}$ by means of the Newtonian potential. The momentum balance may be written in the similarity variables as

$$ R(-\alpha \vec{U} - \beta(\vec{X} \cdot \nabla)\vec{U} + (\vec{U} \cdot \nabla)\vec{U}) = \mu \Delta \vec{U} - \nabla \left( p + \frac{|\vec{B}|^2}{2\mu_0} \right) - \frac{1}{\mu_0} (\vec{B} \cdot \nabla)\vec{B}. \quad (130) $$

Inverting the operator $\mu \Delta$ with the fundamental solution of the Laplacian in $\mathbb{R}^3$ and recalling the decay hypotheses we write the integral Green formalism,

$$ \vec{U}(\vec{X}) = -\frac{1}{4\pi\mu} \int_{\mathbb{R}^3} \frac{\vec{F}(\vec{Y})}{|\vec{X} - \vec{Y}|} d\vec{Y}, \quad (131a) $$

$$ \vec{F}(\vec{Y}) = R \left( -\alpha \vec{U} - \beta(\vec{Y} \cdot \nabla)\vec{U} + (\vec{U} \cdot \nabla)\vec{U})(\vec{Y}) + \frac{1}{\mu_0}(\vec{B} \cdot \nabla)\vec{B}(\vec{Y}) - \nabla(p(\vec{Y}) + \frac{|\vec{B}(\vec{Y})|^2}{2\mu_0}) \right). \quad (131b) $$

The equational block (131) is exact under the stated hypotheses. Iterating the integral notation by substituting the right-hand non-linearities in $\vec{F}$ with their representations in turn produces higher-order convolution closures. Truncating at third order (and isolating the explicit leading-density contribution while leaving the purely magnetic coupling as a remainder) gives the triple-convolution expansion,

$$ \vec{U}(\vec{X}) = -\frac{1}{(4\pi\mu)^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \frac{R(\vec{Y})R(\vec{Z})R(\vec{W})}{|\vec{X} - \vec{Y}||\vec{Y} - \vec{Z}||\vec{Z} - \vec{W}|} \\ \times [-\alpha \vec{U}(\vec{W}) - \beta(\vec{W} \cdot \nabla)\vec{U}(\vec{W}) + (\vec{U}(\vec{W}) \cdot \nabla)\vec{U}(\vec{W})] d\vec{W} d\vec{Z} d\vec{Y} + \mathcal{O}\left(\frac{1}{\mu_0}\right), \quad (132) $$

where the magnetic remainder is given by the double-convolution term,

$$ \mathcal{O}\left(\frac{1}{\mu_0}\right) = \frac{1}{2\mu_0(4\pi\mu)^2} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \frac{(\vec{B}(\vec{Y}) \cdot \nabla)\vec{B}(\vec{Z}) - (\vec{B}(\vec{Z}) \cdot \nabla)\vec{B}(\vec{Y})}{|\vec{X} - \vec{Y}||\vec{Y} - \vec{Z}|} d\vec{Z} d\vec{Y}. \quad (133) $$

Eqq. (131) (132) (133) provide a non-local integral closure of the similarity momentum operator consistent with the spectral and resolvent constructions. The triple integral in (132) explicitly displays the iterated spatial transport of the density-weighted inertial term $R\vec{U}$ and reveals the structure of successive non-local interactions; the $\mathcal{O}(\frac{1}{\mu_0})$ term collects the purely magnetic non-local contributions deferred for separate spectral treatment.
### 8.1.1.2. Iterated Functional Hierarchy and $K_{P,\lambda}$-Cascade

The resolvent integral formulation and the associated kernel operators $K_{P,\lambda}$ and $K_{M,\lambda}$ naturally lead to an iterated integral structure when inserted recursively into the weak formulation of the momentum equation. By iterating the substitution

$$r = -S_{\bar{f},\lambda}^{-1}(\nabla \cdot (R\bar{u})) \quad (134)$$

within the pressure and Lorentz terms and reinserting the resulting non-local expressions back into the functional energy, one obtains a hierarchy of nested resolvent integrals. This provides the mathematical bridge between the non-local kernel formulation and the multi-linear functional structures that emerge in the eigenvalue problem.

The time-dependent functional (8) can now be expressed as a resolvent-induced hierarchy, generated by repeated application of the operator $S_{\bar{f},\lambda}^{-1}$ inside the pressure term. Specifically,

$$\begin{aligned} \mathcal{E}[t] = & \int_{\mathbb{R}^3} \left( \frac{1}{2} \rho |\vec{u}|^2 + \frac{1}{2\mu_0} |\vec{B}|^2 + \rho e(\rho) \right) d\vec{X} \\ & + \sum_{m=1}^{\infty} \frac{1}{m!} \int_0^{\infty} ds_1 \int_0^{\infty} ds_2 \cdots \int_0^{\infty} ds_m \int_{\mathbb{R}^3} d\vec{X}_0 \cdots \int_{\mathbb{R}^3} d\vec{X}_m \\ & \times \prod_{j=1}^{m} e^{-\lambda s_j} J(s_j; \vec{X}_j) [\nabla \cdot (R\bar{u})] (Y(s_j; \vec{X}_j)) \prod_{k=0}^{m-1} K_{P,\lambda}(\vec{X}_k, \vec{X}_{k+1}), \end{aligned} \quad (135)$$

where each integral over $s_j$ corresponds to an iteration of the scalar resolvent inversion, and the nested spatial integrals encode the cumulative non-local interactions generated by repeated substitution.

In Fourier space, the same hierarchy appears as an infinite chain of spectral convolutions,

$$\mathcal{E}[t] = \sum_{m=1}^{\infty} \frac{1}{m!} \int_{\mathbb{R}^3} d\xi_0 \cdots \int_{\mathbb{R}^3} d\xi_m \prod_{j=0}^{m} \hat{\rho}(\xi_j, t) \prod_{k=0}^{m-1} \widehat{K_{P,\lambda}}(\xi_k - \xi_{k+1}) \hat{u}(\xi_{k+1}, t). \quad (136)$$

Incorporating $\mathfrak{F}^s$-weights produces the regularized multi-resolvent expansion,

$$\begin{aligned} \mathcal{E}[t] = & \sum_{m=1}^{\infty} \frac{1}{m!} \int_{\mathbb{R}^3} d\xi_0 \cdots \int_{\mathbb{R}^3} d\xi_m \prod_{j=0}^{m} \Lambda^{-s} F \left\{ \Lambda^s \left[ \rho(\cdot, t) \left( \frac{1}{2} |\vec{u}|^2 + \frac{1}{2\mu_0} |\vec{B}|^2 \right) \right] \right\} (\xi_j) \\ & \times \prod_{k=0}^{m-1} \widehat{K_{P,\lambda}}(\xi_k - \xi_{k+1})^2. \end{aligned} \quad (137)$$

A more explicit description of the iterated structure arises by expanding the scalar resolvent itself as a Neumann-type series, applied to the non-linear transport-dilation operator,

$$S_{\bar{f},\lambda}^{-1} = \sum_{n=0}^{\infty} [S_{\bar{f},\lambda}^{-1} \circ \text{NL}]^n S_{\bar{f},\lambda}^{-1}, \quad (138)$$

where $\text{NL}[\phi] = \nabla \cdot (R\phi)$ is for the non-linear continuity perturbation. Inserting this into the pressure contribution produces the explicit n-th term,

$$\begin{aligned} r^{(n)}(\vec{X}) = & \frac{(-1)^{n+1}}{n!} \int_0^{\infty} ds_1 \cdots \int_0^{\infty} ds_{n+1} \int_{\mathbb{R}^3} d\vec{Z}_1 \cdots \int_{\mathbb{R}^3} d\vec{Z}_{n+1} \\ & \times e^{-\lambda(s_1 + \cdots + s_{n+1})} \prod_{j=1}^{n+1} J(s_j; \vec{Z}_j)^{-1} \nabla_{\vec{Z}_{n+1}} \cdot (R(\vec{Z}_{n+1})\bar{u}(\vec{Z}_{n+1})) \\ & \times \prod_{k=1}^{n} K_{P,\lambda}(\vec{Z}_k, \vec{Z}_{k+1}), \end{aligned} \quad (139)$$

so that the full density perturbation is recovered as $r = \sum_{n=0}^{\infty} r^{(n)}$.
The corresponding kernel composition admits a compact operatorial construct,

$$K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}) = \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_n K_{P,\lambda}(\vec{X}, \vec{Y}_1) K_{P,\lambda}(\vec{Y}_1, \vec{Y}_2) \cdots K_{P,\lambda}(\vec{Y}_n, \vec{Z}), \quad (140)$$

for which the resolvent expansion in kernel form is

$$K_{P,\lambda}^{\text{tot}}(\vec{X}, \vec{Z}) = \sum_{n=0}^{\infty} K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}). \quad (141)$$

The kernel composition propagates directly to the functional $\mathcal{E}[t]$, which can be rewritten as

$$\mathcal{E}[t] = \sum_{n=0}^{\infty} \frac{1}{n!} \int_{\mathbb{R}^3} d\vec{X} \int_{\mathbb{R}^3} d\vec{Z} \vec{u}(\vec{Z}, t) \cdot [K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}) \vec{B}(\vec{Z}, t)] \rho(\vec{X}, t) \\ + \text{magnetic and higher-order coupling terms.} \quad (142)$$

Each order $n$ corresponds to an $(n+1)$-fold convolution in space and an $(n+1)$-fold integral in similarity time, with a hierarchical chain that governs the spectral dynamics of the eigenmodes.

### 8.1.2. Regularized Determinants of Integral Operators

We formulate a single compact integral-operator whose spectral equation is equivalent to the coupled problem and we introduce a Carleman-regularized determinant [78] capturing the spectrum associated with Eq. (116).

Let

$$\mathcal{D}_{\text{loc}}^{\text{mom}}[\vec{v}] = -\frac{1}{2}\vec{v} - \frac{1}{2}(\vec{X} \cdot \nabla)\vec{v} + (\vec{U} \cdot \nabla)\vec{v} + (\vec{v} \cdot \nabla)\vec{U} + \mu \frac{1}{R(\vec{X})} \left[ \Delta \vec{v} + \left(\nabla \ln R(\vec{X})\right) \cdot \nabla \vec{v} \right], \quad (143)$$

where $\mathcal{D}_{\text{loc}}^{\text{mom}}$ is the local differential operator appearing in the reduced momentum equation.

Let

$$a_{\text{ind}}(\lambda)[\vec{h}] = \lambda \vec{h} + \frac{1}{2}\vec{h} + \frac{1}{2}(\vec{X} \cdot \nabla)\vec{h} - (\vec{U} \cdot \nabla)\vec{h} + (\vec{h} \cdot \nabla)\vec{U} - \frac{\nu}{R(\vec{X})} \Delta \vec{h}, \quad (144)$$

being $a_{\text{ind}}(\lambda)$ the linear induction operator acting on $\vec{h}$.

Let

$$w[\vec{u}] = (\vec{u} \cdot \nabla)\vec{H} - (\vec{H} \cdot \nabla)\vec{u}, \quad (145)$$

where $w$ is the bilinear mapping $u \mapsto (\vec{u} \cdot \nabla)\vec{H} - (\vec{H} \cdot \nabla)\vec{u}$.

The induction resolvent (when invertible) is $a_{\text{ind}}(\lambda)^{-1}$ and meets

$$\vec{h} = a_{\text{ind}}(\lambda)^{-1} (w[\vec{u}]). \quad (146)$$

Set $T(\lambda)$ as the composite integral operator acting on $\vec{u}$ by

$$T(\lambda)[\vec{u}](\vec{X}) = (\lambda \mathbb{I} - \mathcal{D}_{\text{loc}}^{\text{mom}})^{-1} \left\{ - \int_{\mathbb{R}^3} K_{P,\lambda}(\vec{X}, \vec{Z}) \cdot \vec{u}(\vec{Z}) d\vec{Z} \right. \\ \qquad \left. + \int_{\mathbb{R}^3} K_{M,\lambda}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z} \right. \\ \qquad \left. + \int_{\mathbb{R}^3} \tilde{K}_M(\vec{X}, \vec{Z}) (a_{\text{ind}}(\lambda)^{-1} w[\vec{u}]) (\vec{Z}) d\vec{Z} \right\}, \quad (147)$$

where the inverse $(\lambda \mathbb{I} - \mathcal{D}_{\text{loc}}^{\text{mom}})^{-1}$ is understood as the resolvent (Green operator) of the local differential operator appearing as left-hand side of the reduced momentum equation, whereas $K_{P,\lambda}$, $K_{M,\lambda}$, and $\tilde{K}_M$ are the kernels of the respective integral operators.

Let

$$f_c(\lambda) = \det_2(\mathbb{I} - T(\lambda)), \quad (148)$$

be the Carleman 2-regularized determinant of $\mathbb{I} - T(\lambda)$ well-defined when $T(\lambda)$ is Hilbert–Schmidt trace-class [59] after the hypotheses on profiles and $\Re\lambda$.

Zeros of $f_c(\lambda)$ correspond to values $\lambda$ where 1 is in the spectrum of $T(\lambda)$.
Note that

$$ \log_{\mathcal{F}_c}(\lambda) = - \sum_{m=1}^{\infty} \frac{1}{m} \mathrm{tr}(T(\lambda)^m) \qquad (149) $$

is valid in the trace-class regime and it gives an explicit computable expansion: the traces are integrals of kernel-diagonals.

If at $\lambda_0$ the operator $T(\lambda_0)$ has an algebraically simple eigenvalue 1 with right eigenfunction $\Phi$ and left eigenfunction $\Psi$ normalized by $\langle\Psi, \Phi\rangle_{L_R^2} = 1$, then the first-order eigenvalue shift under a small perturbation $\delta$ is

$$ \delta\lambda = - \frac{\langle\Psi, \delta T(\lambda_0)\Phi\rangle_{L_R^2}}{\langle\Psi, \partial_\lambda T(\lambda_0)\Phi\rangle_{L_R^2}}. \qquad (150) $$

Let $Q = (\mathbb{I} - T(\lambda_0))^{-1}$ be the reduced resolvent on the $L_R^2$-orthogonal complement of span $\{\Phi\}$. A formal second-order correction is provided by

$$ \delta\lambda^{(2)} = - \frac{\langle\Psi, \delta^2 T(\lambda_0)\Phi\rangle_{L_R^2} + 2\langle\Psi, \delta T(\lambda_0)Q\delta T(\lambda_0)\Phi\rangle_{L_R^2}}{2\langle\Psi, \partial_\lambda T(\lambda_0)\Phi\rangle_{L_R^2}}, \qquad (151) $$

which is the standard second-order Rayleigh-Schrödinger type correction [70] [75] written for the non-linear in $\lambda$ operator $T(\lambda)$. Attention: the formula requires existence of $Q$ on the complementary subspace and is to be interpreted in the usual projected sense.

### 8.1.2.1. Resolvent and Function-operator Nested Integral Cascade

We present a unified operator-functional integral expansion that unfolds the reduced resolvent $(\mathbb{I} - T(\lambda))^{-1}$ in terms of spatial convolution chains built from the kernels appearing in the definition of $T(\lambda)$, and displays the corresponding action on the physical profile $\vec{u}$ as an integral expansion.

Let

$$ T(\lambda)[\vec{u}](\vec{X}) = (\lambda\mathbb{I} - \mathcal{D}_{\mathrm{loc}}^{\mathrm{mom}})^{-1} \left\{ - \int_{\mathbb{R}^3} K_{P,\lambda}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z} + \int_{\mathbb{R}^3} K_{M,\lambda}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z} \right. \\ \left. + \int_{\mathbb{R}^3} \tilde{K}_M(\vec{X}, \vec{Z}) \circ (a_{\mathrm{ind}}(\lambda)^{-1} w[\vec{u}](\vec{Z})) d\vec{Z} \right\}. \qquad (152) $$

We regard each integral kernel appearing in (152) as a building block for the convolution chains which produce the kernels of the powers $T(\lambda)^m$. For $m \ge 1$, the kernel of $T(\lambda)^m$ is the $(m-1)$-fold spatial convolution of the composite integrand in (152); for $m \ge 1$, it is

$$ [T(\lambda)^m](\vec{X}; \vec{Z}) = \int_{\mathbb{R}^3} d\vec{Y}_1 \int_{\mathbb{R}^3} d\vec{Y}_2 \cdots \int_{\mathbb{R}^3} d\vec{Y}_{m-1} K_{T,\lambda}(\vec{X}, \vec{Y}_1; \lambda) K_{T,\lambda}(\vec{Y}_1, \vec{Y}_2; \lambda) \cdots K_{T,\lambda}(\vec{Y}_{m-1}, \vec{Z}; \lambda), \qquad (153) $$

where the composite pointwise kernel $K_{T,\lambda}(\cdot, \cdot; \lambda)$ is the integrand operator in (152); explicitly it is the mapping

$$ K_{T,\lambda}(\vec{X}, \vec{Z}; \lambda) = (\lambda\mathbb{I} - \mathcal{D}_{\mathrm{loc}}^{\mathrm{mom}})^{-1} \left( -K_{P,\lambda}(\vec{X}, \vec{Z}) + K_{M,\lambda}(\vec{X}, \vec{Z}) + \tilde{K}_M(\vec{X}, \vec{Z}) \circ (a_{\mathrm{ind}}(\lambda)^{-1} w[\cdot]) \right). \qquad (154) $$

The Neumann expansion of the reduced resolvent is thereby

$$ (\mathbb{I} - T(\lambda))^{-1} = \sum_{m=0}^{\infty} T(\lambda)^m, \qquad (155) $$

in the operator-norm regime where the series converges.
Acting on a physical profile $\vec{u}$, the resolvent produces the following nested integral cascade,

$$ (\mathbb{I} - T(\lambda))^{-1}[\vec{u}](\vec{X}) = \vec{u}(\vec{X}) + \sum_{m=1}^{\infty} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_m K_{T,\lambda}(\vec{X}, \vec{Y}_1; \lambda) \\ \times K_{T,\lambda}(\vec{Y}_1, \vec{Y}_2; \lambda) \cdots K_{T,\lambda}(\vec{Y}_{m-1}, \vec{Y}_m; \lambda) \vec{u}(\vec{Y}_m). \quad (156) $$

In the trace-class regime the Carleman 2-regularized determinant expansion is the standard trace expansion,

$$ \mathrm{tr}(T(\lambda)^m) = \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_m \mathrm{tr}_{\mathbb{R}^d} \left( \underbrace{K_{T,\lambda}(\vec{Y}_1, \vec{Y}_2; \lambda) K_{T,\lambda}(\vec{Y}_2, \vec{Y}_3; \lambda) \cdots K_{T,\lambda}(\vec{Y}_m, \vec{Y}_1; \lambda)}_{m \text{ matrix factors}} \right), \quad (157) $$

with the finite-dimensional component trace tr$_{R^d}$.

We now combine this operator cascade with a functional Picard procedure for the self-similar profile components. Let $\vec{u}^{(0)}$ be a chosen initial profile and delineate the iterates

$$ \vec{u}^{(0)} = \vec{u}_{\mathrm{init}}, \quad (158) $$

$$ \vec{u}^{(n+1)} = \vec{u}_{\mathrm{init}} + T(\lambda)[\vec{u}^{(n)}], \quad (159) $$

for $n \ge 0$. Each iterate $\vec{u}^{(n)}$ is expressible as a finite sum of integrals by substituting repeatedly the explicit integral of $T(\lambda)$ from (152).

A corresponding iteration for the truncated resolvent is got from

$$ R^{(0)}(\lambda) = \mathbb{I}, \quad (160) $$

$$ R^{(n+1)}(\lambda) = \mathbb{I} + T(\lambda)R^{(n)}(\lambda), \quad (161) $$

so that $R^{(n)}(\lambda)$ has kernel equal to the finite sum $\sum_{m=0}^n [T(\lambda)^m](\cdot; \cdot)$ and converges formally to $(\mathbb{I} - T(\lambda))^{-1}$ as $n \to \infty$ in the regime of convergence.

All of the series above are formal unless the operator $T(\lambda)$ satisfies suitable smallness or Schatten-class hypotheses [72].

Sufficient conditions guaranteeing absolute summability of the integrals include Hilbert–Schmidt bounds on each constituent kernel and an operator-norm bound of the form $\|T(\lambda)\|_{\mathrm{op}} < 1$ on the function space in which the operators act; under such hypotheses the identities (155) (157) and the Picard iterations (159) converge legitimately.

### 8.1.2.2. Kernel-Based Analysis and Spectral Projectors

#### I. Iterated Chain Representation and Fredholm Determinant Series

The fixed-point formulation associated with Eq. (152) admits a fully explicit series expansion as an infinite iterated kernel series. The operator $T(\lambda)$ raised to its $n$-th power has the integral realization,

$$ [T(\lambda)^n[\vec{u}^{(0)}]](\vec{X}) = \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_n K_{T,\lambda}(\vec{X}, \vec{Y}_1) \cdots K_{T,\lambda}(\vec{Y}_{n-1}, \vec{Y}_n) \vec{u}^{(0)}(\vec{Y}_n), \quad (162) $$

and the iterated solution associated with the resolvent formulation can be written as the infinite series,

$$ \vec{u}(\vec{X}) = \sum_{n=0}^{\infty} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_n K_{T,\lambda}(\vec{X}, \vec{Y}_1) \cdots K_{T,\lambda}(\vec{Y}_{n-1}, \vec{Y}_n) \vec{u}^{(0)}(\vec{Y}_n). \quad (163) $$

When $K_{T,\lambda}$ is itself generated by the scalar resolvent $S_{\mathbb{F}\lambda}^{-1}$ through the kernel composition $K_{T,\lambda} = K_{P,\lambda} + K_{M,\lambda}$, each integral factor in (162) expands into a pair of similarity-time integrals and a non-local transport term. Hence

$$ \begin{aligned} [T(\lambda)^n[\vec{u}^{(0)}]](\vec{X}) = & \int_0^\infty ds_1 \cdots \int_0^\infty ds_n \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_n \prod_{j=1}^n e^{-\lambda s_j} J(s_j; \vec{Y}_j)^{-1} \\ & \times K_{P,\lambda}(\vec{X}, \vec{Y}_1) K_{M,\lambda}(\vec{Y}_1, \vec{Y}_2) \cdots K_{P,\lambda}(\vec{Y}_{n-1}, \vec{Y}_n) \vec{u}^{(0)}(\vec{Y}_n). \end{aligned} \quad (164) $$
The Fourier-space counterpart follows directly,

$$T(\lambda)^n[\widehat{\vec{u}}^{(0)}](\xi_0) = \int_{\mathbb{R}^3} d\xi_1 \cdots \int_{\mathbb{R}^3} d\xi_n \widehat{K}_{T,\lambda}(\xi_0 - \xi_1) \cdots \widehat{K}_{T,\lambda}(\xi_{n-1} - \xi_n) \widehat{\vec{u}}^{(0)}(\xi_n). \quad (165)$$

The spectral resolvent itself admits a convergent Neumann series. For the non-homogeneous fixed-point problem,

$$(\mathbb{I} - T(\lambda))[\vec{u}] = \vec{u}^{(0)}, \quad (166)$$

one has

$$\vec{u} = \sum_{n=0}^{\infty} T(\lambda)^n [\vec{u}^{(0)}], \quad (\mathbb{I} - T(\lambda))^{-1} = \sum_{n=0}^{\infty} T(\lambda)^n, \quad (167)$$

whereas the homogeneous spectral problem,

$$(\mathbb{I} - T(\lambda))[\vec{u}] = 0, \quad (168)$$

is treated via the Fredholm condition [22] $\det(\mathbb{I} - T(\lambda)) = 0$, with no $\vec{u}^{(0)}$ present.

Expanding the determinant by the standard trace-log identity yields a multi-integral secular series,

$$\log \det (\mathbb{I} - T(\lambda)) = - \sum_{m=1}^{\infty} \frac{1}{m} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_m \operatorname{tr} [K_{T,\lambda}(\vec{Y}_1, \vec{Y}_2) \cdots K_{T,\lambda}(\vec{Y}_m, \vec{Y}_1)]. \quad (169)$$

The spectral structure of the non-local momentum operator is manifested through an infinite hierarchy of space-time integrals, each order corresponding to an iterated composition of the physical kernels $K_{P,\lambda}$ and $K_{M,\lambda}$ and their Fourier convolutions. The resulting series is the analytic bridge between the kernel formulation and the eigenvalue spectrum.

We come now to the spectral condition $\det(\mathbb{I} - T(\lambda)) = 0$. It admits a deeper analytical expansion that connects the non-local kernel composition to the eigenvalue spectrum through the regularized trace-log hierarchy. Starting from the standard identity,

$$\log \det (\mathbb{I} - T(\lambda)) = - \sum_{m=1}^{\infty} \frac{1}{m} \operatorname{tr} [T(\lambda)^m], \quad (170)$$

each trace term admits a fully expanded kernel chain,

$$\operatorname{tr} [T(\lambda)^m] = \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_m \operatorname{tr} [K_{T,\lambda}(\vec{Y}_1, \vec{Y}_2) K_{T,\lambda}(\vec{Y}_2, \vec{Y}_3) \cdots K_{T,\lambda}(\vec{Y}_m, \vec{Y}_1)]. \quad (171)$$

The expansion organizes into a hierarchy of iterated kernel products, each term encoding the spectral contribution of an m-fold propagation chain. Because the integral domain is cyclic under the trace, one can rewrite the secular series in terms of a cyclic convolution integral,

$$\log \det (\mathbb{I} - T(\lambda)) = - \sum_{m=1}^{\infty} \frac{1}{m} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots d\vec{Y}_m \\ \times \operatorname{tr} [K_{P,\lambda}(\vec{Y}_1, \vec{Y}_2) K_{M,\lambda}(\vec{Y}_2, \vec{Y}_3) \cdots K_{P,\lambda}(\vec{Y}_m, \vec{Y}_1)]. \quad (172)$$

In the Fourier domain, the same hierarchy becomes a spectral convolution series,

$$\log \det (\mathbb{I} - T(\lambda)) = - \sum_{m=1}^{\infty} \frac{1}{m} \int_{\mathbb{R}^3} d\xi_1 \cdots d\xi_m \operatorname{tr} [\widehat{K}_{T,\lambda}(\xi_1 - \xi_2) \cdots \widehat{K}_{T,\lambda}(\xi_m - \xi_1)]. \quad (173)$$

Each term of the hierarchy corresponds to a collective resonance mode, and the vanishing of the Fredholm determinant translates to the existence of a non-trivial eigenmode of the non-local momentum operator. The kernel-chain not only encompasses the full spectral data but also provides the functional analytic connection between the original integral formulation and the spectral theory.
Scholium 8.2 (Higher functional variations and non-local operator differentiation). Beyond the fixed-point formulation, one may inspect the higher functional sensitivity of the non-local integral operator with respect to perturbations of the field $\vec{u}$. Let $\delta\vec{u}$ be an admissible variation in the function space of solutions. The first-order functional derivative of the operator Eq. (166) captures the linearized response,

$$ (\mathbb{I} - T(\lambda))[\delta\vec{u}] - \int_{\mathbb{R}^3} d\vec{Y} \frac{\delta K_{T,\lambda}(\vec{X}, \vec{Y})}{\delta \vec{u}(\vec{Y})} \vec{u}(\vec{Y}) = 0. \quad (174) $$

Iterating the variation, one obtains the m-th order functional derivative of the operator composition,

$$ \frac{\delta^m}{\delta \vec{u}(\vec{Y}_1) \cdots \delta \vec{u}(\vec{Y}_m)} \left[ (\mathbb{I} - T(\lambda))[\vec{u}] \right](\vec{X}) = - \int_{\mathbb{R}^3} d\vec{Z} \frac{\delta^m K_{T,\lambda}(\vec{X}, \vec{Z})}{\delta \vec{u}(\vec{Y}_1) \cdots \delta \vec{u}(\vec{Y}_m)} \vec{u}(\vec{Z}) \\ - \sum_{j=1}^{m} \int_{\mathbb{R}^3} d\vec{Z} \frac{\delta^{m-1} K_{T,\lambda}(\vec{X}, \vec{Z})}{\delta \vec{u}(\vec{Y}_1) \cdots \delta \vec{u}(\vec{Y}_j) \cdots \delta \vec{u}(\vec{Y}_m)} \delta \vec{u}(\vec{Z}), \quad (175) $$

where the hat $\widehat{\delta\vec{u}(\vec{Y}_j)}$ indicates omission of the corresponding derivative.

This expansion is a functional Taylor series for the operator equation in the neighborhood of a solution $\vec{u}$. Each term captures non-linear feedback between the kernel $K_{T,\lambda}$ and the field configuration, and higher-order terms encode the full non-local response structure of the integral problem under perturbations of $\vec{u}$. ◇

## II. Double Composition

The spectral structure associated with the resolvent operator admits an even deeper functional layer if one considers the bilinear composition of two resolvents at distinct spectral parameters. For $\lambda, \mu \in \mathbb{C}$ such that the respective resolvents are defined, the composition reads

$$ (\mathbb{I} - T(\lambda))^{-1} (\mathbb{I} - T(\mu))^{-1} = \sum_{n=0}^{\infty} \sum_{m=0}^{\infty} T(\lambda)^n T(\mu)^m, \quad (176) $$

and each term in the double series possesses a kernel with $n+m$ spatial integrals.

For the action on the initial state $\vec{u}^{(0)}$,

$$ \left[ \left( (\mathbb{I} - T(\lambda))^{-1} (\mathbb{I} - T(\mu))^{-1} [\vec{u}^{(0)}] \right) (\vec{X}) \right. \\ = \sum_{n,m=0}^{\infty} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_{n+m} \left( \underbrace{K_{T,\lambda}(\vec{X}, \vec{Y}_1) K_{T,\lambda}(\vec{Y}_1, \vec{Y}_2) \cdots K_{T,\lambda}(\vec{Y}_{n-1}, \vec{Y}_n)}_{n \text{ } \lambda\text{-factors}} \right. \\ \left. \times \underbrace{K_{T,\mu}(\vec{Y}_n, \vec{Y}_{n+1}) K_{T,\mu}(\vec{Y}_{n+1}, \vec{Y}_{n+2}) \cdots K_{T,\mu}(\vec{Y}_{n+m-1}, \vec{Y}_{n+m})}_{m \text{ } \mu\text{-factors}} \right) \vec{u}^{(0)}(\vec{Y}_{n+m}). \quad (177) $$

The corresponding Fourier-space expression follows from repeated convolution,

$$ \widehat{[\cdot]}(\xi_0) = \sum_{n,m=0}^{\infty} \int_{\mathbb{R}^3} d\xi_1 \cdots \int_{\mathbb{R}^3} d\xi_{n+m} \left( \underbrace{\widehat{K_{T,\lambda}}(\xi_0 - \xi_1) \widehat{K_{T,\lambda}}(\xi_1 - \xi_2) \cdots \widehat{K_{T,\lambda}}(\xi_{n-1} - \xi_n)}_{n \text{ } \lambda\text{-factors}} \right. \\ \left. \times \underbrace{\widehat{K_{T,\mu}}(\xi_n - \xi_{n+1}) \widehat{K_{T,\mu}}(\xi_{n+1} - \xi_{n+2}) \cdots \widehat{K_{T,\mu}}(\xi_{n+m-1} - \xi_{n+m})}_{m \text{ } \mu\text{-factors}} \right) \widehat{\vec{u}^{(0)}}(\xi_{n+m}), \quad (178) $$

$$ \text{where } \widehat{[\cdot]} \stackrel{\text{eqv}}{=} \left[ \left( (\mathbb{I} - T(\lambda))^{-1} (\mathbb{I} - T(\mu))^{-1} [\vec{u}^{(0)}] \right) \right]. $$

This double-resolvent expansion embodies correlations between two distinct spectral layers of the operator $T(\cdot)$ and generates a hierarchical chain of multi-integrals. Each term simultaneously
probes the $\lambda$- and $\mu$-dependence of the spectrum, providing the analytic bridge to bilinear spectral forms and two-point spectral densities.

### III. Spectral Projectors and Contour Integral Sketch

A further refinement of the spectral theory associated with $T(\lambda)$ arises via constructing the spectral projectors through the contour integral of the resolvent. Let $\Gamma_c$ be a simple closed contour in the complex $\lambda$-plane enclosing an isolated eigenvalue $\lambda_*$ of the spectrum of $T(\lambda)$. The corresponding spectral projector is

$$ \mathbb{P}_{\lambda_*} = \frac{1}{2\pi i} \oint_{\Gamma_c} (\mathbb{I} - T(\zeta))^{-1} d\zeta. \qquad (179) $$

In the analytic Fredholm framework, a more accurate form of the same projector includes the derivative of $T(\zeta)$ inside the contour integral,

$$ \mathbb{P}_{\lambda_*} = \frac{1}{2\pi i} \oint_{\Gamma_c} (\mathbb{I} - T(\zeta))^{-1} \frac{dT(\zeta)}{d\zeta} d\zeta, \qquad (180) $$

which holds whenever $(\mathbb{I} - T(\lambda))^{-1}$ is meromorphic in $\lambda$ and $\lambda_*$ is an isolated simple pole of the resolvent. Under this hypothesis the contour integral in (179) extracts precisely the corresponding spectral projector.

Expanding the resolvent under the integral using the Neumann series and inserting the kernel for each power, one reaches a formula for $\mathbb{P}_{\lambda_*}$ as an infinite series of nested space-parameter integrals,

$$ [\mathbb{P}_{\lambda_*} \vec{u}^{(0)}](\vec{X}) = \frac{1}{2\pi i} \oint_{\Gamma_c} d\zeta \sum_{n=0}^{\infty} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_n K_{T,\zeta}(\vec{X}, \vec{Y}_1) \cdots K_{T,\zeta}(\vec{Y}_{n-1}, \vec{Y}_n) \vec{u}^{(0)}(\vec{Y}_n). \quad (181) $$

In Fourier variables this becomes

$$ [\widehat{\mathbb{P}_{\lambda_*} \vec{u}^{(0)}}](\xi_0) = \frac{1}{2\pi i} \oint_{\Gamma_c} d\zeta \sum_{n=0}^{\infty} \int_{\mathbb{R}^3} d\xi_1 \cdots \int_{\mathbb{R}^3} d\xi_n \widehat{K_{T,\zeta}}(\xi_0 - \xi_1) \cdots \widehat{K_{T,\zeta}}(\xi_{n-1} - \xi_n) \widehat{\vec{u}^{(0)}}(\xi_n). \quad (182) $$

This formulation exhibits the spectral projectors as contour integrals of infinite-dimensional integral operators and encodes the local spectral structure of $T(\lambda)$ directly in terms of its kernel chains. Each term of the expansion is a nested multi-integral coupling spatial variables and the spectral parameter $\zeta$, thereby providing a direct analytic construction of eigenmodes and generalized eigenspaces.

**Scholium 8.3 (Asymptotic expansion of the spectral projector).** The contour integral of the projector admits a fully explicit series expansion in terms of derivatives of $T(\lambda)$ at the isolated eigenvalue $\lambda_*$. By expanding $(\mathbb{I} - T(\zeta))^{-1}$ as a Laurent series around $\zeta = \lambda_*$, one gets

$$ \mathbb{P}_{\lambda_*} = \sum_{p=0}^{\infty} \frac{1}{p! 2\pi i} \oint_{\Gamma_c} \frac{d\zeta}{(\zeta - \lambda_*)^{p+1}} \frac{d^p}{d\lambda^p} \left[ (\mathbb{I} - T(\lambda))^{-1} \right]_{\lambda=\lambda_*} . \qquad (183) $$

Each derivative of the inverse operator can be computed recursively from the identity

$$ \frac{d}{d\lambda} (\mathbb{I} - T(\lambda))^{-1} = (\mathbb{I} - T(\lambda))^{-1} \frac{dT(\lambda)}{d\lambda} (\mathbb{I} - T(\lambda))^{-1}, \qquad (184) $$

which implies

$$ \frac{d^p}{d\lambda^p} (\mathbb{I} - T(\lambda))^{-1} = \sum_{r=1}^{p} \sum_{\substack{j_1 + \cdots + j_r = p \\ j_k \ge 1}} (\mathbb{I} - T(\lambda))^{-1} \frac{d^{j_1} T(\lambda)}{d\lambda^{j_1}} (\mathbb{I} - T(\lambda))^{-1} \cdots \frac{d^{j_r} T(\lambda)}{d\lambda^{j_r}} (\mathbb{I} - T(\lambda))^{-1}. \qquad (185) $$

Substituting (185) into (183) it is determined a series involving only the original operator $T(\lambda)$, its derivatives $\frac{d^j T(\lambda)}{d\lambda^j}$, and the resolvent factor $(\mathbb{I} - T(\lambda))^{-1}$ evaluated at $\lambda_*$. The result expresses $\mathbb{P}_{\lambda_*}$ as an infinite operator series elucidating the full local structure of the eigenspace associated with $\lambda_*$.
A complementary perspective expresses the same projector directly in terms of the functional calculus of the operator. When the spectrum of $T(\lambda)$ is discrete in a neighborhood of $\lambda_*$, the contour integral (179) is equivalent to taking the residue of the resolvent operator at the simple pole $\lambda = \lambda_*$,

$$ \mathbb{P}_{\lambda_*} = \mathrm{res}_{\lambda=\lambda_*} \left[ \left( \mathbb{I} - T(\lambda) \right)^{-1} \right], \qquad (186) $$

which connects the integral of $\mathbb{P}_{\lambda_*}$ to the analytic functional calculus and highlights its role as the spectral component of the identity associated with the isolated eigenvalue $\lambda_*$. It also provides a direct bridge between the Laurent expansion (183) and the spectral decomposition of $T(\lambda)$. ◇

## IV. Faà di Bruno Hierarchies

The parametric behaviour of the operator family $(\mathbb{I} - T(\lambda))^{-1}$ as a function of the complex parameter $\lambda$ can be fully described by the classical combinatorial machinery of Faà di Bruno [19] [20], adapted to operator calculus. Since

$$ (\mathbb{I} - T(\lambda))^{-1} = \sum_{n=0}^{\infty} T(\lambda)^n, \qquad (187) $$

we may differentiate term by term and reorganize the resulting series. The first derivative turns out to be

$$ \frac{d}{d\lambda} (\mathbb{I} - T(\lambda))^{-1} = \sum_{n=1}^{\infty} \sum_{k=0}^{n-1} T(\lambda)^k \frac{dT(\lambda)}{d\lambda} T(\lambda)^{n-1-k}, \qquad (188) $$

and acting on the initial state $\vec{u}^{(0)}$ gives

$$ \left[ \frac{d}{d\lambda} (\mathbb{I} - T(\lambda))^{-1} [\vec{u}^{(0)}] \right] (\vec{X}) = \sum_{n=1}^{\infty} \sum_{k=0}^{n-1} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots d\vec{Y}_{n+1} \\ \times K_{T,\lambda}(\vec{X}, \vec{Y}_1) \cdots \frac{\partial K_{T,\lambda}}{\partial \lambda}(\vec{Y}_k, \vec{Y}_{k+1}) \cdots K_{T,\lambda}(\vec{Y}_n, \vec{Y}_{n+1}) \vec{u}^{(0)}(\vec{Y}_{n+1}). \qquad (189) $$

Higher derivatives follow from repeated differentiation. The p-th derivative follows from the recursive non-commutative Faà di Bruno identity,

$$ \frac{d^p}{d\lambda^p} (\mathbb{I} - T(\lambda))^{-1} = \sum_{r=1}^{p} \sum_{\substack{j_1+\cdots+j_r=p \\ j_k \ge 1}} (\mathbb{I} - T(\lambda))^{-1} \frac{d^{j_1} T(\lambda)}{d\lambda^{j_1}} (\mathbb{I} - T(\lambda))^{-1} \cdots \frac{d^{j_r} T(\lambda)}{d\lambda^{j_r}} (\mathbb{I} - T(\lambda))^{-1}. \qquad (190) $$

Each term in (190) translates into a multi-integral involving $\partial_{\lambda}^j K_{T,\lambda}$ for $1 \le j \le p$, interspersed with powers of $K_{T,\lambda}$. Acting on $\vec{u}^{(0)}$,

$$ \left[ \frac{d^p}{d\lambda^p} (\mathbb{I} - T(\lambda))^{-1} [\vec{u}^{(0)}] \right] (\vec{X}) = p! \sum_{m=p}^{\infty} \sum_{z(p)} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots d\vec{Y}_{m+p} \\ \times K_{T,\lambda}(\vec{X}, \vec{Y}_1) \cdots \frac{\partial^{j_1} K_{T,\lambda}}{\partial \lambda^{j_1}}(\vec{Y}_{i_1}, \vec{Y}_{i_1+1}) \cdots \frac{\partial^{j_p} K_{T,\lambda}}{\partial \lambda^{j_p}}(\vec{Y}_{i_p}, \vec{Y}_{i_p+1}) \\ \cdots \vec{u}^{(0)}(\vec{Y}_{m+p}). \qquad (191) $$

The parametric derivative of the Fredholm determinant follows similarly. Differentiating the trace-log identity,

$$ \frac{d}{d\lambda} \log \det (\mathbb{I} - T(\lambda)) = -\mathrm{tr} \left[ (\mathbb{I} - T(\lambda))^{-1} \frac{dT(\lambda)}{d\lambda} \right], \qquad (192) $$

and applying Faà di Bruno combinatorics to higher orders unfolds a hierarchy of trace invariants coupling powers of $(\mathbb{I} - T(\lambda))^{-1}$ with successive derivatives of $T(\lambda)$. Each such invariant admits a nested kernel expansion with iterated $\lambda$-derivatives of $K_{T,\lambda}$.
In Fourier space the first parametric derivative acting on $\vec{u}^{(0)}$ works as

$$ \frac{d}{d\lambda} \widehat{[\cdot]}(\xi_0) = \sum_{n=1}^{\infty} \sum_{k=0}^{n-1} \int_{\mathbb{R}^3} d\xi_1 \cdots d\xi_{n+1} \\ \times \widehat{K_{T,\lambda}}(\xi_0 - \xi_1) \cdots \frac{\partial \widehat{K_{T,\lambda}}}{\partial \lambda}(\xi_k - \xi_{k+1}) \cdots \widehat{K_{T,\lambda}}(\xi_n - \xi_{n+1}) \widehat{u}^{(0)}(\xi_{n+1}), \quad (193) \\ \text{where } \widehat{[\cdot]} \stackrel{\text{eqv}}{=} \left[ (\mathbb{I} - T(\lambda))^{-1} [\widehat{u}^{(0)}] \right]. $$

**Scholium 8.4 (Higher-order derivative expansion of the trace-log).** We can push the differentiation of the trace-log identity to arbitrary order while keeping the operator structure unchanged. By repeated application of the Leibniz rule and the cyclicity of the trace, the p-th derivative of the trace-log satisfies

$$ \frac{d^p}{d\lambda^p} \log \det (\mathbb{I} - T(\lambda)) = - \sum_{r=1}^{p} (-1)^{r-1} (r-1)! \sum_{\substack{j_1+\cdots+j_r=p \\ j_k \ge 1}} \frac{p!}{j_1! \cdots j_r!} \\ \times \operatorname{tr} \left[ \left( \frac{d^{j_1} T(\lambda)}{d\lambda^{j_1}} \right) (\mathbb{I} - T(\lambda))^{-1} \right. \\ \left. \times \left( \frac{d^{j_2} T(\lambda)}{d\lambda^{j_2}} \right) (\mathbb{I} - T(\lambda))^{-1} \cdots \left( \frac{d^{j_r} T(\lambda)}{d\lambda^{j_r}} \right) (\mathbb{I} - T(\lambda))^{-1} \right]. \quad (194) $$

This identity assumes that $(\mathbb{I} - T(\lambda))^{-1}$ is meromorphic in $\lambda$ and that $\lambda_*$ is a simple pole of the resolvent. Here the summation is taken over all $r$-tuples of positive integers $(j_1, \dots, j_r)$ such that $j_1 + \cdots + j_r = p$. Each term consists only of the original operator $T(\lambda)$, its successive derivatives $\frac{d^j T(\lambda)}{d\lambda^j}$, and the resolvent factor $(\mathbb{I} - T(\lambda))^{-1}$, all of which are defined in the original formulation. ◇

## V. Higher-Order Parametric Differentiation

Another layer of spectral analysis is acquired by studying the parametric derivatives of the resolvent with respect to the spectral parameter $\lambda$. These derivatives encode fine geometric information on the spectral distribution and provide access to generalized Green-type functions and higher-order eigenmode couplings.

Starting from the Neumann series expansion,

$$ (\mathbb{I} - T(\lambda))^{-1} = \sum_{n=0}^{\infty} T(\lambda)^n, \quad (195) $$

its first derivative reads

$$ \frac{d}{d\lambda} (\mathbb{I} - T(\lambda))^{-1} = \sum_{n=1}^{\infty} \sum_{k=0}^{n-1} T(\lambda)^k \frac{dT(\lambda)}{d\lambda} T(\lambda)^{n-1-k}. \quad (196) $$

Each term of the double sum admits a kernel with $(n+1)$ spatial integrals. Specifically, for the action on $\vec{u}^{(0)}$,

$$ \left[ \frac{d}{d\lambda} (\mathbb{I} - T(\lambda))^{-1} [\vec{u}^{(0)}] \right] (\vec{X}) = \sum_{n=1}^{\infty} \sum_{k=0}^{n-1} \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_{n+1} \\ \times K_{T,\lambda}(\vec{X}, \vec{Y}_1) \cdots \frac{\partial K_{T,\lambda}}{\partial \lambda}(\vec{Y}_k, \vec{Y}_{k+1}) \cdots K_{T,\lambda}(\vec{Y}_n, \vec{Y}_{n+1}) \vec{u}^{(0)}(\vec{Y}_{n+1}). \quad (197) $$
Higher-order parametric derivatives are generated through repeated differentiation. For the p-th derivative,

$$ \frac{d^p}{d\lambda^p} (\mathbb{I} - T(\lambda))^{-1} = p! \sum_{\substack{n=p \\ k_1+\cdots+k_p=n-p}}^{\infty} \sum_{\substack{j=0 \\ k_j \ge 0}} T(\lambda)^{k_1} \frac{dT(\lambda)}{d\lambda} T(\lambda)^{k_2} \frac{dT(\lambda)}{d\lambda} \cdots \frac{dT(\lambda)}{d\lambda} T(\lambda)^{k_p}. \quad (198) $$

Each derivative term gives rise to a hierarchical chain of nested integrals with $(n+p)$ integration variables, encompassing multipolar spectral information and sensitivity of the solution to variations in $\lambda$.

In Fourier variables, the first derivative can be expressed as

$$ \frac{d}{d\lambda}[\hat{\cdot}](\xi_0) = \sum_{n=1}^{\infty} \sum_{k=0}^{n-1} \int_{\mathbb{R}^3} d\xi_1 \cdots \int_{\mathbb{R}^3} d\xi_{n+1} \\ \times \widehat{K_{T,\lambda}}(\xi_0 - \xi_1) \cdots \frac{\partial \widehat{K_{T,\lambda}}}{\partial \lambda}(\xi_k - \xi_{k+1}) \cdots \widehat{K_{T,\lambda}}(\xi_n - \xi_{n+1}) \widehat{u}^{(0)}(\xi_{n+1}), \quad (199) $$

$$ \text{where } [\cdot] \stackrel{\text{eqv}}{=} \left[ (\mathbb{I} - T(\lambda))^{-1} [\vec{u}^{(0)}] \right]. $$

This parametric derivative expansion provides access to the sensitivity of the spectral structure with respect to the parameter $\lambda$ and plays a fundamental role in perturbative spectral analysis, resonance unfolding, and the construction of higher-order Green operators.

## VI. Vue d'Ensemble

The parametric derivative structure described above admits a reformulation in terms of Volterra integrals [85], contour projectors, and asymptotic kernel expansions. This allows one to relate the resolvent $(\mathbb{I} - T(\lambda))^{-1}$ not only to its differential response with respect to $\lambda$, but also to its complete integral and spectral counterparts.

### A. Volterra-type nested integral expansion

Starting from the iterated kernel cascade and isolating the contribution of $K_{P,\lambda}$ or $K_{M,\lambda}$ at each iteration, one attains

$$ \vec{u}(\vec{X}) = \vec{u}(\vec{X})|_{n=0} + \sum_{n=1}^{\infty} \int_{0}^{\infty} ds_{1} \cdots \int_{0}^{\infty} ds_{n} \int_{\mathbb{R}^{3}} d\vec{Z}_{1} \cdots \int_{\mathbb{R}^{3}} d\vec{Z}_{n} \\ \times \left( \prod_{j=1}^{n} e^{-\lambda s_{j}} J(s_{j}; \vec{Z}_{j-1})^{-1} \right) K_{n}(\vec{Z}_{0}, \vec{Z}_{1}, \dots, \vec{Z}_{n}) \vec{u}(\vec{Z}_{n}), \quad (200) $$

where $\vec{Z}_0 = \vec{X}$ and $K_n$ denotes the ordered product of kernels selected among $K_{P,\lambda}$ and $K_{M,\lambda}$ at each step. This integral chain is the full Neumann expansion of the resolvent in real space.

### B. Projector-type contour formulation

The same series can be recast through the contour integral of the resolvent, where the spectral projector associated with a discrete eigenvalue $\lambda_*$ is fixed by

$$ \mathbb{P}_{\lambda_*} = \frac{1}{2\pi i} \oint_{\Gamma_c} (\mathbb{I} - T(\zeta))^{-1} d\zeta, \quad (201) $$

and, after expanding $(\mathbb{I} - T(\zeta))^{-1}$, one collects

$$ \mathbb{P}_{\lambda_*(\vec{X}, \vec{Z})} = \frac{1}{2\pi i} \oint_{\Gamma_c} \sum_{n=0}^{\infty} \int_{0}^{\infty} ds_1 \cdots \int_{0}^{\infty} ds_n \int_{\mathbb{R}^3} d\vec{Z}_1 \cdots d\vec{Z}_{n-1} \\ \times \left( \prod_{j=1}^{n} e^{-\zeta s_j} J(s_j; \vec{Z}_{j-1})^{-1} \right) K_n(\vec{Z}_0, \vec{Z}_1, \dots, \vec{Z}_n) d\zeta, \quad (202) $$
with $\vec{Z}_0 = \vec{X}$ and $\vec{Z}_n = \vec{Z}$.

This contour formulation isolates the spectral residues linked to $\lambda_*$ and reorganizes the integral hierarchy into a spectral series adapted to eigenfunction reconstruction and analytic continuation.

### C. Large-$\Re\lambda$ asymptotics of the iterated kernels

For $\Re\lambda \gg 1$, the integrals localise near $s_j = 0$, and the kernel structure can be approximated using the Taylor expansions,

$$Y(s; \vec{W}) = \vec{W} + sv(\vec{W}) + \frac{s^2}{2}(\partial v)(\vec{W})[v(\vec{W})] + \dots, \quad (203)$$

$$J(s; \vec{W})^{-1} = 1 - s(\nabla \cdot v)(\vec{W}) + \frac{s^2}{2}[(\nabla \cdot v)^2 - \text{tr}[(\nabla v)^2]](\vec{W}) + \dots, \quad (204)$$

and the Laplace moments,

$$\int_0^\infty e^{-\lambda s} s^k ds = \frac{\Gamma(k+1)}{\lambda^{k+1}} = \frac{k!}{\lambda^{k+1}}, \quad \Re\lambda > 0. \quad (205)$$

One then derives the asymptotic expansion,

$$K_{P,\lambda}^{(n)}(\vec{X}, \vec{Z}) \sim \sum_{m=0}^{\infty} \lambda^{-n-m} A_{n,m}(\vec{X}, \vec{Z}), \quad (206)$$

where each $A_{n,m}$ is a spatial kernel constructed from derivatives of $R$, $\vec{U}$, and $v$. For $n = 1$, the first two terms are

$$K_{P,\lambda}^{(1)}(\vec{X}, \vec{Z}) = \lambda^{-1} \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \nabla_{\vec{Z}} \cdot (R(\vec{Z})\delta(\vec{X} - \vec{Z})) \right\} \\ + \lambda^{-2} \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} (-v(\vec{X}) \cdot \nabla_{\vec{Z}}) \nabla_{\vec{Z}} \cdot (R(\vec{Z})\delta(\vec{X} - \vec{Z})) \right\} + \mathcal{O}(\lambda^{-3}). \quad (207)$$

This asymptotic formulation systematically reconstructs the non-local kernel structure in terms of local differential operators and provides a controlled approximation scheme in the regime $\Re\lambda \gg 1$.

## VII. Pseudo-differential Symbolics and Microlocal Propagation

The integral and spectral constructions admit a complementary formulation in pseudo-differential (microlocal) terms. The objective here is to express the reduced resolvent $(\mathbb{I} - T(\lambda))^{-1}$ and the iterated kernels $K_{T,\lambda}^{(n)}$ through Kohn-Nirenberg symbols [42] and to relate the nested similarity-time integrals to principal symbols obtained by partial Fourier transform in the relative spatial variable.

**Kohn-Nirenberg symbol.** For the kernel $K_{T,\lambda}(\vec{X}, \vec{Z})$, the symbol $t_\lambda(\vec{X}, \Xi)$ is introduced through the partial Fourier transform in the relative variable,

$$t_{\lambda}(\vec{X}, \Xi) = \int_{\mathbb{R}^3} e^{-i\Xi \cdot (\vec{X} - \vec{Z})} K_{T,\lambda}(\vec{X}, \vec{Z}) d\vec{Z}. \quad (208)$$

For any symbol $a(\vec{X}, \Xi)$, we associate the integral operator acting on test functions $\psi$ by the Kohn-Nirenberg quantization formula,

$$(\mathfrak{a}\psi)(\vec{X}) = \frac{1}{(2\pi)^3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} e^{i(\vec{X}-\vec{Y})\cdot\Xi} a(\vec{X}, \Xi)\psi(\vec{Y})d\vec{Y}d\Xi, \quad (209)$$

so that $T(\lambda)$ coincides with the operator given by the symbol $t_\lambda$ in (208).

**Symbol classes and composition.** Let $S_{1,0}^m$ denote the Kohn-Nirenberg symbol class: $a \in S_{1,0}^m$ if for all multi-indices $\alpha, \beta$ there exists $C_{\alpha\beta}$ with

$$|\partial_{\Xi}^{\alpha} \partial_{\vec{X}}^{\beta} a(\vec{X}, \Xi)| \le C_{\alpha\beta} (1 + |\Xi|)^{m-|\alpha|}. \quad (210)$$

The composition of the two quantized operators admits the standard asymptotic expansion,

$$a \sharp b \sim \sum_{\alpha \in \mathbb{N}^3} \frac{(i)^{-|\alpha|}}{\alpha!} \partial_{\Xi}^{\alpha} a(\vec{X}, \Xi) \partial_{\vec{X}}^{\alpha} b(\vec{X}, \Xi), \quad (211)$$
with the remainder estimates in $S_{1,0}^{m_a+m_b-N}$ after truncation at $|\alpha| < N$.

**Symbol-level Neumann series and parametrix.** Let $t_\lambda(\vec{X}, \Xi) \in S_{1,0}^m$ uniformly for $\Re\lambda > \lambda_0$, under the smallness condition on the associated operator norm, in accordance with the Neumann expansion of $(\mathbb{I} - T(\lambda))^{-1}$. Determine the inverse-symbol formal series,

$$s_{\lambda}(\vec{X}, \Xi) \sim \sum_{n=0}^{\infty} t_{\lambda}^{\#n}(\vec{X}, \Xi), \qquad (212)$$

where $t_\lambda^{\#0} = 1$ and $t_\lambda^{\#n} = \underbrace{t_\lambda^{\#} \cdots \# t_\lambda}_{n \text{ factors}}$. Then there is a remainder symbol $r_\lambda \in S^{-\infty}$ such that

$$(1 - t_\lambda)^{\#} s_\lambda = 1 + r_\lambda. \qquad (213)$$

Accordingly the reduced resolvent equals the operator associated with $s_\lambda$ up to a smoothing operator.

**Principal symbol of iterated integral kernels.** Let $K_{T,\lambda}^{(n)}$ signify the $n$-fold spatial composition kernel. Its principal symbol coincides with the principal part of $t_\lambda^{\#n}$. To exhibit the principal contribution of the pressure-type compositions, set for fixed $n$ the transported-density Fourier integral,

$$G_{\lambda}^{(n)}(\vec{X}, \Xi) = \int_{0}^{\infty} \cdots \int_{0}^{\infty} \int_{\mathbb{R}^{3n}} \left[ \prod_{j=1}^{n} e^{-\lambda s_j} J(s_j; \vec{Z}_{j-1})^{-1} R(\vec{Z}_j) \right] \times e^{-i\Xi \cdot (Y(s_n; \vec{Z}_{n-1}) - \vec{Z}_n)} d\vec{Z}_1 \cdots d\vec{Z}_n ds_1 \cdots ds_n, \quad (214)$$

with the convention $\vec{Z}_0 = \vec{X}$. The leading-order principal symbol of $K_{P,\lambda}^{(n)}$ satisfies

$$\sigma_{\text{prin}}(K_{P,\lambda}^{(n)})(\vec{X}, \Xi) = \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \right\} \cdot G_{\lambda}^{(n)}(\vec{X}, \Xi) + \text{lower-order terms in } \Xi. \quad (215)$$

Magnetic compositions $K_{M,\lambda}^{(n)}$ admit the same treatment with the divergence factors replaced by the corresponding magnetic integrands.

**Subprincipal hierarchy and recursive coefficients.** Expand the inverse symbol as an asymptotic series in decreasing orders of $\Xi$,

$$s_{\lambda}(\vec{X}, \Xi) \sim \sum_{j=0}^{\infty} s_{-j}(\vec{X}, \Xi), \qquad (216)$$

whose leading coefficients are given recursively by

$$s_0 = \frac{1}{1 - t_\lambda}, \qquad (217)$$

$$s_{-1} = -s_0 \frac{i}{2} \{t_\lambda, s_0\}, \qquad (218)$$

$$s_{-2} = -s_0 \left[ \frac{i}{2} \{t_\lambda, s_{-1}\} + \frac{1}{8} \{\{t_\lambda, s_0\}\} \right], \qquad (219)$$

where $\{a, b\} = \nabla\Xi a \cdot \nabla_X b - \nabla_X a \cdot \nabla\Xi b$ is the Poisson bracket, and $\{\cdot, \cdot\}$ the iterated bracket of order two.

**Microlocal propagation and wavefront set.** If $(\mathbb{I} - T(\lambda))\vec{u} = \vec{\varphi}$ and $t_\lambda$ is elliptic on a conic region $\Gamma \subset T^*\mathbb{R}^3$, then

$$\mathrm{WF}(\vec{u}) \cap \Gamma = \mathrm{WF}(\vec{\varphi}) \cap \Gamma, \qquad (220)$$

where $\mathrm{WF}(\vec{u})$ is the wavefront set, for the spatial-frequency localization of singularities.
When ellipticity fails, the singularities propagate along the bicharacteristic flow of the Hamiltonian vector field,

$$ \mathcal{H}_{t_{\lambda}} = \nabla_{\Xi} t_{\lambda} \cdot \nabla_{\vec{X}} - \nabla_{\vec{X}} t_{\lambda} \cdot \nabla_{\Xi}, \qquad (221) $$

which refines the description of resonance manifolds and microlocal spectral strata.

**Trace and determinant asymptotics.** Under the above symbol-calculus hypotheses the trace of $T(\lambda)^n$ has the semiclassical-type expansion,

$$ \mathrm{tr}(T(\lambda)^n) \sim (2\pi)^{-3} \int_{\mathbb{R}^3} \int_{\mathbb{R}^3} \mathrm{tr}_{\mathbb{R}^d} (t_{\lambda}^{\sharp n}(\vec{X}, \Xi)) d\Xi d\vec{X} + r_{n,\lambda}, \quad r_{n,\lambda} = \mathcal{O}(\lambda^{-\epsilon}), \quad \epsilon > 0. \qquad (222) $$

The log-det expansion reads

$$ \log \det (\mathbb{I} - T(\lambda)) = - \sum_{n=1}^{\infty} \frac{1}{n} \mathrm{tr}(T(\lambda)^n), \qquad (223) $$

and the leading contributions to the regularized determinant are computable from the integrals of the symbols $t_{\lambda}^{\sharp n}$ in (222).

**Large $\Re\lambda$ localisation.** For $\Re\lambda \gg 1$, the Laplace factors $e^{-\lambda s_j}$ localize the similarity-time integrals near $s_j = 0$. An explicit expansion in inverse powers of $\lambda$ is obtained by expanding transported quantities and evaluating Laplace moments,

$$ t_{\lambda}(\vec{X}, \Xi) \sim \lambda^{-1} t^{(0)}(\vec{X}, \Xi) + \lambda^{-2} t^{(1)}(\vec{X}, \Xi) + \dots, \qquad (224) $$

where $t^{(0)}$ results from replacing each transported factor by its instantaneous value at $s = 0$ and inserting the leading Laplace moment $\lambda^{-1}$.

**Sufficient conditions for validity.** The preceding pseudo-differential description is valid under the following checkable hypotheses;

(i) $K_{T,\lambda}(\vec{X}, \vec{Z})$ is $\mathcal{C}^{\infty}$ off the diagonal $\vec{X} = \vec{Z}$ with uniform derivative bounds so that $t_{\lambda} \in S_{1,0}^m$, for some $m < 0$ and for $\Re\lambda > \lambda_0$;

(ii) the operator associated with $t_{\lambda}$ is bounded on the chosen function space and its operator norm is $< 1$ (this reproduces the smallness hypothesis used earlier);

(iii) for each $n$, the transported-density Fourier integral $G_{\lambda}^{(n)}(\vec{X}, \Xi)$ in (214) qualifies as a symbol in $\Xi$ with the stated growth.

### 8.1.2.3. Radial Reductions, Spherical Decoupling and Mellin-Bessel via Evaluation

We collect a concise set of integral-calculus identities and reductions that are tailored to the spatial integrals appearing in the $T(\lambda)$ expansions.

We first introduce a local shorthand for the precise composite kernel appearing in the definition of $T(\lambda)$,

$$ K_{T,\lambda}(\vec{X}, \vec{Z}) = (\lambda \mathbb{I} - \mathcal{D}_{\mathrm{loc}}^{\mathrm{mom}})^{-1} \left( K_{P,\lambda}(\vec{X}, \vec{Z}) + K_{M,\lambda}(\vec{X}, \vec{Z}) + \tilde{K}_M(\vec{X}, \vec{Z}) \circ (a_{\mathrm{ind}}(\lambda)^{-1} w) \right). \qquad (225) $$

By construction one has $T(\lambda)[\vec{u}](\vec{X}) = \int_{\mathbb{R}^3} K_{T,\lambda}(\vec{X}, \vec{Z}) \vec{u}(\vec{Z}) d\vec{Z}$.

When $K_{T,\lambda}$ is isotropic, or when it depends on $\vec{X}, \vec{Z}$ only through the radii $r = |\vec{X}|, s = |\vec{Z}|$, we write the scalar radial kernel $k_{T,\lambda}(r, s)$ so that

$$ K_{T,\lambda}(\vec{X}, \vec{Z}) = k_{T,\lambda}(r, s) \mathbb{I}_{\mathbb{R}^d}, \quad r = |\vec{X}|, s = |\vec{Z}|. \qquad (226) $$

Switching to spherical coordinates, $\vec{X} = r\hat{X}, \vec{Z} = s\hat{Z}$ and $d\vec{Y} = \rho^2 d\rho d\Omega(\hat{Y})$, one marks out the standard radial reduction for $m$-fold spatial convolution of an isotropic kernel,

$$ \begin{aligned} \left[T(\lambda)^m\right]_{\mathrm{diag}}(r, s) &= \int_{\mathbb{R}^3} d\vec{Y}_1 \cdots \int_{\mathbb{R}^3} d\vec{Y}_{m-1} k_{T,\lambda}(r, \rho_1) k_{T,\lambda}(\rho_1, \rho_2) \cdots k_{T,\lambda}(\rho_{m-1}, s) \\ &= (4\pi)^{m-1} \int_0^\infty \rho_1^2 d\rho_1 \int_0^\infty \rho_2^2 d\rho_2 \cdots \int_0^\infty \rho_{m-1}^2 d\rho_{m-1} \\ &\qquad \times k_{T,\lambda}(r, \rho_1) k_{T,\lambda}(\rho_1, \rho_2) \cdots k_{T,\lambda}(\rho_{m-1}, s), \end{aligned} \qquad (227) $$
where the factor $(4\pi)^{m-1}$ arises from the angular integrals $\int_{S^2} d\Omega(\hat{Y}) = 4\pi$ performed at each intermediate convolution node.

The diagonal traces reduce in consequence to purely radial integrals. For instance, for isotropic $k_{T,\lambda}$,

$$ \mathrm{tr}(T(\lambda)^m) = (4\pi)^{m-1} \int_0^\infty r^2 dr \int_0^\infty \rho_1^2 d\rho_1 \cdots \int_0^\infty \rho_{m-1}^2 d\rho_{m-1} \mathrm{tr}_{\mathbb{R}^d} (k_{T,\lambda}(r, \rho_1) \cdots k_{T,\lambda}(\rho_{m-1}, r)). \quad (228) $$

When the kernel is not isotropic but admits an angular expansion in spherical harmonics, the convolution decouples by angular momentum channels. Write the standard spherical-harmonic decomposition,

$$ K_{T,\lambda}(\vec{X}, \vec{Z}) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} k_{\ell,\lambda}(r, s) Y_{\ell m}(\hat{X}) \overline{Y_{\ell m}(\hat{Z})}, \quad (229) $$

where $Y_{\ell m}$ are the unit-normalized spherical harmonics on $S^2$ and $r = |\vec{X}|, s = |\vec{Z}|$. Orthogonality of $Y_{\ell m}$ implies that the $\ell$-channels decouple under convolution, and one obtains for each $\ell$ an independent radial cascade,

$$ [T(\lambda)^m]_\ell(r, s) = (4\pi)^{m-1} \int_0^\infty \rho_1^2 d\rho_1 \cdots \int_0^\infty \rho_{m-1}^2 d\rho_{m-1} k_{\ell,\lambda}(r, \rho_1) \cdots k_{\ell,\lambda}(\rho_{m-1}, s). \quad (230) $$

A practical evaluation strategy for kernels built from radial Green's functions uses repeated application of the radial resolvent $G_k(r, s; \lambda)$. Suppose $k_{T,\lambda}$ factorizes as a composition of a radial Green's function and a compact multiplier, e.g.

$$ k_{T,\lambda}(r, s) = G_k(r, s; \lambda)m(s), \quad (231) $$

with a local multiplier $m(\rho)$. Then the $m$-fold convolution reduces to repeated radial action of $G_k$ interleaved with $m(\cdot)$ insertions,

$$ k_{T,\lambda}^{(m)}(r, s) = \int_{0}^{\infty} d\rho_{1} \cdots \int_{0}^{\infty} d\rho_{m-1} G_{k}(r, \rho_{1}; \lambda) m(\rho_{1}) G_{k}(\rho_{1}, \rho_{2}; \lambda) \cdots m(\rho_{m-1}) G_{k}(\rho_{m-1}, s; \lambda). \quad (232) $$

When closed forms for $G_k$ in terms of modified Bessel functions $I_\zeta, K_\zeta$ are available, each inner radial integral may be expressed as Mellin-Bessel type integrals [54] which are amenable to contour deformation or to standard integral tables.

Integration by parts is often effective to move radial derivatives off singular weight factors appearing in Green's functions. An identity valid for sufficiently decaying functions at 0 and $\infty$ is

$$ \int_0^\infty \rho^2 \partial_\rho (\Phi(\rho)) \Psi(\rho) d\rho = - \int_0^\infty \Phi(\rho) \partial_\rho (\rho^2 \Psi(\rho)) d\rho, \quad [\rho^2 \Phi(\rho) \Psi(\rho)]_0^\infty = 0, \quad (233) $$

which is used to transfer derivatives appearing inside $G_k$ onto multipliers $m(\rho)$ or onto explicit polynomial weights.

As an illustrative explicit consequence we record the diagonal formulas for the first two trace terms in the isotropic case,

$$ \mathrm{tr}(T(\lambda)) = 4\pi \int_0^\infty r^2 \mathrm{tr}_{\mathbb{R}^d} (k_{T,\lambda}(r, r)) dr, \quad (234) $$

$$ \mathrm{tr}(T(\lambda)^2) = 4\pi \int_0^\infty r^2 dr \int_0^\infty \rho^2 d\rho \mathrm{tr}_{\mathbb{R}^d} (k_{T,\lambda}(r, \rho) k_{T,\lambda}(\rho, r)). \quad (235) $$

For convergence control of the radial chains one may adopt the Hilbert-Schmidt criterion in radial variables. For an isotropic kernel $k_{T,\lambda}$ one has

$$ \|T(\lambda)\|_{\mathrm{HS}}^2 = (4\pi)^2 \int_0^\infty \int_0^\infty r^2 s^2 \|k_{T,\lambda}(r, s)\|_{\mathrm{F}}^2 dr ds, \quad (236) $$

where $\| \cdot \|_{\mathrm{HS}}$ is the Hilbert-Schmidt norm, and $\| \cdot \|_{\mathrm{F}}$ the Frobenius (finite-dimensional) norm [24] on component matrices; absolute summability of the nested integrals follows when $\|T(\lambda)\|_{\mathrm{HS}} < \infty$ and $\|T(\lambda)\|_{\mathrm{op}} < 1$ or under analogous Schatten-class hypotheses.
The identities above provide a practical calculus: reduce spatial convolutions to radial integrals when isotropy permits; otherwise diagonalize angular structure with spherical harmonics and process each ℓ-channel separately; when Green's functions $G_k$ are present evaluate inner integrals by Mellin-Bessel techniques and take advantage of the radial integration by parts to remove singular derivatives.

### 8.1.2.4. Decomposition of the Fredholm Trace-log

At this stage we can provide an explicit integral expansion of the Fredholm trace-log, obtained by substituting the similarity-time resolvent into the kernel-chain trace formula. The result is a multi-integral expression that links the similarity-time resolvent weights $\int_0^\infty e^{-\lambda s} J(s; \cdot)^{-1}(\cdot) ds$ to the secular determinant. Let

$$ \log \det (\mathbb{I} - T(\lambda)) = - \sum_{m=1}^{\infty} \frac{1}{m} \int_{\mathbb{R}^3} d\vec{X}_1 \cdots \int_{\mathbb{R}^3} d\vec{X}_m \operatorname{tr}_{\mathbb{R}^d} [K_{T,\lambda}(\vec{X}_1, \vec{X}_2) \cdots K_{T,\lambda}(\vec{X}_m, \vec{X}_1)] . \quad (237) $$

Substitute the explicit similarity-time integral form of each constituent kernel factor $K_{T,\lambda}$. Recall that $K_{T,\lambda}$ is built from $K_{P,\lambda}$, $K_{M,\lambda}$ and the $\tilde{K}_M$-parts, each containing the resolvent weight $\int_0^\infty e^{-\lambda s} J^{-1}(s) ds$. For every $m \ge 1$,

$$ \begin{aligned} & \operatorname{tr}_{\mathbb{R}^d} [K_{T,\lambda}(\vec{X}_1, \vec{X}_2) \cdots K_{T,\lambda}(\vec{X}_m, \vec{X}_1)] \\ &= \sum_{\sigma \in \{P, M, \tilde{M}\}^m} \int_0^\infty ds_1 \cdots \int_0^\infty ds_m \int_{\mathbb{R}^3} d\vec{Z}_1 \cdots \int_{\mathbb{R}^3} d\vec{Z}_m \prod_{q=1}^m e^{-\lambda s_q} J(s_q; \vec{X}_q)^{-1} \\ & \qquad \times \operatorname{tr}_{\mathbb{R}^d} \left[ \prod_{q=1}^m \mathcal{D}_{\sigma_q}(\vec{X}_q, \vec{Z}_q; s_q) \right] \prod_{q=1}^m \delta(Y(s_q; \vec{X}_q) - \vec{Z}_q). \end{aligned} \quad (238) $$

When $\sigma_q = P$ (pressure-type factor) the local differential density appearing in (238) is the pressure integrand arising from $K_{P,\lambda}$, namely

$$ \mathcal{D}_P(\vec{X}, \vec{Z}; s) = \nabla_{\vec{X}} \left\{ \frac{P'(R(\vec{X}))}{R(\vec{X})} \nabla_{\vec{Z}} \cdot (R(\vec{Z})(\cdot)) \right\} \Bigg|_{(\cdot)=\delta(Y(s;\vec{X})-\vec{Z})}, \quad (239) $$

and when $\sigma_q = M$ or $\tilde{M}$ the density $\mathcal{D}_{\sigma_q}$ is the magnetic integrand for $K_{M,\lambda}$ and $\tilde{K}_M$, with the same resolvent weight under the s-integral. Each $\mathcal{D}_{\sigma_q}$ is a composition of gradients, curls and pointwise multipliers built from $R$, $U$, $H$ and their derivatives.

Combining (237) and (238) we earn the full secular multi-integral,

$$ \begin{aligned} \log \det (\mathbb{I} - T(\lambda)) = & - \sum_{m=1}^{\infty} \frac{1}{m} \sum_{\sigma \in \{P, M, \tilde{M}\}^m} \int_{(\mathbb{R}^3)^m} d\vec{X}_1 \cdots d\vec{X}_m \int_{[0,\infty)^m} ds_1 \cdots ds_m \\ & \times \int_{(\mathbb{R}^3)^m} d\vec{Z}_1 \cdots d\vec{Z}_m \prod_{q=1}^m e^{-\lambda s_q} J(s_q; \vec{X}_q)^{-1} \\ & \times \operatorname{tr}_{\mathbb{R}^d} \left[ \prod_{q=1}^m \mathcal{D}_{\sigma_q}(\vec{X}_q, \vec{Z}_q; s_q) \right] \prod_{q=1}^m \delta(Y(s_q; \vec{X}_q) - \vec{Z}_q). \end{aligned} \quad (240) $$

Eq. (240) displays the m-fold spatial integration, the m-fold similarity-time integration, the distributional concatenation along the characteristics $Y(s; \cdot)$ via the Dirac deltas, and the finite-dimensional trace coupling the matrix-valued integrands.

Which is why the spectral condition $\det(\mathbb{I} - T(\lambda)) = 0$ is equivalent to the vanishing of (240). Truncating the outer sum at $m = M$ it is possible to determine a finite secular approximation that can be used for asymptotic root-finding or for numerical evaluation of spectral moments.
### 8.1.3. Two Supporting Propositions

#### 8.1.3.1. I. On the Magnetic Divergence-Free Condition

**Proposition 8.1** (Propagation of the magnetic constraint in similarity variables). Let $(\vec{U}, \vec{H}, R)$ solve Eqq. (88) (89) and set $v = \nabla \cdot \vec{H}$. Then

$$ -\frac{1}{2}v - \frac{1}{2}(\vec{X} \cdot \nabla)v = \nu \Delta v, \text{ in the stationary similarity regime.} \qquad (241) $$

If $v \in L_{\omega}^2(\mathbb{R}^3)$ with $\omega(\vec{X}) = (1+|\vec{X}|^2)^{-\beta}$, $\beta > \frac{3}{2}$, one concludes that $v = 0$, and any finite-similarity-energy solution satisfies $\nabla \cdot \vec{H} = 0$.

*Proof.* Apply $\nabla \cdot$ to (89). The curl term vanishes, $\nabla \cdot (\nabla \times (\vec{U} \times \vec{H})) \stackrel{\text{eqv}}{=} 0$, which immediately gives the displayed equation for $v$. Test against $v\omega$ and integrate by parts,

$$ \nu \int |\nabla v|^2 \omega = \frac{1}{4} \int v^2 (1 + \vec{X} \cdot \nabla \log \omega) \omega. \qquad (242) $$

For $\omega(\vec{X}) = (1 + |\vec{X}|^2)^{-\beta}$ one has $\vec{X} \cdot \nabla \log \omega = -\frac{2\beta|\vec{X}|^2}{1+|\vec{X}|^2} \le 0$, so the right-hand side $\le \frac{1}{4} \int v^2 \omega$. Thus $\nabla v \stackrel{\text{eqv}}{=} 0$ in $L_{\omega}^2$ and the equation forces $v$ to be constant; the weight excludes non-zero constants, which shows that $v = 0$. $\square$

#### 8.1.3.2. II. On the Maxwell-Transport Cancellation

**Proposition 8.2** (Subtle cancellation aspects). Assume $\nabla \cdot \vec{H} = 0$ and sufficient decay at infinity. Then

$$ \int_{\mathbb{R}^3} ((\nabla \times \vec{H}) \times \vec{H}) \cdot \vec{U} - \int_{\mathbb{R}^3} (\vec{U} \cdot \nabla) \vec{H} \cdot \vec{H} = \frac{1}{2} \int_{\mathbb{R}^3} (\nabla \cdot \vec{U}) |\vec{H}|^2. \qquad (243) $$

In the incompressible limit $\nabla \cdot \vec{U} = 0$ the two non-linear magnetic couplings cancel exactly in the energy balance.

*Proof.* Use $(\nabla \times \vec{H}) \times \vec{H} = (\vec{H} \cdot \nabla)\vec{H} - \nabla \frac{|\vec{H}|^2}{2}$ and $(\vec{U} \cdot \nabla)\frac{|\vec{H}|^2}{2} = (\vec{U} \cdot \nabla)\vec{H} \cdot \vec{H}$. Integration by parts (boundary terms vanish by decay) gives the identity. $\square$

### 9. Leray-Type Self-Similar Ansatz

We adopt a backward self-similar ansatz to analyze finite-time singularities in the compressible MHD system, suitably generalized to account for compressibility and magnetic fields.
Let $T > 0$ be a putative blow-up time. Consider the similarity variables,

$$ \vec{X} = \frac{\vec{x}}{(T-t)^{\beta}}, \quad \tau = -\log(T-t), \qquad (244) $$

and seek self-similar profiles of the form

$$ \vec{u}(\vec{x}, t) = (T - t)^{-\alpha} \vec{U}(\vec{x}), \qquad (245a) $$

$$ \rho(\vec{x}, t) = (T - t)^{-\gamma} R(\vec{x}), \qquad (245b) $$

$$ \vec{B}(\vec{x}, t) = (T - t)^{-\alpha} \vec{B}(\vec{x}), \qquad (245c) $$

$$ p(\vec{x}, t) = (T - t)^{-\gamma_p} p(\vec{x}). \qquad (245d) $$

Alternatively, we consider solutions undergoing finite-time singularities at $t = t^* = 0$. Following the classical approach of Leray, we introduce a backward self-similar ansatz for the compressible
MHD system,

$$ \vec{x} = t^{\beta} \vec{X}, \qquad (246a) $$

$$ \vec{u}(\vec{x}, t) = t^{-\alpha} \vec{U}(\vec{X}), \qquad (246b) $$

$$ \rho(\vec{x}, t) = t^{-\gamma} R(\vec{X}), \qquad (246c) $$

$$ \vec{B}(\vec{x}, t) = t^{-\alpha} \vec{B}(\vec{X}), \qquad (246d) $$

$$ p(\vec{x}, t) = t^{-\gamma_p} p(\vec{X}), \qquad (246e) $$

where $\vec{X} = \frac{\vec{x}}{t^{\beta}}$ are similarity coordinates, and the exponents $\alpha, \beta$, and $\gamma$ control, respectively, the time decay or growth of the velocity, spatial scaling, and density blow-up. This ansatz prescribes homogeneity of the unknowns under parabolic rescaling, preserving the structure of the inviscid compressible MHD system.

Plugging (245) into (3), one writes

$$ \nabla \cdot (R\vec{U}) + \beta \nabla \cdot (\vec{X}R) + \gamma R = 0, \qquad (247a) $$

$$ R\left(\left(\vec{U} + \beta\vec{X}\right) \cdot \nabla\vec{U} + \alpha\vec{U}\right) + \nabla p = \mu \Delta\vec{U} + \mu_{\text{bulk}} \nabla\left(\nabla \cdot \vec{U}\right) + \left(\nabla \times \vec{B}\right) \times \vec{B}, \qquad (247b) $$

$$ (\vec{U} + \beta \vec{X}) \cdot \nabla \vec{B} + \alpha \vec{B} - (\vec{B} \cdot \nabla) \vec{U} + (\nabla \cdot \vec{U}) \vec{B} = \nu \Delta \vec{B}, \qquad (247c) $$

$$ \nabla \cdot \vec{B} = 0, \qquad (247d) $$

$$ p = AR^{\gamma}. \qquad (247e) $$

This ansatz induces a natural homogeneity structure on the compressible MHD system, rendering it invariant under the continuous rescaling group,

$$ G_{\lambda}: \left(\vec{x}, t, \vec{u}, \vec{B}, \rho, p\right) \mapsto \left(\lambda\vec{x}, \lambda^{\frac{1}{\beta}}t, \lambda^{-\alpha}\vec{u}, \lambda^{-\alpha}\vec{B}, \lambda^{-\gamma}\rho, \lambda^{-\gamma_p}p\right). \qquad (248) $$

This ansatz is invariant under a continuous scaling group acting on space-time (Section 11.1).

## 9.1. Geometric Interpretation via Lie Derivatives

The self-similar compressible MHD system admits a natural reformulation in the language of differential geometry. The continuity equation, responsible for mass conservation, exhibits an elegant expression in terms of Lie derivatives, which govern how densities evolve under flows.

The classical continuity equation in Eulerian form,

$$ \frac{\partial R}{\partial t} + \nabla \cdot (R \vec{u}) = 0, \qquad (249) $$

describes the conservation of mass for a density field $R(\vec{x}, t)$ transported by a velocity field $\vec{u}(\vec{x}, t)$. Under the self-similar ansatz, where variables scale with time as

$$ \vec{u}(\vec{x}, t) = \frac{1}{(-t)^{\alpha}} \vec{U} \left( \frac{\vec{x}}{(-t)^{\beta}} \right), \quad R(\vec{x}, t) = \frac{1}{(-t)^{\gamma}} R \left( \frac{\vec{x}}{(-t)^{\beta}} \right), \qquad (250) $$

the continuity equation transforms into an ordinary differential equation in the radial coordinate $r = \frac{|\vec{x}|}{(-t)^{\beta}}$.

Let

$$ \vec{v}(r) = (U(r) - \beta r)\hat{r} \qquad (251) $$

be a vector field that governs the effective advection in similarity coordinates, where the term $-\beta r$ accounts for the contraction induced by scaling.

The radial continuity equation in similarity variables, formally derived in Section 10.2, takes the schematic form

$$ \frac{d}{dr}(RU) + \frac{2}{r}RU + \beta r \frac{dR}{dr} + \gamma R = 0. \qquad (252) $$

This can be interpreted geometrically using the Lie derivative of the density $R$ along the vector field $\vec{v}$,

$$ \mathcal{L}_{\vec{v}} R + (\nabla \cdot \vec{v}) R = -\gamma R. \qquad (253) $$
This expression illustrates how self-similar transport consists of a combination of directional advection by the velocity profile $U(r)$, radial scaling drift induced by $\beta$, and a geometric source term weighted by $\gamma$.

The divergence term $\nabla \cdot \vec{v}$ contributes a geometric stretching or compression, analogous to a volume form deformation. In spherical coordinates one has

$$ \nabla \cdot \vec{v} = \frac{dU}{dr} + \frac{2}{r}U(r) - 3\beta, \qquad (254) $$

reflecting the balance between local velocity gradients and global similarity contraction.

This formulation provides a deeper geometric insight into the mass conservation mechanism in self-similar compressible flows and serves as a bridge to more general tensorial reformulations of the full MHD system, including the momentum and induction equations. The Lie derivative formulation is valid provided the fields are smooth and $R > 0$; it geometrically codifies the conservation of the volume form $RdV$ under the flow generated by $\vec{v}$.

## 9.2. Energetic and Entropic Scalings

To identify admissible values for $(\alpha, \beta, \gamma)$, we insert the ansatz into the compressible MHD equations and demand balance of leading-order terms.

(1) From the continuity scaling one obtains

$$ \partial_t \rho \sim t^{-\gamma-1}, \quad \nabla \cdot (\rho \vec{u}) \sim t^{-\gamma-\alpha-\beta}, \qquad (255) $$

hence equating exponents gives

$$ -\gamma - 1 = -\gamma - \alpha - \beta \implies \alpha + \beta = 1. \qquad (256) $$

(2) For the momentum equation one must choose which terms are intended to balance. The inertial time-derivative term scales as

$$ \rho \partial_t \vec{u} \sim t^{-\gamma-\alpha-1}, \qquad (257) $$

the convective inertial term as

$$ \rho(\vec{u} \cdot \nabla)\vec{u} \sim t^{-\gamma-2\alpha-\beta}, \qquad (258) $$

the pressure gradient as

$$ \nabla p \sim t^{-\gamma_p - \beta}, \qquad (259) $$

and the viscous term as

$$ \mu \Delta \vec{u} \sim t^{-\alpha-2\beta}. \qquad (260) $$

Because of that,

(i) if one enforces time-inertia $\sim$ pressure: $-\gamma - \alpha - 1 = -\gamma_p - \beta$, hence

$$ \gamma_p = \gamma + \alpha + 1 - \beta, \qquad (261) $$

(ii) if one enforces convective inertia $\sim$ pressure: $-\gamma - 2\alpha - \beta = -\gamma_p - \beta$, hence

$$ \gamma_p = \gamma + 2\alpha, \qquad (262) $$

(iii) if one enforces convective inertia $\sim$ viscosity: $-\gamma - 2\alpha - \beta = -\alpha - 2\beta$, hence

$$ \gamma + \alpha = \beta. \qquad (263) $$

Which balance is physically relevant depends on the regime under study: inertial vs. viscous dominated.

(3) The induction scaling, $-\alpha + \beta = -\alpha + \beta$, is consistent (trivially) provided the chosen balances on the momentum side are used; the formal identity of exponents for advection terms holds.

Thus, the system is consistent under the constraints,

$$ \alpha + \beta = 1, \quad \gamma_p = \gamma + \alpha + 1 - \beta, \qquad (264) $$

or, in the convective regime,

$$ \alpha + \beta = 1, \quad \gamma_p = \gamma + 2\alpha, \qquad (265) $$

or, in the viscous regime,

$$ \alpha + \beta = 1, \quad \gamma + \alpha = \beta. \qquad (266) $$
If one assumes a barotropic closure $p = A\rho^{\gamma_p}$, the pressure and density share compatible scaling whenever $\gamma_p = \alpha + \beta$, which follows in the special case $\gamma = 0$, i.e. scale-invariant density.

### 9.3. Entropy and Energy Scaling

The energy integral scales as

$$ \mathcal{E}[t] = \int_{\mathbb{R}^3} \left( \frac{1}{2} \rho |\vec{u}|^2 + \frac{1}{2\mu_0} |\vec{B}|^2 \right) dx \sim t^{-\gamma - 2\alpha + 3\beta}, \qquad (267) $$

while the total entropy behaves as

$$ S(t) = \int_{\mathbb{R}^3} \rho \log \rho dx \sim t^{-\gamma+3\beta}. \qquad (268) $$

These relations inform critical thresholds for anomalous dissipation and finite-time energy blow-up. The case $\gamma = 3\beta$ yields constant entropy, consistent with adiabatic collapse.

### 9.4. Self-Similar Reduction of the Full System

Under the self-similar change $\vec{X} = \frac{\vec{x}}{t^\beta}$ and the ansatz $\vec{u}(\vec{x}, t) = t^{-\alpha}\vec{U}(\vec{X})$, one has, when acting on the profile factors,

$$ \partial_t (t^{-\alpha} \Phi(\vec{X})) = t^{-\alpha-1} \left( -\alpha \Phi(\vec{X}) - \beta \vec{X} \cdot \nabla_{\vec{X}} \Phi(\vec{X}) \right), \qquad (269) $$

so that (as an operator on profile functions)

$$ \partial_t = -\frac{\beta}{t} \vec{X} \cdot \nabla_{\vec{X}} - \frac{\alpha}{t}, \quad \nabla_{\vec{x}} = t^{-\beta} \nabla_{\vec{X}}. \qquad (270) $$

With a Leray-type self-similar ansatz applied to the compressible MHD equations, the full system reduces to a time-independent formulation in the similarity variables $\vec{X} \in \mathbb{R}^3$.

(1) **Continuity | scalar equation for $R(\vec{X})$,**

$$ -\gamma R - \beta \vec{X} \cdot \nabla R + (\vec{U} \cdot \nabla) R + R \nabla \cdot \vec{U} = 0. \qquad (271) $$

This is a first-order transport equation with source, interpreted as a stationarity condition under the renormalized similarity flow derived from inserting the ansatz $\rho(\vec{x}, t) = t^{-\gamma}R(\vec{X})$ into the continuity equation and dividing by the common factor $t^{-\gamma-1}$. In the inviscid limit, this equation admits characteristics aligned with streamlines of $\vec{U} - \beta\vec{X}$, revealing the importance of the so-called “effective velocity” field.

(2) **Momentum | vector equation for $U(\vec{X})$,**

$$ R \left( -\alpha \vec{U} - \beta (\vec{X} \cdot \nabla) \vec{U} + (\vec{U} \cdot \nabla) \vec{U} \right) = \mu \Delta \vec{U} - \nabla \left( p + \frac{|\vec{B}|^2}{2\mu_0} \right) + \frac{1}{\mu_0} (\vec{B} \cdot \nabla) \vec{B}. \qquad (272) $$

This is a compressible generalization of the steady MHD momentum equation in similarity variables. The left-hand side represents the effective convective acceleration adjusted for similarity scaling, and the right-hand side combines viscous diffusion and magnetic stress.

To close the system, we prescribe a barotropic law for the pressure,

$$ p = K_c R^{\gamma_p}, \quad K_c > 0, \quad \gamma_p > 1, \qquad (273) $$

where $K_c = \frac{p_0}{R_0^{\gamma_p}}$ is a proportionality constant that depends on the choice of physical or reference scales.ª This leads to the total pressure

$$ P = p + \frac{|\vec{B}|^2}{2\mu_0} = K_c R^{\gamma_p} + \frac{|\vec{B}|^2}{2\mu_0}. \qquad (274) $$

Consequently, the momentum equation becomes

$$ R \left( -\alpha \vec{U} - \beta (\vec{X} \cdot \nabla) \vec{U} + (\vec{U} \cdot \nabla) \vec{U} \right) = \mu \Delta \vec{U} - \nabla P + \frac{1}{\mu_0} (\vec{B} \cdot \nabla) \vec{B}. \qquad (275) $$

<footnote>ª With $p_0$ and $R_0$ denoting characteristic pressure and density scales, respectively. In non-dimensionalized form, one may set $K = 1$ without loss of generality.</footnote>
The non-linear coupling between $R$, $\vec{U}$, and $\vec{B}$ through the pressure $P$ introduces significant complexity, remarkably in the compressible regime. The system may become overdetermined unless the scaling exponents and closure relations are chosen consistently.

$$ (3) \quad \textbf{Induction} \mid \textbf{vector equation for } \mathbf{B}(\vec{X}), $$

$$ -\alpha \vec{B} - \beta(\vec{X} \cdot \nabla)\vec{B} + (\vec{U} \cdot \nabla)\vec{B} = \nu \Delta \vec{B} + (\vec{B} \cdot \nabla)\vec{U} - \vec{B} \nabla \cdot \vec{U}. \quad (276) $$

This equation governs the evolution of the magnetic field profile $\vec{B}(\vec{X})$ and reflects the balance of advection, magnetic stretching and diffusion in similarity variables. The structure of this equation admits a Lie derivative interpretation,

$$ \mathcal{L}_{\vec{U}-\beta\vec{X}} \vec{B} + \alpha \vec{B} = \nu \Delta \vec{B}, \quad (277) $$

suggesting invariance under the effective flow.

$$ (4) \quad \textbf{Magnetic constraint (scalar) for } \vec{B}(\vec{X}), $$

$$ \nabla \cdot \vec{B} = 0, \quad (278) $$

implying compatibility with Maxwell's equations and the absence of magnetic monopoles.

This similarity system governs the structure of time-invariant profiles $(R, \vec{U}, \vec{B})$ in similarity variables, under the constraint that the full solution $(\rho, \vec{u}, \vec{B})$ exhibits singular behavior at $t = 0$ through scaling.

(5) **Induction in Curl Form.** If the magnetic field satisfies $\nabla \cdot \vec{B} = 0$, we may equivalently write the induction Eq. (276) in curl form,

$$ -\alpha \vec{B} - \beta(\vec{X} \cdot \nabla)\vec{B} = \nabla \times (\vec{U} \times \vec{B}) + \nu \Delta \vec{B}. \quad (279) $$

This expression reveals the geometric origin of magnetic field evolution: the generation of new field lines via the curl of the electromotive force $\vec{U} \times \vec{B}$, along with diffusion via resistivity; in other words, it emphasizes the interplay between the non-linear vortex stretching $\vec{U} \times \vec{B}$ and magnetic diffusion. It also aligns structurally with the Helmholtz decomposition and enables compact variational formulations.

## 9.5. Further Remarks

(1) **Incompressible Limit and Inviscid Collapse.** In the inviscid and incompressible limit,

$$ \nabla \cdot \vec{U} = 0, \quad \rho = \text{const}, \quad \mu = \nu = 0, \quad (280) $$

the system reduces to classical Leray-type similarity equations for incompressible MHD. The current model extends this framework by incorporating compressibility and diffusive operators, significantly enriching the phase space and blow-up scenarios.

(2) **On Inviscid Limit and Scale-Invariance.** In the limit $\mu = \nu = 0$, the system becomes formally scale-invariant under

$$ \vec{x} \mapsto \lambda \vec{x}, \quad t \mapsto \lambda^{\frac{1}{\beta}} t, \quad (281) $$

which preserves the ansatz structure. This invariance motivates the search for self-similar blow-up profiles as attractors or asymptotic states in the inviscid compressible MHD dynamics.

These self-similar solutions serve as a bridge between rigorous fluid dynamics, dynamical systems, and asymptotic analysis.

(3) **Similarity with Incompressible Systems.** Despite the fact that the above system generalizes the classical incompressible Navier-Stokes MHD equations, it preserves structural analogies:

(i) transport-diffusion mechanisms (via $\vec{U}$ and $\Delta$),

(ii) stretching and amplification of $\vec{B}$ by $\nabla \vec{U}$,

(iii) conservation of total energy in the inviscid limit.

These parallels justify further reduction under symmetry assumptions, explored in the next Section.
# 10. Derivation of the Self-Similar ODE System

We now specialize to spherically symmetric solutions in similarity coordinates, assuming that all fields depend only on the radial variable

$$r = |\vec{X}|, \quad \vec{X} = \frac{\vec{x}}{t^{\beta}}. \qquad (282)$$

This radial reduction transforms the self-similar PDE system into a set of coupled scalar ordinary differential equations.

The scalar functions $R(r)$, $U(r)$, $B(r)$ are invariant under the group $G_\lambda$ and belong to equivalence classes of similarity solutions modulo this dilation symmetry. The dynamical system lives effectively in a reduced phase space.

## 10.1. Radial Symmetry and Model Simplification

Assume purely radial fields

$$R(\vec{X}) = R(r), \quad \vec{U}(\vec{X}) = U(r)\hat{r}, \quad \vec{B}(\vec{X}) = B(r)\hat{r}, \qquad (283)$$

where $\hat{r} = \frac{\vec{X}}{|\vec{X}|}$ is the radial unit vector. Angular components vanish, and derivatives reduce to radial derivatives. For a scalar function $\varphi(r)$ and a radial vector field $V(r)\hat{r}$, we recall

$$\nabla \varphi = \frac{d\varphi}{dr} \hat{r}, \qquad (284)$$

$$\Delta \varphi = \frac{d^2 \varphi}{dr^2} + \frac{2}{r} \frac{d\varphi}{dr}, \qquad (285)$$

$$\nabla \cdot (\varphi \hat{r}) = \frac{1}{r^2} \frac{d}{dr} (r^2 \varphi), \qquad (286)$$

$$\Delta(\varphi \hat{r}) = \left( \frac{d^2 \varphi}{dr^2} + \frac{2}{r} \frac{d\varphi}{dr} - \frac{2}{r^2} \varphi \right) \hat{r}. \qquad (287)$$

## 10.2. Final Form of the Self-Similar System

Under the self-similar ansatz, the system of PDEs (271) (272) (276) reduces, under radial symmetry, to the following system of ordinary differential equations for the profiles $R(r)$, $U(r)$, and $B(r)$,

$$-\gamma R - \beta r \dot{R} + U \dot{R} + R \left( \dot{U} + \frac{2}{r} U \right) = 0, \qquad (C)$$

$$R \left( -\alpha U - \beta r \dot{U} + U \dot{U} \right) = \mu \left( \ddot{U} + \frac{2}{r} \dot{U} - \frac{2}{r^2} U \right) - \dot{P} + \frac{1}{\mu_0} B \dot{B}, \qquad (M)$$

$$-\alpha B - \beta r \dot{B} + U \dot{B} = \nu \left( \ddot{B} + \frac{2}{r} \dot{B} - \frac{2}{r^2} B \right) - \frac{2}{r} U B. \qquad (I)$$

Here, overdots denote derivatives, i.e. $\frac{d}{dr}$, with respect to the radial variable $r$, and $P(r)$ denotes the total pressure, determined via the barotropic closure relation.

For reference, the right-hand side expressions of the equations above can be explicitly grouped as follows,

$$\text{RHS}_{\text{self-similar}}^{(1)} = -\gamma R - \beta r \dot{R} + U \dot{R} + R \left( \dot{U} + \frac{2}{r} U \right), \qquad (288)$$

$$\text{RHS}_{\text{self-similar}}^{(2)} = R \left( -\alpha U - \beta r \dot{U} + U \dot{U} \right) - \mu \left( \ddot{U} + \frac{2}{r} \dot{U} - \frac{2}{r^2} U \right) + \dot{P} - \frac{1}{\mu_0} B \dot{B}, \qquad (289)$$

$$\text{RHS}_{\text{self-similar}}^{(3)} = -\alpha B - \beta r \dot{B} + U \dot{B} - \nu \left( \ddot{B} + \frac{2}{r} \dot{B} - \frac{2}{r^2} B \right) + \frac{2}{r} U B. \qquad (290)$$

### 10.2.1. Various Enlightenment

Let us examine the system of equations above more closely.
### 10.2.1.1. Non-local Integral Chain

The momentum, induction, and continuity equations are treated in this order.

(1) We perform a non-local integral reformulation of Eq. (M). We impose regularity at the origin together with the standard radial regularity conditions. In particular, we require

$$U(0) = 0, \qquad (291)$$

$$\lim_{r \to 0} r^2 \dot{U}(r) = 0. \qquad (292)$$

The following equalities are algebraically equivalent to the momentum ODE and expose an explicit Volterra-type nested integral; they will be used to generate a Neumann-style expansion consisting of an iterated chain of integrals,

$$\nu \left( \ddot{U}(r) + \frac{2}{r} \dot{U}(r) - \frac{2}{r^2} U(r) \right) = -\alpha U(r) - \beta r \dot{U}(r) + U(r) \dot{U}(r) - \frac{\dot{P}(r)}{R(r)} + \frac{1}{\mu_0 R(r)} B(r) \dot{B}(r), \qquad (293)$$

$$\frac{d}{dr} (r^2 \dot{U}(r)) = 2U(r) + \frac{r^2}{\nu R(r)} \left[ R(r) (-\alpha U(r) - \beta r \dot{U}(r) + U(r) \dot{U}(r)) - \dot{P}(r) + \frac{1}{\mu_0} B(r) \dot{B}(r) \right], \qquad (294)$$

$$\dot{U}(r) = \frac{1}{r^2} \int_0^r \left\{ 2U(t) + \frac{t^2}{\nu R(t)} \left[ R(t) (-\alpha U(t) - \beta t \dot{U}(t) + U(t) \dot{U}(t)) - \dot{P}(t) + \frac{1}{\mu_0} B(t) \dot{B}(t) \right] \right\} dt, \qquad (295)$$

$$U(r) = \underbrace{\int_0^r \frac{1}{\tau^2} \int_0^\tau \frac{s^2}{\nu R(s)} \left[ R(s) (-\alpha U(s) - \beta s \dot{U}(s) + U(s) \dot{U}(s)) - \dot{P}(s) + \frac{1}{\mu_0} B(s) \dot{B}(s) \right] ds d\tau}_{\text{forcing term}} + \underbrace{\int_0^r \frac{1}{\tau^2} \int_0^\tau 2U(s) ds d\tau}_{\text{self-interaction term}}, \qquad (296)$$

$$r^2 \dot{U}(r) = \int_0^r \frac{t^2}{\nu R(t)} \left[ R(t) (-\alpha U - \beta t \dot{U} + U \dot{U}) - \dot{P} + \frac{1}{\mu_0} B \dot{B} \right] dt + 2 \int_0^r U(t) dt. \qquad (297)$$

The identity (296) is a Volterra integral equation of the second kind. An iterated Neumann-like expansion that produces a chain of nested integrals can be constructed in the following way.

Describe the iteration sequence inside the same notation,

$$U^{(0)}(r) = 0, \qquad (298)$$

$$U^{(1)}(r) = \int_0^r \frac{1}{\tau^2} \int_0^\tau \frac{s^2}{\nu R(s)} \left[ R(s) \left( -\alpha U^{(0)}(s) - \beta s \dot{U}^{(0)}(s) \right) + U^{(0)}(s) \dot{U}^{(0)}(s) \right] - \dot{P}(s) + \frac{1}{\mu_0} B(s) \dot{B}(s) ds d\tau, \qquad (299)$$

$$U^{(n+1)}(r) = \int_0^r \frac{1}{\tau^2} \int_0^\tau 2U^{(n)}(s) ds d\tau + \int_0^r \frac{1}{\tau^2} \int_0^\tau \frac{s^2}{\nu R(s)} \left[ R(s) \left( -\alpha U^{(n)}(s) - \beta s \dot{U}^{(n)}(s) + U^{(n)}(s) \dot{U}^{(n)}(s) \right) \right] ds d\tau, \qquad (300)$$

for $n \ge 1$. All occurrences of $U^{(n)}$ inherit the same dot-differentiation convention, viz. $\dot{U}^{(n)} = \frac{d}{dr} U^{(n)}$.
Unwinding the recursion renders an explicit k-fold nested integral formula,

$$ U^{(k)}(r) = \sum_{m=1}^{k} \underbrace{\int_{0}^{r} \frac{d\tau_1}{\tau_1^2} \int_{0}^{\tau_1} d\tau_2 \int_{0}^{\tau_2} d\tau_3 \cdots \int_{0}^{\tau_{m-1}}}_{m \text{ nested integrals}} \left\{ 2U^{(k-m)}(\tau_m) \right. \\ \left. + \frac{\tau_m^2}{\nu R(\tau_m)} \left[ R(\tau_m) \left( -\alpha U^{(k-m)}(\tau_m) - \beta \tau_m \dot{U}^{(k-m)}(\tau_m) + U^{(k-m)}(\tau_m) \dot{U}^{(k-m)}(\tau_m) \right) \right] \right\} d\tau_m, \quad (301) $$

for any $k \ge 1$.

The full formal Picard sum is

$$ U(r) = \sum_{n=1}^{\infty} U^{(n)}(r), \quad (302) $$

with $U^{(n)}$ given by Eqq. (299) (300). Regularity at the origin fixes the integration constants and selects the regular branch of the solution, so that the single integral constants that would appear in primitive integrations vanish.

(2) The derivation for Eq. (1) proceeds analogously:

$$ \nu \left( \ddot{B}(r) + \frac{2}{r} \dot{B}(r) - \frac{2}{r^2} B(r) \right) = -\alpha B(r) - \beta r \dot{B}(r) + U(r) \dot{B}(r) + \frac{2}{r} U(r) B(r), \quad (303) $$

$$ \frac{d}{dr} \left( r^2 \dot{B}(r) \right) = 2B(r) + \frac{r^2}{\nu} \left[ -\alpha B(r) - \beta r \dot{B}(r) + U(r) \dot{B}(r) + \frac{2}{r} U(r) B(r) \right], \quad (304) $$

$$ \dot{B}(r) = \frac{1}{r^2} \int_0^r \left\{ 2B(t) + \frac{t^2}{\nu} \left[ -\alpha B(t) - \beta t \dot{B}(t) + U(t) \dot{B}(t) + \frac{2}{t} U(t) B(t) \right] \right\} dt, \quad (305) $$

$$ B(r) = \int_0^r \frac{1}{\tau^2} \int_0^\tau \left\{ 2B(s) + \frac{s^2}{\nu} \left[ -\alpha B(s) - \beta s \dot{B}(s) + U(s) \dot{B}(s) + \frac{2}{s} U(s) B(s) \right] \right\} ds d\tau. \quad (306) $$

Then

$$ B^{(0)}(r) = 0, \quad (307) $$

$$ B^{(1)}(r) = \int_{0}^{r} \frac{1}{\tau^2} \int_{0}^{\tau} \frac{s^2}{\nu} \left[ -\alpha B^{(0)}(s) - \beta s \dot{B}^{(0)}(s) + U^{(0)}(s) \dot{B}^{(0)}(s) + \frac{2}{s} U^{(0)}(s) B^{(0)}(s) \right] ds d\tau, \quad (308) $$

$$ B^{(n+1)}(r) = \int_0^r \frac{1}{\tau^2} \int_0^\tau 2B^{(n)}(s) ds d\tau \\ + \int_0^r \frac{1}{\tau^2} \int_0^\tau \frac{s^2}{\nu} \left[ -\alpha B^{(n)}(s) - \beta s \dot{B}^{(n)}(s) + U^{(n)}(s) \dot{B}^{(n)}(s) + \frac{2}{s} U^{(n)}(s) B^{(n)}(s) \right] ds d\tau, \quad (309) $$

for $n \ge 1$.

The nested integral notation denotes $m$ iterated integrations in the obvious order,

$$ \int_0^r \frac{d\tau_1}{\tau_1^2} \int_0^{\tau_1} d\tau_2 \int_0^{\tau_2} d\tau_3 \cdots \int_0^{\tau_{m-1}} d\tau_m. \quad (310) $$
And here one has

$$
\begin{equation}
\begin{split}
B^{(k)}(r) = \sum_{m=1}^{k} & \left( \int_{0}^{r} \frac{d\tau_1}{\tau_1^2} \int_{0}^{\tau_1} d\tau_2 \int_{0}^{\tau_2} d\tau_3 \cdots \int_{0}^{\tau_{m-1}} d\tau_m \right) \left\{ 2B^{(k-m)}(\tau_m) \right. \\
& \left. + \frac{\tau_m^2}{\nu} \left[ -\alpha B^{(k-m)}(\tau_m) - \beta \tau_m \dot{B}^{(k-m)}(\tau_m) + U^{(k-m)}(\tau_m) \dot{B}^{(k-m)}(\tau_m) \right. \right. \\
& \left. \left. + \frac{2}{\tau_m} U^{(k-m)}(\tau_m) B^{(k-m)}(\tau_m) \right] \right\} d\tau_m,
\end{split}
\tag{311}
\end{equation}
$$

and the sum is $B(r) = \sum_{n=1}^{\infty} B^{(n)}(r)$.

(3) The procedure for the ODE continuity is a little different. We rewrite Eq. (C) in an explicit first-order form for $\dot{R}$ and then expose a formal exponential. Under mild regularity at the origin we assume $R(0)$ finite.

Let

$$
-\gamma R(r) - \beta r \dot{R}(r) + U(r) \dot{R}(r) + R(r) \left( \dot{U}(r) + \frac{2}{r} U(r) \right) = 0. \quad (312)
$$

Rearranging the terms in (312) we obtain a linear homogeneous first-order ODE for $R$,

$$
\dot{R}(r) = a(r)R(r), \text{ where } a(r) = \frac{\dot{U}(r) + \frac{2}{r}U(r) - \gamma}{\beta r - U(r)}. \quad (313)
$$

Integrating (313) gives

$$
R(r) = R(0) \exp \left\{ \int_{0}^{r} a(t) dt \right\}. \qquad (314)
$$

Expand the factor $\frac{1}{\beta r - U(r)}$ as a geometric series in $\frac{U}{\beta r}$, and this is a formal expansion which is justified e.g. when $|U(r)| < \theta \beta r$ for some $\theta < 1$ on the integration interval. Thus

$$
\frac{1}{\beta r - U(r)} = \frac{1}{\beta r} \frac{1}{1 - \frac{U(r)}{\beta r}} = \sum_{k=0}^{\infty} \frac{U(r)^k}{\beta^{k+1} r^{k+1}}. \quad (315)
$$

Using (313) and (315) we write $a(r)$ as a formal power series,

$$
a(r) = \left( \dot{U}(r) + \frac{2}{r} U(r) - \gamma \right) \sum_{k=0}^{\infty} \frac{U(r)^k}{\beta^{k+1} r^{k+1}}. \quad (316)
$$

Insert (316) into (314) and expand the exponential in its Taylor series. Specify the mth exponential term for $m \ge 0$ by

$$
T_m(r) = \frac{1}{m!} \left( \int_0^r a(t) dt \right)^m = \frac{1}{m!} \underbrace{\int_0^r dt_1 \int_0^{t_1} dt_2 \cdots \int_0^{t_{m-1}} dt_m}_{m \text{ nested integrals}} \prod_{j=1}^m a(t_j), \quad (317)
$$

with the convention $T_0(r) = 1$.

Substituting (316) into (317) and expanding the product one arrives at

$$
T_m(r) = \frac{1}{m!} \sum_{k_1, \dots, k_m \ge 0} \left( \prod_{j=1}^{m} \beta^{-(k_j+1)} \right) \underbrace{\int_0^r dt_1 \int_0^{t_1} dt_2 \cdots \int_0^{t_{m-1}} dt_m}_{m \text{ nested integrals}} \prod_{j=1}^{m} \left\{ \left( \dot{U}(t_j) + \frac{2}{t_j} U(t_j) - \gamma \right) \frac{U(t_j)^{k_j}}{t_j^{k_j+1}} \right\}, \quad (318)
$$

for any $m \ge 1$.
It is sometimes convenient to rewrite the kernel product in symmetric form as

$$ K_m(t_1, \dots, t_m) = \prod_{j=1}^{m} \left[ \frac{U(t_j)^{k_j}}{t_j^{k_j+1}} \left( \dot{U}(t_j) + \frac{2U(t_j)}{t_j} - \gamma \right) \right], \quad (319) $$

so that

$$ T_m(r) = \frac{1}{m!} \sum_{k_1, \dots, k_m \ge 0} \beta^{-\sum_j (k_j+1)} \int_{0<t_m<\dots<t_1<r} K_m(t_1, \dots, t_m) dt_m \dots dt_1. \quad (320) $$

The full expansion for $R$ is the exponential series built from the $T_m$,

$$ R(r) = R(0) \sum_{m=0}^{\infty} T_m(r), \quad (321) $$

with $T_0(r) = 1$.

One may also articulate the cumulative logarithmic derivative as

$$ \int_0^r a(t) dt = \sum_{k=0}^{\infty} \beta^{-(k+1)} \int_0^r \left( \dot{U}(t) + \frac{2}{t} U(t) - \gamma \right) \frac{U(t)^k}{t^{k+1}} dt, \quad (322) $$

and the entire exponential operator may be written in compact form,

$$ R(r) = R(0) \exp \left\{ \sum_{k=0}^{\infty} \beta^{-(k+1)} \int_{0}^{r} \left( \dot{U}(t) + \frac{2}{t} U(t) - \gamma \right) \frac{U(t)^k}{t^{k+1}} dt \right\}. \quad (323) $$

This is formally equivalent to the series in (318) but often more convenient for analytical manipulations and asymptotic estimates.

### 10.2.1.2. Schur Complement Expansion

Put

$$ R(r, \tau) = R_*(r) + e^{\lambda \tau} \eta(r), \quad U(r, \tau) = U_*(r) + e^{\lambda \tau} v(r), \quad B(r, \tau) = B_*(r) + e^{\lambda \tau} w(r), \quad (324) $$

where $(R_*, U_*, B_*)$ is a stationary self-similar profile solving (C) (M) (I).

The eigenvalue problem for the linearised operator appears as

$$ \lambda \begin{pmatrix} \eta \\ v \\ w \end{pmatrix} = \lambda \begin{pmatrix} \eta \\ v \\ w \end{pmatrix}, \quad (325) $$

with the block operator

$$ \lambda = \begin{pmatrix} \lambda_{RR} & \lambda_{RU} & \lambda_{RB} \\ \lambda_{UR} & \lambda_{UU} & \lambda_{UB} \\ \lambda_{WR} & \lambda_{WU} & \lambda_{WB} \end{pmatrix}, \quad (326) $$

where each entry $\lambda_{ij}$ is a linear differential operator acting on the corresponding component of the perturbation vector.
The explicit expressions for the differential operator entries are

$$ \lambda_{RR}[\eta] = -\gamma\eta - \beta r \dot{\eta} - U_* \dot{\eta} - (\dot{U}_* + \frac{2}{r} U_*)\eta, \quad (327a) $$

$$ \lambda_{RU}[v] = v \dot{R}_* + R_* \dot{v} + \frac{2}{r} R_* v = (R_* \partial_r + \dot{R}_* + \frac{2}{r} R_*)[v], \quad (327b) $$

$$ \lambda_{RB} = 0, \quad (327c) $$

$$ \lambda_{UR}[\eta] = (-\alpha U_* - \beta r \dot{U}_* + U_* \dot{U}_*)\eta - \partial_r (P'(R_*)\eta), \quad (327d) $$

$$ \lambda_{UU}[v] = R_*(-\alpha v - \beta r \dot{v} - U_* \dot{v} - v \dot{U}_*) + \mu \Delta v, \quad (327e) $$

$$ \lambda_{UB}[w] = \frac{1}{\mu_0} (B_* \dot{w} + \dot{B}_* w), \quad (327f) $$

$$ \lambda_{WR}[\eta] = \eta \dot{B}_*, \quad (327g) $$

$$ \lambda_{WU}[v] = v \dot{B}_* + \frac{2}{r} B_* v, \quad (327h) $$

$$ \lambda_{WB}[w] = -\alpha w - \beta r \dot{w} - U_* \dot{w} + \nu \Delta w - \frac{2}{r} U_* w. \quad (327i) $$

To reduce the problem to the velocity component, we construct the Schur complement [76],

$$ Y(\lambda) = \lambda - \lambda_{UU} - \lambda_{UR}(\lambda - \lambda_{RR})^{-1}\lambda_{RU} - \lambda_{UB}(\lambda - \lambda_{WB})^{-1}\lambda_{WU}. \quad (328) $$

As a consequence, the eigenvalue condition reduces to the scalar non-local (generally pseudo-differential) equation,

$$ Y(\lambda)v = 0. \quad (329) $$

For $|\lambda| \gg 1$ one may expand the resolvents in a Neumann series,

$$ (\lambda - \lambda_{RR})^{-1} = \lambda^{-1} (\mathbb{I} + \lambda^{-1} \lambda_{RR} + \lambda^{-2} \lambda_{RR}^2 + \dots), \quad (330) $$

and similarly for $(\lambda - \lambda_{WB})^{-1}$. Inserting this into (329), one arrives at a formal polynomial eigenvalue problem, possibly quadratic or of higher order,

$$ (\lambda - \lambda_{UU} - \frac{1}{\lambda}(\lambda_{UR}\lambda_{RU} + \lambda_{UB}\lambda_{WU}) + \mathcal{O}(\lambda^{-2}))v = 0, \quad (331) $$

which can be reorganised as a quadratic (or higher) spectral problem amenable to Evans-function style techniques. Here, $\mathcal{O}(\lambda^{-2})$ is understood in the sense of Landau notation.

The natural boundary conditions are

$$ v(r) \text{ regular at } r = 0, \quad v \in L^2((0, \infty); \rho_{\text{weight}}(r)dr). \quad (332) $$

We now develop the semiclassical expansion of the Schur operator. Recall

$$ \begin{aligned} Y(\lambda) &= \lambda - \lambda_{UU} - \lambda_{UR}(\lambda - \lambda_{RR})^{-1}\lambda_{RU} - \lambda_{UB}(\lambda - \lambda_{WB})^{-1}\lambda_{WU}, \\ &= \lambda - \lambda_{UU} - \sum_{m \ge 0} \lambda^{-m-1} (\lambda_{UR}\lambda_{RR}^m\lambda_{RU} + \lambda_{UB}\lambda_{WB}^m\lambda_{WU}), \end{aligned} \quad (333) $$

Introducing the semiclassical parameter,

$$ \delta = \lambda^{-1}, \quad (334) $$

we rescale and arrive at

$$ \delta Y\left(\frac{1}{\delta}\right) = \mathbb{I} - \delta \lambda_{UU} - \delta^2 Y_1 - \delta^3 Y_2 + \mathcal{O}(\delta^4), \quad (335) $$

where the coefficients $Y_k$ are given explicitly by

$$ Y_1 = \lambda_{UR}\lambda_{RU} + \lambda_{UB}\lambda_{WU}, \quad (336) $$

$$ Y_2 = - (\lambda_{UR}\lambda_{RR}\lambda_{RU} + \lambda_{UB}\lambda_{WB}\lambda_{WU}), \quad (337) $$

$$ Y_k = \sum_{\substack{m \ge 0 \\ m+1=k}} [(-1)^m \lambda_{UR} \lambda_{RR}^m \lambda_{RU} + (-1)^m \lambda_{UB} \lambda_{WB}^m \lambda_{WU}], \quad k \ge 1. \quad (338) $$

Eq. (335) is the formal semiclassical asymptotic expansion of the Schur operator in powers of $\delta$.
Let $\lambda_{UU}$ admit an isolated simple eigenpair,

$$ \lambda_{UU} v_* = \lambda_* v_*, \quad \lambda_{UU}^* w_* = \bar{\lambda}_* w_*, \quad \langle w_*, v_* \rangle = 1. \qquad (339) $$

Let

$$ \Sigma(\lambda) = \lambda_{UR}(\lambda - \lambda_{RR})^{-1}\lambda_{RU} + \lambda_{UB}(\lambda - \lambda_{WB})^{-1}\lambda_{WU} \qquad (340) $$

be the coupling operator.

The Schur Eq. (329) is equivalent to

$$ (\lambda - \lambda_{UU})v = \Sigma(\lambda)v. \qquad (341) $$

Projecting onto $\text{span}\{v_*\}$ and its complement via Lyapunov-Schmidt reduction [53] [74], one gets the scalar secular equation,

$$ (\lambda - \lambda_*)a = \langle w_*, \Sigma(\lambda)(av_* + z(\lambda, a)) \rangle, \quad z(\lambda, a) = -(\lambda - \lambda_{UU})^{-1}\text{prj}^Q \Sigma(\lambda)(av_* + z), \qquad (342) $$

where

$$ v = av_* + z, \quad z \perp w_*, \qquad (343) $$

and

$$ \text{prj}^1 = v_* \otimes w_*, \quad \text{prj}^Q = \mathbb{I} - \text{prj}^1, \qquad (344) $$

with $\text{prj}^1$ the rank-one projector onto $\text{span}\{v_*\}$ and $\text{prj}^Q$ the complementary projector onto $\ker(w_*) = \{z : \langle w_*, z \rangle = 0\}$.

Neglecting $z$ in (342) gives the leading-order eigenvalue shift,

$$ \Delta\lambda^{(1)} = \langle w_*, \Sigma(\lambda)v_* \rangle = \langle w_*, \lambda_{UR}(\lambda_* - \lambda_{RR})^{-1}\lambda_{RU}v_* \rangle + \langle w_*, \lambda_{UB}(\lambda_* - \lambda_{WB})^{-1}\lambda_{WU}v_* \rangle. \qquad (345) $$

Higher-order corrections are achieved by iterating the equation for $z$ and expanding $\Sigma(\lambda)$ around $\lambda = \lambda_*$. In operator form, the full perturbative series may be written as

$$ \lambda = \lambda_* + \sum_{n \ge 1} \langle w_*, C_n^\lambda(\lambda_*)v_* \rangle, \qquad (346) $$

where

$$ C_1^\lambda(\lambda) = \Sigma(\lambda), \quad C_n^\lambda(\lambda) = \Sigma(\lambda) (R_{\text{prj}^Q}(\lambda)\text{prj}^Q \Sigma(\lambda))^{n-1}, \quad n \ge 2, \qquad (347) $$

with $R_{\text{prj}^Q}(\lambda) = \text{prj}^Q(\lambda - \lambda_{UU})^{-1}\text{prj}^Q$.

Let us try to give a precise mathematical identity to the semiclassical asymptotic expansion of the scalar secular equation.

We expand the coupling operator and solve the scalar secular equation by a systematic semiclassical expansion in the large-$\lambda$ and small-$\delta$ parameters,

$$ \Sigma(\lambda) = \sum_{k \ge 1} \lambda^{-k} Y_k + \mathcal{O}(\lambda^{-N-1}), \quad N \to \infty. \qquad (348) $$

Let

$$ \Delta\lambda \stackrel{\text{eqv}}{=} \lambda - \lambda_*, \quad R_{\text{prj}^Q}^* = R_{\text{prj}^Q}(\lambda_*) = \text{prj}^Q(\lambda_* - \lambda_{UU})^{-1}\text{prj}^Q, \qquad (349) $$

and write the unknown component as $v = av_* + z$ with $z \perp w_*$, solving the complementary equation,

$$ z = - \sum_{n \ge 1} (R_{\text{prj}^Q}(\lambda) \Sigma(\lambda))^n (av_*) . \qquad (350) $$

Inserting this into the projected (scalar) equation gives the exact series,

$$ \Delta\lambda = \sum_{n \ge 1} \langle w_*, \Sigma(\lambda) (R_{\text{prj}^Q}(\lambda) \Sigma(\lambda))^{n-1} v_* \rangle. \qquad (351) $$

Evaluating this expansion asymptotically at $\lambda = \lambda_*$, by expanding $\Sigma$ in powers of $\lambda_*^{-1}$, the asymptotic series will be

$$ \Delta\lambda \sim \sum_{p \ge 1} \lambda_*^{-p} \delta_p, \quad \text{equivalently} \quad \lambda = \lambda_* + \sum_{p \ge 1} \delta^p \delta_p, \qquad (352) $$
with the explicit low-order coefficients,

$$ \delta_1 = \langle w_*, Y_1 v_* \rangle, \qquad (353) $$

$$ \delta_2 = \langle w_*, Y_2 v_* \rangle + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle, \qquad (354) $$

$$ \delta_3 = \langle w_*, Y_3 v_* \rangle + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_2 v_* \rangle + \langle w_*, Y_2 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle \\ + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_1 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle, \qquad (355) $$

$$ \delta_4 = \langle w_*, Y_4 v_* \rangle + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_3 v_* \rangle + \langle w_*, Y_2 R_{\mathrm{prjQ}}^* Y_2 v_* \rangle + \langle w_*, Y_3 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle \\ + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_1 R_{\mathrm{prjQ}}^* Y_2 v_* \rangle + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_2 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle \\ + \langle w_*, Y_2 R_{\mathrm{prjQ}}^* Y_1 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle + \langle w_*, Y_1 R_{\mathrm{prjQ}}^* Y_1 R_{\mathrm{prjQ}}^* Y_1 R_{\mathrm{prjQ}}^* Y_1 v_* \rangle. \qquad (356) $$

More compactly, for each $p \ge 1$,

$$ \delta_p = \sum_{n=1}^{p} \sum_{\substack{k_1+\cdots+k_n=p \\ k_j \ge 1}} \left\langle w_*, Y_{k_1} (R_{\mathrm{prjQ}}^* Y_{k_2}) (R_{\mathrm{prjQ}}^* Y_{k_3}) \cdots (R_{\mathrm{prjQ}}^* Y_{k_n}) v_* \right\rangle. \qquad (357) $$

The above combinatorial formula generates an effective algorithm to compute $\delta_p$ to any desired order in $\lambda_*^{-1}$ once the operators $Y_k$ and the reduced resolvent $R_{\mathrm{prjQ}}^*$ are evaluated on the known eigenpair $(v_*, w_*)$.

Scholium 10.1 (Contour formula and remainder estimate). Assume $\Sigma(\lambda)$ is analytic for $|\lambda - \lambda_*| < \rho$ and that

$$ \| R_{\mathrm{prjQ}}^* \Sigma(\lambda_*) \| < \frac{1}{2}. \qquad (358) $$

The secular equation is equivalent to

$$ \Delta\lambda = F(\lambda), \quad F(\lambda) = \left\langle w_*, \Sigma(\lambda) (\mathbb{I} - R_{\mathrm{prjQ}}(\lambda)\Sigma(\lambda))^{-1} v_* \right\rangle. \qquad (359) $$

By Cauchy projection one has for each $p \ge 1$,

$$ \delta_p = \frac{1}{2\pi i} \oint_{|\zeta|=\rho} \zeta^{-p-1} \left\langle w_*, \Sigma(\lambda_* + \zeta) \left( \mathbb{I} - R_{\mathrm{prjQ}}(\lambda_* + \zeta) \Sigma(\lambda_* + \zeta) \right)^{-1} v_* \right\rangle d\zeta, \qquad (360) $$

provided the integrand is analytic inside the contour.

Truncating the series after $N$ terms,

$$ \Delta\lambda = \sum_{p=1}^{N} \lambda_{*}^{-p} \delta_{p} + R_{N}(\lambda_{*}), \qquad (361) $$

with remainder satisfying for some $C_N > 0$,

$$ |R_N(\lambda_*)| \le C_N \lambda_*^{-(N+1)}. \qquad (362) $$

The scalar equation may be written in Fredholm form,

$$ 0 = \det \left( \mathbb{I} - (\lambda - \lambda_{UU})^{-1} \Sigma(\lambda) \right), \qquad (363) $$

so that $\delta_p$ can also be expressed through cumulant-type expansions involving traces of products of $R_{\mathrm{prjQ}}^*$ and $Y_k$.

Scholium 10.2 (Spectral Schur reduction). We perform a spectral reduction by eliminating the R-block via Schur complement, write the resulting non-local reduced operator in kernel form, give a Green's function [29, § 3] for the RR-resolvent, and record a Fredholm determinant condition together with a first-order perturbation formula.

Let

$$ \langle \varphi, g \rangle_r = \int_0^\infty \varphi(r) g(r) r^2 dr \qquad (364) $$

be the scalar inner product. For vector-valued functions $V = (v, w)^T$ and $W = (\tilde{v}, \tilde{w})^T$ we set

$$ \langle V, W \rangle_r = \int_0^\infty \left( v(r) \tilde{v}(r) + w(r) \tilde{w}(r) \right) r^2 dr \qquad (365) $$
and

$$ (\lambda_{RR} - \lambda \mathbb{I}) [\mathcal{G}(\cdot, s; \lambda)](r) = \delta(r - s), \quad (366) $$

with boundary conditions: $\mathcal{G}(\cdot, s; \lambda)$ is the Green's function of $\lambda_{RR} - \lambda \mathbb{I}$; it is regular at $r = 0$ and decays as $r \to \infty$. The resolvent kernel on a test function $\vec{\phi}$ is

$$ (\lambda_{RR} - \lambda \mathbb{I})^{-1} [\vec{\phi}](r) = \int_{0}^{\infty} \mathcal{G}(r, s; \lambda) \vec{\phi}(s) ds. \quad (367) $$

Set the block vectors and matrices,

$$ V = \begin{pmatrix} v \\ w \end{pmatrix}, \quad \lambda_{VV} = \begin{pmatrix} \lambda_{UU} & \lambda_{UB} \\ \lambda_{WU} & \lambda_{WB} \end{pmatrix}, \quad \lambda_{VR} = \begin{pmatrix} \lambda_{UR} \\ \lambda_{WR} \end{pmatrix}, \quad \lambda_{RV} = \begin{pmatrix} \lambda_{RU} & \lambda_{RB} \end{pmatrix}. \quad (368) $$

For the Schur complement we may write

$$ Y(\lambda) = \lambda_{VV} - \lambda \mathbb{I} - \lambda_{VR} (\lambda_{RR} - \lambda \mathbb{I})^{-1} \lambda_{RV}, \quad (369) $$

so the reduced eigenvalue problem reads

$$ Y(\lambda)V = 0. \quad (370) $$

The non-local coupling in $Y(\lambda)$ can be written in integral kernel form. For $i, j \in \{U, B\}$, let $K_{ij}(r, s; \lambda)$ be the integral kernel of the composition $\lambda_{iR}(\lambda_{RR} - \lambda \mathbb{I})^{-1}\lambda_{Rj}$ through the identity

$$ \begin{aligned} \lambda_{iR}(\lambda_{RR} - \lambda \mathbb{I})^{-1} \lambda_{Rj}[\mathfrak{S}](r) &= \int_{0}^{\infty} K_{ij}(r, s; \lambda) \mathfrak{S}(s) ds \\ &= \int_{0}^{\infty} (\lambda_{iR}[\mathcal{G}(\cdot, s; \lambda)]) (r) (\lambda_{Rj}\mathfrak{S})(s) ds, \end{aligned} \quad (371) $$

for smooth compactly supported $\mathfrak{S}$, which motivates the distributional kernel identification,

$$ K_{ij}(r, s; \lambda) = (\lambda_{iR}[\mathcal{G}(\cdot, s; \lambda)]) (r) c_{Rj}(s), \quad (372) $$

where $c_{Rj}(s)$ denotes the pointwise action of $\lambda_{Rj}$ on the delta at $s$. The reduced operator acts componentwise as

$$ (Y(\lambda)V)_i(r) = (\lambda_{iV} - \lambda \mathbb{I})[V]_i(r) - \sum_j \int_0^\infty K_{ij}(r, s; \lambda) V_j(s) ds. \quad (373) $$

The spectral condition for $\lambda$ reduces to a Fredholm determinant for the compact perturbation of $(\lambda_{VV} - \lambda \mathbb{I})$,

$$ \lambda \text{ is an eigenvalue of } \lambda \iff \det_{L_r^2} (\mathbb{I} - (\lambda_{VV} - \lambda \mathbb{I})^{-1} \lambda_{VR} (\lambda_{RR} - \lambda \mathbb{I})^{-1} \lambda_{RV}) = 0. \quad (374) $$

To make contact with the scalar secular coefficients already displayed, one expands the Fredholm determinant by a trace-log relation.

The starting point is the Laurent expansions [45] of $\Sigma(\lambda)$ and the UU-resolvent,

$$ \Sigma(\lambda) = \sum_{k \ge 1} \lambda^{-k} Y_k, \quad (\lambda - \lambda_{UU})^{-1} = \sum_{n \ge 0} \lambda^{-(n+1)} \lambda_{UU}^n, \quad (375) $$

which upon multiplication give the composite series

$$ (\lambda - \lambda_{UU})^{-1} \Sigma(\lambda) = \sum_{p \ge 2} \lambda^{-p} \sum_{\substack{n \ge 0, k \ge 1 \\ n+k=p-1}} \lambda_{UU}^n Y_k. \quad (376) $$

Inserting this into the determinant, the trace-log expansion is

$$ \log \det (\mathbb{I} - (\lambda - \lambda_{UU})^{-1} \Sigma(\lambda)) = - \sum_{q \ge 1} \frac{1}{q} \operatorname{tr} \{ (\lambda - \lambda_{UU})^{-1} \Sigma(\lambda) \}^q, \quad (377) $$

with $\det_{L_r^2} (\mathbb{I} - (\lambda - \lambda_{UU})^{-1} \Sigma(\lambda)) = 0$.
Collecting terms that contribute to the same power $\lambda^{-p}$, we reorganize the sum over partitions,

$$ \log \det(\cdot) = - \sum_{p \ge 2} \lambda^{-p} \sum_{q \ge 1} \frac{(-1)^{q+1}}{q} \sum_{\substack{p_1+\dots+p_q=p \\ p_j \ge 2}} \text{tr} \left\{ \prod_{j=1}^{q} \left[ \sum_{\substack{n_j \ge 0, k_j \ge 1 \\ n_j+k_j=p_j-1}} \lambda_{UU}^{n_j} Y_{k_j} \right] \right\}. \quad (378) $$

One may equivalently rewrite this in cumulant form, making the analogy with connected diagrams explicit,

$$ \log \det(\cdot) = - \sum_{m \ge 1} \frac{1}{m} \sum_{\substack{k_1, \dots, k_m \ge 1 \\ n_1, \dots, n_m \ge 0}} \lambda^{-\left(\sum_{j=1}^m (n_j + k_j + 1)\right)} \text{tr} \left\{ \prod_{j=1}^m \lambda_{UU}^{n_j} Y_{k_j} \right\}. \quad (379) $$

By expanding the product inside the trace, we can identify the general coefficient of $\lambda^{-p}$ as

$$ C_p = \sum_{m \ge 1} \frac{(-1)^{m+1}}{m} \sum_{\substack{(n_1+k_1+1)+\dots+(n_m+k_m+1)=p \\ k_j \ge 1, n_j \ge 0}} \text{tr}(\lambda_{UU}^{n_1} Y_{k_1} \cdots \lambda_{UU}^{n_m} Y_{k_m}). \quad (380) $$

For clarity, the first few low-order coefficients are explicitly listed,

$$ C_2 = \text{tr}(\lambda_{UU}^0 Y_1) = \text{tr}(Y_1), \quad (381) $$

$$ C_3 = \text{tr}(\lambda_{UU}^1 Y_1) + \text{tr}(\lambda_{UU}^0 Y_2) = \text{tr}(\lambda_{UU} Y_1) + \text{tr}(Y_2), \quad (382) $$

$$ \begin{aligned} C_4 &= \text{tr}(\lambda_{UU}^2 Y_1) + \text{tr}(\lambda_{UU}^1 Y_2) + \text{tr}(\lambda_{UU}^0 Y_3) - \frac{1}{2} \text{tr}(Y_1^2) \\ &= \text{tr}(\lambda_{UU}^2 Y_1) + \text{tr}(\lambda_{UU} Y_2) + \text{tr}(Y_3) - \frac{1}{2} \text{tr}(Y_1^2). \end{aligned} \quad (383) $$

The same $C_p$ coefficients appear in the power expansion of the determinant,

$$ \log \det(\cdot) = - \sum_{p \ge 2} \lambda^{-p} C_p. \quad (384) $$

The rank-one reduction allows to project each operator word onto a scalar secular form,

$$ \text{tr} \left\{ \prod_{j=1}^{m} T_j \right\} \rightarrow \sum_{\text{cyclic}} \langle w_*, T_1 R_{\text{prjQ}}^* T_2 R_{\text{prjQ}}^* \cdots R_{\text{prjQ}}^* T_m v_* \rangle_r, \quad (385) $$

where each $T_j$ is an operator word of the form $\lambda_{UU}^{n_j} Y_{k_j}$ appearing in the traces above (the iterated $R_{\text{prjQ}}^*$-insertions implement the projection onto the range complement as in the main text).

Applying this rule to each trace term, the secular coefficients $\delta_p$ are

$$ \delta_1 = \langle w_*, Y_1 v_* \rangle_r, \quad (386) $$

$$ \delta_2 = \langle w_*, Y_2 v_* \rangle_r + \langle w_*, Y_1 R_{\text{prjQ}}^* Y_1 v_* \rangle_r, \quad (387) $$

$$ \delta_3 = \langle w_*, Y_3 v_* \rangle_r + \langle w_*, \lambda_{UU} Y_1 v_* \rangle_r + \langle w_*, Y_2 R_{\text{prjQ}}^* Y_1 v_* \rangle_r + \langle w_*, Y_1 R_{\text{prjQ}}^* Y_2 v_* \rangle_r, \quad (388) $$

$$ \begin{aligned} \delta_4 = & \langle w_*, Y_4 v_* \rangle_r + \langle w_*, \lambda_{UU} Y_3 v_* \rangle_r + \langle w_*, \lambda_{UU}^2 Y_1 v_* \rangle_r \\ & + \langle w_*, Y_3 R_{\text{prjQ}}^* Y_1 v_* \rangle_r + \langle w_*, Y_2 R_{\text{prjQ}}^* Y_2 v_* \rangle_r \\ & - \frac{1}{2} \langle w_*, Y_1 R_{\text{prjQ}}^* Y_1 v_* \rangle_r + \text{additional cyclic permutations}. \end{aligned} \quad (389) $$

The nested convolutional structure becomes visible: each $C_p$ is a weighted sum over all decompositions of $p$ into $(n_j + k_j + 1)$, with full operator words alternating between $\lambda_{UU}$-powers and $Y_k$-blocks.

At this point, the reduction to the rank-one secular problem is immediate,

$$ \text{tr}(\cdot) \rightarrow \sum_{\text{cyclic}} \langle w_*, (\cdots) v_* \rangle_r, \quad (390) $$

which reproduces the explicit combinatorial structure of the $\delta_p$, for instance,

$$ \delta_1 = \langle w_*, Y_1 v_* \rangle, \quad \delta_2 = \langle w_*, Y_2 v_* \rangle + \langle w_*, Y_1 R_{\text{prjQ}}^* Y_1 v_* \rangle, \quad (391) $$
and analogously for higher orders.

This determinant can therefore be regarded as an Evans-type function suitable for spectral computations.

If the composition $\lambda_{VR}\lambda_{RR}^{-1}\lambda_{RV}$ is bounded, one may expand the RR-resolvent in a Neumann series for $|\lambda| \ll 1$,

$$ (\lambda_{RR} - \lambda \mathbb{I})^{-1} = \lambda_{RR}^{-1} + \lambda \lambda_{RR}^{-2} + \lambda^2 \lambda_{RR}^{-3} + \dots, \quad (392) $$

valid for sufficiently small $|\lambda|$, provided $\|\lambda\lambda_{RR}^{-1}\| < 1$, from which the expansion

$$ Y(\lambda) = Y(0) - \lambda \mathbb{I} - \sum_{m \ge 1} \lambda^{m-1} \lambda_{VR} \lambda_{RR}^{-(m+1)} \lambda_{RV}. \quad (393) $$

Suppose $Y(0)$ has a simple root: let $V_0 \neq 0$ satisfy $Y(0)V_0 = 0$ and let $W_0$ be the adjoint mode for the adjoint of $Y(0)$ with respect to (364). Then the first non-trivial eigenvalue correction is

$$ \lambda_1 = - \frac{\langle [\partial_\lambda Y(\lambda)]_{\lambda=0} V_0, W_0 \rangle_r}{\langle V_0, W_0 \rangle_r} = \frac{\langle V_0 + \lambda_{VR} \lambda_{RR}^{-2} \lambda_{RV} V_0, W_0 \rangle_r}{\langle V_0, W_0 \rangle_r}. \quad (394) $$

Finally, let $\Phi_+(r; \lambda)$ be two independent solutions of

$$ (\lambda_{RR} - \lambda)\Phi = 0, \quad (395) $$

chosen so that $\Phi_-(\cdot; \lambda)$ is regular at $r = 0$ and $\Phi_+(\cdot; \lambda)$ satisfies the appropriate decay condition as $r \to \infty$.

Let

$$ \Gamma_w(s; \lambda) = \Phi_+(s; \lambda) \partial_s \Phi_-(s; \lambda) - \Phi_-(s; \lambda) \partial_s \Phi_+(s; \lambda) \quad (396) $$

be the Wronskian [36].

Then the Green's function admits

$$ \mathcal{G}(r, s; \lambda) = \frac{1}{\Gamma_w(s; \lambda)} \begin{cases} \Phi_-(r; \lambda)\Phi_+(s; \lambda), & r < s, \\ \Phi_+(r; \lambda)\Phi_-(s; \lambda), & r > s. \end{cases} \quad (397) $$

Eqq. (369) (373) (374) (394) are the operative formulæ to

(i) reduce the spectral problem to the V-subspace,

(ii) express the non-local coupling via kernels $K_{ij}$,

(iii) formulate an Evans-type determinant, and

(iv) obtain the leading perturbative eigenvalue correction.

◇

## 10.3. Algebraic Power-Law Profiles

We seek special power-law solutions of the reduced self-similar system by assuming profiles of the form

$$ R(r) = R_0 r^{-s}, \quad U(r) = U_0 r^q, \quad B(r) = B_0 r^m, \quad (398) $$

where $R_0, U_0, B_0 \in \mathbb{R}$ are constants and $s, q, m \in \mathbb{R}$ are scaling exponents. Substitution into the system entails algebraic constraints on these exponents.

(1) **Continuity equation (C),**

$$ -\gamma R - \beta r \dot{R} + U \dot{R} + R \left( \dot{U} + \frac{2}{r} U \right) = 0. \quad (399) $$

We compute

$$ \dot{R} = -s R_0 r^{-s-1}, \quad \dot{U} = q U_0 r^{q-1}. \quad (400) $$

Substituting into the equation

$$ -\gamma R_0 r^{-s} + \beta s R_0 r^{-s} - s R_0 U_0 r^{q-s-1} + q R_0 U_0 r^{q-s-1} + 2 R_0 U_0 r^{q-s-1} = 0, \quad (401) $$

$$ R_0 r^{-s} (-\gamma + \beta s) + R_0 U_0 r^{q-s-1} (-s + q + 2) = 0, \quad (402) $$

to balance powers of $r$, we require both terms to vanish independently, so

$$ \gamma = \beta s, \quad s = q + 2. \quad (403) $$
(2) Momentum equation (M),

$$R(-\alpha U - \beta r \dot{U} + U \ddot{U}) = \mu \left( \ddot{U} + \frac{2}{r} \dot{U} - \frac{2}{r^2} U \right) - \dot{P} + \frac{1}{\mu_0} B \dot{B}. \quad (404)$$

Neglecting pressure gradient for scaling purposes and focusing on leading-order balance, one has

$$\dot{U} = qU_0r^{q-1}, \quad \ddot{U} = q(q-1)U_0r^{q-2}. \quad (405)$$

The left-hand side scales as

$$RU \sim R_0 U_0 r^{q-s}, \quad \text{with factor } (-\alpha + \beta q + q^2). \quad (406)$$

Thus, balance gives

$$\alpha = \beta q + q^2. \quad (407)$$

(3) Induction equation (I),

$$-\alpha B - \beta r \dot{B} + U \ddot{B} = \nu \left( \ddot{B} + \frac{2}{r} \dot{B} - \frac{2}{r^2} B \right) - \frac{2}{r} U B. \quad (408)$$

Computing derivatives

$$\dot{B} = m B_0 r^{m-1}, \quad \ddot{B} = m(m-1) B_0 r^{m-2}. \quad (409)$$

Substituting these expressions into the governing balance and retaining only leading-order terms in $r$, we find

$$\text{LHS} \sim B_0 r^m (-\alpha - \beta m + qm). \quad (410)$$

Matching dominant scaling against the right-hand side then gets

$$\alpha = \beta m + qm + 2. \quad (411)$$

Summary of algebraic constraints. The power-law ansatz leads to the following non-linear system for the scaling exponents,

$$\left\{ \begin{array}{l} \gamma = \beta s, \end{array} \right. \quad (412a)$$

$$\left\{ \begin{array}{l} s = q + 2, \end{array} \right. \quad (412b)$$

$$\left\{ \begin{array}{l} \alpha = \beta q + q^2, \end{array} \right. \quad (412c)$$

$$\left\{ \begin{array}{l} \alpha = \beta m + qm + 2. \end{array} \right. \quad (412d)$$

Solving this system determines admissible exponent triples $(s, q, m)$ that are compatible with self-similar blow-up and asymptotic quiescence. These exact profiles offer valuable insight into the asymptotic structure of phase space and serve as guiding templates for numerical continuation and matched asymptotic expansions.

10.4. Barotropic Closure Conditions

Let the pressure obey a barotropic equation of state,

$$p(r) = AR(r)^{\gamma_p}, \quad A > 0, \quad \gamma_p > 1, \quad (413)$$

so that the total pressure satisfies

$$P(r) = p(r) + \frac{1}{2\mu_0} B(r)^2, \quad (414)$$

$$\dot{P}(r) = A\gamma_p R^{\gamma_p - 1} \dot{R} + \frac{1}{\mu_0} B \dot{B}. \quad (415)$$

This introduces a non-linearity that couples the momentum equation to the evolution of $R$ and $B$.

This algebraic closure can itself be unfolded into a nested integral chain once the continuity and induction expansions are inserted. We derive an explicit nested integral expansion for the total pressure $P(r)$ using the barotropic closure $p(r) = AR(r)^{\gamma_p}$ and the Picard-style induction iterates for $B$.

Let

$$P(r) = AR(r)^{\gamma_p} + \frac{1}{2\mu_0} B(r)^2. \quad (416)$$
We assume the continuity exponential expansion for R, see point (3) on p. 43,

$$R(r) = R(0) \sum_{m=0}^{\infty} T_m(r), \quad T_0(r) = 1, \qquad (417)$$

where the continuity kernel $a(t)$ and the mth exponential term $T_m$ are defined by

$$a(t) = \frac{\dot{U}(t) + \frac{2}{l}U(t) - \gamma}{\beta t - U(t)}, \qquad (418)$$

$$T_m(r) = \frac{1}{m!} \int_0^r dt_1 \int_0^{t_1} dt_2 \cdots \int_0^{t_{m-1}} dt_m \prod_{j=1}^m a(t_j), \quad m \ge 1. \qquad (419)$$

Apply the generalized binomial expansion to $R(r)^{\gamma_p}$ in the form $R(r)^{\gamma_p} = R(0)^{\gamma_p} (1 + X(r))^{\gamma_p}$ with $X(r) = \sum_{m \ge 1} T_m(r)$. The expansion is formal unless $|X(r)| < 1$ holds on the interval of interest.

We use the generalized binomial coefficients $\binom{\gamma_p}{n}$,

$$R(r)^{\gamma_p} = R(0)^{\gamma_p} \sum_{n=0}^{\infty} \binom{\gamma_p}{n} \left( \sum_{m=1}^{\infty} T_m(r) \right)^n, \qquad (420)$$

where

$$\binom{\gamma_p}{n} = \frac{\gamma_p(\gamma_p - 1)\cdots(\gamma_p - n + 1)}{n!}, \quad \binom{\gamma_p}{0} = 1, \qquad (421)$$

for $n \ge 1$.

Expanding the nth power produces a multi-sum over compositions of integers,

$$R(r)^{\gamma_p} = R(0)^{\gamma_p} \sum_{n=0}^{\infty} \binom{\gamma_p}{n} \sum_{m_1, \dots, m_n \ge 1} \prod_{i=1}^{n} T_{m_i}(r). \qquad (422)$$

Substitute the explicit nested-integral of each $T_{m_i}$ from (419) to produce an explicit multi-level formula,

$$R(r)^{\gamma_p} = R(0)^{\gamma_p} \sum_{n=0}^{\infty} \binom{\gamma_p}{n} \sum_{m_1, \dots, m_n \ge 1} \prod_{i=1}^{n} \left( \frac{1}{m_i!} \int_0^r dt_{i,1} \int_0^{t_{i,1}} dt_{i,2} \cdots \int_0^{t_{i,m_i-1}} dt_{i,m_i} \prod_{j=1}^{m_i} a(t_{i,j}) \right). \qquad (423)$$

Next expand the magnetic contribution using the Picard iterates for B. We posit the induction series,

$$B(r) = \sum_{q=1}^{\infty} B^{(q)}(r), \qquad (424)$$

where each $B^{(q)}(r)$ is the q-fold nested integral produced by the induction Picard scheme, see point (2) on p. 42.

Then

$$B(r)^2 = \sum_{p=1}^{\infty} \sum_{q=1}^{\infty} B^{(p)}(r) B^{(q)}(r). \qquad (425)$$

Combining the two expansions and using (416) one achieves the full formal nested-integral expansion for $P(r)$,

$$P(r) = AR(0)^{\gamma_p} \sum_{n=0}^{\infty} \binom{\gamma_p}{n} \sum_{m_1, \dots, m_n \ge 1} \prod_{i=1}^{n} \left( \frac{1}{m_i!} \int_0^r dt_{i,1} \int_0^{t_{i,1}} dt_{i,2} \cdots \int_0^{t_{i,m_i-1}} dt_{i,m_i} \prod_{j=1}^{m_i} a(t_{i,j}) \right) + \frac{1}{2\mu_0} \sum_{p=1}^{\infty} \sum_{q=1}^{\infty} B^{(p)}(r) B^{(q)}(r). \qquad (426)$$
## 10.5. Qualitative Behavior and Invariant Regions

The system (C) (M) (I) preserves key invariant structures

(1) if $R(0) > 0$, then positivity of density is preserved for small $r$,

(2) smooth initial data remain smooth under finite radial evolution,

(3) blow-up in derivatives is prevented near $r = 0$ by analyticity of solutions.

We study these properties more systematically using center manifold theory in the next Section.

# 11. Center Manifold Dynamics

To analyze regularity near the origin, we reformulate the system as a dynamical system in the radial variable $r$. This approach enables a precise classification of solution behavior in the vicinity of $r = 0$ using tools from invariant manifold theory.

## 11.1. Symmetry Groups and Invariant Structures

The similarity variables $(\vec{X}, \vec{U}, \vec{B}, R)$ form a dynamical system invariant under the dilation symmetry,

$$ G_{\lambda}: (\vec{x}, t) \mapsto (\lambda \vec{x}, \lambda^{\frac{1}{3}} t), \quad \lambda > 0, \qquad (427) $$

which induces the following scaling action on the profiles

$$ R(\vec{X}) \mapsto \lambda^{-\gamma} R(\lambda^{-1}\vec{X}), \qquad (428a) $$

$$ \vec{U}(\vec{X}) \mapsto \lambda^{-\alpha} \vec{U}(\lambda^{-1}\vec{X}), \qquad (428b) $$

$$ \vec{B}(\vec{X}) \mapsto \lambda^{-\alpha} \vec{B}, \qquad (428c) $$

$$ P(\vec{X}) \mapsto \lambda^{-\delta} P(\lambda^{-1}\vec{X}), \quad \delta = \gamma_p \beta. \qquad (428d) $$

On the dynamical variable $\vec{Y}$,

$$ G_{\lambda}: \vec{Y}(r) \mapsto \lambda^{-\rho} \vec{Y}(\lambda^{-1}r), \qquad (429) $$

where $\rho = \text{diag}(\gamma, \alpha+1, \alpha+2, \alpha+1, \alpha+2)$. For any solution $\vec{Y}(r)$, the rescaled profile $\vec{Y}_{\lambda}(r) = G_{\lambda}\vec{Y}(r)$ is also a solution. This implies that the vector field $\vec{F}(\vec{Y}, r)$ is in agreement with the quasi-homogeneity condition,

$$ \vec{F}(\lambda^{-\rho}\vec{Y}, \lambda^{-1}r) = \lambda^{-\rho-1}\vec{F}(\vec{Y}, r), \qquad (430) $$

so that $\vec{F}$ is equivariant under the scaling operator $G_{\lambda}$ up to a conformal weight. As a result, the dynamical system admits a conserved dilation current along group orbits,

$$ \frac{d}{dr} \left( \vec{Y} \cdot r \frac{\partial \vec{Y}}{\partial r} \right) = r \left\langle \frac{\partial \vec{Y}}{\partial r}, \frac{\partial \vec{Y}}{\partial r} \right\rangle + \vec{Y} \cdot \frac{\partial \vec{Y}}{\partial r} = \text{tr}(\rho) \|\vec{Y}\|^2 + \text{lower-order terms}, \qquad (431) $$

where the last equality holds when $\vec{F} = -\rho\vec{Y}$ with $\rho$ homogeneous of degree 0, so that $\nabla_{\vec{Y}} \cdot (\rho\vec{Y}) = \text{tr}(\rho)$ is constant along the flow.

The space of self-similar solutions is foliated into $G_{\lambda}$-orbits, each limiting an equivalence class of blow-up profiles. These correspond to invariant trajectories in the reduced ODE phase space, embedded in a $G_{\lambda}$-invariant manifold $M \subset F_{\varepsilon}$, where

$$ F_{\varepsilon} = \left\{ (\vec{X}, \vec{U}(\vec{X}), \vec{B}(\vec{X}), R(\vec{X}), P(\vec{X})) \middle| \vec{X} \in \mathbb{R}^3, \vec{U}, \vec{B} \in C^{\infty}(\mathbb{R}^3; \mathbb{R}^3), R, P \in \mathcal{C}^{\infty}(\mathbb{R}^3; \mathbb{R}_+) \right\} \qquad (432) $$

is the phase space. Invariance is regulated by the infinitesimal generator, a vector field for the symmetry,

$$ \vec{Q} = \underbrace{\vec{X} \cdot \nabla_{\vec{X}}}_{\text{spatial scaling}} + \gamma R \partial_{R} + \alpha (\vec{U} \cdot \nabla_{\vec{U}} + \vec{B} \cdot \nabla_{\vec{B}}) + \delta P \partial_{P}, \qquad (433) $$

satisfying $\mathcal{L}_{\vec{Q}}\omega_s = 0$ for the symplectic form

$$ \omega_s = \int_{\mathbb{R}^3} \left( \delta \vec{U} \wedge \delta \vec{B} + \frac{1}{R} \delta P \wedge \delta R \right) d^3 \vec{X} \qquad (434) $$
of MHD. Factoring out $G_λ$ simplifies the dynamics and enables compactification techniques. Similar approaches appear in [41], where modified Leray-type formulations are applied to the MHD-$α$ model, and in [82], which analyzes asymptotic self-similarity for the ideal compressible MHD system.

## 11.2. Dynamical Reformulation

We assign the notation

$$ \vec{Y}(r) = (R(r), U(r), \dot{U}(r), B(r), \dot{B}(r))^{\mathrm{T}} \in \mathbb{R}^{5}, \qquad (435) $$

the phase vector, so that the self-similar system reduces to a non-autonomous first-order ODE,

$$ \frac{d}{dr} \vec{Y} = \vec{F}(\vec{Y}(r), r), \qquad (436) $$

where $\vec{F}: \mathbb{R}^5 \times (0, \infty) \to \mathbb{R}^5$ is smooth for $r > 0$, but exhibits singularities at $r = 0$ due to radial terms such as $\frac{1}{r}$ and $\frac{1}{r^2}$. For any compact $\Omega \in \mathbb{R}^5$, suppose

(1) $\vec{F}(\cdot, r) \in \mathcal{C}^1(\Omega)$ for each $r > 0$,

(2) $\vec{F}$ is subject to local Lipschitz and polynomial growth,

$$ \|\vec{F}(\vec{Y}, r)\| \le C_M(1 + \|\vec{Y}\|), \quad \forall \vec{Y} \in \Omega, r \in [0, R]. \qquad (437) $$

Then by the Cauchy-Lipschitz or Picard-Lindelöf theorem [63] [50], there exists a unique local solution for initial data $\vec{Y}(0) = \vec{Y}_0 \in \Omega$.

For some $R > 0$, if the solution remains in $\Omega$ for all $r \in [0, R)$, ODE continuation guarantees global extendability beyond $r = R$, unless blow-up occurs,

$$ \sup_{r \in [0, R)} \|\vec{Y}(r)\| < \infty \implies \vec{Y}(r) \text{ extends beyond } R. \qquad (438) $$

Singularity formation corresponds to loss of compactness or divergence of derivatives in finite $r$. Let $F_\varepsilon^c$ be the compactified phase space,

$$ F_{\varepsilon}^{c} = \{ \vec{Y} \in \mathcal{C}^{1}((0, \infty); \mathbb{R}^{5}) \mid \|\vec{Y}\|_{F_{\varepsilon}^{c}} < \infty \}, \qquad (439) $$

with the weighted norm,

$$ \|\vec{Y}\|_{F_{\varepsilon}^{c}} = \sup_{r>0} \left( r^{\kappa} |\vec{Y}(r)| + r^{w} |\dot{\vec{Y}}(r)| \right), \qquad (440) $$

where the weights $(\kappa, w) \in \mathbb{R}^2$ are tuned to reflect asymptotic decay near $r = 0$ and $r \to \infty$; here $w = \kappa + 1$, in order to preserve scaling invariance under $r \mapsto \lambda r$. If $\vec{F}$ is Lipschitz on $F_\varepsilon^c$, one may construct global solutions via contraction mappings in Banach spaces.

**Scholium 11.1 (Weighted coercivity and variational structure).** Suppose the system admits an energy functional $\mathcal{E}[\vec{Y}]$ invariant under radial dilations. Then, selecting $(\kappa, w)$ to saturate the scaling symmetry, we obtain a weighted coercivity inequality,

$$ \mathcal{E}[\vec{Y}] \ge c \|\vec{Y}\|_{F_{\varepsilon}^{c}}^{2} - C_{M}, \qquad (441) $$

for some $c > 0$, uniform on finite-energy profiles. This implies stability under perturbations and facilitates the construction of gradient flows in $F_\varepsilon^c$. ◇

Also, if the vector field $\vec{F}$ derives from a variational principle,

$$ \vec{F} = -\nabla_{F_{\varepsilon}^{c}} \mathcal{E}[\vec{Y}], \qquad (442) $$

then $\vec{Y}(r)$ satisfies a weighted gradient flow, and existence follows from minimizing sequences or energy descent. In such a setting, blow-up corresponds to escape from level sets of $\mathcal{E}$ or loss of compactness in $F_\varepsilon^c$.
### 11.3. Continuation Criterion in $F_g^c$

The local dynamics of the radial self-similar system may be developed as an autonomous evolution in the configuration space $F_g^c$. The following lemma establishes local existence, uniqueness, and a continuation criterion on this conceptual path.

**Lemma 11.1 (Local existence, uniqueness, and continuation in $F_g^c$).** Let $\vec{F}: \mathbb{R}^5 \times (0, \infty) \to \mathbb{R}^5$ be continuous in $r$ and locally Lipschitz in $\vec{Y}$ uniformly on bounded subsets of $F_g^c$. For any initial datum $\vec{Y}(r_0) = \vec{Y}_0$, treat the Cauchy problem,

$$ \frac{d}{dr} \vec{Y}(r) = \vec{F}(\vec{Y}(r), r), \quad \vec{Y}(r_0) = \vec{Y}_0. \qquad (443) $$

Then

(i) (Local existence and uniqueness) There exist radius and Lipschitz constants $\rho_{\text{rad}}$, $c_L > 0$ and $\delta > 0$ such that the problem admits a unique solution,

$$ \vec{Y} \in \mathcal{C}^1([r_0, r_0 + \delta]; \mathbb{R}^5), \quad \sup_{r \in [r_0, r_0 + \delta]} \|\vec{Y}(r) - \vec{Y}_0\|_{F_g^c} \le \rho_{\text{rad}}, \qquad (444) $$

and

$$ \|\vec{F}(\vec{Y}_1, r) - \vec{F}(\vec{Y}_2, r)\|_{F_g^c} \le c_L \|\vec{Y}_1 - \vec{Y}_2\|_{F_g^c} \quad \forall \vec{Y}_1, \vec{Y}_2 \in F_g^c. \qquad (445) $$

(ii) (Continuous dependence) The solution depends continuously on the initial datum $\vec{Y}_0$ in the topology of $F_g^c$.

(iii) (Continuation criterion) Let $[r_0, r_{\text{max}})$ be the maximal forward interval of existence. If $r_{\text{max}} < \infty$ and

$$ \sup_{r \in [r_0, r_{\text{max}})} \|\vec{Y}(r)\|_{F_g^c} < \infty, \qquad (446) $$

then the solution can be extended beyond $r_{\text{max}}$. Consequently,

$$ \limsup_{r \nearrow r_{\text{max}}} \|\vec{Y}(r)\|_{F_g^c} = \infty \qquad (447) $$

is a necessary condition for finite-radius breakdown.

**Proof.** Fix $r_0$ and define the closed ball

$$ \bar{\mathbb{B}} = \{ \vec{Y} \in F_g^c : \|\vec{Y} - \vec{Y}_0\|_{F_g^c} \le \rho_{\text{rad}} \}. \qquad (448) $$

For $\delta > 0$ to be chosen, adopt the Picard map on continuous curves by

$$ (T\vec{Y})(r) = \vec{Y}_0 + \int_{r_0}^{r} \vec{F}(\vec{Y}(s), s) \, ds, \quad r \in [r_0, r_0 + \delta], \qquad (449) $$

where $\int_{r_0}^{r}$ is taken in the Bochner sense in the Banach space $F_g^c$, thereby making $T: \mathcal{C}([r_0, r_0 + \delta]; F_g^c) \to \mathcal{C}([r_0, r_0 + \delta]; F_g^c)$.

Iterating the operator $T$ defines a hierarchy of nested integral chains,

$$ (T^{(2)}\vec{Y})(r) = \vec{Y}_0 + \int_{r_0}^{r} \vec{F}\left(\vec{Y}_0 + \int_{r_0}^{s_1} \vec{F}(\vec{Y}(s_2), s_2) ds_2, s_1\right) ds_1, \qquad (450) $$

$$ (T^{(3)}\vec{Y})(r) = \vec{Y}_0 + \int_{r_0}^{r} \vec{F}\left(\vec{Y}_0 + \int_{r_0}^{s_1} \vec{F}\left(\vec{Y}_0 + \int_{r_0}^{s_2} \vec{F}(\vec{Y}(s_3), s_3) ds_3, s_2\right) ds_2, s_1\right) ds_1, \qquad (451) $$

and so on, producing an $m$-fold Volterra-type structure,

$$ (T^{(m)}\vec{Y})(r) = \vec{Y}_0 + \int_{r_0}^{r} \vec{F}(\vec{Y}(t_m), t_m) \circ \dots \circ \vec{F}(\vec{Y}(t_1), t_1) \, dt_m \dots dt_1, \qquad (452) $$

whose integrand encodes cumulative non-linear feedback of the flow variables. The convergence of $T^{(m)}\vec{Y}$ as $m \to \infty$ corresponds to the local fixed point of the evolution in $F_g^c$.

Let

$$ C_T = \sup_{\substack{\|\vec{Y} - \vec{Y}_0\|_{F_g^c} \le \rho_{\text{rad}} \\ s \in [r_0, r_0 + \delta]}} \|\vec{F}(\vec{Y}, s)\|_{F_g^c} < \infty \qquad (453) $$
be a uniform bound for the vector field $\vec{F}$ on the closed ball $\{\|\vec{Y} - \vec{Y}_0\|_{F_g^c} \le \rho_{\text{rad}}\}$, so that the Picard map remains invariant on its domain.

Picking $\delta$ such that $C_T\delta \le \rho_{\text{rad}}$ and $c_L\delta < 1$, the operator $T$ maps the closed subset

$$ \hat{\Omega}_{\mathbb{B}} = \left\{ \vec{Y} \in \mathcal{C}([r_0, r_0 + \delta]; \mathbb{R}^5) : \sup_{r \in [r_0, r_0 + \delta]} \|\vec{Y}(r) - \vec{Y}_0\|_{F_g^c} \le \rho_{\text{rad}} \right\} \quad (454) $$

into itself and fulfills

$$ \|T\vec{Y} - T\tilde{\vec{Y}}\|_{F_g^c} \le c_L \delta \sup_{r \in [r_0, r_0 + \delta]} \|\vec{Y}(r) - \tilde{\vec{Y}}(r)\|_{F_g^c}. \quad (455) $$

Accordingly $T$ is a strict contraction on $\hat{\Omega}_{\mathbb{B}}$ and admits a unique fixed point by Banach's theorem. This shows (i) and (ii).

To prove (iii), assume $[r_0, r_{\text{max}})$ is maximal with $r_{\text{max}} < \infty$ and $\sup_{r \in [r_0, r_{\text{max}})} \|\vec{Y}(r)\|_{F_g^c} < \infty$. One may restart the local argument at any $r_1 \in (r_{\text{max}} - \varepsilon, r_{\text{max}})$, contradicting maximality. And the solution can fail to continue only if the norm diverges as $r \nearrow r_{\text{max}}$. $\square$

### 11.4. Taylor Expansion at $r = 0$

We seek smooth solutions near the origin by imposing analyticity at $r = 0$ and expanding the profiles as formal even and odd Taylor series,

$$ R(r) = R_0 + R_2 r^2 + R_4 r^4 + \dots, \quad (456) $$

$$ U(r) = U_1 r + U_3 r^3 + U_5 r^5 + \dots, \quad (457) $$

$$ B(r) = B_1 r + B_3 r^3 + B_5 r^5 + \dots, \quad (458) $$

with the regularity conditions

$$ \dot{R}(0) = 0, \quad U(0) = 0, \quad B(0) = 0. \quad (459) $$

We substitute these expansions into the reduced system (C) (M) (I), and isolate each algebraic constraint $A_k$ arising at order $r^k$, for each $k \in \mathbb{N}$.

(1) Order $\mathcal{O}(r^0)$ in (C):

$$ \text{LHS}_{(C)} = -\gamma R_0 + 3R_0 U_1 + \mathcal{O}(r^2). \quad (460) $$

The zeroth-order constraint reads

$$ A_0^{(C)} : U_1 = \frac{\gamma}{3}. \quad (461) $$

(2) Order $\mathcal{O}(r^1)$ in (M):

$$ \text{LHS}_{(M)} = R_0 (-\alpha U_1 - \beta U_1 + U_1^2) r + \mathcal{O}(r^3), \quad (462) $$

$$ \mu \left( \ddot{U} + \frac{2}{r} \dot{U} - \frac{2}{r^2} U \right) = 6\mu U_3 r + \mathcal{O}(r^3), \quad (463) $$

$$ \dot{P} = \frac{d}{dr} \left( K_c R^{\gamma_p} + \frac{1}{2\mu_0} B^2 \right) = \left( 2\gamma_p K_c R_0^{\gamma_p-1} R_2 + \frac{B_1^2}{\mu_0} \right) r + \mathcal{O}(r^3), \quad (464) $$

$$ B\dot{B} = B_1^2 r + \mathcal{O}(r^3). \quad (465) $$

Putting it all together,

$$ R_0 (-\alpha - \beta + U_1) U_1 r = 6\mu U_3 r - 2\gamma_p K_c R_0^{\gamma_p-1} R_2 r + \mathcal{O}(r^3). \quad (466) $$

Oversimplify,

$$ A_1^{(M)} : R_0(-\alpha - \beta + U_1)U_1 = 6\mu U_3 - 2\gamma_p K_c R_0^{\gamma_p-1} R_2. \quad (467) $$

And substituting $U_1 = \frac{\gamma}{3}$ from $A_0^{(C)}$,

$$ R_2 = \frac{1}{2\gamma_p K_c R_0^{\gamma_p-1}} \left[ 6\mu U_3 - R_0 \left(-\alpha - \beta + \frac{\gamma}{3}\right) \cdot \frac{\gamma}{3} \right]. \quad (468) $$
(3) Order $\mathcal{O}(r^1)$ in (I):

$$ \text{LHS}_{(I)} = B_1(-\alpha - \beta + U_1)r + \mathcal{O}(r^3) \quad (469) $$

$$ \nu \left( \ddot{B} + \frac{2}{r}\dot{B} - \frac{2}{r^2}B \right) = 6\nu B_3 r + \mathcal{O}(r^3) \quad (470) $$

$$ -\frac{2}{r}UB = -2U_1B_1r + \mathcal{O}(r^3). \quad (471) $$

From this,

$$ B_1(-\alpha - \beta + U_1) = 6\nu B_3 - 2U_1B_1. \quad (472) $$

Solving for $B_3$,

$$ \mathcal{A}_1^{(1)} : B_3 = \frac{B_1}{6\nu}(-\alpha - \beta + 3U_1) = \frac{B_1}{6\nu}(-\alpha - \beta + \gamma). \quad (473) $$

Summary of leading-order Taylor coefficients,

$$ U_1 = \frac{\gamma}{3}, \quad (474) $$

$$ R_2 = \frac{1}{2\gamma_p K_c R_0^{\gamma_p - 1}} \left[ 6\mu U_3 - R_0 \left( -\alpha - \beta + \frac{\gamma}{3} \right) \cdot \frac{\gamma}{3} \right], \quad (475) $$

$$ B_3 = \frac{B_1}{6\nu}(-\alpha - \beta + \gamma). \quad (476) $$

Each coefficient is determined recursively from lower-order terms. The hierarchy $\mathcal{A}_k$ consists of algebraic constraints comprising the necessary compatibility between the self-similar structure and analyticity at the origin.

## 11.5. Jacobian and Local Spectral Features

Let

$$ J_0 = D_{\vec{Y}} \vec{F} \Big|_{\vec{Y} = \vec{Y}_0(0), r=0} \in \mathbb{R}^{5 \times 5}, \quad (477) $$

be the Jacobian matrix of $\vec{F}$ at the origin, where $J_0$ is a 5 × 5 matrix whose entries are the partial derivatives $\frac{\partial F_i}{\partial Y_j}$ evaluated at the initial configuration $\vec{Y}_0(0)$.

The spectrum of $J_0$,

$$ \text{spec}(J_0) \stackrel{\text{eqv}}{=} \text{spec}(D\vec{F}|_{r=0}) = \{\lambda_i \in \mathbb{C} \mid i = 1, \dots, 5\}, \quad (478) $$

conditions the linear flow near $r=0$. According to the Hartman-Grobman theorem [34] [30],

(1) if $\Re(\lambda_i) \neq 0$ for all $i$, the system is topologically conjugate to its linearization.

(2) if some $\Re(\lambda_i) = 0$, then the dynamics are governed by the center manifold, whose dimension equals the number of such eigenvalues.

## 11.6. Existence of a Local Center Manifold

If $J_0$ has three eigenvalues with $\Re(\lambda_i) = 0$, then there exists a $\mathcal{C}^k$-smooth local center manifold

$$ \mathcal{M}^c(0) \subset \mathbb{R}^5, \dim \mathcal{M}^c = 3, \quad (479) $$

parameterized by $(R_0, U_1, B_1)$. It contains all finite-energy regular solutions near the origin,

$$ R(r) = R_0 + \mathcal{O}(r^2), \quad U(r) = U_1 r + \mathcal{O}(r^3), \quad B(r) = B_1 r + \mathcal{O}(r^3). \quad (480) $$

## 11.7. Autonomous Dynamics and Invariant Geometry

Set $\xi = \log r$, so that

$$ \frac{d\vec{Y}}{d\xi} = \tilde{F}(\vec{Y}) = r\tilde{F}(\vec{Y}, r). \quad (481) $$

The fixed point $\vec{Y}_0 = (R_0, 0, 0, 0, 0)$ admits a center subspace of dimension 3,

$$ D\tilde{F}(\vec{Y}_0) = \tilde{J}_0, \quad \text{spec}(\tilde{J}_0) = \{0, 0, 0, \lambda_-, \lambda_+\}, \quad \Re(\lambda_-) < 0. \quad (482) $$