1 Introduction↩︎

Monte Carlo (MC) integration is central to computational physics, statistical mechanics, and machine learning [1], [2]. Drawing \(M\) i.i.d.samples and forming the empirical mean gives an unbiased estimator with RMSE \(\sim M^{-1/2}\), the asymptotically optimal rate for sample-mean estimators by the CLT and the Cramér–Rao bound. Reaching precision \(\varepsilon\) requires \(M\sim\varepsilon^{-2}\) samples. Quasi-Monte Carlo (QMC) methods [3], [4] improve the prefactor via low-discrepancy sequences but retain the sample-mean structure and its asymptotic variance dependence.

1.0.0.1 Quantum acceleration of MC integration.

Two quantum frameworks have been studied. Quantum amplitude estimation (QAE) [5], [6] estimates an expectation value with \(O(\varepsilon^{-1})\) oracle calls (a quadratic improvement) and has been applied to option pricing [7], [8] and QMC integration [9], though the oracle construction may dominate resource cost on near-term hardware. Szegedy walks [10] accelerate Markov-chain mixing by \(\sqrt{\Delta_c}\) in the spectral gap, underpinning quantum partition-function [11] and simulated-annealing [12] algorithms. Both frameworks accelerate sample-mean estimators and rely on phase estimation or oracle primitives demanding near-fault-tolerant resources.

1.0.0.2 This work.

Coined discrete-time quantum walks (DTQWs) [13], [14] spread ballistically. The Konno limit theorem [15], [16] gives \(\mathrm{Var}(X_T)=\Theta(T^2)\) against the diffusive \(\Theta(T)\), with 2D analogues established by Watabe et al. [17]. We ask whether this faster spreading can be exploited for geometric integration, and identify a class of estimators in which it can. The construction joins (i) the ballistic 2D DTQW position distribution, concentrating probability mass at scale \(\Theta(T)\), and (ii) the Hardy–Huxley lattice asymptotic \(N(R)/R^2 = \pi + O(R^{\gamma-2})\) with \(\gamma\le131/208\) [18], [19]. Taking the \(\alpha\)-quantile \(\hat{R}(\alpha)\) of \(M\) radial measurements and forming \(N(\hat{R})/\hat{R}^2\) yields an estimator of \(\pi\) whose RMSE is dominated by a bias floor \(O(T^{\gamma-2})\) falling with \(T\), not \(M\). For precision \(\varepsilon\) strictly above this floor, a fixed \(\varepsilon\)-independent sample count suffices. This construction applies to integrals reducible to convex-region lattice counts; discrete partition functions, non-convex Boltzmann integrals, and oracle-encoded integrands lie outside its scope (Sec. II G). DTQWs have been realised on IBM superconducting devices, with circuit-depth limits currently restricting reliable runs to short walks [20], [21], and a recent work [22] couples DTQW with QAE for radiation transport, targeting a different problem class than the geometric integrals addressed here. The closest prior quantum estimate of \(\pi\) using the in-circle lattice criterion is that of Gustafson et al. [23], who estimate \(\pi\) via quantum mean estimation, where quantum phase estimation accelerates a sample-mean estimator over a uniform grid. Our construction differs on three structural axes. It is oracle-free in the QAE sense (Sec. 4) and uses no quantum phase estimation, the estimator is a rank statistic rather than a sample mean, and the accuracy-determining quantity is the deterministic lattice-count residual controlled by walk depth rather than the statistical variance controlled by query count.

Reaching precision \(\varepsilon\) above the bias floor requires \(T=\Theta(\varepsilon^{-1/(2-\gamma)})\) and a fixed, \(\varepsilon\)-independent \(M\), giving a sample-count ratio \(M_{\mathrm{MC}}/M_{\mathrm{QW}}=\Theta(\varepsilon^{-2})\) relative to uniform MC. We refer to this \(\varepsilon^{-2}\) ratio as the ballistic acceleration; the corresponding ratio relative to scrambled Sobol QMC is \(\sim\!500\times\) by order-of-magnitude extrapolation (Sec. III A). This is structurally distinct from the \(\varepsilon^{-1}\) oracle-call scaling of QAE and the \(\sqrt{\Delta_c}\) mixing-time scaling of Szegedy walks.

1.0.0.3 Contributions.

Our contributions are six. (i) A sample-complexity bound (Proposition 1) for the disk estimator with explicit error decomposition into bias-floor and sampling terms. (ii) An envelope-form bias-floor identity (Theorem 2) that supplies a parameter-free bias floor at any \((T,\alpha)\) and reproduces the measured floors to within \(1.5\times\). (iii) An asymptotic variance optimum (Proposition 3) yielding a closed-form \(\alpha^\star\) as a function of the Konno limit measure. (iv) A generalised estimator (Proposition 4) extending the framework to convex smooth regions via Krätzel’s lattice asymptotic. (v) A Cavalieri-based estimator (Proposition 5) for smooth-function integrals, together with a corollary for annular super-level sets (Corollary 1). (vi) Numerical experiments on \(\pi\) (\(\sim 7\times 10^4\times\) sample-count reduction at \(T=200\), \(M_{\mathrm{QW}}=100\) relative to uniform MC; the floor is reached already at \(M=50\); \(\sim 500\times\) relative to scrambled Sobol QMC by order-of-magnitude extrapolation), area estimation across eight 2D shapes (reductions of \(10^2\) to \(10^7\times\), Appendix 7), single-batch error propagation across a 2D Gaussian cumulative distribution, and recovery of the 2D quantum-harmonic-oscillator energy ladder \(E_n=\hbar\omega(n+1)\) (\(n=0,\dots,3\)) with a level spacing \(10\times\) tighter than the diffusive random-walk baseline. Every experiment is benchmarked against the classical random walk paired with the identical estimator. The phase-averaged envelope slopes lie in \([-1.506,-1.489]\), with a mean of \(-1.499\approx-3/2\), confirming the Hardy conjecture.

2 Methods↩︎

The theoretical development proceeds in five steps. We first recall why every sample-mean estimator (classical Monte Carlo and quasi-Monte Carlo alike) is bound to the \(M^{-1/2}\) rate by the Cramér–Rao limit (Sec. 2.1). We then introduce the ballistic 2D quantum walk whose position distribution concentrates at the lattice scale \(\Theta(T)\) (Sec. 2.2), and combine it with the Gauss circle lattice asymptotic to build the disk estimator for \(\pi\), whose error splits into a depth-controlled deterministic floor and a sample-controlled statistical term (Sec. 2.3). Two extensions follow from the same identity. The first is a generalised estimator for convex regions via Krätzel’s bound (Sec. 2.4); the second is a Cavalieri estimator for smooth integrals that fixes an entire integral family from one walk through classically precomputed multipliers (Sec. 2.5). Finally we explain why the construction sidesteps the sample-mean barrier. The accuracy-determining quantity is relocated from a statistical term to a deterministic number-theoretic one, and we then delimit precisely which integrals fall in scope (Secs. 2.6–2.7).

2.1 Classical baselines and the Cramér–Rao limit↩︎

We develop a family of lattice-quantile estimators that share the construction \(N_D(\hat{R})/\hat{R}^2\) for some lattice count \(N_D\) and rank statistic \(\hat{R}\). The family has three members: the disk estimator for \(\pi\) (Sec. 2.3), the generalised disk estimator for convex smooth domains (Sec. 2.4), and the Cavalieri estimator for smooth integrals (Sec. 2.5), with an annular-Cavalieri corollary for non-monotonic radial profiles (Cor. 1). All members inherit the same bias-floor structure and Theorem 2 envelope; their differences reside in the choice of lattice count and classical prefactor.

