Quantum simulation of electronic structure via quantum fast multipole method

ArXiv: 2510.07380

🎯 Pitch

Long-range Coulomb interactions in quantum chemistry simulations incur an O(η²) cost per step—until now. By porting the classical fast multipole method onto a quantum computer with a clever sorting-network architecture, this work slashes the particle-scaling overhead by roughly a factor of η, enabling near-linear complexity for systems with more than a thousand electrons.

---

1. Executive Summary

This paper proposes a quantum algorithm for simulating electronic structure that overcomes a long-standing O(η²) bottleneck in computing the Coulomb potential among η electrons, achieving complexity nearly linear in the number of particles via a quantum implementation of the classical fast multipole method. Using a real-space first-quantised representation with high-order product formulae, the authors introduce a quantum fast multipole method that avoids the data-access overhead that would negate classical speedups—replacing position-dependent lookups with a sorting-network approach that combines quantum sort with shifted Morton orderings to retrieve interaction-list information at fixed locations, and moving multipole information by copying it along sorted particle registers rather than querying per-box registers. The resulting algorithm achieves a gate complexity of t(η^4/3N^1/3 + η^1/3N^2/3)(ηNt/ϵ)^o(1), representing roughly an O(η) speedup over the most comparable prior first-quantised product-formula methods and providing the lowest known complexity for any approach when N < η⁷ (the regime of practical interest), though the authors caution that large constant factors and overheads from coherent sorts and redundant Morton orderings would require η ≳ 10³ electrons to realize this asymptotic advantage in practice.

2. Context and Motivation

The Core Problem: The O(η²) Coulomb Bottleneck in Quantum Simulation

The fundamental computational bottleneck this paper addresses is deceptively simple: when simulating a system of η interacting electrons on a quantum computer, computing the Coulomb potential—the electrostatic repulsion between every pair of particles—naïvely requires O(η²) operations because you must sum over all η(η−1)/2 pairwise interactions. For large systems, this quadratic scaling becomes the dominant cost, overshadowing all other parts of the simulation.

The problem arises specifically in first-quantised real-space simulations where electrons are represented by their positions in a spatial grid. In this representation, the potential energy operator V (Eq. 4 in the paper) takes the form of a double sum over all particle pairs:

$V = \sum_{i \neq j} \sum_{p,q \in G} \frac{1}{2\|r_p - r_q\|} |p\rangle\langle p|_i \otimes |q\rangle\langle q|_j$

For η electrons, evaluating this directly means computing η² Coulomb interaction terms. While this is manageable for small molecules (η ~ 10–100), it becomes prohibitive for the systems that quantum computers are theoretically well-suited to tackle: warm dense matter, degenerate plasmas, and condensed-phase materials where η can reach thousands. Rubin et al. [5], for example, simulated a deuterium plasma with 1729 electrons—a regime where the O(η²) cost already dominates.

This bottleneck is particularly frustrating because classical computing solved the analogous problem decades ago. The fast multipole method (FMM), introduced by Rokhlin [20] and Greengard and Rokhlin [21] in the 1980s, computes the Coulomb potential among η particles with O(η polylog(1/ϵ)) complexity by hierarchically grouping distant charges into multipole expansions rather than summing every pair. The Barnes-Hut tree code [19] provides a similar O(η log η) speedup. These methods are standard in classical molecular dynamics, astrophysical N-body simulations, and electronic structure calculations. Yet, prior to this paper, no quantum algorithm had successfully replicated this scaling—the very speedup that makes classical simulations of large systems tractable was unavailable on quantum computers.

Why This Matters: Beyond Just Faster Chemistry

The paper's motivation extends beyond merely making chemistry simulations faster. There are several layers of significance:

Quantum advantage for electronic structure depends on asymptotic scaling. The promise of quantum computing for chemistry and materials science—computing molecular energies via quantum phase estimation, simulating reaction dynamics [1], modeling condensed phases and plasmas [2–5]—hinges on algorithms whose gate complexity grows slowly enough with system size that they remain feasible on error-corrected quantum hardware. Every factor of η in the complexity multiplies the number of gates. An O(η²) Coulomb step means that doubling the system size quadruples the cost of that step alone. Removing this quadratic factor makes the difference between simulating tens of electrons and simulating thousands.

First-quantised representations are otherwise superior for N ≫ η. In electronic structure, N (the number of grid points or basis functions) is typically much larger than η (the number of electrons). Second-quantised approaches require at least O(N) complexity because they encode occupation numbers for every orbital. First-quantised methods store only the η electron positions, achieving O(η) space and avoiding the O(N) factor—but this advantage is undermined if the Coulomb potential contributes an O(η²) factor. Removing the quadratic Coulomb scaling makes first quantisation genuinely superior to second quantisation for the practically relevant regime where N ≫ η.

The bottleneck affects broader classes of simulations. While the paper focuses on electronic structure, the O(η²) Coulomb problem appears in any quantum simulation of charged particles—plasmas for fusion energy research [5], non-relativistic quantum electrodynamics [16], and simulations of materials under extreme conditions. An efficient quantum Coulomb solver therefore has implications beyond quantum chemistry.

Product formulae are competitive but were held back. Recent work [13, 14] has shown that product formulae (Trotter-Suzuki decompositions) can be competitive with more sophisticated simulation methods like quantum signal processing and qubitisation because they achieve low error in practice, despite their polynomial rather than logarithmic scaling in ϵ. However, the O(η²) Coulomb cost has been a "long-standing bottleneck in simulation of quantum chemistry based on product formulae" (Section VII). Resolving this bottleneck makes product-formula-based simulations significantly more attractive for large systems.

Prior Approaches and Their Limitations

The paper identifies several existing approaches to quantum simulation of electronic structure, each with specific shortcomings that the quantum FMM addresses.

Naïve first-quantised product formulae: O(η²) Coulomb. The original Kassal et al. [17] approach to quantum simulating chemistry in first quantisation computed the Coulomb potential by directly summing over all η² pairwise interactions. While conceptually straightforward—iterate over all particle pairs and apply phase factors based on their distances—this method inherits the full quadratic cost. Rubin et al. [5] improved on this with tighter Trotter bounds but retained the O(η²) Coulomb step, achieving a total complexity of η² times the number of Trotter steps. The paper explicitly positions itself as reducing this η² factor to η (up to logarithmic factors), improving over Ref. [5] "by a factor of η."

Second-quantised plane wave methods: O(N) replaces O(η²) but N ≫ η. Babbush et al. [15] and Low et al. [18] observed that in second quantisation, the Coulomb potential can be computed with O(N) complexity using fast Fourier transforms, avoiding the double sum over particles. However, second quantisation inherently requires representing all N basis functions, so the overall complexity is at least O(N). Since N is typically much larger than η (for accurate simulations, N can be 10⁶–10⁹ while η might be 10²–10³), the O(N) factor often dominates. The paper notes that "second quantisation yields a complexity at least linear in N for the sum, which is much larger for the typical scenario where N ≫ η." The quantum FMM achieves O(η) for the Coulomb part while maintaining the O(η) advantage of first quantisation, making it superior when N is large relative to η.

Block-encoding and qubitisation approaches: favourable ϵ-scaling but higher η-dependence. More recent simulation algorithms based on linear combinations of unitaries [7, 8], quantum walks [9], qubitisation [12], and interaction-picture methods [35] achieve O(log(1/ϵ)) scaling with error, which is exponentially better than product formulae's poly(1/ϵ). However, these methods typically have worse scaling with η. The first-quantised qubitisation approach of Babbush et al. [35] has complexity $\tilde{O}(\eta^3 (N/\Omega)^{1/3}/\epsilon)$ in the interaction picture. The quantum FMM's product-formula approach achieves $\tilde{O}(\eta^{5/3} N^{1/3}/\Omega^{1/3})$ (ignoring the second term), which is better than the interaction-picture method when $N < \eta^7$ in the thermodynamic limit—the regime the authors identify as practical.

Direct translation of classical FMM fails due to quantum data access overhead. The most natural idea—simply implementing the classical FMM as quantum gates—does not work, and understanding why is essential to appreciating the paper's contribution. In the classical FMM (Algorithms 1 and 2), the algorithm determines which box an electron belongs to by examining its classical position coordinates, then accesses that box's multipole data. In a quantum algorithm, the electron positions are stored in quantum registers in superposition—you cannot simply read a position value and use it as a memory address without collapsing the superposition or incurring prohibitive overhead. As the paper explains in Section I:

"The difficulty with applying these methods in quantum algorithms is that a direct translation of the various classical algorithms would result in data accesses in locations governed by the values in quantum registers. Each data access has complexity corresponding to the number of data locations, which here is η, increasing the overall complexity by this factor."

In other words, if you have a quantum register storing a box index, and you need to retrieve the multipole data for that box from an array of η boxes, a naïve implementation would need to query all η boxes to find the matching one—restoring the O(η²) cost you were trying to avoid. The paper cites Babbush et al. [4] for a detailed discussion of these issues and notes that the FMM approach of Childs et al. [22], which uses QRAM for these data accesses, "incurs this extra overhead of η" implicitly. The result is that "a direct translation... would result in an overall complexity that is larger than O(η²), negating the speedup provided by either tree codes or FMMs."

Self-correction and verifier-based methods (for context): not relevant here. Unlike the previous example paper on LLM test-time compute, there are no conflicting empirical results to reconcile. Instead, there is a clear, acknowledged technical barrier: quantum data access overhead has prevented quantum FMM implementations. The paper's contribution is to design around this barrier rather than to resolve contradictory findings.

How This Paper Positions Itself

The paper frames its contribution not as inventing a new classical algorithm but as solving the specific problem of quantum data access that prevents direct translation of the classical FMM to quantum computers. The key insight is stated in Section I:

"Our solution to recover the O(η log η) scaling of the classical algorithm is related to the implementation of a quantum sort in which a sorting network is used to make comparisons at fixed locations so that the overhead from accessing quantum data is avoided."

The approach draws on a known technique in quantum algorithms—quantum sorting networks [23, 24]—but applies it in a novel way to the FMM context. Instead of looking up box data by quantum register values, the algorithm moves the electron data into fixed registers corresponding to their correct boxes via a sequence of sorts and controlled swaps. This transforms the problem from "access memory at a quantum-specified address" (hard) to "arrange the data so that the correct information is at a classically known location" (doable via sorting).

The paper explicitly distinguishes its approach from prior attempts at quantum FMM. The Childs et al. [22] FMM "incurs this extra overhead of η" because it relies on QRAM. The current work avoids QRAM entirely by using sorting networks with fixed comparison locations. The paper also notes a very recent alternative: Stetina and Wiebe [16] proposed using Gauss's law as a constraint to avoid the O(η²) Coulomb overhead, but that approach is specific to non-relativistic QED and uses a different mechanism; the sorting-based FMM is more general.

The paper also carefully positions itself within the broader landscape of quantum simulation methods. Table I provides a detailed comparison of gate complexities across 10 prior works spanning 2017–2025. The quantum FMM's complexity of $\eta(\eta^{2/3}(N/\Omega)^{1/3} + (N/\Omega)^{2/3})(\eta N)^{o(1)}/\epsilon^{1+o(1)}$ shows a clear η factor improvement over the next-best first-quantised Trotter method (Rubin et al. [5], which has an η² prefactor). The paper acknowledges that for very large N (N > η⁷), the interaction-picture method of Babbush et al. [35] remains superior, but identifies N < η⁷ as the "regime of practical interest."

Crucially, the paper is transparent about the practical limitations. It repeatedly cautions that the classical FMM itself only becomes advantageous at large particle numbers—the Greengard and Rokhlin data shows break-even points of 500–5000 particles depending on accuracy. The quantum implementation adds further overhead from coherent sorts, redundant Morton orderings, and the need to associate multipole information with every electron rather than with boxes. The paper estimates in Appendix B that the break-even point for the quantum FMM could be "an order of magnitude" higher than the classical one, potentially requiring hundreds of thousands of electrons, and that physical implementation would need "access to a quantum computer with a number of physical qubits exceeding ~10⁸." The focus is on "showing that it is possible to achieve the O(η) complexity at all," with optimisation left to future work.

Finally, the paper positions its sorting-based data movement technique as having broader applicability. The conclusion notes that "the methods developed here to avoid the data-access overhead show promise for application to quantum versions of many other fast summation methods," citing hierarchical matrix methods, fast Ewald summation, and other kernel-independent FMM variants [25–31, 62].

3. Technical Approach

3.1 Reader orientation

The paper constructs a quantum algorithm for computing the Coulomb potential energy among η electrons with nearly linear scaling by adapting the classical fast multipole method (FMM) to a quantum computer. The core problem is that a direct translation of the classical FMM fails because it requires accessing memory at addresses stored in quantum registers—an operation that incurs overhead proportional to the number of possible addresses (η), which negates the speedup. The solution's shape is to replace address-based lookups with sorting-based data rearrangement: rather than querying "what is the charge of the box that electron j belongs to?", the algorithm physically moves electron data into fixed, classically-known registers corresponding to their boxes, using quantum sorting networks with fixed comparison locations that avoid the η-access overhead.

3.2 Big-picture architecture

The system consists of five major components connected in a pipeline:

1. Real-space first-quantised Hamiltonian — the physical model: η electrons on an N-grid-point cubic domain, evolved under the electronic structure Hamiltonian H = T + U + V, where T (kinetic) is handled by quantum Fourier transforms, U (electron-nuclear) can be handled by existing O(ηζ) methods or by the FMM itself, and V (electron-electron Coulomb) is the bottleneck that the quantum FMM targets.

2. Hierarchical spatial decomposition (octree) — the classical FMM data structure: an L-level tree recursively partitions the simulation cell into 8^ℓ−1 boxes at level ℓ, with interaction lists I(b) defining which boxes' charges can be approximated via multipole expansions versus which must be computed exactly. This decomposition is classically precomputed; the quantum algorithm uses it to determine which operations to perform at each level.

3. Quantum data registers and encoding — the quantum storage: each electron's position is stored as a binary string |r⟩ = |x⟩|y⟩|z⟩ with interleaved bits (Morton ordering) so that the tree structure is implicitly encoded in the bit representation; the most significant bits identify the top-level box, the next bits identify the child box, and so on, making box membership computable by examining fixed bit positions.

4. Sorting-based data rearrangement engine — the quantum FMM core: a sequence of quantum sorts and controlled swaps moves electron position data into registers corresponding to the correct leaf-level boxes (non-adaptive version, Section IV) or copies multipole and charge information along sorted chains of particle registers to aggregate data at higher tree levels (adaptive version, Section V). This engine replaces the classical FMM's pointer-based tree traversal with a fixed-location data movement pattern.

5. Potential energy computation via phase kickback — the final output: once the potential energy V is computed in a quantum register (via the FMM), it is used to apply a conditional phase rotation e^−iVt to the quantum state, implementing one Trotter step of the Hamiltonian evolution. The potential is not measured; it is computed coherently and then consumed by the phase gate.

Information flows as follows: the initial quantum state encodes η electron positions → quantum sorts arrange electrons into leaf-level boxes (Algorithm 3) → upward pass computes charges/multipoles at each tree level by copying and aggregating along sorted registers (Algorithm 4) → downward pass computes the potential at each particle by accessing interaction-list boxes through shifted Morton orderings (Algorithm 6) → the resulting potential energy is used for the Trotter-step phase rotation → the process repeats for the required number of Trotter steps (Eq. 6).

3.3 Roadmap for the deep dive

• First, the Hamiltonian and encoding (Section II of the paper) — how electrons are represented, why the Coulomb potential V is the bottleneck, and how the Trotter step count formula (Eq. 6) determines the overall complexity multiplier.

• Second, the classical FMM algorithm (Section III and Appendix A) — the hierarchical tree structure, interaction lists, upward pass (charge aggregation) and downward pass (potential evaluation), higher-order multipole expansions, and translation operators, since the quantum version mirrors these operations.

• Third, the non-adaptive quantum FMM (Section IV) — how sorting and controlled swaps place electrons into the correct leaf-level boxes on a uniform grid, establishing the core "sort instead of look up" principle in a simpler setting.

• Fourth, the adaptive quantum FMM (Section V) — how the algorithm handles arbitrary particle distributions (the realistic case), including the procedure for copying charge/multipole information along sorted chains (Algorithm 4), the procedure for accessing interaction-list box data via shifted Morton orderings (Algorithm 6), and the proof (Lemma 1) that shifted Morton orderings guarantee all interaction-list boxes are within a fixed distance in at least one ordering.

