Let \(n_{\rm gr}\) be the number of grid points per phase-space axis and \(n_t\) the number of time steps. The classical cost of storing/evolving the 6-D distribution scales as

because phase space has six axes. The quantum algorithm instead uses

to output the power spectrum \(P_\nu(k)\) to accuracy \(\epsilon\) — an exponential reduction in the grid dependence.

• QRAM. The external CDM force must be stored in a quantum RAM with \(O(n_{\rm gr}^{3})\) entries — a device that does not yet exist at scale.
• Output only. You get \(P_\nu(k)\), not the full \(f\); reading \(2^n\) amplitudes is impossible, so the “answer” is deliberately a small summary.
• Fault tolerance. The circuits (block-encoding, QAE) assume an error-corrected machine — not today’s NISQ hardware.

Let \(\rho(\mathbf x)\) be the matter density at comoving position \(\mathbf x\) and \(\bar\rho\) its cosmic mean. The dimensionless fluctuation is

By construction \(\langle\delta\rangle=0\); \(\delta>0\) marks an overdensity (a proto-halo or filament), \(\delta=-1\) an empty void, and \(\delta\gg1\) the deeply nonlinear, collapsed regime. \(\delta\) is a single realization of a random field seeded by quantum fluctuations during inflation — so we predict its statistics, not its exact value.

The lowest-order statistic is the two-point correlation \(\xi(r)=\langle\delta(\mathbf x)\delta(\mathbf x+\mathbf r)\rangle\). Its Fourier transform is the power spectrum. Concretely, transform \(\delta\to\tilde\delta(\mathbf k)\); statistical homogeneity (no preferred origin) forces the two-point function in Fourier space to be diagonal, and statistical isotropy (no preferred direction) makes it depend only on \(k=\|\mathbf k\|\):

So \(P(k)\) is the variance of the fluctuations at wavenumber \(k\): how much structure exists at spatial scale \(\sim\!2\pi/k\). Small \(k\) = large scales (superclusters); large \(k\) = small scales (galaxies). A scale-free-ish \(P(k)\) means structure on all scales — the observed cosmic web.

Cosmologists often quote the dimensionless variance per log-\(k\),

the contribution to \(\langle\delta^2\rangle\) from each decade of scale. The neutrino imprint sits at \(k\sim0.01\text{–}1\ h\,\text{Mpc}^{-1}\) — the mildly nonlinear regime where percent-level modelling is both possible and necessary.

Galaxy redshift surveys (SDSS, DESI, Euclid) and weak gravitational lensing directly measure \(P(k)\) (or the closely related lensing/galaxy spectra) across scales, together with the baryon-acoustic-oscillation feature that calibrates distances. Cosmological parameters — matter density, dark-energy equation of state, and crucially the neutrino mass — are inferred by comparing the measured \(P(k)\) to a predicted \(P(k)\) from theory and simulation. A biased prediction biases the inferred \(M_\nu\); hence the demand for a neutrino \(P_\nu(k)\) accurate to the same percent level as the data. That accuracy requirement is the seed of this paper.

CDM particles are heavy and, by assumption, born with negligible thermal velocity dispersion, \(\sigma_v\to0\). At every point their velocity is effectively single-valued, so they behave as a pressureless dust: gravity has no thermal pressure to fight, and they cluster on all scales down to very small \(k^{-1}\). This is why CDM reproduces the observed cosmic web so well — and why it can be treated with a simple fluid or N-body description.

A “hot” species carries large thermal velocities \(\sigma_v\). Random motion acts like a pressure that resists collapse — the kinetic analogue of the Jeans criterion. A perturbation of comoving wavelength \(\lambda\) can only grow if gravity overcomes streaming; the borderline defines the free-streaming scale. Heuristically, a particle covers a comoving distance \(\lambda_{\rm fs}\sim\int v\,dt/a\sim v/(aH)\) before the potential well can turn it around:

Above \(k_{\rm fs}\) the species does not fall into overdensities — it streams straight through and out, damping its own perturbations.

Encode the net effect in the transfer function \(T(k)=\sqrt{P_{\rm with}(k)/P_{\rm without}(k)}\). For a free-streaming species \(T(k)\to1\) at small \(k\) (untouched large scales) and drops toward a suppressed plateau at large \(k\). The observable consequence is a scale-dependent step in the matter power spectrum: full power on large scales, reduced power below \(k_{\rm fs}\). Neutrinos, being only \(\lesssim1\%\) of the matter, produce a partial step of depth \(\Delta P/P\approx-8\,\Omega_\nu/\Omega_m\) (slide 6) rather than the total cut-off a pure-HDM universe would have.

Because the suppression comes from the spread of velocities, you cannot capture it with a single-velocity fluid: you need the full distribution over \(v\). That is exactly why neutrinos demand the Vlasov equation in 6-D phase space (slide 7) — and why they are so much more expensive than CDM.

In the hot early universe, neutrinos were kept in thermal equilibrium by weak interactions. As the universe expanded and cooled, the weak-interaction rate \(\Gamma\propto G_F^2 T^5\) fell faster than the expansion rate \(H\propto T^2\); at \(T\sim 1\ \text{MeV}\) (about one second after the Big Bang) they decoupled — froze out — while still ultra-relativistic (\(T\gg m_\nu\)). They therefore inherited a relativistic Fermi–Dirac momentum distribution, today redshifted to the cosmic neutrino background at \(T_\nu\simeq 1.95\ \text{K}\).

Because they froze out relativistic, their momenta redshift as \(p\propto (1+z)\) but the spread of the distribution never narrows into a cold clump. A characteristic relic speed today is

i.e. hundreds of km/s — comparable to the orbital speeds inside galaxy clusters. That is the operational meaning of “hot”: even with \(m_\nu\sim0.1\ \text{eV}\), a neutrino outruns galactic-scale potential wells and streams away rather than falling in.

The distance a neutrino covers before gravity can turn it around defines a comoving free-streaming wavenumber \(k_{\rm fs}\); below it (larger scales) neutrinos cluster with the CDM, above it (smaller scales) they do not:

This is the physical scale the whole simulation exists to resolve: the location of the suppression feature in \(P(k)\). Getting \(k_{\rm fs}\) and the depth of the dip right is how \(m_\nu\) is read off from data.

Summing over the three mass eigenstates, \(M_\nu=\sum m_\nu\), the present-day neutrino energy density is fixed by the relic number density (\(\approx 336\ \text{cm}^{-3}\)) times the mass, giving \(\Omega_\nu h^2 = M_\nu/93.14\ \text{eV}\). Dividing by the total matter density (paper Eq. 3):

The constant \(93.14\ \text{eV}\) is exactly the summed mass that would make neutrinos all of the dark matter; the measured \(M_\nu\) is far below it, so neutrinos are \(\lesssim 1\%\) of the matter. This one number — the smallness of \(\Omega_\nu/\Omega_m\) — is the physical licence for the linearization in Part 3 (neglecting neutrino self-gravity), which is what makes the quantum algorithm possible at all.

Neutrino oscillation experiments measure only mass-squared differences, \(\Delta m^2_{21}\simeq 7.5\times10^{-5}\ \text{eV}^2\) and \(|\Delta m^2_{31}|\simeq 2.5\times10^{-3}\ \text{eV}^2\), which fix the spacing of the eigenstates but not the overall scale. They imply a floor \(M_\nu\gtrsim 0.06\ \text{eV}\) (normal ordering) or \(\gtrsim 0.10\ \text{eV}\) (inverted). Laboratory endpoint measurements (KATRIN) currently bound a single effective mass to \(\lesssim 0.45\ \text{eV}\).

Cosmology is more sensitive: to leading order the small-scale suppression is

so a percent-level measurement of \(P(k)\) probes \(M_\nu\) at the \(0.01\ \text{eV}\) level — already \(M_\nu<0.12\ \text{eV}\) from Planck+BAO. But extracting this needs theory predictions of matching accuracy, i.e. a neutrino \(P_\nu(k)\) computed without N-body shot noise. That requirement is the doorway to this paper.

One can try to describe the species by velocity moments of \(f\): the density \(\rho=\int f\,d^3v\), the momentum \(\rho\mathbf u=\int \mathbf v f\,d^3v\), the stress \(\Pi=\int \mathbf v\mathbf v f\,d^3v\), and so on. Taking moments of the Vlasov equation gives fluid-like equations — but each moment’s evolution involves the next moment (density needs momentum, momentum needs stress, …). The chain never closes on its own.