For the canonical task of estimating \(\pi\) via the unit disk \(D = \{(x,y) : x^2+y^2 \le 1\}\), dart-throw Monte Carlo forms the sample-mean indicator \[\hat{\pi}_{\mathrm{MC}} = \frac{4}{M}\sum_{i=1}^M \mathbf{1}_D(X_i, Y_i), \qquad (X_i,Y_i)\sim\mathrm{Unif}([-1,1]^2), \label{eq:mc}\tag{1}\] with \(\mathrm{RMSE}(\hat{\pi}_{\mathrm{MC}})=\sqrt{\pi(4-\pi)/M}\approx 1.64/\sqrt M\) [2]. The \(M^{-1/2}\) rate is the optimal scaling for any unbiased sample-mean estimator of an expectation value, by the CLT together with the Cramér–Rao bound [2], [24]: \[\mathrm{RMSE}(\hat{\theta}_M) = \Omega\bigl(M^{-1/2}\bigr) \label{eq:cramer-rao}\tag{2}\] for unbiased sample-mean estimators of an expectation parameter \(\theta=\mathbb{E}[f(X)]\). Quasi-Monte Carlo (QMC) replaces the i.i.d.samples by deterministic low-discrepancy sequences and obtains Koksma–Hlawka rates \(O(M^{-1}\log^2 M)\) in \(d=2\) [3], [4], but retains the sample-mean structure and is therefore still bound by Eq. 2 as far as the expectation-value estimator is concerned. The lattice-quantile estimator we develop below sidesteps this constraint not by improving the sampling distribution but by replacing the estimator structure with a deterministic lattice-counting function of a sample quantile.

2.2 Two-dimensional coined DTQW↩︎

Let \(\mathcal{H} = \mathbb{C}^4 \otimes \ell^2(\mathbb{Z}^2)\) be the Hilbert space of a 2D coined DTQW with four-state coin register and integer position register. The walk operator is \(U = S\,(C \otimes \mathbb{I})\) with coin operator \(C \in U(4)\) and conditional shift \[S\,|c\rangle\,|x,y\rangle = |c\rangle\,|x + d_c, y + e_c\rangle, \label{eq:shift}\tag{3}\] where \((d_c, e_c) = ((-1)^{c_1}, (-1)^{c_0})\) for the binary coin label \(c = (c_1 c_0)_2 \in \{0,1,2,3\}\). We use the 2D Grover diffusion operator \[G = 2|s\rangle\langle s| - \mathbb{I}_4,\qquad |s\rangle = \tfrac{1}{2}\textstyle\sum_{c=0}^{3}|c\rangle, \label{eq:grover}\tag{4}\] the unique coin (up to local phases) respecting the full \(D_4\) symmetry of the square lattice and the standard choice in 2D quantum walks [17]. Two further coins are used to test how stable the bias-floor scaling is under a change of coin: the separable Hadamard coin \(H\otimes H\) and the 4-point Fourier coin \(F\). Their explicit matrices and limit-measure properties are collected in Appendix 6.

The initial state \[|\psi_0\rangle = \tfrac{1}{2}(|0\rangle + i|1\rangle + i|2\rangle - |3\rangle) \otimes |\mathbf{0}\rangle \label{eq:initial}\tag{5}\] is chosen so that the position distribution \(p_T(x,y) = \sum_c |\langle c, (x,y)|U^T|\psi_0\rangle|^2\) carries the full dihedral \(D_4\) symmetry of the square lattice (invariance under the eight reflections and rotations mapping \(\mathbb{Z}^2\) to itself). This \(D_4\) symmetry is what the disk estimator requires: it guarantees that the radial marginal \(r/T\) is well-defined and that the \(\alpha\)-quantile \(\hat{R}(\alpha)\) is isotropic, regardless of whether the underlying support is the velocity disk (Grover, Fourier) or the velocity square (Hadamard). Full continuous rotational invariance holds only in the Grover and Fourier cases; the weaker \(D_4\) symmetry, which holds for all three coins, suffices for everything that follows.

The \(T^2\) scaling of \(\mathrm{Var}(X_T)\) implied by the lemma is the ballistic property of quantum walks, contrasting with the diffusive \(\Theta(T)\) scaling of classical random walks. The disk estimator we develop next uses only the radial marginal \(r/T = \sqrt{v_x^2 + v_y^2}\), for which the relevant question is the location of the radial \(\alpha\)-quantile rather than the precise shape of the support.

2.3 Disk estimator for \(\pi\)↩︎