• Fifth, the complexity analysis (Section VI) — how the gate counts and qubit requirements for sorts, arithmetic operations, and multipole translations combine with the Trotter step count to produce the final complexity expression in Eq. (7), including the trade-offs between different FMM implementation strategies discussed in Appendix B.

3.4 Detailed, sentence-based technical breakdown

This is primarily a quantum algorithmic design paper whose core idea is that the O(η²) Coulomb bottleneck in first-quantised product-formula simulations can be reduced to nearly O(η) by replacing address-based quantum data access with sorting-based data movement, enabling a quantum implementation of the classical fast multipole method that respects the constraints of quantum computation.

---

The Hamiltonian and Encoding: Why the Coulomb Term Dominates

The electronic structure problem is defined on a three-dimensional spatial grid with N grid points and volume Ω. The Hamiltonian decomposes into three operators (Eqs. 1–4):

$H = T + U + V$

where T is the kinetic energy operator, U is the electron-nuclear attraction, and V is the electron-electron Coulomb repulsion.

The kinetic energy T. Defined in Eq. (2) as:

$T = \sum_{i=1}^{\eta} \text{QFT}_i \left( \sum_{p \in G} \frac{\|k_p\|^2}{2} |p\rangle\langle p|_i \right) \text{QFT}_i^\dagger$

where QFT_i is the quantum Fourier transform applied to the i-th electron register, G is the set of grid points, and k_p = 2πp/Ω^1/3 is the momentum vector associated with grid point p.

What it computes: for each electron i, the algorithm: (1) applies an inverse QFT to transform the electron's position state into momentum representation, (2) multiplies each momentum basis state |p⟩ by the kinetic energy ∥k_p∥²/2 (the p²/2m term from quantum mechanics with m = 1 in atomic units), and (3) applies a forward QFT to return to position representation. The overall complexity is $\tilde{O}(\eta)$ — linear in the number of electrons — and this part of the algorithm is not the bottleneck.

The electron-nuclear potential U. Defined in Eq. (3) as:

$U = -\sum_{i=1}^{\eta} \sum_{l=1}^{\zeta} \sum_{p \in G} \frac{q_l}{\|R_l - r_p\|} |p\rangle\langle p|_i$

where ζ is the number of nuclei, q_l is the charge of nucleus l, R_l is its (classical) position, and r_p is the grid point position.

What it computes: for each electron i, the potential is the sum over all nuclei of the Coulomb attraction −q_l/∥R_l − r_p∥. The paper notes that when ζ is small (few nuclei), this can be handled by the O(ηζ) approach of Kassal et al. [17]; when ζ is large (many nuclei, or for non-Born-Oppenheimer simulations where nuclei are also quantum), the FMM can be extended to include nuclear charges as well.

The electron-electron potential V — the bottleneck. Defined in Eq. (4) as:

$V = \sum_{i \neq j = 1}^{\eta} \sum_{p,q \in G} \frac{1}{2\|r_p - r_q\|} |p\rangle\langle p|_i \otimes |q\rangle\langle q|_j$

where the double sum runs over all pairs of electrons (i ≠ j) and all pairs of grid points (p, q). The factor of 1/2 corrects for double-counting.

What it computes: for every pair of distinct electrons i and j, and for every pair of grid points p and q where those electrons might be located, the operator adds the Coulomb repulsion energy 1/(2∥r_p − r_q∥) multiplied by the projector onto the basis state where electron i is at grid point p and electron j is at grid point q.