For CDM the closure is trivial: the velocity is single-valued (cold), so stress and all higher moments vanish and \(\rho+\mathbf u\) suffice. For neutrinos the velocity distribution at a point is broad and skewed, and crucially it develops multiple streams: neutrinos that fell through a potential well from different directions cross, so several distinct \(\mathbf v\) coexist at the same \(\mathbf x\). A single-velocity fluid literally cannot represent this — any truncated-moment (fluid) closure introduces uncontrolled errors exactly on the scales \(k\sim k_{\rm fs}\) we care about.

The only complete, closed description keeps the full distribution \(f(t,\mathbf x,\mathbf v)\) — the number density of neutrinos at position \(\mathbf x\) and velocity \(\mathbf v\). It lives on the 6-D phase space. Its evolution (Vlasov) is closed because it tracks all velocities at once; the observable density is recovered by integrating velocity out (paper Eq. 8):

Completeness costs dimensionality: six axes instead of three. Gridding them is the \(O(n_{\rm gr}^6)\) curse (Part 2); sampling them with particles is the N-body shot-noise problem (Part 2). The phase-space description is non-negotiable for a hot species yet ruinously expensive classically — which is precisely the gap the quantum algorithm sets out to close. Everything from here is about carrying this 6-D object efficiently.

\(f(t,\mathbf x,\mathbf v)\,d^3x\,d^3v\) is the number of neutrinos that, at time \(t\), sit in the phase-space volume \(d^3x\,d^3v\) around the point \((\mathbf x,\mathbf v)\). It is a density on the 6-dimensional space of positions and velocities, not a density in ordinary 3-D space. It is non-negative everywhere, and integrating over all of phase space returns the total count,

Because gravity neither creates nor destroys neutrinos on these scales, \(N_\nu\) is a fixed number for all time — a conservation law the discretization in Part 4 will have to respect exactly.

Everything we can actually measure is a velocity moment of \(f\) — an integral that projects the 6-D object down to a 3-D field. The most important is the zeroth moment, the local density recovered by integrating velocity out (paper Eq. 8),

from which the density contrast \(\delta_\nu\) and ultimately the power spectrum \(P_\nu(k)\) are built. The first moment gives the bulk velocity \(\rho_\nu\mathbf u=\int\mathbf v f\,d^3v\), the second the velocity dispersion (pressure). The reason we cannot stop at \(\rho_\nu\) alone — the reason we keep the full \(f\) — is the moment-closure failure of Part 1: for a hot species the higher moments do not vanish and do not close.

Each neutrino traces a trajectory \(\big(\mathbf x(t),\mathbf v(t)\big)\) — a characteristic of the equation — obeying Newton's law with the gravitational force per unit mass \(\mathbf F\):

The first equation is just the definition of velocity; the second is \(\mathbf F=-\nabla\Phi\) with \(\Phi\) the gravitational potential. The distribution \(f\) is carried, unchanged in value, along each such trajectory — like a dye passively advected by a fluid. The Vlasov equation is precisely the compact statement that \(f\) is conserved along these characteristics, which is what the next slide derives.

On cosmological scales neutrinos interact only gravitationally: the weak-interaction rate is utterly negligible once they have decoupled, so no scattering event ever jumps a neutrino discontinuously in velocity. The only thing that changes \(f\) at a point is the smooth gravitational flow. This is the difference between a Vlasov equation — a first-order transport (Liouville) equation with no right-hand side — and a full Boltzmann equation carrying a collision integral \((\partial_t f)_{\rm coll}\). Collisionlessness is what keeps the equation first-order and linear-in-derivatives, and it is ultimately why the discretized operator can be made antisymmetric (Part 4) and the evolution unitary (Part 5).

"\(f\) is constant along each trajectory" is a statement about the total (convective, Lagrangian) time derivative, following a neutrino as it moves. That derivative vanishes:

Physically this is incompressibility of the flow in phase space: a comoving blob of neutrinos may be stretched and sheared, but its 6-D volume — and hence the density within it — is preserved, because the Hamiltonian flow \((\dot{\mathbf x},\dot{\mathbf v})=(\mathbf v,\mathbf F)\) has zero phase-space divergence, \(\partial_{\mathbf x}\!\cdot\!\mathbf v+\partial_{\mathbf v}\!\cdot\!\mathbf F=0\) (the first term is zero since \(\mathbf v\) is independent of \(\mathbf x\); the second since \(\mathbf F\) depends on \(\mathbf x\), not \(\mathbf v\)).

The convective derivative unpacks into a partial-time term plus advection in each of the six axes. Substituting the equations of motion \(\dot{\mathbf x}=\mathbf v\) and \(\dot{\mathbf v}=\mathbf F\):

Each dot product runs over the three components of that vector, so this single line is really seven scalar terms: one time derivative, three spatial, three velocity. Nothing has been approximated — this is an exact rewriting of the conservation statement.

Read the three terms as three physical processes. \(\partial_t f\): the local change of the distribution at a fixed phase-space point. \(\mathbf v\cdot\partial_{\mathbf x}f\): streaming — particles drift through position at whatever velocity they carry, so a feature is transported in \(\mathbf x\) at rate \(\mathbf v\). \(\mathbf F\cdot\partial_{\mathbf v}f\): forcing — gravity accelerates particles, sliding the distribution through velocity at rate \(\mathbf F\).

These same three terms recur, unchanged in structure, all the way to the quantum circuit. Streaming and forcing become the two operators of Part 3, the two matrix blocks of Part 4, and the two factors of each Trotter step in Part 6; the vanishing right-hand side (collisionlessness) is what will let the generator be made antisymmetric and the evolution unitary. If you understand this boxed equation, the rest of the paper is a sequence of faithful re-encodings of it onto ever-more-exotic hardware.

Gravity is conservative, so the force per unit mass derives from a scalar potential, \(\mathbf F=-\nabla_{\mathbf x}\Phi\), and \(\Phi\) is fixed by the total mass through Poisson's equation,

The mean \(\bar\rho_{\rm tot}\) is subtracted because in an expanding, on-average-homogeneous universe only fluctuations source peculiar gravity — a uniform density produces no net force. So \(\Phi\) responds to the density contrast \(\delta=(\rho-\bar\rho)/\bar\rho\), not to the density itself.

The right-hand side is \(\rho_{\rm tot}=\rho_{\rm CDM}+\rho_{\rm b}+\rho_\nu\) — cold dark matter, baryons, and neutrinos. CDM plus baryons carry essentially all the mass; with \(\Omega_m\simeq0.31\) and the neutrino fraction \(\Omega_\nu/\Omega_m\lesssim1\%\) (next Part), the potential in and around collapsing halos is set overwhelmingly by the cold component. The neutrinos feel this potential fully — they just contribute almost nothing back to it.

If the neutrinos sourced their own \(\Phi\), the force would depend on their own density \(\rho_\nu=\int f\,d^3v\), so the forcing term would read \(\big(\!-\nabla\Phi[\rho_\nu]\big)\cdot\partial_{\mathbf v}f\) — the coefficient of \(\partial_{\mathbf v}f\) is itself built from \(f\). That makes the term quadratic in \(f\), and the coupled Vlasov + Poisson system genuinely nonlinear. This Vlasov–Poisson problem is the hard, self-gravitating case (the plasma-physics quantum-Vlasov literature lives here); its nonlinearity is exactly what obstructs a clean, unconditional quantum speedup, since unitary evolution is linear.

The way out is physical, not mathematical: because neutrinos are under a percent of the matter, the neutrino piece of \(\Phi\) is negligible and the potential is set by CDM alone, external to the neutrino dynamics. Freezing \(\mathbf F\to\mathbf F_{\rm CDM}(t,\mathbf x)\) as a prescribed function severs the feedback loop, turning the quadratic Vlasov–Poisson system into a linear transport equation. That single move — justified quantitatively in Part 3 — is the linchpin that makes the entire quantum algorithm possible.

Discretize each of the six phase-space axes \((x,y,z,v_x,v_y,v_z)\) with \(n_{\rm gr}\) points. Because the axes are independent, the number of grid cells is the product,

The exponent 6 is the whole problem. A modest \(n_{\rm gr}=100\) already gives \(10^{12}\) cells; bump each axis to \(n_{\rm gr}=200\) and the count jumps by \(2^6=64\times\) to \(6.4\times10^{13}\). This is the sixth-order scaling the paper flags as making a Vlasov simulation heavier than the 3-D N-body simulation it would replace.