Given \(M\) projective position measurements \(\{(X_i, Y_i)\}_{i=1}^M\) from the DTQW, define the radial samples \(r_i = \sqrt{X_i^2 + Y_i^2}\) and the \(\alpha\)-quantile \(\hat{R}(\alpha) = r_{(\lceil \alpha M\rceil)}\) (where \(r_{(k)}\) is the \(k\)-th order statistic). The disk estimator is \[\boxed{\;\hat{\pi}(\alpha; T, M) = \frac{N(\hat{R}(\alpha))}{\hat{R}(\alpha)^2}\;} \label{eq:disk}\tag{6}\] with \(N(R) = \#\{(i,j)\in\mathbb{Z}^2 : i^2 + j^2 \le R^2\}\) computed classically.

The estimator combines two ingredients:

Remark (status of the bound). This result is stated as a Proposition rather than a Theorem for two reasons. First, the exponent \(\kappa(\alpha)\) is not independently identifiable from our data: across all configurations tested, the sampling term is already sub-dominant to the bias floor at \(M=5\) (Sec. [sec:results-scaling]), so only its smallness, not its scaling, is empirically constrained. For the Grover coin we give a closed-form derivation in Cor. 2 (Appendix 8) under Conjecture 6, yielding \(\kappa(\alpha)=1\) with an \(\alpha\)-suppressed prefactor; the Hadamard and Fourier cases remain phenomenological. See Appendix 6 for the heuristic Bahadur–Kiefer argument and its caveats for the Grover localised measure, and Lemma 4 for the truncated-BK statement under which the Grover argument applies. Second, the leading \(C_1 T^{\gamma-2}\) term is an envelope, not a pointwise description; the precise behaviour of the measured floor is captured by Theorem 2 below, and explicit upper bounds on \(C_1, C_2\) for Grover are given in Prop. 7.

Remark (empirical validation). Equation ?? is the strongest parameter-free prediction of the framework. Table ¿tbl:tab:floor? compares it to measured floors at \(T\in\{50,100,200,400\}\) and finds agreement to within \(1.5\times\) (mean over 32 seeds at each depth). It also explains two features a monotone power law cannot: the spread of fitted exponents across coins, and the non-monotonic small-\(M\) behaviour reported in Sec. [sec:results-scaling].

Remark (practical relevance of \(\alpha^\star\)). Equation ?? minimises the sampling term of Eq. ?? ; we downgrade it from a Theorem to a Proposition because that term is sub-dominant to the bias floor for all \(M\ge 5\) tested (Sec. [sec:results-scaling]), so the variance optimum is not what governs the measured RMSE in practice. The empirically relevant \(\alpha\)-dependence enters instead through Theorem 2: changing \(\alpha\) moves \(R^\star(\alpha)\) along the oscillating number-theoretic residual, so the \(\alpha\)-sweep of Appendix A (Fig. 5 (c)) should be read as a scan of the floor landscape rather than as variance minimisation. The agreement of its minima with Eq. ?? is therefore partly coincidental; a principled choice of \(\alpha\) would target zeros of the residual (quantile tuning, Sec. [sec:results-scaling]).

2.4 Generalised disk estimator for convex regions↩︎

The framework extends to arbitrary convex smooth domains. Let \(D \subset \mathbb{R}^2\) be a bounded convex region with \(C^2\) boundary. Define \(N_D(L) = \#\{(i,j) \in \mathbb{Z}^2 : (i/L, j/L) \in D\}\).

The curvature condition is genuine: it fails on the axes of the super-ellipse \(x^4+y^4\le1\), where the curvature vanishes. Empirically (Appendix 7) the estimator still performs well there, but with a slightly larger floor than the strictly-curved cases, consistent with the boundary degeneracy contributing a sub-leading lattice-counting correction.

For an ellipse \(D=\{(x,y):x^2/a^2+y^2/b^2\le 1\}\) with \(|D|=\pi a b\), the lattice count is \(N_D(L)=\#\{(i,j):(ib)^2+(ja)^2\le(abL)^2\}\). Boundary regimes that violate the Krätzel hypothesis (vanishing curvature, cusp, polygonal corners, doubly-connected annulus) are treated empirically in Appendix 7.

2.5 Cavalieri estimator for smooth integrals↩︎

For a non-negative smooth integrand \(f: \Omega \to \mathbb{R}_{\ge 0}\) on a bounded region \(\Omega \subset \mathbb{R}^2\), Cavalieri’s principle gives \[I = \iint_\Omega f(x,y)\, dx\, dy = \int_0^{f_{\max}} |S_t|\, dt, \label{eq:cavalieri}\tag{7}\] where \(S_t = \{(x,y) \in \Omega : f(x,y) \ge t\}\) is the super-level set at threshold \(t\). When \(f\) is quasi-concave with rotationally symmetric super-level sets (i.e., \(S_t\) is a disk of radius \(r(t)\)), 7 reduces to \(I = \pi \int_0^{f_{\max}} r^2(t)\, dt = \pi\, Q\), where \(Q\) is computable analytically or by classical numerical quadrature.

The estimator requires only a single batch of \(M\) projective measurements (one batch = \(M\) shots of the state-preparation–evolution–measurement circuit, post-processed into one \(\hat{\pi}\)); the multiplicative constant \(Q\) is precomputed classically. For ellipse-shaped super-level sets, an analogous expression \(\hat{I} = \widehat{\pi a b}(\alpha; T, M) \cdot Q'\) uses the generalised disk estimator ?? . Statistical replicates appearing later in the paper (e.g. trials at \(M=500\) in Figs. 3 and 4) are used only to estimate RMSE; the deployment unit is a single batch.

The disk-shaped hypothesis of Proposition 5 excludes one physically important case that we use in Sec. 3.3: integrands whose super-level sets are annular rather than disk-shaped. This arises whenever the radial profile \(f(r)\) is non-monotonic with an interior maximum, as for the excited-state densities \(r^{2n}e^{-r^2/\sigma^2}\) (\(n\ge1\)) of the harmonic oscillator. We treat this case as a formal corollary.

The proof makes no use of \(C^2\) convexity of the annular boundary because the estimator never counts lattice points inside the annulus: it only uses the disk-estimator output \(\hat{\pi}\) and the precomputed classical \(Q\). The doubly-connected geometry is absorbed into the deterministic constant. We verify the prediction empirically in the annulus entry of Appendix 7 and the energy ladder of Sec. 3.3.

2.6 Origin of the sample-count reduction↩︎

The lattice-quantile estimator’s behaviour relative to Monte Carlo arises from two structural features. First, \(\hat{R}(\alpha)\) is a rank statistic, not a sample mean: by the Bahadur–Kiefer representation [24], [25], its relative fluctuation is the classical \(\Theta(M^{-1/2})\), but its absolute fluctuation carries the ballistic prefactor \(T\) from \(R^\star=\Theta(T)\) (proof of Proposition 1). Second, \(N(R)/R^2\) applied to this rank statistic is a deterministic number-theoretic count whose convergence to \(\pi\) is controlled by the Gauss circle exponent (Lemma 2), a lattice property containing no statistical randomness. The composite \(\hat{\pi}=N(\hat{R})/\hat{R}^2\) is therefore a deterministic function of a rank statistic. The Cramér–Rao bound on sample means 2 does not apply directly: \(\pi\) is not a parameter of the sampling distribution, the dominant error at large \(M\) is the deterministic residual \(O(T^{\gamma-2})\) (Lemma 2), and the bound still governs the residual sampling fluctuation of \(\hat{R}\) at the standard \(M^{-1/2}\) rate. The framework does not evade the \(M^{-1/2}\) statistical law; rather, it relocates the accuracy-determining quantity from a statistical term to a deterministic number-theoretic one controlled by walk depth.

2.6.0.1 The role of ballistic spreading.

The bias-floor exponent \(T^{\gamma-2}\) requires that the walker’s typical position scales as \(\Theta(T)\). If position scales only diffusively as \(\Theta(\sqrt T)\), the same disk estimator paired with a classical 2D random walk produces a bias floor at the shallower rate \(T^{(\gamma-2)/2}\), giving \[\frac{\mathrm{RMSE}^{\mathrm{(QW)}}_{\mathrm{floor}}}{\mathrm{RMSE}^{\mathrm{(RW)}}_{\mathrm{floor}}} = \Theta\bigl(T^{-(2-\gamma)/2}\bigr) \to 0 \quad (T \to \infty). \label{eq:qw-rw-ratio}\tag{8}\] At fixed precision \(\varepsilon\) above the QW bias floor, MC needs \(M=\Theta(\varepsilon^{-2})\) samples while the lattice-quantile estimator needs only \(M=O(1)\), so the sample-count ratio behaves as \(M_{\mathrm{MC}}/M_{\mathrm{QW}}=\Theta(\varepsilon^{-2})\). We refer to this \(\varepsilon^{-2}\) ratio as the framework’s ballistic acceleration; it should be understood as a sample-count comparison at fixed precision rather than a complexity-class separation in the technical sense.

2.7 Scope of applicability↩︎

The framework applies precisely to integrals reducible to lattice point counts in convex smooth regions. Three concrete classes lie within scope. (i) Direct area or volume of convex smooth domains. (ii) Smooth integrals via Cavalieri whose super-level sets are convex. (iii) Quadratic Boltzmann integrals whose super-level sets are ellipsoidal. Problems outside this class (discrete partition functions such as Ising and spin glasses, non-convex Boltzmann integrals, one-dimensional integration) fall outside the framework and remain the domain of complementary methods. The most relevant of these are Szegedy walks for MCMC-based partition function estimation [10]–[12] and amplitude estimation for general oracle-encoded integrands [5], [6].

3 Results↩︎

We organise the numerical results as a progression of increasing physical content, each stage benchmarked against the classical random walk (RW) paired with the identical lattice-quantile estimator so as to isolate the contribution of ballistic spreading. Section 3 A benchmarks \(\pi\) estimation against uniform Monte Carlo, scrambled Sobol QMC, and the RW baseline, and extends to eight 2D shape classes. Section 3 B quantifies how the accuracy of a single \(\hat{\pi}\) propagates across an entire Gaussian cumulative distribution through a classically known multiplier. Section 3 C applies the framework to a spectroscopic task—the energy ladder of the 2D quantum harmonic oscillator—and demonstrates a controlled, analytically predictable error budget across levels. All quantum results use the Grover coin unless otherwise stated, computed on a noiseless statevector emulator of the walk dynamics.

3.1 Benchmarking \(\pi\) estimation across methods↩︎

Figure 1 plots the precision \(\varepsilon\) (RMSE of \(\hat{\pi}\)) achievable at a given sample count \(M\) and fixed walk depth \(T=200\), the central comparison of this work. The two classical sample-mean methods improve only by spending samples: dart-throw MC rides the CLT line \(\varepsilon \propto M^{-1/2}\), from \(\varepsilon=0.18\) at \(M=100\) down to \(0.017\) at \(M=10^4\), and scrambled Sobol QMC sits below it with the faster empirical slope \(\varepsilon\sim M^{-0.76}\). Both lattice-quantile walkers behave qualitatively differently: their precision is flat in \(M\), set by the bias floor rather than the sample count. Both QW estimators (Grover and Fourier) sit at \(\varepsilon\sim 6\times 10^{-4}\) from \(M=50\) onward, while the classical random walk paired with the identical estimator (grey) is flat at \(\varepsilon\approx 4\)–\(7\times10^{-3}\). The practical reading is direct: to reach the QW precision at \(M=100\), MC would need \(M\approx 7\times 10^6\) samples—a sample-count ratio \(M_{\mathrm{MC}}/M_{\mathrm{QW}}\approx 7\times 10^4\) (from \(M_{\mathrm{MC}}^\star = 16\,\mathrm{Var}(\mathbf{1}_D)/\mathrm{RMSE}_{\mathrm{QW}}^2\) with \(\mathrm{Var}(\mathbf{1}_D) = \pi(4-\pi)/16\))—while the diffusive walker cannot reach it at any \(M\). The horizontal QW lines lying an order of magnitude below the grey RW line, both flat, is the visual signature that the advantage comes from ballistic spreading and not from the rank-statistic construction, which both share.

3.1.0.1 Comparison with quasi-Monte Carlo.

Because uniform dart-throw MC is the weakest classical baseline, we also benchmark against scrambled Sobol quasi-Monte Carlo (QMC) [3], which exploits low-discrepancy structure and is widely used in practice. Over 500 independent scrambles, Sobol achieves \(\mathrm{RMSE} \approx 7.8\times 10^{-2}\) at \(M=100\) and \(\approx 2.6\times 10^{-3}\) at \(M=8192\), with an empirical slope of approximately \(M^{-0.76}\) in the range \(M\in[64,8192]\)—better than uniform MC’s \(M^{-1/2}\) but decelerating relative to the theoretical \(O(M^{-1}\log M)\) due to the discontinuous indicator boundary. Extrapolating this fit to the QW bias floor \(\mathrm{RMSE}\approx 6\times 10^{-4}\) gives \(M_{\mathrm{Sobol}}\approx 5\times 10^4\), a factor of \(\sim\!500\times\) more measurements than the QW estimator’s \(M=100\). We note that this extrapolation spans roughly three decades from the QW reference point (\(M=100\)) to the target (\(M\approx5\times10^4\)) and should be treated as an order-of-magnitude estimate. Moreover, since the empirical slope is already decelerating (from \(M^{-0.76}\) to slower) as \(M\) grows—reflecting the degradation of Sobol convergence for discontinuous indicator boundaries at large \(M\)—the actual \(M_{\mathrm{Sobol}}\) needed may be larger than \(5\times 10^4\), making this estimate conservative in favour of Sobol (and thus the QW advantage may be understated). Table ¿tbl:tab:baseline? summarises the four-way comparison.

3.1.0.2 Comparison with the classical random walk.

MC and Sobol differ from the QW estimator in both the sampling distribution and the estimator structure. The cleanest control that isolates the quantum contribution is the classical 2D random walk paired with the same lattice-quantile estimator: identical rank-statistic construction, with ballistic spreading replaced by diffusive. At \(T=200\) this RW baseline is also bias-dominated—its RMSE stays flat at \(\approx 4\)–\(7\times10^{-3}\) for all \(M\in[50,10^4]\) (grey triangles in Fig. 1)—but its floor sits \(\sim 7\)–\(10\times\) above the QW floor, because the diffusive walker reaches only radius \(\Theta(\sqrt T)\) and the lattice-counting error enters at the shallower rate \(T^{(\gamma-2)/2}\) (Eq. 8 ). Empirically the RW floor scales as \(T^{-0.85}\) over \(T\in\{50,100,200,400\}\), so reaching \(\mathrm{RMSE}\approx6\times10^{-4}\) would require \(T\approx 2\times10^3\), an order of magnitude more walk steps than the QW’s \(T=200\). The gap between the grey and blue curves in Fig. 1 is therefore attributable to ballistic spreading alone, not to the estimator construction.

Figure 2 reports the empirical bias-floor scaling for the three coins. The floor is the RMS of an oscillating number-theoretic residual over the quantile neighbourhood (Eq. ?? ), so its value at any single radius depends on the local phase of that oscillation; a naive power-law fit to raw single-radius floors is correspondingly unstable, returning coin-dependent slopes from \(-1.5\) to \(-1.8\) that do not reflect a true asymptotic exponent (panels (b)–(c) make this failure explicit). Evaluating the floor in its phase-averaged form instead yields a clean, coin-independent result, as we now describe.

Table ¿tbl:tab:floor? validates the floor formula directly: the parameter-free RMS prediction of Eq. ?? reproduces the measured floors to within \(1.5\times\) at every depth, whereas the single-point residual \(|N(R^\star)/R^{\star2}-\pi|\) misses by up to \(5\times\) (it can land arbitrarily close to a zero crossing of the oscillation, as happens at \(T=400\)). The oscillatory structure also produces a behaviour that a monotone bias-plus-variance picture forbids and that we observe directly: at \(T=100\), the RMSE at \(M=5\) is \(0.60\times\) the \(M=2\times10^4\) floor (400 trials), because the finite-\(M\) quantile bias shifts \(\hat{R}\) to a more favourable phase of the residual.

Figure 2 (a) reports the bias floor actually delivered by the estimator—the value at the single quantile radius \(\hat{R}\) that one measured walk provides—out to \(T=800\). The measured floors do not trace a clean power law: they oscillate, with each coin’s floor dipping and rising as the depth changes (e.g.the Grover floor falls to \(9.6\times10^{-5}\) at \(T=500\), rises to \(1.56\times10^{-4}\) at \(T=600\), then falls again at \(T=800\)). This is not measurement noise—each point is averaged over \(30\) trials at \(M=2\times10^4\)—but the deterministic consequence of reading a sign-oscillating number-theoretic residual at a single radius, the mechanism dissected in panels (b)–(c). Overlaid as a faint line is the phase-averaged floor of Eq. ?? : the RMS of the same residual over one lattice period around the measured \(\hat{R}\), which is coin-independent and decays at slope \(-1.50\), coincident with the Hardy rate \(T^{-3/2}\) to three significant figures. The measured single-radius floors scatter about this envelope. We emphasise the distinction because the envelope is an analytic trend, not a quantity any single run delivers; the markers are what the estimator returns. The diffusive RW floor, by contrast, sits \(5\)–\(30\times\) higher and decays at the much shallower \(T^{-0.85}\), never approaching the quantum band.

3.1.0.3 Why the data may lie above or below the reference lines.

A note on reading the reference lines. The Hardy and Huxley lines mark slopes, not absolute heights: Huxley’s theorem bounds the single-point discrepancy \(|N(R)-\pi R^2|\le C\,R^{131/208}\) with an unspecified constant, and Hardy’s is a conjectured lower envelope of the amplitude, so the vertical placement of either line is an arbitrary normalisation. The meaningful comparison is therefore the decay exponent, and it is the phase-averaged floor—not the raw single-radius readout—whose slope the theorem speaks to. We make this distinction concrete in the next two panels, because the raw readout is badly misleading.

3.1.0.4 Why the floor must be phase-averaged.

Why the phase average is essential, and not a cosmetic choice, is shown in Fig. 2 (b)–(c). The raw floor at each depth reads the oscillating residual \(E(R)/R^2=N(R)/R^2-\pi\) off at the single radius \(R^\star(\alpha)\) that the \(\alpha\)-quantile selects [panel (b)]. As \(T\) increases, \(R^\star\) lands on a different phase of this oscillation—near a zero crossing at \(T=400\) (residual anomalously small), near a trough at \(T=600\) (residual large)—so the single-radius value does not track the underlying trend. The consequence [panel (c)]: connecting the raw single-radius values gives a slope of \(+3.9\) over \(T=400\)–\(600\)—the wrong sign, a floor that appears to rise with depth—which no genuine decay law can produce. Indeed, fitting raw single-radius floors over \(T\in\{400,500,600,800\}\) in isolation yields wildly coin-dependent slopes (\(-2.1\) Grover, \(-2.7\) Fourier, \(-4.4\) Hadamard); since all three count the same residual \(E(R)\), a true exponent must be coin-independent, so these values cannot be asymptotic. A single-radius readout of a sign-oscillating function distorts, and can even invert, the apparent slope—hence the phase average.

The phase average [purple in panel (c)] removes the distortion: it is monotone, with \(T=400\)–\(600\) slope \(-1.42\) and adjacent slopes \(-1.53,-1.61,-1.39\), no sign reversals, inside the Hardy–Huxley window. We therefore make no claim of an asymptotic exponent from raw single-radius fits—their slope ranges from \(-0.9\) to \(-4.4\) depending on the depth window, the clearest signature that the floor is the RMS of an oscillating residual rather than a clean power law—and claim only the stable, mechanism-supported statement read off panel (a): the phase-averaged envelope decays at the Hardy rate. Equation ?? supplies its exact, classically precomputable value at any \((T,\alpha)\), opening the possibility of quantile tuning—choosing \(\alpha\) so that \(R^\star\) lands near a zero of the residual—which we leave to future work. Konno limit measures, the RMSE-vs-\(M\) comparison at fixed \(T\), and the \(\alpha^\star\) sweep are reported in Appendix A.

The classical 2D random walk paired with the same disk estimator (grey triangles in Fig. 2) yields the slower scaling \(T^{-0.85}\), broadly in agreement with the diffusive prediction \(T^{(\gamma-2)/2} \approx T^{-3/4}\) of Eq. 8 , and its floor remains \(5\)–\(10\times\) above all three quantum coins across the full depth range.

The flatness of the QW curve in Fig. 1 is the operational content of the bias–variance decoupling of Sec. 2.6: to reach any precision \(\varepsilon\) strictly above the bias floor, the lattice-quantile estimator needs only a fixed, \(\varepsilon\)-independent \(M\), whereas the sample-mean class needs \(M\sim\varepsilon^{-2}\). The sample-count gap therefore widens as \(\varepsilon^{-2}\) as the target tightens—reaching \(\sim10^4\times\) at \(\varepsilon=10^{-3}\)—until \(\varepsilon\) approaches the floor, below which no \(M\) suffices and the walk depth \(T\) must instead be increased. This is the direct consequence of the ballistic variance scaling \(\langle r^2\rangle_{\mathrm{QW}}\sim T^2\) versus the diffusive \(\langle r^2\rangle_{\mathrm{RW}}\sim T\) of Lemma 1.

The generalised estimator of Proposition 4 extends the \(\pi\) benchmark to arbitrary bounded regions. We validated it on eight 2D shapes spanning four boundary classes—convex smooth, non-convex smooth, doubly connected, and piecewise smooth with corners—obtaining sample-count reductions of \(10^{2}\)–\(10^{7}\times\) over uniform MC at \(M=100\), \(T=200\); the full table and the boundary-class discussion (including the cardioid cusp and polygonal-corner caveats) are given in Appendix 7. Of particular relevance to what follows is the annulus: its doubly-connected, non-convex super-level geometry is the same class that arises for the radial densities of excited harmonic-oscillator states (Sec. 3.3), and the QW estimator handles it at RMSE \(5.4\times10^{-4}\) (\(2.8\times10^5\times\) over MC), providing the empirical foundation for the energy-ladder reconstruction of Sec. 3 C.

3.2 Error propagation across a Gaussian CDF↩︎

The disk estimator returns \(\hat{\pi}\) from the \(\alpha\)-quantile of a single batch of \(M\) radial measurements. Because the Cavalieri identity (Proposition 5) writes any rotationally symmetric integral as \(\pi\cdot Q\) with \(Q\) precomputed classically, a single quantum-walk run fixes an entire one-parameter family of integrals at no additional quantum cost. We state plainly what this does and does not mean before illustrating it on the 2D Gaussian cumulative distribution \[P(R) = \iint_{x^2+y^2\le R^2}\frac{e^{-r^2/\sigma^2}}{\pi\sigma^2}\,dx\,dy = 1 - e^{-R^2/\sigma^2}, \label{eq:gauss-cdf}\tag{9}\] the 2D analogue of the normal CDF and the ground-state probability of the isotropic harmonic oscillator. Writing \(P_\sigma(R)=\pi\cdot Q_\sigma(R)\) with \(Q_\sigma(R)=(1-e^{-R^2/\sigma^2})/\pi\) a classically tabulated grid, a single walk yields \(\hat{P}_\sigma(R)=\hat{\pi}\cdot Q_\sigma(R)\) across all radii and all widths \(\sigma\) simultaneously: the same measured \(\hat{\pi}\) reconstructs the entire one-parameter family. The shape of every curve resides entirely in the classical prefactor \(Q_\sigma(R)\); the quantum measurement contributes exactly one number, \(\hat{\pi}\), and each \(\hat{P}_\sigma(R)\) is a one-parameter rescaling of a known curve. The reconstruction errors at different radii and widths are therefore perfectly correlated—this follows by algebraic construction: \(\hat{P}_\sigma(R)-P_\sigma(R)=Q_\sigma(R)(\hat{\pi}-\pi)\) exactly since \(Q_\sigma(R)\) is a deterministic scalar, so the trial-by-trial correlation equals \(1.0000\) identically. What this part demonstrates is consequently not curve recovery but the error budget: the single-number accuracy of each walker propagates across the entire family as \(\sigma(\hat{P}_\sigma(R)) = Q_\sigma(R)\,\sigma(\hat{\pi})\), with a deterministic, analytically known multiplier.

Having benchmarked the estimator against the sample-mean paradigm (MC and QMC) in Sec. 3 A, Secs. 3 B and C focus on the comparison that isolates the quantum resource itself: the same lattice-quantile estimator driven by a ballistic quantum walker versus a diffusive classical walker. Figure 3 reconstructs three members of the family (\(\sigma=0.7,1.0,1.5\)) from one walk at \(T=200\), \(M=500\), 50 trials. The QW recovers the sigmoidal shape of every curve visually (panel a), and the absolute error (panel b) separates QW from RW cleanly at each \(\sigma\): the QW reconstruction stays below \(3.3\times10^{-5}\) across the full radial range, while the RW reaches \(4.8\times10^{-4}\)—a factor of \(\sim15\times\) worse, which by construction equals the ratio of the underlying \(\hat{\pi}\) errors of the two walkers at this \((T,M)\) (a finite-\(M\) total-error ratio, somewhat larger than the \(\sim7\)–\(10\times\) asymptotic bias-floor ratio of Sec. 3.1). Because the error of \(\hat{P}_\sigma(R)\) is \(Q_\sigma(R)\cdot\mathrm{RMSE}(\hat{\pi})\), it grows with \(R\) exactly as the classically known prefactor dictates, vanishing at \(R\to0\) where \(Q_\sigma\to0\). This error-propagation structure is exploited more substantively in Sec. 3 C, where the multipliers \(Q_n\) carry the level dependence of a physical spectrum.

3.3 Energy ladder of the 2D quantum harmonic oscillator↩︎

We now apply the framework to a genuine spectroscopic task: reconstructing the low-lying energy ladder of the 2D isotropic harmonic oscillator, \(E_n = \hbar\omega(n+1)\) for \(n=0,1,2,3\) (ground state through third excited state). The lowest radial mode of angular momentum \(|\ell|=n\) has probability density \(|\psi_n(r)|^2 = (\pi\sigma^2 n!)^{-1}(r/\sigma)^{2n}e^{-r^2/\sigma^2}\), which is rotationally symmetric and radially single-peaked. By the virial theorem the energy is \(E_n = m\omega^2\langle r^2\rangle_n\) with \[\langle r^2\rangle_n = \iint r^2\,|\psi_n(r)|^2\,dx\,dy = \pi\cdot Q_n, \label{eq:r2-cavalieri}\tag{10}\] where the integrand \(r^2|\psi_n|^2\) has annular (doubly-connected, non-convex) super-level sets for \(n\ge1\), precisely the boundary class validated in Appendix 7 (Table ¿tbl:tab:shapes?). The Cavalieri prefactors \(Q_n\) are computed classically once; we find \(\pi Q_n = (n+1)\sigma^2\) to machine precision, recovering the exact \(\langle r^2\rangle_n = (n+1)\sigma^2\).

Each \(\hat{E}_n = \hat{\pi}\cdot Q_n\) is estimated over 50 trials at \(T=200\), \(M=500\), for the QW and the RW—again the identical estimator with only the walker exchanged. Figure 4 (a) shows the resulting energy ladder. Both walkers recover the linear spectrum, but the error bars—magnified \(30\times\) for visibility—are markedly tighter for the QW. A weighted linear fit \(E_n = \hbar\omega\,n + E_0\) extracts the level spacing and zero-point energy: \[\begin{align} \text{QW:}\quad & \hbar\omega = 1.0000(2),\quad E_0 = 1.0000(2),\nonumber\\ \text{RW:}\quad & \hbar\omega = 1.001(2),\quad\;\; E_0 = 1.001(2),\nonumber \end{align}\] in units where the true values are \(\hbar\omega = E_0 = 1\). The QW determines the level spacing with a \(10\times\) smaller statistical uncertainty than the RW at identical walk depth and sample count.

A distinctive feature, shown in Fig. 4 (b), is that the statistical error grows linearly with level index: \(\sigma(\hat{E}_n) = Q_n\,\mathrm{RMSE}(\hat{\pi}) \propto (n+1)\), a direct, parameter-free consequence of the error-propagation structure—each walker’s \(\hat{\pi}\) accuracy enters every level through the classically known multiplier \(Q_n\). A linear fit of \(\sigma(\hat{E}_n)\) versus \(n\) gives a slope of \(1.90\times10^{-4}\,\hbar\omega\) per level for the QW and \(1.99\times10^{-3}\,\hbar\omega\) per level for the RW—a ratio of \(10.5\), so the QW error budget widens an order of magnitude more slowly as one climbs the spectrum. The ballistic advantage thus compounds level by level rather than appearing only in a single figure of merit. The framework thus delivers not just a single estimate but a controlled, analytically predictable error budget across the full spectrum.

3.3.0.1 Beyond closed-form targets.

The harmonic ladder is exactly solvable, so the above is a controlled demonstration. To exhibit the intended use case—moments of a known radial density that lack an elementary closed form—we repeat the construction for the anharmonic family \(f_\lambda(r) \propto e^{-r^2/\sigma^2 - \lambda r^4}\), whose second moment \(\langle r^2\rangle_\lambda\) has no elementary expression. The integrand \(r^2 f_\lambda(r)\) again has annular super-level sets, and the annular Cavalieri identity, Eq. ?? , gives \(\langle r^2\rangle_\lambda = \pi\cdot Q_\lambda\) with \(Q_\lambda\) a one-dimensional classical quadrature; we verify \(\pi Q_\lambda\) against direct two-dimensional quadrature to machine precision (\(\langle r^2\rangle_{0.25}=0.638968\), \(\langle r^2\rangle_{1.0}=0.416353\) at \(\sigma=1\)). The propagated accuracies follow immediately from the measured \(\hat{\pi}\) statistics of Fig. 3: the quantum walker determines \(\langle r^2\rangle_\lambda\) to a relative error of \(1.7\times10^{-4}\), the diffusive walker to \(2.1\times10^{-3}\), at identical \(T=200\) and \(M=500\)—the same order-of-magnitude separation as the harmonic ladder, now for a target with no analytic shortcut.

4 Discussion↩︎

The numerical results of Sec. III match the mechanistic picture of Sec. II F. Bias and variance decouple, with the bias floor governed by walk depth \(T\) and the variance term by sample count \(M\). The phase-averaged envelope slopes are \(-1.502\) (Grover), \(-1.489\) (Hadamard), and \(-1.506\) (Fourier) (Fig. 2), averaging \(-1.499\approx -3/2\), which matches the Hardy conjecture to three decimal places. Raw single-radius fits can yield slopes as anomalous as \(+2.7\) or as steep as \(-4.4\) due to the oscillating Gauss residual sampling different phases at successive depths (Fig. 2 (b,c)); these are phase artefacts, not coin-dependent asymptotics, and disappear under phase-averaging. The classical 2D random walk paired with the same disk estimator scales as \(T^{-0.85}\), consistent with the diffusive prediction of Eq. 8 to within \(\sim 13\%\) (computed using the Hardy-conjectured \(\gamma=1/2\), giving \(T^{-3/4}\); the Huxley bound \(\gamma\le131/208\) gives \(T^{-0.685}\), a \(\sim24\%\) discrepancy). The resulting sample-count ratio \(M_{\mathrm{MC}}/M_{\mathrm{QW}} = \Theta(\varepsilon^{-2})\) for any \(\varepsilon\) strictly above the QW bias floor characterises what we call the framework’s ballistic acceleration (defined relative to uniform dart-throw MC; the corresponding ratio relative to scrambled Sobol QMC is \(\sim 500\times\) by order-of-magnitude extrapolation). It is structurally distinct from the \(\Theta(\varepsilon^{-1})\) oracle-call scaling of QAE [5], [6], [9] and from the \(\sqrt{\Delta_c}\) mixing-time scaling of Szegedy walks [10]–[12], [27]. Whether this \(\varepsilon^{-2}\) ratio constitutes a complexity-theoretic separation in any precise sense is a separate question that we do not address here.

The 2D coined DTQW uses \(2 + 2\lceil\log_2(2T+3)\rceil\) qubits with no ancillae; for \(T=200\) this corresponds to \(\sim 20\) qubits, which is within the range of current superconducting devices by qubit count. Gate depth, however, is more demanding: the asymptotic \(O(T \cdot n_{\rm pos}^2)\) cx-count of QFT-based modular adders [28], [29] implies \(\sim\!10^6\)–\(10^7\) cx gates per shot at \(T=200\), far beyond what current hardware can execute reliably. The estimator uses only projective measurement of the position register followed by classical post-processing of \(N(\hat{R})/\hat{R}^2\); it does not require quantum phase estimation, controlled oracles, or controlled rotations of the integrand. We use “oracle-free” in the QAE sense: no unitary \(\mathcal{A}\) that prepares a state encoding the integrand is required, and consequently no controlled \(\mathcal{A}\) inside a QPE register. The cost contrast with amplitude-estimation-based Monte Carlo on near-term hardware [7], [8] is precisely that those primitives have proven costly to compile. Near-term demonstrations of DTQWs suggest a path toward larger walk depths: Razzoli et al. [20] demonstrated circuit-efficient DTQWs on IBM devices beyond the few-step regime, and continued improvements in two-qubit gate fidelity could reduce per-step error accumulation. However, extrapolating from these short-depth demonstrations to the \(T\approx100\)–\(200\) regime required here involves a \(\sim\!15\times\) increase in walk depth, accompanied by a \(\sim\!10^2\times\) increase in gate count from the combined growth in \(T\) and \(n_{\rm pos}\) (from \(\sim 4\) to \(\sim 9\) position qubits, \(T\,n_{\rm pos}^2\) scaling); whether ballistic spreading is preserved under this level of noise accumulation is an open experimental question. We emphasise that strong conclusions about hardware feasibility cannot be drawn without explicit noise-model analysis at these gate depths.

4.0.0.1 Sample count versus total circuit cost.

The sample-count ratios reported throughout this paper compare the number of projective measurements \(M\) required by each method to reach a given RMSE, and do not include quantum circuit execution cost. As noted above, the asymptotic gate-count estimate implies \(\sim\!10^6\)–\(10^7\) cx gates per shot at \(T=200\), compared with \(O(1)\) floating-point operations per MC or QMC sample. The sample-count reduction of \(\sim 100\) versus \(\sim 7\times 10^6\) therefore quantifies an advantage in measurement complexity rather than wall-clock time or total gate count. Quantifying the end-to-end trade-off on fault-tolerant hardware requires reliable circuit synthesis for \(T=200\) and a concrete noise model, both of which are outside the scope of the present paper.

4.0.0.2 From a single number to a spectrum.

The progression of Sec. 3 highlights a structural feature of the lattice-quantile framework that distinguishes it from sample-mean estimators. Because every target integral is written as \(\pi\) times a classically precomputed prefactor (Proposition 5, Corollary 1), one batch of measurements feeds an arbitrary number of derived observables, each inheriting the same \(\hat{\pi}\) bias floor through a known multiplier. Section 3 B makes this explicit on a Gaussian CDF, the curve shape is entirely classical, and the measurement fixes one scalar, and Sec. 3 C turns the same structure into a spectroscopic tool: the 2D harmonic-oscillator energy ladder \(E_n=\hbar\omega(n+1)\) is recovered with a level spacing \(\hbar\omega = 1.0000(2)\), a \(10\times\) tighter determination than the diffusive RW baseline at identical \(T\) and \(M\). The error budget is itself analytically predictable, \(\sigma(\hat{E}_n)=Q_n\,\mathrm{RMSE}(\hat{\pi})\propto(n+1)\), a transparency that sample-mean Monte Carlo, whose error depends on the integrand variance at each \(R\), does not offer. We stress that the harmonic oscillator admits closed-form energies, so Sec. 3 C is a controlled demonstration of error propagation rather than a claim of advantage over analytic methods; the value of the construction lies in problems where the radial density is known but its moments lack a closed form, such as anharmonic or multi-well potentials.

5 Conclusion↩︎

This paper has introduced a lattice-quantile framework joining the ballistic spreading of 2D coined quantum walks (Watabe–Konno limit [15]–[17]) with the Hardy–Huxley solution of the Gauss circle problem [18], [19]. Because the estimator is a deterministic lattice count applied to a rank statistic rather than a sample mean, the Cramér–Rao bound on sample-mean estimators does not directly constrain it. What it does constrain is the residual sampling fluctuation of the rank statistic, which decays as \(M^{-1/2}\) and is sub-dominant to the depth-controlled bias floor in our regime. At walk depth \(T=200\), \(M=50\) measurements suffice to reach \(\mathrm{RMSE}\approx 6\times 10^{-4}\) for \(\pi\), and the Grover coin saturates to \(\mathrm{RMSE}\approx5.0\times 10^{-4}\) at \(M=10^4\). This represents a measurement-count reduction of \(\sim 7\times 10^{4}\times\) over uniform Monte Carlo at \(M_{\mathrm{QW}}=100\) (the floor is already reached at \(M=50\); rising to \(\sim 1.1\times 10^{5}\times\) at the Grover asymptotic floor, \(\mathrm{RMSE}\approx5.0\times10^{-4}\)) and \(\sim 500\times\) over scrambled Sobol QMC by order-of-magnitude extrapolation, while the classical random walk paired with the identical estimator stalls at a \(\sim 10\times\) higher bias floor. The accuracy of a single batch of measurements propagates through classically known multipliers across an entire Gaussian cumulative distribution. The framework recovers the 2D harmonic-oscillator energy ladder with a level spacing \(\hbar\omega=1.0000(2)\), ten times tighter than the diffusive baseline, and the error budget across levels is analytically predictable. These ratios compare measurement counts at fixed precision and do not include quantum circuit execution cost, which we estimate at \(\sim 10^{6}\) to \(10^{7}\) cx gates per shot at \(T=200\).

5.0.0.1 Physical applications.

The lattice-quantile construction targets precisely the class of quantities that reduce to counting integer points inside a smoothly bounded region, which recur throughout physics. The most direct is the density of states: in the free-electron and photon-gas pictures the cumulative number of modes below energy \(E\) is the count of reciprocal-lattice points inside a sphere of radius \(k(E)\), and its derivative is the density of states [30]. The leading volume term is the Weyl law [31], and the correction, governed by the curvature of the energy surface, is exactly the lattice discrepancy our estimator computes, suggesting a route to mode counting for free Bose and Fermi gases and for quantum billiards [32]. The same structure underlies microcanonical phase-space volumes \(\Omega(E)=\#\{\mathbf{n}: \varepsilon(\mathbf{n})\le E\}\), whose logarithm is the Boltzmann entropy; whenever the energy surface \(\varepsilon(\mathbf{n})=E\) is smooth and convex, the generalised estimator of Proposition 4 applies directly, and the anharmonic example of Sec. 3 C is a first step in this direction.

5.0.0.2 Higher dimensions.

The \(d\)-dimensional Gauss-sphere generalisation [26], [33] is particularly favourable: while the planar (\(d=2\)) and \(d=3\) discrepancy exponents remain the famous open problems, for \(d\ge4\) the exponent is rigorously pinned at the volume-shell value by classical results on representations by quadratic forms [33]. A ballistic walk on \(\mathbb{Z}^d\) would therefore inherit a provable, rather than conjectural, bias-floor scaling in four or more dimensions, an unusual situation in which the higher-dimensional problem is theoretically cleaner than the planar one studied here.

5.0.0.3 Hardware and algorithmic outlook.

Although the gate-model circuit depth required for \(T\sim200\) is beyond current superconducting devices, Photonic platforms natively support quantum walks without compiling modular adders: time-multiplexed fibre loops implement DTQWs and have reached over ten steps in two dimensions [34], while integrated waveguide arrays realise continuous-time quantum walks at scale [35]. These platforms offer a complementary route to the walk depths our estimator needs. On the theory side, natural next steps are a noise-threshold analysis, decoherence drives DTQWs toward the diffusive regime [36] and would progressively lift the bias floor toward the classical one, adaptive \(T\)-scheduling, quantile tuning to exploit the oscillatory floor structure of Eq. ?? , and a closer connection to quantum-walk search, where ballistic spreading is the same underlying algorithmic resource [37].

5.0.0.4 Open theoretical questions.

Appendix 8 states formal counterparts of the following four items conditional on a single boundary-density conjecture; rigorous treatment under standard hypotheses is open. (i) A rigorous proof of the inverse-square-root edge behaviour of the Grover Konno radial marginal (Conjecture 6), together with an explicit value of the constant \(A_G\). (ii) Removal of the truncation \(v_0\) in Lemma 4, i.e.a Bahadur–Kiefer-type representation that handles the \(\delta\)-localisation atom of the Grover measure directly. (iii) An explicit numerical value of the Huxley constant \(K_H\) in Prop. 7, computable from the proof of [19] but not given in closed form there. (iv) A \(d\)-dimensional Konno-type weak limit theorem for the Grover coin in \(d\ge 3\), completing the bias-floor scaling of Prop. 8 for that coin. The Hardy line itself remains, of course, conjectural; only the Huxley bound \(\gamma\le 131/208\) is rigorous, and all numerical results assume noiseless evolution.

With these caveats, the lattice-quantile construction offers a third entry in the catalogue of quantum-walk Monte Carlo techniques, complementary to amplitude-estimation and Szegedy-walk approaches. More broadly, harnessing a quantum platform to probe a classical problem in analytic number theory—here the Gauss circle problem—parallels recent quantum approaches to the Riemann zeta zeros [38], suggesting number theory as a structured resource for quantum estimation.

Acknowledgments↩︎

EJK thanks the National Center for Theoretical Sciences and National Yang Ming Chiao Tung University in Taiwan for their support under Grant No. NSTC 114-2112-M-A49-036-MY3.

6 Comprehensive coin operator comparison↩︎

This appendix collects the definitions and supplementary numerical comparisons of the three coin operators. The Grover coin \(G\) is defined in Eq. 4 of the main text. The separable Hadamard coin acts as \(H\otimes H\) on the two coin qubits, \[H\otimes H = \frac{1}{2}\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -1 & 1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & -1 & 1 \end{pmatrix}, \label{eq:hadamard}\tag{11}\] corresponding to two independent 1D Hadamard walks on the \(x\) and \(y\) axes. The 4-point Fourier coin is \[F = \frac{1}{2}\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{pmatrix},\qquad F_{jk}=\tfrac12 e^{2\pi i jk/4}, \label{eq:fourier}\tag{12}\] which produces the most concentrated wavefront near the velocity-disk boundary and hence the smallest variance of the radial \(\alpha\)-quantile (Proposition 3).

6.0.0.1 Limit measure structure.

Figure 5 (a) and Fig. 6 together characterise the Konno limit measure of each coin. The radial CDF (panel a) reveals two qualitatively distinct features: (i) the localisation spike of the Grover coin at \(r \approx 0\), contributing \(\sim 35\%\) of the total probability mass within \(|v| \le 0.05\); and (ii) the sharp boundary edge of all three coins at \(|v| \approx \sqrt{\nu_C^2}\), where the lattice quantile \(\hat{R}(\alpha)\) becomes well-defined for high \(\alpha\).

6.0.0.2 Bahadur–Kiefer caveat for the Grover coin.

The Bahadur–Kiefer representation invoked in Sec. 2.6 assumes the sampling distribution is continuous with positive density at the quantile \(\alpha\). The Hadamard product measure \(\mu_H=\mu_{1\text{D}}\otimes\mu_{1\text{D}}\) satisfies this (Lemma 1(b)). The Fourier coin satisfies the continuity-and-positive-density hypothesis numerically at \(\alpha\ge 0.9\) (Fig. 5 (a)), but lacks a rigorous weak-limit statement (Lemma 1(c)), so the BK invocation for the Fourier coin is on the same numerical footing as the limit measure itself. The Grover coin exhibits a \(\delta\)-function–like localisation at the origin (the central spike in Figs. 5 (a) and 6); for \(\alpha\ge 0.9\) the quantile sits in the smooth, boundary-concentrated component of the measure and the localisation atom does not contribute, so the standard BK representation should be regarded as a heuristic motivating the \(1/M\) variance scaling rather than a strict theorem, with a rigorous truncated-BK lemma left to future work (open question (ii) of Sec. 5). The empirical \(1/M\) scaling observed in Fig. 1 is consistent with this heuristic at the values \(\alpha\ge 0.9\) used in practice.

6.0.0.3 Ballistic constants \(\nu_C^2\).

The variance ratio \(\nu_C^2 = \langle r^2 \rangle / T^2\) is independent of \(T\) in the asymptotic regime (Lemma 1). Empirical values at \(T \in \{50, 100, 200, 400\}\) (Table ¿tbl:tab:coin-summary?) are: Grover \(\nu_G^2 = 0.364\), Hadamard \(\nu_H^2 = 0.586\), and Fourier \(\nu_F^2 = 0.402\). All three are well below the support boundary \(1/2\) predicted by the velocity-disk constraint \(\{|v| \le 1/\sqrt{2}\}\), except for Hadamard whose value exceeds \(1/2\) due to its support extending into the corners of the square \(|v_x|, |v_y| \le 1/\sqrt{2}\) (visible in Fig. 6, middle panel). All values are stable to four significant digits across the four \(T\) values, confirming the asymptotic regime.

6.0.0.4 Practical recommendation.

For the disk estimator at \(T = 200\), \(M \in [50, 10^4]\), all three coins yield comparable RMSE within a factor of two. We recommend the Grover coin for its lattice symmetry and well-studied limit measure [17], the Fourier coin when the steepest bias-floor scaling is desired (e.g., for \(T \ge 400\) where it outperforms the Hardy line), and the Hadamard coin for the simplest experimental implementation as a separable two-qubit walk on \(x\) and \(y\) axes. The Grover coin is used as the default in the main body of this paper.

7 Generalised area estimation across boundary classes↩︎

We validate Proposition 4 on a representative set of 2D shapes spanning four boundary classes: convex smooth (disk, ellipse, super-ellipse), non-convex smooth (rosette, cardioid), doubly-connected (annulus), and piecewise smooth with corners (square, L-shape). Each shape is defined by a region equation \(D = \{(x,y) : \text{condition}\}\) with closed-form true area; the lattice count \(N_D(L)\) is enumerated classically. Table ¿tbl:tab:shapes? reports the RMSE at \(M=100\) for both methods across all eight shapes, alongside the sample count \(M_{\mathrm{MC}}^\star\) that classical Monte Carlo would require to match the QW precision.

7.0.0.1 Beyond-Krätzel regimes.

Several shapes violate the strict hypothesis of Lemma 3. We classify them by which lattice-counting bound the literature supplies in place of \(\gamma_D\):

• Convex piecewise-\(C^2\) boundaries with isolated curvature degeneracies (super-ellipse \(x^4+y^4\le1\), with vanishing curvature on the axes; cardioid, with a cusp at the origin). The disk floor of Lemma 2 is recovered up to a sub-leading correction; empirically the floor sits at \(\sim 2\)–\(3\times\) the strictly-curved disk floor.

• Polygonal boundaries (square, L-shape). Each corner contributes an \(O(1)\) lattice-counting error, giving floor exponent \(T^{-1}\) rather than \(T^{\gamma-2}\) [26], consistent with the \(600\) to \(5800\times\) reductions in the last two rows.

• Doubly-connected annular boundary. The corollary applies via the disk floor, without further hypothesis (Cor. 1).

Fractal boundaries with Hausdorff dimension \(d_H>1\) would be expected to weaken the framework further; we did not test these. The five shapes in (i)–(iii) lie outside the rigorous reach of Proposition 4, and the corresponding entries of Table ¿tbl:tab:shapes? should be read as numerical evidence.

All eight Table ¿tbl:tab:shapes? entries use the Grover coin only. Whether the same range of sample-count ratios is preserved with Hadamard or Fourier coins is plausible but has not been verified here.

8 Supporting theoretical results↩︎

This appendix collects formal statements supporting Sec. II and lists the points at which further work is required.

8.0.0.1 Notation.

We write \(\mu_C\) for the Konno limit measure of coin \(C\in\{G,H,F\}\) and \(g_C(v)\) for its radial marginal density (defined where the measure has an absolutely continuous part). The radial CDF is \(F_C(v)=\int_0^v g_C(u)\,du\) and the population \(\alpha\)-quantile is \(v^\star(\alpha)=F_C^{-1}(\alpha)\).

8.1 Closed-form treatment of \(\kappa(\alpha)\) for the Grover coin↩︎

The exponent \(\kappa(\alpha)\in(0,1]\) entering Proposition 1 parametrises the \(T\)-scaling contribution of \(g_G(v^\star)^{-1}\) inherited from the Bahadur–Kiefer representation of \(\hat{R}(\alpha)\). We make the dependence explicit.

Motivation. The inverse-square-root form is the standard 1D Konno edge behaviour [15], [16] and is expected to carry over to the radial marginal of the 2D Grover walk by the stationary-phase analysis underlying [17]. Watabe et al.proved the weak-limit theorem but did not write the radial-marginal edge expansion in closed form. Verifying the form rigorously and computing \(A_G\) are open problems flagged in Sec. 5.

8.2 Truncated Bahadur–Kiefer representation for atomic Konno measures↩︎

The Grover Konno measure \(\mu_G\) has a localisation atom at the origin of mass \(p_{\mathrm{loc}}\approx 0.35\) (Appendix 6). The standard Bahadur–Kiefer representation assumes a continuous density at the quantile; this fails at the atom. The following truncation argument suffices in the regime used in this paper.

8.3 Explicit upper bounds on \(C_1\) and \(C_2\) for the Grover coin↩︎

8.4 Higher-dimensional generalisation↩︎

For \(d\ge 3\), the lattice-counting residual \(|N_d(R)-V_d R^d|\) is bounded by \(O(R^{\gamma_d})\) (possibly with logarithmic factors) with \(V_d=\pi^{d/2}/\Gamma(d/2+1)\) the unit-ball volume. The state of the art on \(\gamma_d\) is [26], [33], [39]:

• \(d=3\): \(\gamma_3\le 21/16\) (Heath-Brown line [39]); the optimal exponent is open.

• \(d=4\): \(|N_4(R)-V_4 R^4|=O(R^2\log^{2/3}R)\) (Walfisz [33]); \(\gamma_4=2\) up to a logarithmic factor.

• \(d\ge 5\): \(\gamma_d=d-2\) rigorously, with no logarithmic correction, by the representation theory of quadratic forms [33].

References↩︎