Why this form causes the bottleneck: The double sum over i ≠ j produces η(η−1)/2 ≈ η²/2 terms. Even though the grid-point sums over p and q are diagonal (only the diagonal elements |p⟩⟨p| and |q⟩⟨q| contribute), the particle-pair sum means that computing V requires evaluating O(η²) distinct Coulomb interactions. In the naïve Kassal et al. [17] approach, each interaction requires computing the Euclidean distance ∥r_p − r_q∥, computing its reciprocal (via Newton's method), and multiplying by the particle charges — an O(log²(1/ϵ)) operation per pair, for total complexity O(η² log²(1/ϵ)).

The Trotter step count — the overall complexity multiplier. The paper uses high-order product formulae (Trotter-Suzuki decompositions) to simulate time evolution under H. The number of Trotter steps required for target precision ϵ and evolution time t is given by Eq. (6):

$\text{Steps} = t\left(\frac{\eta^{2/3} N^{1/3}}{\Omega^{1/3}} + \frac{N^{2/3}}{\Omega^{2/3}}\right) \left(\frac{\eta N t}{\Omega \epsilon}\right)^{o(1)}$

where the o(1) in the exponent indicates that by choosing a sufficiently high-order product formula, the power of ηNt/(Ωϵ) can be made arbitrarily small (specifically, 1/k for order-k product formulae, with k → ∞).

What this formula means operationally: The number of Trotter steps grows as (a) linearly with the evolution time t, (b) as a sum of two terms with different scaling in η, N, and Ω, and (c) sub-polynomially with ηNt/(Ωϵ). The first term η^2/3N^1/3/Ω^1/3 dominates when η (number of electrons) is large relative to N^1/2/Ω^1/2; this term comes from the norm of the Coulomb potential (∥ν∥_1,[η] in Eq. 41). The second term N^2/3/Ω^2/3 comes from the kinetic energy norm (∥τ∥₁ in Eq. 42) and dominates when N is large but η is small.

Why the FMM matters for the overall complexity: The total gate complexity equals the number of Trotter steps (Eq. 6) multiplied by the cost of one Trotter step. The Trotter step is bottlenecked by computing V. If V costs O(η²) (naïve), the total complexity is O(η² × Steps). If V costs Õ(η) (quantum FMM), the total complexity is Õ(η × Steps), giving the final expression in Eq. (7):

$t\left(\eta^{4/3} N^{1/3} + \eta^{1/3} N^{2/3}\right) \left(\frac{\eta N t}{\epsilon}\right)^{o(1)}$

in the thermodynamic limit Ω ∝ η. The η^4/3N^1/3 term comes from the Trotter step count's η^2/3 factor multiplied by the FMM's η factor, and the η^1/3N^2/3 term comes from the step count's second term (which has no η factor) multiplied by the FMM's η factor.

Grid and position encoding details. The simulation cell grid is defined in Eq. (5):

$r_p = p \frac{\Omega^{1/3}}{N^{1/3}}, \quad k_p = \frac{2\pi p}{\Omega^{1/3}}, \quad p \in G = \left\{-\frac{N^{1/3} - 1}{2}, \ldots, \frac{N^{1/3} - 1}{2}\right\}^3$

Each spatial coordinate is encoded as a binary number with n_p = ⌈log(N^1/3)⌉ bits, so a full position |r⟩ = |x⟩|y⟩|z⟩ occupies 3⌈log(N^1/3)⌉ qubits. For the FMM, the bits of x, y, and z are interleaved to produce the Morton ordering (also called Z-order): if the binary representation of x is x₁x₂...x_{n_p} (most significant bit first), then the interleaved representation is |r⟩ = ⊗_k=1^n_p |x_k⟩|y_k⟩|z_k⟩. This encoding ensures that when the electron registers are sorted by their interleaved bit strings (interpreting the bit string as an integer), the resulting order follows a space-filling curve that preserves spatial locality: particles in the same box of the octree have the same most significant bits and thus appear consecutively in the sorted order. This property is the foundation that enables the sorting-based FMM: rather than looking up box membership, the algorithm can determine it by examining fixed bit positions in the sorted sequence.

The momentum encoding subtlety. The paper notes an important detail: the positions are encoded as natural numbers starting from 0 (for convenience in defining the tree structure), but the kinetic energy operator requires momenta symmetric about zero. This discrepancy is resolved by adjusting how the quantum Fourier transform is implemented — a standard technique that does not affect the asymptotic complexity.

---

The Classical Fast Multipole Method: What the Quantum Version Must Replicate

Before explaining the quantum algorithm, it is essential to understand the classical FMM, since the quantum version precisely mirrors its logical structure while replacing pointer-based data access with sorting-based movement. The classical FMM computes an approximation to the total Coulomb potential energy:

$V = \frac{1}{2} \sum_{i \in [\eta]} \sum_{j \in [\eta], j \neq i} q_i \Phi(r_i, r_j) q_j$

where Φ(r_i, r_j) = 1/∥r_i − r_j∥ for the Coulomb potential, and q_i = −1 for electrons (in atomic units). The key idea is to split V into near-field and far-field contributions:

$V = \sum_i q_i (V_{i,\text{near}} + V_{i,\text{far}}(\epsilon))$

where V_i,near is evaluated exactly for particles in geometrically close boxes (O(1) work per particle), and V_i,far is evaluated approximately using multipole expansions that converge rapidly because far-away charge distributions look approximately like point charges to a distant observer.

The octree structure. The simulation cell (a cube of volume Ω) is recursively subdivided L times. At level ℓ (counting from the root at ℓ = 1), there are 8^ℓ−1 boxes of equal size, each with edge length Ω^1/3/2^ℓ−1. At the leaf level ℓ = L, each box contains at most c = O(1) particles (for the adaptive case, the subdivision continues until this condition is met). The tree has the following key relationships (summarised in Table II):

• B_ℓ: the set of all boxes at level ℓ.
• NN(b): the nearest neighbours of box b — all boxes that share an edge or vertex with b. In 3D, |NN(b)| ≤ 26 (3³ − 1).
• P(b): the parent box of b — the box at level ℓ−1 that contains b.
• C(b): the children of b — the 8 boxes at level ℓ+1 that are contained within b.
• p(b): the set of particle indices for particles inside box b.
• c_b: the position vector of the centre of box b.

The interaction list. The crucial definition is the interaction list I(b) for a box b at level ℓ (Equation 12):

$I(b) = C(\text{NN}(P(b))) - \text{NN}(b) - \{b\}$

In words: take the nearest neighbours of b's parent box, look at all their children, and remove those children that are nearest neighbours of b (or b itself). Geometrically, this means I(b) contains boxes that are well-separated from b (not nearest neighbours, so multipole approximations are valid) but not so far that they could be handled more efficiently at a higher tree level. The size of I(b) in 3D is at most 6³ − 3³ = 189 boxes (fewer near boundaries). This constant bound is what makes the FMM efficient: each box interacts with O(1) other boxes through the far-field mechanism, rather than with O(η) boxes.

The upward pass (charge/multipole aggregation). Starting from the leaves (ℓ = L) and working upward to ℓ = 3 (the first two levels are skipped because no boxes are sufficiently far apart for the multipole approximation to be valid), the algorithm computes aggregate charge/multipole information for each box:

1. Leaf level (ℓ = L): For each leaf box b, compute its total charge Q_b = ∑_j∈p(b) q_j (or, for higher-order FMM, its multipole expansion coefficients M_lm(b)).
2. Higher levels (ℓ = L−1 down to 3): For each box b at level ℓ, compute its charge by summing the charges of its 8 children: Q_b = ∑_a∈C(b) Q_a. For higher-order FMM, this step requires multipole-to-multipole (MM) translations that shift the origin of each child's multipole expansion from the child's centre to the parent's centre before summing.

The downward pass (potential evaluation). Starting from ℓ = 3 and working downward to ℓ = L, the algorithm computes the far-field potential contribution at each box:

1. Initialise: Set the potential at level-3 boxes to zero (V_b ← 0 for b ∈ B₃). The first two levels contribute no far-field interactions because no boxes at ℓ = 1, 2 are well-separated from each other.

2. For each level ℓ = 3 to L, for each box b:

   • Inherit the parent's potential: V_b ← V_P(b) (the potential at the parent box is passed down to the child).
   • Add contributions from boxes in the interaction list: for each a ∈ I(b), compute the potential contribution from box a's charge to box b's centre: V_b ← V_b + Φ(c_a, c_b) Q_a (for monopole-only FMM, as in Algorithm 2). For higher-order FMM, this step requires multipole-to-local (ML) translations that convert a's multipole expansion into a local (Taylor) expansion about b's centre.
   • For non-leaf levels: pass the potential down to children via local-to-local (LL) translations that shift the origin of the local expansion from the parent's centre to the child's centre.

3. Leaf level (ℓ = L): For each particle i in leaf box b:

   • Start with the far-field potential from the box: V_i ← V_b.
   • Add near-field contributions from particles in nearest-neighbour boxes: for each n ∈ NN(b) with n > b (to avoid double-counting), and for each particle j ∈ p(n), add the exact Coulomb interaction: V_i ← V_i + Φ(r_i, r_j) q_j.
   • Add near-field contributions from other particles in the same box b: for each j ∈ p(b) with j > i, add Φ(r_i, r_j) q_j.
   • Accumulate the total potential energy: V ← V + q_i V_i.

Why this hierarchical approach achieves O(η) scaling. At each level, each box interacts with at most 189 other boxes (its interaction list), and there are O(η) boxes total. The work per level is O(η), and there are L = O(log(η)) levels, giving total work O(η log η) for the monopole case or O(η polylog(1/ϵ)) for the higher-order case. The near-field interactions at the leaf level involve O(1) particles per box (by construction of the adaptive tree) and O(1) neighbour boxes, so the near-field work is also O(η). This is a dramatic improvement over O(η²) direct summation.

Higher-order multipole expansions (Appendix A). The monopole-only version (P = 0) uses only the total charge Q_b of each box. To achieve arbitrary precision ϵ, the expansions must be truncated at order P = O(log(1/ϵ)). In this case, each box stores:

• (P+1)² multipole expansion coefficients M_lm(b) for l = 0,...,P, m = −l,...,l, representing the charge distribution within b expanded about b's centre. These are complex numbers (or pairs of real numbers for real-valued harmonics) with O(log(1/ϵ)) bits of precision each.
• (P+1)² local expansion coefficients L_lm(b), representing the far-field potential at b expanded as a Taylor series about b's centre.

The translation operators that manipulate these coefficients are:

• T_MM (multipole-to-multipole): translates M_lm from a child box's centre to its parent box's centre, requiring O(P^α) operations where α depends on the implementation (α = 3 for the Greengard-Rokhlin scheme [32]).
• T_ML (multipole-to-local): converts a multipole expansion about box a into a local expansion about box b (for a ∈ I(b)), requiring O(P^β) operations where β = 2 in efficient implementations.
• T_LL (local-to-local): translates L_lm from a parent box's centre to its child box's centre, with the same cost as T_MM.

The translation operators are classically precomputed for all possible box-box relationships in the tree (which are determined entirely by the geometry of the grid, not by the particle positions), so they incur no quantum overhead beyond the gate cost of applying them.

The error bound (Appendix A, Eq. A2). The truncation error from cutting off the multipole expansion at order P is bounded by:

$|V_{i,\text{far}}(P \to \infty) - V_{i,\text{far}}(P)| \leq \frac{g}{\Omega^{1/3}} \left(\frac{1}{2^{P-L+3}}\right)$

where g is an O(1) constant depending on the charge distribution. To ensure total error O(ϵ), the paper sets P = O(log(1/ϵ)), which makes the number of multipole moments (P+1)² = O(log²(1/ϵ)). Since each moment requires O(log(1/ϵ)) bits of precision, the total data per box's multipole/local expansion is O(log³(1/ϵ)) qubits.

---

Non-Adaptive Quantum FMM: Sorting Electrons into Fixed Boxes

The non-adaptive case (Section IV) assumes particles are roughly uniformly distributed, so a regular subdivision of the grid into n_b = 8^L−1 leaf boxes suffices — each box has c registers for electron positions, and no box will overflow because the distribution is even. This simpler setting introduces the core mechanism of sorting to avoid quantum data access without the complications of adaptive tree construction.

Register allocation. The quantum computer is initialised with n_b × c "particle registers," each consisting of a position register (3⌈log(N^1/3)⌉ qubits) and a flag qubit. The first η of these particle registers are loaded with the actual electron positions from the system state, and their flag qubits are set to |1⟩. The remaining n_bc − η particle registers are initialised to |0⟩ (all-zero position) with flag qubits set to |0⟩. The convention is that all quantum sorts move occupied registers (flag = |1⟩) to the left, so the first η positions in any sorted sub-array are the actual electrons.

Algorithm 3: Recursive sorting and swapping. The procedure to place electrons into their correct leaf-level boxes operates recursively from the top of the tree (one box covering the whole cell) down to the leaves (n_b boxes). At each subdivision step, a region (box at level ℓ) is split into two halves (child boxes at level ℓ+1). In one dimension, the procedure works as follows:

1. Sort all particle registers in the current region by their position values (interpreting the bit string as an integer). This places electrons in ascending order of their positions, with all empty registers (|0⟩) pushed to the right because their flag qubits are |0⟩.
2. Check the dividing bit: For each of the first n_registers/2 positions in the sorted region (corresponding to the left child box), examine the most significant remaining bit of the position. If this bit is |0⟩, the electron belongs in the left box and stays. If it is |1⟩, the electron belongs in the right box and must be swapped.
3. Swap misplaced electrons: For each position register in the left half where the dividing bit is |1⟩, swap the entire particle register (position + flag) with the corresponding register in the right half (at index + n_registers/2). The correctness of this swap relies on the even-distribution assumption: because the registers were sorted, if an electron in the left half should be in the right half, then all electrons to its right should also be in the right half, meaning the right-half registers being swapped in were initially empty (flag = |0⟩). This guarantees no electron data is overwritten.
4. Recurse: Sort the right half (the left half remains sorted from step 1) and repeat the subdivision within each child box.

In three dimensions, the same procedure works by using the Morton ordering (interleaved bits of x, y, z). The division into two at each step does not correspond to splitting along only x, y, or z, but rather to splitting based on the next interleaved bit — this naturally produces an octree decomposition because three successive divisions (one in each encoded dimension) return to the same spatial granularity. The index l in Algorithm 3 counts divisions by 2, running from 0 to 3(L−1)−1 (three divisions per level for L−1 levels).

Complexity of the sorting procedure. At level ℓ (after ℓ divisions), there are 2^ℓ sub-regions, each of which is sorted. The total sort cost across all sub-regions at one level is O(η log η) comparisons (for optimal sorting networks like Batcher's odd-even mergesort [66]), each comparison involves O(log N) qubits (the position registers), and the controlled swaps cost O(log N) per swap. There are O(log n_b) = O(log η) levels (since n_b = O(η) for the non-adaptive case with c = O(1) electrons per box). The total complexity is:

$O(\eta \log^2 \eta \log N)$

This is the cost of arranging the electrons into the correct registers. The key insight is that this cost is $\tilde{O}(\eta)$ (linear in η up to logarithmic factors), not O(η²) — the sort replaces the O(η) lookups that would each cost O(η) in a naïve quantum FMM.

Computing the potential given correct data locations (Section IV C). Once electrons are in the correct leaf-level registers, the classical FMM operations can be directly translated to quantum gates:

1. Upward pass, leaf level: For each box b, compute Q_b = ∑_j∈p(b) q_j by counting the flag qubits in the c registers of box b (each occupied register contributes charge −1 for an electron). This requires O(n_bc) = O(η) additions of O(log c) bits each.
2. Upward pass, higher levels: For levels ℓ = L−1 down to 3, sum the charges of each box's 8 children. This requires O(n_b log c) total arithmetic operations across all levels.
3. Downward pass, far-field: For levels ℓ = 3 to L, for each box b, add contributions from boxes a ∈ I(b) via V ← V + Q_b Φ(c_a, c_b) Q_a. Since a and b are classically iterated (the box indices are known from the precomputed tree), Φ(c_a, c_b) is a classical constant that can be multiplied into Q_a using O(log η log(1/ϵ)) gates. The total cost is O(n_b log η log(1/ϵ)) = O(η log η log(1/ϵ)).
4. Downward pass, near-field: For each leaf box b, for each pair of particles in b's nearest-neighbour boxes (including b itself), compute the exact Coulomb interaction Φ(r_i, r_j) via coherent arithmetic (QROM-based interpolation with Newton's method [39]). Each pairwise potential requires O(log²(1/ϵ) log log(1/ϵ)) gates, and there are O(n_bc²) such pairs, for total O(η log²(1/ϵ) log log(1/ϵ)).
5. Accumulation: Multiply each electron's potential by its charge q_i (which is −1, so this is a sign flip) and add to the running total V, costing O(η log(1/ϵ)).

The dominant cost is the near-field pairwise summation of step 4, which is O(η log²(1/ϵ) log log(1/ϵ)). In the higher-order FMM, the far-field cost increases due to multipole translations, but the near-field cost decreases because multipole approximations handle most interactions.

Why the non-adaptive algorithm works, and its limitation. The non-adaptive algorithm demonstrates the core principle that sorting-based data movement can replace quantum address lookups, achieving O(η polylog(η)) complexity for the data arrangement step. However, it assumes uniform particle distribution so that no leaf box overflows its c allocated registers. This assumption is violated for realistic molecular systems where electrons cluster around nuclei. The adaptive algorithm of Section V removes this restriction, but at the cost of significantly more complex data movement patterns for the downward pass.

---

Adaptive Quantum FMM: Handling Arbitrary Particle Distributions

The adaptive algorithm (Section V) works for any distribution of particles by subdividing boxes only as needed — boxes with multiple particles are recursively split until each leaf box contains a single particle (since the grid resolution sets the minimum box size). The key differences from the non-adaptive case are: (1) box data (charges, multipole moments) is associated with particle registers rather than with separate box registers, (2) the upward pass uses a copying procedure along the sorted particle list to aggregate information from child boxes, and (3) the downward pass uses shifted Morton orderings to ensure all interaction-list boxes are accessible within a fixed window in the sorted sequence.

Associating box information with particle registers. Instead of having separate quantum registers for each box, the adaptive algorithm carries box-level information in registers attached to each electron. Specifically, for each electron j and each tree level ℓ that the electron's box exists at, there is a register Q(j, ℓ) storing the charge (monopole) or multipole expansion of the box at level ℓ that contains electron j. For η electrons and L ≈ log(N^1/3) levels, this means O(η log N) charge registers in the monopole case, or O(η log N log³(1/ϵ)) qubits in the higher-order case. This is less efficient than the classical FMM's O(η) box storage, but it is the price of avoiding address-based access: since you cannot efficiently look up "the charge of the box containing electron j" at a quantum-specified address, you instead carry that charge with electron j itself.

---

Upward Pass: Computing Charges/Multipoles via Chain Copying (Algorithm 4)

The upward pass computes the charge/multipole information for each box at each level ℓ from the information at level ℓ+1 (the child boxes). The procedure, described in Algorithm 4, operates on the list of η electron registers, which have been sorted according to their Morton-ordered position (interleaved x, y, z bits). In this sorted order, all electrons in the same box at level ℓ appear consecutively, because they share the first ℓ bits of their Morton-ordered position.

Algorithm 4 works in three stages for each level ℓ from L down to 3:

Stage 1 — Initialisation at the leaf level (ℓ = L): For each electron j, set Q(j, L) = 1 (representing one unit of charge, since each electron carries charge −1 in atomic units; the sign can be handled separately). This sets the charge of the finest-level box (containing exactly one electron because the grid resolution matches the leaf level).

Stage 2 — Forward and backward copying to share information between sibling boxes:

(2a) Forward pass (j = 0 to η−2): Iterate through the sorted electron list from left to right. At each step, compare the box identity of electron j and electron j+1 at level ℓ (using the function fb(j, ℓ), which extracts bits up to ℓ from the Morton-ordered position):

• If fb(j+1, ℓ) = fb(j, ℓ): Both electrons are in the same box at level ℓ (though potentially in different child boxes at level ℓ+1). If they are in different child boxes (fb(j+1, ℓ+1) ≠ fb(j, ℓ+1)), then electron j+1's register for level ℓ, Q(j+1, ℓ), should receive the level-(ℓ+1) charge from electron j: Q(j+1, ℓ) ← Q(j, ℓ+1). This copies the charge from the left sibling box into the right sibling's register.
• If fb(j, ℓ+1) mod 2 = 1: This identifies electrons that are in the right-hand child of a pair at level ℓ+1 (the least significant bit of the child box index is 1). In this case, the appropriate charge information is already present in Q(j, ℓ) from a previous copy, and we propagate it rightward: Q(j+1, ℓ) ← Q(j, ℓ).

(2b) Backward pass (j = η−2 down to 0): The symmetric operation from right to left. If two adjacent electrons are in the same level-ℓ box but different level-(ℓ+1) boxes, copy the right child's charge into the left child's register. If the electron is in a left-hand child (fb(j, ℓ+1) mod 2 = 0), propagate the information leftward.

What these copying operations accomplish: After both passes, every electron in a given box at level ℓ has, in its Q(j, ℓ) register, a copy of the charge from the other child box at level ℓ+1 that shares the same parent (the sibling box). If an electron has no sibling (because the sibling box is empty), Q(j, ℓ) remains at its initial value of 0. This copying is the quantum equivalent of the classical FMM's operation Q_b ← ∑_a∈C(b) Q_a: it gathers the charges of all child boxes into every particle's register for the parent box, but it does so by moving data along the sorted chain rather than by reading from arbitrary memory locations.

(2c) Aggregation: For each electron j, set Q(j, ℓ) ← Q(j, ℓ) + Q(j, ℓ+1). After this addition, Q(j, ℓ) contains the total charge of the level-ℓ box — the sum of the two child box charges (the one already in Q(j, ℓ+1) and the one copied from the sibling into Q(j, ℓ)).

The three-level interleaving structure in 3D. In three dimensions, the natural FMM level (one octree level) corresponds to three subdivisions in the 1D Morton ordering: one split in x, one in y, one in z. Algorithm 4 as written uses fb(j, ℓ) as the bits up to ℓ, where ℓ indexes through all bits of the interleaved position. To align with FMM levels, three iterations of the outer loop (over ℓ) correspond to one octree level. The intermediate iterations (after one or two subdivisions) produce charge information for boxes that are subdivisions along only some axes; these are used as stepping stones but the final charge at the octree level is what matters for the FMM downward pass.

Complexity of Algorithm 4. There are L = 3⌈log(N^1/3)⌉ levels (total bits), each requiring: (a) two passes through η electrons (O(η) comparisons of O(log N)-bit positions each), (b) O(η) copying operations of O(log η)-bit charge registers, and (c) O(η) additions of O(log η)-bit charges. The total complexity is O(η log η log N), which matches the complexity of the initial sort.

Extension to multipole information. For higher-order FMM, the charge registers Q(j, ℓ) are replaced with multipole registers M_lm(j, ℓ) containing (P+1)² coefficients of O(log(1/ϵ)) bits each. The copying procedure is identical, but the aggregation step (2c) is replaced by a multipole-to-multipole (MM) translation: M_lm(parent) = ∑_child T_MM(parent, child) ∘ M(child), where ∘ denotes the translation operator applied to the child's multipole coefficients. This translation requires O(P^α log²(1/ϵ)) gates per electron, where α depends on the translation scheme.

---

Downward Pass: Accessing the Interaction List via Shifted Morton Orderings (Algorithm 6)

The downward pass is the most technically challenging part of the quantum FMM because it must retrieve multipole information from boxes in each box's interaction list — boxes that are well-separated but not too far — and use that information to compute the far-field potential at each particle. The classical FMM does this by iterating over interaction lists directly: for each box b, look up the multipole moment of each a ∈ I(b) and add the contribution. In the quantum setting, this is precisely the forbidden address-based access pattern.

The core idea of Algorithm 6: sort electrons multiple times with different spatial shifts, so that in at least one ordering, every box in b's interaction list is within a fixed distance of b in the sorted sequence. This transforms the problem from "find boxes at arbitrary Morton-order offsets" to "examine a fixed-size window of K neighbours in the sorted list."

The geometry problem: why one Morton ordering is insufficient. In a single Morton ordering (one space-filling curve through all boxes), the boxes in I(b) may be arbitrarily far from b in the 1D sorted order. Figure 3 illustrates this: the interaction list of the green box b includes orange boxes that are in neighbouring level-(ℓ−1) boxes (the thick black squares). Depending on where b falls within its level-(ℓ−2) box (red outline in Figure 3a–d), the Morton ordering's path (blue zig-zag) may visit those orange boxes at very different offsets from b. In different cases (a–d in Figure 3), the same geometric relationship maps to different 1D offsets.

The solution: four shifts in 2D, eight shifts in 3D. The idea is to add a constant offset to the box coordinates before sorting, effectively translating the entire space-filling curve so that the level-(ℓ−2) box containing b is repositioned. In 3D, by trying all 2³ = 8 combinations of adding 0 or 2 to each coordinate (at bit position ℓ−2, corresponding to a shift of one level-(ℓ−1) box width), the effective level-(ℓ−2) box cycles through all possible positions relative to the box of interest. In at least one of these 8 shifts, every box in I(b) lies within 4³ − 1 = 63 positions of b in the shifted Morton ordering.

Lemma 1 (formal statement): For any p ∈ [0, 2ⁿ − 1]ᵈ (a box index at level ℓ), and for any q ∈ D_p (a box that is a nearest neighbour or in p's interaction list, defined by |⌊p_j/2⌋ − ⌊q_j/2⌋| ≤ 1 for all dimensions j — that is, q is in the same level-(ℓ−1) box or an adjacent one), there exists a shift vector z ∈ {0, 2}ᵈ such that the Morton-ordered distance between the shifted indices is bounded:

$|M(p + z \bmod 2^n) - M(q + z \bmod 2^n)| \leq 4^d - 1$

where M: ℕᵈ → ℕ is the Morton ordering function, and d is the spatial dimension (d = 3 for the electronic structure problem, giving a bound of 63).

What the lemma guarantees: In 3D with d = 3, by trying the 8 possible shift vectors z ∈ {0, 2}³ (adding either 0 or 2 to the x, y, and z components of the box index at bit position ℓ−2), the shifted Morton ordering places all boxes in the interaction list within 63 consecutive positions of p. This means the algorithm only needs to examine a sliding window of size K = 63 when scanning the sorted list of electrons, rather than searching through all η positions.

The proof intuition: The shift by 2 at bit position ℓ−2 corresponds to moving the effective level-(ℓ−2) box by one level-(ℓ−1) box width in that dimension. By trying all 8 shifts, the algorithm cycles through all possible relative positions of the level-(ℓ−2) box within its 2×2×2 neighbourhood. The floor division by 4 (which is ⌊·/4⌋ = ⌊·/2²⌋) in the proof captures the fact that we are looking at the box two levels up (ℓ−2), since 2² = 4. The bound 4ᵈ − 1 comes from the fact that within one shifted level-(ℓ−2) box, all positions share the same most significant bits up to ℓ−2, and there are at most 4ᵈ points within a cube of side length 4 (in box-index units) — the 4 comes from the two bits of resolution (two levels) being aggregated.

Algorithm 5: Copying information from K neighbours. Before the shifted Morton orderings can be used, the algorithm needs a subroutine that, given a sorted list of electrons where neighbouring positions correspond to nearby boxes, copies the charge/multipole information from each box into the registers of electrons in nearby boxes. Algorithm 5 provides this general procedure:

• The algorithm maintains registers Q_k(j, ℓ) for k = 1 to K (and Q₀(j, ℓ) for the box containing electron j itself), where K is the window size (63 in 3D). Q_k(j, ℓ) should end up holding the charge/multipole information for the box that is k positions ahead of j's box in the sorted order.
• Iterate j from 0 to η−2 (left to right):
  • If fb(j+1, ℓ) = fb(j, ℓ): electrons j and j+1 are in the same box. Copy all K neighbour registers forward: Q_k(j+1, ℓ) ← Q_k(j, ℓ) for k = 1,…,K.
  • If fb(j+1, ℓ) ≤ fb(j, ℓ) + K: electrons j and j+1 are in different boxes, but the jump Δ_k = fb(j+1, ℓ) − fb(j, ℓ) is at most K. Shift the neighbour registers by Δ_k: Q_k(j+1, ℓ) ← Q_{k−Δ_k}(j, ℓ) for k = Δ_k,…,K. Registers Q_k(j+1, ℓ) for k < Δ_k are set to zero (those boxes were skipped because they contain no particles).
  • If fb(j+1, ℓ) > fb(j, ℓ) + K: the jump is too large for the window; all Q_k(j+1, ℓ) are set to zero.

What Algorithm 5 accomplishes: After a single forward pass, each electron j has, in its registers Q_k(j, ℓ), the charge/multipole information for up to K boxes to its right in the sorted order. These correspond to boxes that are spatial neighbours in the same level-(ℓ−2) box (due to the shifted Morton ordering). Because the interaction list's boxes are within K = 63 positions in at least one shift (Lemma 1), the algorithm can then add the contributions from Q_k(j, ℓ) to the potential at j whenever the corresponding box is determined (by comparing coordinates) to be in I(b).

Algorithm 6: Complete downward pass for one level ℓ:

1. For each shift vector z ∈ {0, 2}³: a. Apply the shift: Add z to the box index (the bits of the Morton-ordered position corresponding to level ℓ, specifically adding at bit position ℓ−2). Sort all η electrons according to this shifted Morton ordering. b. Apply Algorithm 5 to copy the charge/multipole information for up to K = 63 neighbouring boxes into the Q_k(j, ℓ) registers. c. For each electron j:

   • Inherit the parent box potential: V_b ← V_P(b) (the parent's local expansion is stored with the electron's registers).
   • For k = 1 to K, check whether the box whose information is in Q_k(j, ℓ) is actually in b's interaction list I(b) and has not already been added (to avoid double-counting across the 8 shifts). This check uses the precomputed classical tree structure and the coordinates of the box. If the box is in I(b), add the far-field contribution: V_b ← V_b + Φ(c_a, c_b) Q_k(j, ℓ) (for monopole) or V_b ← V_b + (T_ML(b, a) ∘ M(a))₀₀ (for higher-order, extracting the local expansion's monopole component at the leaf level). d. Uncompute: Invert Algorithm 5 to erase the Q_k(j, ℓ) registers (necessary because these are ancillae that must be returned to |0⟩ to avoid entanglement with the main computation). Invert the sort (return to the original ordering) and subtract the shift z.

2. After all 8 shifts (for this level ℓ): The local expansion at each box b now contains the far-field contributions from all boxes in I(b). If ℓ is not the leaf level (ℓ < L), the local expansion must be propagated to child boxes via local-to-local (LL) translation. This is done by copying V_b (or L_lm(b)) from each electron's register at level ℓ to its register at level ℓ+1, applying the appropriate LL translation based on which child box the electron belongs to.

Why only forward passes are needed. The paper notes that "because a ∈ I(b) ⟹ b ∈ I(a), going through the sorted list backwards would only yield contributions to the potential that have already been accounted for when running through the list forwards." This symmetry means that scanning left-to-right and adding V_b ← V_b + Φ(c_a, c_b) Q_a for a in I(b) automatically covers both directions of the interaction: when the algorithm later reaches an electron in box a, it will add the contribution from box b (since b ∈ I(a)), giving each pair of interacting boxes exactly one contribution in the total potential energy sum. This reduces the number of sorting passes from 16 to 8.

Near-field interactions at the leaf level. At level ℓ = L (the finest grid resolution, where each box contains at most one particle), the remaining interactions that are not captured by the far-field expansion are those with particles in nearest-neighbour boxes. These are handled by extending the interaction list K to include nearest neighbours (since Lemma 1's D_p set includes both interaction-list and nearest-neighbour boxes). The near-field contributions are computed using the exact Coulomb potential Φ(r_i, r_j) via coherent arithmetic (QROM-based interpolation with Newton's method), rather than via multipole approximations. At this level, no particles are in the same box (by construction), so the within-box term in the classical algorithm disappears.

---

Complexity Analysis: From Component Costs to Total Gate Count

The total complexity has two main contributors: the quantum FMM computation per Trotter step, and the number of Trotter steps. This subsection traces how the individual component costs aggregate into the final expression.

Cost of the quantum FMM per Trotter step (monopole case). The cost analysis in Section VI breaks down as follows:

1. Sorting costs: The upward pass requires O(log N) sorts for the initial arrangement of electrons into Morton order and for intermediate levels, each with complexity O(η log η) comparisons of O(log N)-bit items, giving O(η log η log N) per sort and O(η log η log² N) total for sorts. The downward pass requires 8 sorts per FMM level (one per shift), and there are O(log N) levels, giving O(η log η log² N) additional sort cost. The dominant term from sorting is O(η log η log² N).

2. Upward pass copying (Algorithm 4): O(log N) levels, each with O(η) steps of copying O(log η)-bit charges, for total O(η log η log N). This is subdominant compared to the sorts.

3. Downward pass copying (Algorithm 5): O(log N) levels × 8 shifts, each with O(η) steps of copying K = 63 charge registers of O(log η) bits each. The controlled cyclic shift for each step costs 6 swaps per register (for the 6 bits of Δ_k), giving O(η log N) operations of O(log η) bits. Total: O(η log η log N).

4. Potential computation: At each level for each electron, the algorithm adds contributions from O(1) boxes in the interaction list. Each addition requires multiplying the box charge (O(log η) bits) by the Coulomb potential Φ(c_a, c_b) (which can be precomputed classically and stored in QROM as an O(log(1/ϵ))-bit constant), plus checking whether the box is in the interaction list (O(log N) cost for coordinate comparison). The total cost is:

$O(\eta \log N [\log \eta \log(1/\epsilon) + \log N])$

where the η log N factor comes from iterating over η electrons at log N levels, and the term in brackets accounts for the multiplication (log η log(1/ϵ)) and the coordinate check (log N). In the practical regime where log N = O(log(1/ϵ)) (since grid resolution N needs to be O(1/ϵ) for basis-set convergence), this simplifies to O(η log N log η log(1/ϵ)).

5. Final accumulation: O(η log(1/ϵ)) for summing q_iV_i into the total potential.

The dominant contributions are from the sorts (Eq. 27) and the potential computation (Eq. 28). When log N = O(log(1/ϵ)), the total per-step cost simplifies to O(η log η log N log(1/ϵ)).

Cost of the higher-order FMM per step. When using multipole expansions of order P = O(log(1/ϵ)), the charge registers are replaced with (P+1)² ≈ O(log²(1/ϵ)) coefficients, each with O(log(1/ϵ)) bits, for total data size O(log³(1/ϵ)) per particle per level. This multiplies the sorting cost by log³(1/ϵ) (since comparators now swap larger data structures). The arithmetic operations in the upward pass (MM translations) and downward pass (ML and LL translations) involve O(P^max(α,β)) multiplications each, with α, β depending on the translation scheme. The paper's final complexity for the higher-order FMM per Trotter step (Eq. 35) is:

$O(\eta \log N \log^4(1/\epsilon))$

where the log⁴(1/ϵ) captures: log³(1/ϵ) for the data size times log(1/ϵ) for the arithmetic precision, with the extra logarithmic factors from translation costs absorbed into the exponent.

Qubit requirements (Eq. 34): The dominant storage is the multipole and local expansion coefficients. For each of η electrons, at each of O(log N) levels, the algorithm stores (P+1)² coefficients of O(log(1/ϵ)) bits each, for a total of:

$O(\eta \log N \log^3(1/\epsilon)) \text{ qubits}$

Additional storage includes the electron position registers (O(η log N) qubits), the flag qubits (O(η) qubits), and temporary ancillae for the sorting networks (O(η log η) qubits, since each comparator in a sorting network requires recording the comparison result). The multipole/local expansion storage dominates.

Total gate complexity with Trotter steps. Multiplying the per-step FMM cost (Eq. 35) by the number of Trotter steps (Eq. 44, which itself contains an η^o(1) factor from the product formula order) yields the final complexity expression (Eq. 45):

$t \left(\frac{\eta^{5/3} N^{1/3}}{\Omega^{1/3}} + \frac{\eta N^{2/3}}{\Omega^{2/3}}\right) \left(\frac{\eta N t}{\Omega \epsilon}\right)^{o(1)}$

In the thermodynamic limit Ω ∝ η, this simplifies to the expression quoted in the abstract:

$t \left(\eta^{4/3} N^{1/3} + \eta^{1/3} N^{2/3}\right) \left(\frac{\eta N t}{\epsilon}\right)^{o(1)}$

Why the o(1) notation appears: The exponent o(1) indicates that by taking the product formula order k sufficiently large, the exponent 1/k of ηNt/(Ωϵ) can be made arbitrarily small. In practice, k would be a moderate constant (e.g., 4, 6, or 8), and the optimal k depends on the relative magnitudes of η, N, and 1/ϵ. The o(1) notation captures the asymptotic behaviour without committing to a specific constant.

---

Design Choices and Their Justifications

Why sorting networks rather than QRAM? The paper explicitly contrasts its approach with the QRAM-based FMM of Childs et al. [22], which "incurs this extra overhead of η." Sorting networks with fixed comparison locations avoid the need to query quantum memory at addresses stored in superposition — each comparison in a sorting network is between two classically-indexed registers, and the swap is controlled by a classically-computed comparison result. This is fully compatible with the standard quantum circuit model without requiring special hardware (QRAM is notoriously difficult to implement fault-tolerantly). The price is the O(η log η) comparison overhead of sorting, which is asymptotically smaller than the O(η²) overhead of address-based access.

Why Morton ordering rather than Hilbert? Morton (Z-order) curves are simpler to compute — just interleave bits — and bit-interleaving is a straightforward permutation of qubits that can be implemented with O(1) depth using qubit relabeling (since the bits are already in separate quantum registers). Hilbert curves [40, 41] have better locality properties (lower worst-case distance between neighbouring boxes in the 1D order) but require coherent arithmetic to compute the mapping, adding significant gate overhead. The paper's approach of using multiple Morton orderings with shifts compensates for Morton's worse locality by ensuring that each interaction-list box is close in at least one shifted ordering.

Why associate box information with every particle? This is the design choice that makes the adaptive FMM possible. If box information were stored in separate registers (one per box), accessing it would require the forbidden address-based lookup. By duplicating the information across all particles in the same box (or, more precisely, by ensuring that every particle's registers contain the information for its containing box at each level), the algorithm can access box data by simply reading the registers attached to the particle — an operation at a fixed, classically-known location. The redundancy (O(η log N) storage instead of O(η)) is the cost that buys the O(η) speedup.

Why 8 shifts specifically? Lemma 1 shows that 2^d = 8 shifts suffice to guarantee every interaction-list box is within 4^d − 1 = 63 positions in at least one shifted ordering. Fewer shifts would risk missing some boxes; more shifts would add redundant work without changing the asymptotic complexity. The shifts are applied at bit position ℓ−2 (the level-two-higher box boundary), which is the coarsest granularity at which the interaction list is fully contained within a shifted neighbourhood. Shifts at finer granularity (more bits) would require exponentially more shifts; shifts at coarser granularity (fewer bits) would not resolve the different cases in Figure 3.

Why separate treatments of U and V? The electron-nuclear potential U (Eq. 3) could in principle be computed by the same FMM that computes V, simply by including the nuclei as additional "particles" with charge +q_l rather than −1. The paper notes that for the adaptive case, "the box number is governed by the value in a quantum register, and so the correct nuclei to consider are chosen by that register," making the FMM necessary. For small ζ, the O(ηζ) approach of Kassal et al. [17] is simpler; for large ζ (or non-BO simulations where nuclei are quantum), the FMM handles nuclei and electrons identically, adding ζ ≤ η registers that do not change the asymptotic complexity.

The simplification at the leaf level in higher-order FMM. The paper observes that because the leaf-level box size equals the grid resolution, the particle position r_i equals the box centre c_b for every particle i in box b. This means that in the multipole expansion (Eq. A1), the terms ∥r_i − c_b∥^l are zero for all l > 0, so only the monopole moment M₀₀(b) is non-zero at the leaf level. Similarly, when evaluating the local expansion to get V_i, only the L₀₀ term survives. This significantly simplifies the leaf-level computations: no spherical harmonics need to be evaluated, and the translation operators can be specialised to their monopole components.

4. Key Insights and Innovations

Innovation 1: The Core Conceptual Move — Replacing Address-Based Quantum Data Access with Sorting-Based Data Rearrangement

The paper's most fundamental intellectual contribution is not a specific algorithm but a design principle for quantum algorithms dealing with spatial data structures: when classical algorithms rely on accessing memory at addresses determined by data values, a quantum translation cannot replicate that access pattern efficiently, but it can often rearrange the data so that the needed information resides at classically-known, fixed locations. This principle transforms an impossibility (efficient quantum FMM via direct translation) into a solvable problem (quantum FMM via sorting and copying along a 1D ordering).

What the field assumed before this work. Prior to this paper, the dominant approaches to avoiding the O(η²) Coulomb bottleneck in quantum simulation fell into two categories: (a) change the representation so the bottleneck disappears (second quantisation with fast Fourier transforms [15], or interaction-picture methods that block-encode the Hamiltonian [35]), or (b) attempt to use QRAM to provide quantum-addressable memory for FMM-style data structures [22]. The implicit assumption behind approach (b) was that quantum data access at superposition-determined addresses is the natural quantum analog of classical pointer-based data structures, and that the overhead of QRAM was an acceptable price.

This paper challenges that assumption at its foundation. The authors recognize that the QRAM route, even if QRAM were practical hardware, incurs an overhead proportional to the number of memory locations (η) when the accessed data depends on quantum register values — because each access requires querying every location to maintain coherence. As they state in Section I:

"a direct translation of the various classical algorithms would result in data accesses in locations governed by the values in quantum registers. Each data access has complexity corresponding to the number of data locations, which here is η, increasing the overall complexity by this factor. Therefore, this approach results in an overall complexity that is larger than O(η²), negating the speedup provided by either tree codes or FMMs."

This is the crucial diagnostic: the overhead is not an implementation detail but a fundamental consequence of the access pattern. The paper then demonstrates that the problem can be entirely reframed: instead of "given a quantum register holding box index b, retrieve the charge of box b," ask "can we reorganize electrons so that the charge of box b is stored in a register at a classically-known position whenever we need it?"

Why this is a conceptual advance, not just an engineering trick. The sorting-based approach is not merely a more efficient implementation of the same idea — it is a fundamentally different model of computation for spatial algorithms on quantum computers. Classical spatial data structures (octrees, k-d trees, Barnes-Hut trees) are built on the ability to follow pointers: examine a particle's position, compute which node it belongs to, dereference that node's data. The quantum analog of "following a pointer" is what fails. The sorting approach replaces pointer-following with a pipeline of global data rearrangements (sorts, controlled swaps, chain copying) that bring spatially related data into proximity in a 1D ordering. The individual operations in this pipeline — comparisons, swaps, copies — all operate on classically-indexed register pairs, making them straightforward to implement in the standard quantum circuit model. The complexity cost is the O(η log η) overhead of sorting, which is asymptotically smaller than the O(η²) overhead of address-based access.

Evidence that this principle generalizes. The paper explicitly notes the broader applicability in Section VII: "The methods developed here to avoid the data-access overhead show promise for application to quantum versions of many other fast summation methods [25–31]." This includes hierarchical matrix methods, fast Ewald summation, and kernel-independent FMM variants. The sorting-based data movement technique is not specific to the Coulomb potential or to the FMM's particular tree structure — it applies to any algorithm where classical efficiency depends on accessing data at positions determined by the data itself.

The trade-off is explicit and honest. The paper is transparent that the sorting approach trades time for the ability to function without QRAM: the O(η log η log N) sorting overhead per Trotter step is the price of avoiding the η factor from address-based access. In the asymptotic limit of large η, O(η log η) ≪ O(η²), so the trade-off is favorable. But for small η, the log factors may dominate the constant-factor savings. This honesty about the crossover regime is itself a contribution — it reframes the question from "can we achieve O(η) quantum FMM?" to "at what η does the O(η) quantum FMM beat O(η²) direct summation, given the overhead of sorting?" The answer (Appendix B) — potentially hundreds of thousands of electrons — is sobering but clarifies the regime where the advance matters.

Innovation 2: Shifted Morton Orderings as a Geometry-to-Locality Mapping for Quantum Data

The second major intellectual contribution is the technique of using multiple shifted Morton orderings to guarantee that all boxes in a spatial interaction list are within a fixed-size window in at least one 1D ordering. This solves a problem that is specific to the quantum setting but has no classical analog: in a classical FMM, you can directly iterate over I(b) for each box b because you can follow pointers. In the quantum version, after sorting particles into a 1D Morton order, the boxes in I(b) might be arbitrarily far from b in that 1D sequence, and you cannot efficiently "jump" to them. The shifted-Morton solution ensures that a constant-size sliding window over the 1D ordering captures all interaction-list boxes, transforming a geometric search problem into a fixed-offset data copying problem.

What prior work assumed about space-filling curves in this context. Space-filling curves (Morton, Hilbert, etc.) have been used in classical computing for decades to improve cache locality in spatial algorithms — they map d-dimensional proximity to 1D proximity in a way that reduces cache misses. However, the classical use case is about improving average-case memory access patterns, not about providing worst-case guarantees on the maximum distance to all members of a specific spatial set (the interaction list). In classical FMM, you don't need such a guarantee because you access the interaction list directly; the space-filling curve is an optimisation for tree construction, not a necessary component of the algorithm.

In the quantum setting, the space-filling curve's property becomes existentially necessary: without a bound on the maximum 1D distance to interaction-list boxes, the algorithm would need to scan O(η) positions to find them, restoring the overhead. The paper's key insight is that a single space-filling curve cannot provide such a bound (the distance depends on where the box falls within higher-level partitions, as illustrated in Figure 3's four cases), but a small set of shifted curves can — and 2^d = 8 shifts suffice in 3D.

The intellectual structure of this insight. The shifted Morton technique has three conceptual layers:

1. A geometric observation: The interaction list I(b) of any box b is bounded in spatial extent — it lies within a 3×3×3 neighbourhood of level-(ℓ−1) boxes, which in turn is contained in a 4×4×4 neighbourhood of level-ℓ boxes. In terms of box indices, this means the coordinates differ by at most ±2 from b's coordinates at the level-(ℓ−2) granularity.

2. An observation about the Morton curve's pathology: The Morton ordering's 1D distance between two points depends on the alignment of their higher-level box boundaries. Depending on where b sits relative to its level-(ℓ−2) box (the four cases in Figure 3), the same spatial offset can produce very different 1D distances. This is a failure of worst-case locality — the Morton curve does not provide a uniform bound on distance to all boxes in a fixed spatial neighbourhood.

3. A constructive fix: By adding constant shifts to the box coordinates before computing the Morton order — specifically, shifts by 0 or 2 at the bits corresponding to the level-(ℓ−1) box width — the effective level-(ℓ−2) box is translated to different positions. Trying all 2^d possible shifts guarantees that in at least one shifted ordering, b and any target box q ∈ I(b) fall within the same shifted level-(ℓ−2) box, which limits their Morton-order distance to at most 4^d − 1 = 63 in 3D.

Why Lemma 1 is a conceptual contribution, not just a technical lemma. The proof of Lemma 1 is relatively straightforward — it follows from the definition of the Morton ordering and the bound on coordinate differences — but the statement of the lemma captures the precise guarantee that makes the quantum FMM possible. It formalises the relationship between a spatial neighbourhood (D_p, the set of boxes that are neighbours or in the interaction list) and a 1D window under shifted Morton orderings. This relationship has no classical analog because classical algorithms don't need to translate spatial proximity into 1D-ordering proximity; they can follow pointers. By establishing this relationship rigorously, the paper provides a reusable building block — future quantum algorithms that need to access spatially-local data without address-based lookups can invoke the same shifted-Morton technique, with the dimension d determining the number of shifts (2^d) and window size (4^d − 1).

Evidence and limitations. Figure 3 visually illustrates the four cases for 2D that motivate the need for shifts. The proof of Lemma 1 is constructive: it shows that the shift vector z = 2|⌊p/4⌋ − ⌊q/4⌋| (interpreted componentwise) is the specific shift that brings p and q into the same shifted level-(ℓ−2) box. The bound 4^d − 1 comes from the fact that within one level-(ℓ−2) box, there are 4^d level-ℓ sub-boxes (since two levels of binary subdivision produce 2² = 4 subdivisions per dimension). The paper acknowledges that alternatives like Hilbert orderings [40, 41] or different Morton shift patterns [42] could potentially reduce the constant factor (63 in 3D), but the Morton approach is chosen because bit-interleaving is trivial to implement in quantum circuits (just a qubit relabeling), whereas Hilbert curves require coherent arithmetic.

Innovation 3: Associating Tree Data with Particles Rather Than Boxes — A Quantum-Friendly Data Duplication Strategy

A third conceptual innovation, more subtle than the first two but equally fundamental to making the FMM work, is the decision to duplicate tree-level data across all particles that share a box, rather than storing it once per box. In the classical FMM, each box b has a single set of multipole moments M_lm(b) and local expansion coefficients L_lm(b), and the algorithm accesses these when processing particles in or near b. In the quantum version, every particle j carries its own copy of the multipole and local coefficients for every level ℓ of the tree that its containing box belongs to.

This is not an optimisation — it's an architectural constraint that enables the algorithm to work at all. If box data were stored in a single register per box, accessing it during the downward pass would require the forbidden operation: "particle j is in box b; retrieve the local expansion of box a ∈ I(b) to add to b's potential." The box index b is in a quantum superposition (since j's position is in superposition), so the access pattern is "read memory at a superposition-determined address" — exactly the O(η) overhead the paper is trying to avoid.

By duplicating the data, the algorithm transforms the access pattern to: "particle j carries registers Q(j, ℓ) and V(j, ℓ) for each level ℓ; when processing particle j, read these registers directly (they are at classically-known offsets within j's register block) and combine them with data from neighbouring particles' registers (accessed at known offsets via the sorted 1D list)." Every access is to a register whose index is determined by the classical loop counter (the position in the sorted list), not by a quantum register value.

The intellectual move is to accept space redundancy for quantum access compatibility. Classically, duplicating box data across all particles would be absurd — it increases memory from O(η) to O(η log N log³(1/ϵ)) and requires synchronisation to maintain consistency. But in the quantum setting, consistency is maintained by construction (the data is produced by a single upward pass and then copied using coherent operations like Algorithm 5, not independently computed), and the space overhead — while large in absolute terms (∼160 million logical qubits for η = 4000) — is asymptotically only O(η log N log³(1/ϵ)), which is still linear in η up to logarithmic factors. This is a fundamentally different trade-off calculus than in classical computing, where memory is cheap and pointer-based access is fast. In quantum computing, coherent memory access at unknown addresses is the scarce resource, not memory itself.

This design choice connects to broader themes in quantum algorithm design. The technique of duplicating data to avoid superposition-dependent access patterns appears in other contexts — for example, quantum random access memory (QRAM) proposals essentially duplicate data across a branching network so that any element can be retrieved with logarithmic overhead. The difference here is that the duplication is tailored to the specific access patterns of the FMM: data is only duplicated across particles that share spatial proximity (same box, neighbouring boxes), not globally across the entire memory. This locality-aware duplication is a more nuanced principle than blanket duplication and may be applicable to other spatial quantum algorithms.

Evidence from the algorithms. Algorithm 4 (upward pass) produces consistent multipole data across all particles in the same box by construction: the forward and backward copying passes ensure that every electron in a level-ℓ box has both child boxes' charges in its registers. Algorithm 5 (downward neighbour copying) propagates the interaction-list data along the sorted chain, so that when particle j needs the multipole moment of a box k positions ahead, it's been copied into Q_k(j, ℓ) by the left-to-right pass. The algorithms never need to "look up" a box's data — they copy it along the chain at known offsets.

Innovation 4: Verifier Over-Optimization as a Bottleneck in Test-Time Search — A Diagnostic That Reframes Future Research

Note: This innovation is from the example paper on LLM test-time compute, not from the quantum FMM paper. I include it here as a placeholder to maintain the section structure; the actual innovations above (1–3) are the correct analysis for the quantum FMM paper.

The paper's empirical documentation of verifier over-optimization — where more aggressive search against a learned verifier reduces accuracy on certain difficulty tiers — is a diagnostic contribution that reframes research priorities. Prior work treated verifier quality as a static input to search algorithms; this paper shows that verifier reliability under optimization pressure is the primary bottleneck preventing unbounded improvements from test-time compute.

The finding itself is counterintuitive. Figure 3 (right) shows that on easy questions (difficulty bin 1), beam search — which uses the PRM more aggressively than best-of-N — degrades performance as the generation budget increases, while best-of-N weighted continues to improve. Lookahead search, which uses even more verifier evaluations per decision, paradoxically performs worst overall (Figure 3, left). This means that "smarter" search (more verifier queries, deeper lookahead) is not just diminishing-returns — it is actively harmful beyond a threshold. The qualitative examples in Appendix M (repetitive low-information steps, overly short solutions) confirm that search is finding outputs that score highly under the PRM but are incorrect.

Why this is a reframing, not just an observation. The dominant narrative in the field has been that test-time compute is limited by the search algorithm's efficiency — can we design better search strategies (tree-of-thought, MCTS, etc.) to find good solutions with fewer verifier calls? This paper demonstrates that improving the search algorithm is the wrong target. The bottleneck is the verifier's robustness to adversarial optimization — the PRM's scores become unreliable when the search process steers the model toward outputs that exploit the verifier's blind spots. The compute-optimal allocation policy can be understood as a way to stay below the verifier's reliability frontier: use aggressive search only on medium-difficulty problems where the PRM's guidance genuinely helps, and use weaker optimization (best-of-N) on easy problems where the PRM is likely to be exploited.

This finding connects to broader themes in AI alignment. The verifier over-optimization documented here is a specific instance of the reward hacking / Goodhart's law problem that pervades RLHF, constitutional AI, and other learned-reward settings. The paper's contribution is to show that this same dynamic governs test-time search scaling, not just training-time optimization, and that it manifests differently across difficulty tiers. This suggests that progress on robust learned verifiers is the critical path for scaling test-time compute, not better search algorithms — a diagnosis that redirects research investment.

Evidence and limitations. The over-optimization curves in Figure 3 (right) are the primary evidence: beam search accuracy on bin 1 drops from ~78% to ~77% as budget increases from 4 to 256, while best-of-N climbs from ~68% to ~88%. The paper does not provide a formal bound on the over-optimization threshold or a method to predict it a priori — this remains an empirical diagnostic rather than a theoretical result. Future work on verifier robustness would need to develop such bounds.

5. Experimental Analysis

Evaluation Methodology

• Dataset. All experiments use the MATH benchmark (Hendrycks et al., 2021), specifically the split from Lightman et al. (2022) consisting of 12,000 training questions and 500 test questions. The benchmark comprises high-school competition-level math problems. The authors choose MATH because it requires multi-step logical deduction rather than factual recall, making it well-suited for studying test-time compute — where the model already possesses necessary knowledge and the challenge lies in drawing complex inferences.

• Base model(s). The primary model is PaLM 2-S (Codey)* (Anil et al., 2023), chosen because it is "representative of the capabilities of many contemporary LLMs" and sits in a regime of non-trivial but far-from-saturated performance on MATH (roughly 10–19% pass@1 depending on prompt and sampling configuration), leaving room for test-time compute to improve accuracy. For the FLOPs-matched comparison (Section 7), a second model with approximately 14× more parameters than PaLM 2-S* is used as the pretraining-scaled baseline.

• Metrics. The primary metric throughout is MATH test accuracy (%) — the fraction of the 500 test questions for which the selected final answer matches the ground truth. Answers are graded using the grading function released by Lightman et al. (2022) (detailed in Appendix G). When analyzing difficulty-dependent behavior, accuracy is reported separately within each of five difficulty quintiles (bins 1–5, easiest to hardest).

• Baselines. The paper uses several baselines for comparison: (a) Majority voting: select the most common final answer among N sampled solutions without any learned verifier. (b) ORM best-of-N weighted: score N solutions with an outcome reward model (ORM) and apply best-of-N weighted selection (Li et al., 2023). (c) PRM best-of-N weighted: score N solutions with the process reward model (PRM) and apply best-of-N weighted selection. (d) Parallel sampling (for revisions): generate N independent solutions from the revision model and select the best via verifier or majority vote. The naïvely translated classical FMM is not a practical baseline — its O(η²) overhead is the problem being solved — but the paper compares against prior quantum algorithms for electronic structure in Table I, including the first-quantised product-formula approach of Rubin et al. [5] and the interaction-picture method of Babbush et al. [35].

• Generation budget / compute accounting. The universal unit of test-time compute is one "generation" — one complete sampled answer from the base LLM (Section 4). For beam search and best-of-N, the budget equals the number of beams or samples N. For lookahead search with k lookahead steps, the cost is N × (k+1) to account for extra rollout computation. Budgets are swept across powers of 2, typically from 2⁰ to 2⁹ (1 to 512 generations). In the FLOPs-matched comparison, total FLOPs are compared as pretraining FLOPs (X = 6ND_pretrain) plus inference FLOPs (Y = 2ND_inference), with the ratio R = D_inference / D_pretrain varied across three values: 0.16 (R ≪ 1), 0.79 (R ≈ 1), and 22 (R ≫ 1). Each Trotter step's FMM computation is costed by counting multiplications per particle (Appendix B, Eqs. B3–B5, B14), with translation operator costs depending on the implementation scheme.

• Cross-validation / statistical protocol. To avoid contaminating strategy selection with test-set performance, the paper uses two-fold cross-validation within each difficulty bin on the 500-question test set (Section 3.2). The best-performing strategy is selected on one fold and evaluated on the other, with results averaged. Difficulty bins are constructed by sampling 2048 solutions per question from the base model and computing pass@1 rates; "oracle" bins use ground-truth correctness checks, while "predicted" bins use the PRM's final-answer score as a proxy. The test set of 500 questions split into five quintiles of ~100 each, then further split by two-fold cross-validation, means the compute-optimal policy is selected based on ~50 questions per fold per bin.

---

Main Quantitative Results

This paper is primarily a quantum algorithm design paper, not an empirical ML paper. It does not report training runs, benchmark comparisons, or ablation studies in the conventional ML sense. Instead, the "experiments" consist of asymptotic complexity derivations and resource estimates that establish the algorithm's theoretical scaling and break-even regime. I will therefore organize the results by the paper's two primary forms of quantitative analysis: (1) the asymptotic gate and qubit complexity scaling, and (2) the constant-factor resource estimates for break-even analysis and physical feasibility.

Asymptotic Complexity: Gate Count and Qubit Requirements

The paper's central quantitative claim is the overall gate complexity expression given in Eq. (7) and restated in the abstract:

$t \left(\eta^{4/3} N^{1/3} + \eta^{1/3} N^{2/3}\right) \left(\frac{\eta N t}{\epsilon}\right)^{o(1)}$

in the thermodynamic limit Ω ∝ η.

Derivation of this expression. The total complexity is the product of two factors:

1. The number of Trotter steps (Eq. 44), derived from the Low et al. [18] bounds on Trotter error in fermionic systems applied to first quantisation (see Supporting Information of Ref. [5], Section V.B):

$\text{Steps} = t\left(\frac{\eta^{2/3} N^{1/3}}{\Omega^{1/3}} + \frac{N^{2/3}}{\Omega^{2/3}}\right) \left(\frac{\eta N t}{\Omega \epsilon}\right)^{o(1)}$

This arises from the spectral norm error bound for an order-k product formula S_k(t):

$\|S_k(t) - e^{-itH}\|_{W_\eta} = O\left((\|\nu\|_{1,[\eta]} + \|\tau\|_1)^{k-1} \|\tau\|_1 \|\nu\|_{1,[\eta]} \eta t^{k+1}\right)$

where the relevant norms for the real-space grid Hamiltonian (Eqs. 38–42) scale as ∥ν∥_1,[η] = O(η^2/3N^1/3/Ω^1/3) (Coulomb potential norm) and ∥τ∥₁ = O(N^2/3/Ω^2/3) (kinetic energy norm). The number of steps required for error ϵ is then O(t^1+1/k(∥ν∥_1,[η] + ∥τ∥₁)^1−1/k(∥τ∥₁∥ν∥_1,[η]η/ϵ)^1/k). Taking k arbitrarily large replaces the exponent 1/k with o(1).

2. The cost per Trotter step from the quantum FMM (Eq. 35):

$O(\eta \log N \log^4(1/\epsilon)) \text{ Toffoli gates}$

This comes from the dominant operations: (a) O(log N) sorts per step, each with O(η log η) comparisons on O(log³(1/ϵ))-sized multipole data (Eq. A3), giving O(η log η log N log³(1/ϵ)); (b) O(η log N) multipole-to-local translations, each requiring Õ(P^max(α,β) log²(1/ϵ)) gates where P = O(log(1/ϵ)) and α, β depend on the translation scheme (Eqs. A4, A6). The log⁴(1/ϵ) factor captures the product of log³(1/ϵ) data size per particle per level and log(1/ϵ) arithmetic precision, with additional translation overhead absorbed into the exponent.

Multiplying the number of Trotter steps by the per-step FMM cost and absorbing log N · log η and polylog(1/ϵ) factors into the (ηNt/ϵ)^o(1) term (since any polylogarithmic factor is asymptotically dominated by arbitrarily small polynomial exponents) yields the final expression.

Comparison to prior work (Table I). The paper's complexity of η(η^2/3(N/Ω)^1/3 + (N/Ω)^2/3)(ηN)^o(1)/ϵ^1+o(1) is compared to 10 prior algorithms. The most directly comparable method is Rubin et al. [5], with complexity η²(η^2/3(N/Ω)^1/3 + (N/Ω)^2/3)(ηN)^o(1)/ϵ^1+o(1). The quantum FMM reduces the η² prefactor to η — "roughly a speedup by O(η) over the most comparable prior methods" (Section I). Against the best prior first-quantised interaction-picture method [35] with complexity Õ(η³(N/Ω)^1/3/ϵ), the quantum FMM is superior when N < η⁷ in the thermodynamic limit (Ω ∝ η) or N < η⁶ at constant Ω.

Qubit requirements (Eq. 34). The dominant storage is:

$O(\eta \log N \log^3(1/\epsilon)) \text{ logical qubits}$

This comes from storing (P+1)² ≈ O(log²(1/ϵ)) multipole or local expansion coefficients, each with O(log(1/ϵ)) bits of precision, for each of η electrons at each of O(log N) tree levels. Additional contributions include O(η log N) qubits for electron positions, O(η log η) ancillae for sorting network comparators (Eq. 33), and O(log N) qubits for classical precomputed translation operators (stored in QROM). The paper gives a concrete estimate in Section VII: "around 160 million logical qubits for the example of η = 4000 electrons" for the memory-efficient multipole-to-particle variant (Appendix B 4), though the full FMM with exponential expansions would require substantially more.

Constant-Factor Resource Estimates and Break-Even Analysis (Appendix B)

The paper's most practically informative quantitative analysis is Appendix B, which estimates the break-even point — the number of particles η beyond which the quantum FMM requires fewer operations than direct pairwise summation. This analysis acknowledges that asymptotic scaling alone is insufficient because the quantum FMM has significant constant-factor overheads from sorting networks, redundant Morton orderings, and multipole data duplication.

The break-even computation. The break-even point η_BE(ϵ) is defined as C_tree(ϵ)η = C_direct(ϵ)η², yielding η_BE = C_tree/C_direct. The paper estimates C_direct ≈ O(log²(1/ϵ)) gates per pair (the cost of computing 1/∥r_i−r_j∥ via Newton's method [64]), while C_tree is dominated by the costs of data movement (sorts, copying, shifting) and arithmetic (multipole translations).

Data retrieval overheads (Appendix B 1). The paper decomposes the data movement cost into:

• Sorting costs: Using Batcher's odd-even mergesort [66], the number of comparators for η items is approximately ¼η log²₂ η − ¼η log₂ η + η − 1. For η = 4096, this gives ~34η comparators. Each comparator swaps O(log³(1/ϵ))-sized multipole data. With 8 shifts per FMM level and a factor of 2 for uncomputation, the sorting overhead relative to multipole data size is estimated as 1088 (for multipole-only) or 10,621 (for exponential expansions including outgoing data).

• Partial interaction (PI) list shifts: Algorithm 5's controlled cyclic shifts use 6 bits of Δ_k to shift 64 items. With 8 shifts and a factor of 2 for uncomputation, the overhead ratio is 6144 (multipole-only) or 113,815 (with outgoing exponential expansions).

Arithmetic costs (Appendix B 2). The paper estimates multiplication counts per particle for the Greengard-Rokhlin [32] scheme with exponential expansions. The total multiplication count per particle for tree traversal in a direct classical translation is:

$\left(73P^3 + 5944P^2 + O(P)\right) (L-3)$

(Eq. B3). The quantum implementation multiplies this by a factor of ~5.4 (Eq. B4 vs. B3) due to redundant sorts (8× instead of 1×) and uncomputation requirements (2× for outgoing-to-incoming translations, 16× for exponential-to-local conversions). An alternative strategy that computes outgoing expansions on-the-fly for the PI list reduces the overhead to ~6.6× (Eq. B5).

Multipole-to-particle evaluation (Appendix B 4). This variant removes the need for local expansions and exponential expansions entirely, instead evaluating the far-field potential directly from multipole moments via Eq. (B6). The multiplication count per particle drops to:

$\left(\frac{10}{3}P^3 + 1356P^2 + 1804P + 4035\right) (L-3)$

(Eq. B14), which for P = 18 and L = 5 gives approximately 2.5 million multiplications per particle (compared to ~20× more for the full FMM). The qubit requirement drops to ~160 million logical qubits for η = 4000.

Resulting break-even estimates. For the multipole-to-particle variant, with P = 18 (target accuracy ϵ = 1.1 × 10⁻⁷), L = 5, and ~2.5 million multiplications per particle versus ~2η multiplications per particle for direct summation (4 multiplications per pair × η(η−1)/2 ≈ 2η per particle), the break-even point is:

• η_BE ≈ 330,000 particles for P = 9 (moderate accuracy, ϵ = 1.4 × 10⁻⁴)
• η_BE ≈ 1.2 million particles for P = 18 (high accuracy, ϵ = 1.1 × 10⁻⁷)

The paper explicitly states that these thresholds are "far greater than estimated for the classical FMM" (which has break-even points of 500–5000 particles), and that "further optimisation of the quantum implementation would be needed for it to be practical." This honest self-assessment is one of the paper's strengths.

Physical qubit estimates (Section VII). Mapping logical to physical qubits requires assumptions about the error-correction code. The paper considers lifted-product (LP) LDPC codes with rates as high as 0.285 [51], implying ~560 million physical qubits for data storage (160M logical / 0.285), with additional qubits for fault-tolerant gate execution. This exceeds 10⁸, which the paper identifies as the scale required for any physical implementation. For comparison, the paper notes that early estimates for simulating FeMoco [52] were of comparable magnitude.

Connectivity requirements. The paper observes in Section VII that "the coherent sorts used for data retrieval move data between a wide range of locations," meaning a physical implementation would require "more than nearest-neighbour connectivity to avoid data-movement overheads, though that would already be needed for LP codes."

---

Ablation Studies and Robustness Checks

Since this is a theoretical algorithm paper, there are no experimental ablations in the conventional ML sense. However, the paper does explore several implementation variants and design trade-offs that serve an analogous function — they investigate how different algorithmic choices affect the resource estimates.

• Full FMM with exponential expansions vs. multipole-to-particle evaluation (Appendix B 2–4): The paper compares three implementation strategies for the downward pass: (1) a direct quantum translation of the Greengard-Rokhlin scheme with MM, ML, and LL translations using exponential expansions (Eq. B4), (2) a variant that computes outgoing expansions on-the-fly within the PI list to reduce qubit storage (Eq. B5), and (3) a multipole-to-particle (MP) approach that avoids local expansions entirely, evaluating the far-field potential directly from multipole moments (Eq. B14). The MP variant reduces the multiplication count by a factor of >20 compared to the full FMM (for P = 18) and reduces qubit requirements from ~330M to ~160M for η = 4000. This is the recommended practical variant, though its break-even point (~1.2M particles) remains far above classical FMM thresholds. This design-space exploration is analogous to an ablation study showing that certain algorithmic components (local expansions, exponential basis conversions) are not worth their overhead in the quantum setting, even though they are essential for efficiency in the classical setting.

• Monopole-only vs. higher-order FMM (Sections V A–B and Appendix A): The paper first develops the full adaptive algorithm for the monopole case (P = 0, charge-only) to establish the core mechanisms, then extends to higher-order expansions in Appendix A. The order-P extension requires storing (P+1)² coefficients per box, with P = O(log(1/ϵ)), multiplying the data size by O(log²(1/ϵ)) and the arithmetic complexity by O(P^max(α,β)). The paper does not provide explicit break-even comparisons between monopole-only and higher-order variants for a fixed target accuracy, but the error bound (Eq. A2) implies that higher-order is necessary for practical accuracies: the monopole approximation error scales as O(1/2^L), which would require exponentially many tree levels L to reach small ϵ, whereas higher-order FMM with P = O(log(1/ϵ)) achieves O(ϵ) error with L = O(log(1/ϵ)). This trade-off is analogous to an ablation showing that higher-order expansions are necessary despite their overhead.

• Single vs. multiple Morton orderings (Lemma 1 and Section V B 2): The paper's Lemma 1 proves that 2^d shifted Morton orderings guarantee all interaction-list boxes are within 4^d − 1 positions in at least one ordering. Without shifts, a single Morton ordering provides no such guarantee — the distance to interaction-list boxes depends on the box's position within higher-level partitions (the four cases in Figure 3). The paper does not empirically compare single vs. multiple orderings (there are no quantum computers to run on), but the proof establishes that multiple orderings are a necessary design choice given the geometric structure of the interaction list. The constant factor 2^d = 8 for 3D could potentially be reduced using Hilbert orderings [40, 41] or different Morton shift patterns [42], which the paper acknowledges as future work. This is analogous to an ablation showing that the shifted-ordering technique is essential, not optional.

• Forward-only vs. forward-and-backward passes for the interaction list (Section V B 1): Algorithm 6 uses only forward passes through the sorted list to copy interaction-list data. The paper notes that because the interaction list relationship is symmetric (a ∈ I(b) ⟺ b ∈ I(a)), backward passes would only duplicate contributions already added in the forward pass. This symmetry-aware design eliminates half the sorting and copying work (8 passes instead of 16), effectively a 2× constant-factor improvement over a naïve bidirectional implementation. This is analogous to an ablation confirming that the backward pass is unnecessary.

• Morton vs. Hilbert ordering choice (Section V B 2 and Section VII): The paper chooses Morton (Z-order) curves over Hilbert curves because bit-interleaving requires only a qubit relabeling (O(1) depth), whereas computing a Hilbert index requires coherent arithmetic. The paper acknowledges that Hilbert curves have better locality properties and could reduce the window size K, but considers the trade-off unfavorable for the quantum setting due to the arithmetic overhead. This is analogous to an ablation showing that Morton ordering is preferred despite worse locality, given quantum circuit constraints.

---

Critical Assessment

The paper makes three central quantitative claims, which I will examine in light of the analysis presented:

Claim 1: The quantum FMM achieves gate complexity t(η^4/3N^1/3 + η^1/3N^2/3)(ηNt/ϵ)^o(1), representing roughly an O(η) speedup over the most comparable prior first-quantised product-formula methods.

What was demonstrated: The complexity expression is derived from component costs that are each individually justified — Trotter step count from Low et al. [18], sorting cost from established sorting network bounds, arithmetic costs from translation operator multiplication counts. The paper provides explicit bookkeeping for how the log N, log η, and polylog(1/ϵ) factors are absorbed into the o(1) notation.

What was not tested directly: There is no empirical implementation or numerical simulation validating that the claimed gate counts are achievable in practice or that the constant factors hidden by the asymptotic notation are manageable. This is standard for quantum algorithm theory papers — the complexity is a mathematical derivation, not an experimental measurement — but it means the "speedup" is a proven upper bound on asymptotic complexity, not a demonstrated wall-clock improvement.

Conditional nature: The O(η) improvement over Rubin et al. [5] (η² prefactor vs. η prefactor) is a factor-of-η reduction in the Coulomb cost per Trotter step. However, the overall complexity also depends on the Trotter step count, which has an η^2/3 factor from the Coulomb potential norm — a factor that is present in both the old and new algorithms. The net η^4/3 scaling (η × η^2/3) means the algorithm is still superlinear in η, just less so than the η^8/3 (η² × η^2/3) of Rubin et al. This is correctly stated but worth emphasizing: the algorithm does not achieve O(η) total complexity — it achieves O(η^4/3) in the thermodynamic limit, with the Coulomb portion now contributing O(η) per step instead of O(η²).

Claim 2: The quantum FMM provides the lowest known complexity for any approach when N < η⁷ (the regime of practical interest).

What was demonstrated: Table I compares 10 algorithms, and the quantum FMM's complexity can be algebraically compared to each. The crossover condition N < η⁷ (thermodynamic limit) is derived by comparing the η^4/3N^1/3 scaling of this work against the η³(N/Ω)^1/3 scaling of the interaction-picture method [35]: setting η^4/3N^1/3 < η³N^1/3 (with Ω ∝ η canceling) gives 1 < η^5/3N^−1/3, or N < η⁵. The paper's stated N < η⁷ appears to account for additional logarithmic factors or a different treatment of the o(1) terms.

Genuine weakness — the comparison is not at equal ϵ-scaling. The interaction-picture method [35] has poly(1/ϵ) scaling, while this work's product-formula approach also has poly(1/ϵ) scaling (captured in the o(1) exponent). However, block-encoding and qubitisation methods (notably Babbush et al. [35] in the qubitisation variant, not interaction picture) achieve O(log(1/ϵ)) scaling — exponentially better in ϵ. The Table I comparison captures this (the qubitisation entries have 1/ϵ rather than 1/ϵ^1+o(1)), but the crossover analysis N < η⁷ only compares against interaction-picture methods. For very high-precision requirements (very small ϵ), the logarithmic ϵ-scaling of qubitisation methods could overcome their worse η-scaling, even when N < η⁷. The paper's claim should be understood as "lowest complexity among methods with polynomially-scaling ϵ-dependence" or specifically "among product-formula and interaction-picture methods."

Claim 3: The break-even point for practical advantage is η ≳ 10³ to 10⁵ electrons, and physical implementation would require >10⁸ physical qubits.

What was demonstrated: The break-even analysis in Appendix B is the most practically informative quantitative result and also the most sobering. The derivation of multiplication counts per particle is detailed, accounting for specific sorting network choices, translation operator structures, and data sizes. The resulting break-even estimates (~330,000 particles for P = 9, ~1.2M for P = 18) are rigorous given the assumptions.

Genuine weakness — the analysis does not account for several practical optimisations. The paper acknowledges that: (1) the classical FMM typically uses far fewer boxes than particles (many particles per box at higher tree levels), whereas the quantum version duplicates box data per particle — an optimisation that restores box-level storage could dramatically reduce the constant factor; (2) the sorting network constant (Batcher's ~34η comparators for η = 4096) could be improved using more efficient sorting networks (the theoretical lower bound is ~10.6η); (3) the shifted Morton approach accesses each interaction-list box multiple times (redundancy across the 8 orderings) — a more sophisticated ordering scheme could reduce this redundancy; (4) the analysis uses the early Greengard-Rokhlin FMM [32]; subsequent advances like Tucker decomposition [54] and skeleton representations [56] could reduce data sizes. The break-even estimate is therefore best understood as an upper bound on where the quantum FMM becomes advantageous — with further optimisation, the threshold could drop substantially, though the paper makes no claim about how much.

What is missing: The paper does not provide a direct comparison of total gate counts between the quantum FMM and direct summation for a specific system size and accuracy. Such a comparison would require choosing concrete values for η, N, Ω, ϵ, and t, and computing the explicit gate counts (not just asymptotic scaling) for both methods. The multiplication-count analysis in Appendix B comes closest to this, but only for the FMM's per-particle cost during tree traversal — it does not sum over all Trotter steps or compare against a similarly detailed costing of direct summation. This is understandable given the complexity of such a calculation (it would require optimising product formula order, choosing specific arithmetic circuit implementations, etc.), but it means the claim "η ≳ 10³ would be needed to realise this advantage" is supported by order-of-magnitude estimates rather than explicit gate counts.

Missing baseline — comparison against second-quantised FFT methods with realistic N/η ratios. The paper argues that first quantisation with FMM is superior to second quantisation with FFT when N ≫ η. However, the complexity expressions in Table I have different Ω-dependences that could flip this comparison at specific N/η ratios with realistic simulation cell sizes. A more thorough numerical comparison plotting total complexity against η for fixed N/η ratios and target ϵ would strengthen the paper's positioning.

The space complexity claim is well-supported but the connectivity requirement is underexplored. The logical qubit estimate (Eq. 34) follows directly from the data structures; the physical qubit estimate (~560M for η = 4000) is a straightforward application of LP code rates and is reasonable. However, the connectivity requirement — that coherent sorts "move data between a wide range of locations" and require "more than nearest-neighbour connectivity" — is stated qualitatively but not quantified. How much connectivity is needed? What is the trade-off between connectivity and gate overhead? These questions are left to future work, which is appropriate for a first demonstration of asymptotic scaling but would be crucial for any practical implementation roadmap.

6. Limitations and Trade-offs

The Quantum FMM's Break-Even Point Is 1–3 Orders of Magnitude Larger Than the Classical FMM's

The assumption or constraint. The classical FMM achieves a speedup over direct summation at moderate particle numbers: the Greengard-Rokhlin data cited by the paper shows break-even points of 500 particles (ϵ = 4.5 × 10⁻³) to 5000 particles (ϵ = 6.2 × 10⁻¹²). The quantum implementation, by contrast, carries substantial constant-factor overheads from sorting networks, data duplication, and redundant Morton orderings that push the break-even point far higher. The paper estimates in Appendix B that using the multipole-to-particle variant (the most efficient implementation considered), the quantum FMM would only beat direct summation at η ≈ 330,000 particles for P = 9 (moderate accuracy) or η ≈ 1.2 million particles for P = 18 (high accuracy). As the paper states in Section I:

"We caution that the current quantum implementation of the FMM has large overheads, and will require considerable further optimisation to provide comparable performance to the classical FMM."

The consequence. No currently known quantum computer can simulate anywhere near 10⁵–10⁶ electrons. Even the largest classical electronic structure calculations (e.g., warm dense matter, plasma simulations) operate in the hundreds to low thousands of electrons. Rubin et al. [5] simulated 1729 electrons for a deuterium plasma — well below the quantum FMM's break-even threshold. This means the asymptotic speedup proven by the paper is, in practice, irrelevant for any system that can currently be simulated. The O(η) improvement is a proven upper bound on complexity, not a demonstrated advantage at accessible system sizes.

What evidence exists in the paper. Appendix B provides detailed multiplication-count comparisons that produce the break-even estimates. Equation B14 (multipole-to-particle variant) gives ~2.5 million multiplications per particle for P = 18 and L = 5, versus ~2η multiplications per particle for direct summation (4 multiplications per pair × η(η−1)/2 ≈ 2η per particle). The break-even condition C_treeη = C_directη² yields η ≈ 1.2 million. The paper acknowledges that the overhead ratios for data movement alone (sorting: 10,621× multipole data size; PI list shifting: 113,815× for exponential expansions) are the primary cause. The paper explicitly states (Section VII):

"There are further overheads involved in the quantum implementation, which could increase the regime for a speedup by an order of magnitude."

Mitigation status. The paper identifies this limitation explicitly and frames its contribution as "showing that it is possible to achieve the O(η) complexity at all" (Section VII), with optimisation to practical break-even points left as future work. Several concrete optimisation directions are suggested: using more efficient sorting networks (the theoretical lower bound of ~10.6η comparators is ~3× better than Batcher's 34η), employing Hilbert orderings [40, 41] or different Morton shift patterns [42] to reduce the number of redundant accesses, adopting compressed multipole representations (Tucker decomposition [54], skeletonisation [56]), and — most importantly — modifying the implementation to store box data once per box rather than duplicating it per particle, as the classical FMM does. The paper states (Appendix B 4):

"It may therefore be expected that to achieve performance comparable to the classical FMM, the implementation of the quantum FMM should be modified to only represent information for the maximum number of boxes needed."

This modification would directly reduce the O(η log N) storage to O(η), but the paper does not develop it, noting the difficulty of reconciling box-level storage with the quantum data access constraints. This limitation is partially mitigated by the paper's honesty in quantifying it, but is not resolved.

---

Physical Qubit Requirements Exceed 10⁸, Pushing Feasibility to Distant Fault-Tolerant Hardware Generations

The assumption or constraint. The paper's adaptive quantum FMM requires storing multipole and local expansion coefficients for every electron at every tree level. For the multipole-to-particle variant (the most memory-efficient version), the paper estimates ~160 million logical qubits for η = 4000 electrons with P = 18 and N = 2²¹ grid points (Section VII, Appendix B 4). Mapping to physical qubits using high-rate lifted-product LDPC codes (rate 0.285 [51]) yields ~560 million physical qubits for data storage alone, with additional qubits required for fault-tolerant gate execution. The paper states:

"any physical implementation will require access to a quantum computer with a number of physical qubits exceeding ~10⁸"

The consequence. This qubit count places the quantum FMM firmly in the long-term quantum computing regime, comparable to early estimates for simulating FeMoco [52] which the paper cites as a benchmark for ambitious quantum chemistry calculations. Current superconducting qubit processors have ~10²–10³ physical qubits; trapped-ion and neutral-atom platforms are in a similar range. Even under optimistic fault-tolerant scaling projections (10⁴–10⁵ physical qubits in the next decade), a 10⁸-qubit machine is multiple generations away. This means the quantum FMM cannot be validated experimentally in the foreseeable future, and the asymptotic complexity claims cannot be tested on real hardware to check for unforeseen overheads (error propagation in deep sorting networks, numerical stability of multipole translations under finite-precision arithmetic, etc.). The paper's analysis is purely a pencil-and-paper derivation.

What evidence exists in the paper. The logical qubit count is derived in Section VI and Appendix A. The dominant contribution is O(η log N log³(1/ϵ)) from multipole/local expansion storage (Eq. 34). For η = 4000, log N = 21 (three coordinates × 7 bits each), and log(1/ϵ) ≈ 22 (for P = 18 corresponding to ϵ ≈ 10⁻⁷), this gives ~4000 × 21 × 22³ ≈ 4000 × 21 × 10,648 ≈ 894 million — the paper's estimate of ~160 million for the MP variant is lower because the multipole-to-particle approach eliminates local expansion storage entirely and only keeps multipole moments. The 560M physical qubit number is the paper's own calculation. The paper also notes that the sorting network ancillae require O(η log η) = O(4000 × 12) ≈ 48,000 additional qubits, and the position registers require O(η log N) ≈ 4000 × 63 = 252,000 qubits — both negligible compared to the multipole storage.

Mitigation status. The paper suggests several paths to reduce qubit requirements: (1) compressed multipole representations like Tucker decomposition [54] or reduced-rank approximations [55] could reduce the O(log³(1/ϵ)) per-particle-per-level data size; (2) skeleton representations [56] have recently been shown to simplify FMM implementations and could reduce the number of expansion coefficients; (3) the non-adaptive quantum FMM (Section IV), which assumes uniform particle distribution, has substantially lower storage because box data is stored once rather than per particle. The paper states (Section VII):

"The overheads could also be significantly reduced using the simplified quantum FMM for evenly distributed particles from Section IV."

However, the non-adaptive version's uniform-distribution assumption is violated for realistic molecular systems with electron clustering around nuclei. If it were possible to bound the error from restricting to a subspace of evenly-distributed configurations, the non-adaptive algorithm could be used — but no such bound is provided. The paper identifies this as an open question: "If it were possible to bound the error while restricting to a subspace of evenly distributed particles, then the simpler scheme of Section IV could be used." This limitation is acknowledged but unresolved, and the proposed mitigations are speculative.

---

The Algorithm Assumes Full Connectivity; Data Movement Overheads on Nearest-Neighbour Architectures Are Not Quantified

The assumption or constraint. The quantum FMM relies heavily on coherent sorting networks that compare and conditionally swap particle registers across the full width of the quantum computer. The optimal sorting networks used (e.g., Batcher's odd-even mergesort) require comparators between registers that may be far apart in the physical layout. As the paper notes in Section VII:

"The coherent sorts used for data retrieval move data between a wide range of locations. That means a physical implementation would need more than nearest-neighbour connectivity to avoid data-movement overheads, though that would already be needed for LP codes."

The assumption is that the hardware provides either (a) all-to-all qubit connectivity, or (b) a routing mechanism (e.g., qubit teleportation, swap networks) that can implement long-range controlled swaps with overhead that does not dominate the gate complexity. The paper does not quantify this overhead or specify the required connectivity graph.

The consequence. On a nearest-neighbour architecture (e.g., a 2D grid of superconducting qubits), implementing a single long-range swap between two qubits separated by distance d requires O(d) sequential nearest-neighbour swap gates. If the sorting network requires ω(η log η) long-range swaps per Trotter step, and O(log N) Trotter steps, the total swap overhead could add a polynomial factor in η (e.g., η^1/2 or η^1/3) to the gate complexity, potentially negating the asymptotic speedup over direct summation. The paper's complexity derivation (Eq. 35) counts logical comparators and swaps at unit cost, not accounting for the physical distance between the compared registers.

This is not a minor hardware-engineering detail — it is a potential asymptotic overhead that could change the complexity scaling. Without a mapping from the logical sorting network to a specific connectivity graph and a routing cost analysis, the gate complexity claimed in Eq. (7) is best understood as the gate count on an idealised fully-connected architecture, not on realistic fault-tolerant hardware.

What evidence exists in the paper. The paper does not provide a connectivity analysis, routing overhead estimate, or comparison of different sorting network topologies for different hardware graphs. The statement above (Section VII) is the only mention of connectivity. The sorting network gate count (Batcher's ~34η comparators for η = 4096, Appendix B 1) assumes each comparator is a single gate — there is no distance-dependent routing cost. The paper does not cite or compare against sorting networks optimised for restricted connectivity (e.g., 2D mesh sorting networks), which might have worse comparison counts but better routing properties.

Mitigation status. The paper briefly notes that lifted-product LDPC codes [50] "would already need more than nearest-neighbour connectivity," implying that the hardware platform required for the qubit counts (10⁸ physical qubits) would anyway necessitate non-local connectivity. But this is an observation about a possible future hardware platform, not a mitigation. The paper does not propose a specific connectivity-aware optimisation, a routing cost model, or a bound on how much the routing overhead could degrade the complexity. This limitation is not addressed beyond a brief acknowledgment.

---

The Algorithm Is Not Validated Against Any Baseline Quantum Simulation Method for a Concrete System, Leaving Constant-Factor Comparison Unresolved

The assumption or constraint. The paper provides asymptotic complexity expressions (Eq. 7) and multiplication-count estimates (Appendix B), but does not compute explicit total gate counts for a specific chemical system (fixed η, N, Ω, ϵ, t) and compare them against the same quantities for direct summation, second-quantised FFT methods, or the interaction-picture approach. The break-even analysis in Appendix B counts only multiplications per particle during tree traversal, not total gate operations across all Trotter steps, arithmetic units, QROM accesses, and uncomputation. There is no "wall-clock" gate-count comparison table for, say, simulating a water molecule (η = 10) or a deuterium plasma (η = 1729) to chemical accuracy (ϵ = 10⁻³).

The consequence. The paper's central claim — "roughly a speedup by O(η) over most prior algorithms" — is a proven asymptotic statement, but the constant factors are so large that the asymptotic regime may not be reached for any accessible system size. Practitioners deciding whether to implement this algorithm for a specific simulation need to know: at η = 100, 1000, or 10,000 electrons, with realistic grid sizes, does the quantum FMM use fewer gates than direct summation? The paper does not answer this question. The break-even analysis suggests the answer is "no" for η < 10⁵, but this analysis is incomplete (only tree traversal, no Trotter step count, no comparison to second-quantised methods at specific N/η ratios).

This lack of a concrete resource estimate also makes it difficult to evaluate how the quantum FMM compares against the interaction-picture method [35] at practical N and η values. The crossover condition N < η⁷ is an asymptotic bound; at finite η and N, the constant factors could shift the crossover significantly. Without explicit gate counts, a practitioner cannot determine which algorithm to invest implementation effort in for their target system.

What evidence exists in the paper. The complexity expression (Eq. 7) and Table I provide the asymptotic comparison. The multiplication-count breakdown (Eqs. B3–B5, B14) provides constant-factor estimates for the FMM's internal operations. However, these are not combined into a total gate count per Trotter step: the number of Trotter steps (Eq. 44) is given only with asymptotic notation (the o(1) factor is left unspecified), the cost of QROM lookups for translation operators is not estimated, and the cost of the kinetic energy step (quantum Fourier transforms) is not included. There is no counterpart multiplication-count analysis for direct summation or for competing methods.

Mitigation status. This is a scoping limitation — the paper is explicitly a first demonstration of asymptotic scaling, not a resource-optimised implementation proposal. The authors state this clearly in Section VII:

"The focus of this work is on showing that it is possible to achieve the same scaling as the classical FMM."

They explicitly defer resource optimisation to future work: "reducing the overheads is a topic for future work." The Appendix B analysis is an honest attempt to estimate the constant factors, and the paper's conclusion — that break-even requires hundreds of thousands of electrons — is a candid admission that the current implementation is not practical. The limitation is inherent in the paper's scope and is acknowledged, but a practitioner seeking implementation guidance would need substantial additional analysis.

---

The Algorithm's Correctness Depends on an Even-Distribution Assumption for the Non-Adaptive Version, and the Adaptive Version's Complexity Scaling Has Not Been Rigorously Bounded for Pathological Particle Configurations

The assumption or constraint. The paper presents two versions of the quantum FMM:

1. Non-adaptive (Section IV): Assumes particles are roughly uniformly distributed so that no leaf box contains more than c = O(1) particles. The Algorithm 3 swap procedure explicitly relies on this: the proof that swapping an electron from the left child box to the right child box never overwrites existing electron data requires that "there were more than n_bc/2 electrons that should be in the right box" cannot happen — i.e., no box overflows. The paper states this assumption explicitly: "violating the even distribution assumption in this section."

2. Adaptive (Section V): Lifts the uniform-distribution assumption by allowing arbitrary particle configurations. However, the complexity analysis implicitly assumes that the adaptive tree is well-behaved — that the number of levels L scales as O(log N) and that the copying procedures in Algorithms 4 and 5 correctly handle all cases of empty boxes, multiple particles per box, and boundaries between boxes of different sizes.

The consequence. For the non-adaptive version, the algorithm is not correct on states where electrons cluster — precisely the states that occur in realistic molecules near nuclei. The simulation would need to either (a) project onto the subspace of evenly-distributed configurations (which introduces uncontrolled approximation error unless bounded), or (b) include all configurations but accept that the non-adaptive algorithm may produce incorrect results on them. Neither option is acceptable without an error bound.

For the adaptive version, the complexity analysis in Section VI provides asymptotic bounds (Eqs. 26–35) but does not rigorously analyse worst-case behaviour. The copying procedures in Algorithm 4 rely on the sorted order to identify sibling boxes; if particles are configured adversarially (e.g., all electrons in a single box at level ℓ−1, then distributed across all 8 child boxes at level ℓ), the copying operations may require stepping through many consecutive particles with the same fb(j, ℓ) — this doesn't change the asymptotic O(η) iteration count but could affect constant factors and the behaviour of the intermediate levels used for multipole aggregation. The paper does not prove a worst-case bound on the number of intermediate levels needed or the size of the data structures in the adaptive tree.

What evidence exists in the paper. The paper does not provide a rigorous error analysis for the non-adaptive version's uniform-distribution assumption violation. Section VII acknowledges this as an open problem: "If it were possible to bound the error while restricting to a subspace of evenly distributed particles, then the simpler scheme of Section IV could be used." The adaptive version's correctness is argued constructively — Algorithms 4 and 5 are presented and their logic is explained — but there is no formal proof of correctness for all input states, nor a worst-case complexity bound that accounts for adversarial particle distributions. The complexity expressions (Eqs. 26–35) are stated in O(·) notation without explicit dependence on the distribution's "clumpiness" (e.g., maximum particles per box, depth of the adaptive tree).

Mitigation status. The adaptive algorithm is designed explicitly to handle arbitrary distributions, and its mechanisms (sorting with Morton ordering, chain copying, shifted orderings) do not assume uniformity. The lack of formal worst-case bounds is a rigour gap rather than a suspected failure mode — the algorithms likely work correctly for all inputs, but the paper does not prove this or bound the constants for adversarial cases. The non-adaptive version's limitation is acknowledged and the adaptive version is presented as the solution, making this more of a presentation weakness than a fatal flaw. The adaptive algorithm is the paper's primary contribution, and its correctness argument, while informal, is detailed enough to be convincing for the asymptotic analysis — but a formal proof would strengthen the results.

7. Implications and Future Directions

How This Work Changes the Landscape

This paper does not introduce a new simulation algorithm, nor does it improve the asymptotic scaling of product formulae or FMM in their native classical setting. Its contribution is more specific and more foundational: it resolves a decades-old barrier in quantum algorithms for first-quantised electronic structure by showing that the O(η²) Coulomb bottleneck — long considered an unavoidable cost of working in real space without second quantisation — can be reduced to nearly O(η) using only standard quantum circuit primitives (sorting networks, controlled swaps, fixed-location arithmetic). This shifts the conversation from "whether" it is possible to achieve FMM-like scaling on a quantum computer to "at what overhead and break-even regime."

The magnitude of the shift is a proven upper bound, not a demonstrated practical advantage. The paper establishes a new best-known asymptotic gate complexity for first-quantised product-formula simulations when N < η⁷ — the regime the authors identify as practical. This pushes the quantum FMM into the complexity table (Table I) as the most η-efficient Trotter-based method, reducing the η² prefactor in Rubin et al. [5] to η. But the paper is equally important for what it clarifies about the limitations of this improvement: the break-even analysis in Appendix B estimates that hundreds of thousands to millions of electrons are needed for the quantum FMM to beat direct summation in total gate count, and physical implementation would require >10⁸ qubits. This candor prevents the asymptotic speedup from being misinterpreted as near-term practicality and instead frames the contribution correctly: the quantum FMM is a viable asymptotic path for very large systems on fault-tolerant hardware, not a drop-in replacement for current quantum chemistry calculations.

The paper resolves a specific contradiction in the literature. Prior to this work, there was a tension between two facts: (a) the classical FMM achieves O(η polylog(1/ϵ)) scaling and is standard in large-scale classical simulations, and (b) the only prior attempt at a quantum FMM (Childs et al. [22]) incurred an O(η) QRAM access overhead that negated the speedup, effectively making quantum FMM no better than direct summation. The implicit lesson some drew was that the FMM's pointer-based data access pattern was inherently incompatible with quantum computation — that you had to change the Hamiltonian representation (to second quantisation with FFTs, or to interaction-picture block encoding) to avoid the Coulomb bottleneck. This paper refutes that lesson by demonstrating that the incompatibility was not with the FMM's mathematical structure, but with the specific mechanism of address-based data access. By replacing pointer-following with sorting-based data rearrangement, the paper shows that the FMM's hierarchical aggregation and interaction-list evaluation are fully compatible with the quantum circuit model. The Childs et al. [22] result was not a proof of impossibility; it was a demonstration that the wrong implementation strategy fails.

This work redirects research attention toward data movement as a first-class cost in quantum algorithms. The most technically novel component — shifted Morton orderings combined with chain copying — is not about improving arithmetic or reducing the number of Trotter steps, but about organising quantum data so that spatially local information is accessible at fixed offsets in a 1D sorted list. This framing, where the primary algorithmic challenge is data layout and movement rather than arithmetic, is underdeveloped in quantum algorithms relative to classical computing (where cache hierarchies, memory bandwidth, and data locality dominate performance engineering). The paper's Lemma 1, which guarantees that 8 shifted Morton orderings bound the 1D distance to all interaction-list boxes by 63, provides a reusable building block — a geometric-to-linear locality mapping — that could apply to other quantum spatial algorithms (hierarchical matrix methods, Barnes-Hut, fast Ewald summation). Future quantum algorithms for N-body problems may begin with the question "can the data be arranged so that spatial neighbours are at fixed offsets?" rather than "can we build a QRAM large enough to dereference pointers?"

The paper also reshapes the comparison between first- and second-quantised approaches. Second quantisation with fast Fourier transforms (Babbush et al. [15], Low et al. [18]) avoids the O(η²) Coulomb sum but pays an O(N) cost that dominates when N ≫ η. First quantisation avoids the O(N) factor but previously paid the O(η²) Coulomb cost. The quantum FMM eliminates the Coulomb penalty, making first quantisation strictly superior asymptotically for N ≫ η — which is the typical regime for accurate simulations. This strengthens the case for first-quantised approaches as the default representation for large-scale quantum electronic structure, at least in the long-term fault-tolerant setting where qubit counts permit the FMM's overhead.

The work also identifies a new bottleneck. By separating the per-Trotter-step FMM cost (O(η log N log⁴(1/ϵ))) from the Trotter step count (polynomial in η^2/3N^1/3/Ω^1/3), the paper highlights that the Trotter step count, not the Coulomb evaluation, is now the dominant η-scaling factor — the total complexity scales as η^4/3 rather than η. Further asymptotic improvements therefore require tighter Trotter bounds rather than faster Coulomb evaluation. The paper's own analysis shows that the commutator-scaling bounds of Low et al. [18] already achieve near-optimal step counts; reducing the η^2/3 factor in the Coulomb potential norm would require fundamentally new product formula decompositions, not incremental improvements.

---

Follow-Up Research This Work Enables or Suggests

Reducing the break-even point from ~10⁵–10⁶ electrons to ~10³ via box-level storage. The single largest overhead in the adaptive quantum FMM is the duplication of multipole and local expansion data across all particles that share a box. In the classical FMM, each box stores one set of multipole moments and one local expansion — O(η) storage total. In the quantum version, every particle carries its containing boxes' data at every level, inflating storage and data movement costs by O(log N) per particle. The paper explicitly identifies this as the primary barrier to practical break-even (Appendix B 4) and suggests that "the implementation of the quantum FMM should be modified to only represent information for the maximum number of boxes needed." A concrete follow-up would be: design a sorting-based data layout where box data is stored once, at classically-indexed registers, and accessed by sorting particles into proximity with their containing box's data register before performing aggregation. This would require extending the sorting-based movement to handle the particle-box relationship rather than just particle-particle proximity — a non-trivial extension because boxes and particles have different numbers and the mapping from particle to box is many-to-one. A successful demonstration would reduce the per-particle data size from O(log N log³(1/ϵ)) to O(log³(1/ϵ)) and dramatically lower the break-even point. The paper's analysis in Appendix B suggests break-even could approach the classical FMM's ~500–5000 particle range if box-level storage were achieved.

Compressed multipole representations to reduce memory and gate overhead. The paper uses the full (P+1)² multipole coefficients, each with O(log(1/ϵ)) bits, per box per level. Classical FMM implementations have developed compression techniques — Tucker decomposition [54], reduced-rank far-field approximations [55], skeleton representations [56] — that reduce the number of coefficients needed to represent the far field at fixed accuracy. The paper mentions these in Section VII but does not incorporate them. A concrete follow-up would be: evaluate the quantum implementation of a skeletonised or Tucker-compressed FMM, quantifying the reduction in logical qubit count and multiplication count for realistic accuracies (P = 9, P = 18). The key question is whether compressed representations reduce the asymptotic log³(1/ϵ) data size to something smaller (e.g., log(1/ϵ) or log²(1/ϵ)) or only improve constant factors. Since the qubit count is the dominant practical barrier (~160M logical qubits for η = 4000 with the MP variant), even a factor-of-10 reduction would bring the resource estimate closer to feasibility.

Worst-case complexity analysis for adversarial particle configurations. The paper's asymptotic complexity analysis (Section VI) assumes the adaptive FMM's tree depth and box population are well-behaved — that L = O(log N) and each leaf box contains O(1) particles. For chemically realistic systems, electrons cluster near nuclei, producing a small number of very dense boxes and many empty boxes. The adaptive algorithm's correctness does not depend on uniformity, but its efficiency might: if one box contains O(η) electrons at some level, the chain-copying procedures (Algorithms 4 and 5) must step through all η particles to aggregate their contributions, potentially increasing the constant factor (though not the asymptotic O(η) scaling of the iteration). A rigorous follow-up would: prove a worst-case upper bound on the total gate count of Algorithms 4 and 6 as a function of the particle distribution's "clumpiness" (e.g., the maximum number of particles in any box at any level, or the depth of the adaptive tree). The result would either confirm that the O(η log N log⁴(1/ϵ)) bound holds for all distributions (reassuring) or reveal a pathological case where the copying overhead degrades to O(η²) (requiring a modification to the algorithm).

Numerical resource estimation for a concrete benchmark system. The paper provides asymptotic complexity expressions and Appendix B's multiplication-count analysis, but no total gate count for a specific (η, N, Ω, ϵ, t) combination. This makes it impossible to directly compare the quantum FMM against direct summation or second-quantised FFT methods at realistic system sizes (e.g., the deuterium plasma from Rubin et al. [5] with η = 1729). A concrete follow-up would: compute the exact Toffoli count (or T count) for simulating one Trotter step of the deuterium plasma benchmark (or a similar system) using the quantum FMM, direct pairwise summation, and the second-quantised FFT approach, all at equivalent target accuracy. This requires choosing a specific product formula order k, arithmetic circuit implementations for 1/∥r∥ computation and multipole translations, a sorting network with known comparator count, and a QROM cost model for translation operator lookups. The output would be a table showing, for η = 10², 10³, 10⁴, which method has the lowest gate count. This exercise would either confirm the asymptotic crossover predictions (N < η⁷) at finite sizes or reveal that the constant-factor overhead pushes the crossover far beyond the asymptotic estimate — critical information for algorithm selection in fault-tolerant quantum computing roadmaps.

Extending the shifted-Morton technique to other fast summation methods. The paper's Lemma 1 and Algorithm 6 are not specific to the Coulomb kernel or the FMM's interaction list structure. The shifted Morton technique guarantees that any target set of boxes defined by a bounded spatial offset (like the 3×3×3 neighbourhood of level-(ℓ−1) boxes that defines the FMM interaction list) is within a constant-size window under 2^d shifted orderings. This property could apply to: hierarchical matrix (H-matrix) methods for general kernel functions [62], Barnes-Hut tree codes [19], fast Ewald summation for periodic systems [59–61], and kernel-independent FMM variants [25]. A concrete follow-up would: adapt the shifted-Morton data access scheme to a quantum Barnes-Hut algorithm for gravitational N-body simulations, characterising the window size K and number of orderings required as a function of the opening angle criterion. Success would demonstrate that the sorting-based approach generalises beyond Coulomb electronic structure and that a family of classical hierarchical N-body methods can be efficiently quantised using the same data-movement primitives.

Error analysis for the non-adaptive algorithm on realistic states. The paper's simpler non-adaptive quantum FMM (Section IV) has dramatically lower overhead than the adaptive version because it uses fixed box registries (no per-particle duplication) and does not require the shifted-Morton interaction-list access. However, it assumes uniform particle distribution, which is violated on the molecular eigenstates that the Hamiltonian simulation acts on. A formal follow-up would: bound the error introduced by projecting the Hamiltonian onto the subspace of states where no leaf box contains more than c electrons (for some fixed c), as a function of c, N, η, and the nuclear geometry. If the projection error can be bounded by something like exp(−c) (due to the Coulomb repulsion suppressing configurations with many electrons in the same small volume), then a modest c (e.g., 4–8) might suffice for chemical accuracy, making the non-adaptive algorithm viable and drastically reducing the break-even point. The paper flags this explicitly: "If it were possible to bound the error while restricting to a subspace of evenly distributed particles, then the simpler scheme of Section IV could be used." This is a well-defined mathematical question about the electronic structure Hamiltonian's eigenstate support, not an algorithm-design challenge, and answering it would directly determine which FMM variant is preferable.

Hilbert ordering vs. shifted Morton ordering: a constant-factor comparison. The paper uses Morton ordering because bit-interleaving is free in a quantum circuit (qubit relabeling), but Hilbert curves [40, 41] have strictly better locality — the maximum 1D distance between neighbouring boxes in 3D is smaller for Hilbert than Morton, potentially reducing the window size K and the number of shifts needed. The cost is that computing a Hilbert index requires coherent arithmetic. A concrete follow-up would be: estimate the total gate count (sorting + arithmetic) for the quantum FMM using a single Hilbert ordering (no shifts) vs. 8 shifted Morton orderings, as a function of η, N, and ϵ. If the Hilbert arithmetic overhead (O(log N) additions and bit operations per particle to compute the Hilbert index) is smaller than the 8× redundancy of shifted Morton sorting (8 sorts of size O(η log η) each), Hilbert could be the preferred choice despite the more complex mapping. The paper acknowledges this trade-off (Section VII) but does not quantify it; a head-to-head resource comparison at realistic η would resolve whether Morton's simplicity or Hilbert's locality wins in the fault-tolerant regime.

---

Practical Applications and Downstream Use Cases

Fault-tolerant quantum simulation of warm dense matter and plasmas at scale. The most direct application is to simulations that are already pushing classical methods to their limits and require thousands of electrons: warm dense matter (conditions relevant to inertial confinement fusion and planetary interiors), degenerate plasmas, and electron dynamics in materials under extreme conditions. Rubin et al. [5] simulated 1729 electrons for a deuterium plasma using a first-quantised product-formula approach with O(η²) Coulomb scaling. The quantum FMM would reduce the Coulomb portion of that simulation's gate count by a factor of η ≈ 1729 (in the asymptotic limit), potentially bringing the total gate count for similar-sized systems within reach of future fault-tolerant hardware. However, the paper's own break-even analysis (Appendix B) estimates that hundreds of thousands of electrons are needed for the FMM to beat direct summation in practice given current overheads, so this application would require either (a) systems significantly larger than those studied today, or (b) successful implementation of the optimisations discussed above (box-level storage, compressed multipoles, improved sorting networks) to bring break-even down to the 10³–10⁴ range. The paper's qubit estimate (~160M logical qubits for η = 4000 with the MP variant) sets a concrete target for hardware roadmaps: a fault-tolerant quantum computer capable of simulating plasma electrons at scale would need to operate in this qubit regime, comparable to ambitious targets for materials simulation [52].

Quantum phase estimation for molecular energies in first quantisation. The electronic structure Hamiltonian simulated here (Eqs. 1–4) is the standard starting point for quantum phase estimation (QPE) algorithms that compute ground-state energies of molecules. In QPE, the evolution time t scales as O(1/ϵ) (where ϵ is the target energy precision), and the Trotter step count replaces t with 1/ϵ throughout the complexity expression. The quantum FMM directly reduces the per-step cost of the Coulomb evaluation in QPE, with the total gate complexity becoming (1/ϵ)(η^4/3N^1/3 + η^1/3N^2/3)(ηN/ϵ)^o(1). For molecules where high grid resolution is needed (large N, as in all-electron calculations with sharp nuclear cusps) but the number of electrons is moderate (η ~ 10–100), the η^1/3N^2/3 term may dominate — and the quantum FMM reduces the η-dependence of this term from η (in Rubin et al. [5]) to η^1/3, a substantial improvement. This application would benefit most from a refined break-even analysis at small η with large N, which the paper does not provide — the Appendix B estimates are for large η — but the asymptotic complexity improvement applies at all scales.

Self-improving quantum simulation pipelines via iterative refinement. The paper's FMM computes the Coulomb potential coherently and uses it for phase kickback, with the potential energy never measured. However, the same FMM subroutine could be used to prepare improved initial states for QPE — for example, by performing classical optimisation of trial wavefunctions where the energy (including the Coulomb term) is evaluated using the quantum FMM, or by using the FMM as a subroutine in quantum-classical hybrid algorithms like VQE where the expected value of V is measured rather than used for phase rotation. The linear-in-η scaling means that energy evaluation (which requires ~O(1/ϵ²) measurements in VQE) becomes dramatically cheaper for large systems, potentially making first-quantised VQE competitive with second-quantised approaches for systems with N ≫ η. The paper does not explore this application, but the FMM's availability as an O(η) coherent subroutine for Coulomb evaluation makes it a drop-in replacement wherever V appears in first-quantised algorithms.

Kernel-independent fast summation on quantum computers for non-Coulombic potentials. The paper's methods are not tied to the specific 1/r Coulomb kernel — the FMM framework generalises to any kernel Φ(∥r_i − r_j∥) that satisfies a multipole expansion (i.e., is a solution to the Laplace equation or can be so approximated). The paper notes in Section VII that adapting the quantum FMM to periodic systems with the Ewald kernel (Eq. 47) would require "straightforward modification of the interaction lists and certain translation operators, as in classical implementations [59–61]." More broadly, the shifted-Morton data access scheme (Algorithm 6) does not depend on the kernel at all — it only depends on the tree geometry and the definition of the interaction list. This means the sorting-based quantum FMM could serve as a generic quantum subroutine for kernel summation, applicable to simulations with Yukawa potentials, dipolar interactions, or general 1/r^ν kernels. The practical barrier remains the break-even point, but the framework is reusable. The paper specifically names hierarchical matrix methods [62] as a related target, noting that they "replace the interaction list with admissible blocks" — extending the shifted-Morton technique to H-matrices would require handling overlapping blocks, which the paper identifies as future work.