Storing \(f\) needs \(O(n_{\rm gr}^{6})\) numbers; at \(n_{\rm gr}=100\) that is \(10^{12}\) double-precision values, roughly 8 terabytes for a single snapshot — and one must hold and update it every step. Evolving it \(n_t\) time steps costs \(O(n_{\rm gr}^{6}\,n_t)\) operations, and the paper notes the force coefficients alone demand \(O(n_{\rm gr}^{6}\,n_t)\) queries. This exponential-in-dimension blow-up is the textbook curse of dimensionality: adding resolution multiplies cost by a sixth power, so accuracy is bought at a ruinous rate.

Three of the six axes are velocity, and for a hot species the distribution is broad there, so those axes need real resolution too — you cannot cheat them down to one point as you can for cold, single-velocity CDM. State-of-the-art classical Vlasov solvers manage hundreds of thousands of grid points per position axis but only tens per velocity axis, and even that is at the edge of what the largest supercomputers allow. The room to improve accuracy by refining the grid is therefore severely limited — which is precisely the gap a quantum method might open.

A quantum computer stores \(f\) in the amplitudes of its qubits, and \(n\) qubits carry \(2^n\) amplitudes. Encoding all \(N_{\rm gr}=n_{\rm gr}^{6}\) grid values therefore needs only

versus the \(10^{12}\) classical numbers — an exponential compression of storage from terabytes to a few dozen qubits. Memory, though, is only the first hurdle: the algorithm must also perform the evolution and the read-out without ever re-materializing the \(n_{\rm gr}^{6}\) amplitudes, and whether it truly avoids that price is the substance of Parts 6–8.

N-body replaces the continuous \(f\) with \(N_p\) "superparticles" and integrates Newton's equations for each (paper Eq. 2),

Rather than storing values on a fixed 6-D lattice, it samples phase space with a swarm of particles that live only where the matter is — an adaptive, sparse representation that dodges the \(n_{\rm gr}^6\) grid. A common setting takes \(N_p=n_p^3\) with \(n_p\sim O(10^2)\), i.e. of order a million to a billion particles, which is heavy but tractable on a supercomputer.

For cold dark matter the velocity at each point is essentially single-valued, so a handful of particles per cell captures the local state faithfully — the sampling error is small precisely because there is little velocity spread to resolve. This is the standard, well-tested workhorse of structure formation, and it is exactly the tool the paper reuses to precompute the CDM gravity \(\mathbf F_{\rm CDM}\) (Part 3).

To represent a broad velocity distribution, every position must be populated by many particles spanning many velocities. If you use too few, the estimate of \(\rho_\nu=\int f\,d^3v\) is dominated by shot noise — spurious \(\sim\!1/\sqrt{N}\) fluctuations from the discreteness of the sampling. For the fast, low-density neutrino component this noise is severe: the real signal is a percent-level suppression of \(P(k)\), and shot-noise fluctuations of comparable size can swamp it unless an impractically large number of neutrino superparticles is used.

Solving the Vlasov PDE directly carries no particles, hence no shot noise at all — it pays instead with the curse of dimensionality. So the two methods trade one intractability for another: N-body trades storage for noise, direct Vlasov trades noise for storage. The paper's wager is that a fault-tolerant quantum computer can defang the direct method's storage cost (via the exponential amplitude compression), delivering a clean, noise-free neutrino power spectrum that N-body cannot cheaply provide.

From Part 1 (paper Eq. 3), the ratio of the neutrino energy fraction to the total matter fraction is fixed by the summed neutrino mass \(M_\nu=\sum m_\nu\):

With current bounds \(M_\nu\lesssim0.12\,\text{eV}\) (Planck+BAO) and the oscillation floor \(M_\nu\gtrsim0.06\,\text{eV}\), the ratio sits firmly below one percent. So in mean density the neutrinos are a rounding error on the cold matter.

A small mean fraction would not by itself justify neglecting neutrino gravity if neutrinos clumped far more strongly than CDM — the local contribution to \(\Phi\) is what matters near a halo. But the opposite holds: their large thermal velocities keep them smooth. Fully nonlinear simulations find the neutrino overdensity stays of order unity, \(\delta_\nu\sim O(1)\), exactly where the CDM overdensity runs to \(\delta_{\rm CDM}\sim500\text{–}1000\) inside collapsed halos. So neutrinos are subdominant both in the cosmic mean and, more importantly, in the dense regions where gravity is strongest.

The potential is sourced by total matter; split the source into CDM and neutrino pieces. Because \(\rho_\nu\ll\rho_{\rm CDM}\) everywhere in the regime of interest, the neutrino contribution to \(\Phi\) is negligible term-by-term:

Concretely: we keep the force neutrinos feel in full, and discard only the force they exert. The back-reaction of neutrinos onto CDM is dropped as well, which is why the CDM can be evolved on its own beforehand.

Note carefully what the approximation requires: that the neutrino density non-uniformity is small relative to that of CDM — not that the neutrino perturbation itself is small. This is a different and weaker condition than the one behind a standard perturbative expansion (which needs the perturbation small compared to the background). It holds robustly because \(M_\nu/M_{\rm CDM}\lt1\%\), and it would break only for exotic scenarios — much heavier neutrinos, or beyond-Standard-Model hot components that cluster far more strongly. Within the physical regime, this is one of the safest approximations in the whole pipeline.

The CDM distribution is computed separately and in advance, by a standard classical N-body simulation of the cold component alone. From that run one reads off the gravitational force per unit mass felt by a neutrino at every grid point and every time (paper Eq. 5),

This precomputed field is then tabulated — later loaded into a QRAM with \(O(n_{\rm gr}^{3})\) entries (one 3-vector per spatial cell, per time slice) so the quantum circuit can query it.

Crucially, the CDM N-body simulation is the cheap, well-understood half of the problem. CDM is cold, so its velocity spread is tiny and N-body is both accurate and shot-noise-tolerant for it (Part 2). It lives in 3-D configuration space, not 6-D phase space, so it never touches the \(n_{\rm gr}^6\) curse. The paper thus offloads all the hard, nonlinear self-gravitating dynamics onto the component that classical methods already handle well, and reserves the quantum machinery for the one component they handle badly — the neutrinos.

Because \(\mathbf F_{\rm CDM}\) is fixed by CDM alone, it does not depend on the neutrino distribution \(f\): it is a prescribed, known function of \((t,\mathbf x)\) only. The neutrino Vlasov equation therefore becomes a passive transport problem — neutrinos are pushed around by a given external field they do not influence, like tracer particles in a prescribed flow. All the feedback that made the general problem nonlinear has been severed.

Earlier quantum-Vlasov work, rooted in plasma physics, kept the force self-consistent: the charged particles source their own field through Poisson's equation, so \(\mathbf F\) depends on \(f\) and the system is nonlinear. Those methods then need extra machinery (e.g. Carleman linearization) with application conditions and weaker guarantees. Here the force is external and given from the outset — a genuinely different, and far more tractable, setting. This structural simplification is exactly what allows an efficient, unconditional quantum algorithm rather than a conditional one.

Replace \(\mathbf F(t,\mathbf x)\) by the prescribed \(\mathbf F_{\rm CDM}(t,\mathbf x)\) in the Vlasov equation:

Structurally it looks identical to the original Vlasov equation — the same three terms — but one thing has changed decisively: the coefficient \(\mathbf F_{\rm CDM}\) is now a fixed, externally given function, no longer a functional of \(f\).

Check linearity directly. \(\partial_t f\): one power of \(f\). \(\mathbf v\cdot\partial_{\mathbf x}f\): the coefficient \(\mathbf v\) is an independent phase-space coordinate, so this too is one power of \(f\). \(\mathbf F_{\rm CDM}(t,\mathbf x)\cdot\partial_{\mathbf v}f\): the coefficient \(\mathbf F_{\rm CDM}\) is a given function of \((t,\mathbf x)\), not of \(f\), so this is again exactly one power of \(f\). Every term contains \(f\) to the first power and no higher — the definition of a linear equation.

The contrast makes the point sharp. In the full self-gravitating case the force obeyed \(\mathbf F\propto-\nabla\Phi[\rho_\nu]\) with \(\rho_\nu=\int f\,d^3v\), so the forcing term was \(\big(\!-\nabla\Phi[\textstyle\int f\,d^3v]\big)\cdot\partial_{\mathbf v}f\) — two factors of \(f\), a quadratic nonlinearity. Freezing the force to the external \(\mathbf F_{\rm CDM}\) removes the inner \(f\)-dependence, collapsing the quadratic term to a linear one. That is the entire content of the linearization.

A linear PDE with given (possibly time-dependent) coefficients becomes, after discretization, a linear ODE system \(\dot{\mathbf f}=A(t)\,\mathbf f\) (Part 4). Linear ODEs are exactly what map cleanly onto Hamiltonian simulation, because unitary quantum evolution \(i\partial_t|\psi\rangle=H|\psi\rangle\) is itself linear. A surviving nonlinear term would have blocked this route outright — nonlinear quantum-PDE methods must fall back on Carleman linearization, embedding the problem in a larger space with approximation-dependent, conditional guarantees. By contrast this equation is natively linear, so the optimal-complexity Hamiltonian-simulation machinery applies without caveat.

Write the evolution as \(\partial_t f=-\mathcal L f\) and split the generator into two operators,

The two correspond exactly to the two non-trivial terms of the linearized Vlasov equation, and they act on different registers of phase space — one on position, one on velocity — which is what lets them be handled almost independently.

\(\mathcal L_{\rm stream}=\mathbf v\cdot\partial_{\mathbf x}\) is a velocity-weighted spatial derivative: it moves the distribution through position, and the speed of that motion at each point is the local velocity \(\mathbf v\). It is time-independent — the streaming operator never changes during the whole simulation — because \(\mathbf v\) is just a coordinate. It couples the position register to the velocity register (the velocity supplies the coefficient) but differentiates only in position.

\(\mathcal L_{\rm force}=\mathbf F_{\rm CDM}(t,\mathbf x)\cdot\partial_{\mathbf v}\) is an \(\mathbf x\)-dependent velocity derivative: it slides the distribution through velocity, at a rate set by the local gravitational force. It carries all the gravitational information, and it is time-dependent, since \(\mathbf F_{\rm CDM}(t,\mathbf x)\) evolves as the CDM structure grows. This is the operator that will demand the QRAM, because its coefficients are the precomputed force field.

Two properties of this decomposition are what make the quantum implementation efficient. First, sparsity: after central-differencing (Part 4) each derivative couples only nearest-neighbour grid points along one axis, so each operator has \(O(1)\) nonzeros per row — an \(s\)-sparse matrix, which is cheap to block-encode (Part 6). Second, the clean separation onto position vs velocity registers means the two can be Trotterized as separate, well-conditioned factors \(e^{-i\Delta_t H}\approx e^{-i\Delta_t H_{\rm stream}}e^{-i\Delta_t H_{\rm force}}\). Sparsity gives cheap block-encodings; the split gives a clean Trotterization — together they are the reason the evolution avoids the \(n_{\rm gr}^6\) cost.

Confine each position coordinate to \([0,L]\) with periodic boundaries, and each velocity coordinate to \([-V,V]\) with Dirichlet boundaries (\(f=0\) at \(\pm V\)). Physically: the box wraps in space, and no neutrino is faster than \(V\) (choose \(V\) large enough that \(f\) already vanishes there).

with \(i\in\{0,\dots,n_{\rm gr}-1\}\), and analogously for \(y,z\) and \(v,w\). Note \(n_{\rm gr}=2^{m_{\rm gr}}\) is a power of two — so each axis indexes exactly onto \(m_{\rm gr}\) qubits.

A phase-space cell is labelled by six indices \(\mathbf i=(i_x,i_y,i_z,i_u,i_v,i_w)\), flattened to a single integer \(i=\sigma(\mathbf i)\in\{0,\dots,N_{\rm gr}-1\}\) with \(N_{\rm gr}=n_{\rm gr}^{6}\) (paper Eqs. 13–14). The sampled distribution becomes a vector \(\mathbf f\in\mathbb R^{N_{\rm gr}}\) with entries (paper Eq. 17)

with \(\mathbf e_x=(1,0,0)\), and identically for \(y,z\) and \(u,v,w\). A derivative at a point uses only its two nearest neighbours along that axis.

• Sparsity. Each grid point couples to \(O(1)\) neighbours per axis → the resulting matrix has only \(O(1)\) nonzeros per row (\(s\)-sparse). Cheap to block-encode.
• Antisymmetry (previewed). The stencil’s \(+\) neighbour and \(-\) neighbour enter with opposite signs \(\pm 1/(2\Delta)\). This \(\pm\) pairing is precisely what makes the difference operator antisymmetric — proven two slides on.

An upwind stencil would be more stable/positivity-preserving but is not antisymmetric — it would spoil the Hermitian structure and the norm conservation the algorithm relies on. The paper deliberately uses central differences to keep \(A\) antisymmetric, and notes positivity as an open problem to revisit.

Streaming gives the (time-independent) \(A_{x,y,z}\); forcing gives the (time-dependent) \(A_{u,v,w}(t)\).

The \(x\)-streaming term \(v_x\,\partial_x\) becomes a tensor product acting on the \(x\)-position register and the \(u\)-velocity register:

where \(D_{\rm per}\) is the \(n_{\rm gr}\times n_{\rm gr}\) periodic central-difference matrix (\(+1\) on the super-diagonal, \(-1\) on the sub-diagonal, and \(\mp1\) in the wrap-around corners), and \(E_u=\mathrm{diag}\!\big(u_{i}/(2\Delta_x)\big)\) supplies the velocity weight \(v_x\).

The \(u\)-forcing term \(F_{{\rm CDM},x}\,\partial_u\) gives a matrix whose only nonzero entries connect a cell to its \(\pm u\)-neighbours, weighted by the local force:

Again the \(\pm\) sign pairing — the seed of antisymmetry.

Take \(D_{\rm per}\). Its entries are \((D_{\rm per})_{i,i+1}=+1\), \((D_{\rm per})_{i,i-1}=-1\) (indices mod \(n_{\rm gr}\)), zero on the diagonal. Transposing swaps rows and columns: \((D_{\rm per}^{\mathsf T})_{i,i+1}=(D_{\rm per})_{i+1,i}=-1=-(D_{\rm per})_{i,i+1}\). The same holds for every nonzero entry, so \(D_{\rm per}^{\mathsf T}=-D_{\rm per}\). The diagonal weight \(E_u\) is symmetric, and \((X\otimes Y)^{\mathsf T}=X^{\mathsf T}\otimes Y^{\mathsf T}\); so \(A_x^{\mathsf T}=(-D_{\rm per})^{\mathsf T}\!\otimes\dots=-A_x\).

For \(A_u\), the \(\mathbf i\to\mathbf i+\mathbf e_u\) entry is \(-F/2\Delta_v\) and the reverse \(\mathbf i+\mathbf e_u\to\mathbf i\) entry is \(+F/2\Delta_v\) with the same force \(F\) (which depends only on position, common to both). These are negatives of each other ⇒ \((A_u)^{\mathsf T}=-A_u\). Summing antisymmetric matrices preserves antisymmetry, so \(A(t)=\sum_a A_a\) is antisymmetric. \(\blacksquare\)

For \(\dot{\mathbf f}=A\mathbf f\) with \(A^{\mathsf T}=-A\):

So \(\|\mathbf f(t)\|\) is constant: the evolution is a rotation (orthogonal flow). There is no numerical blow-up — a real merit of the central-difference choice (paper, around Eq. 23).

The physical particle count is \(f_{\rm sum}=\sum_{\mathbf i}f_{\mathbf i}\). A short calculation (paper Eq. 24) shows \(\dot f_{\rm sum}\) reduces to boundary terms in the velocity coordinates, which vanish because \(f=0\) at \(\pm V\). Hence the number of neutrinos in the box is conserved — escapes only occur if the box is too small.

These two facts (constant norm, conserved count) are the discrete shadows of Liouville’s theorem. They are also exactly what let us reinterpret the ODE as a Schrödinger equation — the next slide — because unitary quantum evolution also conserves norm. The physics and the quantum formalism agree because both are norm-preserving.

A Hermitian matrix satisfies \(H^\dagger=H\), where \(\dagger\) is conjugate-transpose. With \(A\) real (\(A^*=A\)) and antisymmetric (\(A^{\mathsf T}=-A\)):

So \(H\) is a legitimate quantum Hamiltonian: its eigenvalues are real, its evolution unitary.

If \(A\) had any symmetric part, \(H=iA\) would be non-Hermitian, its evolution non-unitary, and the norm would grow or decay — no quantum Hamiltonian, and the algorithm’s optimal Hamiltonian-simulation machinery would not apply. The whole bridge rests on the antisymmetry proved in slide 4.11–4.12.

Start from \(\dot{\mathbf f}=A\mathbf f\). Since \(A=-iH\),

(Setting \(\hbar=1\).) This is the Schrödinger equation, with the discretized neutrino distribution playing the role of the “wavefunction”.

Nothing about the neutrinos is quantum-mechanical in this equation — it is a classical kinetic PDE. What is true is that its mathematical form is identical to quantum evolution. That formal identity is the entire licence to use a quantum computer: a device built to solve \(i\partial_t|\psi\rangle=H|\psi\rangle\) will solve this too.

For piecewise-constant \(H\) (Part 6), the solution is a product of unitaries (paper Eq. 34),

Implementing each factor \(e^{-i\Delta_t H_{i_t}}\) as a circuit is the task called Hamiltonian simulation.

Store the entries of \(\mathbf f\) as the amplitudes of a quantum state over the grid-index basis \(\{|\mathbf i\rangle\}\):

The basis state \(|\mathbf i\rangle=|i_x\rangle|i_y\rangle|i_z\rangle|i_u\rangle|i_v\rangle|i_w\rangle\) is the binary label of the phase-space cell, using \(\lg N_{\rm gr}=6\,m_{\rm gr}\) qubits.

• Normalization. The state is \(\mathbf f/\|\mathbf f\|\); the overall scale \(\|\mathbf f\|\) is not stored in the state. Because the evolution conserves \(\|\mathbf f\|\) (slide 4.13), this is a harmless constant, tracked classically.
• Readout is limited. You cannot look at the \(2^n\) amplitudes directly; extracting information needs a dedicated procedure (Part 7). Encoding is cheap; decoding is the bottleneck.

The initial state \(|f(0)\rangle\) is prepared by an oracle \(O_{f(0)}\) (paper Eq. 30). For the cosmological initial condition — a Fermi–Dirac velocity distribution modulated by the primordial density perturbation (paper Eqs. 61–64) — this is built from arithmetic circuits + a QRAM holding the initial \(\delta_\nu(0,\mathbf x)\).

\(n\) qubits span a Hilbert space of dimension \(2^{n}\). Encoding \(N_{\rm gr}=n_{\rm gr}^{6}\) grid amplitudes therefore needs only

versus \(N_{\rm gr}=n_{\rm gr}^{6}\) classical numbers. The \(6\)-D curse becomes a logarithmic qubit count — an exponential compression of storage.

Cheap storage would be worthless if evolving or reading the state cost \(O(N_{\rm gr})\). The algorithm’s real work (Parts 6–7) is doing the evolution \(e^{-iHT}\) and the read-out \(P_\nu(k)\) without ever touching all \(2^n\) amplitudes — using sparsity (block-encoding) and amplitude estimation. Only then is the memory win an actual speed win.

Storage: exponential win, clean. Still owed: (i) evolve the state efficiently despite time-dependent \(H\); (ii) load the external force \(\mathbf F_{\rm CDM}\) into the circuit (QRAM); (iii) extract \(P_\nu(k)\) from an un-readable state. These are Parts 6, 6, and 7 respectively.

Given the initial state \(|f(0)\rangle\) and the Hamiltonian \(H(t)=iA(t)\), produce

where \(\mathcal T\) is time-ordering. This is the problem called Hamiltonian simulation — historically the first proposed use of a quantum computer (Feynman, Lloyd).

The external force is known only at discrete times (from the N-body run), so \(H(t)\) is taken piecewise constant on \(n_t\) intervals of length \(\Delta_t=T/n_t\) (paper Assumption 1, Eq. 29). Then time-ordering collapses to a plain product (paper Eq. 34),

So the whole job reduces to implementing one factor \(e^{-i\Delta_t H_{i_t}}\) for a constant Hermitian \(H_{i_t}\), then repeating \(n_t\) times.

Two ways to build \(e^{-i\Delta_t H}\): the elementary Trotter product formula (6.4–6.6), and the near-optimal block-encoding + qubitization route the paper actually uses for its complexity claim (6.7–6.13). We do both, because the Trotter route is what we could run on today’s hardware (Part 10).

Split \(H=H_{\rm stream}+H_{\rm force}\) (the streaming and forcing blocks of Part 3.10). In general \([H_{\rm stream},H_{\rm force}]\neq0\), so \(e^{-iH\tau}\neq e^{-iH_{\rm stream}\tau}e^{-iH_{\rm force}\tau}\) — you cannot exponentiate the pieces separately and multiply.

But for a small time slice \(\tau=t/n\) the error is controlled:

Each factor is easy: \(H_{\rm stream}\) is diagonalized by a Fourier transform on the position registers; \(H_{\rm force}\) by one on the velocity registers. So one Trotter step = a few standard gates.

The error falls as \(1/n\) (first order; symmetric “Strang” splitting gives \(1/n^2\)). More slices → more accurate but a longer circuit. On real hardware more gates means more noise — the exact trade-off our QV-F3/F4 charts measured (Part 10), where the sweet spot was a handful of steps.

Quantum gates are unitary (norm-preserving, invertible). But \(H\) is a general Hermitian matrix, not unitary. To manipulate \(H\) with a circuit we embed it as a sub-block of a larger unitary.

A unitary \(U_H\) on \(a+n\) qubits is an \((\alpha,a)\)-block-encoding of \(H\) if its top-left \(2^n\times2^n\) block equals \(H/\alpha\):

The \(a\) ancilla qubits are the “address” of the block; \(\alpha\ge\|H\|\) is a normalization. Acting with \(U_H\) and post-selecting the ancillas on \(|0\rangle\) applies \(H/\alpha\).

Once you have a block-encoding of \(H\), a whole calculus (quantum signal processing / QSVT) lets you build a block-encoding of almost any function \(g(H)\) — including the one we want, \(e^{-iH\tau}\) — at near-optimal cost. The remaining question is: how do you build \(U_H\) for our \(H\)? Answer: from its sparsity (next slide).

A matrix is \(s\)-sparse if every row and column has at most \(s\) nonzero entries. Our \(A\) (hence \(H\)) is \(s\)-sparse with a small constant \(s\): a derivative couples a grid cell only to its immediate neighbours (slide 4.5). The nonzero positions are fixed by the stencil; the values involve \(v\) and \(\mathbf F_{\rm CDM}\).

A sparse block-encoding is built from two reversible (unitary) look-up circuits:

• Position oracle \(O_A\): given a row \(\mathbf i\) and a neighbour index \(\ell\), returns the column of the \(\ell\)-th nonzero — i.e. “which cells does \(\mathbf i\) couple to.” Computed on the fly from the stencil.
• Value oracle \(O_H\): given \((\mathbf i,\mathbf j)\), returns the matrix entry \(H_{\mathbf i\mathbf j}\). For the streaming block this is arithmetic in \(v\); for the forcing block it needs the external force \(\mathbf F_{\rm CDM}(t,\mathbf x)\) — which must be looked up from a QRAM (next slide).

Given \(O_A,O_H\), a standard construction (paper’s Theorem 3) yields a block-encoding of \(H\) with \(O(1)\) oracle queries and \(O(\mathrm{poly\,log})\) extra gates. Sparsity is exactly why the \(n_{\rm gr}^6\)-dimensional \(H\) never has to be written down — you only ever query a few of its entries.

Given classical data \(\mathcal X=\{x_0,\dots,x_{N-1}\}\), a QRAM implements the coherent look-up

i.e. it returns data for a superposition of addresses at once. Here the data is the precomputed \(\mathbf F_{\rm CDM}(t,\mathbf x)\) on the \(n_{\rm gr}^3\) position grid points.

The QRAM stores \(O(n_{\rm gr}^{3})\) real numbers (position grid only — velocity is handled analytically). That is far smaller than the classical \(O(n_{\rm gr}^{6})\) full phase-space storage — a genuine improvement — but still a large, structured quantum memory. One QRAM is loaded per time step (piecewise-constant force), and the same device can be reused for different neutrino masses.

A scalable QRAM does not yet exist. Circuit-based QRAM proposals need \(O(N)\) gates (though depth can be logarithmic), and error correction for QRAM is an open problem. The paper is explicit: “implementing a QRAM is highly challenging, and we simply assume its availability.” This is caveat #1 from slide 0.6.

From a block-encoding \(U_H\) of \(H\), quantum signal processing (Low–Chuang “qubitization”) constructs a new circuit whose action on the encoded subspace is \(g(H)\) for a chosen polynomial \(g\). Choosing \(g\) to approximate \(e^{-i\alpha\tau x}\) on \([-1,1]\) yields a block-encoding of \(e^{-iH\tau}\).

The number of calls to \(U_H\) needed to simulate for time \(\tau\) is (paper’s Theorem 4)

which is optimal: it matches the “no-fast-forwarding” lower bound (you cannot simulate for time \(\tau\) with sub-linear-in-\(\tau\) work in general). The additive \(\log(1/\epsilon)\) is the exponentially cheap price of accuracy.

The cost scales with \(\|H\|\tau\). For our operator \(\|H\|\sim\) (velocity)/\(\Delta_x\) + (force)/\(\Delta_v\), so \(\|H\|\) grows with the grid resolution \(n_{\rm gr}\). Bounding this product — by choosing the box \(L,V\) sensibly — is what keeps the total query count at \(\widetilde O(n_{\rm gr})\) (Part 8).

By Assumption 1 the force is constant on each of \(n_t\) intervals, so \(H(t)=H_{i_t}\) is constant there (paper Eq. 29). Each interval is simulated independently by qubitization with its own block-encoding (its own QRAM of \(\mathbf F_{\rm CDM}\) for that time), and the results are chained (paper Eq. 34, Fig. 1).

Assembling the \(n_t\) steps, the circuit \(U_{f(T),\epsilon}\) on \((2\lg N_{\rm gr}+5)\) qubits produces \(|f(T)\rangle\) to accuracy \(\epsilon\), querying the force oracles (paper Eq. 32)

The first term is the streaming/forcing cost (\(\|H\|T\)); the second is the overhead of chaining \(n_t\) steps.

Naively each time step needs its own QRAM (\(n_t\) of them). The paper notes you can reuse one QRAM’s qubits, re-loading it per interval, so the qubit count does not multiply by \(n_t\). This is an implementation nicety that keeps the space cost logarithmic.

Each evolution block calls the force oracles \(O_{F_{\rm CDM}}^{i_t}\) (fed by QRAM). The theorem guarantees the output is \(|f(T)\rangle\) up to a small garbage component \(\||\psi_{\rm gar}\rangle\|\le\epsilon\) (paper Eq. 31).

We now hold \(|f(T)\rangle\) — the full evolved neutrino distribution — in \(O(\log n_{\rm gr})\) qubits, built with \(\widetilde O(n_{\rm gr}+n_t)\) oracle queries. But it is unreadable: measuring gives one random grid cell, not \(P_\nu(k)\). Extracting the physics is Part 7.

Integrate \(f\) over velocity to get the neutrino density, then its fractional perturbation (paper Eqs. 7–8, 40–41):

Fourier-transform \(\delta_\nu\) on the grid, \(\tilde\delta^\nu_{\mathbf k}=\frac{1}{n_{\rm gr}^3}\sum_{\mathbf x}\delta^\nu_{\mathbf x}\,e^{i\mathbf k\cdot\mathbf x}\), and define

So the target is the squared magnitude of a single Fourier mode of the neutrino density (then ensemble-averaged). This is a scalar — a tiny summary of the enormous state \(|f(T)\rangle\).

Reading the whole \(f\) (\(2^n\) amplitudes) is impossible in \(\mathrm{poly}(n)\) time — it would destroy the exponential advantage. The algorithm is honest about this: it targets the specific low-dimensional quantity cosmologists actually compare to data, and extracts only that.

Measuring \(|f(T)\rangle=\sum_{\mathbf i}f_{\mathbf i}|\mathbf i\rangle\) in the computational basis returns a single grid cell \(\mathbf i\) with probability \(|f_{\mathbf i}|^2/\|\mathbf f\|^2\). One shot gives one random cell; you learn essentially nothing about a specific Fourier amplitude, and repeating \(2^n\) times to reconstruct \(f\) is exponentially expensive.

The fix: apply a unitary \(W\) that maps the quantity we want onto the amplitude of one specific basis state, then estimate that amplitude. Concretely (paper Eq. 46) \(W\) is a quantum Fourier transform on the position registers and a Hadamard/sum on the velocity registers,

The QFTs turn position into wavenumber \(\mathbf k\); the Hadamards sum over velocities (that is the \(\int d^3v\) that makes the density).

After \(W\), the overlap with the target basis state \(|\Omega_{\mathbf k}\rangle=|\mathbf k\rangle|0\rangle_{\rm vel}\) is

with a known constant \(C\) (paper Eqs. 42, 50). So the Fourier power we want is exactly the probability of finding the system in \(|\Omega_{\mathbf k}\rangle\) — a quantity amplitude estimation can measure efficiently.

Given a circuit that prepares \(|\Phi\rangle=\sqrt{p}\,|\text{good}\rangle+\sqrt{1-p}\,|\text{bad}\rangle\), estimate the amplitude \(p\). Naive sampling needs \(O(1/\epsilon^2)\) repetitions for precision \(\epsilon\) (shot noise). QAE (Brassard et al.; via amplitude amplification + phase estimation) needs only

— a quadratic speedup in the read-out (paper’s Theorem 5).

Set “good” \(=|\Omega_{\mathbf k}\rangle\). The preparation circuit is \(W\,U_{f(T)}\). QAE estimates \(p=|\langle\Omega_{\mathbf k}|W|f(T)\rangle|^2=C|\tilde\delta^\nu_{\mathbf k}|^2\), from which \(|\tilde\delta^\nu_{\mathbf k}|^2\) follows by dividing the known \(C\). To reach accuracy \(\epsilon\) on \(|\tilde\delta|^2\), QAE is run at precision \(\sim C\epsilon\).

Estimating one Fourier mode to accuracy \(\epsilon\) with confidence \(1-\delta\) costs (paper Eq. 43)

i.e. the state-generation cost times the QAE factor \(1/(C\epsilon)\). Still only linear in \(n_{\rm gr}\).

\(P_\nu(k)=\langle|\tilde\delta_{\mathbf k}|^2\rangle\) is an average over random initial conditions from inflation. Classically that means running the whole simulation \(n_{\rm IV}\) times with different seeds and averaging. Quantum-mechanically, you prepare a superposition of the different initial conditions in one extra register and estimate the average with a single QAE — no repeated runs.

Ensemble averaging is intrinsic to the observable, and superposition handles it natively. This is a clean example of quantum parallelism doing real work — one of the paper’s nicer conceptual points, distinct from the raw Hamiltonian-simulation speedup.

A classical grid solver stores the sampled distribution \(f\) on \(N_{\rm gr}=n_{\rm gr}^{6}\) phase-space cells and advances it \(n_t\) time steps; each explicit step touches every cell (a handful of neighbour reads per cell, so an \(O(1)\)-per-cell stencil). The two costs are therefore

The exponent is not an artifact of a clumsy method — it is the dimension of phase space itself. Six axes, each gridded \(n_{\rm gr}\) ways, multiply into a sixth power that no clever data structure removes.

Accurate neutrino work is doubly punishing: the velocity distribution is broad and shifting, so it demands many velocity grid points, and the six axes then multiply. Put in one concrete number: a modest \(n_{\rm gr}=100\) already gives \(N_{\rm gr}=10^{12}\) cells, which at 8 bytes each is \(\sim\!8\) TB for a single snapshot — and an explicit scheme needs at least two buffers. Advancing that for \(n_t\sim10^{2}\text{–}10^{3}\) steps is \(10^{14}\text{–}10^{15}\) cell-updates per run. This is why state-of-the-art classical Vlasov solvers manage only hundreds of position points and tens of velocity points, and why cosmologists retreat to N-body (and swallow its shot noise) instead of solving the PDE directly.

Worse, \(P_\nu(k)\) is an ensemble average over the random initial conditions seeded by inflation, not a single deterministic output. Classically you estimate it by running \(n_{\rm IV}\) independent realizations (different random seeds) and averaging — so the true bill carries a third factor, \(O(n_{\rm gr}^{6}\,n_t\,n_{\rm IV})\). The quantum algorithm folds this ensemble into a single superposed run (slide 7.11), erasing the \(n_{\rm IV}\) factor outright.

So "beating classical" here means beating \(O(n_{\rm gr}^{6})\) in memory, \(O(n_{\rm gr}^{6}n_t)\) in evolution time, and \(O(n_{\rm IV})\) in ensembling — all at once, for the specific deliverable \(P_\nu(k)\) to accuracy \(\epsilon\). Keep this triple in view: the next slide meets it on memory and grid-scaling cleanly, and the slide after that is honest about the fine print.

The simulation domain has to be just large enough that no neutrino escapes during the run of duration \(T\). A particle drifts \(\lesssim VT\) in space and is kicked \(\lesssim F_{\max}T\) in velocity, so the box half-widths \(L\) (space) and \(V\) (velocity) must satisfy

This pins the physical scales to \(T\) and feeds directly into the norm \(\|H\|\), which is what the query count depends on.

Substituting the box scaling into the paper's Theorems 1–2, the query complexity to output \(P_\nu(k)\) to accuracy \(\epsilon\), together with the qubit count, are

The tilde in \(\widetilde O\) hides polylogarithmic factors and the \(1/\epsilon\) amplitude-estimation read-out factor, which are real but sub-dominant to the leading grid dependence.

Optimal (qubitization) Hamiltonian simulation of an \(s\)-sparse \(H\) costs \(\widetilde O(s\,\|H\|\,T)\) queries to the block-encoding. The central-difference stencil is \(O(1)\)-sparse, so \(s\) is a constant. The derivative weight \(1/\Delta\) with \(\Delta\sim L/n_{\rm gr}\) makes \(\|H\|\) grow like \(n_{\rm gr}\), while the box choice keeps \(T\sim L/V\) fixed — so schematically \(\|H\|T\sim n_{\rm gr}\), giving \(\widetilde O(n_{\rm gr})\) for the evolution, plus \(\widetilde O(n_t)\) for stepping the time-dependent force. Contrast the classical \(n_{\rm gr}^{6}\): the sixth power collapses to first power because the state lives in amplitudes, not cells.

The amplitude-encoded state needs only \(n=6\log_2 n_{\rm gr}\) qubits — one \(m_{\rm gr}\)-bit register per axis, already logarithmic in the grid. Qubitization, the sparse-access oracle, and the nested QAE add ancilla registers whose width grows with the target accuracy \(\epsilon\); carrying \(\epsilon\) through the layered primitives raises the logarithm to the \(5/2\) power, \(O(\log^{5/2}(n_{\rm gr}/\epsilon))\). The key point survives the bookkeeping: the qubit count is polylogarithmic in the grid, an exponential compression of the \(n_{\rm gr}^{6}\) classical memory.

Side by side: grid dependence falls from classical \(n_{\rm gr}^{6}\) to quantum \(\widetilde O(n_{\rm gr})\) queries on \(O(\mathrm{polylog})\) qubits — an exponential improvement in memory and a large polynomial improvement in the grid scaling of the run-time. This one sentence is the entire justification the paper is built to earn.

• Memory: \(O(\log n_{\rm gr})\) qubits vs \(O(n_{\rm gr}^6)\) classical numbers — unconditional, a direct consequence of amplitude encoding (\(n=6\log_2 n_{\rm gr}\) qubits hold all \(n_{\rm gr}^6\) amplitudes).
• Grid scaling of evolution: \(\widetilde O(n_{\rm gr})\) queries via optimal Hamiltonian simulation on a sparse, antisymmetric \(H\) — the sixth power is genuinely gone, not hidden.
• Ensemble average: one superposed run replaces \(n_{\rm IV}\) seeded classical runs (slide 7.11).

• QRAM for \(\mathbf F_{\rm CDM}\) (\(O(n_{\rm gr}^{3})\) stored entries): the whole \(\widetilde O(n_{\rm gr})\) count assumes each force query costs \(O(\mathrm{polylog})\). No hardware delivers that at scale today; if QRAM access is not polylogarithmic, the advantage erodes.
• Output is only \(P_\nu(k)\): the \(1/\epsilon\) QAE factor is real, and you cannot cheaply extract the full \(2^n\)-amplitude \(f\). That is a fine match for cosmology (we only ever wanted \(P(k)\)), but it is a genuine restriction, not a free lunch.
• Fault tolerance: block-encoding, qubitization, and QAE are deep, high-precision circuits that assume an error-corrected machine — categorically not runnable on today's NISQ hardware.

Every count above is in the query model: it counts calls to the block-encoding and force oracles, not raw gate count. A single query itself compiles to many gates, and \(\widetilde O\) suppresses polylog and \(1/\epsilon\) factors. None of this cancels the exponential memory win or the polynomial grid-scaling win — but it means "\(\widetilde O(n_{\rm gr})\)" is an asymptotic, oracle-level statement, and the constant prefactors on a fault-tolerant device would be large.

What the paper actually establishes is a rigorous, asymptotic, conditional speedup for a well-posed sub-problem — linearized neutrino Vlasov, output = power spectrum — contingent on fault tolerance and scalable QRAM. That is promising theory with cleanly labelled open problems, not a demonstrated capability. Holding those two apart is the discipline of this whole deck.

Strip the six axes down to one position and one velocity coordinate \((x,u)\), on a \(64\)-point grid (\(n_{\rm gr}=64\)), with spatial box \(L=2\) and velocity range \(V=1\). The externally-supplied CDM force is the simplest nontrivial field, a single spatial sinusoid,

and the initial distribution is uniform in \(x\) and Maxwellian in \(u\) (paper Eq. 70) — a calm, structureless gas that any structure in the final state must have been generated by the force.

Each choice is deliberate. A single Fourier mode \(\sin(Kx)\) means the analytic answer is predictable — a linear equation driven at wavenumber \(K\) should respond at \(K\) and nowhere else — giving a sharp pass/fail test. The Maxwellian initial \(u\)-profile is the honest thermal distribution of a warm species, and the uniform \(x\)-profile guarantees the initial density contrast is exactly zero, so the emergent \(\delta_\nu(x)\) is unambiguously the force's doing.

Assemble the discretized antisymmetric operator \(A\) exactly as in Eqs. 16–22, then evolve classically by the dense matrix exponential,

This \(e^{AT}\) is precisely the operation the quantum circuit would realize by Hamiltonian simulation (\(e^{-iHT}\) with \(H=iA\)); here it is computed on an ordinary laptop, at a size small enough to exponentiate directly, purely to check that the physics is right.

It genuinely exercises the discretization, the antisymmetry \(A^{\mathsf T}=-A\), and the coupled streaming↔forcing dynamics — the entire classical core of the reduction. It does not run any quantum circuit, invoke no block-encoding, no QRAM, no QAE, and demonstrates no speedup. It is a correctness check of the reduction, not of the quantum algorithm — a distinction the next two slides keep sharp.

Three heat-maps of \(f(T,x,u)\) at \(T=0,0.1,0.2\): the horizontal axis is position \(x\), the vertical axis is velocity \(u\), and the colour is the phase-space density of neutrinos at that \((x,u)\). Reading a vertical slice gives the local velocity distribution at a point; a horizontal slice gives the spatial profile at fixed speed.

• \(T=0\): a flat horizontal band — uniform in \(x\), Maxwellian in \(u\). The calm initial gas, structureless in space.
• \(T\gt 0\): the band wrinkles. Streaming \((u\,\partial_x)\) shears it horizontally — fast neutrinos (\(u\gt 0\)) slide right, slow ones (\(u\lt 0\)) left, so the band tilts. Simultaneously the force \((F\,\partial_u)\) pushes density up or down in \(u\), hardest where \(|\sin Kx|\) peaks. The two effects together tilt and fold the band into the characteristic shear.

The evolution is a rigid rotation of the state vector (antisymmetric generator), so phase-space area is preserved and the band cannot simply spread and blur — it must wind. Over cosmological times this same winding is what drives phase mixing and, eventually, the fine filamentary structure of a collisionless species. The toy captures the first \(0.2\) units of that story in a single readable frame.

Throughout the run \(\|\mathbf f\|\) stays constant to machine precision — the discrete conservation law of slide 4.13, now confirmed numerically rather than merely proved. Constancy of the norm is the operational signature that the evolution is orthogonal (a rotation), exactly as the Schrödinger form \(i\partial_t|f\rangle=H|f\rangle\) with Hermitian \(H\) demands. Physics and formalism agree because both are norm-preserving.

Integrate the sheared \(f(T,x,u)\) over velocity to recover the spatial density \(\rho_\nu(x)=\int f(T,x,u)\,du\), then form the contrast \(\delta_\nu(x)=\rho_\nu/\bar\rho_\nu-1\). Starting from a uniform gas (\(\delta_\nu=0\)), a single smooth density wave emerges: neutrinos have been swept toward the convergence points of the force, building an overdensity where \(F\) compresses them and a deficit where it rarefies them.

Take the discrete Fourier transform \(\tilde\delta_\nu(k)\) and plot the power \(|\tilde\delta_\nu(k)|^2\). Only the \(k=K\) mode carries weight — every other mode is zero to numerical precision:

This is a direct consequence of linearity: a linear operator driven at a single wavenumber cannot manufacture other wavenumbers, so a single-mode force can only produce a single-mode density response. The demo confirms the discretization respects that exactly — no spurious mode leakage from the stencil.

The quantity \(|\tilde\delta^\nu_{\mathbf k}|^2\) plotted here is precisely what amplitude estimation would extract on hardware (slide 7.6) — the single Fourier component of the neutrino density that, aggregated over \(k\), is the power spectrum \(P_\nu(k)\). Computing it classically here shows the mechanism produces a well-defined, physically correct number for QAE to estimate.

The entire purpose of the QFT-based read-out operator \(W\) (slide 7.5) is to isolate one Fourier mode of the density. The demo shows the underlying physics really is single-mode for a single-mode force, so \(W\) is measuring a sharp, meaningful quantity rather than an artifact. Read-out design and physics are consistent — the last link in the reduction that the toy can check.

• The discretization is correct: \(A\) is antisymmetric, \(\|\mathbf f\|\) is conserved to machine precision, and the shearing dynamics are physically sensible.
• The forcing drives exactly the matching density mode (\(k=K\) in, \(k=K\) out) — so the full observable pipeline \(f\to\rho_\nu\to\delta_\nu\to\tilde\delta_\nu\) is right, with no spurious mode leakage.
• The reduction "Vlasov \(\to\) linear ODE \(\to\) Schrödinger" is faithful: \(e^{AT}\) is the exact classical shadow of the quantum \(e^{-iHT}\), and it produces the expected physics.

• No quantum circuit was run — no block-encoding, no QRAM, no qubitization, no QAE on any hardware or even an emulator.
• No speedup was demonstrated — this is a \(64\)-point, 1-D toy solved by dense matrix exponential on a classical computer, i.e. the very regime where classical is trivially adequate.
• The full 6-D cosmological problem, and the fault-tolerant resources (QRAM of \(n_{\rm gr}^3\) entries, deep QAE circuits) it would demand, remain entirely untouched.

The demo is a correctness check of the classical core of a quantum-algorithm proposal — genuinely valuable, and honestly presented as such. It confirms the mechanism the quantum circuit would run, but it is a long way from a neutrino power spectrum actually computed on a quantum computer. Closing part of that gap is exactly what our own reproduction (Part 10) set out to do.

We rebuilt the Sec. IV toy from scratch: the antisymmetric operator \(A\) on a \(64\)-point 1-D grid, evolved by the dense exponential \(e^{AT}\). We reproduced both published results — the phase-space shear (Fig. 2) and the single-mode density response (Figs. 3–4) — and confirmed the discrete conservation law numerically, with \(\|\mathbf f\|\) held to \(\sim\!10^{-15}\) (i.e. double-precision round-off). Matching the paper here is the ground truth against which everything quantum is checked.

We then built the actual Trotter circuit for \(e^{-iHT}\) — the gate-level object the paper describes but never runs — splitting the streaming and forcing operators and compiling each into gates. Executed on a noiseless statevector emulator at 4–6 qubits, it converges to the exact \(e^{-iHT}\) as the Trotter step count \(n\) grows, with the error falling like \(\propto 1/n\) exactly as first-order Trotter theory predicts. That is our QV-F3 result: the reduction is not just faithful on paper, it compiles to a circuit that provably approaches the target unitary.

This already goes past what Miyamoto et al. did. The paper stops at the abstract algorithm plus a classical matrix-exponential demo; we constructed and validated a concrete gate circuit for the mechanism, on the smallest honest instance, and verified its convergence. It is the bridge from "a theorem about a circuit" to "a circuit that runs" — and it is what made the hardware run on the next slide possible.

We transpiled the smallest honest instance — the 4-qubit circuit for \(n_{\rm gr}=4\) — and ran it on IonQ Forte-1 (a trapped-ion QPU), at 300 shots per circuit, in a three-way benchmark: exact classical \(\leftrightarrow\) noiseless simulator \(\leftrightarrow\) real hardware. When the device was well-calibrated, the measured output distribution reproduced the target to \(\sim\!90\text{–}95\%\) fidelity. The Vlasov Hamiltonian-simulation kernel genuinely runs on a quantum computer at toy scale — not in principle, in practice.

One result was counter-intuitive and demanded explanation: the shallowest circuit (fewest Trotter steps, fewest gates) was reproducibly the worst, which naïve depth-noise intuition gets backwards — fewer gates should mean less accumulated noise. A depth sweep combined with a depolarizing-noise model resolved it: the shallow-circuit error is coherent and systematic — a large-angle Trotter (algorithmic) error in the few-gate limit — not incoherent gate noise. The two regimes separate cleanly: coherent Trotter error dominates at low depth, hardware noise at high depth, and they cross at an optimum of 2 Trotter steps, the sweet spot.

This is an honest, toy-scale hardware demonstration of the algorithm's core kernel: Hamiltonian simulation of the discretized Vlasov state. It is emphatically not a run of the QRAM, the QAE read-out, or the claimed speedup — none of those touched hardware. It is the faithful smallest slice of a fault-tolerant proposal, executed on a NISQ device with the algorithmic and hardware error budgets cleanly disentangled.

Every row is a downgrade we made deliberately to get something honest onto hardware — and every row is a frontier that still has to be crossed for the speedup to be real.

• QRAM at cosmological scale — storing the CDM force with \(O(n_{\rm gr}^{3})\) entries and querying it in polylog time — does not exist; on hardware we simply hard-coded the toy force.
• Read-out is intrinsically limited to the summary \(P_\nu(k)\), at a \(1/\epsilon\) cost; the full \(2^n\)-amplitude \(f\) is never accessible cheaply.
• Fault tolerance is required for the deep block-encoding and QAE circuits — precisely what today's NISQ machines lack.

On top of these, the discretization's positivity of \(f\) is not guaranteed by central differences and is flagged as future work.

What stands at the end is a rigorous, elegant, conditional quantum-speedup proposal for a genuine cosmological sub-problem — with its core mechanism now demonstrated on real hardware at toy scale, and the road to usefulness mapped clearly and honestly. That combination — real theory, a real (if tiny) hardware run, and unflinching accounting of the gap — is what mastery of this paper looks like.

• Why a hot species needs a 6-D phase-space description (multi-streaming breaks any fluid closure), and how the six gridded axes produce the \(O(n_{\rm gr}^6)\) curse of dimensionality.
• Derive the Vlasov equation from the single premise "\(f\) is constant along trajectories," via the convective derivative and \(\dot{\mathbf x}=\mathbf v,\ \dot{\mathbf v}=\mathbf F\).
• Why neglecting neutrino self-gravity (\(\Omega_\nu/\Omega_m\lt 1\%\)) linearizes the equation by making \(\mathbf F_{\rm CDM}\) an external, prescribed field independent of \(f\).
• Set up the central-difference discretization to the linear ODE \(\dot{\mathbf f}=A(t)\mathbf f\), and say why central (not upwind) differences are chosen.
• Prove \(A^{\mathsf T}=-A\) from the \(\pm\) stencil pairing, hence \(H=iA\) is Hermitian and \(\|\mathbf f\|\) is conserved — the discrete shadow of Liouville's theorem.
• State the Schrödinger recast \(i\partial_t|f\rangle=H|f\rangle\) and the amplitude encoding, and locate the exponential memory win (\(6\log_2 n_{\rm gr}\) qubits for \(n_{\rm gr}^6\) amplitudes).
• Explain block-encoding, sparse-access oracles, and the QRAM that supplies \(\mathbf F_{\rm CDM}\) with \(O(n_{\rm gr}^3)\) entries.
• Explain qubitization (optimal \(e^{-iH\tau}\) on a sparse \(H\)) and how QAE reads out the single mode \(|\tilde\delta^\nu_{\mathbf k}|^2\) that builds \(P_\nu(k)\), at a \(1/\epsilon\) cost.
• Quote the complexity from memory: \(\widetilde O(n_{\rm gr}+n_t)\) queries on \(O(\mathrm{polylog})\) qubits, versus classical \(O(n_{\rm gr}^6 n_t)\) — and sketch where each factor comes from.
• Name the three caveats that make it conditional: QRAM at scale, output-only \(P(k)\) with the \(1/\epsilon\) factor, and fault tolerance for the deep circuits.

Companion decks carry the parts this tutorial only summarized: Deck · Jul 20 (our full reproduction, the IonQ hardware run, and the coherent-vs-depth-noise RCA) and QSim-Lit Survey (where this paper sits in the broader quantum-simulation literature). The original is arXiv:2310.01832, Phys. Rev. Research 6, 013200 (2024).